update benchmarks and doc

JianqiangDing · JianqiangDing · commit c25a6afbbce8 · 2024-01-18T09:12:18.000+08:00
diff --git a/README.md b/README.md
@@ -117,45 +117,41 @@ $$
 
 ```python
 import numpy as np
-
-from pybdr.algorithm import ASB2008CDC
-from pybdr.dynamic_system import NonLinSys
+from pybdr.algorithm import ASB2008CDC, ASB2008CDC
+from pybdr.util.functional import performance_counter, performance_counter_start
 from pybdr.geometry import Zonotope, Interval, Geometry
-from pybdr.geometry.operation import boundary
+from pybdr.geometry.operation import boundary, cvt2
 from pybdr.model import *
 from pybdr.util.visualization import plot
 
-# init dynamic system
-system = NonLinSys(Model(vanderpol, [2, 1]))
-
-# settings for the computation
-options = ASB2008CDC.Options()
-options.t_end = 6.74
-options.step = 0.005
-options.tensor_order = 3
-options.taylor_terms = 4
-options.u = Zonotope.zero(1, 1)
-options.u_trans = np.zeros(1)
-
-z = Zonotope([1.4, 2.4], np.diag([0.17, 0.06]))
-
-# -----------------------------------------------------
-# computing reachable sets with boundary analysis
-options.r0 = boundary(z, 1, Geometry.TYPE.ZONOTOPE)
-# -----------------------------------------------------
-# computing reachable sets without boundary analysis
-# options.r0 = [z]
-# -----------------------------------------------------
-
-# settings for the using geometry
-Zonotope.REDUCE_METHOD = Zonotope.REDUCE_METHOD.GIRARD
-Zonotope.ORDER = 50
-
-# reachable sets computation
-ti, tp, _, _ = ASB2008CDC.reach(system, options)
-
-# visualize the results
-plot(tp, [0, 1])
+if __name__ == '__main__':
+    # settings for the computation
+    options = ASB2008CDC.Options()
+    options.t_end = 6.74
+    options.step = 0.005
+    options.tensor_order = 3
+    options.taylor_terms = 4
+    options.u = Zonotope.zero(1, 1)
+    options.u_trans = np.zeros(1)
+
+    # settings for the using geometry
+    Zonotope.REDUCE_METHOD = Zonotope.REDUCE_METHOD.GIRARD
+    Zonotope.ORDER = 50
+
+    z = Interval([1.23, 2.34], [1.57, 2.46])
+    x0 = cvt2(z, Geometry.TYPE.ZONOTOPE)
+    xs = boundary(z, 1, Geometry.TYPE.ZONOTOPE)
+
+    this_time = performance_counter_start()
+    ri_without_bound, rp_without_bound = ASB2008CDC.reach(vanderpol, [2, 1], options, x0)
+    this_time = performance_counter(this_time, 'reach_without_bound')
+
+    ri_with_bound, rp_with_bound = ASB2008CDC.reach_parallel(vanderpol, [2, 1], options, xs)
+    this_time = performance_counter(this_time, 'reach_with_bound')
+
+    # visualize the results
+    plot(ri_without_bound, [0, 1])
+    plot(ri_with_bound, [0, 1])
 ```
 
 |     With Boundary Analysis (BA)     |       No Boundary Analysis (NBA)       |
@@ -164,16 +160,16 @@ plot(tp, [0, 1])
 
 For large initial sets,
 
