Runge-Kutta: The Sequel

examples/SinGen/mpl.py

@@ -31,7 +31,7 @@ from CBD.simulator import Simulator
 # manager.connect('sin', 'update', lambda d, axis=ax: axis.set_xlim(follow(d[0], 10.0, lower_bound=0.0)))
 
 sim = Simulator(sinGen)
-sim.setVerbose("test.log")
+sim.setVerbose("test.txt")
 # sim.setRealTime()
 sim.setDeltaT(0.1)
 sim.run(20.0)

src/CBD/CBD.py

@@ -36,6 +36,36 @@ class BaseBlock:
         # In wich CBD the baseblock is situated
         self._parent = None
 
+    def addInputPort(self, name):
+        """
+        Adds an input port if there is no port with the given name.
+
+        Args:
+            name (str): The name of the port.
+        """
+        if name not in self.__nameLinks:
+            self.__nameLinks.append(name)
+
+    def addOutputPort(self, name):
+        """
+        Adds an output port if there is no port with the given name.
+
+        Args:
+            name (str): The name of the port.
+        """
+        if name not in self.__signals:
+            self.__signals[name] = []
+
+    def clone(self):
+        """
+        (Deep) copies the current block, ignoring all connections or links
+        that were set on this block.
+        """
+        other = deepcopy(self)
+        other._linksIn.clear()
+        other._parent = None
+        return other
+
     def getBlockName(self):
         """
         Gets the name of the block.
@@ -72,6 +102,9 @@ class BaseBlock:
         """
         return self._linksIn
 
+    def getInputPortNames(self):
+        return self.__nameLinks
+
     def getOutputNameOfInput(self, inputBlock):
         """
         Gets the name of the output port in the :code:`inputBlock` that is linked to this block.
@@ -325,16 +358,37 @@ class CBD(BaseBlock):
         for output_port in output_ports:
             self.addBlock(OutputPortBlock(output_port, self))
 
+    def addInputPort(self, name):
+        BaseBlock.addInputPort(self, name)
+        self.addBlock(InputPortBlock(name, self))
+
+    def addOutputPort(self, name):
+        BaseBlock.addOutputPort(self, name)
+        self.addBlock(OutputPortBlock(name, self))
+
+    def clone(self):
+        other = BaseBlock.clone(self)
+        other.setClock(deepcopy(self.getClock()))
+        other.clearBlocks()
+        for block in self.getBlocks():
+            other.addBlock(block.clone())
+        for block in self.getBlocks():
+            for name_input, link in block.getLinksIn().items():
+                other.addConnection(link.block.getBlockName(), block.getBlockName(), name_input, link.output_port)
+        return other
+
     def getTopCBD(self):
         """
         Finds the highest-level :class:`CBD` instance.
         """
         return self if self._parent is None else self._parent.getTopCBD()
 
-    def flatten(self, parent=None):
+    def flatten(self, parent=None, ignore=None):
         """
         Flatten the CBD parent and call flatten recursively for CBD's in this CBD
         """
+        if ignore is None:
+            ignore = []
         blocksToRemove = []
         blocksToAdd = []
 
@@ -371,8 +425,8 @@ class CBD(BaseBlock):
                                 b._linksIn[portname] = InputLink(wb, "OUT1")
 
         for childBlock in self.__blocks:
-            if isinstance(childBlock, CBD):
-                childBlock.flatten(self)
+            if isinstance(childBlock, CBD) and not isinstance(childBlock, tuple(ignore)):
+                childBlock.flatten(self, ignore)
                 blocksToRemove.append(childBlock)
 
         # Delete blocksToRemove
@@ -403,6 +457,23 @@ class CBD(BaseBlock):
         """
         return self.__blocksDict[name]
 
+    def hasBlock(self, name):
+        """
+        Checks if the CBD has a block with the given name.
+
+        Args:
+            name (str): The name of the block to check.
+        """
+        return name in self.__blocksDict
+
+    def clearBlocks(self):
+        """
+        Clears the block information. Calling this function will
+        "empty" the current block.
+        """
+        self.__blocks.clear()
+        self.__blocksDict.clear()
+
     def getClock(self):
         """
         Gets the current simulation clock.

src/CBD/CBDDraw.py

@@ -1,8 +1,10 @@
 """
 Useful drawing function to easily draw a CBD model in Graphviz.
 """
+import re
+from CBD.lib.std import *
+from CBD.CBD import CBD
 
-from .CBD import *
 def draw(cbd, filename, colors=None):
 	"""
 	Output :class:`CBD` as a graphviz script to filename.
