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          10 : void SketchMapBase::registerKeywords( Keywords& keys ) {
      28          10 :   DimensionalityReductionBase::registerKeywords( keys );
      29          10 :   keys.remove("NLOW_DIM");
      30          40 :   keys.add("compulsory","HIGH_DIM_FUNCTION","as in input action","the parameters of the switching function in the high dimensional space");
      31          40 :   keys.add("compulsory","LOW_DIM_FUNCTION","as in input action","the parameters of the switching function in the low dimensional space");
      32          20 :   keys.add("compulsory","MIXPARAM","0.0","the amount of the pure distances to mix into the stress function");
      33          10 : }
      34             : 
      35           6 : SketchMapBase::SketchMapBase( const ActionOptions& ao ):
      36             :   Action(ao),
      37             :   DimensionalityReductionBase(ao),
      38           6 :   smapbase(NULL),
      39           6 :   normw(0.0)
      40             : {
      41             :   // Check if we have data from a input sketch-map object - we can reuse switching functions wahoo!!
      42           6 :   if( dimredbase ) smapbase = dynamic_cast<SketchMapBase*>( dimredbase );
      43             : 
      44             :   // Read in the switching functions
      45             :   std::string linput,hinput, errors;
      46          12 :   parse("HIGH_DIM_FUNCTION",hinput);
      47           6 :   if( hinput=="as in input action" ) {
      48           1 :     if( !smapbase ) error("high dimensional switching function has not been set - use HIGH_DIM_FUNCTION");
      49           1 :     reuse_hd=true;
      50           1 :     log.printf("  reusing high dimensional filter function defined in previous sketch-map action\n");
      51             :   } else {
      52           5 :     reuse_hd=false;
      53           5 :     highdf.set(hinput,errors);
      54           5 :     if(errors.length()>0) error(errors);
      55          10 :     log.printf("  filter function for dissimilarities in high dimensional space has cutoff %s \n",highdf.description().c_str() );
      56             :   }
      57             : 
      58          12 :   parse("LOW_DIM_FUNCTION",linput);
      59           6 :   if( linput=="as in input action" ) {
      60           1 :     if( !smapbase ) error("low dimensional switching function has not been set - use LOW_DIM_FUNCTION");
      61           1 :     reuse_ld=true;
      62           1 :     log.printf("  reusing low dimensional filter function defined in previous sketch-map action\n");
      63             :   } else {
      64           5 :     reuse_ld=false;
      65           5 :     lowdf.set(linput,errors);
      66           5 :     if(errors.length()>0) error(errors);
      67          10 :     log.printf("  filter function for distances in low dimensionality space has cutoff %s \n",lowdf.description().c_str() );
      68             :   }
      69             : 
      70             :   // Read the mixing parameter
      71           6 :   parse("MIXPARAM",mixparam);
      72           6 :   if( mixparam<0 || mixparam>1 ) error("mixing parameter must be between 0 and 1");
      73           6 :   log.printf("  mixing %f of pure distances with %f of filtered distances \n",mixparam,1.-mixparam);
      74           6 : }
      75             : 
      76           7 : void SketchMapBase::calculateProjections( const Matrix<double>& targets, Matrix<double>& projections ) {
      77           7 :   if( dtargets.size()!=targets.nrows() ) {
      78             :     // These hold data so that we can do stress calculations
      79           6 :     dtargets.resize( targets.nrows() ); ftargets.resize( targets.nrows() ); pweights.resize( targets.nrows() );
      80             :     // Matrices for storing input data
      81             :     transformed.resize( targets.nrows(), targets.ncols() );
      82             :     distances.resize( targets.nrows(), targets.ncols() );
      83             :   }
      84             : 
      85             :   // Stores the weights in an array for faster access, as well as the normalization
      86           7 :   normw=0;
      87        1141 :   for(unsigned i=0; i<targets.nrows() ; ++i) { pweights[i] = getWeight(i); normw+=pweights[i]; }
      88           7 :   normw*=normw;
      89             : 
      90             :   // Transform the high dimensional distances
      91             :   double df; distances=0.; transformed=0.;
      92        1134 :   for(unsigned i=1; i<distances.ncols(); ++i) {
      93      172838 :     for(unsigned j=0; j<i; ++j) {
      94      171711 :       distances(i,j)=distances(j,i)=std::sqrt( targets(i,j) );
      95      171711 :       transformed(i,j)=transformed(j,i)=transformHighDimensionalDistance( distances(i,j), df );
      96             :     }
      97             :   }
      98             :   // And minimse
      99           7 :   minimise( projections );
     100           7 : }
     101             : 
     102       41344 : double SketchMapBase::calculateStress( const std::vector<double>& p, std::vector<double>& d ) {
     103             :   // Zero derivative and stress accumulators
     104      124032 :   for(unsigned i=0; i<p.size(); ++i) d[i]=0.0;
     105       41344 :   double stress=0; std::vector<double> dtmp( p.size() );
     106             :   // Now accumulate total stress on system
     107     6687344 :   for(unsigned i=0; i<ftargets.size(); ++i) {
     108     6646000 :     if( dtargets[i]<epsilon ) continue ;
     109             : 
     110             :     // Calculate distance in low dimensional space
     111     6615785 :     double dd=0;
     112    19847355 :     for(unsigned j=0; j<p.size(); ++j) { dtmp[j]=p[j]-projections(i,j); dd+=dtmp[j]*dtmp[j]; }
     113     6615785 :     dd = std::sqrt(dd);
     114             : 
     115             :     // Now do transformations and calculate differences
     116     6615785 :     double df, fd = transformLowDimensionalDistance( dd, df );
     117     6615785 :     double ddiff = dd - dtargets[i];
     118     6615785 :     double fdiff = fd - ftargets[i];
     119             : 
     120             :     // Calculate derivatives
     121     6615785 :     double pref = 2.*getWeight(i) / dd ;
     122    19847355 :     for(unsigned j=0; j<p.size(); ++j) d[j] += pref*( (1-mixparam)*fdiff*df + mixparam*ddiff )*dtmp[j];
     123             : 
     124             :     // Accumulate the total stress
     125     6615785 :     stress += getWeight(i)*( (1-mixparam)*fdiff*fdiff + mixparam*ddiff*ddiff );
     126             :   }
     127       41344 :   return stress;
     128             : }
     129             : 
     130        1913 : double SketchMapBase::calculateFullStress( const std::vector<double>& p, std::vector<double>& d ) {
     131             :   // Zero derivative and stress accumulators
     132      393213 :   for(unsigned i=0; i<p.size(); ++i) d[i]=0.0;
     133        1913 :   double stress=0; std::vector<double> dtmp( nlow );
     134             : 
     135      195650 :   for(unsigned i=1; i<distances.nrows(); ++i) {
     136      193737 :     double iweight = pweights[i];
     137    14802162 :     for(unsigned j=0; j<i; ++j) {
     138    14608425 :       double jweight =  pweights[j];
     139             :       // Calculate distance in low dimensional space
     140    14608425 :       double dd=0;
     141    43825275 :       for(unsigned k=0; k<nlow; ++k) { dtmp[k]=p[nlow*i+k] - p[nlow*j+k]; dd+=dtmp[k]*dtmp[k]; }
     142    14608425 :       dd = std::sqrt(dd);
     143             : 
     144             :       // Now do transformations and calculate differences
     145    14608425 :       double df, fd = transformLowDimensionalDistance( dd, df );
     146    14608425 :       double ddiff = dd - distances(i,j);
     147    14608425 :       double fdiff = fd - transformed(i,j);;
     148             : 
     149             :       // Calculate derivatives
     150    14608425 :       double pref = 2.*iweight*jweight*( (1-mixparam)*fdiff*df + mixparam*ddiff ) / dd;
     151    43825275 :       for(unsigned k=0; k<nlow; ++k) {
     152    29216850 :         double dterm=pref*dtmp[k]; d[nlow*i+k]+=dterm; d[nlow*j+k]-=dterm;
     153             :       }
     154             : 
     155             :       // Accumulate the total stress
     156    14608425 :       stress += iweight*jweight*( (1-mixparam)*fdiff*fdiff + mixparam*ddiff*ddiff );
     157             :     }
     158             :   }
     159      393213 :   stress /= normw; for (unsigned k=0; k < d.size(); ++k) d[k] /= normw;
     160        1913 :   return stress;
     161             : }
     162             : 
     163             : }
     164             : }
     165             : 
     166             :