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 10421 : PLUMED_REGISTER_ACTION(SketchMapRead,"SKETCHMAP_READ")
67 :
68 2 : void SketchMapRead::registerKeywords( Keywords& keys ) {
69 2 : SketchMapBase::registerKeywords( keys ); keys.remove("USE_OUTPUT_DATA_FROM");
70 4 : keys.add("compulsory","TYPE","OPTIMAL-FAST","the manner in which distances are calculated. More information on the different "
71 : "metrics that are available in PLUMED can be found in the section of the manual on "
72 : "\\ref dists");
73 4 : keys.add("compulsory","REFERENCE","the file containing the sketch-map projection");
74 4 : keys.add("compulsory","PROPERTY","the property to be used in the index. This should be in the REMARK of the reference");
75 4 : keys.addFlag("DISABLE_CHECKS",false,"disable checks on reference input structures.");
76 2 : useCustomisableComponents(keys);
77 2 : }
78 :
79 1 : SketchMapRead::SketchMapRead( const ActionOptions& ao ):
80 : Action(ao),
81 1 : SketchMapBase(ao)
82 : {
83 2 : std::vector<std::string> propnames; parseVector("PROPERTY",propnames);
84 1 : if(propnames.size()==0) error("no properties were specified");
85 1 : nlow=propnames.size();
86 3 : for(unsigned i=0; i<nlow; ++i) {
87 4 : std::size_t dot=propnames[i].find_first_of( getLabel() + "." ); std::string substr=propnames[i].c_str();
88 2 : if( dot!=std::string::npos ) { substr.erase(dot,getLabel().length()+1); }
89 2 : log.printf(",%s", propnames[i].c_str() ); addComponent( substr ); componentIsNotPeriodic( substr );
90 4 : property.insert( std::pair<std::string,std::vector<double> >( propnames[i], std::vector<double>() ) );
91 : }
92 1 : log.printf(" mapped properties are %s ",propnames[0].c_str() );
93 2 : for(unsigned i=1; i<nlow; ++i) log.printf(",%s", propnames[i].c_str() );
94 1 : log.printf("\n");
95 :
96 3 : parse("TYPE",mtype); bool skipchecks; parseFlag("DISABLE_CHECKS",skipchecks);
97 2 : std::string ifilename; parse("REFERENCE",ifilename);
98 1 : FILE* fp=fopen(ifilename.c_str(),"r");
99 1 : if(!fp) error("could not open reference file " + ifilename );
100 :
101 : // Read in the embedding
102 : bool do_read=true; double val, ww, wnorm=0, prop; unsigned nfram=0;
103 85 : while (do_read) {
104 85 : PDB inpdb;
105 : // Read the pdb file
106 170 : do_read=inpdb.readFromFilepointer(fp,plumed.getAtoms().usingNaturalUnits(),0.1/atoms.getUnits().getLength());
107 : // Break if we are done
108 85 : if( !do_read ) break ;
109 : // Check for required properties
110 252 : for(std::map<std::string,std::vector<double> >::iterator it=property.begin(); it!=property.end(); ++it) {
111 168 : if( !inpdb.getArgumentValue( it->first, prop ) ) error("pdb input does not have contain property named " + it->first );
112 168 : it->second.push_back(prop);
113 : }
114 : // And read the frame ( create a measure )
115 84 : myframes.emplace_back( metricRegister().create<ReferenceConfiguration>( mtype, inpdb ) );
116 168 : if( !inpdb.getArgumentValue( "WEIGHT", ww ) ) error("could not find weights in input pdb");
117 84 : weights.push_back( ww ); wnorm += ww; nfram++;
118 : // And create a data collection object
119 84 : analysis::DataCollectionObject new_data; new_data.setAtomNumbersAndArgumentNames( getLabel(), inpdb.getAtomNumbers(), inpdb.getArgumentNames() );
120 84 : new_data.setAtomPositions( inpdb.getPositions() );
121 504 : for(unsigned i=0; i<inpdb.getArgumentNames().size(); ++i) {
122 420 : std::string aname = inpdb.getArgumentNames()[i];
123 420 : if( !inpdb.getArgumentValue( aname, val ) ) error("failed to find argument named " + aname);
124 420 : new_data.setArgument( aname, val );
125 : }
126 84 : data.push_back( new_data );
127 85 : }
128 1 : fclose(fp);
129 : // Finish the setup of the object by getting the arguments and atoms that are required
130 : std::vector<AtomNumber> atoms; std::vector<std::string> args;
131 85 : for(unsigned i=0; i<myframes.size(); ++i) { weights[i] /= wnorm; myframes[i]->getAtomRequests( atoms, skipchecks ); myframes[i]->getArgumentRequests( args, skipchecks ); }
132 1 : requestAtoms( atoms ); std::vector<Value*> req_args; std::vector<std::string> fargs;
133 6 : for(unsigned i=0; i<args.size(); ++i) {
134 : bool found=false;
135 12 : for(std::map<std::string,std::vector<double> >::iterator it=property.begin(); it!=property.end(); ++it) {
136 9 : if( args[i]==it->first ) { found=true; break; }
137 : }
138 5 : if( !found ) { fargs.push_back( args[i] ); }
139 : }
140 1 : interpretArgumentList( fargs, req_args ); mypdb.setArgumentNames( fargs ); requestArguments( req_args );
141 :
142 1 : if(nfram==0 ) error("no reference configurations were specified");
143 1 : log.printf(" found %u configurations in file %s\n",nfram,ifilename.c_str() );
144 3 : }
145 :
146 1 : void SketchMapRead::minimise( Matrix<double>& projections ) {
147 : unsigned j=0;
148 3 : for(std::map<std::string,std::vector<double> >::iterator it=property.begin(); it!=property.end(); ++it) {
149 170 : for(unsigned i=0; i<myframes.size(); ++i) projections(i,j) = it->second[i];
150 2 : j++;
151 : }
152 1 : }
153 :
154 7056 : analysis::DataCollectionObject & SketchMapRead::getStoredData( const unsigned& idata, const bool& calcdist ) {
155 7056 : return data[idata];
156 : }
157 :
158 173 : unsigned SketchMapRead::getNumberOfDataPoints() const {
159 173 : return myframes.size();
160 : }
161 :
162 0 : unsigned SketchMapRead::getDataPointIndexInBase( const unsigned& idata ) const {
163 0 : error("cannot use read in sketch-map to out of sample project data");
164 : return idata;
165 : }
166 :
167 0 : std::vector<Value*> SketchMapRead::getArgumentList() {
168 0 : std::vector<Value*> arglist( ActionWithArguments::getArguments() );
169 0 : for(unsigned i=0; i<nlow; ++i) arglist.push_back( getPntrToComponent(i) );
170 0 : return arglist;
171 : }
172 :
173 : // Highly unsatisfactory solution to problem GAT
174 3486 : double SketchMapRead::getDissimilarity( const unsigned& i, const unsigned& j ) {
175 3486 : plumed_assert( i<myframes.size() && j<myframes.size() );
176 3486 : if( i!=j ) {
177 : double dd;
178 3486 : getStoredData( i, true ).transferDataToPDB( mypdb );
179 3486 : auto myref1=metricRegister().create<ReferenceConfiguration>(mtype, mypdb);
180 3486 : getStoredData( j, true ).transferDataToPDB( mypdb );
181 3486 : auto myref2=metricRegister().create<ReferenceConfiguration>(mtype, mypdb);
182 3486 : dd=distance( getPbc(), ActionWithArguments::getArguments(), myref1.get(), myref2.get(), true );
183 : return dd;
184 3486 : }
185 : return 0.0;
186 : }
187 :
188 168 : double SketchMapRead::getWeight( const unsigned& idata ) {
189 168 : plumed_assert( idata<weights.size() );
190 168 : return weights[idata];
191 : }
192 :
193 : }
194 : }
|
