Merge pull request #313 from godardma/adding_peibos

[functions] final patch for parallelepiped_eval

Author: SimonRohou

src/core/domains/zonotope/codac2_Parallelepiped.cpp

@@ -10,7 +10,6 @@
 #include "codac2_Parallelepiped.h"
 #include "codac2_inversion.h"
 
-using namespace std;
 using namespace codac2;
 
 Parallelepiped::Parallelepiped(const Vector& z_, const Matrix& A_)
@@ -19,7 +18,7 @@ Parallelepiped::Parallelepiped(const Vector& z_, const Matrix& A_)
   assert_release(A.cols() <= z.size() && "too many vectors, you are describing a zonotope");
 }
 
-void generate_vertices(Index i, Index n, const Vector& z, const Matrix& A, vector<Vector>& L_v)
+void generate_vertices(Index i, Index n, const Vector& z, const Matrix& A, std::vector<Vector>& L_v)
 {
   if (i == n)
   {
@@ -34,7 +33,7 @@ void generate_vertices(Index i, Index n, const Vector& z, const Matrix& A, vecto
 
 std::vector<Vector> Parallelepiped::vertices() const
 {
-  vector<Vector> L_v;
+  std::vector<Vector> L_v;
   generate_vertices(0, z.size(),z,A,L_v);
   return L_v;
 }

src/core/domains/zonotope/codac2_Zonotope.cpp

@@ -9,7 +9,6 @@
 
 #include "codac2_Zonotope.h"
 
-using namespace std;
 using namespace codac2;
 
 Zonotope::Zonotope(const Vector& z_, const Matrix& A_)

src/core/peibos/codac2_peibos.cpp

@@ -11,7 +11,6 @@
 #include "codac2_peibos_tools.h"
 #include "codac2_OctaSym_operator.h"
 
-using namespace std;
 using namespace codac2;
 
 namespace codac2
@@ -34,12 +33,12 @@ namespace codac2
     return Parallelepiped(z, A_inf);
   }
 
-  vector<Parallelepiped> PEIBOS(const AnalyticFunction<VectorType>& f, const AnalyticFunction<VectorType>& psi_0, const vector<OctaSym>& Sigma, double epsilon, bool verbose)
+  std::vector<Parallelepiped> PEIBOS(const AnalyticFunction<VectorType>& f, const AnalyticFunction<VectorType>& psi_0, const std::vector<OctaSym>& Sigma, double epsilon, bool verbose)
   {
     return PEIBOS(f, psi_0, Sigma, epsilon, Vector::zero(psi_0.output_size()), verbose);
   }
 
-  vector<Parallelepiped> PEIBOS(const AnalyticFunction<VectorType>& f, const AnalyticFunction<VectorType>& psi_0, const vector<OctaSym>& Sigma, double epsilon, const Vector& offset, bool verbose)
+  std::vector<Parallelepiped> PEIBOS(const AnalyticFunction<VectorType>& f, const AnalyticFunction<VectorType>& psi_0, const std::vector<OctaSym>& Sigma, double epsilon, const Vector& offset, bool verbose)
   {
     Index m = psi_0.input_size();
 
@@ -50,9 +49,9 @@ namespace codac2
 
     clock_t t_start = clock();
 
-    vector<Parallelepiped> output;
+    std::vector<Parallelepiped> output;
 
-    vector<IntervalVector> boxes;
+    std::vector<IntervalVector> boxes;
     double true_eps = split(IntervalVector::constant(m,{-1,1}), epsilon, boxes);
 
     for (const auto& sigma : Sigma)

src/core/peibos/codac2_peibos_tools.cpp

@@ -8,17 +8,15 @@
  */
 
 #include "codac2_peibos_tools.h"
+#include "codac2_template_tools.h"
 #include "codac2_inversion.h"
 #include "codac2_IntvFullPivLU.h"
 
-#include <iostream>
-
-using namespace std;
 using namespace codac2;
 
 namespace codac2
 {
-  double split (const IntervalVector& X, double eps, vector<IntervalVector>& boxes)
+  double split (const IntervalVector& X, double eps, std::vector<IntervalVector>& boxes)
   {
     if (X.max_diam()<=eps)
     {
@@ -58,25 +56,24 @@ namespace codac2
     Index m = A.cols();
     Index n = A.rows();
 
-    IntervalMatrix A_int (A);
-
-    IntvFullPivLU LUdec((IntervalMatrix) A_int.transpose());
+    IntvFullPivLU LUdec((Matrix) A.transpose());
     IntervalMatrix N = LUdec.kernel();
 
-    vector<int> col_indices;
+    std::vector<int> col_indices;
 
     for (int i = 0; i < m; i++)
       if (! (A.col(i).isZero()) )
         col_indices.push_back(i);
 
-    // A_reduced is composed of the non empty vectors of A_int
-    IntervalMatrix A_reduced (n, col_indices.size());
+    // A_reduced is composed of the non empty vectors of A
+
+    Matrix A_reduced (n, col_indices.size());
     for (int i = 0; i < (int) col_indices.size(); i++)
-      A_reduced.col(i) = A_int.col(col_indices[i]);
+      A_reduced.col(i) = A.col(col_indices[i]);
 
     // Construct the new matrix A_tild = [A_reduced | N]
-    IntervalMatrix A_tild (n,n);
-    A_tild << A_reduced,N;
+    Matrix A_tild (n,n);
+    A_tild << A_reduced,N.mid();
 
     IntervalMatrix Q = inverse_enclosure(A_tild.transpose() * A_tild);
 
@@ -87,42 +84,7 @@ namespace codac2
     for (int i = 0; i < m; i++)
       mult(i,i) += e_vec(i);
 
-    // From here we have an IntervalMatrix A_inf, meaning that each generator is an interval vector
-    IntervalMatrix A_inf = A_tild * mult;
-
-    // The initial parallelepiped is <y> = Y*[-1,1]^n
-    Matrix Y = A_inf.smig();
-    
-    // The following is similar to the previous operations. The difference is that here the Matrix Y is square and not interval
-
-    // rho2 is the sum of the radiuses of the circle enclosing each interval generator
-    Interval rho2 (0.);
-
-    for (int i = 0; i < n; i++)
-      rho2 += A_inf.col(i).rad().norm();
-
-    IntervalMatrix M = Y.transpose() * Y;
-
-    double eps = 1e-10 * M.norm().ub();
-    double cond_Y = condition_number(Y);
-
-    if (cond_Y>1e10)
-    {
-      cerr << "Warning: ill-conditioned matrix in the inflate process. The result may not be guaranteed " << endl;
-      M += IntervalMatrix::Identity(n,n)*eps; // if Y is ill-conditioned, we add a small term to M to make it invertible
-    }
-
-    IntervalMatrix Q2 = inverse_enclosure(M);
-    
-    IntervalMatrix Y2 = IntervalMatrix::zero(n,n);
-
-    Interval factor = sqrt(1.0 + eps * cond_Y);
-    for (int i = 0; i < n; i++)
-      if (cond_Y>1e10)
-        Y2.col(i) = Y.col(i)*(1+rho2*sqrt(Q2(i,i))*factor); // sufficient to keep the guarantee ?
-      else
-        Y2.col(i) = Y.col(i)*(1+rho2*sqrt(Q2(i,i)));
-    
-    return Y2.smag();
+    return ((A_tild.template cast<Interval>()) * mult).smag();
   }
+
 }

tests/core/domains/zonotope/codac2_tests_Parallelepiped_eval.cpp

@@ -36,11 +36,10 @@ TEST_CASE("Parallelepiped_eval")
   auto pa = f2.parallelepiped_eval(Interval(-0.1,0.1), Interval(-0.1,0.1));
   auto pb = f2.parallelepiped_eval(1.0, Interval(-1,1));
 
-  
   CHECK(Approx(pa.z,1e-6) == Vector({0,0,0}));
   CHECK(Approx(pa.A,1e-6) == Matrix({{0.14,0,0},{0,0.14,0},{0,0,0.04}}));
   CHECK(Approx(pb.z,1e-6) == Vector({1,0,1}));
-  CHECK(Approx(pb.A,1e-5) == Matrix({{0.894428,0,1.78886},{0,3.00001,0},{1.78886,0,-0.894427}}));
+  CHECK(Approx(pb.A,1e-5) == Matrix({{0.894428,0,1.78886},{0,3,0},{1.78886,0,-0.894427}}));
 
 
   double dx = 0.4;

tests/core/domains/zonotope/codac2_tests_Parallelepiped_eval.py

@@ -38,7 +38,7 @@ def test_parallelepiped_eval(self):
     self.assertTrue(Approx(pa.z,1e-6)==Vector([0,0,0]))
     self.assertTrue(Approx(pa.A,1e-6)==Matrix([[0.14,0,0],[0,0.14,0],[0,0,0.04]]))
     self.assertTrue(Approx(pb.z,1e-6)==Vector([1,0,1]))
-    self.assertTrue(Approx(pb.A,1e-5)==Matrix([[0.894428,0,1.78886],[0,3.00001,0],[1.78886,0,-0.894427]]))
+    self.assertTrue(Approx(pb.A,1e-5)==Matrix([[0.894428,0,1.78886],[0,3,0],[1.78886,0,-0.894427]]))
 
 
     dx = 0.4