-|                                                         System                                                          |                                                                 Code                                                                 |   Reachable Sets (Orange-NBA,Blue-BA)   |
-|:-----------------------------------------------------------------------------------------------------------------------:|:------------------------------------------------------------------------------------------------------------------------------------:|:---------------------------------------:|
-|        [synchronous machine](https://github.com/ASAG-ISCAS/PyBDR/blob/master/pyrat/model/synchronous_machine.py)        | [benchmark_synchronous_machine_cmp.py](https://github.com/ASAG-ISCAS/PyBDR/blob/master/example/benchmark_synchronous_machine_cmp.py) |   ![](doc/imgs/sync_machine_cmp.png)    |
-| [Lotka Volterra model of 2 variables](https://github.com/ASAG-ISCAS/PyBDR/blob/master/pyrat/model/lotka_volterra_2d.py) |   [benchmark_lotka_volterra_2d_cmp.py](https://github.com/ASAG-ISCAS/PyBDR/blob/master/example/benchmark_lotka_volterra_2d_cmp.py)   | ![](doc/imgs/lotka_volterra_2d_cmp.png) |
-|                 [Jet engine](https://github.com/ASAG-ISCAS/PyBDR/blob/master/pyrat/model/jet_engine.py)                 |          [benchmark_jet_engine_cmp.py](https://github.com/ASAG-ISCAS/PyBDR/blob/master/example/benchmark_jet_engine_cmp.py)          |    ![](doc/imgs/jet_engine_cmp.png)     |
+|                                                         System                                                          |                                                                  Code                                                                   |   Reachable Sets (Orange-NBA,Blue-BA)   |
+|:-----------------------------------------------------------------------------------------------------------------------:|:---------------------------------------------------------------------------------------------------------------------------------------:|:---------------------------------------:|
+|        [synchronous machine](https://github.com/ASAG-ISCAS/PyBDR/blob/master/pyrat/model/synchronous_machine.py)        | [benchmark_synchronous_machine_cmp.py](https://github.com/ASAG-ISCAS/PyBDR/blob/master/benchmarks/benchmark_synchronous_machine_cmp.py) |   ![](doc/imgs/sync_machine_cmp.png)    |
+| [Lotka Volterra model of 2 variables](https://github.com/ASAG-ISCAS/PyBDR/blob/master/pyrat/model/lotka_volterra_2d.py) |   [benchmark_lotka_volterra_2d_cmp.py](https://github.com/ASAG-ISCAS/PyBDR/blob/master/benchmarks/benchmark_lotka_volterra_2d_cmp.py)   | ![](doc/imgs/lotka_volterra_2d_cmp.png) |
+|                 [Jet engine](https://github.com/ASAG-ISCAS/PyBDR/blob/master/pyrat/model/jet_engine.py)                 |          [benchmark_jet_engine_cmp.py](https://github.com/ASAG-ISCAS/PyBDR/blob/master/benchmarks/benchmark_jet_engine_cmp.py)          |    ![](doc/imgs/jet_engine_cmp.png)     |
 
 For large time horizons, i.e. consider
 the system [Brusselator](https://github.com/ASAG-ISCAS/PyBDR/blob/master/pyrat/model/brusselator.py)
 > For more details about the following example, please refer to
-> our [code](https://github.com/ASAG-ISCAS/PyBDR/blob/master/example/benchmark_brusselator_cmp.py).
+> our [code](https://github.com/ASAG-ISCAS/PyBDR/blob/master/benchmarks/benchmark_brusselator_cmp.py).
 
 | Time instance | With Boundary Analysis                |        Without Boundary Analysi        |
 |:-------------:|---------------------------------------|:--------------------------------------:|
@@ -238,47 +234,43 @@ $$
 ```python
 import numpy as np
 
-np.seterr(divide='ignore', invalid='ignore')
-import sys
-
-# sys.path.append("./../../") uncomment this line if you need to add path manually
 from pybdr.algorithm import ASB2008CDC
-from pybdr.dynamic_system import NonLinSys
 from pybdr.geometry import Zonotope, Interval, Geometry
 from pybdr.model import *
-from pybdr.util.visualization import plot
-from pybdr.util.functional.neural_ode_generate import neuralODE
-from pybdr.geometry.operation import boundary
+from pybdr.util.visualization import plot, plot_cmp
+from pybdr.geometry.operation import boundary, cvt2
+from pybdr.util.functional import performance_counter_start, performance_counter
 
-# init neural ODE
-system = NonLinSys(Model(neuralODE, [2, 1]))
+if __name__ == '__main__':
+    # settings for the computation
+    options = ASB2008CDC.Options()
+    options.t_end = 1
+    options.step = 0.01
+    options.tensor_order = 2
+    options.taylor_terms = 2
 
-# settings for the computation
-options = ASB2008CDC.Options()
-options.t_end = 1.5
-options.step = 0.01
-options.tensor_order = 2
-options.taylor_terms = 2
-Z = Zonotope([0.5, 0], np.diag([1, 0.5]))
+    options.u = Zonotope([0], np.diag([0]))
+    options.u_trans = options.u.c
 
-# Reachable sets computed with boundary analysis
-# options.r0 = boundary(Z,1,Geometry.TYPE.ZONOTOPE)
+    # settings for the using geometry
+    Zonotope.REDUCE_METHOD = Zonotope.REDUCE_METHOD.GIRARD
+    Zonotope.ORDER = 50
 
-# Reachable sets computed without boundary analysis
-options.r0 = [Z]
+    z = Interval([0, -0.5], [1, 0.5])
+    x0 = cvt2(z, Geometry.TYPE.ZONOTOPE)
+    xs = boundary(z, 2, Geometry.TYPE.ZONOTOPE)
 
-options.u = Zonotope.zero(1, 1)
-options.u_trans = np.zeros(1)
+    print(len(xs))
 
-# settings for the using geometry
-Zonotope.REDUCE_METHOD = Zonotope.REDUCE_METHOD.GIRARD
-Zonotope.ORDER = 50
+    this_time = performance_counter_start()
+    ri_without_bound, rp_without_bound = ASB2008CDC.reach(neural_ode_spiral1, [2, 1], options, x0)
+    this_time = performance_counter(this_time, "reach_without_bound")
 
-# reachable sets computation
-ti, tp, _, _ = ASB2008CDC.reach(system, options)
+    ri_with_bound, rp_with_bound = ASB2008CDC.reach_parallel(neural_ode_spiral1, [2, 1], options, xs)
+    this_time = performance_counter(this_time, "reach_with_bound")
 
-# visualize the results
-plot(tp, [0, 1])
+    # visualize the results
+    plot_cmp([ri_without_bound, ri_with_bound], [0, 1], cs=["#FF5722", "#303F9F"])
 ```
 
 In the following table, we show the reachable computed with boundary analysis and without boundary analysis on different
diff --git a/benchmarks/benchmark_boundary_analysis_vs_naive_partition_cmp.py b/benchmarks/benchmark_boundary_analysis_vs_naive_partition_cmp.py
@@ -2,9 +2,8 @@
 from pybdr.geometry import Interval, Zonotope, Geometry
 from pybdr.geometry.operation import cvt2, partition, boundary
 from pybdr.algorithm import ASB2008CDC
-from pybdr.dynamic_system import NonLinSys
 from pybdr.model import *
-from pybdr.util.visualization import plot, plot_cmp
+from pybdr.util.visualization import plot_cmp
 from pybdr.util.functional import performance_counter, performance_counter_start
 
 if __name__ == '__main__':
@@ -24,35 +23,52 @@
 
     init_set = Interval.identity(2) * 0.5 + 3
 
-    time_start = performance_counter_start()
+    this_time = performance_counter_start()
     _, rp_without_bound = ASB2008CDC.reach(lotka_volterra_2d, [2, 1], options, cvt2(init_set, Geometry.TYPE.ZONOTOPE))
+    this_time = performance_counter(this_time, 'reach_without_bound')
 
     # NAIVE PARTITION
     # --------------------------------------------------------
     # options.r0 = partition(z, 1, Geometry.TYPE.ZONOTOPE)
-    # xs = partition(init_set, 1, Geometry.TYPE.ZONOTOPE)
+    xs_naive_partition_4 = partition(init_set, 1, Geometry.TYPE.ZONOTOPE)
     # 4
     # ASB2008CDCParallel.reach_parallel cost: 20.126129667s MacOS
     # --------------------------------------------------------
-    # xs = partition(init_set, 0.5, Geometry.TYPE.ZONOTOPE)
+    xs_naive_partition_9 = partition(init_set, 0.5, Geometry.TYPE.ZONOTOPE)
     # 9
     # ASB2008CDCParallel.reach_parallel cost: 23.938516459000002s MacOS
     # --------------------------------------------------------
-    # xs = partition(init_set, 0.2, Geometry.TYPE.ZONOTOPE)
+    xs_naive_partition_36 = partition(init_set, 0.2, Geometry.TYPE.ZONOTOPE)
     # 36
     # ASB2008CDCParallel.reach_parallel cost: 65.447113125s MacOS
     # --------------------------------------------------------
 
     #  BOUNDAYR ANALYSIS
     # --------------------------------------------------------
-    xs = boundary(init_set, 1, Geometry.TYPE.ZONOTOPE)
+    xs_with_bound_8 = boundary(init_set, 1, Geometry.TYPE.ZONOTOPE)
     # 8
     # ASB2008CDCParallel.reach_parallel cost: 22.185758250000003s MacOS
 
-    print(len(xs))
+    print(len(xs_naive_partition_4))
+    print(len(xs_naive_partition_9))
+    print(len(xs_naive_partition_36))
+    print(len(xs_with_bound_8))
 
-    _, rp_with_bound = ASB2008CDC.reach_parallel(lotka_volterra_2d, [2, 1], options, xs)
+    this_time = performance_counter(this_time, 'partition and boundary')
 
-    performance_counter(time_start, "ASB2008CDCParallel.reach_parallel")
+    _, rp_naive_partition_4 = ASB2008CDC.reach_parallel(lotka_volterra_2d, [2, 1], options, xs_naive_partition_4)
+    this_time = performance_counter(this_time, 'reach naive partition 4')
 
-    plot_cmp([rp_without_bound, rp_with_bound], [0, 1], cs=["#FF5722", "#303F9F"])
+    _, rp_naive_partition_9 = ASB2008CDC.reach_parallel(lotka_volterra_2d, [2, 1], options, xs_naive_partition_9)
+    this_time = performance_counter(this_time, 'reach naive partition 9')
+
+    _, rp_naive_partition_36 = ASB2008CDC.reach_parallel(lotka_volterra_2d, [2, 1], options, xs_naive_partition_36)
+    this_time = performance_counter(this_time, 'reach naive partition 36')
+
+    _, rp_with_bound_8 = ASB2008CDC.reach_parallel(lotka_volterra_2d, [2, 1], options, xs_with_bound_8)
+    this_time = performance_counter(this_time, 'reach with bound 8')
+
+    plot_cmp([rp_without_bound, rp_naive_partition_4], [0, 1], cs=["#FF5722", "#303F9F"])
+    plot_cmp([rp_without_bound, rp_naive_partition_9], [0, 1], cs=["#FF5722", "#303F9F"])
+    plot_cmp([rp_without_bound, rp_naive_partition_36], [0, 1], cs=["#FF5722", "#303F9F"])
+    plot_cmp([rp_without_bound, rp_with_bound_8], [0, 1], cs=["#FF5722", "#303F9F"])
diff --git a/benchmarks/benchmark_node_spiral1_cmp.py b/benchmarks/benchmark_node_spiral1_cmp.py
@@ -0,0 +1,39 @@
+import numpy as np
+
+from pybdr.algorithm import ASB2008CDC
+from pybdr.geometry import Zonotope, Interval, Geometry
+from pybdr.model import *
+from pybdr.util.visualization import plot, plot_cmp
+from pybdr.geometry.operation import boundary, cvt2
+from pybdr.util.functional import performance_counter_start, performance_counter
+
+if __name__ == '__main__':
+    # settings for the computation
+    options = ASB2008CDC.Options()
+    options.t_end = 1
+    options.step = 0.01
+    options.tensor_order = 2
+    options.taylor_terms = 2
+
+    options.u = Zonotope([0], np.diag([0]))
+    options.u_trans = options.u.c
+
+    # settings for the using geometry
+    Zonotope.REDUCE_METHOD = Zonotope.REDUCE_METHOD.GIRARD
+    Zonotope.ORDER = 50
+
+    z = Interval([0, -0.5], [1, 0.5])
+    x0 = cvt2(z, Geometry.TYPE.ZONOTOPE)
+    xs = boundary(z, 2, Geometry.TYPE.ZONOTOPE)
+
+    print(len(xs))
+
+    this_time = performance_counter_start()
+    ri_without_bound, rp_without_bound = ASB2008CDC.reach(neural_ode_spiral1, [2, 1], options, x0)
+    this_time = performance_counter(this_time, "reach_without_bound")
+
+    ri_with_bound, rp_with_bound = ASB2008CDC.reach_parallel(neural_ode_spiral1, [2, 1], options, xs)
+    this_time = performance_counter(this_time, "reach_with_bound")
+
+    # visualize the results
+    plot_cmp([ri_without_bound, ri_with_bound], [0, 1], cs=["#FF5722", "#303F9F"])
diff --git a/benchmarks/benchmark_vanderpol_cmp.py b/benchmarks/benchmark_vanderpol_cmp.py
@@ -0,0 +1,36 @@
+import numpy as np
+from pybdr.algorithm import ASB2008CDC, ASB2008CDC
+from pybdr.util.functional import performance_counter, performance_counter_start
+from pybdr.geometry import Zonotope, Interval, Geometry
+from pybdr.geometry.operation import boundary, cvt2
+from pybdr.model import *
+from pybdr.util.visualization import plot
+
+if __name__ == '__main__':
+    # settings for the computation
+    options = ASB2008CDC.Options()
+    options.t_end = 6.74
+    options.step = 0.005
+    options.tensor_order = 3
+    options.taylor_terms = 4
+    options.u = Zonotope.zero(1, 1)
+    options.u_trans = np.zeros(1)
+
+    # settings for the using geometry
+    Zonotope.REDUCE_METHOD = Zonotope.REDUCE_METHOD.GIRARD
+    Zonotope.ORDER = 50
+
+    z = Interval([1.23, 2.34], [1.57, 2.46])
+    x0 = cvt2(z, Geometry.TYPE.ZONOTOPE)
+    xs = boundary(z, 1, Geometry.TYPE.ZONOTOPE)
+
+    this_time = performance_counter_start()
+    ri_without_bound, rp_without_bound = ASB2008CDC.reach(vanderpol, [2, 1], options, x0)
+    this_time = performance_counter(this_time, 'reach_without_bound')
+
+    ri_with_bound, rp_with_bound = ASB2008CDC.reach_parallel(vanderpol, [2, 1], options, xs)
+    this_time = performance_counter(this_time, 'reach_with_bound')
+
+    # visualize the results
+    plot(ri_without_bound, [0, 1])
+    plot(ri_with_bound, [0, 1])
diff --git a/pybdr/algorithm/asb2008cdc.py b/pybdr/algorithm/asb2008cdc.py
@@ -53,7 +53,7 @@ def linearize(dyn: Callable, dims, r: Geometry.Base, opt: Options):
         a = sys.evaluate((opt.lin_err_x, opt.lin_err_u), "numpy", 1, 0)
         b = sys.evaluate((opt.lin_err_x, opt.lin_err_u), "numpy", 1, 1)
         assert not (np.any(np.isnan(a))) or np.any(np.isnan(b))
-        lin_sys = LinSys(xa=a)
+        lin_sys = LinSys(xa=a, ub=None)
         lin_opt = ALK2011HSCC.Options()
         lin_opt.step = opt.step
         lin_opt.taylor_terms = opt.taylor_terms
diff --git a/pybdr/dynamic_system/continuous_system/linear_system.py b/pybdr/dynamic_system/continuous_system/linear_system.py
@@ -6,18 +6,21 @@
 
 
 class LinSys:
-    def __init__(self, xa, ub):
+    def __init__(self, xa, ub=None):
         self._xa = np.atleast_2d(xa)
-        self._ub = np.atleast_2d(ub)
-        assert self._xa.shape[1] == self._ub.shape[0]
+        if ub is not None:
+            self._ub = np.atleast_2d(ub)
+            assert self._xa.shape[1] == self._ub.shape[0]
+        else:
+            self._ub = None
 
     @property
     def dim(self):
         return self._xa.shape[1]
 
     @property
     def type(self):
-        return 'linear_simple'
+        return 'linear'
 
     @property
     def xa(self):
diff --git a/pybdr/dynamic_system/continuous_system/nonlinear_system.py b/pybdr/dynamic_system/continuous_system/nonlinear_system.py
@@ -1,5 +1,4 @@
 from __future__ import annotations
-import numpy as np
 from pybdr.model import Model
 
 
diff --git a/pybdr/util/visualization/plot.py b/pybdr/util/visualization/plot.py
@@ -95,7 +95,7 @@ def __2d_plot(
                 raise NotImplementedError
 
     ax.autoscale_view()
-    ax.axis("equal")
+    # ax.axis("equal")
     ax.set_xlabel("x" + str(dims[0]))
     ax.set_ylabel("x" + str(dims[1]))
 

Original file line number	Diff line number	Diff line change
`@@ -1,5 +1,4 @@`
`1`	`1`	`from __future__ import annotations`
`2`		`-import numpy as np`
`3`	`2`	`from pybdr.model import Model`
`4`	`3`
`5`	`4`