Line data Source code
1 : /*******************************************************************************
2 : * *
3 : * Cosmos:(C)oncept et (O)utils (S)tatistique pour les (Mo)deles *
4 : * (S)tochastiques *
5 : * *
6 : * Copyright (C) 2009-2012 LSV & LACL *
7 : * Authors: Paolo Ballarini Benoît Barbot & Hilal Djafri *
8 : * Website: http://www.lsv.ens-cachan.fr/Software/cosmos *
9 : * *
10 : * This program is free software; you can redistribute it and/or modify *
11 : * it under the terms of the GNU General Public License as published by *
12 : * the Free Software Foundation; either version 3 of the License, or *
13 : * (at your option) any later version. *
14 : * *
15 : * This program is distributed in the hope that it will be useful, *
16 : * but WITHOUT ANY WARRANTY; without even the implied warranty of *
17 : * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
18 : * GNU General Public License for more details. *
19 : * *
20 : * You should have received a copy of the GNU General Public License along *
21 : * with this program; if not, write to the Free Software Foundation, Inc., *
22 : * 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. *
23 : * file timeGen.cpp created by Benoit Barbot on 25/01/12. *
24 : *******************************************************************************
25 : */
26 :
27 :
28 :
29 : #include <boost/math/special_functions/prime.hpp>
30 : #include <float.h>
31 : #include <array>
32 : #include <tuple>
33 : #include <limits>
34 : #include <math.h>
35 : #include <random>
36 :
37 : #include "timeGen.hpp"
38 :
39 : /*
40 : Copy from Boost: boost/math/tools/roots.hpp
41 : // (C) Copyright John Maddock 2006.
42 : // Use, modification and distribution are subject to the
43 : // Boost Software License, Version 1.0. (See accompanying file
44 : // LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
45 : */
46 : template <class F>
47 0 : double boost_newton_raphson_iterate(F f, double guess, double min, double max, int max_iter)
48 : {
49 0 : int digits = 16;
50 0 : double f0(0), f1;
51 0 : double result = guess;
52 :
53 0 : double factor = static_cast<double>(ldexp(1.0, 1 - digits));
54 0 : double delta = 1;
55 0 : double delta1 = std::numeric_limits<double>::max();
56 0 : double delta2 = std::numeric_limits<double>::max();
57 :
58 0 : unsigned int count(max_iter);
59 :
60 0 : do{
61 0 : delta2 = delta1;
62 0 : delta1 = delta;
63 0 : std::tie(f0, f1) = f(result);
64 0 : if(0 == f0)
65 0 : break;
66 0 : if(f1 == 0)
67 : {
68 0 : std::cerr << "Zero derivative in newton raphson iterate" << std::endl;
69 0 : exit(EXIT_FAILURE);
70 : }
71 : else
72 : {
73 0 : delta = f0 / f1;
74 : }
75 0 : if(fabs(delta * 2) > fabs(delta2))
76 : {
77 : // last two steps haven't converged, try bisection:
78 0 : delta = (delta > 0) ? (result - min) / 2 : (result - max) / 2;
79 : }
80 0 : guess = result;
81 0 : result -= delta;
82 0 : if(result <= min)
83 : {
84 0 : delta = 0.5F * (guess - min);
85 0 : result = guess - delta;
86 0 : if((result == min) || (result == max))
87 : break;
88 : }
89 0 : else if(result >= max)
90 : {
91 0 : delta = 0.5F * (guess - max);
92 0 : result = guess - delta;
93 0 : if((result == min) || (result == max))
94 : break;
95 : }
96 : // update brackets:
97 0 : if(delta > 0)
98 0 : max = guess;
99 : else
100 0 : min = guess;
101 0 : }while(--count && (fabs(result * factor) < fabs(delta)));
102 :
103 0 : max_iter -= count;
104 :
105 0 : return result;
106 : }
107 :
108 :
109 : using namespace std;
110 :
111 26 : VanDerCorputState::VanDerCorputState(){
112 26 : it = 20;
113 26 : baseId = 0;
114 26 : }
115 :
116 26 : void VanDerCorputState::seed(unsigned long seed){
117 26 : it = seed % 10000223 /* large prime */ ;
118 26 : }
119 :
120 0 : double VanDerCorputState::sample(){
121 0 : long n = it;
122 0 : unsigned int base = boost::math::prime(baseId);
123 0 : baseId = baseId+1;
124 : //(5437*baseId + 3671) % 9973; // A simple random (I choose the parameters randomly) linear congruantial pseudo random number generator.
125 0 : double vdc = 0, denom = 1;
126 0 : while (n){
127 0 : vdc += fmod(n, base) / (denom *= base);
128 0 : n /= base; // note: conversion from 'double' to 'int'
129 : }
130 0 : if(P.verbose > 4)cerr << "Sampling VDC; base: "<< base << "\tn: "<< it << " ->" << vdc << endl;
131 :
132 0 : return vdc;
133 : }
134 :
135 0 : void VanDerCorputState::newTrajectory(){
136 0 : baseId = 0;
137 0 : it++;
138 0 : }
139 :
140 :
141 26 : void timeGen::initRandomGenerator(unsigned int seed){
142 26 : RandomNumber.seed(seed);
143 26 : vanDerCorputSampler.seed(RandomNumber.operator()());
144 26 : kroneckerSampler.seed(seed);
145 26 : }
146 :
147 :
148 26 : KroneckerState::KroneckerState(){
149 26 : it = 0;
150 26 : baseId = 1;
151 26 : phi_d = 1.0 / compute_gamma(P.horizon);
152 26 : if(P.verbose > 4)cerr << "Initializing Kronecker Sampler: horizon: "<< P.horizon << "\tgamma: "<< 1.0/phi_d << endl;
153 :
154 26 : }
155 :
156 26 : double unif01ofint(unsigned long v){
157 26 : double res= 1.0;
158 1660 : while( v != 0){
159 817 : if(v %2){
160 428 : res /= 2.0;
161 : }else{
162 389 : res = (0.5 + res)/2.0;
163 : }
164 817 : v /=2;
165 : }
166 26 : return res;
167 : }
168 :
169 26 : void KroneckerState::seed(unsigned long seed){
170 : //double * tmp = reinterpret_cast<double*>(&seed);
171 26 : seed_val = unif01ofint(seed);
172 :
173 26 : }
174 :
175 0 : double KroneckerState::sample(){
176 0 : double alpha = pow(phi_d,baseId);
177 : //double alpha = sqrt(boost::math::prime(baseId-1));
178 :
179 : double d;
180 0 : double sample = modf( seed_val + (double)it * alpha, &d);
181 0 : if(P.verbose > 4)cerr << "Sampling Kronecker; seed: "<< seed_val << " base: "<< baseId << "\tn: "<< it << "\tbasealpha: " << phi_d << "\talpha: "<< alpha << "->" << sample << endl;
182 0 : baseId = baseId+1;
183 0 : return sample;
184 : }
185 :
186 26 : double KroneckerState::compute_gamma(unsigned long d){
187 26 : double x=1.0;
188 806 : for(int i=0; i< 30; i++){
189 780 : x = x-(pow(x,d+1)-x-1)/((d+1)*pow(x,d)-1);
190 : }
191 26 : return x;
192 : }
193 :
194 0 : void KroneckerState::newTrajectory(){
195 0 : baseId = 1;
196 0 : it++;
197 0 : }
198 :
199 0 : double timeGen::sampleQuasiRandom(size_t){
200 0 : double gentime =0.0;
201 0 : switch (P.randomGen) {
202 : case VDC:
203 0 : gentime = vanDerCorputSampler.sample();
204 0 : break;
205 : case Kronecker:
206 0 : gentime = vanDerCorputSampler.sample();
207 0 : break;
208 : case MT19937:
209 0 : gentime=0.0;
210 : }
211 : //cerr << gentime << "\t";
212 0 : return gentime;
213 : }
214 :
215 216102 : void timeGen::reset(){
216 216102 : switch (P.randomGen) {
217 : case VDC:
218 0 : vanDerCorputSampler.newTrajectory();
219 0 : break;
220 : case Kronecker:
221 0 : kroneckerSampler.newTrajectory();
222 : default: {}
223 : }
224 : //cerr << endl;
225 216102 : }
226 :
227 0 : string timeGen::string_of_dist(DistributionType d,const array<double,PARAM_TBL_SIZE> ¶m, const CustomDistr &cd)const{
228 : //use for debuging
229 0 : switch (d) {
230 : case NORMAL:
231 0 : return "Normal("+ to_string(get<0>(param)) +","+ to_string(get<1>(param))+")";
232 : case GAMMA:
233 0 : return "Gamma("+ to_string(get<0>(param)) +","+ to_string(get<1>(param))+")";
234 : case UNIFORM:
235 0 : return "Uniform("+ to_string(get<0>(param)) +","+ to_string(get<1>(param))+")";
236 : case EXPONENTIAL:
237 0 : return "Exponential("+ to_string(get<0>(param))+")";
238 : case DETERMINISTIC:
239 0 : return "Deterministic("+ to_string(get<0>(param))+")";
240 : case LOGNORMAL:
241 0 : return "LogNormal";
242 : case TRIANGLE:
243 0 : return "Triangle";
244 : case GEOMETRIC:
245 0 : return "Geometric";
246 : case ERLANG:
247 0 : return "Erlang("+ to_string(get<0>(param)) +","+ to_string(get<1>(param))+")";
248 : case PARETO:
249 0 : return "Pareto("+ to_string(get<0>(param)) +","+ to_string(get<1>(param))+")";
250 : case WEIBULL:
251 0 : return "Weibull("+ to_string(get<0>(param)) +","+ to_string(get<1>(param))+")";
252 : case DISCRETEUNIF:
253 0 : return "DiscreteUnif("+ to_string(get<0>(param)) +","+ to_string(get<1>(param))+")";
254 : case MASSACTION:
255 0 : return "MassAction("+ to_string(get<0>(param))+")";
256 : case IMMEDIATE:
257 : case IMDT:
258 0 : return "Immediate";
259 :
260 : // return "Userdefine("+ to_string(param[0]) +","+ to_string(param[1])+","+ to_string(param[2])+","+ to_string(param[3])+")";
261 : case DISCRETEUSERDEFINE:
262 0 : return "DiscreteUserDefine("+ to_string(get<0>(param)) +","+ to_string(get<1>(param))+")";
263 : case USERDEFINE:
264 : case USERDEFINEPOLYNOMIAL: {
265 0 : string s = "UserdefinePolynomial("+ to_string(get<0>(param)) +","+ to_string(get<1>(param))+","+ to_string(get<2>(param))+","+ to_string(get<3>(param))+")";
266 0 : return s + cd.print_poly( (int)get<0>(param));
267 : }
268 : case SYNC:
269 0 : return "Synchronisation";
270 : default:
271 0 : return "Unknown distribution";
272 : }
273 :
274 : }
275 :
276 :
277 : /**
278 : * Call the random generator to generate fire time.
279 : * @param distribution is the type of distribution
280 : * @param param is a vector of parameters of the distribution.
281 : */
282 83919280 : double timeGen::GenerateTime(DistributionType distribution,size_t trid, const array<double,PARAM_TBL_SIZE> ¶m, const CustomDistr &cd) {
283 : //cerr << "sampling " << string_of_dist(distribution,param) << endl;;
284 83919280 : switch (distribution) {
285 : case UNIFORM:
286 : {//UNIF
287 811745 : if(get<0>(param) == get<1>(param)) return get<0>(param);
288 811745 : if(P.randomGen != MT19937){
289 0 : return get<0>(param) + (get<1>(param)-get<0>(param))*sampleQuasiRandom(trid);
290 : }else{
291 811745 : std::uniform_real_distribution<> unif(get<0>(param), get<1>(param));
292 811745 : return unif(RandomNumber);
293 : }
294 : }
295 :
296 : case MASSACTION:
297 : case EXPONENTIAL:
298 : {//EXP
299 : //Exponential distribution is the only marking dependent parameters. Check of validity is required
300 :
301 : //-------------- Rare Event -----------------
302 : //cout << "rate:" << param[0] << endl;
303 71161694 : if(get<0>(param) <= 0) { return 1e200; };
304 : //------------- /Rare Event -----------------
305 :
306 70854499 : if (get<0>(param) <= 0) {
307 0 : cout << "Exponential ditribution should be with rate > 0 not "
308 0 : << get<0>(param) << "\n End of Simulation" << endl;
309 0 : exit(1);
310 : }
311 :
312 70854499 : if(P.randomGen != MT19937){
313 0 : return -log(
314 0 : 1.0 - sampleQuasiRandom(trid)) / get<0>(param);
315 : }else{
316 70854499 : std::exponential_distribution<> exp(get<0>(param));
317 70854499 : return exp(RandomNumber);
318 : }
319 :
320 : }
321 :
322 : case SYNC:
323 : case PLAYER1:
324 : case IMDT:
325 : case IMMEDIATE:
326 0 : return 0;
327 :
328 : case DETERMINISTIC:
329 : {//DETERMINISTIC
330 11945841 : return get<0>(param);
331 : }
332 :
333 : case NORMAL:
334 : {
335 0 : std::normal_distribution<> normal(get<0>(param), get<1>(param));
336 0 : return fmax(0.0,normal(RandomNumber));
337 : }
338 :
339 : case LOGNORMAL:
340 : {//LogNormal
341 0 : double m = get<0>(param); //location
342 0 : double s = get<1>(param); //scale
343 : //double mean = exp(m + s*s/2);
344 : //double variance = (exp(s*s)-1)*exp(2*m+s*s);
345 0 : std::lognormal_distribution<> lognormal(m, s);
346 0 : return lognormal(RandomNumber);
347 : }
348 :
349 : case TRIANGLE:
350 : {//Triangle
351 0 : std::uniform_real_distribution<> unif(0, 1);
352 0 : double p = unif(RandomNumber);
353 0 : double a = get<0>(param); //lower
354 0 : double b = get<2>(param); //upper
355 0 : double c = get<1>(param); //mode
356 0 : double p0 = (c-a)/(b-a); //inflection point
357 0 : if(p < p0)return sqrt((b-a)*(c-a)*p)+a;
358 0 : if(p > p0)return b-sqrt((b-a)*(b-c)*(1.0-p));
359 0 : if(p ==p0)return c; //Prob of this case is zero
360 : }
361 : case GEOMETRIC:
362 : {//GEOMETRIC
363 0 : std::uniform_real_distribution<> unif(0, 1);
364 0 : double p = unif(RandomNumber);
365 0 : return (get<1>(param) * (ceil(log(p) / log(1 - get<0>(param)))));
366 : /*if (p >= param[0]){
367 : return param[1];
368 : } else {
369 : return param[1] * ceil(log(p / param[0]) / log(1 - param[0]) + 1);
370 : }*/
371 : }
372 : case ERLANG:
373 : {//ERLANG
374 0 : std::uniform_real_distribution<> unif(0, 1);
375 0 : double sum = 0.0;
376 0 : for (int i = 0; i < get<0>(param); i++)
377 0 : sum -= log(unif(RandomNumber));
378 0 : return sum / get<1>(param);
379 : }
380 : case PARETO:
381 : {
382 0 : std::uniform_real_distribution<> unif(0, 1);
383 0 : return get<0>(param) / pow(unif(RandomNumber),1/get<1>(param)) ;
384 : }
385 :
386 : case WEIBULL:
387 : {
388 0 : std::weibull_distribution<> weibull(get<0>(param), get<1>(param));
389 0 : return weibull(RandomNumber);
390 : }
391 :
392 :
393 : case GAMMA:
394 : {//GAMMA
395 0 : std::gamma_distribution<> gamma(get<0>(param));
396 0 : return get<1>(param) * gamma(RandomNumber);
397 : }
398 : case DISCRETEUNIF:
399 : {//DISCRETEUNIF
400 0 : std::uniform_int_distribution<> unif((int)get<0>(param), (int)get<1>(param));
401 0 : return unif(RandomNumber);
402 : break;
403 : }
404 :
405 : case DISCRETEUSERDEFINE:
406 : {
407 0 : std::uniform_int_distribution<> unif(0, 100000);
408 0 : unsigned int i= unif(RandomNumber);
409 0 : return cd.userDefineDiscreteDistr(param,i);
410 : }
411 :
412 : case USERDEFINEPOLYNOMIAL:
413 : case USERDEFINE:
414 : {
415 : double gentime;
416 0 : if(P.randomGen != MT19937){
417 0 : gentime = sampleQuasiRandom(trid);
418 : }else{
419 0 : std::uniform_real_distribution<> unif(0.0, 1.0);
420 0 : gentime = unif(RandomNumber);
421 : }
422 :
423 : //cerr << gentime << "\t";
424 :
425 0 : double lower = cd.userDefineLowerBound(param);
426 0 : double upper = cd.userDefineUpperBound(param);
427 :
428 :
429 : /*//Isotrop Hack
430 : std::uniform_real_distribution<> unif(lower, upper);
431 : gentime = unif(RandomNumber);
432 : //gentime = sampleQuasiRandom(trid);
433 : //gentime = (upper-lower)*gentime + lower;
434 : cerr << gentime << "\t";
435 : //cerr << cd.userDefineCDF(param,gentime) << "\t";
436 : return gentime;
437 : //Isotrop Hack
438 : */
439 :
440 : //cerr << "sample(" << gentime << ",[" << lower << "," << upper << "]):" <<endl;
441 0 : double initialpt = (lower+upper)/ 2.0;
442 0 : double pt = boost_newton_raphson_iterate([&](double x){
443 0 : const auto cdf = cd.userDefineCDF(param,x);
444 0 : const auto pdf = cd.userDefinePDF(param,x);
445 : // cerr << "it:" << x << endl;
446 0 : return make_tuple(cdf-gentime, pdf);
447 0 : }, initialpt, lower, upper, 100);
448 : //cerr << pt << "\t";
449 0 : return pt;
450 : break;
451 :
452 : }
453 :
454 :
455 : }
456 0 : return DBL_MIN;
457 87 : }
|
