From b3aa51a4efca327975bf64b458fc9de070310177 Mon Sep 17 00:00:00 2001
From: Tim Rheinfels <tim.rheinfels@fau.de>
Date: Wed, 24 May 2023 11:32:02 +0200
Subject: [PATCH] evaluation: add tools to create statistics for simulation
 runs

---
 src/evaluation/statistics.py                  | 205 ++++++++++++++++++
 src/test/evaluation/statistics/__init__.py    |   0
 .../statistics/compute_statistics.py          | 169 +++++++++++++++
 .../evaluation/statistics/count_or_zero.py    |  69 ++++++
 src/test/evaluation/statistics/max_or_zero.py |  60 +++++
 .../evaluation/statistics/mean_or_zero.py     |  60 +++++
 src/test/evaluation/statistics/std_or_zero.py |  60 +++++
 7 files changed, 623 insertions(+)
 create mode 100644 src/evaluation/statistics.py
 create mode 100644 src/test/evaluation/statistics/__init__.py
 create mode 100644 src/test/evaluation/statistics/compute_statistics.py
 create mode 100644 src/test/evaluation/statistics/count_or_zero.py
 create mode 100644 src/test/evaluation/statistics/max_or_zero.py
 create mode 100644 src/test/evaluation/statistics/mean_or_zero.py
 create mode 100644 src/test/evaluation/statistics/std_or_zero.py

