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 stateSpace.cpp created by Benoit Barbot on 01/12/11. *
24 : *******************************************************************************
25 : */
26 :
27 : #include "stateSpace.hpp"
28 : #include <deque>
29 : #include <set>
30 : #include <iostream>
31 : #include <fstream>
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 :
37 : #include "sparse_io.hpp"
38 :
39 : namespace boostmat = boost::numeric::ublas;
40 : //using namespace boost::numeric::ublas;
41 :
42 4 : stateSpace::stateSpace(){
43 4 : nbState=0;
44 4 : nbTrans=0;
45 4 : maxRate=0;
46 4 : }
47 :
48 4253771 : int stateSpace::findHash(const vector<int>* vect)const{
49 4253771 : const auto it = S.find (vect);
50 4253771 : if (it != S.end ())
51 4253771 : return(it->second); // found
52 : else
53 0 : return(-1);
54 : }
55 :
56 4252229 : double stateSpace::getMu(int state)const{
57 4252229 : return (*muvect)[state];
58 : }
59 :
60 602 : void stateSpace::add_state(vector<int> v){
61 602 : vector<int>* v2 = new vector<int>(v);
62 602 : S[v2] = (int)nbState;
63 : //for(int i = 0 ; i< v2->size(); i++)cerr << (*v2)[i]<< ",";
64 : //cerr << endl;
65 602 : findstate->push_back(*v2);
66 602 : nbState++;
67 602 : if((nbState % 100000) ==0)cerr << "Number of states :" <<nbState<< endl;
68 :
69 602 : }
70 :
71 3 : void stateSpace::exploreStateSpace(){
72 : // apply a Dijkstra algorithm on the product of the SPN an the LHA to produce
73 : // the state space. The list of state is store in the hash table S and
74 : // the transition list is stored in transitionList.
75 :
76 3 : N.reset();
77 3 : cerr << "Exploring state space" << endl;
78 :
79 6 : stack<vector<int> ,vector<vector<int> > > toBeExplore;
80 6 : auto init = N.getState().getVector();
81 :
82 3 : A.reset(N.getState());
83 3 : init.push_back(A.CurrentLocation);
84 3 : toBeExplore.push(init);
85 3 : findstate = new vector<vector<int> >(0);
86 3 : add_state(init);
87 :
88 6 : vector<size_t> immTrans;
89 12 : for (size_t t = 0; t < N.tr; t++)
90 9 : if( N.Transition[t].DistTypeIndex == IMMEDIATE || N.Transition[t].DistTypeIndex == DETERMINISTIC )
91 0 : immTrans.push_back(t);
92 :
93 1207 : while (!toBeExplore.empty()) {
94 1204 : vector<int> place = toBeExplore.top();
95 602 : toBeExplore.pop();
96 :
97 1204 : vector<int> currentstate = place;
98 :
99 602 : int lhaloc = place.back();
100 602 : A.CurrentLocation = lhaloc;
101 602 : place.pop_back();
102 602 : N.setState(place);
103 :
104 602 : bool existImm = false;
105 :
106 602 : for (size_t t = 0; t < immTrans.size()&& !existImm ; t++)
107 0 : for(let b : N.Transition[immTrans[t]].bindingList)
108 0 : if (N.IsEnabled(immTrans[t],b)){
109 0 : existImm = true;
110 0 : break;
111 : }
112 1204 : const auto savMark = N.getState();
113 :
114 2408 : for (size_t t = 0; t < N.tr; t++)
115 3612 : for(let b : N.Transition[t].bindingList){
116 :
117 1806 : N.setState(savMark);
118 1806 : A.CurrentLocation = lhaloc;
119 :
120 4988 : if (N.IsEnabled(t,b) &&
121 1591 : (!existImm || N.Transition[t].DistTypeIndex== IMMEDIATE
122 0 : || N.Transition[t].DistTypeIndex== DETERMINISTIC )) {
123 :
124 1591 : N.fire(t,b,0.0);
125 1591 : N.symmetrize();
126 :
127 1591 : int SE = A.synchroniseWith(t,N.getState(),b);
128 :
129 1591 : if (SE > -1) {
130 : //cerr << "transition:" << *it << endl;
131 3084 : vector<int> marking = N.getState().getVector();
132 :
133 1542 : nbTrans++;
134 1542 : marking.push_back( A.CurrentLocation );
135 : //vector<double> Param = N.GetDistParameters(*it);
136 : //transitionsList.push( make_pair(make_pair(currentstate, marking),Param[0] ));
137 :
138 1542 : auto its = S.find (&marking);
139 1542 : if (its == S.end ()){
140 :
141 : /*
142 : //N.Marking.printHeader(cerr);
143 : cerr << "state:"<< nbState << " -> ";
144 : for( let i : marking)cerr << i << " ";
145 : //N.Marking.print(cerr, 0.0);
146 : cerr << endl;
147 : */
148 :
149 :
150 599 : toBeExplore.push(marking);
151 599 : add_state(marking);
152 : } //else { cerr << " -> " << its->second << endl;}
153 : }
154 : //N.unfire(t,b);
155 : }
156 :
157 : }
158 :
159 : }
160 3 : cerr << nbState << " states found" << endl
161 6 : << nbTrans << " transitions found" << endl;
162 3 : }
163 :
164 3 : void stateSpace::buildTransitionMatrix()
165 : {
166 3 : cerr << "Building transition matrix" << endl;
167 :
168 : // transform the transition list into a sparse transition probability matrix
169 3 : transitionsMatrix = new boost::numeric::ublas::compressed_matrix<double>(nbState, nbState, nbTrans);
170 3 : auto &mat = *transitionsMatrix;
171 :
172 3 : cerr << "Exploring graph" << endl;
173 :
174 605 : for (size_t i=0; i<nbState; i++) {
175 1204 : vector<int> place = (*findstate)[i];
176 :
177 602 : int lhaloc = place.back();
178 602 : place.pop_back();
179 602 : N.setState(place);
180 :
181 : /*cerr << "state:";
182 : for (auto it2=currentstate.begin(); it2!= currentstate.end() ; it2++) {
183 : cerr << *it2 << ":";
184 : }
185 : cerr << endl;*/
186 :
187 602 : mat(i,i) = 1.0;
188 2408 : for (size_t t = 0; t < N.tr; t++){
189 : //Loop over binding here
190 3612 : abstractBinding b;
191 1806 : if (N.IsEnabled(t,b)) {
192 1591 : A.CurrentLocation = lhaloc;
193 1591 : N.fire(t,b,0.0);
194 3182 : vector<int> marking = N.getState().getVector();
195 1591 : int SE = A.synchroniseWith( t , N.getState(),b);
196 1591 : if (SE > -1) {
197 1542 : N.unfire(t,b);
198 1542 : marking.push_back( A.CurrentLocation );
199 1542 : N.GetDistParameters(t,b);
200 :
201 1542 : int j = findHash(&marking);
202 1542 : mat (i,j) = N.ParamDistr[0];
203 :
204 : /*cerr << "->state:";
205 : for (auto it2=marking.begin(); it2!= marking.end() ; it2++) {
206 : cerr << *it2 << ":";
207 : }
208 : cerr << N.ParamDistr[0] << endl;*/
209 49 : }else N.unfire(t,b);
210 : }
211 : }
212 :
213 :
214 : }
215 :
216 : /* Quick fix to redo */
217 :
218 : /*
219 : cerr << "uniformize to 1" << endl;
220 : for (it1_t it1 = mat.begin1(); it1 != mat.end1(); it1++)
221 : {
222 : double sum = 0.0;
223 : for (it2_t it2 = it1.begin(); it2 != it1.end(); it2++){
224 : //cerr << "non null:" << it2.index1() << ":" << it2.index2() << endl;
225 : if(it2.index1()!= it2.index2())sum += *it2;
226 : }
227 : for (it2_t it2 = it1.begin(); it2 != it1.end(); it2++){
228 : //cerr << "non null:" << it2.index1() << ":" << it2.index2() << endl;
229 : *it2 /= sum;
230 : }
231 : mat(it1.index1(),it1.index1())= 1.0;
232 : }
233 : */
234 :
235 :
236 3 : cerr << "Adding self loop" << endl;
237 : // Add self loop to ensure that mat is a probability matrix.
238 : // If the model is a CTMC the value on diagonal are wrong.
239 :
240 605 : for (auto it1 = mat.begin1(); it1 != mat.end1(); it1++)
241 : {
242 602 : double sum = 1.0;
243 2746 : for (auto it2 = it1.begin(); it2 != it1.end(); it2++){
244 : //cerr << "non null:" << it2.index1() << ":" << it2.index2() << endl;
245 2144 : if(it2.index1()!= it2.index2())sum -= *it2;
246 : }
247 602 : mat(it1.index1(),it1.index1())= sum;
248 : }
249 :
250 3 : cerr << " copying" << endl;
251 :
252 :
253 3 : finalVector = new boost::numeric::ublas::vector<double> (nbState);
254 605 : for(hash_state::iterator it=S.begin(); it!=S.end() ; it++){
255 602 : A.CurrentLocation = it->first->back();
256 602 : if(A.isFinal()){
257 17 : (*finalVector)(it->second)=1.0;
258 : //cerr << "final:" << it->second << endl;
259 : }else {
260 585 : (*finalVector)(it->second)=0.0;
261 : }
262 :
263 : }
264 :
265 3 : }
266 :
267 : /*double stateSpace::maxRate(){
268 : double t = 0.0;
269 : for(boost::numeric::ublas::compressed_matrix<double>::iterator1 it
270 : = transitionsMatrix->begin1(); it!= transitionsMatrix->end1(); it++){
271 : for(boost::numeric::ublas::compressed_matrix<double>::iterator2 it2
272 : = it.begin(); it2!= it.end(); it2++){
273 : if(it.index1() != it2.index2())t = max(t,*it2);
274 : };
275 : };
276 : return(t);
277 : }*/
278 :
279 0 : double stateSpace::uniformizeMatrix(){
280 : //First Compute infinitesimal generator
281 : //replace all value on the diagonal by opposite of the sum
282 :
283 0 : double lambda = 0.0;
284 0 : for (auto it1 = transitionsMatrix->begin1(); it1 != transitionsMatrix->end1(); it1++)
285 : {
286 0 : double sum = 0.0;
287 0 : for (auto it2 = it1.begin(); it2 != it1.end(); it2++){
288 : //cerr << "non null:" << it2.index1() << ":" << it2.index2() << endl;
289 0 : if(it2.index1()!= it2.index2())sum += *it2;
290 : }
291 0 : lambda = max(lambda ,sum);
292 0 : (*transitionsMatrix)(it1.index1(),it1.index1())= -sum;
293 : }
294 : // Divide each coefficient of the matrix by lambda
295 : // and add 1 on the diagonal
296 :
297 0 : for (auto it1 = transitionsMatrix->begin1(); it1 != transitionsMatrix->end1(); it1++)
298 : {
299 0 : for (auto it2 = it1.begin(); it2 != it1.end(); it2++){
300 : //cerr << "non null:" << it2.index1() << ":" << it2.index2() << endl;
301 0 : *it2 /= lambda;
302 0 : if(it2.index1()== it2.index2())*it2 +=1.0;
303 : }
304 : }
305 0 : maxRate = lambda;
306 0 : return lambda;
307 : }
308 :
309 0 : void stateSpace::printP(){
310 0 : cerr << "Probability transition matrix:" << endl;
311 0 : for(size_t i=0; i< transitionsMatrix->size1() ; i++){
312 0 : for(size_t j = 0;j< transitionsMatrix->size2() ; j++)
313 0 : cerr << (*transitionsMatrix)(i,j) << "\t";
314 0 : cerr << endl;
315 : }
316 0 : cerr << endl << "Final Vector"<< endl;
317 0 : for(size_t i=0; i< finalVector->size() ; i++){
318 0 : cerr << (*finalVector)(i);
319 0 : cerr << endl;
320 : }
321 0 : }
322 :
323 1 : void stateSpace::outputMat(){
324 1 : cerr << "Exporting the transition matrix" << endl;
325 :
326 2 : fstream outputFile;
327 1 : outputFile.open("matrixFile",fstream::out);
328 :
329 1 : outputFile << boostmat::io::sparse(*transitionsMatrix);
330 1 : outputFile << *finalVector << endl;
331 :
332 292 : for(hash_state::iterator it= S.begin() ; it != S.end(); it++){
333 291 : outputFile << "(";
334 582 : vector<int> vect = *(*it).first;
335 873 : for(size_t i=0; i< N.Msimpletab.size();i++){
336 582 : if(i>0)outputFile << ",";
337 582 : outputFile << vect[N.Msimpletab[i]];
338 : };
339 :
340 : /*for(int i =0; i< vect.size()-1; i++){
341 : if(i>0)outputFile << ",";
342 : outputFile << vect[i];
343 : }*/
344 291 : outputFile << ")=";
345 291 : outputFile << (*it).second << endl;
346 : }
347 :
348 1 : outputFile.close();
349 1 : }
350 :
351 2 : void stateSpace::outputPrism(){
352 2 : cerr << "Exporting the model for Prism" << endl;
353 :
354 4 : fstream outputFile;
355 2 : outputFile.open("prismStates.sta",fstream::out);
356 :
357 2 : outputFile << "(" ;
358 7 : for(size_t i=0; i< N.Place.size();i++){
359 5 : outputFile << N.Place[i].label ;
360 5 : outputFile << ",";
361 : };
362 2 : outputFile << "automata)" << endl;
363 :
364 313 : for(size_t it=0 ; it < findstate->size(); it++){
365 311 : outputFile << it << ":(";
366 622 : vector<int> vect = (*findstate)[it];
367 1224 : for(size_t i=0; i< N.Place.size();i++){
368 913 : outputFile << vect[i];
369 913 : outputFile << ",";
370 : };
371 311 : outputFile << vect[vect.size()-1];
372 :
373 311 : outputFile << ")" << endl;
374 : }
375 :
376 2 : outputFile.close();
377 :
378 4 : fstream outputMatrixFile;
379 2 : outputMatrixFile.open("prismMatrix.tra",fstream::out);
380 2 : outputMatrixFile << nbState << " " << nbTrans << endl;
381 :
382 313 : for (auto it1 = transitionsMatrix->begin1(); it1 != transitionsMatrix->end1(); it1++)
383 : {
384 1410 : for (auto it2 = it1.begin(); it2 != it1.end(); it2++){
385 1099 : if( *it2 >= (1e-16)){
386 907 : outputMatrixFile << it2.index1() << " " << it2.index2() << " " << *it2 << endl;
387 : }
388 : }
389 : }
390 2 : outputMatrixFile.close();
391 :
392 4 : fstream outputProperty;
393 2 : outputProperty.open("prismProperty.ctl",fstream::out);
394 2 : outputProperty << "P=? [ (true) U (";
395 2 : bool first =true;
396 2 : int lhaloc = A.CurrentLocation;
397 8 : for(int it = 0 ;it < (int)A.NbLoc; ++it){
398 6 : A.CurrentLocation = it;
399 6 : if(A.isFinal()){
400 2 : if (first){first=false;}else{outputProperty << "|";}
401 2 : outputProperty << "(automata = " << it << ")";
402 : }
403 : }
404 2 : A.CurrentLocation = lhaloc;
405 2 : outputProperty << ")]";
406 2 : outputProperty.close();
407 :
408 4 : fstream outputLabel;
409 2 : outputLabel.open("prismLabel.lbl",fstream::out);
410 2 : outputLabel << "0='init' 1='deadlock'\n0: 0";
411 2 : outputLabel.close();
412 :
413 2 : }
414 :
415 2 : void stateSpace::launchPrism(string prismPath){
416 2 : cerr << "Starting Prism"<< endl;
417 2 : string cmd = prismPath + " -gs -maxiters 1000000000 -ctmc -importtrans prismMatrix.tra -importstates prismStates.sta -importlabels prismLabel.lbl -v -cuddmaxmem 10000000 prismProperty.ctl > prismOutput";
418 2 : if(0 != system(cmd.c_str())){
419 2 : cerr << "Fail to launch prism" << endl;
420 2 : exit(EXIT_FAILURE);
421 : }
422 0 : cerr << "Prism finish" << endl;
423 0 : }
424 :
425 0 : void stateSpace::importPrism(){
426 0 : cerr << "Importing Prism result" << endl;
427 0 : string line;
428 0 : size_t poseq =1;
429 0 : string pos;
430 0 : string prob;
431 0 : ifstream myfile ("prismOutput");
432 0 : if (myfile.is_open())
433 : {
434 0 : do{
435 0 : getline (myfile,line);
436 0 : }while(myfile.good() &&
437 0 : line != "Probabilities (non-zero only) for all states:");
438 :
439 0 : muvect = new vector<double> (nbState,0.0);
440 : //int n=0;
441 0 : while ( myfile.good() && poseq>0)
442 : {
443 0 : getline (myfile,line);
444 : //cerr << line << endl;
445 0 : poseq = line.find("=");
446 :
447 0 : if(poseq != string::npos){
448 : //cerr << line << endl;
449 0 : size_t si = 1+line.find("(",0);
450 0 : pos = line.substr(si,poseq-1-si);
451 : //cerr << "pos:" << pos << endl;
452 0 : prob = line.substr(poseq+1,line.size());
453 :
454 0 : vector<int> vect;
455 0 : size_t it = 0;
456 : //cerr << "v:";
457 0 : while(it < pos.length()){
458 0 : size_t it2 = pos.find(",",it);
459 0 : if(it2 == string::npos ) it2 = pos.length();
460 : //cerr << "test:" << it<< ":" << it2 << endl;
461 0 : vect.push_back(atoi((pos.substr(it,it2-it)).c_str() ));
462 : //cerr << atoi((pos.substr(it,it2-it)).c_str() ) << ",";
463 0 : it = it2+1;
464 : }
465 : //cerr << endl;
466 :
467 0 : int state = findHash(&vect);
468 : //cerr << "state" << state << ":";
469 0 : (*muvect)[state] = atof(prob.c_str());
470 : //cerr << atof(prob.c_str());
471 : //muvect->push_back(atof(prob.c_str()));
472 : //S[new vector<int>(vect)] = n;
473 : //n++;
474 :
475 :
476 : }
477 : }
478 0 : myfile.close();
479 : //nbState = n;
480 : }
481 0 : }
482 :
483 1 : void stateSpace::outputTmpLumpingFun(){
484 1 : cerr << "Exporting the temporary lumping function" << endl;
485 2 : fstream outputlumptmp;
486 :
487 1 : outputlumptmp.open("lumpingfunTmp.cpp",fstream::out);
488 1 : outputlumptmp << "#include \"markingImpl.hpp\"" << endl << endl;
489 :
490 3 : for(size_t i=0; i< N.Msimpletab.size();i++){
491 2 : int j = N.Msimpletab[i];
492 2 : outputlumptmp << "const int reducePL_" << N.Place[j].label.substr(1,N.Place[j].label.length()-1 ) << " = " << i << ";" << endl;
493 : };
494 :
495 1 : outputlumptmp << "void SPN::print_state(const vector<int> &vect){" << endl;
496 3 : for(size_t i=0; i< N.Msimpletab.size();i++){
497 2 : int j = N.Msimpletab[i];
498 2 : outputlumptmp << "\tcerr << \"" << N.Place[j].label.substr(1,N.Place[j].label.length()-1 ) << " = \" << vect[reducePL_" << N.Place[j].label.substr(1,N.Place[j].label.length()-1 ) << "] << endl;" << endl;
499 : };
500 1 : outputlumptmp << "}" << endl;
501 :
502 1 : outputlumptmp << "bool SPN::precondition(const abstractMarking &Marking){return true;}" << endl;
503 :
504 1 : outputlumptmp << endl << "void SPN::lumpingFun(const abstractMarking &Marking,vector<int> &vect){" << endl;
505 3 : for(size_t i=0; i< N.Msimpletab.size();i++){
506 2 : int j = N.Msimpletab[i];
507 2 : outputlumptmp << "\tvect[reducePL_" << N.Place[j].label.substr(1,N.Place[j].label.length()-1 ) << "] = Marking.P->_PL_" << N.Place[j].label.substr(1,N.Place[j].label.length()-1 ) << "; //To Complete" << endl;
508 : };
509 1 : outputlumptmp << "}" << endl;
510 1 : outputlumptmp.close();
511 :
512 1 : }
513 :
514 0 : void stateSpace::outputVect(){
515 0 : cerr << "Exporting the probability vector" << endl;
516 :
517 0 : fstream outputFile;
518 0 : outputFile.open("muFile",fstream::out);
519 0 : outputFile.precision(15);
520 :
521 0 : outputFile << "[" << muvect->size() << "](";
522 0 : for (size_t i =0; i<muvect->size(); i++) {
523 0 : if(i>0)outputFile << ",";
524 0 : outputFile << (*muvect)[i];
525 : }
526 0 : outputFile << ")" << endl;
527 :
528 0 : for(hash_state::iterator it= S.begin() ; it != S.end(); it++){
529 0 : outputFile << "(";
530 0 : vector<int> vect = *(*it).first;
531 0 : for(size_t i=0; i< N.Msimpletab.size();i++){
532 0 : if(i>0)outputFile << ",";
533 0 : outputFile << vect[N.Msimpletab[i]];
534 : };
535 :
536 : /*for(int i =0; i< vect.size()-1; i++){
537 : if(i>0)outputFile << ",";
538 : outputFile << vect[i];
539 : }*/
540 0 : outputFile << ")=";
541 0 : outputFile << (*it).second << endl;
542 : }
543 :
544 0 : outputFile.close();
545 0 : }
546 :
547 0 : double stateSpace::returnPrismResult(){
548 0 : return (*muvect)[0];
549 : }
550 :
551 1 : void stateSpace::inputVect(){
552 1 : cerr<< "Start reading muFile" << endl;
553 2 : ifstream inputFile("muFile",fstream::in);
554 :
555 1 : if(!inputFile.good()){
556 0 : cerr << "Fail to open muFile at" <<endl;
557 0 : exit(EXIT_FAILURE);
558 : }
559 :
560 2 : boostmat::vector<double> v1;
561 1 : inputFile >> v1;
562 1 : nbState = v1.size();
563 1 : muvect = new vector<double>(nbState);
564 292 : for(size_t i=0; i< nbState; i++){
565 291 : (*muvect)[i] = v1 (i);
566 : }
567 :
568 2 : string line;
569 : size_t poseq;
570 2 : string pos;
571 2 : string stateid;
572 587 : while ( inputFile.good() )
573 : {
574 293 : getline (inputFile,line);
575 : //cerr << line << endl;
576 293 : poseq = line.find("=");
577 :
578 293 : if(poseq != string::npos){
579 291 : pos = line.substr(1,poseq-2);
580 291 : stateid = line.substr(poseq+1,line.size());
581 :
582 582 : vector<int> vect;
583 291 : size_t it = 0;
584 1455 : while(it < pos.length()){
585 582 : size_t it2 = pos.find(",",it);
586 582 : if(it2 == string::npos) it2 = pos.length();
587 582 : vect.push_back(atoi((pos.substr(it,it2-it)).c_str() ));
588 582 : it = it2+1;
589 : }
590 :
591 291 : S[new vector<int>(vect)] = (int)atoi(stateid.c_str());
592 :
593 : }
594 : }
595 :
596 1 : inputFile.close();
597 1 : if(S.empty()){
598 0 : cerr << "muFile empty" << endl;
599 0 : exit(EXIT_FAILURE);
600 : }
601 1 : cerr<< "Finished reading muFile with "<< nbState << " states" << endl;
602 1 : }
603 :
604 0 : void stateSpace::inputMat(){
605 0 : fstream inputFile("matrixFile",fstream::in);
606 :
607 0 : if(!inputFile.good()){
608 0 : cerr << "Fail to open matrixFile"<<endl;
609 0 : exit(EXIT_FAILURE);
610 : return;
611 : }
612 :
613 : /*boostmat::matrix<double> m1;
614 : inputFile >> m1;
615 : nbState = m1.size1();*/
616 0 : boostmat::compressed_matrix<double , boostmat::row_major> m;
617 0 : inputFile >> boostmat::io::sparse(m);
618 :
619 : /*for (unsigned i = 0; i < nbState; ++ i)
620 : for (unsigned j = 0; j < nbState; ++ j)
621 : if(m1 (i,j) != 0.) m (i, j) = m1 (i,j);*/
622 0 : transitionsMatrix = new boostmat::compressed_matrix<double>(m);
623 :
624 : //cerr << *transitionsMatrix << endl;
625 :
626 0 : finalVector = new boostmat::vector<double>();
627 0 : inputFile >> (*finalVector);
628 :
629 :
630 0 : string line;
631 : size_t poseq;
632 0 : string pos;
633 0 : string stateid;
634 0 : while ( inputFile.good() )
635 : {
636 0 : getline (inputFile,line);
637 : //cerr << line << endl;
638 0 : poseq = line.find("=");
639 :
640 0 : if(poseq != string::npos ){
641 0 : pos = line.substr(1,poseq-2);
642 0 : stateid = line.substr(poseq+1,line.size());
643 :
644 0 : vector<int> vect;
645 0 : size_t it = 0;
646 0 : while(it < pos.length()){
647 0 : size_t it2 = pos.find(",",it);
648 0 : if(it2 == string::npos) it2 = pos.length();
649 0 : vect.push_back(atoi((pos.substr(it,it2-it)).c_str() ));
650 0 : it = it2+1;
651 : }
652 :
653 0 : S[new vector<int>(vect)] = (int)atoi(stateid.c_str());
654 :
655 : }
656 : }
657 :
658 0 : inputFile.close();
659 :
660 0 : cerr << "DTMC size:" << finalVector->size() << endl;
661 :
662 12 : }
|
