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