Line data    Source code

       1             : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
       2             :    Copyright (c) 2015-2023 The plumed team
       3             :    (see the PEOPLE file at the root of the distribution for a list of names)
       4             : 
       5             :    See http://www.plumed.org for more information.
       6             : 
       7             :    This file is part of plumed, version 2.
       8             : 
       9             :    plumed is free software: you can redistribute it and/or modify
      10             :    it under the terms of the GNU Lesser General Public License as published by
      11             :    the Free Software Foundation, either version 3 of the License, or
      12             :    (at your option) any later version.
      13             : 
      14             :    plumed is distributed in the hope that it will be useful,
      15             :    but WITHOUT ANY WARRANTY; without even the implied warranty of
      16             :    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
      17             :    GNU Lesser General Public License for more details.
      18             : 
      19             :    You should have received a copy of the GNU Lesser General Public License
      20             :    along with plumed.  If not, see <http://www.gnu.org/licenses/>.
      21             : +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
      22             : #include "SketchMapBase.h"
      23             : 
      24             : namespace PLMD {
      25             : namespace dimred {
      26             : 
      27          22 : void SketchMapBase::registerKeywords( Keywords& keys ) {
      28          22 :   DimensionalityReductionBase::registerKeywords( keys );
      29          22 :   keys.remove("NLOW_DIM");
      30          44 :   keys.add("compulsory","HIGH_DIM_FUNCTION","as in input action","the parameters of the switching function in the high dimensional space");
      31          44 :   keys.add("compulsory","LOW_DIM_FUNCTION","as in input action","the parameters of the switching function in the low dimensional space");
      32          44 :   keys.add("compulsory","MIXPARAM","0.0","the amount of the pure distances to mix into the stress function");
      33          22 : }
      34             : 
      35           6 : SketchMapBase::SketchMapBase( const ActionOptions& ao ):
      36             :   Action(ao),
      37             :   DimensionalityReductionBase(ao),
      38           6 :   smapbase(NULL),
      39           6 :   normw(0.0) {
      40             :   // Check if we have data from a input sketch-map object - we can reuse switching functions wahoo!!
      41           6 :   if( dimredbase ) {
      42           5 :     smapbase = dynamic_cast<SketchMapBase*>( dimredbase );
      43             :   }
      44             : 
      45             :   // Read in the switching functions
      46             :   std::string linput,hinput, errors;
      47          12 :   parse("HIGH_DIM_FUNCTION",hinput);
      48           6 :   if( hinput=="as in input action" ) {
      49           1 :     if( !smapbase ) {
      50           0 :       error("high dimensional switching function has not been set - use HIGH_DIM_FUNCTION");
      51             :     }
      52           1 :     reuse_hd=true;
      53           1 :     log.printf("  reusing high dimensional filter function defined in previous sketch-map action\n");
      54             :   } else {
      55           5 :     reuse_hd=false;
      56           5 :     highdf.set(hinput,errors);
      57           5 :     if(errors.length()>0) {
      58           0 :       error(errors);
      59             :     }
      60          10 :     log.printf("  filter function for dissimilarities in high dimensional space has cutoff %s \n",highdf.description().c_str() );
      61             :   }
      62             : 
      63          12 :   parse("LOW_DIM_FUNCTION",linput);
      64           6 :   if( linput=="as in input action" ) {
      65           1 :     if( !smapbase ) {
      66           0 :       error("low dimensional switching function has not been set - use LOW_DIM_FUNCTION");
      67             :     }
      68           1 :     reuse_ld=true;
      69           1 :     log.printf("  reusing low dimensional filter function defined in previous sketch-map action\n");
      70             :   } else {
      71           5 :     reuse_ld=false;
      72           5 :     lowdf.set(linput,errors);
      73           5 :     if(errors.length()>0) {
      74           0 :       error(errors);
      75             :     }
      76          10 :     log.printf("  filter function for distances in low dimensionality space has cutoff %s \n",lowdf.description().c_str() );
      77             :   }
      78             : 
      79             :   // Read the mixing parameter
      80           6 :   parse("MIXPARAM",mixparam);
      81           6 :   if( mixparam<0 || mixparam>1 ) {
      82           0 :     error("mixing parameter must be between 0 and 1");
      83             :   }
      84           6 :   log.printf("  mixing %f of pure distances with %f of filtered distances \n",mixparam,1.-mixparam);
      85           6 : }
      86             : 
      87           7 : void SketchMapBase::calculateProjections( const Matrix<double>& targets, Matrix<double>& projections ) {
      88           7 :   if( dtargets.size()!=targets.nrows() ) {
      89             :     // These hold data so that we can do stress calculations
      90           6 :     dtargets.resize( targets.nrows() );
      91           6 :     ftargets.resize( targets.nrows() );
      92           6 :     pweights.resize( targets.nrows() );
      93             :     // Matrices for storing input data
      94             :     transformed.resize( targets.nrows(), targets.ncols() );
      95             :     distances.resize( targets.nrows(), targets.ncols() );
      96             :   }
      97             : 
      98             :   // Stores the weights in an array for faster access, as well as the normalization
      99           7 :   normw=0;
     100        1141 :   for(unsigned i=0; i<targets.nrows() ; ++i) {
     101        1134 :     pweights[i] = getWeight(i);
     102        1134 :     normw+=pweights[i];
     103             :   }
     104           7 :   normw*=normw;
     105             : 
     106             :   // Transform the high dimensional distances
     107             :   double df;
     108             :   distances=0.;
     109             :   transformed=0.;
     110        1127 :   for(unsigned i=1; i<distances.ncols(); ++i) {
     111      172838 :     for(unsigned j=0; j<i; ++j) {
     112      171711 :       distances(i,j)=distances(j,i)=std::sqrt( targets(i,j) );
     113      171711 :       transformed(i,j)=transformed(j,i)=transformHighDimensionalDistance( distances(i,j), df );
     114             :     }
     115             :   }
     116             :   // And minimse
     117           7 :   minimise( projections );
     118           7 : }
     119             : 
     120       41344 : double SketchMapBase::calculateStress( const std::vector<double>& p, std::vector<double>& d ) {
     121             :   // Zero derivative and stress accumulators
     122      124032 :   for(unsigned i=0; i<p.size(); ++i) {
     123       82688 :     d[i]=0.0;
     124             :   }
     125             :   double stress=0;
     126       41344 :   std::vector<double> dtmp( p.size() );
     127             :   // Now accumulate total stress on system
     128     6687344 :   for(unsigned i=0; i<ftargets.size(); ++i) {
     129     6646000 :     if( dtargets[i]<epsilon ) {
     130       30215 :       continue ;
     131             :     }
     132             : 
     133             :     // Calculate distance in low dimensional space
     134     6615785 :     double dd=0;
     135    19847355 :     for(unsigned j=0; j<p.size(); ++j) {
     136    13231570 :       dtmp[j]=p[j]-projections(i,j);
     137    13231570 :       dd+=dtmp[j]*dtmp[j];
     138             :     }
     139     6615785 :     dd = std::sqrt(dd);
     140             : 
     141             :     // Now do transformations and calculate differences
     142     6615785 :     double df, fd = transformLowDimensionalDistance( dd, df );
     143     6615785 :     double ddiff = dd - dtargets[i];
     144     6615785 :     double fdiff = fd - ftargets[i];
     145             : 
     146             :     // Calculate derivatives
     147     6615785 :     double pref = 2.*getWeight(i) / dd ;
     148    19847355 :     for(unsigned j=0; j<p.size(); ++j) {
     149    13231570 :       d[j] += pref*( (1-mixparam)*fdiff*df + mixparam*ddiff )*dtmp[j];
     150             :     }
     151             : 
     152             :     // Accumulate the total stress
     153     6615785 :     stress += getWeight(i)*( (1-mixparam)*fdiff*fdiff + mixparam*ddiff*ddiff );
     154             :   }
     155       41344 :   return stress;
     156             : }
     157             : 
     158        1913 : double SketchMapBase::calculateFullStress( const std::vector<double>& p, std::vector<double>& d ) {
     159             :   // Zero derivative and stress accumulators
     160      393213 :   for(unsigned i=0; i<p.size(); ++i) {
     161      391300 :     d[i]=0.0;
     162             :   }
     163             :   double stress=0;
     164        1913 :   std::vector<double> dtmp( nlow );
     165             : 
     166      195650 :   for(unsigned i=1; i<distances.nrows(); ++i) {
     167      193737 :     double iweight = pweights[i];
     168    14802162 :     for(unsigned j=0; j<i; ++j) {
     169    14608425 :       double jweight =  pweights[j];
     170             :       // Calculate distance in low dimensional space
     171    14608425 :       double dd=0;
     172    43825275 :       for(unsigned k=0; k<nlow; ++k) {
     173    29216850 :         dtmp[k]=p[nlow*i+k] - p[nlow*j+k];
     174    29216850 :         dd+=dtmp[k]*dtmp[k];
     175             :       }
     176    14608425 :       dd = std::sqrt(dd);
     177             : 
     178             :       // Now do transformations and calculate differences
     179    14608425 :       double df, fd = transformLowDimensionalDistance( dd, df );
     180    14608425 :       double ddiff = dd - distances(i,j);
     181    14608425 :       double fdiff = fd - transformed(i,j);;
     182             : 
     183             :       // Calculate derivatives
     184    14608425 :       double pref = 2.*iweight*jweight*( (1-mixparam)*fdiff*df + mixparam*ddiff ) / dd;
     185    43825275 :       for(unsigned k=0; k<nlow; ++k) {
     186    29216850 :         double dterm=pref*dtmp[k];
     187    29216850 :         d[nlow*i+k]+=dterm;
     188    29216850 :         d[nlow*j+k]-=dterm;
     189             :       }
     190             : 
     191             :       // Accumulate the total stress
     192    14608425 :       stress += iweight*jweight*( (1-mixparam)*fdiff*fdiff + mixparam*ddiff*ddiff );
     193             :     }
     194             :   }
     195        1913 :   stress /= normw;
     196      393213 :   for (unsigned k=0; k < d.size(); ++k) {
     197      391300 :     d[k] /= normw;
     198             :   }
     199        1913 :   return stress;
     200             : }
     201             : 
     202             : }
     203             : }
     204             : 
     205             :