Skip to content

Commit 980a1f4

Browse files
TaoranWuTaoranWu
committed
add codac implementation and benchmark
1 parent 61b081e commit 980a1f4

File tree

6 files changed

+305
-7
lines changed

6 files changed

+305
-7
lines changed

‎benchmarks/benchmark_brusselator_nonconvex.py‎

Lines changed: 45 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,45 @@
1+
import numpy as np
2+
import codac
3+
from pybdr.algorithm import ASB2008CDC
4+
from pybdr.geometry import Zonotope, Interval, Polytope, Geometry
5+
from pybdr.geometry.operation import cvt2
6+
from pybdr.model import *
7+
from pybdr.util.visualization import plot_cmp, plot
8+
from pybdr.util.functional import extract_boundary
9+
from pybdr.util.functional import performance_counter, performance_counter_start
10+
11+
if __name__ == '__main__':
12+
# settings for the computation
13+
options = ASB2008CDC.Options()
14+
options.t_end = 10.0
15+
options.step = 0.05
16+
options.tensor_order = 2
17+
options.taylor_terms = 4
18+
options.u = Zonotope.zero(1, 1)
19+
options.u_trans = np.zeros(1)
20+
21+
# settings for the using geometry
22+
Zonotope.REDUCE_METHOD = Zonotope.REDUCE_METHOD.GIRARD
23+
Zonotope.ORDER = 50
24+
25+
x = codac.VectorVar(2)
26+
f = codac.AnalyticFunction([x], [(0.31) * 1 + (2.57) * x[0] + (-0.00) * x[1] + (-4.58) * x[0] ** 2 + (-2.71) * x[0] * x[1] + (13.63) * x[1] ** 2 + (-9.76) * x[0] ** 3 + (-1.94) * x[0] ** 2 * x[1] + (-6.82) * x[0] * x[1] ** 2 + (0.75) * x[1] ** 3 + (-16.03) * x[0] ** 4 + (7.56) * x[0] ** 3 * x[1] + (15.39) * x[0] ** 2 * x[1] ** 2 + (6.70) * x[0] * x[1] ** 3 + (-83.93) * x[1] ** 4])
27+
init_box = Interval([-0.5, -0.5], [0.5, 0.5])
28+
29+
this_time = performance_counter_start()
30+
xs = extract_boundary(init_interval=init_box, init_set=f, eps=0.04)
31+
ri_with_bound, rp_with_bound = ASB2008CDC.reach_parallel(brusselator, [2, 1], options, xs)
32+
this_time = performance_counter(this_time, 'reach_without_bound')
33+
34+
# visualize the results
35+
plot(ri_with_bound, [0, 1], c="#303F9F")
36+
37+
boundary_boxes = []
38+
for i in xs:
39+
boundary_boxes.append(cvt2(i, Geometry.TYPE.INTERVAL))
40+
constraint = ("(0.31)*1 + (2.57)*x + (-0.00)*y + (-4.58)*x^2 + (-2.71)*x y + (13.63)*y^2 + (-9.76)*x^3 + "
41+
"(-1.94)*x^2 y + (-6.82)*x y^2 + (0.75)*y^3 + (-16.03)*x^4 + (7.56)*x^3 y + (15.39)*x^2 y^2 + "
42+
"(6.70)*x y^3 + (-83.93)*y^4")
43+
constraint = constraint.replace(" y", "*y")
44+
plot(boundary_boxes, [0, 1], xlim=[-0.5, 0.5], ylim=[-0.5, 0.5], init_set=constraint)
45+

‎benchmarks/benchmark_brusselator_polytope.py‎

