Add spec test for pumping example

src/main/scala/symsim/examples/concrete/pumping/Pump.scala

@@ -1,26 +1,21 @@
 package symsim
 package examples.concrete.pumping
-/** A case study of control of a water extraction at Marbjerg (operated by HOFOR).
-  *
-  * Source: Andreas Holck Høeg-Petersen. Reinforcement Learning for Controlling
-  * Groundwater Extraction. MSc thesis at the IT University of Copenhagen.
-  * August 2021. (esp. Chapters 2 and 3)
-  */
 
 import symsim.concrete.Randomized
+import scala.math
 
 
 /** Constants as fit in the model of AHHP (reference above).
   *
   * They are called a, b, c, and d parameters in his report.
   */ 
 val c1: Double = 1.0 / 12
-val c2: Double = 0.15 / (c1*80*12*24*365)
-val c3: Double = 0.1 / (12*24*365)
+val c2: Double = 0.15 / (c1 * 80 * 12 * 24 * 365)
+val c3: Double = 0.1 / (12 * 24 * 365)
 val c4: Double = 0.05
 
 /** Size of the history window for the head-level */
-val k: Int = 5
+val MemoryWindowSize: Int = 5
 
 
 /** Complete state of a pump, both observable and unobservable. */
@@ -33,6 +28,14 @@ case class PumpState (
   w:   Double,        // Water level in the ground deposits (unobservable)
   phm: List[Double]): // Past head means (a sliding window of history) 
 
+  require (f >= FlowMin,  s"Flow is too low: $f")
+  require (f <= FlowMax,  s"Flow is too high: $f")
+  require (h <= HeadMax,  s"Head is too high: $h")
+  require (hm >= HeadMin, s"Head mean is too low: $hm")
+  require (hm <= HeadMax, s"Head mean is too high: $hm")
+  require (w >= WaterMin, s"Water level is too low: $w")
+  require (w <= WaterMax, s"Water level is too high: $w")
+
   override def toString: String =
     s"[flow=$f, head=$h, head mean=$hm, tank level=$tl, time=$t, "
       + s"water=$w, past head means=$phm]"
