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 SimulatorRE.cpp created by Benoit Barbot on 09/11/11. *
24 : *******************************************************************************
25 : */
26 :
27 : #include "EventsQueue.hpp"
28 : #include <iostream>
29 : #include <set>
30 : #include <math.h>
31 : #include <float.h>
32 : #include <time.h>
33 : #include "marking.hpp"
34 : //#include "RareEvent.hpp"
35 :
36 : #include "SimulatorRE.hpp"
37 :
38 : using namespace std;
39 :
40 : template<class DEDS>
41 2457038 : void generateEvent(double ctime,Event& E,size_t Id,const abstractBinding& b, timeGen &TG, DEDS & N){
42 2457038 : double t = ctime;
43 2457038 : if (N.Transition[Id].DistTypeIndex != IMMEDIATE) {
44 2457038 : N.getParams(Id, b);
45 2457038 : t += TG.GenerateTime(N.Transition[Id].DistTypeIndex, Id, N.ParamDistr, N.customDistr);
46 :
47 2457038 : N.Rate_Table[Id] = N.ParamDistr[0];
48 2457038 : N.Origine_Rate_Table[Id] = N.ParamDistr[1];
49 2457038 : N.Rate_Sum = N.Rate_Sum + N.ParamDistr[0];
50 2457038 : N.Origine_Rate_Sum = N.Origine_Rate_Sum + N.ParamDistr[1];
51 :
52 : }
53 2457038 : double w=0.0;
54 2457038 : if (N.Transition[Id].DistTypeIndex > 2) {
55 2457038 : N.ParamDistr[0]= N.GetWeight(Id,b);
56 2457038 : w = TG.GenerateTime(EXPONENTIAL, Id, N.ParamDistr, N.customDistr);
57 : //vector<double> wParam(1, N.GetWeight(Id));
58 : //w = GenerateTime(2, wParam);
59 : }
60 :
61 2457038 : E.transition = Id;
62 2457038 : E.time = t;
63 2457038 : E.priority = N.GetPriority(Id,b);
64 2457038 : E.weight = w;
65 2457038 : E.binding = b;
66 :
67 2457038 : }
68 :
69 :
70 : template <class S, class DEDS>
71 1 : SimulatorREBase<S, DEDS>::SimulatorREBase(DEDS& N,LHA_orig<typename DEDS::stateType>& A):SimulatorBase<S, EventsQueue, DEDS>(N,A){};
72 :
73 : template <class S>
74 1 : SPNBaseRE<S>::SPNBaseRE(bool doubleIS):SPNBase<S,EventsQueue>(0),doubleIS_mode(doubleIS){
75 1 : rareEventEnabled=false;
76 1 : Rate_Sum = 0;
77 1 : Origine_Rate_Sum = 0;
78 1 : Rate_Table = vector<double>(this->tr,0.0);
79 1 : Origine_Rate_Table = vector<double>(this->tr,0.0);
80 1 : }
81 :
82 : template <class S>
83 1 : void SPNBaseRE<S>::initialize( stateSpace *muprob){
84 1 : this->muprob=muprob;
85 1 : }
86 :
87 : template <class S, class DEDS>
88 1 : void SimulatorREBase<S, DEDS>::initVect(){
89 1 : muprob = new stateSpace();
90 1 : muprob->inputVect();
91 1 : this->N.initialize(muprob);
92 1 : }
93 :
94 : template<class S>
95 6000 : void SPNBaseRE<S>::InitialEventsQueue(EventsQueue &EQ,timeGen &TG) {
96 6000 : Rate_Sum = 0;
97 6000 : Origine_Rate_Sum = 0;
98 6000 : SPNBase<S,EventsQueue>::initialEventsQueue(EQ,TG);
99 6000 : }
100 :
101 : template <class S, class DEDS>
102 5347 : void SimulatorREBase<S, DEDS>::returnResultTrue(){
103 10694 : this->A.getFinalValues(
104 5347 : this->N.getState(),
105 : this->Result.quantR,
106 : this->Result.qualR);
107 5347 : this->Result.accept = true;
108 10694 : for(size_t i = 0; i< this->A.FormulaVal.size() ; i++)
109 5347 : this->Result.quantR[i] *= this->A.Likelihood;
110 5347 : if(P.verbose>3)cerr << "---------------\n TRUE: Likelyhood: "<< this->A.Likelihood <<" \n------\n";
111 5347 : }
112 :
113 : template <class S>
114 655847 : void SPNBaseRE<S>::update(double ctime,size_t t, const abstractBinding& b,EventsQueue &EQ,timeGen &TG){
115 : //If rareevent not require yet call the parent function
116 :
117 655847 : if(!rareEventEnabled){
118 0 : if(precondition(this->Marking)){
119 0 : rareEventEnabled = true;
120 : //A.Likelihood = 1.0;
121 : }else{
122 0 : SPNBase<S,EventsQueue>::update(ctime, t, b,EQ, TG);
123 0 : return;
124 : }
125 : }
126 :
127 1311694 : Event F;
128 :
129 655847 : Rate_Sum = 0;
130 655847 : Origine_Rate_Sum = 0;
131 :
132 : //Run over all transitions
133 3279235 : for (const auto &tr : this->Transition) {
134 2623388 : if(tr.Id != this->tr - 1){
135 3935082 : for(const auto &bindex : tr.bindingList ){
136 1967541 : if(this->IsEnabled(tr.Id, bindex)){
137 1783191 : if (EQ.isScheduled(tr.Id, bindex.id())) {
138 1735498 : generateEvent(ctime,F, tr.Id ,bindex, TG, *this );
139 1735498 : EQ.replace(F);
140 : } else {
141 47693 : generateEvent(ctime,F, tr.Id ,bindex, TG, *this );
142 47693 : EQ.insert(F);
143 : }
144 : }else{
145 184350 : if(EQ.isScheduled(tr.Id, bindex.id()))
146 41936 : EQ.remove(tr.Id,bindex.id());
147 : }
148 : }
149 : }
150 : }
151 :
152 1311694 : abstractBinding bpuit;
153 655847 : generateEvent(ctime,F, (this->tr-1), bpuit,TG, *this);
154 655847 : if(!doubleIS_mode){
155 655847 : EQ.replace(F);
156 : }
157 :
158 : /*
159 : //In Debug mode check that transition are scheduled iff they are enabled
160 : for (const auto &tr : N.Transition) {
161 : for(const auto &bindex : tr.bindingList){
162 : if (N.IsEnabled(tr.Id, bindex) !=
163 : EQ->isScheduled(tr.Id, bindex.id())){
164 : cerr << "N.IsEnabled(" << tr.label << ",";
165 : bindex.print();
166 : cerr <<")" << endl;
167 : if(EQ->isScheduled(tr.Id, bindex.id())){
168 : cerr << "Scheduled and not enabled!"<< endl;
169 : }else{
170 : cerr << "Enabled and not scheduled!" << endl;
171 : }
172 : assert(N.IsEnabled(tr.Id, bindex) ==
173 : EQ->isScheduled(tr.Id, bindex.id()));
174 : }
175 : }
176 : }
177 : */
178 :
179 : };
180 :
181 : template<class S, class DEDS>
182 661194 : void SimulatorREBase<S, DEDS>::updateLikelihood(size_t E1_transitionNum){
183 :
184 661194 : if(P.verbose>4){
185 0 : cerr << "initialised?:\t" << E1_transitionNum << "\t" << this->A.Likelihood << endl;
186 0 : cerr << this->N.Rate_Sum << "\t" << this->N.Origine_Rate_Sum << "\t[";
187 0 : for(let cv:this->N.Rate_Table)cerr << cv << ", ";
188 0 : cerr << "]\t[";
189 0 : for(let cv:this->N.Origine_Rate_Table)cerr << cv << ", ";
190 0 : cerr << "]" << endl;
191 : }
192 :
193 661194 : if(this->N.doubleIS_mode){
194 0 : this->A.Likelihood = this->A.Likelihood *
195 0 : (this->N.Origine_Rate_Table[E1_transitionNum] / this->N.Origine_Rate_Sum) *
196 0 : ((this->N.Rate_Sum-this->N.Rate_Table[this->N.tr-1]) / this->N.Rate_Table[E1_transitionNum]);
197 : }else{
198 1983582 : this->A.Likelihood = this->A.Likelihood *
199 : //
200 : //(N.Origine_Rate_Table[E1_transitionNum] / 1.0) *
201 1322388 : (this->N.Origine_Rate_Table[E1_transitionNum] / this->N.Origine_Rate_Sum) *
202 661194 : (this->N.Rate_Sum / this->N.Rate_Table[E1_transitionNum]);
203 : }
204 661194 : }
205 :
206 : template<class S, class DEDS>
207 661847 : bool SimulatorREBase<S, DEDS>::transitionSink(size_t i){
208 661847 : return (i==this->N.tr-1);
209 : }
210 :
211 : /*
212 : Useless for the moment
213 : void SimulatorRE::GenerateDummyEvent(Event& E, size_t Id) {
214 : E.transition = Id;
215 : E.time = 0.0;
216 : E.priority = N.GetPriority(Id);
217 : E.weight = 0.0;
218 : }*/
219 :
220 : template<class S, class DEDS>
221 6000 : void SimulatorREBase<S, DEDS>::reset(){
222 6000 : this->A.Likelihood=1.0;
223 6000 : this->N.reset();
224 6000 : this->A.reset(this->N.getState());
225 6000 : this->EQ->reset();
226 6000 : }
227 :
228 : /**
229 : * Simulate a whole trajectory in the system. Result is store in SimOutput
230 : */
231 : template<class S, class DEDS>
232 6000 : void SimulatorREBase<S, DEDS>::SimulateSinglePath() {
233 6000 : this->N.rareEventEnabled = this->N.precondition(this->N.getState());
234 6000 : static_cast<S*>(this)->reset();
235 6000 : this->N.InitialEventsQueue(*(this->EQ),*this);
236 :
237 6000 : if(this->logtrace.is_open())this->logtrace << "New Path"<< endl;
238 :
239 6000 : bool continueb = true;
240 1329694 : while ((!this->EQ->isEmpty()) && continueb ) {
241 : //cerr << "continue path"<< endl;
242 661847 : if(this->logtrace.is_open()){
243 0 : this->logtrace << this->A.CurrentTime << "\t";
244 0 : this->N.getState().print(this->logtrace,0.0);
245 0 : this->A.printState(this->logtrace);
246 0 : this->logtrace << endl;
247 : }
248 661847 : if(P.verbose>3){
249 : //Print marking and location of the automata
250 : //Usefull to track a simulation
251 0 : this->N.getState().printHeader(cerr);
252 0 : this->A.printHeader(cerr);
253 0 : cerr << endl;
254 0 : this->N.getState().print(cerr,0.0);
255 0 : this->A.printState(cerr);
256 0 : cerr << "\t" << this->A.Likelihood;
257 0 : cerr << endl;
258 0 : if(P.verbose>4)this->EQ->view(this->N.getTransitionLabels());
259 0 : if(P.verbose==6)this->interactiveSimulation();
260 : }
261 :
262 661847 : continueb = static_cast<S*>(this)->SimulateOneStep();
263 661847 : if(!this->N.rareEventEnabled)this->N.rareEventEnabled = this->N.precondition(this->N.getState());
264 : }
265 : //cerr << "finish path"<< endl;
266 6000 : }
267 :
268 : template <class S>
269 2457038 : void SPNBaseRE<S>::getParams(size_t Id,const abstractBinding& b){
270 : //If rareevent not require yet call the parent function
271 2457038 : if(!rareEventEnabled){
272 0 : this->GetDistParameters(Id,b);
273 0 : return;
274 : }
275 2457038 : this->GetDistParameters(Id,b);
276 2457038 : this->ParamDistr[1]= this->ParamDistr[0];
277 2457038 : this->ParamDistr[0]= ComputeDistr( Id, b, this->ParamDistr[0]);
278 : //N.ParamDistr[0]= N.ParamDistr[1]; /////////////////////////////////////////////////to remove
279 : }
280 :
281 : template <class S>
282 4252229 : double SPNBaseRE<S>::mu(){
283 :
284 8504458 : vector<int> vect (muprob->S.begin()->first->size(),0);
285 :
286 4252229 : lumpingFun(this->Marking,vect);
287 : //cerr << "test(";
288 4252229 : int i = muprob->findHash(&vect);
289 4252229 : if(i<0 || P.verbose>4){
290 0 : cerr << "state:(";
291 0 : for (size_t j =0; j < vect.size(); j++) {
292 0 : cerr << vect[j] << ",";
293 : }
294 0 : cerr << ") ->" << i << endl;
295 0 : this->Marking.printHeader(cerr);
296 0 : cerr << endl;
297 0 : this->Marking.print(cerr,0.0);
298 0 : cerr << endl;
299 0 : print_state(vect);
300 0 : if(i<0)exit(EXIT_FAILURE);
301 : }
302 4252229 : if(P.verbose>3)cerr << "muValue: " << muprob->getMu(i) << endl;
303 8504458 : return muprob->getMu(i);
304 : }
305 :
306 : template <class S>
307 2457038 : double SPNBaseRE<S>::ComputeDistr(size_t t , const abstractBinding& b, double origin_rate){
308 2457038 : if(P.verbose>4)cerr << "trans: " << this->Transition[t].label << " mu origine:";
309 2457038 : double mux = mu();
310 2457038 : if( mux==0.0 || mux==1.0) return(origin_rate);
311 :
312 2457038 : if(t== this->tr-1){
313 661847 : if(Origine_Rate_Sum >= Rate_Sum){
314 358569 : return( Origine_Rate_Sum - Rate_Sum );
315 : }else{
316 303278 : if(P.verbose>3 && (Origine_Rate_Sum < 0.99*Rate_Sum)){
317 0 : cerr << "Reduce model does not guarantee variance" << endl;
318 0 : cerr << "Initial sum of rate: " << Origine_Rate_Sum << " Reduce one: " << Rate_Sum << " difference: " << Origine_Rate_Sum - Rate_Sum << endl ;
319 : //exit(EXIT_FAILURE);
320 : }
321 :
322 303278 : return 0.0 ;};
323 : };
324 1795191 : if(P.verbose>4)cerr << "mu target : ";
325 : double distr;
326 1795191 : this->fire(t,b,0.0);
327 1795191 : distr = origin_rate *( mu() / mux);
328 1795191 : this->unfire(t,b);
329 1795191 : if(P.verbose>4)cerr <<endl;
330 1795191 : return(distr);
331 : }
|
