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