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