Prevent division by zero in Z3

6 files changed, 61 insertions, 38 deletions

diff --git a/Assignment2/src/DTMC.dcl b/Assignment2/src/DTMC.dcl
index 04243c4..05145f0 100644
--- a/Assignment2/src/DTMC.dcl
+++ b/Assignment2/src/DTMC.dcl
@@ -9,6 +9,7 @@ from StdOverloaded import class toString
 :: *DTMC =
 	{ nr_states  :: !Int
 	, parameters :: ![String]
+	, unzero     :: ![String]
 	, states     :: !*{Maybe State}
 	}
 
diff --git a/Assignment2/src/DTMC.icl b/Assignment2/src/DTMC.icl
index bd6f9e1..87a9130 100644
--- a/Assignment2/src/DTMC.icl
+++ b/Assignment2/src/DTMC.icl
@@ -60,8 +60,13 @@ stateElimination :: !*DTMC -> *DTMC
 stateElimination dtmc=:{nr_states}
 # (ps,dtmc) = paths dtmc
 = case [p \\ p=:[_:_:_:_] <- ps] of
-	[]    -> {dtmc & states=prune (nr_states-1) dtmc.states}
-	[p:_] -> stateElimination {dtmc & states=uncurry elim (penultimate_trans p) dtmc.states}
+	[]
+		= {dtmc & states=prune (nr_states-1) dtmc.states}
+	[p:_]
+		# (unzs,states) = uncurry elim (penultimate_trans p) dtmc.states
+		# dtmc & states = states
+		# dtmc & unzero = unzs ++ dtmc.unzero
+		= stateElimination dtmc
 		with
 			penultimate_trans [x,y,_] = (x,y)
 			penultimate_trans [_:xs] = penultimate_trans xs
@@ -76,55 +81,52 @@ where
 		remove :: !Int !*{Maybe State} -> *{Maybe State}
 		remove i ss = {ss & [i]=Nothing}
 
-	elim :: !Int !Int !*{Maybe State} -> *{Maybe State}
+	elim :: !Int !Int !*{Maybe State} -> *(![String], !*{Maybe State})
 	elim s1 s ss
 	# (Just s1s,ss) = ss![s1]
 	# prob = 'M'.get s s1s.transitions
 
 	# (succ,ss) = ss![s]
-	| isNothing succ = ss
+	| isNothing succ = ([], ss)
 	# succ = fromJust succ
-	| succ.init || isAbsorbing succ = ss
+	| succ.init || isAbsorbing succ = ([], ss)
 	# succ = succ.transitions
 
 	# self_prob = fromMaybe "0.0" $ 'M'.get s succ
 	# succ = 'M'.toList succ
-	| s1 == s
-		# ss = seqSt (uncurry $ update_eq self_prob) succ ss
-		# (Just s1s,ss) = ss![s1]
-		# ss & [s1] = Just {s1s & transitions='M'.del s s1s.transitions}
-		= ss
-	| otherwise
-		# ss = seqSt (uncurry $ update_ne self_prob) succ ss
-		# (Just s1s,ss) = ss![s1]
-		# ss & [s1] = Just {s1s & transitions='M'.del s s1s.transitions}
-		= ss
+
+	# (unzs,ss) = mapSt (uncurry $ if (s1==s) update_eq update_ne self_prob) succ ss
+	# (Just s1s,ss) = ss![s1]
+	# ss & [s1] = Just {s1s & transitions='M'.del s s1s.transitions}
+	= (catMaybes unzs,ss)
 	where
-		update_eq :: !String !Int !String !*{Maybe State} -> *{Maybe State}
+		update_eq :: !String !Int !String !*{Maybe State} -> *(!Maybe String, !*{Maybe State})
 		update_eq self_prob s2 prob2 ss
-		| s == s2 = ss
+		| s == s2 = (Nothing, ss)
 		# (Just s1s,ss) = ss![s1]
 		# ss & [s1] = Just
 			{ s1s
 			& transitions = 'M'.put s2
-				("(" +++ prob2 +++ ")/(1.0-(" +++ self_prob +++ "))")
+				("(" +++ prob2 +++ ")/" +++ denominator)
 				s1s.transitions
 			}
