Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2015-2019 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 "Function.h"
23 : #include "ActionRegister.h"
24 :
25 : using namespace std;
26 :
27 : namespace PLMD {
28 : namespace function {
29 :
30 : //+PLUMEDOC FUNCTION STATS
31 : /*
32 : Calculates statistical properties of a set of collective variables with respect to a set of reference values.
33 : In particular it calculates and store as components the sum of the squared deviations, the correlation, the
34 : slope and the intercept of a linear fit.
35 :
36 : The reference values can be either provided as values using PARAMETERS or using value without derivatives
37 : from other actions using PARARG (for example using experimental values from collective variables such as
38 : \ref CS2BACKBONE, \ref RDC, \ref NOE, \ref PRE).
39 :
40 : \par Examples
41 :
42 : The following input tells plumed to print the distance between three couple of atoms
43 : and compare them with three reference distances.
44 :
45 : \plumedfile
46 : d1: DISTANCE ATOMS=10,50
47 : d2: DISTANCE ATOMS=1,100
48 : d3: DISTANCE ATOMS=45,75
49 : st: STATS ARG=d1,d2,d3 PARAMETERS=1.5,4.0,2.0
50 : PRINT ARG=d1,d2,d3,st.*
51 : \endplumedfile
52 :
53 : */
54 : //+ENDPLUMEDOC
55 :
56 :
57 93 : class Stats :
58 : public Function
59 : {
60 : std::vector<double> parameters;
61 : bool sqdonly;
62 : bool components;
63 : bool upperd;
64 : public:
65 : explicit Stats(const ActionOptions&);
66 : void calculate();
67 : static void registerKeywords(Keywords& keys);
68 : };
69 :
70 :
71 6483 : PLUMED_REGISTER_ACTION(Stats,"STATS")
72 :
73 32 : void Stats::registerKeywords(Keywords& keys) {
74 32 : Function::registerKeywords(keys);
75 32 : useCustomisableComponents(keys);
76 64 : keys.use("ARG");
77 128 : keys.add("optional","PARARG","the input for this action is the scalar output from one or more other actions without derivatives.");
78 128 : keys.add("optional","PARAMETERS","the parameters of the arguments in your function");
79 96 : keys.addFlag("SQDEVSUM",false,"calculates only SQDEVSUM");
80 96 : keys.addFlag("SQDEV",false,"calculates and store the SQDEV as components");
81 96 : keys.addFlag("UPPERDISTS",false,"calculates and store the SQDEV as components");
82 128 : keys.addOutputComponent("sqdevsum","default","the sum of the squared deviations between arguments and parameters");
83 128 : keys.addOutputComponent("corr","default","the correlation between arguments and parameters");
84 128 : keys.addOutputComponent("slope","default","the slope of a linear fit between arguments and parameters");
85 128 : keys.addOutputComponent("intercept","default","the intercept of a linear fit between arguments and parameters");
86 128 : keys.addOutputComponent("sqd","SQDEV","the squared deviations between arguments and parameters");
87 32 : }
88 :
89 31 : Stats::Stats(const ActionOptions&ao):
90 : Action(ao),
91 : Function(ao),
92 : sqdonly(false),
93 : components(false),
94 62 : upperd(false)
95 : {
96 62 : parseVector("PARAMETERS",parameters);
97 31 : if(parameters.size()!=static_cast<unsigned>(getNumberOfArguments())&&!parameters.empty())
98 0 : error("Size of PARAMETERS array should be either 0 or the same as of the number of arguments in ARG1");
99 :
100 : vector<Value*> arg2;
101 62 : parseArgumentList("PARARG",arg2);
102 :
103 31 : if(!arg2.empty()) {
104 14 : if(parameters.size()>0) error("It is not possible to use PARARG and PARAMETERS together");
105 14 : if(arg2.size()!=getNumberOfArguments()) error("Size of PARARG array should be the same as number for arguments in ARG");
106 17722 : for(unsigned i=0; i<arg2.size(); i++) {
107 17694 : parameters.push_back(arg2[i]->get());
108 11796 : if(arg2[i]->hasDerivatives()==true) error("PARARG can only accept arguments without derivatives");
109 : }
110 : }
111 :
112 31 : if(parameters.size()!=getNumberOfArguments())
113 0 : error("PARARG or PARAMETERS arrays should include the same number of elements as the arguments in ARG");
114 :
115 31 : if(getNumberOfArguments()<2) error("STATS need at least two arguments to be used");
116 :
117 62 : parseFlag("SQDEVSUM",sqdonly);
118 62 : parseFlag("SQDEV",components);
119 62 : parseFlag("UPPERDISTS",upperd);
120 :
121 31 : if(sqdonly&&components) error("You cannot used SQDEVSUM and SQDEV at the sametime");
122 :
123 31 : if(components) sqdonly = true;
124 :
125 45 : if(!arg2.empty()) log.printf(" using %zu parameters from inactive actions:", arg2.size());
126 34 : else log.printf(" using %zu parameters:", arg2.size());
127 17969 : for(unsigned i=0; i<parameters.size(); i++) log.printf(" %f",parameters[i]);
128 31 : log.printf("\n");
129 :
130 31 : if(sqdonly) {
131 17 : if(components) {
132 168 : for(unsigned i=0; i<parameters.size(); i++) {
133 48 : std::string num; Tools::convert(i,num);
134 96 : addComponentWithDerivatives("sqd_"+num);
135 96 : componentIsNotPeriodic("sqd_"+num);
136 : }
137 : } else {
138 10 : addComponentWithDerivatives("sqdevsum");
139 10 : componentIsNotPeriodic("sqdevsum");
140 : }
141 : } else {
142 28 : addComponentWithDerivatives("sqdevsum");
143 28 : componentIsNotPeriodic("sqdevsum");
144 28 : addComponentWithDerivatives("corr");
145 28 : componentIsNotPeriodic("corr");
146 28 : addComponentWithDerivatives("slope");
147 28 : componentIsNotPeriodic("slope");
148 28 : addComponentWithDerivatives("intercept");
149 28 : componentIsNotPeriodic("intercept");
150 : }
151 :
152 :
153 31 : checkRead();
154 31 : }
155 :
156 122 : void Stats::calculate()
157 : {
158 122 : if(sqdonly) {
159 :
160 : double nsqd = 0.;
161 : Value* val;
162 106 : if(!components) val=getPntrToComponent("sqdevsum");
163 469 : for(unsigned i=0; i<parameters.size(); ++i) {
164 121 : double dev = getArgument(i)-parameters[i];
165 121 : if(upperd&&dev<0) dev=0.;
166 121 : if(components) {
167 0 : val=getPntrToComponent(i);
168 0 : val->set(dev*dev);
169 : } else {
170 121 : nsqd += dev*dev;
171 : }
172 121 : setDerivative(val,i,2.*dev);
173 : }
174 53 : if(!components) val->set(nsqd);
175 :
176 : } else {
177 :
178 : double scx=0., scx2=0., scy=0., scy2=0., scxy=0.;
179 :
180 18621 : for(unsigned i=0; i<parameters.size(); ++i) {
181 : const double tmpx=getArgument(i);
182 6161 : const double tmpy=parameters[i];
183 6161 : scx += tmpx;
184 6161 : scx2 += tmpx*tmpx;
185 6161 : scy += tmpy;
186 6161 : scy2 += tmpy*tmpy;
187 6161 : scxy += tmpx*tmpy;
188 : }
189 :
190 69 : const double ns = parameters.size();
191 :
192 69 : const double num = ns*scxy - scx*scy;
193 69 : const double idev2x = 1./(ns*scx2-scx*scx);
194 69 : const double idevx = sqrt(idev2x);
195 69 : const double idevy = 1./sqrt(ns*scy2-scy*scy);
196 :
197 : /* sd */
198 69 : const double nsqd = scx2 + scy2 - 2.*scxy;
199 : /* correlation */
200 69 : const double correlation = num * idevx * idevy;
201 : /* slope and intercept */
202 69 : const double slope = num * idev2x;
203 69 : const double inter = (scy - slope * scx)/ns;
204 :
205 138 : Value* valuea=getPntrToComponent("sqdevsum");
206 138 : Value* valueb=getPntrToComponent("corr");
207 138 : Value* valuec=getPntrToComponent("slope");
208 138 : Value* valued=getPntrToComponent("intercept");
209 :
210 : valuea->set(nsqd);
211 : valueb->set(correlation);
212 : valuec->set(slope);
213 : valued->set(inter);
214 :
215 : /* derivatives */
216 18621 : for(unsigned i=0; i<parameters.size(); ++i) {
217 6161 : const double common_d1 = (ns*parameters[i]-scy)*idevx;
218 6161 : const double common_d2 = num*(ns*getArgument(i)-scx)*idev2x*idevx;
219 6161 : const double common_d3 = common_d1 - common_d2;
220 :
221 : /* sqdevsum */
222 6161 : const double sq_der = 2.*(getArgument(i)-parameters[i]);
223 : /* correlation */
224 6161 : const double co_der = common_d3*idevy;
225 : /* slope */
226 6161 : const double sl_der = (common_d1-2.*common_d2)*idevx;
227 : /* intercept */
228 6161 : const double int_der = -(slope+ scx*sl_der)/ns;
229 :
230 : setDerivative(valuea,i,sq_der);
231 : setDerivative(valueb,i,co_der);
232 : setDerivative(valuec,i,sl_der);
233 : setDerivative(valued,i,int_der);
234 : }
235 :
236 : }
237 122 : }
238 :
239 : }
240 4839 : }
241 :
242 :
|
