Jc/remove general domain (#419)

* make periodic_domain working in UniformGrid solver

* remove GeneralDomain

* remove comments

* remove comments from hierarchic

* Fix format

* Update src/domain/grid_domain/nested_domain.jl

Co-authored-by: Benoît Legat <[email protected]>

* Update src/domain/grid_domain/nested_domain.jl

Co-authored-by: Benoît Legat <[email protected]>

---------

Co-authored-by: Benoît Legat <[email protected]>

Author: JulienCalbert

docs/src/examples/solvers/Hierarchical-abstraction.jl

@@ -22,18 +22,17 @@ include(joinpath(dirname(dirname(pathof(Dionysos))), "problems", "simple_problem
 ## specific functions
 function post_image(abstract_system, concrete_system, xpos, u)
     Xdom = abstract_system.Xdom
-    x = DO.get_coord_by_pos(Xdom.grid, xpos)
+    x = DO.get_coord_by_pos(Xdom, xpos)
     Fx = concrete_system.f_eval(x, u)
-    r = Xdom.grid.h / 2.0 + concrete_system.measnoise
+    r = DO.get_grid(Xdom).h / 2.0 + concrete_system.measnoise
     Fr = r
 
-    rectI = DO.get_pos_lims_outer(Xdom.grid, UT.HyperRectangle(Fx .- Fr, Fx .+ Fr))
+    rectI = DO.get_pos_lims_outer(DO.get_grid(Xdom), UT.HyperRectangle(Fx .- Fr, Fx .+ Fr))
     ypos_iter = Iterators.product(DO._ranges(rectI)...)
     over_approx = []
     allin = true
     for ypos in ypos_iter
-        ypos = DO.set_in_period_pos(Xdom, ypos)
-        if !(ypos in Xdom)
+        if !(ypos in abstract_system)
             allin = false
             break
         end
@@ -44,17 +43,16 @@ function post_image(abstract_system, concrete_system, xpos, u)
 end
 
 function pre_image(abstract_system, concrete_system, xpos, u)
-    grid = abstract_system.Xdom.grid
-    x = DO.get_coord_by_pos(grid, xpos)
+    Xdom = abstract_system.Xdom
+    x = DO.get_coord_by_pos(Xdom, xpos)
     potential = Int[]
     x_prev = concrete_system.f_backward(x, u)
-    xpos_cell = DO.get_pos_by_coord(grid, x_prev)
+    xpos_cell = DO.get_pos_by_coord(Xdom, x_prev)
     n = 2
     for i in (-n):n
         for j in (-n):n
             x_n = (xpos_cell[1] + i, xpos_cell[2] + j)
-            x_n = DO.set_in_period_pos(abstract_system.Xdom, x_n)
-            if x_n in abstract_system.Xdom
+            if x_n in abstract_system
                 cell = SY.get_state_by_xpos(abstract_system, x_n)[1]
                 if !(cell in potential)
                     push!(potential, cell)
@@ -87,9 +85,9 @@ minimum_transition_cost(symmodel, contsys, source, target) = 1.0
 concrete_problem = SimpleProblem.problem(;
     rectX = UT.HyperRectangle(SVector(0.0, 0.0), SVector(60.0, 60.0)),
     obstacles = [UT.HyperRectangle(SVector(22.0, 21.0), SVector(25.0, 32.0))],
-    periodic = Int[],
-    periods = [30.0, 30.0],
-    T0 = [0.0, 0.0],
+    periodic = SVector{0, Int}(),
+    periods = SVector{0, Float64}(),
+    T0 = SVector{0, Float64}(),
     rectU = UT.HyperRectangle(SVector(-2.0, -2.0), SVector(2.0, 2.0)),
     Uobstacles = [UT.HyperRectangle(SVector(-0.5, -0.5), SVector(0.5, 0.5))],
     _I_ = UT.HyperRectangle(SVector(6.5, 6.5), SVector(7.5, 7.5)),
@@ -103,8 +101,8 @@ concrete_problem = SimpleProblem.problem(;
 concrete_system = concrete_problem.system;
 
 # Local optimizer parameters
-hx_local = [0.5, 0.5]
-hx_heuristic = [1.0, 1.0]
+hx_local = SVector(0.5, 0.5)
+hx_heuristic = SVector(1.0, 1.0)
 u0 = SVector(0.0, 0.0)
 hu = SVector(0.5, 0.5)
 Ugrid = DO.GridFree(u0, hu)
@@ -124,7 +122,7 @@ AB.LazyAbstraction.set_optimizer_parameters!(
 )
 
 # Global optimizer parameters
-hx_global = [10.0, 10.0]
+hx_global = SVector(10.0, 10.0)
 u0 = SVector(0.0, 0.0)
 hu = SVector(0.5, 0.5)
 Ugrid = DO.GridFree(u0, hu)
@@ -168,7 +166,7 @@ println("Cost:\t $(cost)")
 
 # ## Display the results
 # # Display the specifications, domains and trajectory
-fig = plot(; aspect_ratio = :equal);
+fig1 = plot(; aspect_ratio = :equal);
 
 #We display the concrete domain
 plot!(concrete_system.X; color = :grey, opacity = 0.5, label = "");
@@ -184,7 +182,7 @@ plot!(concrete_problem.target_set; dims = [1, 2], color = :red, opacity = 0.8);
 plot!(cost_control_trajectory; ms = 0.5)
 
 # # Display the lazy abstraction 
-fig = plot(; aspect_ratio = :equal);
+fig2 = plot(; aspect_ratio = :equal);
 
 plot!(
     optimizer.hierarchical_problem;

docs/src/examples/solvers/Lazy-abstraction-reachability.jl

@@ -34,18 +34,17 @@ include(joinpath(dirname(dirname(pathof(Dionysos))), "problems", "simple_problem
 ## specific functions
 function post_image(abstract_system, concrete_system, xpos, u)
     Xdom = abstract_system.Xdom
-    x = DO.get_coord_by_pos(Xdom.grid, xpos)
+    x = DO.get_coord_by_pos(Xdom, xpos)
     Fx = concrete_system.f_eval(x, u)
-    r = Xdom.grid.h / 2.0 + concrete_system.measnoise
+    r = DO.get_grid(Xdom).h / 2.0 + concrete_system.measnoise
     Fr = r
 
-    rectI = DO.get_pos_lims_outer(Xdom.grid, UT.HyperRectangle(Fx .- Fr, Fx .+ Fr))
+    rectI = DO.get_pos_lims_outer(DO.get_grid(Xdom), UT.HyperRectangle(Fx .- Fr, Fx .+ Fr))
     ypos_iter = Iterators.product(DO._ranges(rectI)...)
     over_approx = []
     allin = true
     for ypos in ypos_iter
-        ypos = DO.set_in_period_pos(Xdom, ypos)
-        if !(ypos in Xdom)
+        if !(ypos in abstract_system)
             allin = false
             break
         end
@@ -56,17 +55,16 @@ function post_image(abstract_system, concrete_system, xpos, u)
 end
 
 function pre_image(abstract_system, concrete_system, xpos, u)
-    grid = abstract_system.Xdom.grid
-    x = DO.get_coord_by_pos(grid, xpos)
+    Xdom = abstract_system.Xdom
+    x = DO.get_coord_by_pos(Xdom, xpos)
     potential = Int[]
     x_prev = concrete_system.f_backward(x, u)
-    xpos_cell = DO.get_pos_by_coord(grid, x_prev)
+    xpos_cell = DO.get_pos_by_coord(Xdom, x_prev)
     n = 2
     for i in (-n):n
         for j in (-n):n
             x_n = (xpos_cell[1] + i, xpos_cell[2] + j)
-            x_n = DO.set_in_period_pos(abstract_system.Xdom, x_n)
-            if x_n in abstract_system.Xdom
+            if x_n in abstract_system
                 cell = SY.get_state_by_xpos(abstract_system, x_n)[1]
                 if !(cell in potential)
                     push!(potential, cell)
@@ -98,11 +96,11 @@ minimum_transition_cost(symmodel, contsys, source, target) = 1.0
 concrete_problem = SimpleProblem.problem()
 concrete_system = concrete_problem.system
 
-hx = [0.5, 0.5]
+hx = SVector(0.5, 0.5)
 u0 = SVector(0.0, 0.0)
 hu = SVector(0.5, 0.5)
 Ugrid = DO.GridFree(u0, hu)
-hx_heuristic = [1.0, 1.0] * 1.5
+hx_heuristic = SVector(1.5, 1.5)
 maxIter = 100
 
 optimizer = MOI.instantiate(AB.LazyAbstraction.Optimizer)

docs/src/examples/solvers/Path planning.jl

@@ -1,4 +1,3 @@
-using Test     #src
 # # Example: Path planning problem solved by [Uniform grid abstraction](https://github.com/dionysos-dev/Dionysos.jl/blob/master/docs/src/manual/manual.md#solvers).
 #
 #md # [![Binder](https://mybinder.org/badge_logo.svg)](@__BINDER_ROOT_URL__/generated/Path planning.ipynb)
@@ -97,13 +96,13 @@ set_attribute(
 set_attribute(model, "efficient", true)
 
 x0 = SVector(0.0, 0.0, 0.0);
-h = SVector(0.2, 0.2, 0.2);
-set_attribute(model, "state_grid", Dionysos.Domain.GridFree(x0, h))
+hx = SVector(0.2, 0.2, 0.2);
+set_attribute(model, "state_grid", Dionysos.Domain.GridFree(x0, hx))
 
 # Definition of the grid of the input-space on which the abstraction is based (origin `u0` and input-space discretization `h`):
 u0 = SVector(0.0, 0.0);
-h = SVector(0.3, 0.3);
-set_attribute(model, "input_grid", Dionysos.Domain.GridFree(u0, h))
+hu = SVector(0.3, 0.3);
+set_attribute(model, "input_grid", Dionysos.Domain.GridFree(u0, hu))
 
 optimize!(model);
 
@@ -124,8 +123,6 @@ println("Time to solve the abstract problem: $(abstract_problem_time)")
 total_time = MOI.get(model, MOI.RawOptimizerAttribute("solve_time_sec"))
 println("Total time: $(total_time)")
 
-@test length(abstract_controller.data) == 19400 #src
-
 # ### Trajectory display
 # We choose a stopping criterion `reached` and the maximal number of steps `nsteps` for the sampled system, i.e. the total elapsed time: `nstep`*`tstep`
 # as well as the true initial state `x0` which is contained in the initial state-space `_I_` defined previously.

docs/src/examples/utils/Grid.jl

@@ -67,16 +67,6 @@ function plot_deformed_grid_with_DomainList(f, fi)
     return plot!(dom; show = true, color = :grey, efficient = false)
 end
 
-function plot_deformed_grid_with_GeneralDomain(f, fi)
-    X = UT.HyperRectangle(SVector(0.0, 0.0), SVector(30.0, 10.0))
-    obstacle = UT.HyperRectangle(SVector(10.0, 10.0), SVector(15.0, 15.0))
-    hx = [2.5, 2.5]
-    d = DO.RectangularObstacles(X, [obstacle])
-    dom = DO.GeneralDomainList(hx; elems = d, f = f, fi = fi, fit = true)
-    plot(; aspect_ratio = :equal)
-    return plot!(dom; show = true, color = :grey, efficient = false)
-end
-
 rect = UT.HyperRectangle(SVector(0.0, 0.0), SVector(2.0, 2.0))
 shape = UT.DeformedRectangle(rect, f2)
 plot(; aspect_ratio = :equal)
@@ -85,53 +75,66 @@ plot!(shape; color = :red, efficient = false, label = "Deformed")
 
 # Display some deformed Grids
 plot_deformed_grid_with_DomainList(f1, fi1)
-plot_deformed_grid_with_GeneralDomain(f1, fi1)
 
 plot_deformed_grid_with_DomainList(f2, fi2)
-plot_deformed_grid_with_GeneralDomain(f2, fi2)
 
 plot_deformed_grid_with_DomainList(f3, fi3)
-plot_deformed_grid_with_GeneralDomain(f3, fi3)
 
 plot_deformed_grid_with_DomainList(f4, fi4)
-plot_deformed_grid_with_GeneralDomain(f4, fi4)
 
 f, fi = build_f_rotation(π / 3.0)
 plot_deformed_grid_with_DomainList(f, fi)
 
-# ### Nested classical grid
+# ### Nested domain
+# --- Base domain and obstacle ---
 X = UT.HyperRectangle(SVector(0.0, 0.0), SVector(30.0, 30.0))
 obstacle = UT.HyperRectangle(SVector(15.0, 15.0), SVector(20.0, 20.0))
-hx = [3.0, 1.0] * 2.0
-periodic = Int[1]
-periods = [30.0, 30.0]
-T0 = [0.0, 0.0]
-d = DO.RectangularObstacles(X, [obstacle])
-dom = DO.GeneralDomainList(
-    hx;
-    elems = d,
-    periodic = periodic,
-    periods = periods,
-    T0 = T0,
-    fit = true,
-)
-Ndomain = DO.NestedDomain(dom)
-DO.cut_pos!(Ndomain, (2, 2), 1)
-DO.cut_pos!(Ndomain, (2, 3), 1)
-DO.cut_pos!(Ndomain, (4, 4), 2)
-
-fig = plot(; aspect_ratio = :equal, legend = false);
+
+# Grid settings
+hx = SVector(6.0, 2.0)
+periodic_dims = SVector(1)
+periods = SVector(30.0)
+start = SVector(0.0)
+
+# Create a LazySetMinus: free space = X minus obstacle
+free_space = UT.LazySetMinus(X, UT.LazyUnionSetArray([obstacle]))
+
+# Create periodic grid domain
+domain = DO.PeriodicDomainList(periodic_dims, periods, start, hx)
+# domain = DO.DomainList(hx)
+
+DO.add_set!(domain, free_space, DO.INNER)
+
+# Create the NestedDomain
+Ndomain = DO.NestedDomain(domain)
+
+# --- Refinement ---
+# Refine certain cells
+div = 3
+DO.cut_pos!(Ndomain, (2, 2), 1; div = div)
+DO.cut_pos!(Ndomain, (2, 3), 1; div = div)
+DO.cut_pos!(Ndomain, (4, 4), 2; div = div)
+
+# --- Plotting ---
+fig = plot(; aspect_ratio = :equal, legend = false)
+plot!(free_space; color = :blue, label = "Base domain", efficient = false)
 plot!(Ndomain; color = :grey, efficient = false)
 
-# ### Ellipsoidal partition based on a grid
+# ----------------------------------------------------------------------
+# 2D Example: Ellipsoidal partition using a grid
+
 x0 = SVector(0.0, 0.0)
 n_step = 2
 h = SVector(1.0 / n_step, 1.0 / n_step)
 P = 0.5 * diagm((h ./ 2) .^ (-2))
-rectX = UT.HyperRectangle(SVector(-2, -2), SVector(2, 2))
 
-grid = DO.GridEllipsoidalRectangular(x0, h, P)
-domain = DO.DomainList(grid)
-DO.add_set!(domain, rectX, DO.OUTER)
-plot(; aspect_ratio = :equal);
-plot!(domain; color = :grey, opacity = 0.5, efficient = false)
+rectX = UT.HyperRectangle(SVector(-2.0, -2.0), SVector(2.0, 2.0))
+
+# Create an ellipsoidal grid domain
+ellip_grid = DO.GridEllipsoidalRectangular(x0, h, P)
+ellip_domain = DO.DomainList(ellip_grid)
+DO.add_set!(ellip_domain, rectX, DO.OUTER)
+
+# Plot
+fig2 = plot(; aspect_ratio = :equal)
+plot!(ellip_domain; color = :grey, opacity = 0.5, efficient = false)

docs/src/reference/Domain.md

@@ -23,14 +23,32 @@ Dionysos.Domain.get_volume
 ```@docs
 Dionysos.Domain.DomainType
 Dionysos.Domain.GridDomainType
+```
+
+```@docs
 Dionysos.Domain.merge_to_hyperrectangles_pos
 Dionysos.Domain.merge_to_hyperrectangles_real
+```
+
+## Concrete domains 
+
+### Classical domains
+```@docs
 Dionysos.Domain.DomainList
+Dionysos.Domain.CustomList
+```
+
+### Periodic domain 
+```@docs
 Dionysos.Domain.PeriodicDomainList
+Dionysos.Domain.is_periodic
 Dionysos.Domain.has_same_periodicity
 Dionysos.Domain.wrap_pos
 Dionysos.Domain.wrap_coord
-Dionysos.Domain.GeneralDomainList
-Dionysos.Domain.RectangularObstacles
-Dionysos.Domain.CustomList
+Dionysos.Domain._make_periodic_index_map
+```
+
+### Nested domain 
+```@docs
+Dionysos.Domain.NestedDomain
 ```

docs/src/reference/System.md

@@ -114,6 +114,7 @@ Dionysos.System.AffineController
 
 ## Trajectories 
 ```@docs
+Dionysos.System.wrap_coord
 Dionysos.System.DiscreteTrajectory
 Dionysos.System.ContinuousTrajectory
 Dionysos.System.HybridTrajectory

docs/src/reference/Utils.md

@@ -46,6 +46,7 @@ Dionysos.Utils.Tree
 ## Geometric shapes
 
 ```@docs
+Dionysos.Utils.set_in_period 
 Dionysos.Utils.HyperRectangle
 Dionysos.Utils.DeformedRectangle
 ```

problems/simple_problem.jl

@@ -66,9 +66,9 @@ end
 function problem(;
     rectX = UT.HyperRectangle(SVector(0.0, 0.0), SVector(30.0, 30.0)),
     obstacles = [UT.HyperRectangle(SVector(15.0, 15.0), SVector(20.0, 20.0))],
-    periodic = Int[],
-    periods = [30.0, 30.0],
-    T0 = [0.0, 0.0],
+    periodic = SVector{0, Int}(),
+    periods = SVector{0, Float64}(),
+    T0 = SVector{0, Float64}(),
     rectU = UT.HyperRectangle(SVector(-2.0, -2.0), SVector(2.0, 2.0)),
     Uobstacles = [UT.HyperRectangle(SVector(-0.5, -0.5), SVector(0.5, 0.5))],
     _I_ = UT.HyperRectangle(SVector(5.0, 5.0), SVector(6.0, 6.0)),

src/domain/domain.jl

@@ -20,7 +20,6 @@ include("grid_domain/grid.jl")
 include("grid_domain/grid_domain.jl")
 include("grid_domain/domain_list.jl")
 include("grid_domain/periodic_domain.jl")
-include("grid_domain/general_domain.jl")
 include("grid_domain/nested_domain.jl")
 
 include("continuous_domain.jl")