@@ -42,9 +45,9 @@ end PumpState
 
 /** Discretize cut points for flow, head, and tank variables./
   *
-  * We observe to the first point on the list that is lower than the actual
-  * value of the variable.  So we 'round down'. Values below the last element
-  * in the cut-point list are rounded up (to the last value).
+  * We discretize to the first point on the list that is lower than the actual
+  * value of the variable.  So we 'round down'. Values below the last element in
+  * the cut-point list are rounded up (to the last value).
   *
   * For now, we have preconditions about low values in discretization so the
   * getOrElse (rounding up) should never happen.
@@ -57,24 +60,24 @@ def closest (value: Double) (cutPoints: List[Double]): Double =
     .find { value >= _ }
     .getOrElse (cutPoints.last)
 
-val HEAD_MIN: Double = 7.0
-val HEAD_MAX: Double = 10.84
-val FLOW_MIN: Double = 0.0
-val FLOW_MAX: Double = 120.0
-val TANK_MIN: Double = 0.0
-val TANK_MAX: Double = 2000.0
-val WATER_MIN: Double = 0.0
-val WATER_MAX: Double = 9.11
+val HeadMin  = 7.0
+val HeadMax  = 10.84
+val FlowMin  = 0.0
+val FlowMax  = 120.0
+val TankMin  = 0.0
+val TankMax  = 2000.0
+val WaterMin = 0.0
+val WaterMax = 9.11
 
-val flowCutPoints = List (FLOW_MAX, 115.0, 110.0, 105.0, 100.0, 95.0,
-  90.0, 85.0, 80.0, 75.0, 70.0, 65.0, 60.0, 55.0, 50.0, FLOW_MIN)
+val FlowCutPoints = List (FlowMax, 115.0, 110.0, 105.0, 100.0, 95.0,
+  90.0, 85.0, 80.0, 75.0, 70.0, 65.0, 60.0, 55.0, 50.0, FlowMin)
 
-val headCutPoints = List (HEAD_MAX, 10.68, 10.52, 10.36, 10.2, 10.04,
+val HeadCutPoints = List (HeadMax, 10.68, 10.52, 10.36, 10.2, 10.04,
   9.88, 9.72, 9.56, 9.4, 9.24, 9.08, 8.92, 8.76, 8.6, 8.44, 8.28, 
-  8.12, 7.96, 7.8, 7.64, 7.48, 7.32, 7.16, HEAD_MIN)
+  8.12, 7.96, 7.8, 7.64, 7.48, 7.32, 7.16, HeadMin)
 
-val tankCutPoints = List (TANK_MAX, 1800.0, 1600.0, 1400.0, 
-  1200.0, 1000.0, 800.0, 600.0, 400.0, 200.0, TANK_MIN)
+val TankCutPoints = List (TankMax, 1800.0, 1600.0, 1400.0, 
+  1200.0, 1000.0, 800.0, 600.0, 400.0, 200.0, TankMin)
 
 /** The (discrete) state, observable by the pump controller. */
 case class ObservablePumpState (
@@ -83,6 +86,11 @@ case class ObservablePumpState (
   hm: Double,  // Discretized, derived head mean
   tl: Double): // Discretized observable tank level
 
+  require (FlowCutPoints.contains (f), s"Wrong discrete flow: $f")
+  require (HeadCutPoints.contains (h), s"Wrong discrete head level: $h")
+  require (HeadCutPoints.contains (hm), s"Wrong discrete mean head level:  $hm")
+  require (TankCutPoints.contains (tl), s"Wrong discrete tank level: $tl")
+
   override def toString: String =
     s"(flow=$f, head=$h, head mean=$hm, tank level=$tl)"
 
@@ -93,70 +101,66 @@ type PumpAction = Double
 type PumpReward = Double
 
 
+/** A case study of control of a water extraction at Marbjerg (operated by HOFOR).
+  *
+  * Source: Andreas Holck Høeg-Petersen. Reinforcement Learning for Controlling
+  * Groundwater Extraction. MSc thesis at the IT University of Copenhagen.
+  * August 2021. (esp. Chapters 2 and 3)
+  */
 object Pump extends 
-  Agent[PumpState, ObservablePumpState, PumpAction, PumpReward, Randomized],
-  Episodic:
-
-  /** An upper bound on episode duration. Used only in testing */
-  val TimeHorizon: Int = 20000
-
+  Agent[PumpState, ObservablePumpState, PumpAction, PumpReward, Randomized]:
 
   /** If isFinal is true in a state, then the episode is complete */
   def isFinal (s: PumpState): Boolean =
-    s.t >= 4000 || s.h < HEAD_MIN || s.tl > TANK_MAX || s.tl < TANK_MIN
+    s.t >= 4000 || s.h < HeadMin || s.tl > TankMax || s.tl < TankMin
 
 
   def observe (s: PumpState): ObservablePumpState =
-    require (s.tl >= TANK_MIN, s"s.tl = ${s.tl} >= $TANK_MIN")
-
-    val df = closest (s.f) (flowCutPoints)
-    val dh = closest (s.h) (headCutPoints)
-    val dhm = closest (s.hm) (headCutPoints)
-    val dtl = closest (s.tl) (tankCutPoints)
-
+    require (s.tl >= TankMin, s"s.tl = ${s.tl} >= $TankMin")
+    val df = closest (s.f) (FlowCutPoints)
+    val dh = closest (s.h) (HeadCutPoints)
+    val dhm = closest (s.hm) (HeadCutPoints)
+    val dtl = closest (s.tl) (TankCutPoints)
     ObservablePumpState (df, dh, dhm, dtl)
 
 
   def reward (source: PumpState) (target: PumpState) (a: PumpAction): Double =
-    if target.h < HEAD_MIN || target.tl < TANK_MIN || target.tl > TANK_MAX 
-    then -9999
+    if isFinal (target) then -9999.0
     else 
-      val flowReward = if target.f != source.f then -0.5 else 0
+      val flowReward = if target.f != source.f then -0.5 else 0.0
       val headCost = (1 + (target.h - target.hm).abs) * (1 + (target.h - target.hm).abs)
       flowReward - headCost
 
 
-  val amp: Double = 0.41852857594808646
-  val freq: Double = 0.0000597030105413
-  val phase: Double = -6347.109214682171
+  val Amp   = 0.41852857594808646
+  val Freq  = 0.0000597030105413
+  val Phase = -6347.109214682171
 
 
   override def step (s: PumpState) (a: PumpAction)
     : Randomized[(PumpState, PumpReward)] =
     require (instances.enumAction.membersAscending.contains (a))
     for
-      nf   <- Randomized.gaussian (0.0, 1.0)
-      f1   =  a + nf
-      nd   <- Randomized.gaussian (0.1, 0.01)
-      cd   <- getDemand (s.t%24 + 1)
-      d    =  cd + nd
+      f1  <- Randomized.gaussian (a, 1.0)
+      cd  <- demand (s.t % 24 + 1)
+      d   <- Randomized.gaussian (cd + 0.1, 0.01)
       tl1  = s.tl + c1 * (f1 - d)
-      h1   = s.h + c4 * (s.w + (c1*f1 / Math.PI))
-      nw   <- Randomized.gaussian (0.0, 1.0)
-      w1   = s.w - c2*(c1*f1) + c3 + 
-             (amp*Math.sin (2*Math.PI*(s.t + phase) / freq)) + nw
-      hm1  = (1.0 / k)*s.phm.sum
-      phm1 = (s.hm:: s.phm).slice (0, k)
+      h1   = s.h + c4 * (s.w + (c1 * f1 / math.Pi))
+      nw  <- Randomized.gaussian (0.0, 1.0)
+      w1   = s.w - c2 * c1 * f1 + c3 + 
+             (Amp * math.sin (2 * math.Pi * (s.t + Phase) / Freq)) + nw
+      hm1  = (1.0/MemoryWindowSize) * s.phm.sum
+      phm1 = (s.hm:: s.phm).slice (0, MemoryWindowSize)
       s1   = PumpState (f1, h1, hm1, tl1, s.t + 1, w1, phm1)
       pr   = reward (s) (s1) (a)
     yield (s1, pr)
 
 
-  def getDemand (t: Int): Randomized[Double] =
+  def demand (t: Int): Randomized[Double] =
     require (t >= 0 && t <= 24)
     if t < 5 then Randomized.between (5.0, 15.0)
-    if t < 12 then Randomized.between (15.0, 45.0)
-    if t < 22 then Randomized.between (20.0, 38.0)
+    else if t < 12 then Randomized.between (15.0, 45.0)
+    else if t < 22 then Randomized.between (20.0, 38.0)
     else Randomized.between (5.0, 20.0)
 
 
@@ -196,10 +200,10 @@ object PumpInstances
 
   given enumState: BoundedEnumerable[ObservablePumpState] =
     val ss: List[ObservablePumpState] = for
-      f  <- flowCutPoints
-      h  <- headCutPoints 
-      hm <- headCutPoints 
-      tl <- tankCutPoints
+      f  <- FlowCutPoints
+      h  <- HeadCutPoints 
+      hm <- HeadCutPoints 
+      tl <- TankCutPoints
     yield ObservablePumpState (f, h, hm, tl)
     BoundedEnumerableFromList (ss*)
 
@@ -213,13 +217,13 @@ object PumpInstances
     concrete.Randomized.canTestInRandomized
 
   lazy val genPumpState: Gen[PumpState] = for
-    f   <- Gen.choose[Double] (FLOW_MIN, FLOW_MAX)
-    h   <- Gen.choose[Double] (HEAD_MIN, HEAD_MAX)
-    hm  <- Gen.choose[Double] (HEAD_MIN, HEAD_MAX)
-    tl  <- Gen.choose[Double] (TANK_MIN, TANK_MAX)
+    f   <- Gen.choose[Double] (FlowMin, FlowMax)
+    h   <- Gen.choose[Double] (HeadMin, HeadMax)
+    hm  <- Gen.choose[Double] (HeadMin, HeadMax)
+    tl  <- Gen.choose[Double] (TankMin, TankMax)
     t   <- Gen.choose[Int] (0, 24)
-    w   <- Gen.choose[Double] (WATER_MIN, WATER_MAX)
-    phm <- Gen.listOfN (5, Gen.choose (HEAD_MIN, HEAD_MAX))
+    w   <- Gen.choose[Double] (WaterMin, WaterMax)
+    phm <- Gen.listOfN (5, Gen.choose (HeadMin, HeadMax))
   yield PumpState (f, h, hm, tl, t, w, phm)
 
   given arbitraryState: Arbitrary[PumpState] = Arbitrary (genPumpState)

src/main/scala/symsim/laws/AgentLaws.scala

@@ -22,7 +22,8 @@ case class AgentLaws[State, ObservableState, Action, Reward, Scheduler[_]]
 
   import agent.instances.given
 
-  val observableStates = agent.instances.enumState.membersAscending
+  val observableStates = 
+    agent.instances.enumState.membersAscending.toSet
 
   val laws: RuleSet = SimpleRuleSet (
     "agent",

src/test/scala/symsim/examples/concrete/pumping/Experiments.scala

@@ -2,10 +2,9 @@ package symsim
 package examples.concrete.pumping
 
 class Experiments
-  extends ExperimentSpec[PumpState,ObservablePumpState,PumpAction]:
+  extends ExperimentSpec[PumpState, ObservablePumpState, PumpAction]:
 
   // Import evidence that states and actions can be enumerated
-  import Pump.*
   import Pump.instances.given
 
   val sarsa = symsim.concrete.ConcreteSarsa (

src/test/scala/symsim/examples/concrete/pumping/PumpIsAgentSpec.scala

@@ -1,11 +1,10 @@
 package symsim
 package examples.concrete.pumping
 
-import laws.AgentLaws
-import laws.EpisodicLaws
+import symsim.laws.AgentLaws
+import symsim.laws.EpisodicLaws
 
 class PumpIsAgentSpec
   extends SymSimSpec:
 
   checkAll ("concrete.pumping.Pump is an Agent", AgentLaws (Pump).laws)
-  checkAll ("concrete.pumping.Pump is Episodic", EpisodicLaws (Pump).laws)

src/test/scala/symsim/examples/concrete/pumping/PumpSpec.scala

@@ -0,0 +1,71 @@
+package symsim
+package examples.concrete.pumping
+
+import symsim.concrete.Randomized.given
+import symsim.CanTestIn.given
+import symsim.CanTestIn.*
+
+import org.scalacheck.{Arbitrary, Gen, Prop}
+import org.scalacheck.Arbitrary.*
+import org.scalacheck.Prop.forAll
+import org.scalacheck.Prop.*
+import examples.concrete.pumping.*
+
+// To eliminate the warning on Pumping, until scalacheck makes it open
+import scala.language.adhocExtensions
+
+
+class PumpSpec
+  extends org.scalacheck.Properties ("Pumping"):
+
+    // Generators of test data
+    val flow = Gen.choose[Double] (FlowMin, FlowMax)
+    val head = Gen.choose[Double] (HeadMin, HeadMax)
+    val head_mean = Gen.choose[Double] (HeadMin, HeadMax)
+    val tank = Gen.choose[Double] (TankMin, TankMax)
+    val time = Gen.choose[Int] (0, 24)
+    val water = Gen.choose[Double] (WaterMin, WaterMax)
+    val past_head_mean = Gen.listOfN (5, Gen.choose (HeadMin, HeadMax))
+    val actions = Gen.oneOf (Pump.instances.enumAction.membersAscending)
+    val obsStates = Gen.oneOf (Pump.instances.allObservableStates)
+
+
+    // Tests
+
+    property ("The tank level never can be negative") = {
+      forAll (flow, head, head_mean, tank, time, water, past_head_mean, actions) {
+        (f, h, hm, tl, t, w, phm, a) =>
+        for (s1, r) <- Pump.step (PumpState (f, h, hm, tl, t, w, phm)) (a)
+        yield s1.tl >= TankMin
+      }
+    }
+
+    property ("There is an observable state for each valid state in the environment") = {
+      forAll (flow, head, head_mean, tank, time, water, past_head_mean, actions) {
+        (f, h, hm, tl, t, w, phm, a) =>
+          val (s1, r) = Pump.step (PumpState (f, h, hm, tl, t, w, phm)) (a).head
+          exists (obsStates) {s => s == Pump.observe (s1)}
+      }
+    }
+
+    property ("The time can be more than TankMax (comparing is not correct)") = {
+      exists (time) { t =>
+        t > TankMax
+      }
+    }
+
+    property ("Test getDemand") = {
+      forAll (time) { t1 =>
+        exists (time) {t2 => Pump.demand (t1) != Pump.demand (t2)}
+      }
+    }
+
+    property ("There is an action for each state which results no overflow") = {
+      forAll (flow, head, head_mean, tank, time, water, past_head_mean) {
+        (f, h, hm, tl, t, w, phm) =>
+          exists(actions) { a =>
+            for (s1, r) <- Pump.step(PumpState(f, h, hm, tl, t, w, phm))(a)
+              yield (s1.tl <= TankMax)
+          }
+      }
+    }