diff --git a/src/ellipsoid.py b/src/ellipsoid.py
new file mode 100644
index 0000000000000000000000000000000000000000..43841807cb39537025b5b7cf32dd91b295ace418
--- /dev/null
+++ b/src/ellipsoid.py
@@ -0,0 +1,261 @@
+## This file is part of the simulative evaluation for the qronos observer abstractions.
+## Copyright (C) 2022-2023 Tim Rheinfels <tim.rheinfels@fau.de>
+## See https://gitlab.cs.fau.de/qronos-state-abstractions/simulation
+##
+## Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met:
+##
+## 1. Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer.
+##
+## 2. Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following disclaimer in the documentation and/or other materials provided with the distribution.
+##
+## 3. Neither the name of the copyright holder nor the names of its contributors may be used to endorse or promote products derived from this software without specific prior written permission.
+##
+## THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS “AS IS” AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+
+###
+### @file ellipsoid.py
+###
+### @brief Provides functionality for ellipsoidal sets
+###
+### @author Tim Rheinfels <tim.rheinfels@fau.de>
+###
+
+import json
+import matplotlib.pyplot as plot
+import numpy as np
+import scipy as sp
+import scipy.linalg
+
+import analysis
+from util import Serializable
+
+###
+### @brief Class representing a (possibly degenerate) ellipsoidal set
+###
+### @note The paper \cite ECRTS23-Rheinfels avoids the need to distinguish between
+### an ellipsoid's shape matrix @f$ P^{-1} @f$ and its inverse used to compute
+### the norms, i.e. @f$ \norm{x}_P^2 = x^T P x @f$. For the sake of brevity, it
+### refers to the latter as shape matrix. In contrast, this class (and the rest
+### of this package) uses the explicit distinction pointed out above.
+###
+class Ellipsoid(Serializable):
+
+ ###
+ ### @brief Constructor
+ ###
+ ### @param Pi Shape matrix @f$ P^{-1} \in \R^{n \times m}, P^{-1} \succeq 0 @f$ or None for zero-dimensional ellipsoid
+ ### @param c Center point @f$ c \in \R^n @f$ or None for centered ellipsoid
+ ###
+ def __init__(self, Pi=None, c=None):
+ # Promote matrices
+ if Pi is None:
+ Pi = np.zeros((0, 0))
+ elif isinstance(Pi, float) or isinstance(Pi, int):
+ Pi = np.array([[Pi]])
+
+ if c is None:
+ c = np.zeros(Pi.shape[1])
+ elif isinstance(c, float) or isinstance(c, int):
+ c = np.array([c])
+
+ # Check dimension
+ assert(c is None or Pi.shape[0] == len(c))
+
+ self._n = Pi.shape[0]
+ self._Pi = Pi
+
+ # Based on the dimension...
+ if self._n != 0:
+ # ... compute (pseudo-) inverse of Pi and its Cholesky (/LDL) decomposition
+ assert(analysis.is_spsd(Pi))
+ if analysis.is_spd(Pi):
+ self._P = np.linalg.inv(self._Pi)
+ self._Pi_sqrt = np.linalg.cholesky(self._Pi)
+ self._degenerate = False
+ else:
+ self._P = np.linalg.pinv(Pi)
+ L, D, perm = sp.linalg.ldl(self._Pi)
+ self._Pi_sqrt = L @ np.sqrt(D)
+ self._degenerate = True
+ else:
+ # ... or forget about it
+ self._P = self._Pi
+ self._Pi_sqrt = self._Pi
+ self._degenerate = False
+
+ self._c = c.flatten()
+ self._centered = (not np.any(self._c))
+
+ ###
+ ### @brief Getter for the dimension
+ ###
+ ### @returns Dimension @f$ n \in \N @f$
+ ###
+ @property
+ def n(self):
+ return self._n
+
+ ###
+ ### @brief Getter for the degenerate flag
+ ###
+ ### @returns Either
+ ### - True if ellipsoid is degenerate, i.e., @f$ P \succeq 0 \wedge P \nsucc 0 @f$ or
+ ### - False if ellipsoid is non-degenerate, i.e., @f$ P \succ 0 @f$
+ ###
+ @property
+ def degenerate(self):
+ return self._degenerate
+
+ ###
+ ### @brief Getter for the centered flag
+ ###
+ ### @returns Either
+ ### - True if ellipsoid is centered, i.e., @f$ c = 0 @f$ or
+ ### - False if ellipsoid is not centered, i.e., @f$ c \neq 0 @f$
+ ###
+ @property
+ def centered(self):
+ return self._centered
+
+ ###
+ ### @brief Getter for the shape matrix
+ ###
+ ### @returns Shape matrix @f$ P^{-1} \in \R^{n \times n}, P^{-1} \succeq 0 @f$
+ ###
+ @property
+ def Pi(self):
+ return self._Pi
+
+ ###
+ ### @brief Getter for the shape matrix's inverse or pseudo inverse in the degenerate case
+ ###
+ ### @returns Either
+ ### - Shape matrix's inverse @f$ P \in \R^{n \times n}, P \succ 0 @f$ in the non-degenerate case or
+ ### - shape matrix's pseudo inverse @f$ P^\dagger \in \R^{n \times n}, P^\dagger \succeq 0 @f$ in the degenerate case
+ ###
+ @property
+ def P(self):
+ return self._P
+
+ ###
+ ### @brief Getter for the shape matrix's Cholesky or semi-definite triangular decomposition in the degenerate case
+ ###
+ ### @returns Either
+ ### - Shape matrix's Cholesky decomposition @f$ P^{-\frac{T}{2}} \succ 0 @f$ in the non-degenerate case or
+ ### - shape matrix decomposed as @f$ L \sqrt{D} \succeq 0 @f$ with the LDL decomposition @f$ P^{-1} = LDL^T @f$
+ ###
+ @property
+ def Pi_sqrt(self):
+ return self._Pi_sqrt
+
+ ###
+ ### @brief Getter for the center position
+ ###
+ ### @returns Center position @f$ c \in \R^n @f$
+ ###
+ @property
+ def c(self):
+ return self._c
+
+ ###
+ ### @brief Computes the ellipsoid norm of a point @p x
+ ###
+ ### @param x Point @f$ x \in \R^n @f$ compute norm for
+ ###
+ ### @returns Ellipsoid norm of @p x as @f$ \norm{x}_P := \sqrt{(x - c)^T P (x - c)} @f$
+ ###
+ ### @note - For values less than one, @p x lies in the interior of the ellipsoid,
+ ### - for values equal to one, @p x lies on the surface of the ellipsoid, and
+ ### - for values greater than one, @p x lies outside of the ellipsoid.
+ ###
+ def norm(self, x):
+ assert(not self.degenerate)
+ xc = x - self.c
+ return np.sqrt(np.dot(xc, self.P @ xc))
+
+ ###
+ ### @brief Checks if a point @p x is an element of the ellipsoid
+ ###
+ ### @param x Point @f$ x \in \R^n @f$ to check
+ ###
+ ### @returns Either
+ ### - True if @p x is an element of the ellipsoid or
+ ### - False if @p x is not an element of the ellipsoid
+ ###
+ def contains(self, x):
+ return (self.norm(x) <= 1.0)
+
+ ###
+ ### @brief Uniformly draws a point from the ellipsoid's surface
+ ###
+ ### This functions uniformly draws a point from the ellipsoid's surface, i.e.
+ ###
+ ### @f[ x \sim \U(\partial \E(P, c)) @f]
+ ###
+ ### by first sampling uniformly from the unit sphere @f$ \B @f$ and
+ ### then applying the affine transform onto the ellipsoid's surface:
+ ###
+ ### With @f$ y \sim \U(\B) @f$ let @f$ x = P^{\frac{1}{2}} y + c \sim \U(\partial \E(P, c)). @f$
+ ###
+ ### @param rs RandomState to use
+ ###
+ ### @returns Point @f$ x @f$ uniformly sampled from the ellipsoid's surface
+ ###
+ def draw_from_surface(self, rs):
+ assert(isinstance(rs, np.random.Generator))
+ # Zero dimensions and completely degenerate ellipsoids are boring...
+ if self.n == 0 or np.all(self.Pi == 0):
+ return self.c
+
+ # Continue sampling (0,1)-normally distributed iid vectors...
+ while True:
+ y = rs.normal(0, 1, self.n)
+ n_y = np.linalg.norm(y)
+ # ... until norm is sufficiently large
+ if n_y > 1e-3:
+ # Division by norm yields sample on the unit sphere
+ y /= n_y
+
+ # Then apply the ellipsoidal transform
+ return self.Pi_sqrt @ y + self.c
+
+ ###
+ ### @brief Uniformly draws a point from the ellipsoid's volume
+ ###
+ ### This functions uniformly draws a point from the ellipsoid's volume, i.e.
+ ###
+ ### @f[ x \sim \U(\E(P, c)) @f]
+ ###
+ ### by first sampling uniformly from its surface and then applying a
+ ### scaling around the center such that the result is uniformly distributed
+ ### inside of the ellipsoid's volume:
+ ###
+ ### With @f$ y \sim \U(\partial \E(P, c)) @f$ and @f$ r \sim \U([0, 1]) @f$ let @f$ x = r^{\frac{1}{n}} (y-c) + c \sim \U(\E(P, c)). @f$
+ ###
+ ### @param rs RandomState to use
+ ###
+ ### @returns Point @f$ x @f$ uniformly sampled from the ellipsoid's volume
+ ###
+ def draw_from_volume(self, rs):
+ assert(isinstance(rs, np.random.Generator))
+ if self.n == 0 or np.all(self.Pi == 0):
+ return self.c
+ x = self.draw_from_surface(rs)
+ r = rs.uniform(0, 1) ** (1.0 / self.n)
+ return r*(x-self.c) + self.c
+
+ ###
+ ### @brief Converts the ellipsoid's data into a dictionary
+ ###
+ ### @returns Dictionary characterizing the ellipsoid
+ ###
+ ### @see util.Serializable
+ ###
+ def to_dict(self):
+ return {
+ 'Pi': self.Pi,
+ 'P': self.P,
+ 'c': self.c,
+ 'centered': self.centered,
+ 'degenerate': self.degenerate,
+ }
diff --git a/src/test/ellipsoid/__init__.py b/src/test/ellipsoid/__init__.py
new file mode 100644
index 0000000000000000000000000000000000000000..e69de29bb2d1d6434b8b29ae775ad8c2e48c5391
diff --git a/src/test/ellipsoid/constructor.py b/src/test/ellipsoid/constructor.py
new file mode 100644
index 0000000000000000000000000000000000000000..ac726a0f333807ee298006f5b73374671540ba6f
--- /dev/null
+++ b/src/test/ellipsoid/constructor.py
@@ -0,0 +1,151 @@
+## This file is part of the simulative evaluation for the qronos observer abstractions.
+## Copyright (C) 2022-2023 Tim Rheinfels <tim.rheinfels@fau.de>
+## See https://gitlab.cs.fau.de/qronos-state-abstractions/simulation
+##
+## Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met:
+##
+## 1. Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer.
+##
+## 2. Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following disclaimer in the documentation and/or other materials provided with the distribution.
+##
+## 3. Neither the name of the copyright holder nor the names of its contributors may be used to endorse or promote products derived from this software without specific prior written permission.
+##
+## THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS “AS IS” AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+
+###
+### @file test/ellipsoid/constructor.py
+###
+### @brief Provides unit tests for @ref ellipsoid.Ellipsoid.__init__
+###
+### @author Tim Rheinfels <tim.rheinfels@fau.de>
+###
+
+import numpy as np
+import scipy as sp
+import scipy.linalg
+import unittest
+
+import analysis
+from ellipsoid import Ellipsoid
+
+###
+### @brief Encapsulates the unit tests for @ref ellipsoid.Ellipsoid.__init__
+###
+class Test(unittest.TestCase):
+
+ ###
+ ### @brief Sets up the same random number generator for all tests
+ ###
+ def setUp(self):
+ self.rng = np.random.default_rng(0)
+
+ ###
+ ### @brief Checks that empty parameters produce a correctly initialized zero-dimensional ellipsoid
+ ###
+ def test_empty(self):
+ ellipsoid = Ellipsoid()
+ self.assertEqual(ellipsoid.n, 0)
+ self.assertFalse(ellipsoid.degenerate)
+ self.assertTrue(ellipsoid.centered)
+ self.assertTrue(np.allclose(ellipsoid.Pi, np.zeros((0, 0))))
+ self.assertTrue(np.allclose(ellipsoid.P, np.zeros((0,))))
+ self.assertTrue(np.allclose(ellipsoid.Pi_sqrt, np.zeros((0,))))
+ self.assertTrue(np.allclose(ellipsoid.c, np.zeros((0,))))
+
+ ###
+ ### @brief Draws 1000 sets of scalars and checks that the various combinations yield a correctly initialized ellipsoid
+ ###
+ def test_scalar(self):
+ for i in range(1000):
+ Pi = np.abs(self.rng.normal())
+ if self.rng.integers(0, 100) >= 90:
+ Pi = 0.0
+ c = self.rng.normal()
+ if self.rng.integers(0, 100) >= 90:
+ c = None
+ elif self.rng.integers(0, 50) >= 40:
+ c = 0.0
+ ellipsoid = Ellipsoid(Pi, c)
+ self.assertEqual(ellipsoid.n, 1)
+ self.assertEqual(ellipsoid.degenerate, np.isclose(Pi, 0.0))
+ if c is None:
+ self.assertTrue(ellipsoid.centered)
+ else:
+ self.assertEqual(ellipsoid.centered, np.allclose(c, 0.0))
+ self.assertEqual(float(ellipsoid.c), float(c))
+ self.assertEqual(float(ellipsoid.Pi), float(Pi))
+ self.assertAlmostEqual(float(ellipsoid.P), (Pi > 0.0) and 1.0/Pi or 0.0)
+ self.assertAlmostEqual(float(ellipsoid.Pi_sqrt), np.sqrt(float(Pi)))
+
+ ###
+ ### @brief Draws 1000 non-symmetric shape matrices of dimensions 2-20 and checks that an exception is raised
+ ###
+ def test_unsymmetric(self):
+ for i in range(1000):
+ n = self.rng.integers(2, 20)
+ Pi = self.rng.normal(size=(n, n))
+ Pi = Pi + Pi.T
+ Pi[0, n-1] = 1-Pi[0, n-1]
+ c = self.rng.normal(size=n)
+ if self.rng.integers(0, 100) >= 90:
+ c = None
+ elif self.rng.integers(0, 50) >= 40:
+ c = 0.0
+
+ with self.assertRaises(Exception):
+ Ellipsoid(Pi, c)
+
+ ###
+ ### @brief Draws 1000 spsd matrices of dimensions 2-20 with actual zero eigenvalues and checks that the ellipsoid is correctly initialized
+ ###
+ def test_spsd(self):
+ for i in range(1000):
+ n = self.rng.integers(1, 20)
+ Pi = self.rng.uniform(-1, +1, size=(n, n))
+ Pi = (Pi + Pi.T) / 2.0 + n * np.eye(n)
+ Pi[0, :] = 0.0
+ Pi[:, 0] = 0.0
+ c = self.rng.normal(size=n)
+ if self.rng.integers(0, 100) >= 90:
+ c = None
+ elif self.rng.integers(0, 50) >= 40:
+ c = np.zeros(n)
+
+ ellipsoid = Ellipsoid(Pi, c)
+ self.assertEqual(ellipsoid.n, n)
+ self.assertTrue(ellipsoid.degenerate)
+ if c is None:
+ self.assertTrue(ellipsoid.centered)
+ else:
+ self.assertEqual(ellipsoid.centered, np.allclose(c, 0.0))
+ self.assertTrue(np.allclose(ellipsoid.c, c))
+ self.assertTrue(np.allclose(ellipsoid.Pi, Pi))
+ self.assertTrue(np.allclose(ellipsoid.P, np.linalg.pinv(Pi)))
+ L, D, _ = sp.linalg.ldl(Pi)
+ self.assertTrue(np.allclose(ellipsoid.Pi_sqrt, L @ np.sqrt(D)))
+
+ ###
+ ### @brief Draws 1000 spd matrices of dimensions 2-20 with actual zero eigenvalues and checks that the ellipsoid is correctly initialized
+ ###
+ def test_spd(self):
+ for i in range(1000):
+ n = self.rng.integers(1, 20)
+ Pi = self.rng.uniform(-1, +1, size=(n, n))
+ Pi = (Pi + Pi.T) / 2.0 + n * np.eye(n)
+ c = self.rng.normal(size=n)
+ if self.rng.integers(0, 100) >= 90:
+ c = None
+ elif self.rng.integers(0, 50) >= 40:
+ c = np.zeros(n)
+
+ ellipsoid = Ellipsoid(Pi, c)
+ self.assertEqual(ellipsoid.n, n)
+ self.assertFalse(ellipsoid.degenerate)
+ if c is None:
+ self.assertTrue(ellipsoid.centered)
+ else:
+ self.assertEqual(ellipsoid.centered, np.allclose(c, 0.0))
+ self.assertTrue(np.allclose(ellipsoid.c, c))
+ self.assertTrue(np.allclose(ellipsoid.Pi, Pi))
+ self.assertTrue(np.allclose(ellipsoid.P, np.linalg.inv(Pi)))
+ self.assertTrue(np.allclose(ellipsoid.Pi_sqrt, np.linalg.cholesky(Pi)))
diff --git a/src/test/ellipsoid/contains.py b/src/test/ellipsoid/contains.py
new file mode 100644
index 0000000000000000000000000000000000000000..3df70177f832639585fa59fa70d21f2e63e0e2a2
--- /dev/null
+++ b/src/test/ellipsoid/contains.py
@@ -0,0 +1,71 @@
+## This file is part of the simulative evaluation for the qronos observer abstractions.
+## Copyright (C) 2022-2023 Tim Rheinfels <tim.rheinfels@fau.de>
+## See https://gitlab.cs.fau.de/qronos-state-abstractions/simulation
+##
+## Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met:
+##
+## 1. Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer.
+##
+## 2. Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following disclaimer in the documentation and/or other materials provided with the distribution.
+##
+## 3. Neither the name of the copyright holder nor the names of its contributors may be used to endorse or promote products derived from this software without specific prior written permission.
+##
+## THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS “AS IS” AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+
+###
+### @file test/ellipsoid/contains.py
+###
+### @brief Provides unit tests for @ref ellipsoid.Ellipsoid.contains
+###
+### @author Tim Rheinfels <tim.rheinfels@fau.de>
+###
+
+import numpy as np
+import scipy as sp
+import scipy.linalg
+import unittest
+
+import analysis
+from ellipsoid import Ellipsoid
+
+###
+### @brief Encapsulates the unit tests for @ref ellipsoid.Ellipsoid.contains
+###
+class Test(unittest.TestCase):
+
+ ###
+ ### @brief Sets up the same random number generator for all tests
+ ###
+ def setUp(self):
+ self.rng = np.random.default_rng(0)
+
+ ###
+ ### @brief Draws 1000 random ellipsoids of dimensions 1-20 and checks that their center is contained
+ ###
+ def test_random_center(self):
+ for i in range(1000):
+ n = self.rng.integers(1, 20)
+ Pi = np.eye(n)
+ Pi = (Pi + Pi.T) / 2.0 + n * np.eye(n)
+ c = self.rng.normal(size=n)
+
+ E = Ellipsoid(Pi, c)
+ self.assertTrue( E.contains(c) )
+
+ ###
+ ### @brief Draws 100 points from the surfaces of 1000 random ellipsoids of dimensions 1-20 each and checks that nudging the points in and our by a small factor yields the correct results for contains
+ ###
+ def test_random_nudge(self):
+ for i in range(1000):
+ n = self.rng.integers(1, 20)
+ Pi = np.eye(n)
+ Pi = (Pi + Pi.T) / 2.0 + n * np.eye(n)
+ c = self.rng.normal(size=n)
+
+ E = Ellipsoid(Pi, c)
+
+ for j in range(100):
+ x = E.draw_from_surface(self.rng)
+ f = self.rng.uniform(1e-3, 0.999)
+ self.assertTrue( E.contains((x-c)*f+c) )
+ self.assertFalse( E.contains((x-c)/f+c) )
diff --git a/src/test/ellipsoid/draw_from_surface.py b/src/test/ellipsoid/draw_from_surface.py
new file mode 100644
index 0000000000000000000000000000000000000000..dac99e53752e357ad73df66c5199cdec1824d449
--- /dev/null
+++ b/src/test/ellipsoid/draw_from_surface.py
@@ -0,0 +1,109 @@
+## This file is part of the simulative evaluation for the qronos observer abstractions.
+## Copyright (C) 2022-2023 Tim Rheinfels <tim.rheinfels@fau.de>
+## See https://gitlab.cs.fau.de/qronos-state-abstractions/simulation
+##
+## Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met:
+##
+## 1. Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer.
+##
+## 2. Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following disclaimer in the documentation and/or other materials provided with the distribution.
+##
+## 3. Neither the name of the copyright holder nor the names of its contributors may be used to endorse or promote products derived from this software without specific prior written permission.
+##
+## THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS “AS IS” AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+
+###
+### @file test/ellipsoid/draw_from_surface.py
+###
+### @brief Provides unit tests for @ref ellipsoid.Ellipsoid.draw_from_surface
+###
+### @author Tim Rheinfels <tim.rheinfels@fau.de>
+###
+
+import numpy as np
+import scipy as sp
+import scipy.linalg
+import scipy.stats
+import unittest
+
+import analysis
+from ellipsoid import Ellipsoid
+
+###
+### @brief Encapsulates the unit tests for @ref ellipsoid.Ellipsoid.draw_from_surface
+###
+class Test(unittest.TestCase):
+
+ ###
+ ### @brief Sets up the same random number generator for all tests
+ ###
+ def setUp(self):
+ self.rng = np.random.default_rng(0)
+
+ ###
+ ### @brief Draws 1000 points from surface of the empty ellipsoid and checks that the resulting vectors are empty themselves
+ ###
+ def test_empty(self):
+ Pi = np.zeros((0, 0))
+ E = Ellipsoid(Pi)
+
+ for i in range(1000):
+ self.assertEqual(E.draw_from_surface(self.rng).size, 0)
+
+ ###
+ ### @brief Draws 100 points from the surface of 1000 degenerate ellipsoids of dimension 2-10 each and checks their distance from the center transformed to the zero-space has negligible length
+ ###
+ def test_random_degenerate(self):
+ for i in range(1000):
+ n = self.rng.integers(2, 10)
+
+ # Draw random orthonormal basis for eigenvectors
+ V = sp.stats.ortho_group.rvs(n, random_state=self.rng)
+
+ # Draw random and positive eivenvalus
+ L = self.rng.uniform(1, 1e4, size=n)
+
+ # Draw directions of degeneracy by
+ # randomly setting the corresponding
+ # eigenvalues to zero
+ idx = []
+ for i in range(n):
+ if self.rng.integers(0, 100) >= 50:
+ idx.append(i)
+
+ idx = np.array(idx, dtype=int)
+
+ L[idx] = 1e-11
+
+ # Compose ellipsoid data
+ Pi = V.T @ np.diag(L) @ V
+ c = self.rng.normal(size=n)
+
+ E = Ellipsoid(Pi, c)
+
+ for j in range(100):
+ x = E.draw_from_surface(self.rng)
+
+ # Transform into zero-space
+ x_z = V[idx, :] @ (x - c)
+
+ # Check if zero
+ self.assertLess(np.linalg.norm(x_z), 1e-5)
+
+ ###
+ ### @brief Draws 100 points from the surface of 1000 non-degenerate ellipsoids of dimension 2-10 each and validates the result using @f$ (x-c)^T P (x-c) = 1 @f$
+ ###
+ def test_random_non_degenerate(self):
+ for i in range(1000):
+ n = self.rng.integers(1, 20)
+ Pi = np.eye(n)
+ Pi = (Pi + Pi.T) / 2.0 + n * np.eye(n)
+ c = self.rng.normal(size=n)
+
+ E = Ellipsoid(Pi, c)
+
+ P = np.linalg.inv(Pi)
+
+ for j in range(100):
+ x = E.draw_from_surface(self.rng)
+ self.assertAlmostEqual( np.dot((x-c), P@(x-c)), 1.0 )
diff --git a/src/test/ellipsoid/draw_from_volume.py b/src/test/ellipsoid/draw_from_volume.py
new file mode 100644
index 0000000000000000000000000000000000000000..ba20a348bbd5114a112ccd80811ad74892952b3e
--- /dev/null
+++ b/src/test/ellipsoid/draw_from_volume.py
@@ -0,0 +1,109 @@
+## This file is part of the simulative evaluation for the qronos observer abstractions.
+## Copyright (C) 2022-2023 Tim Rheinfels <tim.rheinfels@fau.de>
+## See https://gitlab.cs.fau.de/qronos-state-abstractions/simulation
+##
+## Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met:
+##
+## 1. Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer.
+##
+## 2. Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following disclaimer in the documentation and/or other materials provided with the distribution.
+##
+## 3. Neither the name of the copyright holder nor the names of its contributors may be used to endorse or promote products derived from this software without specific prior written permission.
+##
+## THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS “AS IS” AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+
+###
+### @file test/ellipsoid/draw_from_volume.py
+###
+### @brief Provides unit tests for @ref ellipsoid.Ellipsoid.draw_from_volume
+###
+### @author Tim Rheinfels <tim.rheinfels@fau.de>
+###
+
+import numpy as np
+import scipy as sp
+import scipy.linalg
+import scipy.stats
+import unittest
+
+import analysis
+from ellipsoid import Ellipsoid
+
+###
+### @brief Encapsulates the unit tests for @ref ellipsoid.Ellipsoid.draw_from_volume
+###
+class Test(unittest.TestCase):
+
+ ###
+ ### @brief Sets up the same random number generator for all tests
+ ###
+ def setUp(self):
+ self.rng = np.random.default_rng(0)
+
+ ###
+ ### @brief Draws 1000 points from volume of the empty ellipsoid and checks that the resulting vectors are empty themselves
+ ###
+ def test_empty(self):
+ Pi = np.zeros((0, 0))
+ E = Ellipsoid(Pi)
+
+ for i in range(1000):
+ self.assertEqual(E.draw_from_volume(self.rng).size, 0)
+
+ ###
+ ### @brief Draws 100 points from the volume of 1000 degenerate ellipsoids of dimension 2-10 each and checks their distance from the center transformed to the zero-space has negligible length
+ ###
+ def test_random_degenerate(self):
+ for i in range(1000):
+ n = self.rng.integers(2, 10)
+
+ # Draw random orthonormal basis for eigenvectors
+ V = sp.stats.ortho_group.rvs(n, random_state=self.rng)
+
+ # Draw random and positive eivenvalus
+ L = self.rng.uniform(1, 1e4, size=n)
+
+ # Draw directions of degeneracy by
+ # randomly setting the corresponding
+ # eigenvalues to zero
+ idx = []
+ for i in range(n):
+ if self.rng.integers(0, 100) >= 50:
+ idx.append(i)
+
+ idx = np.array(idx, dtype=int)
+
+ L[idx] = 1e-11
+
+ # Compose ellipsoid data
+ Pi = V.T @ np.diag(L) @ V
+ c = self.rng.normal(size=n)
+
+ E = Ellipsoid(Pi, c)
+
+ for j in range(100):
+ x = E.draw_from_volume(self.rng)
+
+ # Transform into zero-space
+ x_z = V[idx, :] @ (x - c)
+
+ # Check if zero
+ self.assertLess(np.linalg.norm(x_z), 1e-5)
+
+ ###
+ ### @brief Draws 100 points from the volume of 1000 non-degenerate ellipsoids of dimension 2-10 each and validates the result using @f$ (x-c)^T P (x-c) \leq 1 @f$
+ ###
+ def test_random(self):
+ for i in range(1000):
+ n = self.rng.integers(1, 20)
+ Pi = np.eye(n)
+ Pi = (Pi + Pi.T) / 2.0 + n * np.eye(n)
+ c = self.rng.normal(size=n)
+
+ E = Ellipsoid(Pi, c)
+
+ P = np.linalg.inv(Pi)
+
+ for j in range(100):
+ x = E.draw_from_volume(self.rng)
+ self.assertTrue( np.dot((x-c), P@(x-c)) <= 1.0 )
diff --git a/src/test/ellipsoid/norm.py b/src/test/ellipsoid/norm.py
new file mode 100644
index 0000000000000000000000000000000000000000..61b3b6dcf6d224316ae298e31df93a9392218556
--- /dev/null
+++ b/src/test/ellipsoid/norm.py
@@ -0,0 +1,74 @@
+## This file is part of the simulative evaluation for the qronos observer abstractions.
+## Copyright (C) 2022-2023 Tim Rheinfels <tim.rheinfels@fau.de>
+## See https://gitlab.cs.fau.de/qronos-state-abstractions/simulation
+##
+## Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met:
+##
+## 1. Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer.
+##
+## 2. Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following disclaimer in the documentation and/or other materials provided with the distribution.
+##
+## 3. Neither the name of the copyright holder nor the names of its contributors may be used to endorse or promote products derived from this software without specific prior written permission.
+##
+## THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS “AS IS” AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+
+###
+### @file test/ellipsoid/norm.py
+###
+### @brief Provides unit tests for @ref ellipsoid.Ellipsoid.norm
+###
+### @author Tim Rheinfels <tim.rheinfels@fau.de>
+###
+
+import numpy as np
+import scipy as sp
+import scipy.linalg
+import unittest
+
+import analysis
+from ellipsoid import Ellipsoid
+
+###
+### @brief Encapsulates the unit tests for @ref ellipsoid.Ellipsoid.norm
+###
+class Test(unittest.TestCase):
+
+ ###
+ ### @brief Sets up the same random number generator for all tests
+ ###
+ def setUp(self):
+ self.rng = np.random.default_rng(0)
+
+ ###
+ ### @brief Draws 100 normally distributed points for 1000 random spherical ellipsoids of dimension 1-20 each and checks that their ellipsoidal norm coincides with the analytical result
+ ###
+ def test_sphere(self):
+ for i in range(1000):
+ n = self.rng.integers(1, 20)
+ r = self.rng.uniform(1e-3, 1e+3)
+ Pi = r**2 * np.eye(n)
+ c = self.rng.normal(size=n)
+
+ E = Ellipsoid(Pi, c)
+
+ for j in range(100):
+ x = self.rng.normal(size=n)
+ self.assertAlmostEqual( np.linalg.norm(x-c) / r, E.norm(x) )
+
+ ###
+ ### @brief Draws 100 normally distributed points for 1000 non-degenerate ellipsoids of dimension 1-10 each, computes their ellipsoidal norm and validates the result using @f$ \norm{P^{\frac{T}{2}} (x-c)} @f$
+ ###
+ def test_random(self):
+ for i in range(1000):
+ n = self.rng.integers(1, 20)
+ Pi = np.eye(n)
+ Pi = (Pi + Pi.T) / 2.0 + n * np.eye(n)
+ c = self.rng.normal(size=n)
+
+ E = Ellipsoid(Pi, c)
+
+ P_sqrt_T = np.linalg.cholesky(np.linalg.inv(Pi)).T
+
+ for j in range(100):
+ x = self.rng.normal(size=n)
+ self.assertAlmostEqual( np.linalg.norm(P_sqrt_T@(x-c)), E.norm(x) )