Started on variable step size

src/CBD/CBD.py

@@ -1,5 +1,6 @@
 from .util import enum
 from collections import namedtuple
+from copy import deepcopy
 
 InputLink = namedtuple("InputLink", ["block", "output_port"])
 Signal = namedtuple("Signal", ["time", "value"])
@@ -119,7 +120,7 @@ class BaseBlock:
         """
         name_output = "OUT1" if name_output is None else name_output
         assert name_output in self.__signals.keys()
-        return self.__signals[name_output] if name_output is not None else self.__signals["OUT1"]
+        return self.__signals[name_output]
 
     def getSignals(self):
         return self.__signals
@@ -225,16 +226,28 @@ class BaseBlock:
                 repr += "  On input " + key + ": " + in_block.getPath() + ":" + in_block.getBlockType() + "\n"
         return repr
 
+    def _rewind(self):
+        """
+        Rewinds the CBD model to the previous iteration.
+
+        Danger:
+            Normally, this function should only be used by the internal mechanisms
+            of the CBD simulator, **not** by a user. Using this function without a
+            full understanding of the simulator may result in undefined behaviour.
+        """
+        for signal in self.getSignals().keys():
+            self.__signals[signal].pop()
+
 class InputPortBlock(BaseBlock):
     """
     The input port of a CBD
     """
     def __init__(self, block_name, parent):
         BaseBlock.__init__(self, block_name, [], ["OUT1"])
-        self.parent = parent
+        self._parent = parent
 
     def compute(self, curIteration):
-        self.appendToSignal(self.parent.getInputSignal(curIteration, self.getBlockName()).value)
+        self.appendToSignal(self._parent.getInputSignal(curIteration, self.getBlockName()).value)
 
 
 class OutputPortBlock(BaseBlock):
@@ -243,7 +256,7 @@ class OutputPortBlock(BaseBlock):
     """
     def __init__(self, block_name, parent):
         BaseBlock.__init__(self, block_name, ["IN1"], ["OUT1"])
-        self.parent = parent
+        self._parent = parent
 
     def compute(self, curIteration):
         self.appendToSignal(self.getInputSignal(curIteration, "IN1").value)
@@ -270,7 +283,7 @@ class SequenceBlock(BaseBlock):
         self.__sequence = sequence
 
     def compute(self, curIteration):
-        if len(self.__sequence) < curIteration:
+        if len(self.__sequence) <= curIteration:
             raise ValueError("Sequence is not long enough")
         self.appendToSignal(self.__sequence[curIteration])
 
@@ -508,5 +521,15 @@ class CBD(BaseBlock):
         assert portBlock is not None
         return portBlock.getSignal("OUT1")
 
+    def getSignals(self):
+        res = {}
+        for port in super().getSignals().keys():
+            res[port] = self.getSignal(port)
+        return res
+
     def compute(self, curIteration):
         pass
+
+    def _rewind(self):
+        for block in self.getBlocks():
+            block._rewind()

src/CBD/lib/io.py

@@ -99,7 +99,7 @@ class ReadCSV(BaseBlock):
 		self.repeat = repeat
 		self.time_col = time_col
 		self.__read_file(file_name, **kwargs)
-		BaseBlock.__init__(self, block_name, [], list(self.data.keys()))
+		BaseBlock.__init__(self, block_name, [], [x for x in self.data.keys() if x != time_col])
 		self.hold = hold
 		self.index = -1
 
@@ -153,7 +153,7 @@ class ReadCSV(BaseBlock):
 				return
 			if self.index < 0:
 				self.index = 0
-		if self.index < 0:
+		if self.index < 0 or self.index >= L:
 			next_time = last_time = 0.0
 		else:
 			next_time = last_time = self.data[self.time_col][self.index]
@@ -168,6 +168,10 @@ class ReadCSV(BaseBlock):
 			val = Interpolation.interpolate(a, b, time, self.hold)
 			self.appendToSignal(val, col)
 
+	def _rewind(self):
+		BaseBlock._rewind(self)
+		self.index -= 1
+
 
 class WriteCSV(BaseBlock):
 	"""

src/CBD/lib/std.py

@@ -294,14 +294,18 @@ class DelayBlock(BaseBlock):
 
 class TimeBlock(BaseBlock):
 	"""
