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 numericalSolverBB.cpp created by Benoit Barbot on 06/01/12. *
24 : *******************************************************************************
25 : */
26 :
27 : #include "numSolverBB.hpp"
28 :
29 0 : void numSolverBB::initVect(int nT){
30 0 : T=nT;
31 0 : u=T;
32 0 : l = (int)sqrt(T);
33 : time_t start, endt;
34 0 : time(&start);
35 :
36 0 : minT = -1;
37 :
38 0 : circularvect = new vector< boostmat::vector<double> > (l+1, boostmat::zero_vector<double> (finalVector->size()));
39 0 : checkPoint = new vector< boostmat::vector<double> > (T/l+1, *finalVector);
40 :
41 0 : time(&endt);
42 0 : cerr << "time for allocation:" << difftime(endt, start) << endl;
43 :
44 0 : lastCP = (T/l)*l;
45 0 : boostmat::vector<double> itervect = *finalVector;
46 0 : boostmat::vector<double> itervect2 = boostmat::zero_vector<double> (finalVector->size());
47 : //We suppose here that the initial state is the first of the vector
48 0 : if(itervect(0) != 0)minT=0;
49 :
50 :
51 0 : for(int i=1; i<=lastCP ; i++){
52 0 : itervect2.clear();
53 0 : sparseProd(&itervect2, &itervect, transitionsMatrix);
54 0 : itervect = itervect2;
55 :
56 0 : if((itervect (0) != 0) && (minT== -1))minT=i;
57 : //cerr << "i:" << i << "v0:\t" << itervect (0) << "mint:\t" << minT<< endl;
58 :
59 :
60 : //itervect = boostmat::prod ((*transitionsMatrix), itervect);
61 0 : if( i % l ==0) (*checkPoint)[i/l]= itervect;
62 : }
63 :
64 0 : time(&endt);
65 0 : cerr << "time for precalculation:" << difftime(endt, start) << endl;
66 :
67 0 : cerr << "Starting the simulation" << endl;
68 :
69 0 : }
70 :
71 0 : void numSolverBB::reset(){
72 : time_t start, endt;
73 0 : time(&start);
74 :
75 0 : (*circularvect)[0]=(*checkPoint)[lastCP/l];
76 :
77 : /*for(int i=1; i<=T-lastCP ; i++){
78 : (*circularvect)[i] = boostmat::prod ((*transitionsMatrix), (*circularvect)[i-1]);
79 : }*/
80 0 : for(int i=1; i<=T-lastCP ; i++){
81 0 : ((*circularvect)[i]).clear();
82 0 : sparseProd(&((*circularvect)[i]), &((*circularvect)[i-1]), transitionsMatrix);
83 : //cerr << "i:" << i << "v0:\t" << (*circularvect)[i] (0) << "mint:\t" << minT<< endl;
84 0 : if(((*circularvect)[i] (0) != 0.0) && (minT== -1))minT=i+lastCP;
85 :
86 : //cerr << "itervect:" << i << ":"<< (*circularvect)[i] << endl;
87 : }
88 0 : time(&endt);
89 0 : cerr << "time for reset:" << difftime(endt, start) << endl << endl;
90 :
91 :
92 0 : u=T;
93 0 : matOffset = T-lastCP;
94 0 : }
95 :
96 :
97 0 : void numSolverBB::stepVect(){
98 : //cerr << "step:" << u << endl;
99 0 : if(matOffset > 0 || u==0){
100 0 : matOffset--;
101 0 : u--;
102 : }else{
103 0 : u--;
104 0 : (*circularvect)[l]=(*circularvect)[0];
105 0 : (*circularvect)[0]=(*checkPoint)[u/l];
106 0 : for(int i=1; i<l ; i++){
107 0 : (*circularvect)[i].clear();
108 0 : sparseProd(&((*circularvect)[i]), &((*circularvect)[i-1]), transitionsMatrix);
109 : //(*circularvect)[i] = boostmat::prod ((*transitionsMatrix), (*circularvect)[i-1]);
110 : }
111 0 : matOffset=l-1;
112 : }
113 0 : }
114 :
115 0 : void numSolverBB::previousVect(){
116 0 : u++;
117 0 : matOffset++;
118 0 : }
119 :
120 :
121 0 : void numSolverBB::printState(){
122 0 : cerr << u;
123 0 : }
124 :
125 0 : int numSolverBB::currentRound(){
126 0 : return T-u;
127 3 : }
128 :
129 :
|