-		= ss
+		= (Just $ "(not (= 0 " +++ parenthesise (parseInfix denominator) +++ "))", ss)
+		where denominator = "(1.0-(" +++ self_prob +++ "))"
 
-		update_ne :: !String !Int !String !*{Maybe State} -> *{Maybe State}
+		update_ne :: !String !Int !String !*{Maybe State} -> *(!Maybe String, !*{Maybe State})
 		update_ne self_prob s2 prob2 ss
-		| s == s2 = ss
+		| s == s2 = (Nothing, ss)
 		# (Just s1s,ss) = ss![s1]
 		# ps1s2 = fromMaybe "0.0" $ 'M'.get s2 s1s.transitions
 		# ps1s = fromMaybe "0.0" $ 'M'.get s s1s.transitions
 		# ss & [s1] = Just
 			{ s1s
 			& transitions = 'M'.put s2
-				(ps1s2 +++ "+(" +++ ps1s +++ ")*(" +++ prob2 +++ ")/(1.0-(" +++ self_prob +++ "))")
+				(ps1s2 +++ "+(" +++ ps1s +++ ")*(" +++ prob2 +++ ")/" +++ denominator)
 				s1s.transitions
 			}
-		= ss
+		= (Just $ "(not (= 0 " +++ parenthesise (parseInfix denominator) +++ "))", ss)
+		where denominator = "(1.0-(" +++ self_prob +++ "))"
 
 isAbsorbing :: !State -> Bool
 isAbsorbing s = all ((==) s.state_id) [to \\ (to,_) <- 'M'.toList s.transitions]
@@ -141,6 +143,7 @@ parseDTMCFromLines lines
 =
 	{ nr_states  = toInt nr_states
 	, parameters = []
+	, unzero     = []
 	, states     = {Just s \\ s <- states}
 	}
 where
@@ -294,7 +297,7 @@ assertProperty (op, prob, target) dtmc z3 w
 # formula = 'M'.get target init.transitions
 | isNothing formula = abort "Cannot reach property target from initial state\n"
 # formula = fromJust formula
-# formula = "( " <+ op <+ " " <+ prob <+ " " <+ parenthesise (parseInfix formula) <+ ")"
+# formula = "(" <+ op <+ " " <+ prob <+ " " <+ parenthesise (parseInfix formula) <+ ")"
 # w = addAssert z3 formula w
 = (dtmc, w)
 
@@ -322,20 +325,21 @@ Start w
 # props = map (parseProperty o trim) $ fromOk props
 
 # dtmc = addVars dtmc
-# dtmc = stateElimination dtmc
 # (config,dtmc) = repairConfig SumSquaredErrors dtmc
+# dtmc = stateElimination dtmc
 
 # (z3,w) = startZ3 w
 # (dtmc,w) = z3AddVars dtmc z3 w
 # (dtmc,w) = assertProperties props dtmc z3 w
 # w = seqSt (addAssert z3) config.totality_restrictions w
 # w = seqSt (addAssert z3) config.transition_restrictions w
+# w = seqSt (addAssert z3) dtmc.unzero w
 # w = addVariable z3 "Real" "cost" w
 # w = addAssert z3 ("(= cost " +++ config.minimalisation_goal +++ ")") w
-//# w = addMinimize z3 "cost" w // TODO: not supported by Z3!
+//# w = addMinimize z3 "cost" w // NOTE: z3 fails with this, but you can approach a bound by asserting an upper bound on cost
 
 # (sat,w) = checkSat z3 w
-| sat
+| not (sat=:Unsat)
 	# (model,w) = getModel z3 w
 	= model
 = "Could not satisfy the property\n"
diff --git a/Assignment2/src/Z3.dcl b/Assignment2/src/Z3.dcl
index 0f636f2..5cff893 100644
--- a/Assignment2/src/Z3.dcl
+++ b/Assignment2/src/Z3.dcl
@@ -4,11 +4,16 @@ from System.Process import :: ProcessHandle, :: ProcessIO
 
 :: Z3
 
