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 "core/ActionRegister.h"
23 : #include "SketchMapBase.h"
24 : #include "tools/ConjugateGradient.h"
25 : #include "gridtools/GridSearch.h"
26 :
27 : //+PLUMEDOC DIMRED SKETCHMAP_POINTWISE
28 : /*
29 : Optimize the sketch-map stress function using a pointwise global optimization algorithm.
30 :
31 : \par Examples
32 :
33 : */
34 : //+ENDPLUMEDOC
35 :
36 : namespace PLMD {
37 : namespace dimred {
38 :
39 : class SketchMapPointwise : public SketchMapBase {
40 : private:
41 : unsigned ncycles;
42 : double cgtol, gbuf;
43 : std::vector<unsigned> npoints, nfgrid;
44 : public:
45 : static void registerKeywords( Keywords& keys );
46 : explicit SketchMapPointwise( const ActionOptions& ao );
47 : void minimise( Matrix<double>& ) override;
48 : };
49 :
50 10421 : PLUMED_REGISTER_ACTION(SketchMapPointwise,"SKETCHMAP_POINTWISE")
51 :
52 2 : void SketchMapPointwise::registerKeywords( Keywords& keys ) {
53 2 : SketchMapBase::registerKeywords( keys );
54 4 : keys.add("compulsory","NCYCLES","5","the number of cycles of global optimization to attempt");
55 4 : keys.add("compulsory","CGTOL","1E-6","the tolerance for the conjugate gradient minimization");
56 4 : keys.add("compulsory","BUFFER","1.1","grid extent for search is (max projection - minimum projection) multiplied by this value");
57 4 : keys.add("compulsory","CGRID_SIZE","10","number of points to use in each grid direction");
58 4 : keys.add("compulsory","FGRID_SIZE","0","interpolate the grid onto this number of points -- only works in 2D");
59 2 : }
60 :
61 1 : SketchMapPointwise::SketchMapPointwise( const ActionOptions& ao ):
62 : Action(ao),
63 : SketchMapBase(ao),
64 2 : npoints(nlow),
65 1 : nfgrid(nlow)
66 : {
67 3 : parseVector("CGRID_SIZE",npoints); parse("BUFFER",gbuf);
68 1 : if( npoints.size()!=nlow ) error("vector giving number of grid point in each direction has wrong size");
69 2 : parse("NCYCLES",ncycles); parse("CGTOL",cgtol);
70 :
71 2 : parseVector("FGRID_SIZE",nfgrid);
72 1 : if( nfgrid[0]!=0 && nlow!=2 ) error("interpolation only works in two dimensions");
73 :
74 1 : log.printf(" doing %u cycles of global optimization sweeps\n",ncycles);
75 1 : log.printf(" using coarse grid of points that is %u",npoints[0]);
76 1 : log.printf(" and that is %f larger than the difference between the position of the minimum and maximum projection \n",gbuf);
77 2 : for(unsigned j=1; j<npoints.size(); ++j) log.printf(" by %u",npoints[j]);
78 1 : if( nfgrid[0]>0 ) {
79 1 : log.printf(" interpolating stress onto grid of points that is %u",nfgrid[0]);
80 2 : for(unsigned j=1; j<nfgrid.size(); ++j) log.printf(" by %u",nfgrid[j]);
81 1 : log.printf("\n");
82 : }
83 1 : log.printf(" tolerance for conjugate gradient algorithm equals %f \n",cgtol);
84 1 : }
85 :
86 1 : void SketchMapPointwise::minimise( Matrix<double>& projections ) {
87 1 : std::vector<double> gmin( nlow ), gmax( nlow ), mypoint( nlow );
88 :
89 : // Find the extent of the grid
90 3 : for(unsigned j=0; j<nlow; ++j) gmin[j]=gmax[j]=projections(0,j);
91 100 : for(unsigned i=1; i<getNumberOfDataPoints(); ++i) {
92 297 : for(unsigned j=0; j<nlow; ++j) {
93 198 : if( projections(i,j) < gmin[j] ) gmin[j] = projections(i,j);
94 198 : if( projections(i,j) > gmax[j] ) gmax[j] = projections(i,j);
95 : }
96 : }
97 3 : for(unsigned j=0; j<nlow; ++j) {
98 2 : double gbuffer = 0.5*gbuf*( gmax[j]-gmin[j] ) - 0.5*( gmax[j]- gmin[j] );
99 2 : gmin[j]-=gbuffer; gmax[j]+=gbuffer;
100 : }
101 :
102 : // And do the search
103 : ConjugateGradient<SketchMapPointwise> mycgminimise( this );
104 1 : gridtools::GridSearch<SketchMapPointwise> mygridsearch( gmin, gmax, npoints, nfgrid, this );
105 : // Run multiple loops over all projections
106 3 : for(unsigned i=0; i<ncycles; ++i) {
107 202 : for(unsigned j=0; j<getNumberOfDataPoints(); ++j) {
108 : // Setup target distances and target functions for calculate stress
109 20200 : for(unsigned k=0; k<getNumberOfDataPoints(); ++k) setTargetDistance( k, distances(j,k) );
110 :
111 : // Find current projection of jth point
112 600 : for(unsigned k=0; k<mypoint.size(); ++k) mypoint[k]=projections(j,k);
113 : // Minimise using grid search
114 200 : bool moved=mygridsearch.minimise( mypoint, &SketchMapPointwise::calculateStress );
115 200 : if( moved ) {
116 : // Reassign the new projection
117 66 : for(unsigned k=0; k<mypoint.size(); ++k) projections(j,k)=mypoint[k];
118 : // Minimise output using conjugate gradient
119 22 : mycgminimise.minimise( cgtol, projections.getVector(), &SketchMapPointwise::calculateFullStress );
120 : }
121 : }
122 : }
123 1 : }
124 :
125 : }
126 : }
|
