add more test units; refactor filter function; add IV support part1

Author: rssalessio

README.md

@@ -1,4 +1,4 @@
-# PythonVRFT Library - Version 0.0.4
+# PythonVRFT Library - Version 0.0.5
 VRFT Adaptive Control Library written in Python. Aim of this library is to provide an implementation of the VRFT (Virtual Reference Feedback Tuning) algorithm.
 
 You can find the package also at the following [link](https://pypi.org/project/pythonvrft/)
@@ -25,7 +25,7 @@ pip install .
 ``` 
 
 ## Examples
-Examples are located in the examples/ folder. At the moment only 1 example is available.
+Examples are located in the examples/ folder. At the moment 2 examples are available. One without instrumental variables and one using instrumental variables.
 
 ## Tests
 To execute tests run the following command
@@ -34,11 +34,11 @@ python -m unittest
 ``` 
 
 ## Changelog
-- [**DONE - V0.0.2**][26.03.2017] Implement the basic VRFT algorithm (1 DOF. offline, linear controller, controller expressed as scalar product theta*f(z))
-- [**DONE - V0.0.3**][05.01.2020] Code refactoring and conversion to Python 3; Removed support for Python Control library.
-- [**In Progress - V0.0.4**][07.01.2020-] Add Documentation and Latex formulas
+- [**DONE**][V. 0.0.2][26.03.2017] Implement the basic VRFT algorithm (1 DOF. offline, linear controller, controller expressed as scalar product theta*f(z))
+- [**DONE**][V. 0.0.3][05.01.2020] Code refactoring and conversion to Python 3; Removed support for Python Control library.
+- [**DONE**][v. 0.0.5][07.01.2020] Add Instrumental Variables (IVs) Support
+- [**In Progress**][-][07.01.2020-] Add Documentation and Latex formulas
 - [**TODO**] Add MIMO Support
-- [**TODO**] Add IV Support
 - [**TODO**] Generalize to other kind of controllers (e.g., neural nets)
 
 ## Citations

setup.py

@@ -6,7 +6,7 @@
 
 
 setup(name = 'pythonvrft',
-    version = '0.0.4',
+    version = '0.0.5',
     description = 'VRFT Python Library',
     long_description = long_description,
     long_description_content_type = 'text/markdown',

test/test_iddata.py

@@ -1,3 +1,8 @@
+# test_iddata.py - Unittest for the iddata object
+#
+# Code author: [Alessio Russo - alessior@kth.se]
+# Last update: 07th January 2020, by alessior@kth.se
+#
 # Copyright [2017-2020] [Alessio Russo - alessior@kth.se]  
 # This file is part of PythonVRFT.
 # PythonVRFT is free software: you can redistribute it and/or modify
@@ -10,33 +15,32 @@
 # You should have received a copy of the GNU General Public License
 # along with PythonVRFT.  If not, see <http://www.gnu.org/licenses/>.
 #
-# Code author: [Alessio Russo - alessior@kth.se]
-# Last update: 06th January 2020, by alessior@kth.se
-#
+
 
 from unittest import TestCase
 from vrft.iddata import iddata
 import numpy as np
+import scipy.signal as scipysig
 
 class TestIDData(TestCase):
     def test_type(self):
         a = iddata(0.0, 0.0, 0.0)
         with self.assertRaises(ValueError):
-            a.checkData()
+            a.check()
 
         a =  iddata(0.0, [1], 0.0)
         with self.assertRaises(ValueError):
-            a.checkData()
+            a.check()
 
         a = iddata(np.zeros(10), 1, 0.0)
         with self.assertRaises(ValueError):
-            a.checkData()
+            a.check()
 
         a = iddata([0 for i in range(10)], [0 for i in range(10)], 1.0)
-        self.assertTrue(a.checkData())
+        self.assertTrue(a.check())
 
         a = iddata(np.zeros(10), np.zeros(10), 1.0)
-        self.assertTrue(a.checkData())
+        self.assertTrue(a.check())
 
     def test_size(self):
         a = iddata(np.zeros(10), np.zeros(10), 0.0)
@@ -51,26 +55,124 @@ def test_size(self):
 
         a = iddata(np.zeros(10), np.zeros(9), 0.0)
         with self.assertRaises(ValueError):
-            a.checkData()
+            a.check()
 
         a = iddata(np.zeros(8), np.zeros(9), 0.0)
         with self.assertRaises(ValueError):
-            a.checkData()
+            a.check()
 
 
     def test_sampling_time(self):
         a = iddata(np.zeros(10), np.zeros(10), 0.0)
         with self.assertRaises(ValueError):
-            a.checkData()
+            a.check()
 
         a = iddata(np.zeros(10), np.zeros(10), 1e-9)
         with self.assertRaises(ValueError):
-            a.checkData()
+            a.check()
 
         a = iddata(np.zeros(10), np.zeros(10), -0.1)
         with self.assertRaises(ValueError):
-            a.checkData()
+            a.check()
+
+        a = iddata(np.zeros(10), np.zeros(10), 0.1)
+        self.assertTrue(a.check())
 
+    def test_copy(self):
         a = iddata(np.zeros(10), np.zeros(10), 0.1)
-        self.assertTrue(a.checkData())
+        b = a.copy()
+        self.assertTrue(a.check())
+        self.assertTrue(b.check())
 
+        self.assertTrue(np.all(a.y == b.y))
+        self.assertTrue(np.all(a.u == b.u))
+        self.assertTrue(np.all(a.y0 == b.y0))
+        self.assertTrue(a.ts == b.ts)
+
+    def test_filter(self):
+        a = iddata(np.zeros(10), np.zeros(10), 0.1)
+        L = scipysig.dlti([1], [1], dt=0.1)
+        b = a.copy()
+        a.filter(L)
+        self.assertTrue(np.all(a.y == b.y))
+        self.assertTrue(np.all(a.u == b.u))
+        self.assertTrue(np.all(a.y0 == b.y0))
+        self.assertTrue(a.ts == b.ts)
+
+    def test_split(self):
+        n = 10
+        a = iddata(np.random.normal(size=n), np.random.normal(size=n), 0.1)
+
+        b, c = a.split()
+
+        self.assertTrue(b.y.size == c.y.size)
+        self.assertTrue(b.u.size == c.u.size)
+        self.assertTrue(b.y0 == c.y0)
+        self.assertTrue(b.ts == c.ts)
+        self.assertTrue(b.ts == a.ts)
+        self.assertTrue(np.all(b.y == a.y[:n//2]))
+        self.assertTrue(np.all(b.u == a.u[:n//2]))
+        self.assertTrue(b.y0 == a.y0)
+
+        self.assertTrue(np.all(c.y == a.y[n//2:]))
+        self.assertTrue(np.all(c.u == a.u[n//2:]))
+        self.assertTrue(c.y0 == a.y0)
+
+        y0 = [-1, 2]
+        a = iddata(np.random.normal(size=n), np.random.normal(size=n), 0.1, y0)
+        n0 = len(y0)
+        n1 = (n + n0) // 2
+        b, c = a.split()
+
+        self.assertTrue(b.y.size == c.y.size)
+        self.assertTrue(b.u.size == c.u.size)
+        self.assertTrue(b.ts == c.ts)
+        self.assertTrue(b.ts == a.ts)
+        self.assertTrue(np.all(b.y == a.y[:n1 - n0]))
+        self.assertTrue(np.all(b.u == a.u[:n1 - n0]))
+        self.assertTrue(np.all(b.y0 == a.y0))
+
+        self.assertTrue(np.all(c.y == a.y[n1:]))
+        self.assertTrue(np.all(c.u == a.u[n1:]))
+        self.assertTrue(np.all(c.y0 == a.y[n1 - n0:n1]))
+
+
+        y0 = [-1, 2]
+        n = 9
+        a = iddata(np.random.normal(size=n), np.random.normal(size=n), 0.1, y0)
+        n0 = len(y0)
+        n -= 1
+        n1 = (n + n0) // 2
+        b, c = a.split()
+
+        self.assertTrue(b.y.size == c.y.size)
+        self.assertTrue(b.u.size == c.u.size)
+        self.assertTrue(b.ts == c.ts)
+        self.assertTrue(b.ts == a.ts)
+        self.assertTrue(np.all(b.y == a.y[:n1 - n0]))
+        self.assertTrue(np.all(b.u == a.u[:n1 - n0]))
+        self.assertTrue(np.all(b.y0 == a.y0))
+
+        self.assertTrue(np.all(c.y == a.y[n1:n]))
+        self.assertTrue(np.all(c.u == a.u[n1:n]))
+        self.assertTrue(np.all(c.y0 == a.y[n1 - n0:n1]))
+
+        y0 = [-1]
+        n = 10
+        a = iddata(np.random.normal(size=n), np.random.normal(size=n), 0.1, y0)
+        n0 = len(y0)
+        n -= 1
+        n1 = (n + n0) // 2
+        b, c = a.split()
+
+        self.assertTrue(b.y.size == c.y.size)
+        self.assertTrue(b.u.size == c.u.size)
+        self.assertTrue(b.ts == c.ts)
+        self.assertTrue(b.ts == a.ts)
+        self.assertTrue(np.all(b.y == a.y[:n1 - n0]))
+        self.assertTrue(np.all(b.u == a.u[:n1 - n0]))
+        self.assertTrue(np.all(b.y0 == a.y0))
+
+        self.assertTrue(np.all(c.y == a.y[n1:n]))
+        self.assertTrue(np.all(c.u == a.u[n1:n]))
+        self.assertTrue(np.all(c.y0 == a.y[n1 - n0:n1]))

test/test_reference.py

@@ -1,4 +1,8 @@
-# Copyright [2017-2020] [Alessio Russo - alessior@kth.se]  
+# test_reference.py - Unittest for virtual reference algorithm
+#
+# Code author: [Alessio Russo - alessior@kth.se]
+# Last update: 07th January 2020, by alessior@kth.se
+#
 # This file is part of PythonVRFT.
 # PythonVRFT is free software: you can redistribute it and/or modify
 # it under the terms of the GNU General Public License as published by
@@ -10,9 +14,6 @@
 # You should have received a copy of the GNU General Public License
 # along with PythonVRFT.  If not, see <http://www.gnu.org/licenses/>.
 #
-# Code author: [Alessio Russo - alessior@kth.se]
-# Last update: 06th January 2020, by alessior@kth.se
-#
 
 from unittest import TestCase
 from vrft.iddata import *

test/test_utils.py

@@ -1,3 +1,8 @@
+# test_utils.py - Unittest for utilities
+#
+# Code author: [Alessio Russo - alessior@kth.se]
+# Last update: 07th January 2020, by alessior@kth.se
+#
 # Copyright [2017-2020] [Alessio Russo - alessior@kth.se]  
 # This file is part of PythonVRFT.
 # PythonVRFT is free software: you can redistribute it and/or modify
@@ -10,14 +15,13 @@
 # You should have received a copy of the GNU General Public License
 # along with PythonVRFT.  If not, see <http://www.gnu.org/licenses/>.
 #
-# Code author: [Alessio Russo - alessior@kth.se]
-# Last update: 06th January 2020, by alessior@kth.se
-#
+
 
 from unittest import TestCase
 from vrft.utils import *
 from vrft.extended_tf import ExtendedTF
 from vrft.vrft_algo import virtualReference
+from vrft.iddata import iddata
 import numpy as np
 import scipy.signal as scipysig
 
@@ -33,7 +37,7 @@ def test_deconvolve(self):
         y = y[:, 0]
         data = iddata(y, u, t_step, [0])
         r1, _ = virtualReference(data, sys.num, sys.den)
-        r2 = deconvolve_signal(sys, data.y, data.ts)
+        r2 = deconvolve_signal(sys, data.y)
         self.assertTrue(np.linalg.norm(r2-r1[:r2.size], np.infty) <  1e-3)
 
 

test/test_vrft.py

@@ -1,3 +1,8 @@
+# test_vrft.py - Unittest for VRFT
+#
+# Code author: [Alessio Russo - alessior@kth.se]
+# Last update: 07th January 2020, by alessior@kth.se
+#
 # Copyright [2017-2020] [Alessio Russo - alessior@kth.se]  
 # This file is part of PythonVRFT.
 # PythonVRFT is free software: you can redistribute it and/or modify
@@ -10,9 +15,7 @@
 # You should have received a copy of the GNU General Public License
 # along with PythonVRFT.  If not, see <http://www.gnu.org/licenses/>.
 #
-# Code author: [Alessio Russo - alessior@kth.se]
-# Last update: 06th January 2020, by alessior@kth.se
-#
+
 
 from unittest import TestCase
 from vrft.iddata import *
@@ -37,13 +40,48 @@ def test_vrft(self):
         data = iddata(y,u,t_step,[0])
 
         refModel = ExtendedTF([0.2], [1, -0.8], dt=t_step)
+        prefilter = refModel * (1-refModel)
+
+        control = [ExtendedTF([1], [1,0], dt=t_step),
+                ExtendedTF([1], [1,0,0], dt=t_step),
+                ExtendedTF([1], [1,0,0,0], dt=t_step),
+                ExtendedTF([1, 0], [1,1], dt=t_step)]
+
+        theta1, _, loss1, _ = compute_vrft(data, refModel, control, prefilter)
+        theta2, _, loss2, _ = compute_vrft([data], refModel, control, prefilter)
+        theta3, _, loss3, _ = compute_vrft([data, data], refModel, control, prefilter)
+        
+        self.assertTrue(np.isclose(loss1, loss2))
+        self.assertTrue(np.isclose(loss1, loss3))
+        self.assertTrue(np.linalg.norm(theta1-theta2)<1e-15)
+        self.assertTrue(np.linalg.norm(theta1-theta3)<1e-15)
+
+    def test_iv(self):
+        t_start = 0
+        t_end = 10
+        t_step = 1e-2
+        t = np.arange(t_start, t_end, t_step)
+        u = np.ones(len(t)).tolist()
+
+        num = [0.1]
+        den = [1, -0.9]
+        sys = scipysig.TransferFunction(num, den, dt=t_step)
+        t, y = scipysig.dlsim(sys, u, t)
+        y = y.flatten() + np.random.normal(size=t.size)
+        data = iddata(y,u,t_step,[0])
+
+        refModel = ExtendedTF([0.2], [1, -0.8], dt=t_step)
+        prefilter = refModel * (1-refModel)
 
         control = [ExtendedTF([1], [1,0], dt=t_step),
                 ExtendedTF([1], [1,0,0], dt=t_step),
                 ExtendedTF([1], [1,0,0,0], dt=t_step),
                 ExtendedTF([1, 0], [1,1], dt=t_step)]
 
-        theta, _, _, loss, _ = compute_vrft(data, refModel, control, refModel * (1-refModel))
+        ##import pdb
+        #pdb.set_trace()
+        
+        # theta, _, loss, _ = compute_vrft(data, refModel, control, prefilter, iv=True)
 
 
 

vrft/extended_tf.py

@@ -76,7 +76,7 @@ def __add__(self,other):
         if type(other) in [int, float]:
             return ExtendedTF(polyadd(self.num, self.den*other), self.den, dt=self._dt)
         if type(other) in [TransFun, type(self)]:
-            if np.all(self.den == other.den):
+            if len(self.den) == len(other.den) and np.all(self.den == other.den):
                 numer = polyadd(self.num, other.num)
                 denom = self.den
             else:
@@ -88,7 +88,7 @@ def __sub__(self, other):
         if type(other) in [int, float]:
             return ExtendedTF(polyadd(self.num, -self.den*other), self.den, dt=self._dt)
         if type(other) in [TransFun, type(self)]:
-            if np.all(self.den == other.den):
+            if len(self.den) == len(other.den) and np.all(self.den == other.den):
                 numer = polyadd(self.num, -other.num)
                 denom = self.den
             else:
@@ -100,7 +100,7 @@ def __rsub__(self, other):
         if type(other) in [int, float]:
             return ExtendedTF(polyadd(-self.num, self.den*other), self.den, dt=self._dt)
         if type(other) in [TransFun, type(self)]:
-            if np.all(self.den == other.den):
+            if len(self.den) == len(other.den) and np.all(self.den == other.den):
                 numer = polyadd(self.num, -other.num)
                 denom = self.den
             else: