1 /* ----------------------------------------------------------------------
2     This is the
3 
4     ██╗     ██╗ ██████╗  ██████╗  ██████╗ ██╗  ██╗████████╗███████╗
5     ██║     ██║██╔════╝ ██╔════╝ ██╔════╝ ██║  ██║╚══██╔══╝██╔════╝
6     ██║     ██║██║  ███╗██║  ███╗██║  ███╗███████║   ██║   ███████╗
7     ██║     ██║██║   ██║██║   ██║██║   ██║██╔══██║   ██║   ╚════██║
8     ███████╗██║╚██████╔╝╚██████╔╝╚██████╔╝██║  ██║   ██║   ███████║
9     ╚══════╝╚═╝ ╚═════╝  ╚═════╝  ╚═════╝ ╚═╝  ╚═╝   ╚═╝   ╚══════╝®
10 
11     DEM simulation engine, released by
12     DCS Computing Gmbh, Linz, Austria
13     http://www.dcs-computing.com, office@dcs-computing.com
14 
15     LIGGGHTS® is part of CFDEM®project:
16     http://www.liggghts.com | http://www.cfdem.com
17 
18     Core developer and main author:
19     Christoph Kloss, christoph.kloss@dcs-computing.com
20 
21     LIGGGHTS® is open-source, distributed under the terms of the GNU Public
22     License, version 2 or later. It is distributed in the hope that it will
23     be useful, but WITHOUT ANY WARRANTY; without even the implied warranty
24     of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. You should have
25     received a copy of the GNU General Public License along with LIGGGHTS®.
26     If not, see http://www.gnu.org/licenses . See also top-level README
27     and LICENSE files.
28 
29     LIGGGHTS® and CFDEM® are registered trade marks of DCS Computing GmbH,
30     the producer of the LIGGGHTS® software and the CFDEM®coupling software
31     See http://www.cfdem.com/terms-trademark-policy for details.
32 
33 -------------------------------------------------------------------------
34     Contributing author and copyright for this file:
35     (if not contributing author is listed, this file has been contributed
36     by the core developer)
37 
38     Copyright 2012-     DCS Computing GmbH, Linz
39     Copyright 2009-2012 JKU Linz
40 ------------------------------------------------------------------------- */
41 
42 #ifndef LMP_PROBABILITY_DISTRIBUTION_H
43 #define LMP_PROBABILITY_DISTRIBUTION_H
44 
45 #include <cmath>
46 #include <stdio.h>
47 #include <string.h>
48 #include "random_park.h"
49 #include "error.h"
50 #include "pointers.h"
51 
52 enum{RANDOM_CONSTANT,RANDOM_UNIFORM,RANDOM_GAUSSIAN,RANDOM_LOGNORMAL};
53 
54 namespace LMP_PROBABILITY_NS {
55   class PDF
56   {
57     public:
PDF(LAMMPS_NS::Error * error)58       PDF(LAMMPS_NS::Error *error)
59       {
60           mu_ = sigma_ = min_ = max_ = 0.;
61           h1_ = h2_ = 0.;
62           this->error = error;
63       }
~PDF()64       ~PDF(){}
65 
66       int rand_style_;
67 
68       double mu_,sigma_;
69       double min_,max_;
70 
71       // helper
72       double h1_,h2_;
73 
74       LAMMPS_NS::Error *error;
75 
rand_style()76       inline int rand_style()
77       { return rand_style_; }
78 
set_min_max(double min,double max)79       inline void set_min_max(double min,double max)
80       {
81           min_ = min;
82           max_ = max;
83       }
84 
set_params(double)85       template<int RAND_STYLE> void set_params(double)
86       { error->all(FLERR,"Faulty usage of Probability::set_params"); }
87 
set_params(double,double)88       template<int RAND_STYLE> void set_params(double,double)
89       { error->all(FLERR,"Faulty usage of Probability::set_params"); }
90   };
91 
pdf_max(PDF * pdf)92   inline double pdf_max(PDF *pdf)
93   {
94       return pdf->max_;
95   }
96 
pdf_min(PDF * pdf)97   inline double pdf_min(PDF *pdf)
98   {
99       return pdf->min_;
100   }
101 
expectancy_value(PDF * pdf)102   template <int RAND_STYLE> inline double expectancy_value(PDF *pdf)
103   {
104       pdf->error->all(FLERR,"Faulty usage of Probability::expectancy");
105       return 0.;
106   }
107 
cubic_expectancy_value(PDF * pdf)108   template <int RAND_STYLE> inline double cubic_expectancy_value(PDF *pdf)
109   {
110       pdf->error->all(FLERR,"Faulty usage of Probability::volume_expectancy");
111       return 0.;
112   }
113 
rand_value(PDF * pdf,LAMMPS_NS::RanPark * rp)114   template <int RAND_STYLE> inline double rand_value(PDF *pdf,LAMMPS_NS::RanPark *rp)
115   {
116       pdf->error->all(FLERR,"Faulty usage of Probability::rand");
117       return 0.;
118   }
119 
120   //------------------------------------------------------------------------------
121   // CONSTANT
122   //------------------------------------------------------------------------------
123 
124   template<> inline void PDF::set_params<RANDOM_CONSTANT>(double val)
125   {
126       rand_style_ = RANDOM_CONSTANT;
127       mu_ = val;
128       set_min_max(mu_,mu_);
129   }
130 
131   template<> inline double cubic_expectancy_value<RANDOM_CONSTANT>(PDF *pdf)
132   {
133 
134       return pdf->mu_*pdf->mu_*pdf->mu_;
135   }
136 
137   template<> inline double expectancy_value<RANDOM_CONSTANT>(PDF *pdf)
138   {
139 
140       return pdf->mu_;
141   }
142 
143   template<> inline double rand_value<RANDOM_CONSTANT>(PDF *pdf,LAMMPS_NS::RanPark *rp)
144   {
145       return pdf->mu_;
146   }
147 
148   //------------------------------------------------------------------------------
149   // UNIFORM
150   //------------------------------------------------------------------------------
151 
152   template<> inline void PDF::set_params<RANDOM_UNIFORM>(double min, double max)
153   {
154       rand_style_ = RANDOM_UNIFORM;
155       set_min_max(min,max);
156       h1_ = 2./(1./(min_*min_)-1./(max_*max_));
157       h2_ = h1_/(2.*min_*min_);
158   }
159 
160   template<> inline double cubic_expectancy_value<RANDOM_UNIFORM>(PDF *pdf)
161   {
162       return 0.25*(pdf->max_*pdf->max_*pdf->max_+
163                    pdf->max_*pdf->max_*pdf->min_+
164                    pdf->max_*pdf->min_*pdf->min_+
165                    pdf->min_*pdf->min_*pdf->min_);
166   }
167 
168   template<> inline double expectancy_value<RANDOM_UNIFORM>(PDF *pdf)
169   {
170       return sqrt(pdf->h1_/(2.*(pdf->h2_-0.5)));
171   }
172 
173   template<> inline double rand_value<RANDOM_UNIFORM>(PDF *pdf,LAMMPS_NS::RanPark *rp)
174   {
175       double rn =  rp->uniform();
176       return sqrt(pdf->h1_/(2.*(pdf->h2_-rn)));
177   }
178 
179   //------------------------------------------------------------------------------
180   // GAUSSIAN
181   //------------------------------------------------------------------------------
182 
183   template<> inline void PDF::set_params<RANDOM_GAUSSIAN>(double mu, double sigma)
184   {
185       rand_style_ = RANDOM_GAUSSIAN;
186       mu_ = mu;
187       sigma_ = sigma;
188 
189       // set min-max to +- 3 sigma (99.73% of all values)
190       set_min_max(mu_-3.*sigma_, mu_+3.*sigma_);
191 
192       if(min_ < 0.)
193          error->all(FLERR,"Probablity distribution: mu-3*sigma < 0, please increase mu or decrease sigma");
194   }
195 
196   template<> inline double cubic_expectancy_value<RANDOM_GAUSSIAN>(PDF *pdf)
197   {
198       return pdf->mu_*(pdf->mu_*pdf->mu_+3*pdf->sigma_*pdf->sigma_);
199   }
200 
201   template<> inline double expectancy_value<RANDOM_GAUSSIAN>(PDF *pdf)
202   {
203       return pdf->mu_;
204   }
205 
206   template<> inline double rand_value<RANDOM_GAUSSIAN>(PDF *pdf,LAMMPS_NS::RanPark *rp)
207   {
208       double value;
209       do
210       {
211           value = pdf->mu_ + rp->gaussian() * pdf->sigma_;
212       } while (value < pdf->min_ || value > pdf->max_);
213       return value;
214   }
215 
216   //------------------------------------------------------------------------------
217   // LOGNORMAL
218   //------------------------------------------------------------------------------
219 
220   template<> inline void PDF::set_params<RANDOM_LOGNORMAL>(double mu, double sigma)
221   {
222       error->all(FLERR,"lognormal distribution currently deactivated");
223 
224       rand_style_ = RANDOM_LOGNORMAL;
225       mu_ = mu;
226       sigma_ = sigma;
227 
228       // also here, take +- 3 sigma as min/max
229       // change in expectancy considered negligable
230       double min =  exp(mu_ - 3. * sigma_);
231       double max =  exp(mu_ + 3. * sigma_);
232       set_min_max(min, max);
233       if(min_ < 0.)
234             error->all(FLERR,"Probablity distribution: exp(mu-3*sigma) < 0, please increase mu or decrease sigma");
235   }
236 
237   template<> inline double cubic_expectancy_value<RANDOM_LOGNORMAL>(PDF *pdf)
238   {
239       return exp(3.*pdf->mu_+4.5*pdf->sigma_*pdf->sigma_);
240   }
241 
242   template<> inline double expectancy_value<RANDOM_LOGNORMAL>(PDF *pdf)
243   {
244       return exp(pdf->mu_ + 0.5 * pdf->sigma_ * pdf->sigma_);
245   }
246 
247   template<> inline double rand_value<RANDOM_LOGNORMAL>(PDF *pdf,LAMMPS_NS::RanPark *rp)
248   {
249      double value;
250      do
251      {
252          value = exp(pdf->mu_ + rp->gaussian() * pdf->sigma_);
253      } while (value < pdf->min_ || value > pdf->max_);
254      return value;
255   }
256 
257   //------------------------------------------------------------------------------
258   // MASTER FUNCTIONS
259   //------------------------------------------------------------------------------
260 
expectancy(PDF * pdf)261   inline double expectancy(PDF *pdf)
262   {
263       if(pdf->rand_style_ == RANDOM_CONSTANT) return expectancy_value<RANDOM_CONSTANT>(pdf);
264       else if(pdf->rand_style_ == RANDOM_UNIFORM) return expectancy_value<RANDOM_UNIFORM>(pdf);
265       else if(pdf->rand_style_ == RANDOM_GAUSSIAN) return expectancy_value<RANDOM_GAUSSIAN>(pdf);
266       else if(pdf->rand_style_ == RANDOM_LOGNORMAL) return expectancy_value<RANDOM_LOGNORMAL>(pdf);
267       else pdf->error->all(FLERR,"Faulty implemantation in Probability::expectancy");
268       return 0.;
269   }
270 
cubic_expectancy(PDF * pdf)271   inline double cubic_expectancy(PDF *pdf)
272   {
273       if(pdf->rand_style_ == RANDOM_CONSTANT) return cubic_expectancy_value<RANDOM_CONSTANT>(pdf);
274       else if(pdf->rand_style_ == RANDOM_UNIFORM) return cubic_expectancy_value<RANDOM_UNIFORM>(pdf);
275       else if(pdf->rand_style_ == RANDOM_GAUSSIAN) return cubic_expectancy_value<RANDOM_GAUSSIAN>(pdf);
276       else if(pdf->rand_style_ == RANDOM_LOGNORMAL) return cubic_expectancy_value<RANDOM_LOGNORMAL>(pdf);
277       else pdf->error->all(FLERR,"Faulty implemantation in Probability::expectancy");
278       return 0.;
279   }
280 
rand(PDF * pdf,LAMMPS_NS::RanPark * rp)281   inline double rand(PDF *pdf,LAMMPS_NS::RanPark *rp)
282   {
283       if(pdf->rand_style_ == RANDOM_CONSTANT) return rand_value<RANDOM_CONSTANT>(pdf,rp);
284       else if(pdf->rand_style_ == RANDOM_UNIFORM) return rand_value<RANDOM_UNIFORM>(pdf,rp);
285       else if(pdf->rand_style_ == RANDOM_GAUSSIAN) return rand_value<RANDOM_GAUSSIAN>(pdf,rp);
286       else if(pdf->rand_style_ == RANDOM_LOGNORMAL) return rand_value<RANDOM_LOGNORMAL>(pdf,rp);
287       else pdf->error->all(FLERR,"Faulty implemantation in Probability::rand");
288       return 0.;
289   }
290 
291 };
292 
293 #endif
294