diff --git a/src/evaluation/statistics.py b/src/evaluation/statistics.py
new file mode 100644
index 0000000..7f9f1a0
--- /dev/null
+++ b/src/evaluation/statistics.py
@@ -0,0 +1,205 @@
+## 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  evaluation/statistics.py
+###
+###  @brief  Provides the functionality for computing the evaluation statistics from a @ref closed_loop.ClosedLoopSimulation
+###
+###  @author  Tim Rheinfels  <tim.rheinfels@fau.de>
+###
+
+import logging
+import matplotlib.pyplot as plot
+import numpy as np
+
+from simulation.closed_loop import ClosedLoopSimulation
+
+###
+###  @brief  Computes the mean of @p x handling empty vectors correctly
+###
+###  @param  x  Vector of real values
+###
+###  @returns  Either
+###      - empirical mean of @p x if @p x is not empty or
+###      - 0 if @p x is empty
+###
+def mean_or_zero(x):
+    if len(x) == 0:
+        return 0.0
+    return np.mean(x)
+
+###
+###  @brief  Computes the standard deviation of @p x handling empty vectors correctly
+###
+###  @param  x  Vector of real values
+###
+###  @returns  Either
+###      - empirical standard deviation of @p x if @p x is not empty or
+###      - 0 if @p x is empty
+###
+def std_or_zero(x):
+    if len(x) == 0:
+        return 0.0
+    return np.std(x)
+
+###
+###  @brief  Computes the maximum of @p x handling empty vectors correctly
+###
+###  @param  x  Vector of real values
+###
+###  @returns  Either
+###      - empirical maximum of @p x if @p x is not empty or
+###      - 0 if @p x is empty
+###
+def max_or_zero(x):
+    if len(x) == 0:
+        return 0.0
+    return np.max(x)
+
+###
+###  @brief  Computes the number of non-zero entries of @p x handling empty vectors correctly
+###
+###  @param  x  Vector of real values
+###
+###  @returns  Either
+###      - number of non-zero entries of @p x if @p x is not empty or
+###      - 0 if @p x is empty
+###
+def count_or_zero(x):
+    if len(x) == 0:
+        return 0
+    return np.count_nonzero(x)
+
+###
+###  @brief  Computes the evaluation statistics resulting from and @p experiment
+###
+###  @param  experiment  A @ref closed_loop.ClosedLoopSimulation after it was run
+###
+###  @returns  Dictionary with five keys
+###
+###     The former are given by
+###       - 'norm_x_P': System's state under the analysis norm, i.e., @f$ \norm{x}_P @f$,
+###       - 'v_sys': System under the safety norm, i.e., @f$ v_{sys} := \norm{s}_S @f$,
+###       - 'v': Abstraction value @f$ v @f$, and
+###       - 'v_over_v_sys': Quotient quantifying the abstraction's overapproximation, i.e., @f$ \frac{v}{v_{sys}} @f$
+###
+###     where each is again a dictionary with the keys
+###       - mean: emprical mean,
+###       - std: empirical standard deviation,
+###       - max: empirical maximum,
+###       - violations: number of timesteps for which the quantity was greater than 1,
+###
+###     each of which contain a 3x3 array of values (dimension 0: switching no, yes, combined, dimension 1: disturbance no, yes, combined)
+###     splitting up the statistics in parts where switching / disturbances are absent or present.
+###
+###     The last key 'sigma' is again a dictionary with statistics about the switching signal:
+###       - 'counts': Absolute number of occurences for each symbol in the alphabet
+###       - 'ratios': Relative number of occurences for each symbol in the alphabet
+###
+###     where either key is a three-dimensional tensor. The first two dimensions are structured the same as the
+###     tableaus above (dimension 0: switching, dimension 1: disturbance) whereas the last dimension indicates
+###     the mode for which the statistics were generated. As an example, counts[1, 2, 0] denotes the absolute
+###     number of occurences for mode @f$ \sigma @f$ when switching is enabled whether a disturbance is present or not.
+###
+def compute_statistics(experiment):
+    assert(isinstance(experiment, ClosedLoopSimulation))
+
+    # Dim 0: Switching
+    # Dim 1: Disturbance
+    # Dim 2: Switching Alphabet
+    n_x = experiment.system_model_simulation.system_model.n_x
+    n_Sigma = experiment.system_model_simulation.system_model.n_Sigma
+    counts = np.zeros((3, 3, n_Sigma))
+    norm_x_P = np.array([[None] * 3] * 3)
+    v_sys = np.array([[None] * 3] * 3)
+    v = np.array([[None] * 3] * 3)
+    v_over_v_sys = np.array([[None] * 3] * 3)
+
+    timesteps = 0
+    for i in range(experiment.run_count):
+        run = experiment.get_run(i)
+
+        # Extract signals, match lengths, and convert to numpy arrays
+        t = np.array(run['sigma']['time_data'])
+        timesteps_run = t.size
+        timesteps += timesteps_run
+        sigma_run = np.concatenate(run['sigma']['signal_data'][:timesteps])
+        x_run = np.concatenate(run['x']['signal_data'][:timesteps])
+        v_sys_run = np.concatenate(run['v_sys']['signal_data'][:timesteps])
+        v_run = np.concatenate(run['v']['signal_data'][:timesteps])
+
+        # Compute x under p norm
+        norm_x_P_run = np.zeros(t.size)
+        for j in range(t.size):
+            norm_x_P_run[j] = experiment.state_abstraction_simulation.abstraction.E_P.norm(x_run[j])
+
+        # Create masks
+        masks_switching = (
+            (t < experiment.switching_sequence_simulation.T_start) | (t > experiment.switching_sequence_simulation.T_end),
+            (t >= experiment.switching_sequence_simulation.T_start) & (t <= experiment.switching_sequence_simulation.T_end),
+            np.array([True] * len(t)),
+        )
+
+        masks_disturbance = (
+            (t < experiment.disturbance_simulation.T_start) | (t > experiment.disturbance_simulation.T_end),
+            (t >= experiment.disturbance_simulation.T_start) & (t <= experiment.disturbance_simulation.T_end),
+            np.array([True] * len(t)),
+        )
+
+        # Iterate over masks
+        for j, mask_switching in enumerate(masks_switching):
+            for k, mask_disturbance in enumerate(masks_disturbance):
+
+                # Combine masks and get relevant indices
+                idx = np.where(mask_switching & mask_disturbance)[0]
+
+                # Iterate over alphabet to collect switching statistics
+                for sigma in range(n_Sigma):
+                    # Count the occurences of sigma == l given the masks
+                    counts[j, k, sigma] += np.count_nonzero(sigma_run[idx] == sigma)
+
+                # Gather relevant v_sys values
+                v_sys[j, k] = v_sys_run[idx]
+                v[j, k] = v_run[idx]
+                norm_x_P[j, k] = norm_x_P_run[idx]
+                v_over_v_sys[j, k] = v[j, k] / v_sys[j, k]
+
+    # Compute switching ratios
+    ratios = counts.astype(float) / float(timesteps)
+
+    # Arrange statistics nicely in a dictionary
+    statistics = {
+        'sigma': {
+            'counts': counts,
+            'ratios': ratios,
+        }
+    }
+
+    for signal_name, signal in zip(('norm_x_P', 'v_sys', 'v', 'v_over_v_sys'), (norm_x_P, v_sys, v, v_over_v_sys)):
+        statistics[signal_name] = {
+            'mean': np.zeros((3, 3)),
+            'std': np.zeros((3, 3)),
+            'max': np.zeros((3, 3)),
+            'violations': np.zeros((3, 3)),
+        }
+
+        for j in range(3):
+            for k in range(3):
+                statistics[signal_name]['mean'][j, k] = mean_or_zero(signal[j, k])
+                statistics[signal_name]['std'][j, k] = std_or_zero(signal[j, k])
+                statistics[signal_name]['max'][j, k] = max_or_zero(signal[j, k])
+                statistics[signal_name]['violations'][j, k] = count_or_zero(signal[j, k] > 1.0)
+
+    return statistics
diff --git a/src/test/evaluation/statistics/__init__.py b/src/test/evaluation/statistics/__init__.py
new file mode 100644
index 0000000..e69de29
diff --git a/src/test/evaluation/statistics/compute_statistics.py b/src/test/evaluation/statistics/compute_statistics.py
new file mode 100644
index 0000000..6e60e66
--- /dev/null
+++ b/src/test/evaluation/statistics/compute_statistics.py
@@ -0,0 +1,169 @@
+## 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/evaluation/statistics/compute_statistics.py
+###
+###  @brief  Provides unit tests for @ref statistics.compute_statistics
+###
+###  @author  Tim Rheinfels  <tim.rheinfels@fau.de>
+###
+
+import numpy as np
+import unittest
+
+from ellipsoid import Ellipsoid
+
+from state_abstraction import StateAbstraction
+from system_model import SystemModel
+
+from simulation.closed_loop import ClosedLoopSimulation
+from simulation.disturbance import DisturbanceSimulation
+from simulation.state_abstraction import StateAbstractionSimulation
+from simulation.switching_sequence import MultinomialSwitchingSequenceSimulation
+from simulation.system_model import SystemModelSimulation
+
+import evaluation.statistics as statistics
+
+runs = 5
+timesteps = 12
+
+###
+###  @brief  Encapsulates the unit tests for @ref statistics.compute_statistics
+###
+class Test(unittest.TestCase):
+
+    ###
+    ###  @brief  Sets up the same random number generator and test data for all tests
+    ###
+    def setUp(self):
+        self.rng = np.random.default_rng(0)
+
+        self.system_model = SystemModel(
+            'test_system',
+            [SystemModel.Mode(
+                'test_mode_%d' % i,
+                A,
+                1.0,
+                Ellipsoid(1.0),
+                1.0,
+                1.0,
+                Ellipsoid(1.0),
+            ) for i, A in enumerate((0.5, 1.0, 2.0))],
+            Ellipsoid(1.0),
+            None,
+            Ellipsoid(100.0),
+        )
+
+        E_P = Ellipsoid(1.0)
+
+        self.state_abstraction = StateAbstraction('test_abstraction', self.system_model, E_P)
+
+        self.disturbance_simulation = DisturbanceSimulation(
+            DisturbanceSimulation.Strategy.Zero,
+            self.system_model,
+            self.state_abstraction,
+            3.0,
+            9.0,
+        )
+
+        self.experiment = ClosedLoopSimulation.configure_with_mk_switching(
+            'test-simulation',
+            self.system_model,
+            self.state_abstraction,
+            self.disturbance_simulation,
+            [1, 0],
+            1,
+            1,
+            9.0,
+            12.0,
+        )
+
+    ###
+    ###  @brief  Checks that an empty vector yields a zero count
+    ###
+    def test_invalid(self):
+        for experiment in (None, 0.1, 0, -1, ClosedLoopSimulation.simulate):
+            with self.assertRaises(Exception):
+                statistics.compute_statistics( experiment )
+
+    ###
+    ###  @brief  Checks that the output format is correcty when inserting valid simulation data
+    ###
+    def test_output_format(self):
+        self.experiment.simulate(runs, timesteps)
+
+        stats = statistics.compute_statistics(self.experiment)
+
+        self.assertTrue( isinstance(stats, dict) )
+        self.assertEqual( set(stats.keys()), set(('sigma', 'norm_x_P', 'v_sys', 'v', 'v_over_v_sys')) )
+
+        self.assertTrue( isinstance(stats['sigma'], dict) )
+        self.assertEqual( set(stats['sigma'].keys()), set(('counts', 'ratios')) )
+        for stat in ('counts', 'ratios'):
+            self.assertTrue( isinstance( stats['sigma'][stat], np.ndarray ) )
+            self.assertEqual( stats['sigma'][stat].ndim, 3 )
+            self.assertEqual( stats['sigma'][stat].shape[0], 3 )
+            self.assertEqual( stats['sigma'][stat].shape[1], 3 )
+            self.assertEqual( stats['sigma'][stat].shape[2], self.system_model.n_Sigma )
+
+        for signal in ('norm_x_P', 'v_sys', 'v', 'v_over_v_sys'):
+            self.assertTrue( isinstance(stats[signal], dict) )
+            self.assertEqual( set(stats[signal].keys()), set(('mean', 'std', 'max', 'violations')) )
+
+            for stat in ('mean', 'std', 'max', 'violations'):
+                self.assertTrue( isinstance( stats[signal][stat], np.ndarray ) )
+                self.assertEqual( stats[signal][stat].ndim, 2 )
+                self.assertEqual( stats[signal][stat].shape[0], 3 )
+                self.assertEqual( stats[signal][stat].shape[1], 3 )
+
+    ###
+    ###  @brief  Checks that the sums of switching counts match across the tableau
+    ###
+    def test_switching_counts_sum(self):
+        self.experiment.simulate(runs, timesteps)
+
+        stats = statistics.compute_statistics(self.experiment)
+
+        summed = np.sum(stats['sigma']['counts'], axis=2)
+
+        self.assertEqual( summed[0, 0] + summed[0, 1], summed[0, 2] )
+        self.assertEqual( summed[1, 0] + summed[1, 1], summed[1, 2] )
+        self.assertEqual( summed[2, 0] + summed[2, 1], summed[2, 2] )
+
+        self.assertEqual( summed[0, 0] + summed[1, 0], summed[2, 0] )
+        self.assertEqual( summed[0, 1] + summed[1, 1], summed[2, 1] )
+        self.assertEqual( summed[0, 2] + summed[1, 2], summed[2, 2] )
+
+        self.assertEqual( summed[2, 2], runs * (timesteps+1) )
+
+    ###
+    ###  @brief  Checks that the sums of switching ratios match across the tableau
+    ###
+    def test_switching_ratios_sum(self):
+        self.experiment.simulate(runs, timesteps)
+
+        stats = statistics.compute_statistics(self.experiment)
+
+        summed = np.sum(stats['sigma']['ratios'], axis=2)
+
+        self.assertAlmostEqual( summed[0, 0] + summed[0, 1], summed[0, 2] )
+        self.assertAlmostEqual( summed[1, 0] + summed[1, 1], summed[1, 2] )
+        self.assertAlmostEqual( summed[2, 0] + summed[2, 1], summed[2, 2] )
+
+        self.assertAlmostEqual( summed[0, 0] + summed[1, 0], summed[2, 0] )
+        self.assertAlmostEqual( summed[0, 1] + summed[1, 1], summed[2, 1] )
+        self.assertAlmostEqual( summed[0, 2] + summed[1, 2], summed[2, 2] )
+
+        self.assertAlmostEqual( summed[2, 2], 1.0 )
diff --git a/src/test/evaluation/statistics/count_or_zero.py b/src/test/evaluation/statistics/count_or_zero.py
new file mode 100644
index 0000000..bf6c755
--- /dev/null
+++ b/src/test/evaluation/statistics/count_or_zero.py
@@ -0,0 +1,69 @@
+## 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/evaluation/statistics/count_or_zero.py
+###
+###  @brief  Provides unit tests for @ref statistics.count_or_zero
+###
+###  @author  Tim Rheinfels  <tim.rheinfels@fau.de>
+###
+
+import numpy as np
+import unittest
+
+import evaluation.statistics as statistics
+
+###
+###  @brief  Encapsulates the unit tests for @ref statistics.count_or_zero
+###
+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 an empty vector yields a zero count
+    ###
+    def test_empty(self):
+        for vector in ([], np.zeros(0)):
+            self.assertEqual( statistics.count_or_zero(vector), 0.0 )
+
+    ###
+    ###  @brief  Checks that the count of a single 0 is 0
+    ###
+    def test_zero(self):
+        self.assertEqual( statistics.count_or_zero([0]), 0 )
+
+    ###
+    ###  @brief  Draws 1000 nonzero random scalars and checks that they equal their count
+    ###
+    def test_nonzero_scalar(self):
+        for i in range(1000):
+            while True:
+                x = self.rng.normal()
+                if np.abs(x) >= 1e-4: #pragma: no cover
+                    break
+            self.assertEqual( statistics.count_or_zero([x]), 1 )
+
+    ###
+    ###  @brief  Draws 100 random scalars 1000 times and checks that the count is computed correctly
+    ###
+    def test_random(self):
+        for i in range(1000):
+            x = self.rng.normal(size=(100,)) * self.rng.binomial(1, 0.5, size=(100,))
+            self.assertEqual( statistics.count_or_zero(x), np.count_nonzero(x) )
diff --git a/src/test/evaluation/statistics/max_or_zero.py b/src/test/evaluation/statistics/max_or_zero.py
new file mode 100644
index 0000000..c09d87b
--- /dev/null
+++ b/src/test/evaluation/statistics/max_or_zero.py
@@ -0,0 +1,60 @@
+## 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/evaluation/statistics/max_or_zero.py
+###
+###  @brief  Provides unit tests for @ref statistics.max_or_zero
+###
+###  @author  Tim Rheinfels  <tim.rheinfels@fau.de>
+###
+
+import numpy as np
+import unittest
+
+import evaluation.statistics as statistics
+
+###
+###  @brief  Encapsulates the unit tests for @ref statistics.max_or_zero
+###
+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 an empty vector yields a zero max
+    ###
+    def test_empty(self):
+        for vector in ([], np.zeros(0)):
+            self.assertEqual( statistics.max_or_zero(vector), 0.0 )
+
+    ###
+    ###  @brief  Draws 1000 random scalars and checks that they equal their max
+    ###
+    def test_scalar(self):
+        for i in range(1000):
+            x = self.rng.normal()
+            self.assertEqual( statistics.max_or_zero([x]), x )
+
+    ###
+    ###  @brief  Draws 100 random scalars 1000 times and checks that the max is computed correctly
+    ###
+    def test_random(self):
+        for i in range(1000):
+            x = self.rng.normal(size=(100,))
+            self.assertEqual( statistics.max_or_zero(x), np.max(x) )
diff --git a/src/test/evaluation/statistics/mean_or_zero.py b/src/test/evaluation/statistics/mean_or_zero.py
new file mode 100644
index 0000000..53a404f
--- /dev/null
+++ b/src/test/evaluation/statistics/mean_or_zero.py
@@ -0,0 +1,60 @@
+## 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/evaluation/statistics/mean_or_zero.py
+###
+###  @brief  Provides unit tests for @ref statistics.mean_or_zero
+###
+###  @author  Tim Rheinfels  <tim.rheinfels@fau.de>
+###
+
+import numpy as np
+import unittest
+
+import evaluation.statistics as statistics
+
+###
+###  @brief  Encapsulates the unit tests for @ref statistics.mean_or_zero
+###
+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 an empty vector yields a zero mean
+    ###
+    def test_empty(self):
+        for vector in ([], np.zeros(0)):
+            self.assertEqual( statistics.mean_or_zero(vector), 0.0 )
+
+    ###
+    ###  @brief  Draws 1000 random scalars and checks that they equal their mean
+    ###
+    def test_scalar(self):
+        for i in range(1000):
+            x = self.rng.normal()
+            self.assertEqual( statistics.mean_or_zero([x]), x )
+
+    ###
+    ###  @brief  Draws 100 random scalars 1000 times and checks that the mean is computed correctly
+    ###
+    def test_random(self):
+        for i in range(1000):
+            x = self.rng.normal(size=(100,))
+            self.assertEqual( statistics.mean_or_zero(x), np.mean(x) )
diff --git a/src/test/evaluation/statistics/std_or_zero.py b/src/test/evaluation/statistics/std_or_zero.py
new file mode 100644
index 0000000..11f163f
--- /dev/null
+++ b/src/test/evaluation/statistics/std_or_zero.py
@@ -0,0 +1,60 @@
+## 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/evaluation/statistics/std_or_zero.py
+###
+###  @brief  Provides unit tests for @ref statistics.std_or_zero
+###
+###  @author  Tim Rheinfels  <tim.rheinfels@fau.de>
+###
+
+import numpy as np
+import unittest
+
+import evaluation.statistics as statistics
+
+###
+###  @brief  Encapsulates the unit tests for @ref statistics.std_or_zero
+###
+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 an empty vector yields a zero std
+    ###
+    def test_empty(self):
+        for vector in ([], np.zeros(0)):
+            self.assertEqual( statistics.std_or_zero(vector), 0.0 )
+
+    ###
+    ###  @brief  Draws 1000 random scalars and checks that their standard deviation equals zero
+    ###
+    def test_scalar(self):
+        for i in range(1000):
+            x = self.rng.normal()
+            self.assertEqual( statistics.std_or_zero([x]), 0.0 )
+
+    ###
+    ###  @brief  Draws 100 random scalars 1000 times and checks that the std is computed correctly
+    ###
+    def test_random(self):
+        for i in range(1000):
+            x = self.rng.normal(size=(100,))
+            self.assertEqual( statistics.std_or_zero(x), np.std(x) )
-- 
GitLab