+:: Satisfaction
+	= Sat
+	| Unsat
+	| Unknown
+
 startZ3 :: !*World -> (!Z3, !*World)
 addAssert :: !Z3 !String !*World -> *World
 addMinimize :: !Z3 !String !*World -> *World
 addVariable :: !Z3 !String !String !*World -> *World
-checkSat :: !Z3 !*World -> (!Bool, !*World)
+checkSat :: !Z3 !*World -> (!Satisfaction, !*World)
 
 // Note: getModel terminates z3
 getModel :: !Z3 !*World -> (!String, !*World)
diff --git a/Assignment2/src/Z3.icl b/Assignment2/src/Z3.icl
index ca869a1..4db349e 100644
--- a/Assignment2/src/Z3.icl
+++ b/Assignment2/src/Z3.icl
@@ -1,5 +1,7 @@
 implementation module Z3
 
+import StdList
+import StdMisc
 import StdString
 
 import System.Process
@@ -43,13 +45,17 @@ addVariable z3 type s w
 # (_, w) = writePipe (trace_id $ "(declare-const " <+ s <+ " " <+ type <+ ")\n") z3.Z3.stdIn w
 = w
 
-checkSat :: !Z3 !*World -> (!Bool, !*World)
+checkSat :: !Z3 !*World -> (!Satisfaction, !*World)
 checkSat z3 w
 # (_, w) = writePipe "(check-sat)\n" z3.Z3.stdIn w
 # (out, w) = readPipeBlocking z3.Z3.stdOut w
 # out = fromOk out
 # out = trace_id out
-= (endsWith "sat\n" out, w)
+= (case last $ init $ split "\n" out of
+	"sat"     -> Sat
+	"unsat"   -> Unsat
+	"unknown" -> Unknown
+	_         -> abort "Z3 failed\n", w)
 
 getModel :: !Z3 !*World -> (!String, !*World)
 getModel z3 w
diff --git a/Assignment2/src/die.prism b/Assignment2/src/die.prism
index 45442d6..06de40a 100644
--- a/Assignment2/src/die.prism
+++ b/Assignment2/src/die.prism
@@ -1,16 +1,18 @@
 dtmc
 
+const double p = 0.4;
+
 module die
 	s : [0..7] init 0;
 	d : [0..6] init 0;
 
-	[] s=0 -> 0.5 : (s'=1) + 0.5 : (s'=2);
-	[] s=1 -> 0.5 : (s'=3) + 0.5 : (s'=4);
-	[] s=2 -> 0.5 : (s'=5) + 0.5 : (s'=6);
-	[] s=3 -> 0.5 : (s'=1) + 0.5 : (s'=7) & (d'=1);
-	[] s=4 -> 0.5 : (s'=7) & (d'=2) + 0.5 : (s'=7) & (d'=3);
-	[] s=5 -> 0.5 : (s'=7) & (d'=4) + 0.5 : (s'=7) & (d'=5);
-	[] s=6 -> 0.5 : (s'=2) + 0.5 : (s'=7) & (d'=6);
+	[] s=0 -> p : (s'=1) + (1-p) : (s'=2);
+	[] s=1 -> p : (s'=3) + (1-p) : (s'=4);
+	[] s=2 -> (1-p) : (s'=5) + p : (s'=6);
+	[] s=3 -> p : (s'=1) + (1-p) : (s'=7) & (d'=1);
+	[] s=4 -> p : (s'=7) & (d'=2) + (1-p) : (s'=7) & (d'=3);
+	[] s=5 -> p : (s'=7) & (d'=4) + (1-p) : (s'=7) & (d'=5);
+	[] s=6 -> p : (s'=2) + (1-p) : (s'=7) & (d'=6);
 	[] s=7 -> 1 : (s'=7);
 endmodule
 
diff --git a/Assignment2/src/die.properties b/Assignment2/src/die.properties
index f507f66..ce7e591 100644
--- a/Assignment2/src/die.properties
+++ b/Assignment2/src/die.properties
@@ -1 +1,6 @@
-(EQ, 0.5, 7)
+(LT, 0.166666666, 7)
+(LT, 0.166666666, 8)
+(LT, 0.166666666, 9)
+(LT, 0.166666666, 10)
+(LT, 0.166666666, 11)
+(LT, 0.166666666, 12)