Currently debugging the computation of transient probabilities in CTMCs.

Former-commit-id: 6671e0205d

examples/ctmc/cluster/cluster.sm

@@ -0,0 +1,116 @@
+// Workstation cluster [HHK00]
+// dxp/gxn 11/01/00
+
+ctmc
+
+const int N; // Number of workstations in each cluster
+const int left_mx = N; // Number of work stations in left cluster
+const int right_mx = N; // Number of work stations in right cluster
+
+// Failure rates
+const double ws_fail = 1/500; // Single workstation: average time to fail = 500 hrs
+const double switch_fail = 1/4000; // Switch: average time to fail = 4000 hrs
+const double line_fail = 1/5000; // Backbone: average time to fail = 5000 hrs
+
+// Left cluster
+module Left 
+	
+	left_n : [0..left_mx] init left_mx; // Number of workstations operational
+	left : bool; // Being repaired?
+	
+	[startLeft] !left & (left_n<left_mx) -> 1 : (left'=true);
+	[repairLeft] left & (left_n<left_mx) -> 1 : (left'=false) & (left_n'=left_n+1);
+	[] (left_n>0) -> ws_fail*left_n : (left_n'=left_n-1);
+	
+endmodule
+
+// Right cluster
+module Right = Left[left_n=right_n,
+                    left=right,
+                    left_mx=right_mx,
+                    startLeft=startRight,
+                    repairLeft=repairRight ]
+endmodule
+
+// Repair unit
+module Repairman
+	
+	r : bool; // Repairing?
+	
+	[startLeft]    !r -> 10 : (r'=true); // Inspect Left 
+	[startRight]   !r -> 10 : (r'=true); // Inspect Right 
+	[startToLeft]  !r -> 10 : (r'=true); // Inspect ToLeft
+	[startToRight] !r -> 10 : (r'=true); // Inspect ToRight 
+	[startLine]    !r -> 10 : (r'=true); // Inspect Line 
+	
+	[repairLeft]    r -> 2     : (r'=false); // Repair Left 
+	[repairRight]   r -> 2     : (r'=false); // Repair Right
+	[repairToLeft]  r -> 0.25  : (r'=false); // Repair ToLeft
+	[repairToRight] r -> 0.25  : (r'=false); // Repair ToRight
+	[repairLine]    r -> 0.125 : (r'=false); // Repair Line
+	
+endmodule
+
+// Line/backbone
+module Line 
+	
+	line :   bool; // Being repaired?
+	line_n : bool init true; // Working?
+	
+	[startLine] !line & !line_n -> 1 : (line'=true);
+	[repairLine] line & !line_n -> 1 : (line'=false) & (line_n'=true);
+	[] line_n -> line_fail : (line_n'=false);
+	
+endmodule
+
+// Left switch
+module ToLeft = Line[line=toleft,
+                     line_n=toleft_n,
+                     line_fail=switch_fail,
+                     startLine=startToLeft,
+                     repairLine=repairToLeft ]
+endmodule
+
+// Right switch
+module ToRight = Line[line=toright,
+                      line_n=toright_n,
+                      line_fail=switch_fail,
+                      startLine=startToRight,
+                      repairLine=repairToRight ]
+endmodule
+
+// Formulas + labels
+
+// Minimum QoS requires 3/4 connected workstations operational
+const int k = floor(0.75*N);
+// left_operational_i : left_n>=i & toleft_n
+// right_operational_i : right_n>=i & toright_n
+// operational_i : (left_n+right_n)>=i & toleft_n & line_n & toright_n
+// minimum_k : left_operational_k | right_operational_k | operational_k
+formula minimum = (left_n>=k & toleft_n) | 
+                  (right_n>=k & toright_n) | 
+                  ((left_n+right_n)>=k & toleft_n & line_n & toright_n);
+label "minimum" = (left_n>=k & toleft_n) | (right_n>=k & toright_n) | ((left_n+right_n)>=k & toleft_n & line_n & toright_n);
+// premium = minimum_N
+label "premium" = (left_n>=left_mx & toleft_n) | (right_n>=right_mx & toright_n) | ((left_n+right_n)>=left_mx & toleft_n & line_n & toright_n);
+
+// Reward structures
+
+// Percentage of operational workstations stations
+rewards "percent_op"
+	true : 100*(left_n+right_n)/(2*N);
+endrewards
+
+// Time that the system is not delivering at least minimum QoS
+rewards "time_not_min"
+	!minimum : 1; 
+endrewards
+
+// Number of repairs
+rewards "num_repairs"
+	[repairLeft]    true : 1;
+	[repairRight]   true : 1;
+	[repairToLeft]  true : 1;
+	[repairToRight] true : 1;
+	[repairLine]    true : 1;
+endrewards

examples/ctmc/embedded/embedded.sm

@@ -0,0 +1,151 @@
+ctmc
+
+// constants
+const int MAX_COUNT;
+const int MIN_SENSORS = 2;
+const int MIN_ACTUATORS = 1;
+
+// rates
+const double lambda_p = 1/(365*24*60*60); // 1 year
+const double lambda_s = 1/(30*24*60*60); // 1 month
+const double lambda_a = 1/(2*30*24*60*60); // 2 months
+const double tau = 1/60; // 1 min
+const double delta_f = 1/(24*60*60); // 1 day
+const double delta_r = 1/30; // 30 secs
+
+// sensors
+module sensors
+	
+	s : [0..3] init 3; // number of sensors working
+	
+	[] s>1 -> s*lambda_s : (s'=s-1); // failure of a single sensor
+	
+endmodule
+
+// input processor
+// (takes data from sensors and passes onto main processor)
+module proci
+	
+	i : [0..2] init 2; // 2=ok, 1=transient fault, 0=failed
+	
+	[] i>0 & s>=MIN_SENSORS -> lambda_p : (i'=0); // failure of processor
+	[] i=2 & s>=MIN_SENSORS -> delta_f : (i'=1); // transient fault
+	[input_reboot] i=1 & s>=MIN_SENSORS -> delta_r : (i'=2); // reboot after transient fault
+	
+endmodule
+
+// actuators
+module actuators
+	
+	a : [0..2] init 2; // number of actuators working
+	
+	[] a>0 -> a*lambda_a : (a'=a-1); // failure of a single actuator
+	
+endmodule
+
+// output processor
+// (receives instructions from main processor and passes onto actuators)
+module proco = proci [ i=o, s=a, input_reboot=output_reboot, MIN_SENSORS=MIN_ACTUATORS ] endmodule
+
+// main processor
+// (takes data from proci, processes it, and passes instructions to proco)
+module procm
+	
+	m : [0..1] init 1; // 1=ok, 0=failed
+	count : [0..MAX_COUNT+1] init 0; // number of consecutive skipped cycles
+	
+	// failure of processor
+	[] m=1 -> lambda_p : (m'=0);
+	// processing completed before timer expires - reset skipped cycle counter
+	[timeout]  comp -> tau : (count'=0);
+	// processing not completed before timer expires - increment skipped cycle counter
+	[timeout] !comp -> tau : (count'=min(count+1, MAX_COUNT+1));
+	
+endmodule
+
+// connecting bus
+module bus
+	
+	// flags
+	// main processor has processed data from input processor
+	// and sent corresponding instructions to output processor (since last timeout)
+	comp : bool init true; 
+	// input processor has data ready to send
+	reqi : bool init true; 
+	// output processor has instructions ready to be processed
+	reqo : bool init false;
+	
+	// input processor reboots
+	[input_reboot]  true -> 1 :
+	// performs a computation if has already done so or
+	// it is up and ouput clear (i.e. nothing waiting)
+	(comp'=(comp | (m=1 & !reqo))) 
+	// up therefore something to process
+	& (reqi'=true)
+	// something to process if not functioning and either
+	// there is something already pending
+	// or the main processor sends a request
+	& (reqo'=!(o=2 & a>=1) & (reqo | m=1));
+	
+	// output processor reboots
+	[output_reboot] true -> 1 :
+	// performs a computation if it has already or
+	// something waiting and is up
+	// (can be processes as the output has come up and cleared pending requests)
+	(comp'=(comp | (reqi & m=1)))
+	// something to process it they are up or
+	// there was already something and the main processor acts
+	// (output now up must be due to main processor being down)
+	& (reqi'=(i=2 & s>=2) | (reqi & m=0))
+	// output and actuators up therefore nothing can be pending
+	& (reqo'=false);
+	
+	// main processor times out
+	[timeout] true -> 1 :
+	// performs a computation if it is up something was pending
+	// and nothing is waiting for the output
+	(comp'=(reqi & !reqo & m=1))
+	// something to process if up or
+	// already something and main process cannot act 
+	// (down or outputs pending)
+	& (reqi'=(i=2 & s>=2) | (reqi & (reqo | m=0)))
+	// something to process if they are not functioning and 
+	// either something is already pending
+	// or the main processor acts
+	& (reqo'=!(o=2 & a>=1) & (reqo | (reqi & m=1)));
+	
+endmodule
+
+
+// the system is down
+formula down = (i=2&s<MIN_SENSORS)|(count=MAX_COUNT+1)|(o=2&a<MIN_ACTUATORS)|(m=0);
+// transient failure has occured but the system is not down
+formula danger = !down & (i=1 | o=1);
+// the system is operational
+formula up = !down & !danger;
+
+
+// reward structures
+
+rewards "up"
+	up : 1/3600;
+endrewards
+
+rewards "danger"
+	danger : 1/3600;
+endrewards
+rewards "down"
+	down : 1/3600;
+endrewards
+
+//labels
+// causes of failues
+label "fail_sensors" = i=2&s<MIN_SENSORS; // sensors have failed
+label "fail_actuators" = o=2&a<MIN_ACTUATORS; // actuators have failed
+label "fail_io" = count=MAX_COUNT+1; // IO has failed
+label "fail_main" = m=0; // ,main processor has failed
+
+// system status
+label "down" = (i=2&s<MIN_SENSORS)|(count=MAX_COUNT+1)|(o=2&a<MIN_ACTUATORS)|(m=0); // system has shutdown
+label "danger" = !down & (i=1 | o=1); // transient fault has occured
+label "up" = !down & !danger; 

examples/ctmc/embedded/embedded_debug.sm

@@ -0,0 +1,33 @@
+ctmc
+
+// constants
+const int MAX_COUNT;
+const int MIN_SENSORS = 2;
+const int MIN_ACTUATORS = 1;
+
+// rates
+const double lambda_p = 1/(365*24*60*60); // 1 year
+const double lambda_s = 1/(30*24*60*60); // 1 month
+const double lambda_a = 1/(2*30*24*60*60); // 2 months
+const double tau = 1/60; // 1 min
+const double delta_f = 1/(24*60*60); // 1 day
+const double delta_r = 1/30; // 30 secs
+
+// sensors
+module sensors
+	
+	s : [0..3] init 3; // number of sensors working
+	
+	[] s>1 -> s*lambda_s : (s'=s-1); // failure of a single sensor
+	
+endmodule
+
+// input processor
+// (takes data from sensors and passes onto main processor)
+module proci
+	
+	i : [0..2] init 2; // 2=ok, 1=transient fault, 0=failed
+	
+	[] i>0 & s>=MIN_SENSORS -> lambda_p : (i'=0); // failure of processor
+	
+endmodule

examples/ctmc/polling/polling2.sm

@@ -0,0 +1,51 @@
+// polling example [IT90]
+// gxn/dxp 26/01/00
+
+ctmc
+
+const int N = 2;
+
+const double mu		= 1;
+const double gamma	= 200;
+const double lambda	= mu/N;
+
+module server
+	
+	s : [1..2]; // station
+	a : [0..1]; // action: 0=polling, 1=serving
+	
+	[loop1a] (s=1)&(a=0) -> gamma	: (s'=s+1);
+	[loop1b] (s=1)&(a=0) -> gamma	: (a'=1);
+	[serve1] (s=1)&(a=1) -> mu		: (s'=s+1)&(a'=0);
+	
+	[loop2a] (s=2)&(a=0) -> gamma	: (s'=1);
+	[loop2b] (s=2)&(a=0) -> gamma	: (a'=1);
+	[serve2] (s=2)&(a=1) -> mu		: (s'=1)&(a'=0);
+	
+endmodule
+
+module station1
+	
+	s1 : [0..1]; // state of station: 0=empty, 1=full
+	
+	[loop1a] (s1=0) -> 1 : (s1'=0);
+	[]       (s1=0) -> lambda : (s1'=1);
+	[loop1b] (s1=1) -> 1 : (s1'=1);
+	[serve1] (s1=1) -> 1 : (s1'=0);
+	
+endmodule
+
+// construct further stations through renaming
+
+module station2 = station1 [ s1=s2, loop1a=loop2a, loop1b=loop2b, serve1=serve2 ] endmodule
+// (cumulative) rewards
+
+// expected time station 1 is waiting to be served
+rewards "waiting"
+	s1=1 & !(s=1 & a=1) : 1;
+endrewards
+
+// expected number of times station 1 is served
+rewards "served"
+	[serve1] true : 1;
+endrewards

examples/ctmc/polling/polling5.sm

@@ -0,0 +1,66 @@
+// polling example [IT90]
+// gxn/dxp 26/01/00
+
+ctmc
+
+const int N = 5;
+
+const double mu		= 1;
+const double gamma	= 200;
+const double lambda	= mu/N;
+
+module server
+	
+	s : [1..5]; // station
+	a : [0..1]; // action: 0=polling, 1=serving
+	
+	[loop1a] (s=1)&(a=0) -> gamma	: (s'=s+1);
+	[loop1b] (s=1)&(a=0) -> gamma	: (a'=1);
+	[serve1] (s=1)&(a=1) -> mu		: (s'=s+1)&(a'=0);
+	
+	[loop2a] (s=2)&(a=0) -> gamma	: (s'=s+1);
+	[loop2b] (s=2)&(a=0) -> gamma	: (a'=1);
+	[serve2] (s=2)&(a=1) -> mu		: (s'=s+1)&(a'=0);
+	
+	[loop3a] (s=3)&(a=0) -> gamma	: (s'=s+1);
+	[loop3b] (s=3)&(a=0) -> gamma	: (a'=1);
+	[serve3] (s=3)&(a=1) -> mu		: (s'=s+1)&(a'=0);
+	
+	[loop4a] (s=4)&(a=0) -> gamma	: (s'=s+1);
+	[loop4b] (s=4)&(a=0) -> gamma	: (a'=1);
+	[serve4] (s=4)&(a=1) -> mu		: (s'=s+1)&(a'=0);
+	
+	[loop5a] (s=5)&(a=0) -> gamma	: (s'=1);
+	[loop5b] (s=5)&(a=0) -> gamma	: (a'=1);
+	[serve5] (s=5)&(a=1) -> mu		: (s'=1)&(a'=0);
+	
+endmodule
+
+module station1
+	
+	s1 : [0..1]; // state of station: 0=empty, 1=full
+	
+	[loop1a] (s1=0) -> 1 : (s1'=0);
+	[]       (s1=0) -> lambda : (s1'=1);
+	[loop1b] (s1=1) -> 1 : (s1'=1);
+	[serve1] (s1=1) -> 1 : (s1'=0);
+	
+endmodule
+
+// construct further stations through renaming
+
+module station2 = station1 [ s1=s2, loop1a=loop2a, loop1b=loop2b, serve1=serve2 ] endmodule
+module station3 = station1 [ s1=s3, loop1a=loop3a, loop1b=loop3b, serve1=serve3 ] endmodule
+module station4 = station1 [ s1=s4, loop1a=loop4a, loop1b=loop4b, serve1=serve4 ] endmodule
+module station5 = station1 [ s1=s5, loop1a=loop5a, loop1b=loop5b, serve1=serve5 ] endmodule
+// (cumulative) rewards
+
+// expected time station 1 is waiting to be served
+rewards "waiting"
+	s1=1 & !(s=1 & a=1) : 1;
+endrewards
+
+// expected number of times station 1 is served
+rewards "served"
+	[serve1] true : 1;
+endrewards

examples/ctmc/tiny/tiny.sm

@@ -1,11 +1,11 @@
 ctmc
 
 module one
-  s : [0 .. 3] init 1;
+  s : [0 .. 3] init 0;
   
   [] s<3 -> 3/2 : (s'=s+1);
   [] s>0 -> 3 : (s'=s-1);
 endmodule
 
 label "empty" = s=0;
-label "full" = s=3;
+label "full" = s=3;

src/logic/BinaryPathFormula.cpp

@@ -14,6 +14,10 @@ namespace storm {
             return this->getLeftSubformula().isPctlStateFormula() && this->getRightSubformula().isPctlStateFormula();
         }
         
+        bool BinaryPathFormula::isCslPathFormula() const {
+            return this->getLeftSubformula().isCslStateFormula() && this->getRightSubformula().isCslStateFormula();
+        }
+        
         bool BinaryPathFormula::isLtlFormula() const {
             return this->getLeftSubformula().isLtlFormula() && this->getRightSubformula().isLtlFormula();
         }

src/logic/BinaryPathFormula.h

@@ -18,6 +18,7 @@ namespace storm {
             virtual bool isBinaryPathFormula() const override;
             
             virtual bool isPctlPathFormula() const override;
+            virtual bool isCslPathFormula() const override;
             virtual bool isLtlFormula() const override;
             virtual bool containsBoundedUntilFormula() const override;
             virtual bool containsNextFormula() const override;

src/logic/BinaryStateFormula.cpp

@@ -13,6 +13,10 @@ namespace storm {
         bool BinaryStateFormula::isPctlStateFormula() const {
             return this->getLeftSubformula().isPctlStateFormula() && this->getRightSubformula().isPctlStateFormula();
         }
+
+        bool BinaryStateFormula::isCslStateFormula() const {
+            return this->getLeftSubformula().isCslStateFormula() && this->getRightSubformula().isCslStateFormula();
+        }
         
         bool BinaryStateFormula::isLtlFormula() const {
             return this->getLeftSubformula().isLtlFormula() && this->getRightSubformula().isLtlFormula();

src/logic/BinaryStateFormula.h

@@ -16,6 +16,7 @@ namespace storm {
             virtual bool isBinaryStateFormula() const override;
 
             virtual bool isPctlStateFormula() const override;
+            virtual bool isCslStateFormula() const override;
             virtual bool isLtlFormula() const override;
             virtual bool containsBoundedUntilFormula() const override;
             virtual bool containsNextFormula() const override;

src/logic/ProbabilityOperatorFormula.cpp

@@ -26,6 +26,10 @@ namespace storm {
             return this->getSubformula().isPctlPathFormula();
         }
         
+        bool ProbabilityOperatorFormula::isCslStateFormula() const {
+            return this->getSubformula().isCslPathFormula();
+        }
+        
         bool ProbabilityOperatorFormula::isPltlFormula() const {
             return this->getSubformula().isLtlFormula();
         }

src/logic/ProbabilityOperatorFormula.h

@@ -18,6 +18,7 @@ namespace storm {
             }
             
             virtual bool isPctlStateFormula() const override;
+            virtual bool isCslStateFormula() const override;
             virtual bool isPltlFormula() const override;
             virtual bool containsProbabilityOperator() const override;
             virtual bool containsNestedProbabilityOperators() const override;

src/modelchecker/csl/SparseCtmcCslModelChecker.cpp

@@ -36,14 +36,23 @@ namespace storm {
         
         template<class ValueType>
         std::unique_ptr<CheckResult> SparseCtmcCslModelChecker<ValueType>::computeBoundedUntilProbabilities(storm::logic::BoundedUntilFormula const& pathFormula, bool qualitative, boost::optional<storm::logic::OptimalityType> const& optimalityType) {
-            STORM_LOG_THROW(pathFormula.isIntervalBounded(), storm::exceptions::InvalidPropertyException, "Cannot treat non-interval bounded until.");
-            
             std::unique_ptr<CheckResult> leftResultPointer = this->check(pathFormula.getLeftSubformula());
             std::unique_ptr<CheckResult> rightResultPointer = this->check(pathFormula.getRightSubformula());
             ExplicitQualitativeCheckResult const& leftResult = leftResultPointer->asExplicitQualitativeCheckResult();;
             ExplicitQualitativeCheckResult const& rightResult = rightResultPointer->asExplicitQualitativeCheckResult();
-            std::pair<double, double> const& intervalBounds =  pathFormula.getIntervalBounds();
-            std::unique_ptr<CheckResult> result = std::unique_ptr<CheckResult>(new ExplicitQuantitativeCheckResult<ValueType>(this->computeBoundedUntilProbabilitiesHelper(leftResult.getTruthValuesVector(), rightResult.getTruthValuesVector(), this->getModel().getExitRateVector(), qualitative, intervalBounds.first, intervalBounds.second)));
+            double lowerBound = 0;
+            double upperBound = 0;
+            if (pathFormula.isIntervalBounded()) {
+                std::pair<double, double> const& intervalBounds =  pathFormula.getIntervalBounds();
+                lowerBound = intervalBounds.first;
+                upperBound = intervalBounds.second;
+            } else {
+                upperBound = pathFormula.getUpperBound();
+            }
+            
+            std::cout << "initial: " << this->getModel().getInitialStates() << std::endl;
+            
+            std::unique_ptr<CheckResult> result = std::unique_ptr<CheckResult>(new ExplicitQuantitativeCheckResult<ValueType>(this->computeBoundedUntilProbabilitiesHelper(leftResult.getTruthValuesVector(), rightResult.getTruthValuesVector(), this->getModel().getExitRateVector(), qualitative, lowerBound, upperBound)));
             
             return result;
         }
@@ -98,6 +107,7 @@ namespace storm {
                     storm::utility::vector::setVectorValues<ValueType>(result, psiStates, storm::utility::one<ValueType>());
                 } else {
                     if (comparator.isZero(lowerBound)) {
+                        std::cout << "case [0, t]" << std::endl;
                         // In this case, the interval is of the form [0, t].
                         // Note that this excludes [0, inf] since this is untimed reachability and we considered this case earlier.
 
@@ -106,41 +116,57 @@ namespace storm {
                         for (auto const& state : statesWithProbabilityGreater0NonPsi) {
                             uniformizationRate = std::max(uniformizationRate, exitRates[state]);
                         }
+                        uniformizationRate *= 1.02;
                         STORM_LOG_THROW(uniformizationRate > 0, storm::exceptions::InvalidStateException, "The uniformization rate must be positive.");
                         
                         // Compute the uniformized matrix.
-                        storm::storage::SparseMatrix<ValueType> uniformizedMatrix = this->computeUniformizedMatrix(this->getModel().getTransitionMatrix(), storm::storage::BitVector(this->getModel().getNumberOfStates(), true), psiStates, uniformizationRate, exitRates);
+                        storm::storage::SparseMatrix<ValueType> uniformizedMatrix = this->computeUniformizedMatrix(this->getModel().getTransitionMatrix(), statesWithProbabilityGreater0, psiStates, uniformizationRate, exitRates);
+//                        storm::storage::SparseMatrix<ValueType> uniformizedMatrix = this->computeUniformizedMatrix(this->getModel().getTransitionMatrix(), storm::storage::BitVector(this->getModel().getNumberOfStates(), true), psiStates, uniformizationRate, exitRates);
 
                         // Finally compute the transient probabilities.
                         std::vector<ValueType> psiStateValues(statesWithProbabilityGreater0.getNumberOfSetBits(), storm::utility::zero<ValueType>());
                         storm::utility::vector::setVectorValues(psiStateValues, psiStates % statesWithProbabilityGreater0, storm::utility::one<ValueType>());
+//                        storm::utility::vector::setVectorValues(psiStateValues, psiStates, storm::utility::one<ValueType>());
+                        
                         std::vector<ValueType> subresult = this->computeTransientProbabilities(uniformizedMatrix, uniformizationRate * upperBound, psiStateValues, *this->linearEquationSolver);
                         storm::utility::vector::setVectorValues(result, statesWithProbabilityGreater0, subresult);
+//                        result = this->computeTransientProbabilities(uniformizedMatrix, uniformizationRate * upperBound, psiStateValues, *this->linearEquationSolver);
                     } else if (comparator.isInfinity(upperBound)) {
+                        std::cout << "case [t, inf]" << std::endl;
                         // In this case, the interval is of the form [t, inf] with t != 0.
                         
                         // Start by computing the (unbounded) reachability probabilities of reaching psi states while
                         // staying in phi states.
                         result = this->computeUntilProbabilitiesHelper(this->getModel().getTransitionMatrix(), backwardTransitions, phiStates, psiStates, qualitative, *this->linearEquationSolver);
                         
+                        storm::storage::BitVector absorbingStates = ~phiStates;
                         ValueType uniformizationRate = 0;
-                        for (auto const& state : statesWithProbabilityGreater0) {
+                        for (auto const& state : ~absorbingStates) {
                             uniformizationRate = std::max(uniformizationRate, exitRates[state]);
                         }
+                        uniformizationRate *= 1.02;
+                        STORM_LOG_THROW(uniformizationRate > 0, storm::exceptions::InvalidStateException, "The uniformization rate must be positive.");
+                        
+                        // Set the result to zero for all states that are known to violate phi.
+                        storm::utility::vector::setVectorValues(result, absorbingStates, storm::utility::zero<ValueType>());
+                        
+                        // FIXME: optimization: just consider the states that can reach a state with positive entry in the result.
                         
                         // Compute the uniformized matrix.
-                        storm::storage::SparseMatrix<ValueType> uniformizedMatrix = this->computeUniformizedMatrix(this->getModel().getTransitionMatrix(), statesWithProbabilityGreater0, storm::storage::BitVector(statesWithProbabilityGreater0.getNumberOfSetBits()), uniformizationRate, exitRates);
+                        storm::storage::SparseMatrix<ValueType> uniformizedMatrix = this->computeUniformizedMatrix(this->getModel().getTransitionMatrix(), storm::storage::BitVector(this->getModel().getNumberOfStates(), true), absorbingStates, uniformizationRate, exitRates);
                         
                         // Finally compute the transient probabilities.
                         result = this->computeTransientProbabilities(uniformizedMatrix, uniformizationRate * lowerBound, result, *this->linearEquationSolver);
                     } else {
                         // In this case, the interval is of the form [t, t'] with t != 0 and t' != inf.
+                        std::cout << "case [t, t']" << std::endl;
                         
                         // Find the maximal rate of all 'maybe' states to take it as the uniformization rate.
                         ValueType uniformizationRate = 0;
                         for (auto const& state : statesWithProbabilityGreater0NonPsi) {
                             uniformizationRate = std::max(uniformizationRate, exitRates[state]);
                         }
+                        uniformizationRate *= 1.02;
                         STORM_LOG_THROW(uniformizationRate > 0, storm::exceptions::InvalidStateException, "The uniformization rate must be positive.");
                         
                         // Compute the (first) uniformized matrix.
@@ -151,7 +177,7 @@ namespace storm {
                         storm::utility::vector::setVectorValues(psiStateValues, psiStates % statesWithProbabilityGreater0, storm::utility::one<ValueType>());
                         std::vector<ValueType> subresult = this->computeTransientProbabilities(uniformizedMatrix, uniformizationRate * (upperBound - lowerBound), psiStateValues, *this->linearEquationSolver);
                         storm::utility::vector::setVectorValues(result, statesWithProbabilityGreater0, subresult);
-
+                        
                         // Then compute the transient probabilities of being in such a state after t time units. For this,
                         // we must re-uniformize the CTMC, so we need to compute the second uniformized matrix.
                         storm::storage::BitVector absorbingStates = ~phiStates;
@@ -159,6 +185,12 @@ namespace storm {
                         for (auto const& state : ~absorbingStates) {
                             uniformizationRate = std::max(uniformizationRate, exitRates[state]);
                         }
+                        uniformizationRate *= 1.02;
+
+                        // Set the result to zero for all states that are known to violate phi.
+                        storm::utility::vector::setVectorValues(result, absorbingStates, storm::utility::zero<ValueType>());
+                        
+                        // FIXME: optimization: just consider the states that can reach a state with positive entry in the result.
                         
                         // Finally, we compute the second set of transient probabilities.
                         uniformizedMatrix = this->computeUniformizedMatrix(this->getModel().getTransitionMatrix(), storm::storage::BitVector(this->getModel().getNumberOfStates(), true), absorbingStates, uniformizationRate, exitRates);
@@ -167,6 +199,11 @@ namespace storm {
                 }
             }
             
+            std::cout << "final result" << std::endl;
+            for (uint_fast64_t index = 0; index < result.size(); ++index) {
+                std::cout << "result[" << index << "] = " << result[index] << std::endl;
+            }
+            
             return result;
         }
         
@@ -177,6 +214,7 @@ namespace storm {
             storm::storage::SparseMatrix<ValueType> uniformizedMatrix = transitionMatrix.getSubmatrix(false, maybeStates, maybeStates, true);
             
             // Make the appropriate states absorbing.
+            // FIXME: optimization: do not make absorbing, but rather add a fixed vector. This also prevents the "wrong" 0.999... for target states.
             uniformizedMatrix.makeRowsAbsorbing(absorbingStates % maybeStates);
             
             // Now we need to perform the actual uniformization. That is, all entries need to be divided by
@@ -189,7 +227,7 @@ namespace storm {
                         if (absorbingStates.get(state)) {
                             // Nothing to do here, since the state has already been made absorbing.
                         } else {
-                            element.setValue(-(exitRates[state] + element.getValue()) / uniformizationRate + storm::utility::one<ValueType>());
+                            element.setValue(-(exitRates[state] - element.getValue()) / uniformizationRate + storm::utility::one<ValueType>());
                         }
                     } else {
                         element.setValue(element.getValue() / uniformizationRate);
@@ -197,6 +235,8 @@ namespace storm {
                 }
                 ++currentRow;
             }
+            
+            return uniformizedMatrix;
         }
         
         template<class ValueType>
@@ -209,30 +249,63 @@ namespace storm {
                 element /= std::get<2>(foxGlynnResult);
             }
             
+            std::cout << "fox glynn:" << std::endl;
+            std::cout << "left: " << std::get<0>(foxGlynnResult) << std::endl;
+            std::cout << "right: " << std::get<1>(foxGlynnResult) << std::endl;
+            std::cout << "total: " << std::get<2>(foxGlynnResult) << std::endl;
+            int i = 0;
+            for (auto const& weight : std::get<3>(foxGlynnResult)) {
+                std::cout << "weight[" << i << "]: " << weight << std::endl;
+                ++i;
+            }
+            
             // Initialize result.
             std::vector<ValueType> result;
-            if (std::get<0>(foxGlynnResult) == 0) {
+            uint_fast64_t startingIteration = std::get<0>(foxGlynnResult);
+            if (startingIteration == 0) {
                 result = values;
                 storm::utility::vector::scaleVectorInPlace(result, std::get<3>(foxGlynnResult)[0]);
+                ++startingIteration;
             } else {
                 result = std::vector<ValueType>(values.size());
             }
             std::vector<ValueType> multiplicationResult(result.size());
             
-            // Perform the matrix-vector multiplications (without adding)
+            std::cout << "starting vector:" << std::endl;
+            for (int i = 0; i < result.size(); ++i) {
+                std::cout << i << ": " << result[i] << std::endl;
+            }
+            
+            // Perform the matrix-vector multiplications (without adding).
             if (std::get<0>(foxGlynnResult) > 1) {
-                linearEquationSolver.performMatrixVectorMultiplication(uniformizedMatrix, values, &result, std::get<0>(foxGlynnResult) - 1, &multiplicationResult);
+                linearEquationSolver.performMatrixVectorMultiplication(uniformizedMatrix, values, nullptr, std::get<0>(foxGlynnResult) - 1, &multiplicationResult);
+            }
+            
+            std::cout << "vector after initial mult phase:" << std::endl;
+            for (int i = 0; i < result.size(); ++i) {
+                std::cout << i << ": " << result[i] << std::endl;
             }
             
             // For the indices that fall in between the truncation points, we need to perform the matrix-vector
             // multiplication, scale and add the result.
             ValueType weight = 0;
-            std::function<ValueType(ValueType const&, ValueType const&)> addAndScale = [&weight] (ValueType const& a, ValueType const& b) { return a + weight * b; };
-            for (uint_fast64_t index = std::get<0>(foxGlynnResult); index <= std::get<1>(foxGlynnResult); ++index) {
+            std::function<ValueType(ValueType const&, ValueType const&)> addAndScale = [&weight] (ValueType const& a, ValueType const& b) { std::cout << "computing " << a << " + " << weight << " * " << b << " = " << (a + weight * b) << std::endl; return a + weight * b; };
+            for (uint_fast64_t index = startingIteration; index <= std::get<1>(foxGlynnResult); ++index) {
                 linearEquationSolver.performMatrixVectorMultiplication(uniformizedMatrix, values, nullptr, 1, &multiplicationResult);
                 
-                ValueType weight = std::get<3>(foxGlynnResult)[index - std::get<0>(foxGlynnResult)];
+                weight = std::get<3>(foxGlynnResult)[index - std::get<0>(foxGlynnResult)];
+                std::cout << "setting weight for index " << index << "(" << (index - std::get<0>(foxGlynnResult)) << ") " << " to " << std::get<3>(foxGlynnResult)[index - std::get<0>(foxGlynnResult)] << std::endl;
                 storm::utility::vector::applyPointwiseInPlace(result, values, addAndScale);
+
+                std::cout << "values after index: " << index << std::endl;
+                for (int i = 0; i < values.size(); ++i) {
+                    std::cout << i << ": " << values[i] << std::endl;
+                }
+
+                std::cout << "result after index: " << index << std::endl;
+                for (int i = 0; i < result.size(); ++i) {
+                    std::cout << i << ": " << result[i] << std::endl;
+                }
             }
             
             return result;

src/parser/FormulaParser.cpp

@@ -44,7 +44,7 @@ namespace storm {
             notStateFormula = (-unaryBooleanOperator_ >> atomicStateFormula)[qi::_val = phoenix::bind(&FormulaParser::createUnaryBooleanStateFormula, phoenix::ref(*this), qi::_2, qi::_1)];
             notStateFormula.name("negation formula");
             
-            eventuallyFormula = (qi::lit("F") >> -(qi::lit("<=") >> qi::uint_) >> pathFormulaWithoutUntil)[qi::_val = phoenix::bind(&FormulaParser::createEventuallyFormula, phoenix::ref(*this), qi::_1, qi::_2)];
+            eventuallyFormula = (qi::lit("F") >> -timeBound >> pathFormulaWithoutUntil)[qi::_val = phoenix::bind(&FormulaParser::createEventuallyFormula, phoenix::ref(*this), qi::_1, qi::_2)];
             eventuallyFormula.name("eventually formula");
             
             globallyFormula = (qi::lit("G") >> pathFormulaWithoutUntil)[qi::_val = phoenix::bind(&FormulaParser::createGloballyFormula, phoenix::ref(*this), qi::_1)];
@@ -56,12 +56,15 @@ namespace storm {
             pathFormulaWithoutUntil = eventuallyFormula | globallyFormula | nextFormula | stateFormula;
             pathFormulaWithoutUntil.name("path formula");
             
-            untilFormula = pathFormulaWithoutUntil[qi::_val = qi::_1] >> *(qi::lit("U") >> -(qi::lit("<=") >> qi::uint_) >> pathFormulaWithoutUntil)[qi::_val = phoenix::bind(&FormulaParser::createUntilFormula, phoenix::ref(*this), qi::_val, qi::_1, qi::_2)];
+            untilFormula = pathFormulaWithoutUntil[qi::_val = qi::_1] >> *(qi::lit("U") >> -timeBound >> pathFormulaWithoutUntil)[qi::_val = phoenix::bind(&FormulaParser::createUntilFormula, phoenix::ref(*this), qi::_val, qi::_1, qi::_2)];
             untilFormula.name("until formula");
             
             conditionalFormula = untilFormula[qi::_val = qi::_1] >> *(qi::lit("||") >> untilFormula)[qi::_val = phoenix::bind(&FormulaParser::createConditionalFormula, phoenix::ref(*this), qi::_val, qi::_1)];
             conditionalFormula.name("conditional formula");
             
+            timeBound = (qi::lit("[") > qi::double_ > qi::lit(",") > qi::double_ > qi::lit("]"))[qi::_val = phoenix::construct<std::pair<double, double>>(qi::_1, qi::_2)] | (qi::lit("<=") >> strict_double)[qi::_val = phoenix::construct<std::pair<double, double>>(0, qi::_1)] | (qi::lit("<=") >  qi::uint_)[qi::_val = qi::_1];
+            timeBound.name("time bound");
+            
             pathFormula = conditionalFormula;
             pathFormula.name("path formula");
             
@@ -193,9 +196,14 @@ namespace storm {
             return std::shared_ptr<storm::logic::Formula>(new storm::logic::AtomicLabelFormula(label));
         }
         
-        std::shared_ptr<storm::logic::Formula> FormulaParser::createEventuallyFormula(boost::optional<unsigned> const& stepBound, std::shared_ptr<storm::logic::Formula> const& subformula) const {
-            if (stepBound) {
-                return std::shared_ptr<storm::logic::Formula>(new storm::logic::BoundedUntilFormula(createBooleanLiteralFormula(true), subformula, static_cast<uint_fast64_t>(stepBound.get())));
+        std::shared_ptr<storm::logic::Formula> FormulaParser::createEventuallyFormula(boost::optional<boost::variant<std::pair<double, double>, uint_fast64_t>> const& timeBound, std::shared_ptr<storm::logic::Formula> const& subformula) const {
+            if (timeBound) {
+                if (timeBound.get().which() == 0) {
+                    std::pair<double, double> const& bounds = boost::get<std::pair<double, double>>(timeBound.get());
+                    return std::shared_ptr<storm::logic::Formula>(new storm::logic::BoundedUntilFormula(createBooleanLiteralFormula(true), subformula, bounds.first, bounds.second));
+                } else {
+                    return std::shared_ptr<storm::logic::Formula>(new storm::logic::BoundedUntilFormula(createBooleanLiteralFormula(true), subformula, static_cast<uint_fast64_t>(boost::get<uint_fast64_t>(timeBound.get()))));
+                }
             } else {
                 return std::shared_ptr<storm::logic::Formula>(new storm::logic::EventuallyFormula(subformula));
             }
@@ -209,9 +217,14 @@ namespace storm {
             return std::shared_ptr<storm::logic::Formula>(new storm::logic::NextFormula(subformula));
         }
         
-        std::shared_ptr<storm::logic::Formula> FormulaParser::createUntilFormula(std::shared_ptr<storm::logic::Formula> const& leftSubformula, boost::optional<unsigned> const& stepBound, std::shared_ptr<storm::logic::Formula> const& rightSubformula) {
-            if (stepBound) {
-                return std::shared_ptr<storm::logic::Formula>(new storm::logic::BoundedUntilFormula(leftSubformula, rightSubformula, static_cast<uint_fast64_t>(stepBound.get())));
+        std::shared_ptr<storm::logic::Formula> FormulaParser::createUntilFormula(std::shared_ptr<storm::logic::Formula> const& leftSubformula, boost::optional<boost::variant<std::pair<double, double>, uint_fast64_t>> const& timeBound, std::shared_ptr<storm::logic::Formula> const& rightSubformula) {
+            if (timeBound) {
+                if (timeBound.get().which() == 0) {
+                    std::pair<double, double> const& bounds = boost::get<std::pair<double, double>>(timeBound.get());
+                    return std::shared_ptr<storm::logic::Formula>(new storm::logic::BoundedUntilFormula(leftSubformula, rightSubformula, bounds.first, bounds.second));
+                } else {
+                    return std::shared_ptr<storm::logic::Formula>(new storm::logic::BoundedUntilFormula(leftSubformula, rightSubformula, static_cast<uint_fast64_t>(boost::get<uint_fast64_t>(timeBound.get()))));
+                }
             } else {
                 return std::shared_ptr<storm::logic::Formula>(new storm::logic::UntilFormula(leftSubformula, rightSubformula));
             }

src/parser/FormulaParser.h

@@ -145,14 +145,17 @@ namespace storm {
             qi::rule<Iterator, std::shared_ptr<storm::logic::Formula>(), Skipper> eventuallyFormula;
             qi::rule<Iterator, std::shared_ptr<storm::logic::Formula>(), Skipper> nextFormula;
             qi::rule<Iterator, std::shared_ptr<storm::logic::Formula>(), Skipper> globallyFormula;
-            qi::rule<Iterator, std::shared_ptr<storm::logic::Formula>(), Skipper> boundedUntilFormula;
             qi::rule<Iterator, std::shared_ptr<storm::logic::Formula>(), Skipper> untilFormula;
+            qi::rule<Iterator, boost::variant<std::pair<double, double>, uint_fast64_t>(), Skipper> timeBound;
             
             qi::rule<Iterator, std::shared_ptr<storm::logic::Formula>(), Skipper> rewardPathFormula;
             qi::rule<Iterator, std::shared_ptr<storm::logic::Formula>(), Skipper> cumulativeRewardFormula;
             qi::rule<Iterator, std::shared_ptr<storm::logic::Formula>(), Skipper> reachabilityRewardFormula;
             qi::rule<Iterator, std::shared_ptr<storm::logic::Formula>(), Skipper> instantaneousRewardFormula;
             
+            // Parser that is used to recognize doubles only (as opposed to Spirit's double_ parser).
+            boost::spirit::qi::real_parser<double, boost::spirit::qi::strict_real_policies<double>> strict_double;
+            
             // Methods that actually create the expression objects.
             std::shared_ptr<storm::logic::Formula> createInstantaneousRewardFormula(unsigned stepCount) const;
             std::shared_ptr<storm::logic::Formula> createCumulativeRewardFormula(unsigned stepBound) const;
@@ -160,10 +163,10 @@ namespace storm {
             std::shared_ptr<storm::logic::Formula> createAtomicExpressionFormula(storm::expressions::Expression const& expression) const;
             std::shared_ptr<storm::logic::Formula> createBooleanLiteralFormula(bool literal) const;
             std::shared_ptr<storm::logic::Formula> createAtomicLabelFormula(std::string const& label) const;
-            std::shared_ptr<storm::logic::Formula> createEventuallyFormula(boost::optional<unsigned> const& stepBound, std::shared_ptr<storm::logic::Formula> const& subformula) const;
+            std::shared_ptr<storm::logic::Formula> createEventuallyFormula(boost::optional<boost::variant<std::pair<double, double>, uint_fast64_t>> const& timeBound, std::shared_ptr<storm::logic::Formula> const& subformula) const;
             std::shared_ptr<storm::logic::Formula> createGloballyFormula(std::shared_ptr<storm::logic::Formula> const& subformula) const;
             std::shared_ptr<storm::logic::Formula> createNextFormula(std::shared_ptr<storm::logic::Formula> const& subformula) const;
-            std::shared_ptr<storm::logic::Formula> createUntilFormula(std::shared_ptr<storm::logic::Formula> const& leftSubformula, boost::optional<unsigned> const& stepBound, std::shared_ptr<storm::logic::Formula> const& rightSubformula);
+            std::shared_ptr<storm::logic::Formula> createUntilFormula(std::shared_ptr<storm::logic::Formula> const& leftSubformula, boost::optional<boost::variant<std::pair<double, double>, uint_fast64_t>> const& timeBound, std::shared_ptr<storm::logic::Formula> const& rightSubformula);
             std::shared_ptr<storm::logic::Formula> createConditionalFormula(std::shared_ptr<storm::logic::Formula> const& leftSubformula, std::shared_ptr<storm::logic::Formula> const& rightSubformula) const;
             std::shared_ptr<storm::logic::Formula> createLongRunAverageOperatorFormula(std::tuple<boost::optional<storm::logic::OptimalityType>, boost::optional<storm::logic::ComparisonType>, boost::optional<double>> const& operatorInformation, std::shared_ptr<storm::logic::Formula> const& subformula) const;
             std::shared_ptr<storm::logic::Formula> createRewardOperatorFormula(std::tuple<boost::optional<storm::logic::OptimalityType>, boost::optional<storm::logic::ComparisonType>, boost::optional<double>> const& operatorInformation, std::shared_ptr<storm::logic::Formula> const& subformula) const;

src/settings/modules/GeneralSettings.cpp

@@ -73,7 +73,7 @@ namespace storm {
                 this->addOption(storm::settings::OptionBuilder(moduleName, counterexampleOptionName, false, "Generates a counterexample for the given PRCTL formulas if not satisfied by the model")
                                 .addArgument(storm::settings::ArgumentBuilder::createStringArgument("filename", "The name of the file to which the counterexample is to be written.").setDefaultValueString("-").setIsOptional(true).build()).setShortName(counterexampleOptionShortName).build());
                 this->addOption(storm::settings::OptionBuilder(moduleName, bisimulationOptionName, false, "Sets whether to perform bisimulation minimization.").setShortName(bisimulationOptionShortName).build());
-                this->addOption(storm::settings::OptionBuilder(moduleName, transitionRewardsOptionName, "", "If given, the transition rewards are read from this file and added to the explicit model. Note that this requires the model to be given as an explicit model (i.e., via --" + explicitOptionName + ").")
+                this->addOption(storm::settings::OptionBuilder(moduleName, transitionRewardsOptionName, false, "If given, the transition rewards are read from this file and added to the explicit model. Note that this requires the model to be given as an explicit model (i.e., via --" + explicitOptionName + ").")
                                 .addArgument(storm::settings::ArgumentBuilder::createStringArgument("filename", "The file from which to read the transition rewards.").addValidationFunctionString(storm::settings::ArgumentValidators::existingReadableFileValidator()).build()).build());
                 this->addOption(storm::settings::OptionBuilder(moduleName, stateRewardsOptionName, false, "If given, the state rewards are read from this file and added to the explicit model. Note that this requires the model to be given as an explicit model (i.e., via --" + explicitOptionName + ").")
                                 .addArgument(storm::settings::ArgumentBuilder::createStringArgument("filename", "The file from which to read the state rewards.").addValidationFunctionString(storm::settings::ArgumentValidators::existingReadableFileValidator()).build()).build());

src/storage/SparseMatrix.cpp

@@ -239,7 +239,7 @@ namespace storm {
         template<typename ValueType>
         SparseMatrix<ValueType>::SparseMatrix(index_type columnCount, std::vector<index_type> const& rowIndications, std::vector<MatrixEntry<index_type, ValueType>> const& columnsAndValues, std::vector<index_type> const& rowGroupIndices) : rowCount(rowIndications.size() - 1), columnCount(columnCount), entryCount(columnsAndValues.size()), nonzeroEntryCount(0), columnsAndValues(columnsAndValues), rowIndications(rowIndications), rowGroupIndices(rowGroupIndices) {
             for (auto const& element : *this) {
-                if (!comparator.isZero(element.getValue())) {
+                if (element.getValue() != storm::utility::zero<ValueType>()) {
                     ++this->nonzeroEntryCount;
                 }
             }
@@ -248,7 +248,7 @@ namespace storm {
         template<typename ValueType>
         SparseMatrix<ValueType>::SparseMatrix(index_type columnCount, std::vector<index_type>&& rowIndications, std::vector<MatrixEntry<index_type, ValueType>>&& columnsAndValues, std::vector<index_type>&& rowGroupIndices) : rowCount(rowIndications.size() - 1), columnCount(columnCount), entryCount(columnsAndValues.size()), nonzeroEntryCount(0), columnsAndValues(std::move(columnsAndValues)), rowIndications(std::move(rowIndications)), rowGroupIndices(std::move(rowGroupIndices)) {
             for (auto const& element : *this) {
-                if (!comparator.isZero(element.getValue())) {
+                if (element.getValue() != storm::utility::zero<ValueType>()) {
                     ++this->nonzeroEntryCount;
                 }
             }
@@ -613,7 +613,7 @@ namespace storm {
             // First, we need to count how many entries each column has.
             for (index_type group = 0; group < columnCount; ++group) {
                 for (auto const& transition : joinGroups ? this->getRowGroup(group) : this->getRow(group)) {
-                    if (!comparator.isZero(transition.getValue())) {
+                    if (transition.getValue() != storm::utility::zero<ValueType>()) {
                         ++rowIndications[transition.getColumn() + 1];
                     }
                 }
@@ -632,7 +632,7 @@ namespace storm {
             // Now we are ready to actually fill in the values of the transposed matrix.
             for (index_type group = 0; group < columnCount; ++group) {
                 for (auto const& transition : joinGroups ? this->getRowGroup(group) : this->getRow(group)) {
-                    if (!comparator.isZero(transition.getValue())) {
+                    if (transition.getValue() != storm::utility::zero<ValueType>()) {
                         columnsAndValues[nextIndices[transition.getColumn()]] = std::make_pair(group, transition.getValue());
                         nextIndices[transition.getColumn()]++;
                     }

src/storage/SparseMatrix.h

@@ -839,9 +839,6 @@ namespace storm {
             
             // A vector indicating the row groups of the matrix.
             std::vector<index_type> rowGroupIndices;
-            
-            // A comparator that can be used to check whether some values satisfy some conditions.
-            storm::utility::ConstantsComparator<value_type> comparator;
         };
 
     } // namespace storage

src/utility/numerical.h

@@ -14,9 +14,11 @@ namespace storm {
                 storm::utility::ConstantsComparator<ValueType> comparator;
                 STORM_LOG_THROW(!comparator.isZero(lambda), storm::exceptions::InvalidArgumentException, "Error in Fox-Glynn algorithm: lambda must not be zero.");
                 
+                std::cout << "calling Fox-Glynn with " << lambda << ", " << overflow << ", " << underflow << ", " << accuracy << std::endl;
+                
                 // This code is a modified version of the one in PRISM. According to their implementation, for lambda
                 // smaller than 400, we compute the result using the naive method.
-                if (lambda < 25) {
+                if (lambda < 400) {
                     ValueType eToPowerMinusLambda = std::exp(-lambda);
                     ValueType targetValue = (1 - accuracy) / eToPowerMinusLambda;
                     std::vector<ValueType> weights;
@@ -55,16 +57,16 @@ namespace storm {
                         ValueType aLambda = 0, bLambda = 0, sqrtLambda = 0;
                         if (m < 400) {
                             sqrtLambda = std::sqrt(400.0);
-                            aLambda = (1.0 + 1.0 / 400.0) * exp(0.0625) * sqrt2;
-                            bLambda = (1.0 + 1.0 / 400.0) * exp(0.125 / 400.0);
+                            aLambda = (1.0 + 1.0 / 400.0) * std::exp(0.0625) * sqrt2;
+                            bLambda = (1.0 + 1.0 / 400.0) * std::exp(0.125 / 400.0);
                         } else {
                             sqrtLambda = std::sqrt(lambda);
-                            aLambda = (1.0 + 1.0 / lambda) * exp(0.0625) * sqrt2;
-                            bLambda = (1.0 + 1.0 / lambda) * exp(0.125 / lambda);
+                            aLambda = (1.0 + 1.0 / lambda) * std::exp(0.0625) * sqrt2;
+                            bLambda = (1.0 + 1.0 / lambda) * std::exp(0.125 / lambda);
                         }
                         
                         // Use Corollary 1 from the paper to find the right truncation point.
-                        uint_fast64_t k = 3;
+                        uint_fast64_t k = 4;
                         
                         ValueType dkl = 1.0 / (1 - std::exp(-(2.0 / 9.0) * (k * sqrt2 * sqrtLambda + 1.5)));
                         
@@ -77,22 +79,27 @@ namespace storm {
                             dkl = 1.0 / (1 - std::exp(-(2.0 / 9.0)*(k * sqrt2 * sqrtLambda + 1.5)));
                         }
                         
+                        std::cout << "k for upper: " << k << std::endl;
+                        
                         // Left hand side of the equation in Corollary 1.
                         rightTruncationPoint = static_cast<uint_fast64_t>(std::ceil((m + k * sqrt2 * sqrtLambda + 1.5)));
                         
                         // Use Corollary 2 to find left truncation point.
-                        k = 3;
+                        k = 4;
                         
                         // Right hand side of the equation in Corollary 2.
-                        while ((accuracy / 2.0) < ((bLambda * exp(-(k*k / 2.0))) / (k * sqrt2 * sqrtpi))) {
+                        while ((accuracy / 2.0) < ((bLambda * std::exp(-(k*k / 2.0))) / (k * sqrt2 * sqrtpi))) {
+                            std::cout << "k=" << k << " produces: " << ((bLambda * std::exp(-(k*k / 2.0))) / (k * sqrt2 * sqrtpi)) << std::endl;
                             ++k;
                         }
+                        std::cout << "k=" << k << " produces: " << ((bLambda * std::exp(-(k*k / 2.0))) / (k * sqrt2 * sqrtpi)) << std::endl;
+                        std::cout << "k for lower: " << k << std::endl;
                         
                         // Left hand side of the equation in Corollary 2.
-                        leftTruncationPoint = std::max(static_cast<uint_fast64_t>(std::trunc(m - k * sqrtLambda - 1.5)), uint_fast64_t(0));
+                        leftTruncationPoint = static_cast<uint_fast64_t>(std::max(std::trunc(m - k * sqrtLambda - 1.5), storm::utility::zero<ValueType>()));
                         
                         // Check for underflow.
-                        STORM_LOG_THROW(static_cast<uint_fast64_t>(std::trunc((m - k * sqrtLambda - 1.5))) >= 0, storm::exceptions::OutOfRangeException, "Error in Fox-Glynn algorithm: Underflow of left truncation point.");
+                        STORM_LOG_THROW(std::trunc((m - k * sqrtLambda - 1.5)) >= 0, storm::exceptions::OutOfRangeException, "Error in Fox-Glynn algorithm: Underflow of left truncation point.");
                     }
                     
                     std::vector<ValueType> weights(rightTruncationPoint - leftTruncationPoint + 1);