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/ActionWithVector.h"
23 : #include "core/ActionRegister.h"
24 : #include "tools/ConjugateGradient.h"
25 : #include "tools/SwitchingFunction.h"
26 : #include "tools/OpenMP.h"
27 : #include "tools/Random.h"
28 :
29 : namespace PLMD {
30 : namespace dimred {
31 :
32 : //+PLUMEDOC DIMRED PROJECT_POINTS
33 : /*
34 : Find the projection of a point in a low dimensional space by matching the (transformed) distance between it and a series of reference configurations that were input
35 :
36 : \par Examples
37 :
38 : */
39 : //+ENDPLUMEDOC
40 :
41 : class ProjectPoints : public ActionWithVector {
42 : private:
43 : double cgtol;
44 : unsigned dimout;
45 : mutable std::vector<unsigned> rowstart;
46 : std::vector<SwitchingFunction> switchingFunction;
47 : ConjugateGradient<ProjectPoints> myminimiser;
48 : void getProjection( const unsigned& current, std::vector<double>& point ) const ;
49 : public:
50 : static void registerKeywords( Keywords& keys );
51 : ProjectPoints( const ActionOptions& );
52 0 : unsigned getNumberOfDerivatives() {
53 0 : return 0;
54 : }
55 : void prepare() override ;
56 : void performTask( const unsigned& current, MultiValue& myvals ) const override ;
57 : double calculateStress( const std::vector<double>& pp, std::vector<double>& der );
58 : void calculate() override ;
59 168 : void apply() override {}
60 : };
61 :
62 : PLUMED_REGISTER_ACTION(ProjectPoints,"PROJECT_POINTS")
63 :
64 6 : void ProjectPoints::registerKeywords( Keywords& keys ) {
65 6 : ActionWithVector::registerKeywords( keys );
66 12 : keys.addInputKeyword("compulsory","ARG","vector","the projections of the landmark points");
67 12 : keys.addInputKeyword("numbered","TARGET","vector/matrix","the matrix of target quantities that you would like to match");
68 6 : keys.add("numbered","FUNC","a function that is applied on the distances between the points in the low dimensional space");
69 12 : keys.addInputKeyword("numbered","WEIGHTS","vector","the matrix with the weights of the target quantities");
70 6 : keys.add("compulsory","CGTOL","1E-6","the tolerance for the conjugate gradient minimization");
71 12 : keys.addOutputComponent("coord","default","scalar/vector","the coordinates of the points in the low dimensional space");
72 6 : }
73 :
74 :
75 2 : ProjectPoints::ProjectPoints( const ActionOptions& ao ) :
76 : Action(ao),
77 : ActionWithVector(ao),
78 2 : rowstart(OpenMP::getNumThreads()),
79 4 : myminimiser( this ) {
80 2 : dimout = getNumberOfArguments();
81 2 : unsigned nvals=getPntrToArgument(0)->getNumberOfValues();
82 6 : for(unsigned i=0; i<getNumberOfArguments(); ++i) {
83 4 : if( nvals!=getPntrToArgument(i)->getNumberOfValues() ) {
84 0 : error("mismatch between numbers of projections");
85 : }
86 : }
87 2 : std::vector<Value*> args( getArguments() ), target, weights;
88 : std::string sfd, errors;
89 : unsigned ntoproj=0;
90 : // Read in target "distances" and target weights
91 2 : for(unsigned i=1;; ++i) {
92 4 : target.resize(0);
93 8 : if( !parseArgumentList("TARGET",i,target) ) {
94 : break;
95 : }
96 : std::string inum;
97 2 : Tools::convert( i, inum );
98 2 : if( target.size()!=1 ) {
99 0 : error("should only be one value in input to TARGET" + inum );
100 : }
101 2 : if( (target[0]->getRank()!=1 && target[0]->getRank()!=2) || target[0]->hasDerivatives() ) {
102 0 : error("input to TARGET" + inum + " keyword should be a vector/matrix");
103 : }
104 2 : if( target[0]->getShape()[0]!=nvals ) {
105 0 : error("number of rows in target matrix should match number of input coordinates");
106 : }
107 2 : if( i==1 && target[0]->getRank()==1 ) {
108 : ntoproj = 1;
109 1 : } else if( ntoproj==1 && target[0]->getRank()!=1 ) {
110 0 : error("mismatch between numbers of target distances");
111 1 : } else if( i==1 ) {
112 1 : ntoproj = target[0]->getShape()[1];
113 0 : } else if( ntoproj!=target[0]->getShape()[1] ) {
114 0 : error("mismatch between numbers of target distances");
115 : }
116 4 : if( !parseArgumentList("WEIGHTS",i,weights) ) {
117 0 : error("missing WEIGHTS" + inum + " keyword in input");
118 : }
119 2 : if( weights.size()!=1 ) {
120 0 : error("should only be one value in input to WEIGHTS" + inum );
121 : }
122 2 : if( weights[0]->getRank()!=1 || weights[0]->hasDerivatives() ) {
123 0 : error("input to WEIGHTS" + inum + " keyword should be a vector");
124 : }
125 2 : if( weights[0]->getShape()[0]!=nvals ) {
126 0 : error("number of weights should match number of input coordinates");
127 : }
128 2 : target[0]->buildDataStore();
129 2 : weights[0]->buildDataStore();
130 2 : args.push_back( target[0] );
131 2 : args.push_back( weights[0] );
132 2 : bool has_sf = parseNumbered("FUNC",i,sfd);
133 2 : switchingFunction.push_back( SwitchingFunction() );
134 2 : if( !has_sf ) {
135 0 : switchingFunction[i-1].set( "CUSTOM FUNC=1-sqrt(x2) R_0=1.0", errors );
136 : } else {
137 2 : switchingFunction[i-1].set( sfd, errors );
138 2 : if( errors.length()!=0 ) {
139 0 : error("problem reading switching function description " + errors);
140 : }
141 : }
142 2 : log.printf(" %sth term seeks to match tranformed distances with those in matrix %s \n", inum.c_str(), target[0]->getName().c_str() );
143 2 : 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() );
144 2 : 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() );
145 2 : }
146 2 : std::vector<unsigned> shape(1);
147 2 : shape[0]=ntoproj;
148 2 : if( ntoproj==1 ) {
149 1 : shape.resize(0);
150 : }
151 6 : for(unsigned i=0; i<dimout; ++i) {
152 : std::string num;
153 4 : Tools::convert( i+1, num );
154 4 : addComponent( "coord-" + num, shape );
155 8 : componentIsNotPeriodic( "coord-" + num );
156 : }
157 : // Create a list of tasks to perform
158 2 : parse("CGTOL",cgtol);
159 2 : log.printf(" tolerance for conjugate gradient algorithm equals %f \n",cgtol);
160 2 : requestArguments( args );
161 2 : checkRead();
162 2 : }
163 :
164 169 : void ProjectPoints::prepare() {
165 169 : if( getPntrToComponent(0)->getRank()==0 ) {
166 168 : return;
167 : }
168 :
169 1 : std::vector<unsigned> shape(1);
170 1 : shape[0] = getPntrToArgument(dimout)->getShape()[0];
171 3 : for(unsigned i=0; i<dimout; ++i) {
172 2 : if( getPntrToComponent(i)->getShape()[0]!=shape[0] ) {
173 2 : getPntrToComponent(i)->setShape(shape);
174 : }
175 : }
176 : }
177 :
178 24180 : double ProjectPoints::calculateStress( const std::vector<double>& pp, std::vector<double>& der ) {
179 24180 : unsigned nmatrices = ( getNumberOfArguments() - dimout ) / 2;
180 : double stress=0;
181 :
182 24180 : unsigned t=OpenMP::getThreadNum();
183 24180 : std::vector<double> dtmp( pp.size() );
184 24180 : unsigned nland = getPntrToArgument(0)->getShape()[0];
185 4889418 : for(unsigned i=0; i<nland; ++i) {
186 : // Calculate distance in low dimensional space
187 : double dd2 = 0;
188 14595714 : for(unsigned k=0; k<pp.size(); ++k) {
189 9730476 : dtmp[k] = pp[k] - getPntrToArgument(k)->get(i);
190 9730476 : dd2 += dtmp[k]*dtmp[k];
191 : }
192 :
193 9730476 : for(unsigned k=0; k<nmatrices; ++k ) {
194 : // Now do transformations and calculate differences
195 4865238 : double df, fd = 1. - switchingFunction[k].calculateSqr( dd2, df );
196 : // Get the weight for this connection
197 4865238 : double weight = getPntrToArgument( dimout + 2*k + 1 )->get( i );
198 : // Get the difference for the connection
199 4865238 : double fdiff = fd - getPntrToArgument( dimout + 2*k )->get( rowstart[t]+i );
200 : // Calculate derivatives
201 4865238 : double pref = -2.*weight*fdiff*df;
202 14595714 : for(unsigned n=0; n<pp.size(); ++n) {
203 9730476 : der[n]+=pref*dtmp[n];
204 : }
205 : // Accumulate the total stress
206 4865238 : stress += weight*fdiff*fdiff;
207 : }
208 : }
209 24180 : return stress;
210 : }
211 :
212 668 : void ProjectPoints::getProjection( const unsigned& current, std::vector<double>& point ) const {
213 668 : Value* targ = getPntrToArgument( dimout );
214 668 : unsigned nland = getPntrToArgument(0)->getShape()[0];
215 668 : unsigned base = current;
216 668 : if( targ->getRank()==2 ) {
217 500 : base = current*targ->getShape()[1];
218 : }
219 : unsigned closest=0;
220 668 : double mindist = targ->get( base );
221 139112 : for(unsigned i=1; i<nland; ++i) {
222 138444 : double dist = targ->get( base + i );
223 138444 : if( dist<mindist ) {
224 : mindist=dist;
225 : closest=i;
226 : }
227 : }
228 : // Put the initial guess near to the closest landmark -- may wish to use grid here again Sandip??
229 668 : Random random;
230 668 : random.setSeed(-1234);
231 2004 : for(unsigned j=0; j<dimout; ++j) {
232 1336 : point[j] = getPntrToArgument(j)->get(closest) + (random.RandU01() - 0.5)*0.01;
233 : }
234 : // And do the optimisation
235 668 : rowstart[OpenMP::getThreadNum()]=current;
236 668 : if( targ->getRank()==2 ) {
237 500 : rowstart[OpenMP::getThreadNum()] = current*targ->getShape()[1];
238 : }
239 668 : myminimiser.minimise( cgtol, point, &ProjectPoints::calculateStress );
240 668 : }
241 :
242 500 : void ProjectPoints::performTask( const unsigned& current, MultiValue& myvals ) const {
243 500 : std::vector<double> point( dimout );
244 500 : getProjection( current, point );
245 1500 : for(unsigned j=0; j<dimout; ++j) {
246 1000 : myvals.setValue( getConstPntrToComponent(j)->getPositionInStream(), point[j] );
247 : }
248 500 : }
249 :
250 169 : void ProjectPoints::calculate() {
251 169 : if( getPntrToComponent(0)->getRank()==0 ) {
252 168 : std::vector<double> point( dimout );
253 168 : getProjection( 0, point );
254 504 : for(unsigned i=0; i<dimout; ++i) {
255 336 : getPntrToComponent(i)->set(point[i]);
256 : }
257 : } else {
258 1 : runAllTasks();
259 : }
260 169 : }
261 :
262 : }
263 : }
|