Lines changed: 50 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,50 @@
1+
import numpy as np
2+
import codac
3+
from pybdr.algorithm import ASB2008CDC
4+
from pybdr.geometry import Zonotope, Interval, Polytope, Geometry
5+
from pybdr.geometry.operation import cvt2
6+
from pybdr.model import *
7+
from pybdr.util.visualization import plot_cmp, plot
8+
from pybdr.util.functional import extract_boundary
9+
from pybdr.util.functional import performance_counter, performance_counter_start
10+
11+
if __name__ == '__main__':
12+
# settings for the computation
13+
options = ASB2008CDC.Options()
14+
options.t_end = 10.0
15+
options.step = 0.05
16+
options.tensor_order = 2
17+
options.taylor_terms = 4
18+
options.u = Zonotope.zero(1, 1)
19+
options.u_trans = np.zeros(1)
20+
21+
# settings for the using geometry
22+
Zonotope.REDUCE_METHOD = Zonotope.REDUCE_METHOD.GIRARD
23+
Zonotope.ORDER = 50
24+
25+
A = np.array([[ 0.84680084, 0.53191008],
26+
[ 0.17672503, 0.98426026],
27+
[-0.68443591, 0.72907303],
28+
[-0.99997959, -0.00638965],
29+
[-0.92101281, -0.3895323 ],
30+
[-0.39691115, -0.91785704],
31+
[ 0.60973787, -0.79260314],
32+
[ 0.93952577, -0.34247821]])
33+
b = np.array([0.40646441, 0.43238234, 0.45970287, 0.41445799,
34+
0.40973516, 0.44057138, 0.38044951, 0.45097237], dtype=float)
35+
init_set = Polytope(A, b)
36+
init_box = codac.IntervalVector([[-0.5, 0.5], [-0.5, 0.5]])
37+
38+
this_time = performance_counter_start()
39+
xs = extract_boundary(init_box, init_set, eps=0.04)
40+
ri_with_bound, rp_with_bound = ASB2008CDC.reach_parallel(brusselator, [2, 1], options, xs)
41+
this_time = performance_counter(this_time, 'reach_without_bound')
42+
43+
# visualize the results
44+
plot(ri_with_bound, [0, 1], c="#303F9F")
45+
46+
boundary_boxes = []
47+
for i in xs:
48+
boundary_boxes.append(cvt2(i, Geometry.TYPE.INTERVAL))
49+
plot(boundary_boxes, [0, 1], init_set=init_set)
50+

‎benchmarks/benchmark_brusselator_zonotope.py‎

Lines changed: 44 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,44 @@
1+
import numpy as np
2+
import codac
3+
from pybdr.algorithm import ASB2008CDC
4+
from pybdr.geometry import Zonotope, Interval, Polytope, Geometry
5+
from pybdr.geometry.operation import cvt2
6+
from pybdr.model import *
7+
from pybdr.util.visualization import plot_cmp, plot
8+
from pybdr.util.functional import extract_boundary
9+
from pybdr.util.functional import performance_counter, performance_counter_start
10+
11+
if __name__ == '__main__':
12+
# settings for the computation
13+
options = ASB2008CDC.Options()
14+
options.t_end = 10.0
15+
options.step = 0.05
16+
options.tensor_order = 2
17+
options.taylor_terms = 4
18+
options.u = Zonotope.zero(1, 1)
19+
options.u_trans = np.zeros(1)
20+
21+
# settings for the using geometry
22+
Zonotope.REDUCE_METHOD = Zonotope.REDUCE_METHOD.GIRARD
23+
Zonotope.ORDER = 50
24+
25+
data = np.array(
26+
[[0.0, 0.24, -0.10, 0.03, -0.08], # x
27+
[0.0, 0.05, 0.18, -0.12, -0.04]],
28+
dtype=float,
29+
)
30+
init_set = Zonotope(data[:, 0], data[:, 1:])
31+
init_box = codac.IntervalVector([[-0.5, 0.5], [-0.5, 0.5]])
32+
33+
this_time = performance_counter_start()
34+
xs = extract_boundary(init_box, init_set, eps=0.04)
35+
ri_with_bound, rp_with_bound = ASB2008CDC.reach_parallel(brusselator, [2, 1], options, xs)
36+
this_time = performance_counter(this_time, 'reach_without_bound')
37+
38+
# visualize the results
39+
plot(ri_with_bound, [0, 1], c="#303F9F")
40+
41+
boundary_boxes = []
42+
for i in xs:
43+
boundary_boxes.append(cvt2(i, Geometry.TYPE.INTERVAL))
44+
plot(boundary_boxes, [0, 1], init_set=init_set)

‎pybdr/util/functional/__init__.py‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -3,6 +3,7 @@
33
from .kd_tree import *
44
from .realpaver_wrapper import RealPaver
55
from .simulator import Simulator
6+
from .codac_wrapper import extract_boundary
67

