From 2a1487dc390049fc438f7c1c11057d6e782e28a5 Mon Sep 17 00:00:00 2001
From: Timo Philipp Gros <timopgros@gmail.com>
Date: Wed, 29 Nov 2017 20:05:25 +0100
Subject: [PATCH] back to copied version of foxglynn, leaving too small values

---
 .../csl/helper/SparseMarkovAutomatonCslHelper.cpp     | 11 ++++++-----
 .../csl/helper/SparseMarkovAutomatonCslHelper.h       |  5 +++--
 2 files changed, 9 insertions(+), 7 deletions(-)

diff --git a/src/storm/modelchecker/csl/helper/SparseMarkovAutomatonCslHelper.cpp b/src/storm/modelchecker/csl/helper/SparseMarkovAutomatonCslHelper.cpp
index b4c7b0526..23a8accde 100644
--- a/src/storm/modelchecker/csl/helper/SparseMarkovAutomatonCslHelper.cpp
+++ b/src/storm/modelchecker/csl/helper/SparseMarkovAutomatonCslHelper.cpp
@@ -267,7 +267,7 @@ namespace storm {
                         // Vd
                         res = storm::utility::zero<ValueType>();
                         for (uint64_t i = k ; i<N ; i++){
-                            if (i < poisson.size()){
+                            if (poisson[i]>1e-300){
                                 ValueType between = poisson[i];
                                 res+=between;
                             }
@@ -652,6 +652,8 @@ namespace storm {
                     }
 
 
+                    std::vector<double> poisson = foxGlynnProb(lambda*T, N+1, epsilon*kappa);
+
                     // (4) define vectors/matrices
                     std::vector<ValueType> init(numberOfStates, -1);
                     vd = std::vector<std::vector<ValueType>>(N + 1, init);
@@ -669,13 +671,13 @@ namespace storm {
                         for (uint64_t k = N; k <= N; k--) {
                             calculateUnifPlusVector(k, i, 0, lambda, probSize, relReachability, dir, unifVectors,
                                                     fullTransitionMatrix, markovianStates, psiStates, solver, logfile,
-                                                    std::get<3>(foxGlynnResult));
+                                                    poisson);
                             calculateUnifPlusVector(k, i, 2, lambda, probSize, relReachability, dir, unifVectors,
                                                     fullTransitionMatrix, markovianStates, psiStates, solver, logfile,
-                                                    std::get<3>(foxGlynnResult));
+                                                    poisson);
                             calculateVu(relReachability, dir, k, i, 1, lambda, probSize, unifVectors,
                                         fullTransitionMatrix, markovianStates, psiStates, solver, logfile,
-                                        std::get<3>(foxGlynnResult));
+                                        poisson);
                             //also use iteration to keep maxNorm of vd and vup to date, so the loop-condition is easy to prove
                             ValueType diff = std::abs(unifVectors[0][k][i] - unifVectors[1][k][i]);
                             maxNorm = std::max(maxNorm, diff);
@@ -700,7 +702,6 @@ namespace storm {
 
                 logfile.close();
                 return unifVectors[0][0];
-
             }
 
             template <typename ValueType, typename std::enable_if<storm::NumberTraits<ValueType>::SupportsExponential, int>::type>
diff --git a/src/storm/modelchecker/csl/helper/SparseMarkovAutomatonCslHelper.h b/src/storm/modelchecker/csl/helper/SparseMarkovAutomatonCslHelper.h
index ab036c8c3..4ce22a17b 100644
--- a/src/storm/modelchecker/csl/helper/SparseMarkovAutomatonCslHelper.h
+++ b/src/storm/modelchecker/csl/helper/SparseMarkovAutomatonCslHelper.h
@@ -106,16 +106,17 @@ namespace storm {
                     std::vector<double> poisson_p(N,0.0);
                     unsigned long iter_num;
 
+                    std::cout << "fg left " << fg->left << " fh right " << fg->right <<"\n";
                     //for(int i=fg->left; i<=fg->right; i++) {
                     for (int i = 0; i<N ; i++){
                         poisson_p[i] = fg->weights[i-fg->left]/fg->total_weight;
                         sumOfPoissonProbs+=poisson_p[i];
                     }
 
-                    /*for(int i=fg->left-1; i>=0; i--) {
+                    for(int i=fg->left-1; i>=0; i--) {
                         poisson_p[i] = poisson_p[i+1]*((i+1)/(lambdaT));
                         sumOfPoissonProbs+=poisson_p[i];
-                    }*/
+                    }
 
                     iter_num = fg->right;