@@ -26,6 +28,9 @@ digraph graphname {
 	if colors is None:
 		colors = {}
 
+	def formatName(name):
+		return re.sub("[^a-zA-Z_0-9]", "", name)
+
 	def writeBlock(block):
 		"""
 		Writes a block to graphviz.
@@ -46,7 +51,7 @@ digraph graphname {
 		if block.getBlockName() in colors:
 			col = ", color=\"{0}\", fontcolor=\"{0}\"".format(colors[block.getBlockName()])
 
-		write("{b} [label=\"{lbl}\"{shape}{col}];\n".format(b=block.getBlockName(),
+		write("{b} [label=\"{lbl}\"{shape}{col}];\n".format(b=formatName(block.getBlockName()),
 			lbl=label,
 			shape=shape,
 			col=col))
@@ -55,16 +60,16 @@ digraph graphname {
 	for block in cbd.getBlocks():
 		writeBlock(block)
 		for (name, other) in  block.getLinksIn().items():
-			label = ""
+			label = name
 
-			if not name.startswith("IN"):
-				label=name
+			# if not name.startswith("IN"):
+			# 	label=name
 
 			if not other.output_port.startswith("OUT"):
 				label = label + " / " + other.output_port
 
-			write("{a} -> {b} [label=\"{lbl}\"];\n".format(a=other.block.getBlockName(),
-				b=block.getBlockName(),
+			write("{a} -> {b} [label=\"{lbl}\"];\n".format(a=formatName(other.block.getBlockName()),
+				b=formatName(block.getBlockName()),
 				lbl=label))
 
 	write("\n}")

src/CBD/lib/std.py

@@ -8,7 +8,7 @@ import math
 __all__ = ['ConstantBlock', 'NegatorBlock', 'InverterBlock', 'AdderBlock', 'ProductBlock', 'ModuloBlock',
            'RootBlock', 'AbsBlock', 'IntBlock', 'ClampBlock', 'GenericBlock', 'MultiplexerBlock', 'LessThanBlock',
            'EqualsBlock', 'LessThanOrEqualsBlock', 'NotBlock', 'OrBlock', 'AndBlock', 'DelayBlock', 'TimeBlock',
-           'LoggingBlock', 'AddOneBlock', 'DerivatorBlock', 'IntegratorBlock']
+           'LoggingBlock', 'AddOneBlock', 'DerivatorBlock', 'IntegratorBlock', 'DeltaTBlock']
 
 class ConstantBlock(BaseBlock):
 	"""
@@ -60,14 +60,28 @@ class InverterBlock(BaseBlock):
 
 class AdderBlock(BaseBlock):
 	"""
-	The adderblock will add the 2 inputs
+	The adderblock will add all the inputs.
+
+	Args:
+		block_name (str):       The name of the block.
+		numberOfInputs (int):   The amount of input ports to set.
 	"""
-	def __init__(self, block_name):
-		BaseBlock.__init__(self, block_name, ["IN1", "IN2"], ["OUT1"])
+	def __init__(self, block_name, numberOfInputs=2):
+		BaseBlock.__init__(self, block_name, ["IN%d" % (x+1) for x in range(numberOfInputs)], ["OUT1"])
+		self.__numberOfInputs = numberOfInputs
 
 	def compute(self, curIteration):
 		# TO IMPLEMENT
-		self.appendToSignal(self.getInputSignal(curIteration, "IN1").value + self.getInputSignal(curIteration, "IN2").value)
+		result = 0
+		for i in range(1, self.__numberOfInputs+1):
+			result += self.getInputSignal(curIteration, "IN%d"%i).value
+		self.appendToSignal(result)
+
+	def getNumberOfInputs(self):
+		"""
+		Gets the total number of input ports.
+		"""
+		return self.__numberOfInputs
 
 
 class ProductBlock(BaseBlock):

src/CBD/preprocessing/rungekutta.py

@@ -0,0 +1,216 @@
+"""
+This module contains all the logic for Runge-Kutta preprocessing.
+"""
+from CBD.CBD import CBD
+from CBD.lib.std import *
+
+class RKPreprocessor:
+	def __init__(self, tableau):
+		self._tableau = tableau
+
+	def preprocess(self, original):
+		# 1. Detect all IVPs and group them in their own blocks
+		model = original.clone()
+		model.flatten(ignore=[IntegratorBlock])
+		blocks = model.getBlocks()
+		IVP = []
+		for block in blocks:
+			if isinstance(block, IntegratorBlock):
+				i = len(IVP)
+				ivp = CBD("IVP-%d" % i, ["time"], ["OUT1"])
+				collection = self.collect(block, ["IN1"], [IntegratorBlock, DelayBlock, TimeBlock])
+				for child in collection:
+					ivp.addBlock(child.clone())
+				for child in collection:
+					for name_input, link in child.getLinksIn().items():
+						lbn = link.block.getBlockName()
+						lop = link.output_port
+						if not ivp.hasBlock(lbn):
+							if link.block.getBlockType() == "TimeBlock":
+								lbn = 'time'
+							else:
+								lbn = name_input + "-" + child.getBlockName()
+								ivp.addInputPort(lbn)
+							lop = None
+						ivp.addConnection(lbn, child.getBlockName(), name_input, lop)
+				fin = block.getBlockConnectedToInput("IN1")
+				ivp.addConnection(fin.block.getBlockName(), "OUT1", None, fin.output_port)
+				IVP.append((block, ivp))
+
+		# 2. Foreach IVP: create an RK-model, based on the given tableau
+		RKs = []
+		for block, ivp in IVP:
+			RKs.append((block, self.create_RK(ivp)))
+
+		# 3. Substitute the RK-model collection in the original CBD
+		# TODO: multiple outputs
+		# TODO: multiple IVP?
+		outputs = original.getSignals().keys()
+		new_model = CBD(original.getBlockName(), [], outputs)
+		new_model.addBlock(TimeBlock("time"))
+		new_model.addBlock(DeltaTBlock("delta_t"))
+		for integ, rk in RKs:
+			rkname = rk.getBlockName()
+			new_model.addBlock(rk)
+
+			# TODO: more complex IC?
+			ic = integ.getBlockConnectedToInput("IC").block
+			new_model.addBlock(ic.clone())
+			new_model.addConnection(ic.getBlockName(), rkname, "IC")
+
+			new_model.addConnection("delta_t", rkname, "h")
+			new_model.addConnection("time", rkname, "t")
+			for outp in outputs:
+				conn = original.getBlockByName(outp).getBlockConnectedToInput("IN1").block
+				if conn.getBlockName() == integ.getBlockName():
+					new_model.addConnection(rkname, outp)
+		return new_model
+
+
+	def collect(self, start, sport=None, finish=None):
+		if finish is None:
+			finish = []
+		collection = [x[1].block for x in start.getLinksIn().items() if \
+		              ((sport is not None and x[0] in sport) or (sport is None)) and not isinstance(x[1].block, tuple(finish))]
+		n_collection = [x.getBlockName() for x in collection]
+		for block in collection:
+			ccoll = self.collect(block, None, finish)
+			for child in ccoll:
+				cname = child.getBlockName()
+				if cname not in n_collection:
+					n_collection.append(cname)
+					collection.append(child)
+		return collection
+
+	def create_RK(self, f):
+		RK = CBD("RK", ["h", "t", "IC"], ["OUT1"])
+		fy = [x for x in f.getInputPortNames() if x != 'time']
+		weights = self._tableau.getWeights()[0]
+		RK.addBlock(AdderBlock("RKSum", len(weights)))
+		RK.addBlock(DelayBlock("delay"))
+		for i in range(len(weights)):
+			j = i + 1
+			RK.addBlock(self.create_K(j, f.clone()))
+			RK.addBlock(ProductBlock("Mult_%d" % j))
+			RK.addBlock(ConstantBlock("B_%d" % j, weights[i]))
+			RK.addConnection("h", "RK-K_%d" % j, "h")
+			RK.addConnection("t", "RK-K_%d" % j, "t")
+			for y in fy:
+				RK.addConnection("delay", "RK-K_%d" % j, y)
+			RK.addConnection("B_%d" % j, "Mult_%d" % j)
+			RK.addConnection("RK-K_%d" % j, "Mult_%d" % j)
+			RK.addConnection("Mult_%d" % j, "RKSum")
+			for s in range(i):
+				RK.addConnection("RK-K_%d" % (s+1), "RK-K_%d" % j, "k_%d" % (s+1))
+
+		# Initial Conditions
+		RK.addBlock(NegatorBlock("neg"))
+		RK.addBlock(AdderBlock("ICSum"))
+		RK.addConnection("RKSum", "neg")
+		RK.addConnection("neg", "ICSum")
+		RK.addConnection("IC", "ICSum")
+		RK.addConnection("ICSum", "delay", "IC")
+
+		# Loop
+		RK.addBlock(AdderBlock("YSum"))
+		RK.addConnection("delay", "YSum")
+		RK.addConnection("YSum", "delay", "IN1")
+		RK.addConnection("RKSum", "YSum")
+		RK.addConnection("delay", "OUT1")
+
+		return RK
+
+	def create_K(self, s, f):
+		input_ports = ["h", "t"] + ["k_%d" % (i+1) for i in range(s-1)]
+		fy = [x for x in f.getInputPortNames() if x != "time"]
+		input_ports += fy
+		K = CBD("RK-K_%d" % s, input_ports, ["OUT1"])
+		K.addBlock(f)
+
+		# Time parameter
+		K.addBlock(ConstantBlock("C", self._tableau.getNodes()[s-1]))
+		K.addBlock(ProductBlock("CMult"))
+		K.addBlock(AdderBlock("CSum"))
+		K.addConnection("h", "CMult")
+		K.addConnection("C", "CMult")
+		K.addConnection("t", "CSum")
+		K.addConnection("CMult", "CSum")
+		K.addConnection("CSum", f.getBlockName(), "time")
+
+		# Y parameters
+		if s - 1 > 0:
+			K.addBlock(AdderBlock("KSum", s - 1))
+			for i in range(s-1):
+				j = i + 1
+				K.addBlock(ConstantBlock("A_%d" % j, self._tableau.getA(s-1, j)))
+				K.addBlock(ProductBlock("Mult_%d" % j))
+				K.addConnection("A_%d" % j, "Mult_%d" % j)
+				K.addConnection("k_%d" % j, "Mult_%d" % j)
+				K.addConnection("Mult_%d" % j, "KSum")
+			for y in fy:
+				K.addInputPort(y)
+				K.addBlock(AdderBlock("YSum-%s" % y))
+				K.addConnection(y, "YSum-%s" % y)
+				K.addConnection("KSum", "YSum-%s" % y)
+				K.addConnection("YSum-%s" % y, f.getBlockName(), y)
+		else:
+			for y in fy:
+				K.addInputPort(y)
+				K.addConnection(y, f.getBlockName(), y)
+
+		# Finishing Up
+		K.addBlock(ProductBlock("FMult"))
+		K.addConnection("h", "FMult")
+		K.addConnection(f.getBlockName(), "FMult")
+		K.addConnection("FMult", "OUT1")
+
+		return K
+
+
+if __name__ == '__main__':
+	from CBD.stepsize import ButcherTableau as BT
+	from CBD.CBDDraw import draw
+	DELTA_T = 0.1
+
+	class Test(CBD):
+		def __init__(self, block_name):
+			CBD.__init__(self, block_name, input_ports=[], output_ports=['y'])
+
+			# Create the Blocks
+			self.addBlock(IntegratorBlock("int"))
+			self.addBlock(ConstantBlock("IC", value=(0)))
+			self.addBlock(ProductBlock("mult"))
+			self.addBlock(AdderBlock("sum"))
+			self.addBlock(ConstantBlock("one", value=(1)))
+			self.addBlock(ConstantBlock("time", value=(DELTA_T)))
+
+			# Create the Connections
+			self.addConnection("IC", "int", output_port_name='OUT1', input_port_name='IC')
+			self.addConnection("int", "mult", output_port_name='OUT1', input_port_name='IN1')
+			self.addConnection("int", "mult", output_port_name='OUT1', input_port_name='IN2')
+			self.addConnection("int", "y", output_port_name='OUT1')
+			self.addConnection("mult", "sum", output_port_name='OUT1', input_port_name='IN2')
+			self.addConnection("one", "sum", output_port_name='OUT1', input_port_name='IN1')
+			self.addConnection("sum", "int", output_port_name='OUT1', input_port_name='IN1')
+			self.addConnection("time", "int", output_port_name='OUT1', input_port_name='delta_t')
+
+	prep = RKPreprocessor(BT.Heun())
+	model = prep.preprocess(Test("Test"))
+	draw(model.findBlock("RK.RK-K_2")[0], "test.dot")
+	# model = Test("Test")
+
+	# from CBD.simulator import Simulator
+	# sim = Simulator(model)
+	# sim.setDeltaT(0.1)
+	# sim.run(1.4)
+	#
+	# s = model.getSignal("y")
+	# L = len(s)
+	#
+	# print("+------------+------------+")
+	# print("|    TIME    |    VALUE   |")
+	# print("+------------+------------+")
+	# for i in range(L):
+	# 	t, v = s[i]
+	# 	print(f"| {t:10.7f} | {v:10.7f} |")
+	# print("+------------+------------+")

src/CBD/preprocessing/test.dot

@@ -0,0 +1,37 @@
+
+digraph graphname {
+ h [label="InputPortBlock (h)"];
+t [label="InputPortBlock (t)"];
+k_1 [label="InputPortBlock (k_1)"];
+IN1mult [label="InputPortBlock (IN1-mult)"];
+IN2mult [label="InputPortBlock (IN2-mult)"];
+OUT1 [label="OutputPortBlock (OUT1)"];
+FMult -> OUT1 [label="IN1"];
+IVP0 [label="CBD (IVP-0)",shape=Msquare];
+CSum -> IVP0 [label="time"];
+YSumIN1mult -> IVP0 [label="IN1-mult"];
+YSumIN2mult -> IVP0 [label="IN2-mult"];
+C [label="ConstantBlock (C)\n1"];
+CMult [label="ProductBlock (CMult)"];
+h -> CMult [label="IN1"];
+C -> CMult [label="IN2"];
+CSum [label="AdderBlock (CSum)"];
+t -> CSum [label="IN1"];
+CMult -> CSum [label="IN2"];
+KSum [label="AdderBlock (KSum)"];
+Mult_1 -> KSum [label="IN1"];
+A_1 [label="ConstantBlock (A_1)\n1"];
+Mult_1 [label="ProductBlock (Mult_1)"];
+A_1 -> Mult_1 [label="IN1"];
+k_1 -> Mult_1 [label="IN2"];
+YSumIN1mult [label="AdderBlock (YSum-IN1-mult)"];
+IN1mult -> YSumIN1mult [label="IN1"];
+KSum -> YSumIN1mult [label="IN2"];
+YSumIN2mult [label="AdderBlock (YSum-IN2-mult)"];
+IN2mult -> YSumIN2mult [label="IN1"];
+KSum -> YSumIN2mult [label="IN2"];
+FMult [label="ProductBlock (FMult)"];
+h -> FMult [label="IN1"];
+IVP0 -> FMult [label="IN2"];
+
+}

src/CBD/stepsize.py

@@ -279,6 +279,16 @@ class ButcherTableau:
 		"""
 		return self._weights
 
+	def getA(self, i, j):
+		"""
+		Obtains an element from the Runge-Kutta matrix.
+
+		Args:
+			i (int):    The row (1-indexed).
+			j (int):    The column (1-indexed).
+		"""
+		return self._matrix[i - 1][j - 1]
+
 	@staticmethod
 	def Heun():
 		r"""

src/CBD/tracers/color.py

@@ -38,7 +38,7 @@ class COLOR:
 		Args:
 			text (str): The string to remove color from.
 		"""
-		return re.sub(r"\[\d+m", "", text)
+		return re.sub(r"\033\[\d+m", "", text)
 
 	@staticmethod
 	def rainbow(text, colors=None):