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 numericalSolver.cpp created by Benoit Barbot on 05/12/11. *
24 : *******************************************************************************
25 : */
26 :
27 : #include "numericalSolver.hpp"
28 :
29 : #include <iostream>
30 : #include <fstream>
31 : #include <time.h>
32 : #include <boost/numeric/ublas/matrix_sparse.hpp>
33 : #include <boost/numeric/ublas/matrix.hpp>
34 : #include <boost/numeric/ublas/matrix_expression.hpp>
35 : #include <boost/numeric/ublas/io.hpp>
36 : #include <boost/numeric/ublas/operation_sparse.hpp>
37 : #include <boost/numeric/ublas/operation.hpp>
38 : #include <boost/numeric/ublas/vector.hpp>
39 :
40 : #define BOOST_UBLAS_NDEBUG
41 :
42 : using namespace std;
43 :
44 :
45 0 : numericalSolver::numericalSolver(){
46 0 : cerr << "Initialising the simulator" << endl;
47 0 : cerr << "Reading Matrix" << endl;
48 : time_t start, endt;
49 0 : time(&start);
50 0 : inputMat();
51 0 : time(&endt);
52 0 : cerr << "time for input:" << endt-start << endl;
53 0 : }
54 :
55 :
56 0 : void numericalSolver::sparseProd(boostmat::vector<double> *result,boostmat::vector<double> *vect, boostmat::compressed_matrix<double> *mat){
57 : //*result = boostmat::zero_vector<double> (vect->size());
58 :
59 : //boostmat::axpy_prod((const boostmat::compressed_matrix<double,boostmat::row_major>) *mat,(const boostmat::vector<double>) *vect,*result,true);
60 :
61 :
62 : //boostmat::axpy_prod(*mat, *vect, *result, boostmat::row_major );
63 :
64 : //boostmat::sparse_prod<double>(mat, vect);
65 :
66 : //boostmat::sparse_prod(mat, vect, result);
67 :
68 :
69 :
70 0 : for(boostmat::compressed_matrix<double>::iterator1 it = mat->begin1(); it!= mat->end1(); it++){
71 0 : for(boostmat::compressed_matrix<double>::iterator2 it2 = it.begin(); it2!= it.end(); it2++){
72 : //cerr << "iteration: " << it.index1() << ":" << it2.index2() << endl;
73 0 : (*result) (it.index1()) += ((*vect) (it2.index2()) * (*it2));
74 : };
75 : };
76 0 : }
77 :
78 0 : void numericalSolver::initVect(int nT){
79 : time_t start, endt;
80 0 : time(&start);
81 0 : minT=-1;
82 :
83 0 : T=nT;
84 0 : circularvect = new vector< boostmat::vector<double> > (nT+1, boostmat::zero_vector<double> (finalVector->size()));
85 0 : (*circularvect)[0] = *finalVector;
86 0 : time(&endt);
87 0 : cerr << "time for allocation:" << difftime(endt, start) << endl;
88 :
89 : //cerr << "itervect:" << 0 << ":"<< (*circularvect)[0] << endl;
90 :
91 : //We suppose here that the initial state is the first of the vector
92 0 : if( !((*circularvect)[0] (0) == 0))minT=0;
93 :
94 :
95 : //int test=0;
96 0 : boostmat::vector<double> itervect = boostmat::vector<double>(*finalVector);
97 : //boostmat::vector<double> itervect2(finalVector->size());
98 0 : for(int i=1; i<=nT ; i++){
99 :
100 0 : sparseProd(&((*circularvect)[i]), &((*circularvect)[i-1]), transitionsMatrix);
101 :
102 0 : if(((*circularvect)[i] (0) != 0) && (minT== -1))minT=i;
103 : //cerr << "i:" << i<< "\t"<<(*circularvect)[i] (0) << endl;
104 :
105 :
106 : //boostmat::sparse_prod(transitionsMatrix, (*circularvect)[i] , (*circularvect)[i-1] );
107 :
108 :
109 : //cerr << "itervect:" << i << ":"<< (*circularvect)[i] << endl;
110 : }
111 0 : nbVect = nT;
112 0 : matOffset = nT;
113 0 : time(&endt);
114 0 : cerr << "time for precalculation:" << difftime(endt, start) << endl;
115 :
116 0 : cerr << "Starting the simulation" << endl;
117 0 : }
118 :
119 0 : int numericalSolver::currentRound(){
120 0 : return T-matOffset;
121 : }
122 :
123 0 : void numericalSolver::switchOff(){
124 0 : matOffset=-1;
125 0 : }
126 :
127 0 : void numericalSolver::reset(){
128 0 : matOffset=T;
129 : //cerr << (*circularvect)[matOffset] << endl;
130 :
131 0 : }
132 :
133 0 : boostmat::vector<double> numericalSolver::getVect(){
134 0 : if (matOffset<0) {
135 0 : return boostmat::vector<double>(1,0.0);
136 : }
137 0 : return (*circularvect)[matOffset];
138 : }
139 :
140 0 : double numericalSolver::getMu(int i)const{
141 : //cerr << "indice" << matOffset << endl;
142 0 : if (matOffset>=0) return (*circularvect)[matOffset][i];
143 0 : else return 1.0;
144 : }
145 :
146 0 : void numericalSolver::stepVect(){
147 0 : matOffset--;
148 : //cerr << " # "<< matOffset <<" # ";
149 0 : }
150 :
151 0 : void numericalSolver::previousVect(){
152 0 : matOffset++;
153 0 : }
154 :
155 0 : void numericalSolver::printState(){
156 0 : cerr << matOffset;
157 3 : }
158 :
159 :
160 :
161 :
162 :
163 :
164 :
165 :
166 :
|
