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_GENERALIZED_EXTREME_VALUE
32 : /*
33 : Generalized extreme value distribution (static).
34 :
35 : Employ a target distribution given by a
36 : [generalized extreme value distribution](https://en.wikipedia.org/wiki/Generalized_extreme_value_distribution)
37 : that is defined as
38 : \f[
39 : p(s) =
40 : \frac{1}{\sigma} \, t(s)^{\xi+1} \, e^{-t(s)},
41 : \f]
42 : where
43 : \f[
44 : t(s) =
45 : \begin{cases}
46 : \left( 1 + \xi \left( \frac{s-\mu}{\sigma} \right) \right)^{-1/\xi} & \mathrm{if\ }\xi \neq 0 \\
47 : \exp\left(- \frac{s-\mu}{\sigma} \right) & \mathrm{if\ } \xi = 0
48 : \end{cases},
49 : \f]
50 : and \f$\mu\f$ is the location parameter which approximately determines the location of the
51 : maximum of the distribution, \f$\sigma>0\f$ is the scale parameter that determines the
52 : broadness of the distribution, and \f$\xi\f$ is the shape parameter that determines
53 : the tail behavior of the distribution. For \f$\xi=0\f$, \f$\xi>0\f$, and \f$\xi<0\f$
54 : the Gumbel, Frechet, and Weibull families of distributions are obtained, respectively.
55 :
56 : The location parameter \f$\mu\f$ is given using the LOCATION keyword, the scale parameter \f$\sigma\f$
57 : using the SCALE keyword, and the shape parameter \f$\xi\f$ using the SHAPE
58 : keyword.
59 :
60 : This target distribution action is only defined for one dimension, for multiple dimensions
61 : it should be used in combination with \ref TD_PRODUCT_DISTRIBUTION action.
62 :
63 : \par Examples
64 :
65 : Generalized extreme value distribution with \f$\mu=0.0\f$, \f$\sigma=2.0\f$, and \f$\xi=0.0\f$ (Gumbel distribution)
66 : \plumedfile
67 : td: TD_GENERALIZED_EXTREME_VALUE LOCATION=0.0 SCALE=2.0 SHAPE=0.0
68 : \endplumedfile
69 :
70 :
71 : Generalized extreme value distribution with \f$\mu=-5.0\f$, \f$\sigma=1.0\f$, and \f$\xi=0.5\f$ (Frechet distribution)
72 : \plumedfile
73 : td: TD_GENERALIZED_EXTREME_VALUE LOCATION=-5.0 SCALE=1.0 SHAPE=0.5
74 : \endplumedfile
75 :
76 :
77 : Generalized extreme value distribution with \f$\mu=5.0\f$, \f$\sigma=2.0\f$, and \f$\xi=-0.5\f$ (Weibull distribution)
78 : \plumedfile
79 : td: TD_GENERALIZED_EXTREME_VALUE LOCATION=5.0 SCALE=1.0 SHAPE=-0.5
80 : \endplumedfile
81 :
82 :
83 : The generalized extreme value distribution is only defined for one dimension so for multiple
84 : dimensions we have to use it in combination with the \ref TD_PRODUCT_DISTRIBUTION action as shown in
85 : the following example where we have a Generalized extreme value distribution for argument 1
86 : and uniform distribution for argument 2
87 : \plumedfile
88 : td_gev: TD_GENERALIZED_EXTREME_VALUE LOCATION=-5.0 SCALE=1.0 SHAPE=0.5
89 :
90 : td_uni: TD_UNIFORM
91 :
92 : td_pd: TD_PRODUCT_DISTRIBUTION DISTRIBUTIONS=td_gev,td_uni
93 : \endplumedfile
94 :
95 :
96 : */
97 : //+ENDPLUMEDOC
98 :
99 15 : class TD_GeneralizedExtremeValue: public TargetDistribution {
100 : std::vector<double> center_;
101 : std::vector<double> scale_;
102 : std::vector<double> shape_;
103 : std::vector<double> normalization_;
104 : double GEVdiagonal(const std::vector<double>&, const std::vector<double>&, const std::vector<double>&, const std::vector<double>&, const std::vector<double>&) const;
105 : public:
106 : static void registerKeywords(Keywords&);
107 : explicit TD_GeneralizedExtremeValue(const ActionOptions& ao);
108 : double getValue(const std::vector<double>&) const;
109 : };
110 :
111 :
112 6457 : PLUMED_REGISTER_ACTION(TD_GeneralizedExtremeValue,"TD_GENERALIZED_EXTREME_VALUE")
113 :
114 :
115 6 : void TD_GeneralizedExtremeValue::registerKeywords(Keywords& keys) {
116 6 : TargetDistribution::registerKeywords(keys);
117 24 : keys.add("compulsory","LOCATION","The \\f$\\mu\\f$ parameter of the generalized extreme value distribution.");
118 24 : keys.add("compulsory","SCALE","The \\f$\\sigma\\f$ parameter for the generalized extreme value distribution given as a postive number.");
119 24 : keys.add("compulsory","SHAPE","The \\f$\\xi\\f$ parameter for the generalized extreme value distribution.");
120 12 : keys.use("WELLTEMPERED_FACTOR");
121 12 : keys.use("SHIFT_TO_ZERO");
122 12 : keys.use("NORMALIZE");
123 6 : }
124 :
125 :
126 5 : TD_GeneralizedExtremeValue::TD_GeneralizedExtremeValue(const ActionOptions& ao):
127 : PLUMED_VES_TARGETDISTRIBUTION_INIT(ao),
128 : center_(0),
129 : scale_(0),
130 : shape_(0),
131 5 : normalization_(0)
132 : {
133 10 : parseVector("LOCATION",center_);
134 10 : parseVector("SCALE",scale_);
135 10 : parseVector("SHAPE",shape_);
136 :
137 5 : setDimension(center_.size());
138 5 : 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.");}
139 5 : if(scale_.size()!=getDimension()) {plumed_merror(getName()+": the SCALE keyword does not match the given dimension in MINIMA");}
140 5 : if(shape_.size()!=getDimension()) {plumed_merror(getName()+": the SHAPE keyword does not match the given dimension in MINIMA");}
141 :
142 5 : normalization_.resize(getDimension());
143 15 : for(unsigned int k=0; k<getDimension(); k++) {
144 10 : if(scale_[k]<0.0) {plumed_merror(getName()+": the value given for the scale parameter in SCALE should be larger than 0.0");}
145 5 : normalization_[k] = 1.0/scale_[k];
146 : }
147 5 : checkRead();
148 5 : }
149 :
150 :
151 1605 : double TD_GeneralizedExtremeValue::getValue(const std::vector<double>& argument) const {
152 1605 : return GEVdiagonal(argument,center_,scale_,shape_,normalization_);
153 : }
154 :
155 :
156 1605 : double TD_GeneralizedExtremeValue::GEVdiagonal(const std::vector<double>& argument, const std::vector<double>& center, const std::vector<double>& scale, const std::vector<double>& shape, const std::vector<double>& normalization) const {
157 : double value = 1.0;
158 7215 : for(unsigned int k=0; k<argument.size(); k++) {
159 4815 : double arg=(argument[k]-center[k])/scale[k];
160 : double tx;
161 1605 : if(shape_[k]!=0.0) {
162 2106 : if( shape_[k]>0 && argument[k] <= (center[k]-scale[k]/shape[k]) ) {return 0.0;}
163 1916 : if( shape_[k]<0 && argument[k] > (center[k]-scale[k]/shape[k]) ) {return 0.0;}
164 1134 : tx = pow( (1.0+arg*shape[k]), -1.0/shape[k] );
165 : }
166 : else {
167 201 : tx = exp(-arg);
168 : }
169 2670 : value *= normalization[k] * pow(tx,shape[k]+1.0) * exp(-tx);
170 : }
171 : return value;
172 : }
173 :
174 :
175 :
176 : }
177 4839 : }
|
