Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2016-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 "reference/ReferenceConfiguration.h"
25 : #include "reference/MetricRegister.h"
26 : #include "core/PlumedMain.h"
27 : #include "core/Atoms.h"
28 : #include "tools/PDB.h"
29 :
30 : //+PLUMEDOC DIMRED SKETCHMAP_READ
31 : /*
32 : Read in a sketch-map projection from an input file
33 :
34 : \par Examples
35 :
36 : */
37 : //+ENDPLUMEDOC
38 :
39 : namespace PLMD {
40 : namespace dimred {
41 :
42 : class SketchMapRead : public SketchMapBase {
43 : private:
44 : std::string mtype;
45 : PDB mypdb;
46 : std::vector<Value*> all_values;
47 : std::vector<double> weights;
48 : /// The list of properties in the property map
49 : std::map<std::string,std::vector<double> > property;
50 : /// The data collection objects we have
51 : std::vector<analysis::DataCollectionObject> data;
52 : /// The frames that we are using to calculate distances
53 : std::vector<std::unique_ptr<ReferenceConfiguration> > myframes;
54 : public:
55 : static void registerKeywords( Keywords& keys );
56 : explicit SketchMapRead( const ActionOptions& ao );
57 : void minimise( Matrix<double>& ) override;
58 : analysis::DataCollectionObject& getStoredData( const unsigned& idata, const bool& calcdist ) override;
59 : unsigned getNumberOfDataPoints() const override;
60 : std::vector<Value*> getArgumentList() override;
61 : unsigned getDataPointIndexInBase( const unsigned& idata ) const override;
62 : double getDissimilarity( const unsigned& i, const unsigned& j ) override;
63 : double getWeight( const unsigned& idata ) override;
64 : };
65 :
66 13787 : PLUMED_REGISTER_ACTION(SketchMapRead,"SKETCHMAP_READ")
67 :
68 5 : void SketchMapRead::registerKeywords( Keywords& keys ) {
69 5 : SketchMapBase::registerKeywords( keys );
70 5 : keys.remove("USE_OUTPUT_DATA_FROM");
71 10 : keys.add("compulsory","TYPE","OPTIMAL-FAST","the manner in which distances are calculated. More information on the different "
72 : "metrics that are available in PLUMED can be found in the section of the manual on "
73 : "\\ref dists");
74 10 : keys.add("compulsory","REFERENCE","the file containing the sketch-map projection");
75 10 : keys.add("compulsory","PROPERTY","the property to be used in the index. This should be in the REMARK of the reference");
76 10 : keys.addFlag("DISABLE_CHECKS",false,"disable checks on reference input structures.");
77 5 : useCustomisableComponents(keys);
78 5 : }
79 :
80 1 : SketchMapRead::SketchMapRead( const ActionOptions& ao ):
81 : Action(ao),
82 1 : SketchMapBase(ao) {
83 : std::vector<std::string> propnames;
84 2 : parseVector("PROPERTY",propnames);
85 1 : if(propnames.size()==0) {
86 0 : error("no properties were specified");
87 : }
88 1 : nlow=propnames.size();
89 3 : for(unsigned i=0; i<nlow; ++i) {
90 4 : std::size_t dot=propnames[i].find_first_of( getLabel() + "." );
91 2 : std::string substr=propnames[i].c_str();
92 2 : if( dot!=std::string::npos ) {
93 2 : substr.erase(dot,getLabel().length()+1);
94 : }
95 2 : log.printf(",%s", propnames[i].c_str() );
96 2 : addComponent( substr );
97 2 : componentIsNotPeriodic( substr );
98 4 : property.insert( std::pair<std::string,std::vector<double> >( propnames[i], std::vector<double>() ) );
99 : }
100 1 : log.printf(" mapped properties are %s ",propnames[0].c_str() );
101 2 : for(unsigned i=1; i<nlow; ++i) {
102 1 : log.printf(",%s", propnames[i].c_str() );
103 : }
104 1 : log.printf("\n");
105 :
106 1 : parse("TYPE",mtype);
107 : bool skipchecks;
108 2 : parseFlag("DISABLE_CHECKS",skipchecks);
109 : std::string ifilename;
110 2 : parse("REFERENCE",ifilename);
111 1 : FILE* fp=this->fopen(ifilename.c_str(),"r");
112 1 : if(!fp) {
113 0 : error("could not open reference file " + ifilename );
114 : }
115 : // call fclose when fp goes out of scope
116 1 : auto deleter=[this](FILE* f) {
117 1 : this->fclose(f);
118 1 : };
119 : std::unique_ptr<FILE,decltype(deleter)> fp_deleter(fp,deleter);
120 :
121 :
122 : // Read in the embedding
123 : bool do_read=true;
124 : double val, ww, wnorm=0, prop;
125 : unsigned nfram=0;
126 85 : while (do_read) {
127 85 : PDB inpdb;
128 : // Read the pdb file
129 170 : do_read=inpdb.readFromFilepointer(fp,plumed.getAtoms().usingNaturalUnits(),0.1/atoms.getUnits().getLength());
130 : // Break if we are done
131 85 : if( !do_read ) {
132 : break ;
133 : }
134 : // Check for required properties
135 252 : for(std::map<std::string,std::vector<double> >::iterator it=property.begin(); it!=property.end(); ++it) {
136 168 : if( !inpdb.getArgumentValue( it->first, prop ) ) {
137 0 : error("pdb input does not have contain property named " + it->first );
138 : }
139 168 : it->second.push_back(prop);
140 : }
141 : // And read the frame ( create a measure )
142 84 : myframes.emplace_back( metricRegister().create<ReferenceConfiguration>( mtype, inpdb ) );
143 168 : if( !inpdb.getArgumentValue( "WEIGHT", ww ) ) {
144 0 : error("could not find weights in input pdb");
145 : }
146 84 : weights.push_back( ww );
147 84 : wnorm += ww;
148 84 : nfram++;
149 : // And create a data collection object
150 : analysis::DataCollectionObject new_data;
151 84 : new_data.setAtomNumbersAndArgumentNames( getLabel(), inpdb.getAtomNumbers(), inpdb.getArgumentNames() );
152 84 : new_data.setAtomPositions( inpdb.getPositions() );
153 504 : for(unsigned i=0; i<inpdb.getArgumentNames().size(); ++i) {
154 420 : std::string aname = inpdb.getArgumentNames()[i];
155 420 : if( !inpdb.getArgumentValue( aname, val ) ) {
156 0 : error("failed to find argument named " + aname);
157 : }
158 420 : new_data.setArgument( aname, val );
159 : }
160 84 : data.push_back( new_data );
161 85 : }
162 : // Finish the setup of the object by getting the arguments and atoms that are required
163 : std::vector<AtomNumber> atoms;
164 : std::vector<std::string> args;
165 85 : for(unsigned i=0; i<myframes.size(); ++i) {
166 84 : weights[i] /= wnorm;
167 84 : myframes[i]->getAtomRequests( atoms, skipchecks );
168 84 : myframes[i]->getArgumentRequests( args, skipchecks );
169 : }
170 1 : requestAtoms( atoms );
171 : std::vector<Value*> req_args;
172 : std::vector<std::string> fargs;
173 6 : for(unsigned i=0; i<args.size(); ++i) {
174 : bool found=false;
175 12 : for(std::map<std::string,std::vector<double> >::iterator it=property.begin(); it!=property.end(); ++it) {
176 9 : if( args[i]==it->first ) {
177 : found=true;
178 : break;
179 : }
180 : }
181 5 : if( !found ) {
182 3 : fargs.push_back( args[i] );
183 : }
184 : }
185 1 : interpretArgumentList( fargs, req_args );
186 1 : mypdb.setArgumentNames( fargs );
187 1 : requestArguments( req_args );
188 :
189 1 : if(nfram==0 ) {
190 0 : error("no reference configurations were specified");
191 : }
192 1 : log.printf(" found %u configurations in file %s\n",nfram,ifilename.c_str() );
193 4 : }
194 :
195 1 : void SketchMapRead::minimise( Matrix<double>& projections ) {
196 : unsigned j=0;
197 3 : for(std::map<std::string,std::vector<double> >::iterator it=property.begin(); it!=property.end(); ++it) {
198 170 : for(unsigned i=0; i<myframes.size(); ++i) {
199 168 : projections(i,j) = it->second[i];
200 : }
201 2 : j++;
202 : }
203 1 : }
204 :
205 7056 : analysis::DataCollectionObject & SketchMapRead::getStoredData( const unsigned& idata, const bool& calcdist ) {
206 7056 : return data[idata];
207 : }
208 :
209 173 : unsigned SketchMapRead::getNumberOfDataPoints() const {
210 173 : return myframes.size();
211 : }
212 :
213 0 : unsigned SketchMapRead::getDataPointIndexInBase( const unsigned& idata ) const {
214 0 : error("cannot use read in sketch-map to out of sample project data");
215 : return idata;
216 : }
217 :
218 0 : std::vector<Value*> SketchMapRead::getArgumentList() {
219 0 : std::vector<Value*> arglist( ActionWithArguments::getArguments() );
220 0 : for(unsigned i=0; i<nlow; ++i) {
221 0 : arglist.push_back( getPntrToComponent(i) );
222 : }
223 0 : return arglist;
224 : }
225 :
226 : // Highly unsatisfactory solution to problem GAT
227 3486 : double SketchMapRead::getDissimilarity( const unsigned& i, const unsigned& j ) {
228 3486 : plumed_assert( i<myframes.size() && j<myframes.size() );
229 3486 : if( i!=j ) {
230 : double dd;
231 3486 : getStoredData( i, true ).transferDataToPDB( mypdb );
232 3486 : auto myref1=metricRegister().create<ReferenceConfiguration>(mtype, mypdb);
233 3486 : getStoredData( j, true ).transferDataToPDB( mypdb );
234 3486 : auto myref2=metricRegister().create<ReferenceConfiguration>(mtype, mypdb);
235 3486 : dd=distance( getPbc(), ActionWithArguments::getArguments(), myref1.get(), myref2.get(), true );
236 : return dd;
237 3486 : }
238 : return 0.0;
239 : }
240 :
241 168 : double SketchMapRead::getWeight( const unsigned& idata ) {
242 168 : plumed_assert( idata<weights.size() );
243 168 : return weights[idata];
244 : }
245 :
246 : }
247 : }
|
