From 25a7b6c71a652b067af212cd2cecde7e1748e1ee Mon Sep 17 00:00:00 2001
From: Timo Philipp Gros <timopgros@gmail.com>
Date: Sun, 19 Nov 2017 12:20:54 +0100
Subject: [PATCH] implemented trajans alg to identify prob Cycles

---
 .../helper/SparseMarkovAutomatonCslHelper.cpp | 74 ++++++++++++++++++-
 .../helper/SparseMarkovAutomatonCslHelper.h   | 15 ++++
 2 files changed, 88 insertions(+), 1 deletion(-)

diff --git a/src/storm/modelchecker/csl/helper/SparseMarkovAutomatonCslHelper.cpp b/src/storm/modelchecker/csl/helper/SparseMarkovAutomatonCslHelper.cpp
index 8c3d33fc2..88d470dce 100644
--- a/src/storm/modelchecker/csl/helper/SparseMarkovAutomatonCslHelper.cpp
+++ b/src/storm/modelchecker/csl/helper/SparseMarkovAutomatonCslHelper.cpp
@@ -337,10 +337,82 @@ namespace storm {
                 logfile << " res = " << res << "\n";
             }
 
+            template <typename ValueType, typename std::enable_if<storm::NumberTraits<ValueType>::SupportsExponential, int>::type=0>
+            int SparseMarkovAutomatonCslHelper::trajans(storm::storage::SparseMatrix<ValueType> const& transitionMatrix, uint64_t node, std::vector<uint64_t >& disc, std::vector<uint64_t >& finish, uint64_t* counter) {
+                auto const &rowGroupIndice = transitionMatrix.getRowGroupIndices();
+
+                disc[node] = *counter;
+                finish[node] = *counter;
+                (*counter)+=1;
+
+                auto from = rowGroupIndice[node];
+                auto to = rowGroupIndice[node+1];
+
+                for(uint64_t i =from; i<to ; i++ ) {
+                    for(auto element : transitionMatrix.getRow(i)){
+                        if (disc[element.getColumn()]==0){
+                            uint64_t back = trajans(transitionMatrix,element.getColumn(),disc,finish, counter);
+                            finish[node]=std::min(finish[node], back);
+                        } else {
+                            finish[node]=std::min(finish[node], disc[element.getColumn()]);
+                        }
+                    }
+                }
+
+                return finish[node];
+            }
+
+
+            template <typename ValueType, typename std::enable_if<storm::NumberTraits<ValueType>::SupportsExponential, int>::type=0>
+            storm::storage::BitVector SparseMarkovAutomatonCslHelper::identifyProbCycles(storm::storage::SparseMatrix<ValueType> const& transitionMatrix, storm::storage::BitVector const& markovianStates, storm::storage::BitVector const& psiStates){
+
+                storm::storage::BitVector const& probabilisticStates = ~markovianStates;
+                storm::storage::BitVector const& probabilisticNonGoalStates = ~markovianStates & ~psiStates;
+
+                storm::storage::SparseMatrix<ValueType> const& probMatrix = transitionMatrix.getSubmatrix(true, probabilisticNonGoalStates, probabilisticNonGoalStates);
+                uint64_t  probSize = probMatrix.getRowGroupCount();
+                std::vector<uint64_t> disc(probSize, 0), finish(probSize, 0);
+
+                uint64_t counter =1;
+
+                for (uint64_t i =0; i<probSize; i++){
+                    if (disc[i]==0) {
+                        trajans(probMatrix, i, disc, finish, &counter);
+                    }
+                }
+
+                std::cout << "giving the circle vectors \n";
+                for (uint64_t i = 0 ; i<probSize; i++){
+                    std::cout << disc[i] << '\t' << finish[i] << "\n";
+                }
+
+
+                storm::storage::BitVector cycleStates(markovianStates.size(), false);
+                for (int i = 0 ; i< finish.size() ; i++){
+                    auto f = finish[i];
+                    for (int j =i; j<finish.size() ; j++){
+                        if (finish[j]==f){
+
+                            cycleStates.set(transformIndice(probabilisticNonGoalStates,i),true);
+                            cycleStates.set(transformIndice(probabilisticNonGoalStates,j),true);
+                        }
+                    }
+                }
+
+                return cycleStates;
+            }
+
+
             template <typename ValueType, typename std::enable_if<storm::NumberTraits<ValueType>::SupportsExponential, int>::type>
             std::vector<ValueType> SparseMarkovAutomatonCslHelper::unifPlus( std::pair<double, double> const& boundsPair, std::vector<ValueType> const& exitRateVector, storm::storage::SparseMatrix<ValueType> const& transitionMatrix, storm::storage::BitVector const& markovianStates, storm::storage::BitVector const& psiStates){
                 STORM_LOG_TRACE("Using UnifPlus to compute bounded until probabilities.");
-                
+
+                storm::storage::BitVector const& probabilisticNonGoalStates = ~markovianStates & ~psiStates;
+                auto cycleStates = identifyProbCycles(transitionMatrix, markovianStates, psiStates);
+
+                for (int i=0 ; i < cycleStates.size() ; i++ ){
+                    std::cout << cycleStates[i] <<  "\t" << probabilisticNonGoalStates[i] <<"\n";
+                }
                 std::ofstream logfile("U+logfile.txt", std::ios::app);
                 ValueType maxNorm = 0;
 
diff --git a/src/storm/modelchecker/csl/helper/SparseMarkovAutomatonCslHelper.h b/src/storm/modelchecker/csl/helper/SparseMarkovAutomatonCslHelper.h
index ff27180bd..559912c2e 100644
--- a/src/storm/modelchecker/csl/helper/SparseMarkovAutomatonCslHelper.h
+++ b/src/storm/modelchecker/csl/helper/SparseMarkovAutomatonCslHelper.h
@@ -56,7 +56,20 @@ namespace storm {
                 static std::vector<ValueType> computeReachabilityTimes(OptimizationDirection dir, storm::storage::SparseMatrix<ValueType> const& transitionMatrix, storm::storage::SparseMatrix<ValueType> const& backwardTransitions, std::vector<ValueType> const& exitRateVector, storm::storage::BitVector const& markovianStates, storm::storage::BitVector const& psiStates, storm::solver::MinMaxLinearEquationSolverFactory<ValueType> const& minMaxLinearEquationSolverFactory);
                 
             private:
+                static uint64_t transformIndice(storm::storage::BitVector const& subset, uint64_t fakeId){
+                    uint64_t id =0;
+                    uint64_t counter =0;
+                    while(counter<=fakeId){
+                        if(subset[id]){
+                            counter++;
+                        }
+                        id++;
+                    }
+                    return id-1;
+                }
 
+                template <typename ValueType, typename std::enable_if<storm::NumberTraits<ValueType>::SupportsExponential, int>::type=0>
+                static storm::storage::BitVector identifyProbCycles(storm::storage::SparseMatrix<ValueType> const& TransitionMatrix, storm::storage::BitVector const& markovianStates, storm::storage::BitVector const& psiStates);
                 /*!
                  * Computes the poission-distribution
                  *
@@ -69,6 +82,8 @@ namespace storm {
                 template <typename ValueType>
                 static ValueType poisson(ValueType lambda, uint64_t i);
 
+                template <typename ValueType, typename std::enable_if<storm::NumberTraits<ValueType>::SupportsExponential, int>::type=0>
+                static int trajans(storm::storage::SparseMatrix<ValueType> const& TransitionMatrix, uint64_t node, std::vector<uint64_t>& disc, std::vector<uint64_t>& finish, uint64_t * counter);
 
                 /*!
                  * Computes vd vector according to UnifPlus