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