Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2016-2018 The VES code team
3 : (see the PEOPLE-VES file at the root of this folder for a list of names)
4 :
5 : See http://www.ves-code.org for more information.
6 :
7 : This file is part of VES code module.
8 :
9 : The VES code module 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 : The VES code module 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 the VES code module. If not, see <http://www.gnu.org/licenses/>.
21 : +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
22 :
23 : #include "TargetDistribution.h"
24 :
25 : #include "core/ActionRegister.h"
26 :
27 :
28 : namespace PLMD {
29 : namespace ves {
30 :
31 : //+PLUMEDOC VES_TARGETDIST TD_CHISQUARED
32 : /*
33 : Chi-squared distribution (static).
34 :
35 : Employ a target distribution given by a
36 : [chi-squared distribution](https://en.wikipedia.org/wiki/Chi-squared_distribution)
37 : that is defined as
38 : \f[
39 : p(s) =
40 : \frac
41 : {1}
42 : {\sigma \, 2^{\frac{k}{2}} \, \Gamma\left(\frac{k}{2}\right) }
43 : \, \left(\frac{s-a}{\sigma}\right)^{\frac{k}{2}-1} \, \exp\left(- \frac{1}{2}
44 : \left(\frac{s-a}{\sigma}\right) \right),
45 : \f]
46 : where \f$a\f$ is the minimum of the distribution that is defined on the interval \f$[a,\infty)\f$,
47 : the parameter \f$k\f$ (given as a postive integer larger than 2) determines how far
48 : the peak of the distribution is from the minimum (known as the "degrees of freedom"),
49 : and the parameter \f$\sigma>0\f$ determines the broadness of the distribution.
50 :
51 : The minimum \f$a\f$ is given using the MINIMUM keyword, the parameter \f$k\f$ is given
52 : using the KAPPA keyword, and the parameter \f$\sigma\f$ is given using the SIGMA keyword.
53 :
54 : This target distribution action is only defined for one dimension, for multiple dimensions
55 : it should be used in combination with the \ref TD_PRODUCT_DISTRIBUTION action.
56 :
57 : \par Examples
58 :
59 : Chi-squared distribution with \f$a=-10.0\f$, \f$\sigma=2.0\f$, and \f$k=2\f$
60 : \plumedfile
61 : td: TD_CHISQUARED MINIMUM=-10.0 SIGMA=2.0 KAPPA=2
62 : \endplumedfile
63 :
64 : The Chi-squared distribution is only defined for one dimension so for multiple
65 : dimensions we have to use it in combination with the \ref TD_PRODUCT_DISTRIBUTION action as shown in
66 : the following example where we have a Chi-squared distribution for argument 1
67 : and uniform distribution for argument 2
68 : \plumedfile
69 : td_chisq: TD_CHISQUARED MINIMUM=10.0 SIGMA=2.0 KAPPA=2
70 :
71 : td_uni: TD_UNIFORM
72 :
73 : td_pd: TD_PRODUCT_DISTRIBUTION DISTRIBUTIONS=td_chisq,td_uni
74 : \endplumedfile
75 :
76 : */
77 : //+ENDPLUMEDOC
78 :
79 27 : class TD_ChiSquared: public TargetDistribution {
80 : std::vector<double> minima_;
81 : std::vector<double> sigma_;
82 : std::vector<double> kappa_;
83 : std::vector<double> normalization_;
84 : public:
85 : static void registerKeywords(Keywords&);
86 : explicit TD_ChiSquared(const ActionOptions& ao);
87 : double getValue(const std::vector<double>&) const;
88 : };
89 :
90 :
91 6461 : PLUMED_REGISTER_ACTION(TD_ChiSquared,"TD_CHISQUARED")
92 :
93 :
94 10 : void TD_ChiSquared::registerKeywords(Keywords& keys) {
95 10 : TargetDistribution::registerKeywords(keys);
96 40 : keys.add("compulsory","MINIMUM","The minimum of the chi-squared distribution.");
97 40 : keys.add("compulsory","SIGMA","The \\f$\\sigma\\f$ parameter of the chi-squared distribution given as a postive number.");
98 40 : keys.add("compulsory","KAPPA","The \\f$k\\f$ parameter of the chi-squared distribution given as postive integer larger than 2.");
99 20 : keys.use("WELLTEMPERED_FACTOR");
100 20 : keys.use("SHIFT_TO_ZERO");
101 20 : keys.use("NORMALIZE");
102 10 : }
103 :
104 :
105 9 : TD_ChiSquared::TD_ChiSquared(const ActionOptions& ao):
106 : PLUMED_VES_TARGETDISTRIBUTION_INIT(ao),
107 : minima_(0),
108 : sigma_(0),
109 : kappa_(0),
110 9 : normalization_(0)
111 : {
112 18 : parseVector("MINIMUM",minima_);
113 18 : parseVector("SIGMA",sigma_);
114 45 : for(unsigned int k=0; k<sigma_.size(); k++) {
115 9 : if(sigma_[k] < 0.0) {plumed_merror(getName()+": the value given in SIGMA should be postive.");}
116 : }
117 :
118 9 : std::vector<unsigned int> kappa_int(0);
119 18 : parseVector("KAPPA",kappa_int);
120 9 : if(kappa_int.size()==0) {plumed_merror(getName()+": some problem with KAPPA keyword, should given as postive integer larger than 2");}
121 9 : kappa_.resize(kappa_int.size());
122 45 : for(unsigned int k=0; k<kappa_int.size(); k++) {
123 9 : if(kappa_int[k] < 2) {plumed_merror(getName()+": KAPPA should be an integer 2 or higher");}
124 9 : kappa_[k] = static_cast<double>(kappa_int[k]);
125 : }
126 :
127 9 : setDimension(minima_.size());
128 9 : if(getDimension()>1) {plumed_merror(getName()+": only defined for one dimension, for multiple dimensions it should be used in combination with the TD_PRODUCT_DISTRIBUTION action.");}
129 9 : if(sigma_.size()!=getDimension()) {plumed_merror(getName()+": the SIGMA keyword does not match the given dimension in MINIMUM");}
130 9 : if(kappa_.size()!=getDimension()) {plumed_merror(getName()+": the KAPPA keyword does not match the given dimension in MINIMUM");}
131 :
132 9 : normalization_.resize(getDimension());
133 27 : for(unsigned int k=0; k<getDimension(); k++) {
134 45 : normalization_[k] = 1.0/(pow(2.0,0.5*kappa_[k])*tgamma(0.5*kappa_[k])*sigma_[k]);
135 : }
136 9 : checkRead();
137 9 : }
138 :
139 :
140 1509 : double TD_ChiSquared::getValue(const std::vector<double>& argument) const {
141 : double value = 1.0;
142 7545 : for(unsigned int k=0; k<argument.size(); k++) {
143 4527 : double arg=(argument[k]-minima_[k])/sigma_[k];
144 1509 : if(arg<0.0) {plumed_merror(getName()+": the chi-squared istribution is not defined for values less that ones given in MINIMUM");}
145 3018 : value *= normalization_[k] * pow(arg,0.5*kappa_[k]-1.0) * exp(-0.5*arg);
146 : }
147 1509 : return value;
148 : }
149 :
150 :
151 : }
152 4839 : }
|