78
__all__ = [
89
"is_empty",

‎pybdr/util/functional/codac_wrapper.py‎

Lines changed: 129 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,129 @@
1+
from __future__ import annotations
2+
import codac
3+
4+
5+
def extract_boundary(init_interval, init_set, eps: float = 0.04,
6+
simplify_result: bool = True, simplify_eps: float = 0.005):
7+
from pybdr.geometry import Geometry, Interval, Polytope, Zonotope
8+
from pybdr.geometry.operation.convert import cvt2
9+
10+
if isinstance(init_set, codac.AnalyticFunction):
11+
return extract_boundary_analytic_function(init_interval, init_set, eps, simplify_result, simplify_eps)
12+
elif isinstance(init_set, Polytope) or isinstance(init_set, Zonotope):
13+
if isinstance(init_set, Zonotope):
14+
init_set = cvt2(init_set, Geometry.TYPE.POLYTOPE)
15+
return extract_boundary_polytope(init_interval, init_set, eps)
16+
else:
17+
raise NotImplementedError()
18+
19+
20+
def extract_boundary_polytope(init_interval, init_set, eps: float = 0.04):
21+
from pybdr.geometry import Geometry, Polytope, Interval, Zonotope
22+
from pybdr.geometry.operation.convert import cvt2
23+
assert isinstance(init_set, Polytope)
24+
n = init_set.a.shape[1]
25+
m = init_set.a.shape[0]
26+
27+
x = codac.VectorVar(n)
28+
29+
funcs = [sum(init_set.a[i, j] * x[j] for j in range(n)) - init_set.b[i] for i in range(m)]
30+
g = codac.AnalyticFunction([x], funcs)
31+
32+
Y = codac.IntervalVector([[float("-inf"), 0.0]] * m)
33+
ctc = codac.CtcInverse(g, Y)
34+
35+
boxes = []
36+
pave_boxes(init_interval, ctc, eps, boxes)
37+
38+
boundary_boxes = []
39+
for box in boxes:
40+
if box.is_empty():
41+
continue
42+
vals = g.eval(box)
43+
for i in range(m):
44+
if vals[i].lb() <= 0.0 <= vals[i].ub():
45+
boundary_boxes.append(box)
46+
break
47+
48+
x_init = []
49+
for i in boundary_boxes:
50+
if not i.is_empty():
51+
i_interval = Interval(i.lb(), i.ub())
52+
# x_init.append(i_interval)
53+
x_init.append(cvt2(i_interval, Geometry.TYPE.ZONOTOPE))
54+
return x_init
55+
56+
57+
def extract_boundary_analytic_function(init_interval, init_set, eps: float = 0.04,
58+
simplify_result: bool = True, simplify_eps: float = 0.005):
59+
from pybdr.geometry import Geometry, Interval
60+
from pybdr.geometry.operation.convert import cvt2
61+
62+
ctc = codac.CtcInverse(init_set, [0.0])
63+
# ctc = codac.CtcInverse(init_set_func, codac.IntervalVector([codac.Interval(0, 1e9)]))
64+
65+
init_sub_interval_list = []
66+
for i in range(init_interval.shape[0]):
67+
init_sub_interval_list.append(codac.Interval(init_interval[i].inf[0], init_interval[i].sup[0]))
68+
init_domain = codac.IntervalVector(init_sub_interval_list)
69+
result = []
70+
pave_boxes(init_domain, ctc, eps, result)
71+
72+
if simplify_result:
73+
isolated_boxes = extract_isolated_boxes(result)
74+
75+
for i in isolated_boxes:
76+
boxes_i = []
77+
i_copy = codac.IntervalVector(i)
78+
pave_boxes(i_copy, ctc, simplify_eps, boxes_i)
79+
80+
result.remove(i)
81+
result += boxes_i
82+
83+
x_init = []
84+
for i in result:
85+
if not i.is_empty():
86+
i_interval = Interval(i.lb(), i.ub())
87+
# x_init.append(i_interval)
88+
x_init.append(cvt2(i_interval, Geometry.TYPE.ZONOTOPE))
89+
return x_init
90+
91+
92+
def pave_boxes(box, ctc: codac.CtcInverse, eps, result):
93+
ctc.contract(box)
94+
if box.is_empty():
95+
return
96+
97+
if box.max_diam() <= eps:
98+
result.append(box)
99+
return
100+
101+
i = box.max_diam_index()
102+
left, right = box.bisect(i)
103+
pave_boxes(left, ctc, eps, result)
104+
pave_boxes(right, ctc, eps, result)
105+
106+
107+
def boxes_touch(b1: codac.IntervalVector, b2: codac.IntervalVector, tol=1e-15):
108+
n = b1.size()
109+
for i in range(n):
110+
if b1[i].ub() < b2[i].lb() - tol or b2[i].ub() < b1[i].lb() - tol:
111+
return False
112+
return True
113+
114+
115+
def extract_isolated_boxes(boxes, tol=1e-12):
116+
isolated = []
117+
for i, b in enumerate(boxes):
118+
connected = False
119+
for j, b2 in enumerate(boxes):
120+
if i == j:
121+
continue
122+
if boxes_touch(b, b2, tol):
123+
connected = True
124+
break
125+
if not connected:
126+
isolated.append(b)
127+
return isolated
128+
129+

‎pybdr/util/visualization/plot.py‎

Lines changed: 36 additions & 7 deletions
Original file line numberDiff line numberDiff line change
@@ -23,7 +23,7 @@ def __2d_add_interval(ax, i: "Interval", dims, color, filled):
2323
closed=True,
2424
alpha=1,
2525
fill=filled,
26-
linewidth=5,
26+
linewidth=3,
2727
edgecolor=color,
2828
facecolor=color,
2929
)
@@ -67,11 +67,39 @@ def __2d_plot(
6767
ylim=None,
6868
c=None,
6969
filled=False,
70+
init_set=None
7071
):
7172
assert len(dims) == 2
7273
px = 1 / plt.rcParams["figure.dpi"]
7374
fig, ax = plt.subplots(figsize=(width * px, height * px), layout="constrained")
7475