-	Outputs the current time of the simulation
+	Outputs the current time of the simulation.
+
+	:Ports:
+		- :code:`OUT1`:     The real simulation time
+		- :code`relative`:  The simulation time since simulation start
 	"""
 	def __init__(self, block_name):
-		BaseBlock.__init__(self, block_name, [], ["real", "relative"])
+		BaseBlock.__init__(self, block_name, [], ["OUT1", "relative"])
 
 	def compute(self, curIteration):
 		time = self.getClock().getTime()
-		self.appendToSignal(time, "real")
+		self.appendToSignal(time, "OUT1")
 		self.appendToSignal(time - self.getClock().getStartTime(), "relative")
 
 

src/CBD/simulator.py

@@ -33,6 +33,9 @@ class Clock:
 		"""
 		return self.__time
 
+	def setTime(self, time):
+		self.__time = time
+
 	def getStartTime(self):
 		"""
 		Gets the starting simulation time.
@@ -45,6 +48,18 @@ class Clock:
 		"""
 		self.__time = self.__time + self.__delta_t
 
+	def _rewind(self):
+		"""
+        Rewinds the simulation clock to the previous iteration, assuming the delta
+        has not been changed.
+
+        Danger:
+            Normally, this function should only be used by the internal mechanisms
+            of the CBD simulator, **not** by a user. Using this function without a
+            full understanding of the simulator may result in undefined behaviour.
+        """
+		self.__time = self.__time - self.__delta_t
+
 	def setDeltaT(self, new_delta_t):
 		"""
 		Sets the delta in-between timesteps.
@@ -92,6 +107,9 @@ class Simulator:
 		# simulation data [dep graph, strong components, curIt]
 		self.__sim_data = [None, None, 0]
 
+		self.__stepsize_backend = Fixed(self.__deltaT)
+		self.__deltas = []
+
 		self.__threading_backend = None
 		self.__threading_backend_subsystem = Platform.PYTHON
 		self.__threading_backend_args = []
@@ -112,6 +130,9 @@ class Simulator:
 		# TODO: make this variable, given more solver implementations
 		self.__solver = GaussianJordanLinearSolver(self.__logger)
 
+	def setBackend(self, back):
+		self.__stepsize_backend = back
+
 	def run(self, term_time=None, start_time=0.0):
 		"""
 		Simulates the model.
@@ -248,6 +269,7 @@ class Simulator:
 		clock = self.getClock()
 		if clock is not None:
 			clock.setDeltaT(delta_t)
+		self.__stepsize_backend.delta_t = delta_t
 
 	def getDeltaT(self):
 		"""
@@ -262,6 +284,9 @@ class Simulator:
 		"""
 		return self.__deltaT
 
+	def setSimData(self, data):
+		self.__sim_data = data
+
 	def setRealTime(self, enabled=True, scale=1.0):
 		"""
 		Makes the simulation run in (scaled) real time.
@@ -416,6 +441,37 @@ class Simulator:
 		"""
 		self.__threading_backend.step(time)
 
+	def do_single_step(self):
+		curIt = self.__sim_data[2]
+		# Efficiency reasons: dep graph only changes at these times
+		#   in the given set of library blocks.
+		# TODO: Must be set to "every time" instead.
+		if curIt < 2 or self.__sim_data[0] is None:
+			self.__sim_data[0] = createDepGraph(self.model, curIt)
+			self.__sim_data[1] = self.__sim_data[0].getStrongComponents(curIt)
+		self.__computeBlocks(self.__sim_data[1], self.__sim_data[0], self.__sim_data[2])
+		self.__sim_data[2] += 1
+
+	def update_clock(self):
+		if self.__sim_data[2] > 0:
+			new_dt = self.__stepsize_backend.getNextStepSize(self)
+			self.setDeltaT(new_dt)
+			self.__deltas.append(new_dt)
+			self.getClock().step()
+
+	def step_back(self):
+		"""
+        Rewinds the simulator to the previous iteration.
+
+        Danger:
+            Normally, this function should only be used by the internal mechanisms
+            of the CBD simulator, **not** by a user. Using this function without a
+            full understanding of the simulator may result in undefined behaviour.
+        """
+		self.getClock()._rewind()
+		self.model._rewind()
+		self.__sim_data[2] -= 1
+
 	def getDurationLog(self):
 		"""
 		Get the list of timings for every iteration.
@@ -425,6 +481,9 @@ class Simulator:
 		"""
 		return self.__duration_log
 
