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 numericalSolverSH.cpp created by Benoit Barbot on 06/01/12. *
24 : *******************************************************************************
25 : */
26 :
27 : #include <time.h>
28 : #include "numSolverSH.hpp"
29 :
30 0 : bool numSolverSH::readbit(int a,int b){
31 0 : int c = a & (1 << b);
32 0 : return (0 != c);
33 : }
34 :
35 :
36 0 : void numSolverSH::initVect(int nT){
37 : time_t start, endt;
38 0 : time(&start);
39 :
40 0 : T=nT;
41 0 : u=T;
42 0 : l = (int)log2(T);
43 :
44 : //cerr << "T: " <<T << " l: " << l << endl;
45 :
46 :
47 : //powTVect = new vector< boostmat::vector<double> > (l+1, *finalVector);
48 0 : lastOne = new vector< boostmat::vector<double> > (l+2, *finalVector);
49 0 : ktable = new vector< double > (l+2,0);
50 0 : lastPowT= 1<<l;
51 :
52 : //cerr << "lastPow" << lastPowT << " log(T) " << l << endl;
53 :
54 0 : boostmat::vector<double> itervect = *finalVector;
55 : boostmat::vector<double> itervect2=
56 0 : boostmat::zero_vector<double> (finalVector->size());
57 :
58 0 : time(&endt);
59 0 : cerr << "time for allocation:" << difftime(endt, start) << endl;
60 :
61 : //We suppose here that the initial state is the first of the vector
62 0 : if(itervect (0) != 0.0)minT=0;
63 :
64 0 : for(int i=1; i<=lastPowT ; i++){
65 : //cerr << "currPow " << currPow << " newPow " << nextPow << endl;
66 0 : itervect2.clear();
67 0 : sparseProd(&itervect2,&itervect, transitionsMatrix);
68 0 : itervect=itervect2;
69 :
70 0 : if((itervect (0) != 0.0) && (minT== -1))minT=i;
71 : //cerr << "i: " <<i << " v: " <<itervect << endl;
72 : //itervect = boostmat::prod ((*transitionsMatrix), itervect);
73 : }
74 0 : (*lastOne)[l]= itervect;
75 0 : (*ktable)[l] = lastPowT;
76 0 : (*ktable)[l+1] = 0;
77 :
78 0 : time(&endt);
79 0 : cerr << "time for precalculation:" << difftime(endt, start) << endl;
80 :
81 0 : cerr << "Starting the simulation" << endl;
82 : //cerr << "finish init" << endl;
83 0 : }
84 :
85 0 : void numSolverSH::compPow(int kp,int up){
86 : //cerr << "compPow("<< kp << "," << up << ")" << endl;
87 :
88 0 : int k = kp-1;
89 0 : boostmat::vector<double> itervect = (*lastOne)[kp];
90 : boostmat::vector<double> itervect2 =
91 0 : boostmat::zero_vector<double> (finalVector->size());
92 :
93 0 : while (k>0 && readbit(up, k)==0){
94 : //cerr << itervect<< endl;
95 0 : (*lastOne)[k]=itervect;
96 0 : (*ktable)[k] =(*ktable)[kp];
97 : //cerr << "k--" << endl;
98 0 : k--;
99 : }
100 :
101 0 : int m = 1<<k; // (m = 2^k
102 0 : for (int i = (int)(*ktable)[kp]+1; i<=up; i++) {
103 : //cerr << "i: " << i << " k: " << k << " m: " << m << endl;
104 0 : itervect2.clear();
105 0 : sparseProd(&itervect2,&itervect, transitionsMatrix);
106 0 : itervect=itervect2;
107 :
108 0 : if((itervect (0) != 0.0) && (minT== -1))minT=i;
109 :
110 : //itervect = boostmat::prod ((*transitionsMatrix), itervect);
111 0 : m--;
112 0 : if(m==0){
113 0 : (*lastOne)[k]=itervect;
114 0 : (*ktable)[k] = i;
115 0 : k--;
116 : //cerr << "k--1" << endl;
117 0 : while (k>0 && readbit(up, k)==0){
118 0 : (*lastOne)[k]=itervect;
119 0 : (*ktable)[k] =i;
120 : //cerr << "k--" << endl;
121 0 : k--;
122 : }
123 :
124 0 : m= 1<<k;
125 : }
126 :
127 : }
128 : //cerr << "finish";
129 : //cerr << itervect << endl;
130 0 : previous_vect = current_vect;
131 0 : current_vect=itervect;
132 0 : is_previous=false;
133 : /*current_vect= boostmat::zero_vector<double> (itervect.size());
134 : sparseProd(¤t_vect, &itervect, transitionsMatrix);*/
135 : //current_vect = boostmat::prod ((*transitionsMatrix), itervect);
136 : //cerr << "finish reset" << endl;
137 :
138 :
139 0 : }
140 :
141 0 : void numSolverSH::reset(){
142 : time_t start, endt;
143 0 : time(&start);
144 :
145 0 : u=T;
146 0 : compPow(l, T);
147 :
148 0 : time(&endt);
149 0 : cerr << "time for reset:" << difftime(endt, start) << endl << endl;
150 :
151 : //cerr << "finish reset" << endl;
152 :
153 :
154 0 : }
155 :
156 0 : void numSolverSH::switchOff(){
157 0 : u= -1;
158 0 : }
159 :
160 0 : void numSolverSH::stepVect(){
161 : //cerr << "step vect u:" << u << "->" << u-1 << "(";
162 0 : u--;
163 0 : if(is_previous){ is_previous=false;}
164 : else {
165 0 : if(u==0){
166 0 : previous_vect = current_vect;
167 0 : current_vect = *finalVector;
168 0 : }else if(u>0){
169 0 : int kp = (int)log2((u^(u+1))+1);
170 : //cerr << "u: " << u<< ": kp: " << kp << endl;
171 0 : compPow(kp, u);
172 : }
173 : }
174 :
175 : //cerr << ")" << endl;
176 0 : }
177 :
178 0 : void numSolverSH::previousVect(){
179 0 : u++;
180 0 : if(is_previous){
181 0 : cerr << "fail to compute previous vect" << endl;
182 : }else{
183 0 : is_previous=true;
184 : }
185 0 : }
186 :
187 0 : boostmat::vector<double> numSolverSH::getVect(){
188 0 : if(is_previous){
189 0 : return previous_vect;
190 0 : }else return current_vect;
191 : }
192 :
193 0 : double numSolverSH::getMu(int i)const{
194 0 : if(u<0)return 1.0;
195 :
196 0 : if(is_previous){
197 0 : return previous_vect[i];
198 0 : }else return current_vect[i];
199 : }
200 :
201 0 : int numSolverSH::currentRound(){
202 0 : return T-u;
203 : }
204 :
205 :
206 0 : void numSolverSH::printState(){
207 0 : cerr << u;
208 3 : }
|