76+
if init_set is not None:
77+
if isinstance(init_set, Polytope) or isinstance(init_set, Zonotope) or isinstance(init_set, Interval):
78+
if init_set.type == Geometry.TYPE.INTERVAL:
79+
__2d_add_interval(ax, init_set, dims, "#F4B083", "True")
80+
elif init_set.type == Geometry.TYPE.POLYTOPE:
81+
__2d_add_polytope(ax, init_set, dims, "#F4B083", "True")
82+
elif init_set.type == Geometry.TYPE.ZONOTOPE:
83+
__2d_add_zonotope(ax, init_set, dims, "#F4B083", "True")
84+
elif isinstance(init_set, str):
85+
assert xlim is not None
86+
assert ylim is not None
87+
mesh_num = 400
88+
import sympy
89+
x, y = sympy.symbols('x y')
90+
poly_expr = sympy.sympify(init_set.replace('^', '**'))
91+
f = sympy.lambdify((x, y), poly_expr, 'numpy')
92+
93+
X = np.linspace(xlim[0], xlim[1], mesh_num)
94+
Y = np.linspace(ylim[0], ylim[1], mesh_num)
95+
XX, YY = np.meshgrid(X, Y)
96+
ZZ = f(XX, YY)
97+
98+
ax.contourf(XX, YY, ZZ, levels=[0, ZZ.max()], colors=['#F4B083'], alpha=1)
99+
# ax.contour(XX, YY, ZZ, levels=[0], colors='#ffcccc', linewidths=3.0)
100+
else:
101+
raise NotImplementedError
102+
75103
for i in range(len(objs)):
76104
this_color = plt.cm.turbo(i / len(objs)) if c is None else c
77105
geos = (
@@ -94,8 +122,8 @@ def __2d_plot(
94122
else:
95123
raise NotImplementedError
96124

97-
ax.autoscale_view()
98-
# ax.axis("equal")
125+
# ax.autoscale_view()
126+
ax.axis("equal")
99127
ax.set_xlabel("x" + str(dims[0]))
100128
ax.set_ylabel("x" + str(dims[1]))
101129

@@ -120,9 +148,10 @@ def plot(
120148
ylim=None,
121149
c=None,
122150
filled=False,
151+
init_set=None
123152
):
124153
if mod == "2d":
125-
return __2d_plot(objs, dims, width, height, xlim, ylim, c, filled)
154+
return __2d_plot(objs, dims, width, height, xlim, ylim, c, filled, init_set)
126155
elif mod == "3d":
127156
return __3d_plot(objs, dims, width, height)
128157
else:
@@ -162,8 +191,8 @@ def __2d_plot_cmp(collections, dims, width, height, xlim, ylim, cs, filled, show
162191
else:
163192
raise NotImplementedError
164193

165-
ax.autoscale_view()
166-
# ax.axis("equal")
194+
# ax.autoscale_view()
195+
ax.axis("equal")
167196
ax.set_xlabel("x" + str(dims[0]))
168197
ax.set_ylabel("x" + str(dims[1]))
169198

@@ -193,7 +222,7 @@ def plot_cmp(
193222
ylim=None,
194223
cs=None,
195224
filled=False,
196-
show=False,
225+
show=True,
197226
save_file_name=None
198227
):
199228
if mod == "2d":

0 commit comments

Comments
 (0)