diff --git a/doc/manual/index.rst b/doc/manual/index.rst index 29fe3a84a..50cc9462f 100644 --- a/doc/manual/index.rst +++ b/doc/manual/index.rst @@ -283,6 +283,7 @@ User manual * SepNot * SepAction * Analytic separators + * :ref:`sec-ctc-analytic-sepimage` * SepInverse * SepTransform * Geometrical separators diff --git a/doc/manual/manual/contractors/analytic/index.rst b/doc/manual/manual/contractors/analytic/index.rst index d0e781ade..65b9a5d3f 100644 --- a/doc/manual/manual/contractors/analytic/index.rst +++ b/doc/manual/manual/contractors/analytic/index.rst @@ -4,5 +4,6 @@ Analytic contractors .. toctree:: Directed operators + sepimage.rst ctcinverse.rst - CtcInverseNotIn \ No newline at end of file + CtcInverseNotIn diff --git a/doc/manual/manual/contractors/analytic/sep_image.png b/doc/manual/manual/contractors/analytic/sep_image.png new file mode 100644 index 000000000..a11e6806f Binary files /dev/null and b/doc/manual/manual/contractors/analytic/sep_image.png differ diff --git a/doc/manual/manual/contractors/analytic/sep_image_fine.png b/doc/manual/manual/contractors/analytic/sep_image_fine.png new file mode 100644 index 000000000..81258cf78 Binary files /dev/null and b/doc/manual/manual/contractors/analytic/sep_image_fine.png differ diff --git a/doc/manual/manual/contractors/analytic/sepimage.rst b/doc/manual/manual/contractors/analytic/sepimage.rst new file mode 100644 index 000000000..371b31e92 --- /dev/null +++ b/doc/manual/manual/contractors/analytic/sepimage.rst @@ -0,0 +1,176 @@ +.. _sec-ctc-analytic-sepimage: + +The SepImage separator +====================== + + Main author: `Maël Godard `_ + +Definition +---------- + +Consider a function :math:`\mathbf{f}:\mathbb{R}^n\to \mathbb{R}^p`. +The ``SepImage`` separator allows one to handle constraints of the form :math:`\mathbf{y}\in\mathbf{f}(\mathbb{X})` +by separating input boxes :math:`[\mathbf{y}]\in\mathbb{IR}^p`. + +Construction and basic usage +---------------------------- + +The ``SepImage`` relies on a boundary approach to compute the image set. The boundary of the initial set needs to be +covered by a :ref:`gnomonic atlas `. To compute the image of the boundary, +the initial box :math:`\left[-1,1\right]^m` can be subdivided up to a size :math:`\epsilon` and the image of each resulting box +is computed. Note that if one chooses to take :math:`\epsilon=2`, only one computation will be done per chart of the atlas. + +Once this image of the boundary has been computed, the ``SepImage`` needs to a way to characterize if a point is inside or +outside of the image set. To do so, it looks for an antecedent of the point in the initial set. An additionnal contractor on +the initial set is then required. + +The typical workflow is: + +1. Define analytic variables (scalar, vector, matrix) associated with the domain of the function. +2. Build an :class:`~codac2.AnalyticFunction`. +3. Define a gnomonic atlas on the boundary of the initial set +4. Define a contractor on this same initial set +5. Instantiate ``SepImage`` with the atlas, the function, the resolution :math:`\epsilon` and the contractor. +6. Contract an input box :math:`[\mathbf{y}]` or pave the separator in the image space. + +Example +------- + +Consider that we want to construct a separator on the image of the unit disk by the function + +.. math:: + \mathbf{f}(\mathbf{x})= + \left( + \begin{array}{c} + 3(x_{1}+1)\\ + x_{2}+ 0.5\sin (3 x_{1}) + \end{array} + \right) + +This function can be constructed in codac as follows + +.. tabs:: + + .. group-tab:: Python + + .. literalinclude:: src.py + :language: py + :start-after: [sepimage-1-beg] + :end-before: [sepimage-1-end] + :dedent: 4 + + .. group-tab:: C++ + + .. literalinclude:: src.cpp + :language: c++ + :start-after: [sepimage-1-beg] + :end-before: [sepimage-1-end] + :dedent: 4 + +We then need a contractor for the unit disk, and a gnomonic atlas for its boundary (the unit circle). +For the atlas, the image of :math:`\left[-1,1\right]` by the function + +.. math:: + \psi_{0} = + \left( + \begin{array}{c} + \cos\left(x\cdot\frac{\pi}{2}\right)\\ + \sin\left(x\cdot\frac{\pi}{2}\right) + \end{array} + \right) + +is half of the unit cicle. The other half can be obtained with a rotation of :math:`\pi` rad. Such atlas +is constucted in codac as + +.. tabs:: + + .. group-tab:: Python + + .. literalinclude:: src.py + :language: py + :start-after: [sepimage-2-beg] + :end-before: [sepimage-2-end] + :dedent: 4 + + .. group-tab:: C++ + + .. literalinclude:: src.cpp + :language: c++ + :start-after: [sepimage-2-beg] + :end-before: [sepimage-2-end] + :dedent: 4 + +In this example, the constraint on the initial set can be seen as a distance constraint. It can be treated as an inversion as follows : + +.. tabs:: + + .. group-tab:: Python + + .. literalinclude:: src.py + :language: py + :start-after: [sepimage-3-beg] + :end-before: [sepimage-3-end] + :dedent: 4 + + .. group-tab:: C++ + + .. literalinclude:: src.cpp + :language: c++ + :start-after: [sepimage-3-beg] + :end-before: [sepimage-3-end] + :dedent: 4 + +The separator can then be constructed and used, for example with a paver to get both an inner and an outer approximation of the image set. +The resulting paving is showed in the next figure. + +.. tabs:: + + .. group-tab:: Python + + .. literalinclude:: src.py + :language: py + :start-after: [sepimage-4-beg] + :end-before: [sepimage-4-end] + :dedent: 4 + + .. group-tab:: C++ + + .. literalinclude:: src.cpp + :language: c++ + :start-after: [sepimage-4-beg] + :end-before: [sepimage-4-end] + :dedent: 4 + +.. figure:: ./sep_image.png + :width: 400px + :align: center + +With a lower resolution, the performances of the separator improve as shown in the following figure. + +.. tabs:: + + .. group-tab:: Python + + .. literalinclude:: src.py + :language: py + :start-after: [sepimage-5-beg] + :end-before: [sepimage-5-end] + :dedent: 4 + + .. group-tab:: C++ + + .. literalinclude:: src.cpp + :language: c++ + :start-after: [sepimage-5-beg] + :end-before: [sepimage-5-end] + :dedent: 4 + +.. figure:: ./sep_image_fine.png + :width: 400px + :align: center + +A more complex example is available on the public github repository. It treats the topic of the explored area, +which is a classical problem in robotics. + +* `Python version `_ +* `C++ version `_ diff --git a/doc/manual/manual/contractors/analytic/src.cpp b/doc/manual/manual/contractors/analytic/src.cpp index 2b8438005..adbffa9de 100644 --- a/doc/manual/manual/contractors/analytic/src.cpp +++ b/doc/manual/manual/contractors/analytic/src.cpp @@ -12,6 +12,7 @@ #include #include #include +#include #include #include @@ -86,4 +87,40 @@ TEST_CASE("CtcInverse - manual") CHECK(c.fnc().input_size() == 2); CHECK(c.fnc().output_size() == 1); } +} + +TEST_CASE("SepImage - manual") +{ + { + // [sepimage-1-beg] + VectorVar y (2); + AnalyticFunction f ({y},{3.*(y[0]+1),y[1]+0.5*sin(3.*y[0])}); + // [sepimage-1-end] + + // [sepimage-2-beg] + // {psi0,Sigma} is a gnomonic atlas of the unit circle + VectorVar X(1); + AnalyticFunction psi0 ({X},{cos(X[0]*PI/2.),sin(X[0]*PI/2.)}); + + OctaSym id ({1, 2}); + OctaSym s ({-1, -2}); + + vector Sigma ({id,s}); + // [sepimage-2-end] + + // [sepimage-3-beg] + AnalyticFunction h ({y},sqrt(sqr(y[0])+sqr(y[1]))); + CtcInverse ctc_in (h,Interval(0,1)); + // [sepimage-3-end] + + // [sepimage-4-beg] + SepImage sep1 (f,psi0,Sigma,0.125,ctc_in); + DefaultFigure::pave({{-0.5,6.5},{-1.5,1.5}},sep1,0.05); + // [sepimage-4-end] + + // [sepimage-5-beg] + SepImage sep2 (f,psi0,Sigma,0.0625,ctc_in); + DefaultFigure::pave({{-0.5,6.5},{-1.5,1.5}},sep2,0.05); + // [sepimage-5-end] + } } \ No newline at end of file diff --git a/doc/manual/manual/contractors/analytic/src.py b/doc/manual/manual/contractors/analytic/src.py index 60505c4dd..2631ba291 100644 --- a/doc/manual/manual/contractors/analytic/src.py +++ b/doc/manual/manual/contractors/analytic/src.py @@ -106,5 +106,38 @@ def tests_CtcInverse_manual(test): assert c.fnc().output_size() == 1 # [ctcinv-7-end] +def tests_SepImage_manual(test): + + # [sepimage-1-beg] + y = VectorVar(2) + f = AnalyticFunction([y], [3.*(y[0]+1),y[1]+0.5*sin(3.*y[0])]) + # [sepimage-1-end] + + # [sepimage-2-beg] + # {psi0,Sigma} is a gnomonic atlas of the box [-1,1]^2 + X = VectorVar(1) + psi0 = AnalyticFunction([X],[cos(X[0]*PI/2.),sin(X[0]*PI/2.)]) + + id = OctaSym([1,2]) + s = OctaSym([-1,-2]) + + Sigma = [id,s] + # [sepimage-2-end] + + # [sepimage-3-beg] + h = AnalyticFunction([y],sqrt(sqr(y[0])+sqr(y[1]))) + ctc_in = CtcInverse(h,Interval(0,1)) + # [sepimage-3-end] + + # [sepimage-4-beg] + sep = SepImage(f,psi0,Sigma,0.125,ctc_in) + DefaultFigure.pave([[-0.5,6.5],[-1.5,1.5]],sep,0.05) + # [sepimage-4-end] + + # [sepimage-5-beg] + sep = SepImage(f,psi0,Sigma,0.0625,ctc_in) + DefaultFigure.pave([[-0.5,6.5],[-1.5,1.5]],sep,0.05) + # [sepimage-5-end] + if __name__ == '__main__': unittest.main() \ No newline at end of file diff --git a/doc/manual/manual/contractors/index.rst b/doc/manual/manual/contractors/index.rst index 9b2537821..8e6eea3c8 100644 --- a/doc/manual/manual/contractors/index.rst +++ b/doc/manual/manual/contractors/index.rst @@ -5,6 +5,7 @@ Contractors, separators .. toctree:: CtcInter + SepImage CtcInverse CtcLohner CtcDist @@ -77,6 +78,9 @@ Overview of contractors and separators * - ``CtcAction`` - ``SepAction`` + * - — + - :ref:`SepImage ` + * - :ref:`CtcInverse ` - ``SepInverse`` diff --git a/doc/manual/manual/extensions/capd/peibos_capd.rst b/doc/manual/manual/extensions/capd/peibos_capd.rst index 3fae692b3..453a7303f 100644 --- a/doc/manual/manual/extensions/capd/peibos_capd.rst +++ b/doc/manual/manual/extensions/capd/peibos_capd.rst @@ -15,7 +15,7 @@ Considering a dynamical system :math:`\dot{\mathbf{x}}=\gamma(\mathbf{x})`, the Gnomonic atlas -------------- -To handle the boundary of the initial set :math:`\mathbb{X}_0`, the PEIBOS tool relies on a gnomonic atlas. See :ref:`subsec-functions-peibos-gnomonic-atals`. +To handle the boundary of the initial set :math:`\mathbb{X}_0`, the PEIBOS tool relies on a gnomonic atlas. See :ref:`subsec-functions-peibos-gnomonic-atlas`. Use --- diff --git a/doc/manual/manual/functions/parallelepiped/parallelepiped_eval.rst b/doc/manual/manual/functions/parallelepiped/parallelepiped_eval.rst index 92a58eb57..f56f29a40 100644 --- a/doc/manual/manual/functions/parallelepiped/parallelepiped_eval.rst +++ b/doc/manual/manual/functions/parallelepiped/parallelepiped_eval.rst @@ -9,7 +9,7 @@ Parallelepiped inclusion function Use case -------- -Consider a function :math:`\mathbf{f}:\mathbb{R}^n\to\mathbb{R}^m`. In the case where :math:`0 + +using namespace std; +using namespace codac2; + +int main() +{ + set_nb_threads(max_threads()); + + // {psi0,Sigma} is a gnomonic atlas of the box [-1,1]^2 + VectorVar X(1); + AnalyticFunction psi0 ({X},{X[0],1}); + + OctaSym id ({1, 2}); + OctaSym s ({-2, 1}); + + vector Sigma ({id,s,s*s,s.invert()}); + + // In this box, we consider that the x-axis is the width of the linear sensor + // and the y-axis is the time + VectorVar y (2); + auto L = 0.1*y[0]; + auto t = 1.2*y[1]; + + // We construct the trajectory of the robot + auto traj = vec(pow(t,3)-t, 1-sqr(t)); + // We need its derivative to compute the orthogonal to the trajectory (for the sensor) + auto dtraj = vec(3*sqr(t)-1, -2*t); + auto dtraj_norm = sqrt(sqr(dtraj[0])+sqr(dtraj[1])); + + // The image of the box [-1,1]^2 by f is the swept area + AnalyticFunction f ({y},{traj[0]-L*dtraj[1]/dtraj_norm, traj[1]+L*dtraj[0]/dtraj_norm}); + + // For the SepImage, we need a contractor on the initial set (here a simple box) + IntervalVector X0 = IntervalVector::constant(2,{-1,1}); + CtcWrapper ctc_in (X0); + + // Separator on the area seen by a robot + SepImage sep (f,psi0,Sigma,0.01,ctc_in); + + // Visualizing the separator + IntervalVector Y0 ({{-0.8,0.8},{-0.7,1.3}}); + + DefaultFigure::pave(Y0,sep,0.01); +} \ No newline at end of file diff --git a/examples/04_explored_area/main_peibos.py b/examples/04_explored_area/main_peibos.py new file mode 100644 index 000000000..d5b55bd57 --- /dev/null +++ b/examples/04_explored_area/main_peibos.py @@ -0,0 +1,42 @@ +from codac import * +import numpy as np + +if __name__=="__main__": + set_nb_threads(max_threads()) + + # {psi0,Sigma} is a gnomonic atlas of the box [-1,1]^2 + X = VectorVar(1) + psi0 = AnalyticFunction([X],[X[0],1]) + + id = OctaSym([1,2]) + s = OctaSym([-2,1]) + + Sigma = [id,s,s*s,s.invert()] + + # In this box, we consider that the x-axis is the width of the linear sensor + # and the y-axis is the time + y = VectorVar(2) + L = 0.1*y[0] + t = 1.2*y[1] + + # We construct the trajectory of the robot + traj = vec(pow(t,3)-t,1-sqr(t)) + # We need its derivative to compute the orthogonal to the trajectory (for the sensor) + dtraj = vec(3*sqr(t)-1, -2*t) + dtraj_norm = sqrt(sqr(dtraj[0])+sqr(dtraj[1])) + + # The image of the box [-1,1]^2 by f is the swept area + f = AnalyticFunction([y], [traj[0]-L*dtraj[1]/dtraj_norm, traj[1]+L*dtraj[0]/dtraj_norm]) + + # For the SepImage, we need a contractor on the initial set (here a simple box) + X0 = IntervalVector.constant(2,[-1,1]) + ctc_in = CtcWrapper(X0) + + # Separator on the area seen by a robot + sep = SepImage(f,psi0,Sigma,0.01,ctc_in) + + # Visualizing the separator + Y0 = IntervalVector([[-0.8,0.8],[-0.7,1.3]]) + + DefaultFigure.pave(Y0,sep,0.01) + diff --git a/python/src/core/CMakeLists.txt b/python/src/core/CMakeLists.txt index 25221ecb7..c8a7c7a43 100644 --- a/python/src/core/CMakeLists.txt +++ b/python/src/core/CMakeLists.txt @@ -107,6 +107,7 @@ separators/codac2_py_SepCross.cpp separators/codac2_py_SepCtcBoundary.cpp separators/codac2_py_SepCtcPair.cpp + separators/codac2_py_SepImage.cpp separators/codac2_py_SepInter.cpp separators/codac2_py_SepInverse.cpp separators/codac2_py_SepNot.cpp diff --git a/python/src/core/codac2_py_core.cpp b/python/src/core/codac2_py_core.cpp index 3b0f73931..f0cefd778 100644 --- a/python/src/core/codac2_py_core.cpp +++ b/python/src/core/codac2_py_core.cpp @@ -145,6 +145,7 @@ void export_SepChi(py::module& m, py::class_& pysep); void export_SepCross(py::module& m, py::class_& pysep); void export_SepCtcBoundary(py::module& m, py::class_& pysep); void export_SepCtcPair(py::module& m, py::class_& pysep); +void export_SepImage(py::module& m, py::class_& pysep); void export_SepInter(py::module& m, py::class_& sep); void export_SepInverse(py::module& m, py::class_& sep); void export_SepNot(py::module& m, py::class_& sep); @@ -330,6 +331,7 @@ PYBIND11_MODULE(_core, m) export_SepCross(m,py_sep); export_SepCtcBoundary(m,py_sep); export_SepCtcPair(m,py_sep); + export_SepImage(m,py_sep); export_SepInter(m,py_sep); export_SepInverse(m,py_sep); export_SepNot(m,py_sep); diff --git a/python/src/core/separators/codac2_py_SepImage.cpp b/python/src/core/separators/codac2_py_SepImage.cpp new file mode 100644 index 000000000..62b3ca5d9 --- /dev/null +++ b/python/src/core/separators/codac2_py_SepImage.cpp @@ -0,0 +1,42 @@ +/** + * Codac binding (core) + * ---------------------------------------------------------------------------- + * \date 2026 + * \author Maël Godard + * \copyright Copyright 2026 Codac Team + * \license GNU Lesser General Public License (LGPL) + */ + +#include +#include +#include +#include +#include +#include +#include "codac2_py_Sep.h" +#include "codac2_py_SepImage_docs.h" // Generated file from Doxygen XML (doxygen2docstring.py): + +using namespace std; +using namespace codac2; +namespace py = pybind11; +using namespace pybind11::literals; + +void export_SepImage(py::module& m, py::class_& pysep) +{ + py::class_ exported(m, "SepImage", pysep, SEPIMAGE_MAIN); + exported + + .def(py::init( + [](const AnalyticFunction& f, const AnalyticFunction& psi_0, const std::vector& Sigma, double epsilon, const CtcBase& ctc_init, double epsilon_pave) + { + return std::make_unique(f,psi_0, Sigma, epsilon, ctc_init.copy(), epsilon_pave); + }), + SEPIMAGE_SEPIMAGE_CONST_ANALYTICFUNCTION_VECTORTYPE_REF_CONST_ANALYTICFUNCTION_VECTORTYPE_REF_CONST_VECTOR_OCTASYM_REF_DOUBLE_CONST_C_REF_DOUBLE, + "f"_a, "psi_0"_a, "Sigma"_a, "epsilon"_a, "ctc_init"_a, "epsilon_pave"_a = -1) + + .def("separate", &SepImage::separate, + BOXPAIR_SEPIMAGE_SEPARATE_CONST_INTERVALVECTOR_REF_CONST, + "x"_a) + + ; +} \ No newline at end of file diff --git a/src/core/CMakeLists.txt b/src/core/CMakeLists.txt index 4ad18aaa0..1573630ce 100644 --- a/src/core/CMakeLists.txt +++ b/src/core/CMakeLists.txt @@ -237,6 +237,8 @@ ${CMAKE_CURRENT_SOURCE_DIR}/separators/codac2_SepCtcPair.h ${CMAKE_CURRENT_SOURCE_DIR}/separators/codac2_SepEllipse.cpp ${CMAKE_CURRENT_SOURCE_DIR}/separators/codac2_SepEllipse.h + ${CMAKE_CURRENT_SOURCE_DIR}/separators/codac2_SepImage.cpp + ${CMAKE_CURRENT_SOURCE_DIR}/separators/codac2_SepImage.h ${CMAKE_CURRENT_SOURCE_DIR}/separators/codac2_SepInter.cpp ${CMAKE_CURRENT_SOURCE_DIR}/separators/codac2_SepInter.h ${CMAKE_CURRENT_SOURCE_DIR}/separators/codac2_SepInverse.h diff --git a/src/core/functions/analytic/codac2_AnalyticFunction.h b/src/core/functions/analytic/codac2_AnalyticFunction.h index 0fe6fd645..a8928b619 100644 --- a/src/core/functions/analytic/codac2_AnalyticFunction.h +++ b/src/core/functions/analytic/codac2_AnalyticFunction.h @@ -214,7 +214,7 @@ namespace codac2 Parallelepiped parallelepiped_eval(const Args&... x) const { this->check_valid_inputs(x...); - assert_release(this->input_size() < this->output_size() && + assert_release(this->input_size() <= this->output_size() && "Parallelepiped evaluation requires more outputs than inputs."); assert_release(this->input_size() > 0 && "Parallelepiped evaluation requires at least one input."); diff --git a/src/core/separators/codac2_SepImage.cpp b/src/core/separators/codac2_SepImage.cpp new file mode 100644 index 000000000..6a5da8153 --- /dev/null +++ b/src/core/separators/codac2_SepImage.cpp @@ -0,0 +1,18 @@ +/** + * codac2_SepImage.cpp + * ---------------------------------------------------------------------------- + * \date 2026 + * \author Maël Godard + * \copyright Copyright 2024 Codac Team + * \license GNU Lesser General Public License (LGPL) + */ + +#include "codac2_SepImage.h" +#include "codac2_CtcPolar.h" + +using namespace codac2; + +BoxPair SepImage::separate(const IntervalVector& x) const +{ + return _sep.front()->separate(x); +} \ No newline at end of file diff --git a/src/core/separators/codac2_SepImage.h b/src/core/separators/codac2_SepImage.h new file mode 100644 index 000000000..dd8ce939d --- /dev/null +++ b/src/core/separators/codac2_SepImage.h @@ -0,0 +1,232 @@ +/** + * \file codac2_SepImage.h + * ---------------------------------------------------------------------------- + * \date 2026 + * \author Maël Godard + * \copyright Copyright 2025 Codac Team + * \license GNU Lesser General Public License (LGPL) + */ + +#pragma once + +#include "codac2_Sep.h" +#include "codac2_Collection.h" +#include "codac2_Paving.h" +#include "codac2_PavingNode.h" +#include "codac2_Subpaving.h" +#include "codac2_AnalyticFunction.h" +#include "codac2_OctaSym.h" +#include "codac2_peibos.h" +#include "codac2_CtcUnion.h" +#include "codac2_CtcWrapper.h" +#include "codac2_CtcIdentity.h" +#include "codac2_SepCtcPair.h" +#include "codac2_SepInverse.h" +#include "codac2_SepWrapper.h" +#include "codac2_SepInter.h" +#include "codac2_pave.h" +#include + +namespace codac2 +{ + /** + * \brief A separator on the direct image of a set by a function. It relies on a boundary approach + */ + class SepImage : public Sep + { + public: + + /** + * \brief Constructor of the separator + * + * \param f The analytic function \f$\mathbf{f}:\mathbb{R}^n\rightarrow\mathbb{R}^p,p\geq n\f$ + * \param psi_0 The transformation function \f$\psi_0:\mathbb{R}^m\rightarrow\mathbb{R}^n\f$ to construct the atlas + * \param Sigma The set of symmetry operators \f$\sigma\f$ to construct the atlas + * \param epsilon The maximum diameter of the boxes to split \f$[-1,1]^m\f$ before computing the parallelepiped inclusions + * \param ctc_init The contractor on the initial set + * \param epsilon_pave An optionnal parameter to custom the resolution of the paving step + */ + template + requires IsCtcBaseOrPtr + SepImage(const AnalyticFunction& f, const AnalyticFunction& psi_0, const std::vector& Sigma, double epsilon, const C& ctc_init, double epsilon_pave = -1) + :Sep(f.output_size()) + { + // Using PEIBOS to get the boundary of the image set + int dim = f.output_size(); + auto v_par = PEIBOS(f,psi_0,Sigma,epsilon); + + IntervalVector X0 (dim); + ctc_init.contract(X0); + + // Contractor on the boundary + CtcUnion ctc_union (dim); + for (const auto& par : v_par) + ctc_union |= CtcWrapper(par); + + SepCtcPair sep_boundary(CtcIdentity(dim), ctc_union); + + IntervalVector Y0 (dim); + ctc_union.contract(Y0); + Y0.inflate(Y0.min_diam()/10.); + + // paving resolution, can be tuned + if (epsilon_pave == -1) + epsilon_pave = Y0.min_diam()/50.; + + // Initial paving, with only an outer aproximation + PavingInOut p = pave(Y0,sep_boundary,epsilon_pave); + + // Determination of the inner connected susets + auto v_cs = p.connected_subsets(PavingInOut::outer_complem); + std::list> cs_to_color; + for (const auto& cs: v_cs) + { + IntervalVector max_box = IntervalVector::empty(dim); + for (const auto& box : cs.boxes()) + if (box.volume()>max_box.volume()) + max_box=box; + + CtcInverse ctc_inv (f, max_box.mid()); + IntervalVector X0_copy (X0); + bool could_contract = false; + + ctc_inv.contract(X0_copy); + ctc_init.contract(X0_copy); + + if (X0_copy != X0) + { + could_contract = true; + if (!X0_copy.is_empty()) + cs_to_color.push_back(std::make_shared(cs)); + } + + // if the CtcInverse was not able to conclude directly, we subdivide the initial box + + if (!could_contract) + { + std::vector boxes; + split(X0_copy,0.125,boxes); + for (auto box: boxes) + { + IntervalVector box_copy (box); + ctc_inv.contract(box); + ctc_init.contract(box); + if (box != box_copy) + { + could_contract = true; + if (!box.is_empty()) + { + cs_to_color.push_back(std::make_shared(cs)); + break; + } + } + } + } + + if (!could_contract) + printf("warning, the SepImage was not able to conclude on the inner regions\n"); + } + + // Turning the previously identified connected subsets as inner (displayed green) + auto visitor = [&](std::shared_ptr> n) + { + std::list boxes_to_color; + for (const auto& cs : cs_to_color) + for (const auto& box : cs->boxes()) + boxes_to_color.push_back(box); + + IntervalVector h = n->hull(); + IntervalVector u = n->unknown(); + + for (const auto& bi : boxes_to_color) + { + IntervalVector bih = bi & h; + if (!bih.is_empty()) + { + if (std::get<0>(n->boxes()).is_empty()) + { + n->set_boxes(std::make_tuple(std::get<1>(n->boxes()),std::get<0>(n->boxes()))); + return false; + } + + std::list lbi; + + for (const auto& bj : boxes_to_color) + { + IntervalVector bjh = bj & h; + if (!bjh.is_empty()) + lbi.push_back(bjh); + } + + IntervalVector prev_x_in = IntervalVector::empty(dim); + + while (h != prev_x_in) + { + prev_x_in = IntervalVector(h); + + for (const auto& li : lbi) + { + std::list d = h.diff(li); + + if (!d.empty()) + { + + IntervalVector hi = IntervalVector::empty(dim); + for (const auto& di : d) + hi |= (h & di); + h &= hi; + + } + } + } + + n->set_boxes(std::make_tuple(hull(lbi)|u, h)); + return true; + } + } + return true; + }; + + // Use the visitor to color the inner regions + p.tree()->visit(visitor); + + // Create a separator on this paving + SepWrapper sep_paving(p); + + // By construction, we know that we have the whole image set. + // What is outside of the paving box Y0 is for sure ouside of it. + SepInter sep_inter = sep_paving&Y0; + + _sep.push_back(sep_inter); + } + + /** + * \brief Constructor of the separator + * + * \param f The analytic function \f$\mathbf{f}:\mathbb{R}^n\rightarrow\mathbb{R}^p,p\geq n\f$ + * \param psi_0 The transformation function \f$\psi_0:\mathbb{R}^m\rightarrow\mathbb{R}^n\f$ to construct the atlas + * \param Sigma The set of symmetry operators \f$\sigma\f$ to construct the atlas + * \param epsilon The maximum diameter of the boxes to split \f$[-1,1]^m\f$ before computing the parallelepiped inclusions + * \param ctc_init A shared pointer to the contractor on the initial set + * \param epsilon_pave An optionnal parameter to custom the resolution of the paving step + */ + template + requires IsCtcBaseOrPtr + SepImage(const AnalyticFunction& f, const AnalyticFunction& psi_0, const std::vector& Sigma, double epsilon, const std::shared_ptr& ctc_init, double epsilon_pave = -1) + : SepImage(f, psi_0, Sigma, epsilon, *ctc_init, epsilon_pave) + { } + + /** + * \brief Separates the box. + * + * \param x The box to separate. + * + * \return A box pair according to separator arithmetics. + */ + BoxPair separate(const IntervalVector& x) const; + + protected: + + Collection _sep; + }; +} \ No newline at end of file diff --git a/src/core/separators/codac2_SepTest.h b/src/core/separators/codac2_SepTest.h index bae989c9e..48fb6cc71 100644 --- a/src/core/separators/codac2_SepTest.h +++ b/src/core/separators/codac2_SepTest.h @@ -39,6 +39,8 @@ namespace codac2 * \brief Separates the box. Both the inner and outer are either empty or the box itself. * * \param x The box to separate. + * + * \return A box pair according to separator arithmetics. */ BoxPair separate(const IntervalVector& x) const;