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 13787 : PLUMED_REGISTER_ACTION(SketchMapPointwise,"SKETCHMAP_POINTWISE")
51 :
52 5 : void SketchMapPointwise::registerKeywords( Keywords& keys ) {
53 5 : SketchMapBase::registerKeywords( keys );
54 10 : keys.add("compulsory","NCYCLES","5","the number of cycles of global optimization to attempt");
55 10 : keys.add("compulsory","CGTOL","1E-6","the tolerance for the conjugate gradient minimization");
56 10 : keys.add("compulsory","BUFFER","1.1","grid extent for search is (max projection - minimum projection) multiplied by this value");
57 10 : keys.add("compulsory","CGRID_SIZE","10","number of points to use in each grid direction");
58 10 : keys.add("compulsory","FGRID_SIZE","0","interpolate the grid onto this number of points -- only works in 2D");
59 5 : }
60 :
61 1 : SketchMapPointwise::SketchMapPointwise( const ActionOptions& ao ):
62 : Action(ao),
63 : SketchMapBase(ao),
64 2 : npoints(nlow),
65 1 : nfgrid(nlow) {
66 1 : parseVector("CGRID_SIZE",npoints);
67 2 : parse("BUFFER",gbuf);
68 1 : if( npoints.size()!=nlow ) {
69 0 : error("vector giving number of grid point in each direction has wrong size");
70 : }
71 1 : parse("NCYCLES",ncycles);
72 1 : parse("CGTOL",cgtol);
73 :
74 2 : parseVector("FGRID_SIZE",nfgrid);
75 1 : if( nfgrid[0]!=0 && nlow!=2 ) {
76 0 : error("interpolation only works in two dimensions");
77 : }
78 :
79 1 : log.printf(" doing %u cycles of global optimization sweeps\n",ncycles);
80 1 : log.printf(" using coarse grid of points that is %u",npoints[0]);
81 1 : log.printf(" and that is %f larger than the difference between the position of the minimum and maximum projection \n",gbuf);
82 2 : for(unsigned j=1; j<npoints.size(); ++j) {
83 1 : log.printf(" by %u",npoints[j]);
84 : }
85 1 : if( nfgrid[0]>0 ) {
86 1 : log.printf(" interpolating stress onto grid of points that is %u",nfgrid[0]);
87 2 : for(unsigned j=1; j<nfgrid.size(); ++j) {
88 1 : log.printf(" by %u",nfgrid[j]);
89 : }
90 1 : log.printf("\n");
91 : }
92 1 : log.printf(" tolerance for conjugate gradient algorithm equals %f \n",cgtol);
93 1 : }
94 :
95 1 : void SketchMapPointwise::minimise( Matrix<double>& projections ) {
96 1 : std::vector<double> gmin( nlow ), gmax( nlow ), mypoint( nlow );
97 :
98 : // Find the extent of the grid
99 3 : for(unsigned j=0; j<nlow; ++j) {
100 2 : gmin[j]=gmax[j]=projections(0,j);
101 : }
102 100 : for(unsigned i=1; i<getNumberOfDataPoints(); ++i) {
103 297 : for(unsigned j=0; j<nlow; ++j) {
104 198 : if( projections(i,j) < gmin[j] ) {
105 7 : gmin[j] = projections(i,j);
106 : }
107 198 : if( projections(i,j) > gmax[j] ) {
108 8 : gmax[j] = projections(i,j);
109 : }
110 : }
111 : }
112 3 : for(unsigned j=0; j<nlow; ++j) {
113 2 : double gbuffer = 0.5*gbuf*( gmax[j]-gmin[j] ) - 0.5*( gmax[j]- gmin[j] );
114 2 : gmin[j]-=gbuffer;
115 2 : gmax[j]+=gbuffer;
116 : }
117 :
118 : // And do the search
119 : ConjugateGradient<SketchMapPointwise> mycgminimise( this );
120 1 : gridtools::GridSearch<SketchMapPointwise> mygridsearch( gmin, gmax, npoints, nfgrid, this );
121 : // Run multiple loops over all projections
122 3 : for(unsigned i=0; i<ncycles; ++i) {
123 202 : for(unsigned j=0; j<getNumberOfDataPoints(); ++j) {
124 : // Setup target distances and target functions for calculate stress
125 20200 : for(unsigned k=0; k<getNumberOfDataPoints(); ++k) {
126 20000 : setTargetDistance( k, distances(j,k) );
127 : }
128 :
129 : // Find current projection of jth point
130 600 : for(unsigned k=0; k<mypoint.size(); ++k) {
131 400 : mypoint[k]=projections(j,k);
132 : }
133 : // Minimise using grid search
134 200 : bool moved=mygridsearch.minimise( mypoint, &SketchMapPointwise::calculateStress );
135 200 : if( moved ) {
136 : // Reassign the new projection
137 66 : for(unsigned k=0; k<mypoint.size(); ++k) {
138 44 : projections(j,k)=mypoint[k];
139 : }
140 : // Minimise output using conjugate gradient
141 22 : mycgminimise.minimise( cgtol, projections.getVector(), &SketchMapPointwise::calculateFullStress );
142 : }
143 : }
144 : }
145 1 : }
146 :
147 : }
148 : }
|
