Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2016-2020 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 "core/ActionWithValue.h"
23 : #include "core/ActionWithArguments.h"
24 : #include "core/ActionRegister.h"
25 : #include "tools/ConjugateGradient.h"
26 : #include "tools/SwitchingFunction.h"
27 : #include "gridtools/GridSearch.h"
28 : #include "SMACOF.h"
29 :
30 : namespace PLMD {
31 : namespace dimred {
32 :
33 : //+PLUMEDOC DIMRED ARRANGE_POINTS
34 : /*
35 : Arrange points in a low dimensional space so that the (transformed) distances between points in the low dimensional space match the dissimilarities provided in an input matrix.
36 :
37 : \par Examples
38 :
39 : */
40 : //+ENDPLUMEDOC
41 :
42 : class ArrangePoints :
43 : public ActionWithValue,
44 : public ActionWithArguments {
45 : private:
46 : unsigned dimout, maxiter, ncycles, current_index;
47 : double cgtol, gbuf;
48 : std::vector<unsigned> npoints, nfgrid;
49 : std::vector<double> mypos;
50 : double smacof_tol, smacof_reg;
51 : int dist_target;
52 : enum {conjgrad,pointwise,smacof} mintype;
53 : std::vector<SwitchingFunction> switchingFunction;
54 : void checkInputMatrix( const std::string& key, const unsigned& nvals, const std::vector<Value*>& mat ) const ;
55 : double recalculateSmacofWeights( const std::vector<double>& p, SMACOF& mysmacof ) const ;
56 : protected:
57 : double calculateStress( const std::vector<double>& p, std::vector<double>& d );
58 : double calculateFullStress( const std::vector<double>& p, std::vector<double>& d );
59 : public:
60 : static void registerKeywords( Keywords& keys );
61 : ArrangePoints( const ActionOptions& );
62 0 : unsigned getNumberOfDerivatives() override { return 0; }
63 : void prepare() override ;
64 : void calculate() override ;
65 : virtual void optimize( std::vector<double>& pos );
66 2 : void apply() override {}
67 : };
68 :
69 : PLUMED_REGISTER_ACTION(ArrangePoints,"ARRANGE_POINTS")
70 :
71 12 : void ArrangePoints::registerKeywords( Keywords& keys ) {
72 12 : Action::registerKeywords( keys ); ActionWithValue::registerKeywords( keys );
73 12 : ActionWithArguments::registerKeywords( keys ); keys.use("ARG");
74 24 : keys.add("numbered","TARGET","the matrix of target quantities that you would like to match");
75 24 : keys.add("numbered","FUNC","a function that is applied on the distances between the points in the low dimensional space");
76 24 : keys.add("numbered","WEIGHTS","the matrix with the weights of the target quantities");
77 24 : keys.add("compulsory","MINTYPE","conjgrad","the method to use for the minimisation");
78 24 : keys.add("compulsory","MAXITER","1000","maximum number of optimization cycles for optimisation algorithms");
79 24 : keys.add("compulsory","CGTOL","1E-6","the tolerance for the conjugate gradient minimization");
80 24 : keys.add("compulsory","NCYCLES","5","the number of cycles of global optimization to attempt");
81 24 : keys.add("compulsory","BUFFER","1.1","grid extent for search is (max projection - minimum projection) multiplied by this value");
82 24 : keys.add("compulsory","CGRID_SIZE","10","number of points to use in each grid direction");
83 24 : keys.add("compulsory","FGRID_SIZE","0","interpolate the grid onto this number of points -- only works in 2D");
84 24 : keys.add("compulsory","SMACTOL","1E-4","the tolerance for the smacof algorithm");
85 24 : keys.add("compulsory","SMACREG","0.001","this is used to ensure that we don't divide by zero when updating weights for SMACOF algorithm");
86 24 : keys.addOutputComponent("coord","default","the coordinates of the points in the low dimensional space");
87 12 : }
88 :
89 :
90 5 : ArrangePoints::ArrangePoints( const ActionOptions& ao ) :
91 : Action(ao),
92 : ActionWithValue(ao),
93 : ActionWithArguments(ao),
94 5 : current_index(0),
95 5 : dist_target(-1)
96 : {
97 5 : dimout = getNumberOfArguments();
98 5 : std::vector<unsigned> shape(1); shape[0]=getPntrToArgument(0)->getNumberOfValues();
99 15 : for(unsigned i=0; i<getNumberOfArguments(); ++i) {
100 10 : if( shape[0]!=getPntrToArgument(i)->getNumberOfValues() ) error("mismatch between sizes of input coordinates");
101 10 : std::string num; Tools::convert( i+1, num ); addComponent( "coord-" + num, shape );
102 20 : componentIsNotPeriodic( "coord-" + num ); getPntrToArgument(i)->buildDataStore();
103 : }
104 5 : std::vector<Value*> args( getArguments() ), target, weights; std::string sfd, errors;
105 : // Read in target "distances" and target weights
106 5 : for(unsigned i=1;; ++i) {
107 22 : target.resize(0); if( !parseArgumentList("TARGET",i,target) ) break;
108 6 : std::string inum; Tools::convert( i, inum ); checkInputMatrix( "TARGET" + inum, shape[0], target );
109 12 : if( !parseArgumentList("WEIGHTS",i,weights) ) error("missing WEIGHTS" + inum + " keyword in input");
110 12 : checkInputMatrix( "WEIGHTS" + inum, shape[0], weights );
111 6 : args.push_back( target[0] ); args.push_back( weights[0] );
112 12 : bool has_sf = parseNumbered("FUNC",i,sfd); switchingFunction.push_back( SwitchingFunction() );
113 6 : if( !has_sf ) {
114 2 : switchingFunction[i-1].set( "CUSTOM FUNC=1-sqrt(x2) R_0=1.0", errors ); dist_target=i-1;
115 : } else {
116 5 : switchingFunction[i-1].set( sfd, errors );
117 5 : if( errors.length()!=0 ) error("problem reading switching function description " + errors);
118 : }
119 6 : log.printf(" %sth term seeks to match tranformed distances with those in matrix %s \n", inum.c_str(), target[0]->getName().c_str() );
120 6 : log.printf(" in %sth term distances are transformed by 1-switching function with r_0=%s \n", inum.c_str(), switchingFunction[i-1].description().c_str() );
121 6 : log.printf(" in %sth term weights of matrix elements in stress function are given by %s \n", inum.c_str(), weights[0]->getName().c_str() );
122 6 : }
123 10 : std::string mtype; parse("MINTYPE",mtype);
124 5 : if( mtype=="conjgrad" ) {
125 3 : mintype=conjgrad;
126 3 : log.printf(" minimimising stress function using conjugate gradients\n");
127 2 : } else if( mtype=="pointwise") {
128 1 : mintype=pointwise;
129 1 : log.printf(" minimimising stress function using pointwise global optimisation\n");
130 1 : npoints.resize(dimout); nfgrid.resize(dimout);
131 3 : parseVector("CGRID_SIZE",npoints); parse("BUFFER",gbuf); parse("NCYCLES",ncycles);
132 2 : parseVector("FGRID_SIZE",nfgrid); if( nfgrid[0]!=0 && dimout!=2 ) error("interpolation only works in two dimensions");
133 1 : log.printf(" doing %u cycles of global optimization sweeps\n",ncycles);
134 1 : log.printf(" using coarse grid of points that is %u",npoints[0]);
135 2 : for(unsigned j=1; j<npoints.size(); ++j) log.printf(" by %u",npoints[j]);
136 1 : log.printf("\n grid is %f times larger than the difference between the position of the minimum and maximum projection \n",gbuf);
137 1 : if( nfgrid[0]>0 ) {
138 1 : log.printf(" interpolating stress onto grid of points that is %u",nfgrid[0]);
139 2 : for(unsigned j=1; j<nfgrid.size(); ++j) log.printf(" by %u",nfgrid[j]);
140 1 : log.printf("\n");
141 : }
142 1 : } else if( mtype=="smacof" ) {
143 1 : mintype=smacof; if( dist_target<0 ) error("one of targets must be distances in order to use smacof");
144 1 : log.printf(" minimising stress fucntion using smacof\n");
145 2 : parse("SMACTOL",smacof_tol); parse("SMACREG",smacof_reg);
146 1 : log.printf(" tolerance for smacof algorithms equals %f \n", smacof_tol);
147 1 : log.printf(" using %f as regularisation parameter for weights in smacof algorithm\n", smacof_reg);
148 0 : } else error("invalid MINTYPE");
149 5 : if( mintype!=smacof) {
150 8 : parse("CGTOL",cgtol); log.printf(" tolerance for conjugate gradient algorithm equals %f \n",cgtol);
151 : }
152 10 : parse("MAXITER",maxiter); log.printf(" maximum number of iterations for minimimization algorithms equals %d \n", maxiter );
153 5 : requestArguments( args ); checkRead();
154 5 : }
155 :
156 12 : void ArrangePoints::checkInputMatrix( const std::string& key, const unsigned& nvals, const std::vector<Value*>& mat ) const {
157 12 : mat[0]->buildDataStore();
158 12 : if( mat.size()!=1 ) error("should only be one value in input to " + key );
159 12 : if( mat[0]->getRank()!=2 || mat[0]->hasDerivatives() ) error("input to " + key + " keyword should be a matrix");
160 12 : if( mat[0]->getShape()[0]!=nvals || mat[0]->getShape()[1]!=nvals ) error("input to " + key + " keyword has the wrong size");
161 12 : }
162 :
163 24600 : double ArrangePoints::calculateStress( const std::vector<double>& p, std::vector<double>& d ) {
164 73800 : double stress=0; for(unsigned i=0; i<p.size(); ++i) d[i]=0.0;
165 24600 : std::vector<double> dtmp(dimout);
166 24600 : std::vector<unsigned> shape( getPntrToArgument( dimout )->getShape() );
167 24600 : unsigned targi=shape[0]*current_index;
168 24600 : unsigned nmatrices = ( getNumberOfArguments() - dimout ) / 2;
169 2484600 : for(unsigned i=0; i<shape[0]; ++i) {
170 2460000 : if( i==current_index ) continue ;
171 : // Calculate distance in low dimensional space
172 7306200 : double dd2=0; for(unsigned k=0; k<dimout; ++k) { dtmp[k]=p[k] - mypos[dimout*i+k]; dd2+=dtmp[k]*dtmp[k]; }
173 :
174 4870800 : for(unsigned k=0; k<nmatrices; ++k ) {
175 : // Now do transformations and calculate differences
176 2435400 : double df, fd = 1. - switchingFunction[k].calculateSqr( dd2, df );
177 : // Get the weight for this connection
178 : double weight = 0;
179 245975400 : for(unsigned j=0; j<shape[0]; ++j) weight += getPntrToArgument( dimout + 2*k + 1 )->get( shape[0]*i+j );
180 : // Get the difference for the connection
181 2435400 : double fdiff = fd - getPntrToArgument( dimout + 2*k )->get( targi+i );
182 : // Calculate derivatives
183 2435400 : double pref = -2.*weight*fdiff*df;
184 7306200 : for(unsigned n=0; n<dimout; ++n) d[n] += pref*dtmp[n];
185 : // Accumulate the total stress
186 2435400 : stress += weight*fdiff*fdiff;
187 : }
188 : }
189 24600 : return stress;
190 : }
191 :
192 1932 : double ArrangePoints::calculateFullStress( const std::vector<double>& p, std::vector<double>& d ) {
193 : // Zero derivative and stress accumulators
194 387932 : for(unsigned i=0; i<p.size(); ++i) d[i]=0.0;
195 1932 : double stress=0; std::vector<double> dtmp( dimout );
196 :
197 1932 : unsigned nmatrices = ( getNumberOfArguments() - dimout ) / 2;
198 1932 : std::vector<unsigned> shape( getPntrToArgument( dimout )->getShape() );
199 193000 : for(unsigned i=1; i<shape[0]; ++i) {
200 14134568 : for(unsigned j=0; j<i; ++j) {
201 : // Calculate distance in low dimensional space
202 41830500 : double dd2=0; for(unsigned k=0; k<dimout; ++k) { dtmp[k]=p[dimout*i+k] - p[dimout*j+k]; dd2+=dtmp[k]*dtmp[k]; }
203 :
204 27887000 : for(unsigned k=0; k<nmatrices; ++k ) {
205 : // Now do transformations and calculate differences
206 13943500 : double df, fd = 1. - switchingFunction[k].calculateSqr( dd2, df );
207 : // Get the weight for this connection
208 13943500 : double weight = getPntrToArgument( dimout + 2*k + 1 )->get( shape[0]*i+j );
209 : // Get the difference for the connection
210 13943500 : double fdiff = fd - getPntrToArgument( dimout + 2*k )->get( shape[0]*i+j );
211 : // Calculate derivatives
212 13943500 : double pref = -2.*weight*fdiff*df;
213 41830500 : for(unsigned n=0; n<dimout; ++n) { double dterm=pref*dtmp[n]; d[dimout*i+n]+=dterm; d[dimout*j+n]-=dterm; }
214 : // Accumulate the total stress
215 13943500 : stress += weight*fdiff*fdiff;
216 : }
217 : }
218 : }
219 1932 : return stress;
220 : }
221 :
222 2 : double ArrangePoints::recalculateSmacofWeights( const std::vector<double>& p, SMACOF& mysmacof ) const {
223 : double stress=0, totalWeight=0;
224 2 : unsigned nmatrices = ( getNumberOfArguments() - dimout ) / 2;
225 2 : std::vector<unsigned> shape( getPntrToArgument( dimout )->getShape() );
226 1000 : for(unsigned i=1; i<shape[0]; ++i) {
227 250498 : for(unsigned j=0; j<i; ++j) {
228 : // Calculate distance in low dimensional space
229 748500 : double dd2=0; for(unsigned k=0; k<dimout; ++k) { double dtmp=p[dimout*i+k] - p[dimout*j+k]; dd2+=dtmp*dtmp; }
230 : // Calculate difference between target difference and true difference
231 249500 : double wval=0, dd1 = sqrt(dd2); double diff = mysmacof.getDistance(i,j) - dd1;
232 :
233 748500 : for(unsigned k=0; k<nmatrices; ++k ) {
234 : // Don't need to do anything for distances we are matching
235 499000 : if( k==dist_target ) continue;
236 : // Now do transformations and calculate differences
237 249500 : double df, fd = 1. - switchingFunction[k].calculateSqr( dd2, df );
238 : // Get the weight for this connection
239 249500 : double weight = getPntrToArgument( dimout + 2*k + 1 )->get( shape[0]*i+j );
240 : // Get the difference for the connection
241 249500 : double fdiff = getPntrToArgument( dimout + 2*k )->get( shape[0]*i+j ) - fd;
242 : // Now set the weight if difference in distance is larger than regularisation parameter
243 249500 : if( fabs(diff)>smacof_reg ) wval -= weight*fdiff*df*dd1 / diff;
244 : // And the total stress and weights
245 249500 : stress += weight*fdiff*fdiff; totalWeight += weight;
246 : }
247 249500 : mysmacof.setWeight( j, i, wval ); mysmacof.setWeight( i, j, wval );
248 : }
249 : }
250 2 : return stress / totalWeight;
251 : }
252 :
253 7 : void ArrangePoints::optimize( std::vector<double>& pos ) {
254 : ConjugateGradient<ArrangePoints> mycgminimise( this );
255 7 : if( mintype==conjgrad ) {
256 5 : mycgminimise.minimise( cgtol, pos, &ArrangePoints::calculateFullStress );
257 2 : } else if( mintype==pointwise ) {
258 1 : unsigned nvals=getPntrToArgument( dimout )->getShape()[0];
259 1 : std::vector<double> gmin( dimout ), gmax( dimout ), mypoint( dimout );
260 : // Find the extent of the grid
261 3 : for(unsigned j=0; j<dimout; ++j) gmin[j]=gmax[j]=pos[j];
262 100 : for(unsigned i=1; i<nvals; ++i) {
263 297 : for(unsigned j=0; j<dimout; ++j) {
264 198 : if( pos[dimout*i+j] < gmin[j] ) gmin[j] = pos[dimout*i+j];
265 198 : if( pos[dimout*i+j] > gmax[j] ) gmax[j] = pos[dimout*i+j];
266 : }
267 : }
268 3 : for(unsigned j=0; j<dimout; ++j) {
269 2 : double gbuffer = 0.5*gbuf*( gmax[j]-gmin[j] ) - 0.5*( gmax[j]- gmin[j] );
270 2 : gmin[j]-=gbuffer; gmax[j]+=gbuffer;
271 : }
272 202 : mypos.resize( pos.size() ); for(unsigned i=0; i<mypos.size(); ++i) mypos[i] = pos[i];
273 1 : gridtools::GridSearch<ArrangePoints> mygridsearch( gmin, gmax, npoints, nfgrid, this );
274 : // Run multiple loops over all projections
275 3 : for(unsigned i=0; i<ncycles; ++i) {
276 202 : for(unsigned j=0; j<nvals; ++j) {
277 : // Setup target distances and target functions for calculate stress
278 200 : current_index=j;
279 :
280 : // Find current projection of jth point
281 600 : for(unsigned k=0; k<dimout; ++k) mypoint[k]=mypos[j*dimout+k];
282 : // Minimise using grid search
283 200 : bool moved=mygridsearch.minimise( mypoint, &ArrangePoints::calculateStress );
284 200 : if( moved ) {
285 : // Reassign the new projection
286 66 : for(unsigned k=0; k<dimout; ++k) mypos[dimout*j+k]=mypoint[k];
287 : // Minimise output using conjugate gradient
288 22 : mycgminimise.minimise( cgtol, mypos, &ArrangePoints::calculateFullStress );
289 : }
290 : }
291 402 : for(unsigned i=0; i<mypos.size(); ++i) pos[i] = mypos[i];
292 : }
293 2 : } else if( mintype==smacof ) {
294 1 : SMACOF mysmacof( getPntrToArgument( dimout + 2*dist_target) ); double stress = recalculateSmacofWeights( pos, mysmacof );
295 :
296 1 : for(unsigned i=0; i<maxiter; ++i) {
297 : // Optimise using smacof and current weights
298 1 : mysmacof.optimize( smacof_tol, maxiter, pos );
299 : // Recalculate weights matrix and sigma
300 1 : double newsig = recalculateSmacofWeights( pos, mysmacof );
301 : // Test whether or not the algorithm has converged
302 1 : if( fabs( newsig - stress )<smacof_tol ) break;
303 : // Make initial sigma into new sigma so that the value of new sigma is used every time so that the error can be reduced
304 : stress=newsig;
305 : }
306 : }
307 7 : }
308 :
309 8 : void ArrangePoints::prepare() {
310 : // Make sure all the components are the right size
311 8 : std::vector<unsigned> shape(1,getPntrToArgument( dimout )->getShape()[0]);
312 24 : for(unsigned j=0; j<dimout; ++j) {
313 16 : if( getPntrToComponent(j)->getShape()[0]!=shape[0] ) getPntrToComponent(j)->setShape( shape );
314 : }
315 8 : }
316 :
317 7 : void ArrangePoints::calculate() {
318 : // Retrive the initial value
319 7 : unsigned nvals = getPntrToArgument( dimout )->getShape()[0];
320 7 : std::vector<double> pos( dimout*nvals );
321 1107 : for(unsigned i=0; i<nvals; ++i) {
322 3300 : for(unsigned j=0; j<dimout; ++j) pos[ dimout*i + j ] = getPntrToArgument(j)->get(i);
323 : }
324 : // Do the optimization
325 7 : optimize( pos );
326 : // And set the final values
327 1107 : for(unsigned i=0; i<nvals; ++i) {
328 3300 : for(unsigned j=0; j<dimout; ++j) getPntrToComponent(j)->set( i, pos[dimout*i+j] );
329 : }
330 7 : }
331 :
332 : }
333 : }
|
