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 SimulatorContinuousBounded.cpp created by Benoit Barbot on 31/01/12. *
24 : *******************************************************************************
25 : */
26 :
27 :
28 : #include <iostream>
29 : #include "SimulatorContinuousBounded.hpp"
30 : #include "float.h"
31 : #include "math.h"
32 : #include <sys/resource.h>
33 : #include <algorithm>
34 : //#include <boost/math/distributions/normal.hpp>
35 :
36 : #include <boost/math/distributions/normal.hpp>
37 : #include <boost/math/distributions/binomial.hpp>
38 : #include <limits>
39 :
40 : template<class DEDS>
41 0 : SimulatorContinuousBounded<DEDS>::SimulatorContinuousBounded(DEDS &N,LHA_orig<typename DEDS::stateType>& A,int m,double e,int js):SimulatorBoundedREBase<SimulatorContinuousBounded,DEDS>(N,A,m){
42 0 : epsilon = e;
43 0 : if(js>0){
44 0 : jumpsize = js;
45 0 : singleIS =false;
46 : }else{
47 0 : jumpsize = -js;
48 0 : singleIS = true;
49 : }
50 : //boost::math::normal norm;
51 : //Normal_quantile = quantile(norm, 0.5 + 0.99 / 2.0);
52 0 : }
53 :
54 : template<class DEDS>
55 0 : void SimulatorContinuousBounded<DEDS>::initVectCo(double t){
56 :
57 0 : this->lambda = this->numSolv->uniformizeMatrix();
58 0 : cerr << "lambda:" << this->lambda<< endl;
59 0 : fg = unique_ptr<FoxGlynn>(new FoxGlynn((this->lambda * t), DBL_EPSILON , 1/DBL_EPSILON ,epsilon));
60 0 : if (fg->isValid){
61 0 : cerr << "fox_glyn:" << fg->left << "," << fg->right << " Total weigts:"<< fg->total_weight<< endl;
62 : /*for(int i = 0; i<= fg->right - fg->left; i++){
63 : cerr << fg->left+i << " " << fg->weights[i]/ fg->total_weight << endl;
64 : }*/
65 : }
66 : //fg->right = 10;
67 0 : double lambda2= this->lambda;
68 0 : this->initVect(fg->right+1);
69 0 : this->lambda = lambda2;
70 0 : }
71 :
72 : template<class DEDS>
73 0 : BatchR SimulatorContinuousBounded<DEDS>::RunBatch(){
74 : //cerr << "test(";
75 0 : this->numSolv->reset();
76 :
77 : //cerr << ")" << endl;
78 0 : int left = (int)max(this->numSolv->getMinT(),fg->left);
79 0 : int Nmax = (int)ceil(((double)fg->right - left)/jumpsize) ;
80 : //cerr << "minT:\t" << numSolv->getMinT() << endl;
81 0 : cerr << "Number of batch:\t" << Nmax+1 << "\tleft: " << left << "\tright: " << fg->right << "\tjumpsize:" << jumpsize << endl;
82 : //cerr << "left:\t" << left << endl;
83 :
84 0 : int n =-1;
85 :
86 :
87 0 : BatchR batchResult(2*(Nmax+1),0);
88 :
89 :
90 : //Copy and sum Poisson Weight
91 0 : for(int i =0; i<= Nmax; i++)
92 0 : batchResult.Min[2*i] = 0.0;
93 0 : for(int i=left; i<= fg->right; i++){
94 0 : int j = (int)ceil((double)( i - left)/jumpsize);
95 0 : batchResult.Min[2*j] += fg->weights[i - fg->left];
96 : }
97 0 : for(int i =0; i<= Nmax; i++)
98 0 : batchResult.Min[2*i] /= fg->total_weight;
99 :
100 0 : list<simulationState<DEDS, EventsQueue >> statevect((Nmax+1)*this->BatchSize);
101 :
102 0 : int c =0;
103 0 : if(P.verbose>=1)cerr << "new round:"<< n << "\tremaining trajectories: "<< statevect.size() << endl;
104 0 : for (auto it= statevect.begin(); it != statevect.end() ; it++) {
105 0 : this->N.Origine_Rate_Table = vector<double>(this->N.tr,0.0);
106 0 : this->N.Rate_Table = vector<double>(this->N.tr,0.0);
107 :
108 0 : bool continueb = false;
109 : //we try to find a trajectory reaching the precondition.
110 0 : while(!continueb){
111 : //we build a new trajectory from the initial state.
112 0 : continueb =true;
113 0 : this->EQ = new EventsQueue(this->N.getNbTransition());
114 0 : this->reset();
115 : //We simulate until either the condition is satisfied or the trajectory reach a deadend.
116 0 : while(!this->N.precondition(this->N.getState()) && continueb){
117 0 : continueb = this->SimulateOneStep();
118 : }
119 : }
120 :
121 0 : if (singleIS) {
122 0 : it->maxStep = fg->right;
123 : }else{
124 0 : it->maxStep = (int)fmin(fg->right , left + (Nmax - (c/this->BatchSize))*jumpsize);
125 : }
126 :
127 : //cerr << "new path:\t" << it->maxStep << endl;
128 :
129 0 : c++;
130 :
131 0 : it->saveState(&(this->N),&(this->A),&(this->EQ));
132 : }
133 :
134 : //cout << "new batch" << endl;
135 0 : while (!statevect.empty()) {
136 0 : this->numSolv->stepVect();
137 0 : if(P.verbose>=3)cerr << this->numSolv->getVect() << endl;
138 0 : n++;
139 0 : if(P.verbose>=1)cerr << "new round:"<< n << "\tremaining trajectories: "<< statevect.size() << endl;
140 :
141 0 : for (auto it= statevect.begin(); it != statevect.end() ; it++) {
142 0 : if(it->maxStep >= fg->right -n){
143 :
144 : //cerr << "vect:\t" << it->maxStep;
145 0 : it->loadState(&(this->N),&(this->A),&(this->EQ));
146 :
147 :
148 : //cerr << A.Likelihood << endl;
149 :
150 : //cerr << "mu:\t" << mu() << " -> ";
151 :
152 0 : if (it->maxStep == fg->right -n) {
153 : //We first need to initialise the trajectory
154 0 : this->N.InitialEventsQueue(*(this->EQ),*this);
155 0 : if(P.verbose>=2)
156 : //cerr << "new Path: " << it->maxStep << "\tmuinit: " << mu() << endl;
157 0 : it->saveState(&(this->N),&(this->A),&(this->EQ));
158 : } else {
159 :
160 0 : bool continueb = this->SimulateOneStep();
161 : //cerr << "\t" << mu() << endl;
162 :
163 0 : if((!this->EQ->isEmpty()) && continueb) {
164 0 : it->saveState(&(this->N),&(this->A),&(this->EQ));
165 : } else {
166 0 : int i = (int)ceil((double)(it->maxStep- left)/jumpsize) ;
167 0 : int i2 =(int)fmax(0.0,ceil((double)(n - left)/jumpsize));
168 : //cerr << "Successfull trajectory: "<< it->maxStep << " -> " << i << endl;
169 0 : if (this->Result.accept ) {
170 :
171 0 : if (this->Result.quantR[0] * (1 - this->Result.quantR[0]) != 0) batchResult.IsBernoulli[0] = false;
172 :
173 0 : batchResult.Isucc+=1;
174 0 : if(singleIS)for(int j = i2; j<= Nmax; j++){
175 0 : batchResult.Mean[2*j]+=1;
176 0 : batchResult.Mean[2*j+1] += this->Result.quantR[0];
177 0 : batchResult.M2[2*j+1] += pow(this->Result.quantR[0] , 2);
178 0 : batchResult.M3[2*j+1] += pow(this->Result.quantR[0] , 3);
179 0 : batchResult.M4[2*j+1] += pow(this->Result.quantR[0] , 4);
180 0 : batchResult.Min[2*j+1] = fmin(batchResult.Min[2*j+1],this->Result.quantR[0]);
181 0 : batchResult.Max[2*j+1] = fmax(batchResult.Max[2*j+1],this->Result.quantR[0]);
182 : } else {
183 0 : batchResult.Mean[2*i]+=1;
184 0 : batchResult.Mean[2*i+1] += this->Result.quantR[0];
185 0 : batchResult.M2[2*i+1] += pow(this->Result.quantR[0] , 2);
186 0 : batchResult.M3[2*i+1] += pow(this->Result.quantR[0] , 3);
187 0 : batchResult.M4[2*i+1] += pow(this->Result.quantR[0] , 4);
188 0 : batchResult.Min[2*i+1] = fmin(batchResult.Min[2*i+1],this->Result.quantR[0]);
189 0 : batchResult.Max[2*i+1] = fmax(batchResult.Max[2*i+1],this->Result.quantR[0]);
190 : }
191 :
192 : }
193 0 : if(singleIS){for(int j = 0; j<= Nmax; j++)batchResult.M2[2*j]+=1;}
194 0 : else batchResult.M2[2*i]+=1;
195 0 : batchResult.I++;
196 :
197 :
198 0 : delete this->EQ;
199 0 : it = statevect.erase(it);
200 :
201 0 : it--;
202 : }
203 : }
204 :
205 : }
206 :
207 : }
208 :
209 : }
210 :
211 0 : return (batchResult);
212 : }
213 :
|
