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 SimulatorBoundedRE.cpp created by Benoit Barbot on 06/12/11. *
24 : *******************************************************************************
25 : */
26 :
27 : #include <list>
28 : #include <sys/resource.h>
29 :
30 :
31 : #include "SimulatorBoundedRE.hpp"
32 :
33 : #include "numSolverBB.hpp"
34 : #include "numSolverSH.hpp"
35 :
36 : using namespace std;
37 :
38 : template<class S,class DEDS>
39 0 : SimulatorBoundedREBase<S,DEDS>::SimulatorBoundedREBase(DEDS& N,LHA_orig<typename DEDS::stateType>& A,int m):SimulatorREBase<S,DEDS>(N,A){
40 0 : switch (m) {
41 : case 1:
42 0 : numSolv = new numericalSolver();
43 0 : break;
44 :
45 : case 2:
46 0 : numSolv = new numSolverBB();
47 0 : break;
48 :
49 : case 3:
50 0 : numSolv = new numSolverSH();
51 0 : break;
52 : }
53 0 : this->muprob = numSolv;
54 0 : delete this->EQ;
55 0 : this->N.initialize( this->muprob);
56 : //numSolv->initVect(T);
57 0 : }
58 :
59 : template<class S,class DEDS>
60 0 : void SimulatorBoundedREBase<S,DEDS>::initVect(int T){
61 0 : lambda = numSolv->uniformizeMatrix();
62 0 : cerr << "lambda:" << lambda<< endl;
63 0 : numSolv->initVect(T);
64 0 : }
65 :
66 : template<class S,class DEDS>
67 0 : BatchR SimulatorBoundedREBase<S,DEDS>::RunBatch(){
68 :
69 : //cerr << "test(";
70 0 : numSolv->reset();
71 : //cerr << ")" << endl;
72 :
73 0 : BatchR batchResult(1,0);
74 :
75 0 : list<simulationState<DEDS, EventsQueue> > statevect(this->BatchSize);
76 : //delete EQ;
77 :
78 0 : if(P.verbose>=1){
79 0 : cerr << "Initial round:";
80 0 : numSolv->printState();
81 0 : cerr << "\tremaining trajectories: "<< statevect.size() << "\tInit Prob:";
82 0 : cerr << numSolv->getVect()[0] << endl;
83 : }
84 0 : for (auto &it : statevect) {
85 0 : this->N.Origine_Rate_Table = vector<double>(this->N.tr,0.0);
86 0 : this->N.Rate_Table = vector<double>(this->N.tr,0.0);
87 0 : this->EQ = new EventsQueue(this->N.getNbTransition());
88 0 : this->reset();
89 :
90 0 : this->N.initialEventsQueue(*(this->EQ),*this);
91 :
92 : //AE = A.GetEnabled_A_Edges( N.Marking);
93 :
94 0 : it.saveState(&(this->N),&(this->A),&(this->EQ));
95 : }
96 :
97 : //cout << "new batch" << endl;
98 0 : while (!statevect.empty()) {
99 0 : numSolv->stepVect();
100 0 : if(P.verbose>=1){
101 0 : cerr << "new round:";
102 0 : numSolv->printState();
103 0 : cerr << "\tremaining trajectories: "<< statevect.size() << "\tInit Prob:";
104 0 : cerr << numSolv->getVect()[0] << endl;
105 : }
106 :
107 0 : for (auto it= statevect.begin(); it != statevect.end() ; it++) {
108 :
109 0 : it->loadState(&(this->N),&(this->A),&(this->EQ));
110 :
111 :
112 : //cerr << A.Likelihood << endl;
113 : //cerr << "mu:\t" << mu() << " ->\t";
114 0 : bool continueb = this->SimulateOneStep();
115 : //cerr << mu() << endl;
116 0 : if(numSolv->getVect().size() <= 1){
117 0 : continueb=false;
118 0 : this->Result.accept=false;
119 : }
120 :
121 0 : if((!this->EQ->isEmpty()) && continueb) {
122 0 : it->saveState(&(this->N),&(this->A),&(this->EQ));
123 : } else {
124 0 : batchResult.addSim(this->Result);
125 0 : delete this->EQ;
126 0 : it = statevect.erase(it);
127 0 : it--;
128 :
129 : //log the result
130 0 : if (this->Result.accept && this->logResult){
131 0 : for(size_t i=0; i<this->Result.quantR.size();i++){
132 0 : if (i>0)this->logvalue << "\t";
133 0 : this->logvalue << this->Result.quantR[i];
134 : }
135 0 : this->logvalue << endl;
136 : }
137 : }
138 :
139 :
140 :
141 : }
142 : }
143 :
144 : //cerr << "test" << endl;
145 : //exit(0);
146 :
147 : rusage ruse;
148 0 : getrusage(RUSAGE_SELF, &ruse);
149 0 : cerr <<endl << endl << "Total Time: "<< ruse.ru_utime.tv_sec + ruse.ru_utime.tv_usec / 1000000.
150 0 : << "\tTotal Memory: " << ruse.ru_maxrss << "ko" << endl << endl;
151 :
152 0 : return (batchResult);
153 : }
154 :
155 :
156 : template <class S>
157 0 : double SPNBaseBoundedRE<S>::mu(){
158 :
159 0 : vector<int> vect (this->muprob->S.begin()->first->size(),0);
160 :
161 0 : static_cast< S* >(this)->lumpingFun(this->Marking,vect);
162 0 : int stateN = this->muprob->findHash(&vect);
163 :
164 0 : if(stateN<0){
165 : //cerr << numSolv->getVect()<< endl
166 0 : cerr << "statevect(";
167 0 : for(size_t i =0 ; i<vect.size() ; i++)cerr << vect[i]<< ",";
168 0 : cerr << ")" << endl<<"state not found" << endl;
169 0 : this->print_state(vect);
170 0 : return 0.0;
171 : //exit(EXIT_FAILURE);
172 : }
173 :
174 0 : return(this->muprob->getMu(stateN));
175 : }
176 :
177 :
178 : template <class S>
179 0 : void SPNBaseBoundedRE<S>::update(double ctime,size_t,const abstractBinding&,EventsQueue &EQ, timeGen &TG){
180 0 : Event F;
181 : //check if the current transition is still enabled
182 :
183 0 : this->Rate_Sum = 0;
184 0 : this->Origine_Rate_Sum = 0;
185 :
186 : //Run over all transition
187 0 : for (size_t it = 0; it < SPN::tr-2; it++) {
188 0 : for(auto bindex = this->Transition[it].bindingList.begin() ;
189 0 : bindex != this->Transition[it].bindingList.end() ; ++bindex){
190 0 : if(this->IsEnabled(it, *bindex)){
191 0 : if (EQ.isScheduled(it, bindex->id())) {
192 0 : generateEvent(ctime,F, it ,*bindex, TG, *(static_cast<S *>(this)) );
193 0 : EQ.replace(F);
194 : } else {
195 0 : generateEvent(ctime,F, it ,*bindex, TG, *(static_cast<S *>(this)) );
196 0 : EQ.insert(F);
197 : }
198 : }else{
199 0 : if(EQ.isScheduled(it, bindex->id()))
200 0 : EQ.remove(it,bindex->id());
201 : }
202 : }
203 : }
204 :
205 0 : abstractBinding bpuit;
206 0 : generateEvent(ctime,F, (SPN::tr-2),bpuit, TG, *(static_cast<S *>(this)));
207 0 : if(! this->doubleIS_mode){
208 0 : EQ.replace(F);
209 : }
210 :
211 0 : generateEvent(ctime,F, (SPN::tr-1),bpuit, TG, *(static_cast<S *>(this)));
212 0 : if(! this->doubleIS_mode){
213 0 : EQ.replace(F);
214 : }
215 :
216 0 : };
217 :
218 :
219 : template <class S>
220 0 : void SPNBaseBoundedRE<S>::getParams(size_t Id,const abstractBinding& b){
221 :
222 0 : static_cast<S*>(this)->GetDistParameters(Id,b);
223 0 : double origin_rate = this->ParamDistr[0];
224 0 : if(Id== SPN::tr-2){
225 0 : origin_rate = this->muprob->maxRate - this->Origine_Rate_Sum;
226 : //cerr << "lambda:\t" << lambda << "\tselfloop:\t" << origin_rate << endl;
227 : }
228 0 : this->ParamDistr[0]= static_cast<S*>(this)->ComputeDistr( Id, b,origin_rate);
229 0 : this->ParamDistr[1]= origin_rate;
230 0 : }
231 :
232 :
233 : template <class S>
234 0 : double SPNBaseBoundedRE<S>::ComputeDistr(size_t t ,const abstractBinding& b, double origin_rate ){
235 :
236 : //cerr << endl<< "mux" << endl;
237 0 : double mux = mu();
238 :
239 0 : if(t== SPN::tr-1){
240 0 : if (mux==0.0)return 1E200;
241 :
242 0 : if(P.verbose>3){
243 0 : this->Marking.printHeader(cerr);cerr << endl;
244 0 : this->Marking.print(cerr,0.0);cerr << endl;
245 0 : vector<int> vect (this->muprob->S.begin()->first->size(),0);
246 0 : static_cast<S*>(this)->lumpingFun(this->Marking,vect);
247 0 : static_cast<S*>(this)->print_state(vect);
248 : }
249 :
250 0 : if(this->Origine_Rate_Sum >= this->Rate_Sum){
251 0 : if(P.verbose>3 )cerr << "trans:sink distr: "<< this->Origine_Rate_Sum - this->Rate_Sum << " origine_rate:" << this->Origine_Rate_Sum <<" Rate: " << this->Rate_Sum << endl;
252 : //cerr << "strange !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!" << endl;
253 0 : return( (this->Origine_Rate_Sum - this->Rate_Sum) );
254 : }else{
255 0 : if(P.verbose>3 && (this->Origine_Rate_Sum < 0.99* this->Rate_Sum)){
256 0 : cerr << "Reduce model does not guarantee variance" << endl;
257 0 : cerr << "Initial sum of rate: " << this->Origine_Rate_Sum << " Reduce one: " << this->Rate_Sum << " difference: " << this->Origine_Rate_Sum - this->Rate_Sum << endl ;
258 : //exit(EXIT_FAILURE);
259 : }
260 : //cerr << "trans:sink distr: 0 " << endl;
261 0 : return 0.0 ;};
262 : };
263 0 : if( mux==0.0 || mux==1.0) return(origin_rate);
264 :
265 : double distr;
266 0 : SPN::fire(t,b,0.0);
267 0 : static_cast<numericalSolver*>(this->muprob)->stepVect();
268 0 : distr = origin_rate *( mu() / mux);
269 0 : if(P.verbose>3 )cerr << "trans: " << this->Transition[t].label << "\tdistr: "<< distr << "\torigin Rate: "<< origin_rate << "\tmu: " << mu()<< "\tmu prec: " << mux << endl;
270 :
271 0 : static_cast<numericalSolver*>(this->muprob)->previousVect();
272 0 : SPN::unfire(t,b);
273 0 : return(distr);
274 : }
|
