With peibos by SimonRohou · Pull Request #250 · codac-team/codac

SimonRohou · 2025-05-15T07:30:01Z

PR juste pour pouvoir ajouter des commentaires, l'interface ne semble pas le permettre autrement.

…ibos

SimonRohou · 2025-05-15T07:30:24Z

@@ -0,0 +1,85 @@
+// Example from the publication:
+// https://www.ensta-bretagne.fr/jaulin/paper_centeredActa.pdf


À retirer / changer

SimonRohou · 2025-05-15T07:31:18Z

+  VectorVar X_2d(1);
+  AnalyticFunction psi0_2d ({X_2d},{cos(X_2d[0]*PI/4.-PI/2),sin(X_2d[0]*PI/4.-PI/2)});
+
+  vector<vector<int>> generators_2d ({{1,2},


Les générators sont-ils une matrice d'entiers ?

On peut le voir comme ça, mais je préfère le voir comme un vecteur dont chaque élément est un vector qui sert à créer une OctaSym. A partir des symmétries définies par l'utilisateur je "complète" avec la fonction generate_symmetries.
Dans les réflexions en cours je me demande si il ne vaut pas mieux demander à l'utilisateur de fournir directement la liste de symmétries pour qu'il puisse faire ce qu'il veut, quitte à le laisser appeler la fonction generate_symmetries par lui-même... à discuter

SimonRohou · 2025-05-15T07:32:03Z

+
+  Figure2D figure_2d ("Henon Map", GraphicOutput::VIBES);
+  figure_2d.set_window_properties({25,50},{500,500});
+  figure_2d.set_axes(axis(0,{-1.4,2.2}), axis(1,{-0.4,0.3}));


D'un point de vue cosmétique il est possible d'écrire set_axes({0,{-1.4,2.2}}, {1,{-0.4,0.3}})

SimonRohou · 2025-05-15T07:32:58Z

+
+  // 3D example of the PEIBOS algorithm
+  VectorVar y_3d(3);
+  AnalyticFunction f_3d({y_3d},{sqr(y_3d[0])-sqr(y_3d[1])+y_3d[0],2*y_3d[0]*y_3d[1]+y_3d[1],y_3d[2]});


Pour info, sqr(y_3d[0]) peut désormais aussi s'écrire y_3d[0]^2 (et également y_3d[0]**2 en Python).

SimonRohou · 2025-05-15T07:34:29Z

+  for p in v_par_nd:
+    for vertice in p.vertices():
+      figure_2d_nd_xy.draw_point(vertice, [Color.red(), Color.red(0.5)])
+      figure_2d_nd_zy.draw_point(vertice, [Color.red(), Color.red(0.5)])


Duplicate ?

Pas tout à fait, ce sont deux "vues 2D", chacun représente un plan (respectivement XY et ZY):
figure_2d_nd_xy.set_axes(axis(0,{-1.2,1.2}), axis(1,{-1.2,1.2}));
figure_2d_nd_zy.set_axes(axis(2,{-1.2,1.2}), axis(1,{-1.2,1.2}));

SimonRohou · 2025-05-15T07:34:46Z

@@ -0,0 +1,61 @@
+// Example from the publication:
+// https://www.ensta-bretagne.fr/jaulin/paper_centeredActa.pdf


SimonRohou · 2025-05-15T07:38:50Z

+  .def("vertices", &Parallelepiped::vertices,
+    VECTOR_VECTOR_PARALLELEPIPED_VERTICES_CONST)
+
+  .def_readwrite("z", &Parallelepiped::z)


La doc des attributs peut aussi être ajoutée

SimonRohou · 2025-05-15T07:39:22Z

+void export_peibos(py::module& m)
+{
+  m.def("PEIBOS", 
+    [](const py::object& f, const py::object& psi_0, vector<vector<int>> generators, double epsilon)


vector<vector<int>> en const reference

SimonRohou · 2025-05-15T07:40:03Z

+  m.def("PEIBOS", 
+    [](const py::object& f, const py::object& psi_0, vector<vector<int>> generators, double epsilon)
+    {
+      return PEIBOS(cast<AnalyticFunction<VectorType>>(f), cast<AnalyticFunction<VectorType>>(psi_0), generators, epsilon);


Si les générateurs sont des indices de coefficients, il faudrait les shifter pour coller avec les conventions Matlab (premier indice à 1).

chaque générateur est utilisé tel quel (i.e. un vecteur d'entier) pour créer une OctaSym, donc pas de problème d'indice

SimonRohou · 2025-05-15T07:40:50Z

+      return PEIBOS(cast<AnalyticFunction<VectorType>>(f), cast<AnalyticFunction<VectorType>>(psi_0), generators, epsilon, offset);
+    },
+    VECTOR_PARALLELEPIPED_PEIBOS_CONST_ANALYTICFUNCTION_VECTORTYPE_REF_CONST_ANALYTICFUNCTION_VECTORTYPE_REF_CONST_VECTOR_VECTOR_INT_REF_DOUBLE_CONST_VECTOR_REF,
+    "f"_a, "psi_0"_a, "generators"_a, "epsilon"_a, "offset"_a);


Valeur par défaut à Vector::zero(..) pour offset ?

pas vraiment possible, d'où le doublon avec une fonction avec et une sans offset : la taille du vecteur d'offset dépend de la taille de sortie de la fonction psi_0. De fait comme je n'ai pas le droit d'écrire
Vector::Zero(psi_0.output_size())
comme valeur par défaut (puisque psi_0 est lui même un argument) je n'ai pas trouvé d'autre solution que d'avoir deux fonctions... peut-être qu'il existe une astuce pour contourner le problème ?

SimonRohou · 2025-05-15T07:41:53Z

+Parallelepiped::Parallelepiped(const Vector& z_, const Matrix& A_)
+  : z(z_), A(A_)
+{
+  assert(z.size() == A.rows());


Je l'aurais plutôt mis en assert_release

SimonRohou · 2025-05-15T07:42:14Z

+  Vector v6 = z + A.col(0) + A.col(1) - A.col(2);
+  Vector v7 = z + A.col(0) - A.col(1) - A.col(2);
+  Vector v8 = z - A.col(0) - A.col(1) - A.col(2);
+  vertices.push_back(v1);


vertices.push_back(z + A.col(0) + A.col(1) + A.col(2));

SimonRohou · 2025-05-15T07:42:50Z

+
+vector<Vector> Parallelepiped::vertices() const
+{
+  assert(A.rows() == 3);


Je l'aurais aussi mis en assert_release pour donner des indications à l'utilisateur en cas de mauvais usage.

SimonRohou · 2025-05-15T07:44:37Z

+    return (a.cross(b)).norm()/((b-a).norm());
+  }
+
+  bool contains (const vector<OctaSym>& symmetries, const OctaSym& symmetry, const AnalyticFunction<VectorType>& psi_0)


Peut-être des noms d'arguments plus différenciés ? symmetries, symmetry

SimonRohou · 2025-05-15T07:45:10Z

+    IntervalVector test_box =  Interval(-1.,1.)*IntervalVector::Ones(psi_0.input_size());
+    IntervalVector psi_0_x = psi_0.eval(EvalMode::NATURAL,test_box);
+
+    for (OctaSym s : symmetries)


const OctaSym& s

SimonRohou · 2025-05-15T07:45:48Z

+  bool contains (const vector<OctaSym>& symmetries, const OctaSym& symmetry, const AnalyticFunction<VectorType>& psi_0)
+  {
+    IntervalVector test_box =  Interval(-1.,1.)*IntervalVector::Ones(psi_0.input_size());
+    IntervalVector psi_0_x = psi_0.eval(EvalMode::NATURAL,test_box);


Pourquoi restreindre à EvalMode::NATURAL ? Tu as peut-être perçu un bug avec le mode d'évaluation par défaut ?

J'ai "peur" que l'utilisation de la forme centrée puisse donner deux résultats très légèrement différents pour deux symmétries qui devraient être les mêmes (par exemple tourner à gauche de 90° deux fois ou tourner à droite de 90° deux fois). De plus ici je ne cherche pas un "bon résultat", juste à faire un test d'unicité, donc l'évaluation naturelle est suffisante.

SimonRohou · 2025-05-15T07:47:21Z

+    return false;
+  }
+
+  vector<OctaSym> generate_symmetries (const vector<vector<int>>& generators, const AnalyticFunction<VectorType>& psi_0)


Peut-être intéressant de remplacer const vector<vector<int>> par une Matrix voire une Matrix<int>.

Voir mon commentaire précédent sur ce que représente ce vecteur de vecteur, mais dans l'absolu on pourrait effectivement utiliser une Matrix.

SimonRohou · 2025-05-15T07:49:26Z

+    for (int i = 0; i < ((int) generators.size()); i++)
+    {
+      OctaSym symmetry = OctaSym(generators[i]);
+      symmetries.push_back(symmetry);
+    }


Juste pour info :

for(const auto& g : generators) symmetries.push_back(OctaSym(g)); // voire même symmetries.push_back({ g })

SimonRohou · 2025-05-15T07:51:32Z

+    auto xc = X.mid();
+
+    IntervalVector dX=X-xc;
+    IntervalMatrix JJg_punc=JJf_punc*IntervalMatrix(symmetry.permutation_matrix())*psi_0.diff(xc);


Pour info, tu peux écrire symmetry.permutation_matrix().template cast<Interval>() au lieu de IntervalMatrix(symmetry.permutation_matrix()). C'est probablement plus efficace pour Eigen (mais rien de fondamental).

SimonRohou · 2025-05-15T07:54:59Z

+    auto xc = X.mid();
+
+    IntervalVector dX=X-xc;
+    IntervalMatrix JJg_punc=JJf_punc*IntervalMatrix(symmetry.permutation_matrix())*psi_0.diff(xc);
+
+    IntervalMatrix JJg=JJf*IntervalMatrix(symmetry.permutation_matrix())*psi_0.diff(X);
+
+    IntervalVector E = (JJg - JJg_punc)*dX;
+    Interval N = sqr(E[0]) + sqr(E[1]);
+
+    return std::sqrt(N.ub());


Tu peux aussi choisir de ne pas typer explicitement dX, JJg_punc, JJg, E, N, avec auto, pour laisser Eigen faire ses optimisations par expressions templates. https://eigen.tuxfamily.org/index.php?title=Expression_templates

SimonRohou · 2025-05-15T07:56:12Z

+    }
+  }
+
+  Matrix inflate_flat_parallelepiped(const Matrix& Jz, double epsilon, double rho)


Fonctionne uniquement pour les flat parallelepipeds ?

C'est une appellation un peu abusive effectivement. Au départ cette fonction prenait un parallelepipède 2D et le "gonflait" en 3D d'une valeur rho. Maintenant on part d'un parallélépipède m-D et on le gonfle en n-D (m<n), il faudrait peut-être changer le nom de la fonction pour être plus explicite.

SimonRohou · 2025-05-15T07:57:32Z

+  vector<Parallelepiped> PEIBOS(const AnalyticFunction<VectorType>& f, const AnalyticFunction<VectorType>& psi_0, const vector<vector<int>>& generators , double epsilon)
+  {
+    return PEIBOS(f, psi_0, generators, epsilon, Vector::Zero(psi_0.output_size()));
+  }


Tu pourrais ne faire qu'une fonction avec une valeur par défaut pour l'argument offset.

Cf mon commentaire sur le binding python de cette même fonction

…ibos

godardma added 8 commits May 7, 2025 23:54

first version of peibos in codac

5be5460

[peibos] corrected python binding

cf75011

[peibos] n-D implementation

6171810

[unsupported] update for non_capd compilation

afc5fbb

[unsupported] update for non_capd compilation

2fbfa8c

[peibos] the const-correctness commit

931db23

[peibos] n-D implementation for CAPD

ec5c2ca

Merge branch 'codac2_dev' of github.com:codac-team/codac into with_pe…

ad4edda

…ibos

SimonRohou commented May 15, 2025

View reviewed changes

godardma added 29 commits June 30, 2025 10:00

minor update

d756ed1

Merge branch 'codac2_dev' of github.com:codac-team/codac into with_pe…

2c80aa6

…ibos

[peibos] minor modifications

3ca6a08

[peibos] fix of uninitialized matrix

a7d6ac0

[peibos] updated example

04374bd

[peibos] updated example

3cbdaa4

Merge branch 'codac2_dev' of github.com:godardma/codac into with_peibos

b4b8c00

Merge branch 'codac2_dev' of github.com:godardma/codac into with_peibos

5416a13

Merge branch 'codac2_dev' of github.com:codac-team/codac into with_pe…

56cb26a

…ibos

[peibos] updated example

535ce69

[peibos] Zonotopes base and Parallelepiped projection

14898c0

[peibos] fixed ill-formatted assert

f603b06

[python] added attribute doc

c90b84e

[peibos] added full example in matlab

a382978

[peibos_capd] multiple timesteps allowed

448c94a

[peibos_capd] map for multiple timesteps implementation

85c32f8

[peibos_capd] Implementation in python + example

9401837

[actions] updated capd setup

ae7a71c

[actions] updated capd setup

480ab6c

[peibos_capd] Example in matlab

9bcf9ca

[peibos_capd] minor update

0e9dd4c

[peibos_capd] minimal capd include

360dd7d

[peibos_capd] fix

f937d95

[peibos_capd] fix

54fbd47

Merge branch 'codac2_dev' of github.com:codac-team/codac into with_pe…

51560a5

…ibos

[peibos] removed automatic generate_symmetries

bfafd69

[peibos] removed automatic generate_symmetries

0672375

Merge branch 'codac2_dev' of github.com:codac-team/codac into with_pe…

f795a07

…ibos

minor update

e5008a3

SimonRohou closed this Oct 21, 2025

		@@ -0,0 +1,85 @@
		// Example from the publication:
		// https://www.ensta-bretagne.fr/jaulin/paper_centeredActa.pdf

		@@ -0,0 +1,61 @@
		// Example from the publication:
		// https://www.ensta-bretagne.fr/jaulin/paper_centeredActa.pdf

Uh oh!

Conversation

SimonRohou commented May 15, 2025

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

2 participants