+	def getDeltaLog(self):
+		return self.__deltas
+
 	def __realtimeWait(self):
 		"""
 		Wait until next realtime event.
@@ -467,16 +526,10 @@ class Simulator:
 				self.__finish()
 				break
 
+			self.update_clock()
+
 			before = time.time()
-			# Efficiency reasons: dep graph only changes at these times
-			#   in the given set of library blocks.
-			# TODO: Must be set to "every time" instead.
-			curIt = self.__sim_data[2]
-			if curIt < 2:
-				self.__sim_data[0] = createDepGraph(self.model, curIt)
-				self.__sim_data[1] = self.__sim_data[0].getStrongComponents(curIt)
-			self.__step(*self.__sim_data)
-			self.__sim_data[2] += 1
+			self.do_single_step()
 			self.__duration_log.append(time.time() - before)
 
 			if self.__threading_backend_subsystem == Platform.GLA:
@@ -487,19 +540,6 @@ class Simulator:
 				# Next event has been scheduled, kill this process
 				break
 
-	def __step(self, depGraph, sortedGraph, curIteration):
-		"""
-		Executes a single step on all blocks that must be computed.
-
-		Args:
-			depGraph:           A dependency graph.
-			sortedGraph:        The set of strong components.
-			curIteration (int): Current simulation iteration.
-		"""
-		self.__computeBlocks(sortedGraph, depGraph, curIteration)
-		self.setDeltaT(self.getDeltaT())
-		self.getClock().step()
-
 	def __computeBlocks(self, sortedGraph, depGraph, curIteration):
 		"""
 		Compute the new state of the model.
@@ -605,3 +645,5 @@ class Simulator:
 			raise ValueError("Invalid signal '%s' in Simulator." % name)
 		for evt in self.__events[name]:
 			evt(*args)
+
+from .stepsize import *

src/CBD/stepsize.py

@@ -0,0 +1,81 @@
+
+from .simulator import Simulator
+from copy import deepcopy
+
+class StepSize:
+	def __init__(self, delta_t):
+		self.delta_t = delta_t
+
+	def getNextStepSize(self, sim):
+		raise NotImplementedError()
+
+
+class Fixed(StepSize):
+	def getNextStepSize(self, sim):
+		return self.delta_t
+
+
+class Euler2(StepSize):
+	def __init__(self, delta_t, epsilon):
+		StepSize.__init__(self, delta_t)
+		self.epsilon = epsilon
+
+	def getNextStepSize(self, sim):
+		sim.setBackend(Fixed(self.delta_t))
+		sim.update_clock()
+		sim.do_single_step()
+		A1 = deepcopy(sim.model.getSignals())
+
+		sim.step_back()
+
+		sim.setDeltaT(self.delta_t / 2)
+		sim.update_clock()
+		sim.do_single_step()
+		A2 = deepcopy(sim.model.getSignals())
+
+		sim.step_back()
+		sim.setBackend(self)
+		sim.setDeltaT(self.delta_t)
+
+
+		r_max = float('-inf')
+		for port in A1.keys():
+			r_max = max(r_max, abs(self.inp(A1, port) - A2[port][-1].value))
+		r_max /= self.delta_t
+		if r_max > self.epsilon:
+			self.delta_t *= .9 * self.epsilon / r_max
+			self.delta_t = round(self.delta_t, 6)
+			# self.delta_t /= 2
+			return self.getNextStepSize(sim)
+
+		# No error => can we go larger?
+		sim.setBackend(Fixed(self.delta_t * 2))
+		sim.update_clock()
+		sim.do_single_step()
+		A3 = deepcopy(sim.model.getSignals())
+
+		sim.step_back()
+		sim.setBackend(self)
+		sim.setDeltaT(self.delta_t)
+
+		r_max = float('-inf')
+		for port in A1.keys():
+			r_max = max(r_max, abs(self.inp(A3, port) - A1[port][-1].value))
+		r_max /= self.delta_t
+		if r_max > self.epsilon:
+			return self.delta_t
+		if r_max < 1e-6:
+			self.delta_t *= 2
+			self.delta_t = round(self.delta_t, 6)
+		else:
+			self.delta_t *= .9 * self.epsilon / r_max
+			self.delta_t = round(self.delta_t, 6)
+		# self.delta_t /= 2
+		return self.getNextStepSize(sim)
+
+	@staticmethod
+	def inp(A1, port):
+		return A1[port][-2].value - (A1[port][-2].value - A1[port][-1].value) / 2
+
+
+class RungeKuttaFehlberg(StepSize): pass