From 5ba2ee0b9b25726a20435d6bad26fc01f64a90d8 Mon Sep 17 00:00:00 2001 From: SimonRohou Date: Fri, 8 Sep 2023 11:53:48 +0200 Subject: [PATCH 1/8] [core] templated Interval/Vector/Matrix --- src/core/2/domains/interval/codac2_Interval.h | 4 +- .../domains/interval/codac2_IntervalMatrix.h | 273 +++++++++++++++++- .../domains/interval/codac2_IntervalVector.h | 180 ++---------- src/core/2/domains/paving/codac2_Paving.h | 10 + src/core/2/variables/codac2_Matrix.h | 21 +- src/core/2/variables/codac2_Vector.h | 26 +- 6 files changed, 344 insertions(+), 170 deletions(-) diff --git a/src/core/2/domains/interval/codac2_Interval.h b/src/core/2/domains/interval/codac2_Interval.h index c87c3e56c..c829cb992 100644 --- a/src/core/2/domains/interval/codac2_Interval.h +++ b/src/core/2/domains/interval/codac2_Interval.h @@ -35,7 +35,7 @@ namespace Eigen IsComplex = 0, IsInteger = 0, IsSigned = 1, - RequireInitialization = 1, + RequireInitialization = 0, ReadCost = 1, AddCost = 3, MulCost = 3 @@ -49,7 +49,7 @@ namespace codac2 inline const Interval& real(const Interval& x) { return x; } inline Interval imag(const Interval&) { return 0.; } inline Interval abs(const Interval& x) { return ibex::abs(x); } - inline Interval abs2(const Interval& x) { return x*x; } + inline Interval abs2(const Interval& x) { return ibex::sqr(x); } } // namespace codac diff --git a/src/core/2/domains/interval/codac2_IntervalMatrix.h b/src/core/2/domains/interval/codac2_IntervalMatrix.h index d65a69fa1..dbcba4c66 100644 --- a/src/core/2/domains/interval/codac2_IntervalMatrix.h +++ b/src/core/2/domains/interval/codac2_IntervalMatrix.h @@ -13,18 +13,285 @@ #define __CODAC2_INTERVALMATRIX_H__ #include -#include #include #include #include +#include namespace codac2 { - template - class IntervalMatrix_ : public Eigen::Matrix + using Eigen::Dynamic; + + template + class IntervalMatrix_ : public Eigen::Matrix + { + public: + + IntervalMatrix_() + : Eigen::Matrix() + { + + } + + IntervalMatrix_(size_t nb_rows, size_t nb_cols) + : Eigen::Matrix(nb_rows, nb_cols) + { + assert(R == Dynamic || R == nb_rows); + assert(C == Dynamic || C == nb_cols); + } + + template + IntervalMatrix_(const Eigen::MatrixBase& other) + : Eigen::Matrix(other) + { } + + // This method allows you to assign Eigen expressions to IntervalMatrix_ + template + IntervalMatrix_& operator=(const Eigen::MatrixBase& other) + { + this->Eigen::Matrix::operator=(other); + return *this; + } + + constexpr size_t size() const + { + return this->Eigen::Matrix::size(); + } + + bool is_empty() const + { + for(size_t i = 0 ; i < size() ; i++) + if((this->data()+i)->is_empty()) + return true; + return false; + } + + bool intersects(const IntervalMatrix_& x) const + { + for(size_t i = 0 ; i < size() ; i++) + if(!(this->data()+i)->intersects(*(x.data()+i))) + return false; + return true; + } + + size_t thinnest_diam_index() const + { + return extr_diam_index(true); + } + + size_t largest_diam_index() const + { + return extr_diam_index(false); + } + + size_t extr_diam_index(bool min) const + { + // This code originates from the ibex-lib + // See: ibex_TemplateVector.h + // Author: Gilles Chabert + + double d = min ? POS_INFINITY : -1; // -1 to be sure that even a 0-diameter interval can be selected + size_t selected_index = -1; + bool unbounded = false; + assert(!this->is_empty() && "Diameter of an empty IntervalVector is undefined"); + + size_t i; + + for(i = 0 ; i < this->size() ; i++) + { + if((this->data()+i)->is_unbounded()) + { + unbounded = true; + if(!min) break; + } + else + { + double w = (this->data()+i)->diam(); + if(min ? wd) + { + selected_index = i; + d = w; + } + } + } + + if(min && selected_index == -1) + { + assert(unbounded); + // the selected interval is the first one. + i = 0; + } + + // The unbounded intervals are not considered if we look for the minimal diameter + // and some bounded intervals have been found (selected_index!=-1) + if(unbounded && (!min || selected_index==-1)) + { + double pt = min ? NEG_INFINITY : POS_INFINITY; // keep the point less/most distant from +oo (we normalize if the lower bound is -oo) + selected_index = i; + for(; i < this->size() ; i++) + { + if((this->data()+i)->lb() == NEG_INFINITY) + { + if((this->data()+i)->ub() == POS_INFINITY) + { + if(!min) + { + selected_index = i; + break; + } + } + if((min && (-(this->data()+i)->ub() > pt)) || (!min && (-(this->data()+i)->ub() < pt))) + { + selected_index = i; + pt = -(this->data()+i)->ub(); + } + } + else if((this->data()+i)->ub() == POS_INFINITY) + { + if((min && ((this->data()+i)->lb() > pt)) || (!min && ((this->data()+i)->lb() < pt))) + { + selected_index = i; + pt = (this->data()+i)->lb(); + } + } + } + } + + return selected_index; + } + + auto bisect(float ratio = 0.49) const + { + assert(Interval(0,1).interior_contains(ratio)); + size_t i = largest_diam_index(); + assert((this->data()+i)->is_bisectable()); + + auto p = std::make_pair(*this,*this); + auto pi = (this->data()+i)->bisect(ratio); + *(p.first.data()+i) = pi.first; + *(p.second.data()+i) = pi.second; + return p; + } + + double volume() const + { + if(this->is_empty()) + return 0.; + double v = 0.; + for(size_t i = 0 ; i < this->size() ; i++) + { + if((this->data()+i)->is_unbounded()) return POS_INFINITY; + if((this->data()+i)->is_degenerated()) return 0.; + v += log((this->data()+i)->diam()); + } + return exp(v); + } + + Matrix_ lb() const + { + assert(!this->is_empty()); // todo: use nan instead of assert? + Matrix_ lb; + for(size_t i = 0 ; i < this->size() ; i++) + *(lb.data()+i) = (this->data()+i)->lb(); + return lb; + } + + Matrix_ ub() const + { + assert(!this->is_empty()); // todo: use nan instead of assert? + Matrix_ ub; + for(size_t i = 0 ; i < this->size() ; i++) + *(ub.data()+i) = (this->data()+i)->ub(); + return ub; + } + + Matrix_ mid() const + { + assert(!this->is_empty()); // todo: use nan instead of assert? + Matrix_ m; + for(size_t i = 0 ; i < this->size() ; i++) + *(m.data()+i) = (this->data()+i)->mid(); + return m; + } + + void set_empty() + { + for(size_t i = 0 ; i < this->size() ; i++) + (this->data()+i)->set_empty(); + } + + static IntervalMatrix_ empty_set(size_t nb_rows = R, size_t nb_cols = C) + { + IntervalMatrix_ x(nb_rows, nb_cols); + x.set_empty(); + return x; + } + + auto& inflate(double r) + { + assert(r >= 0.); + for(size_t i = 0 ; i < this->size() ; i++) + (this->data()+i)->inflate(r); + return *this; + } + + auto& operator&=(const IntervalMatrix_& x) + { + if(!this->is_empty()) + { + if(x.is_empty()) + this->set_empty(); + else + for(size_t i = 0 ; i < this->size() ; i++) + *(this->data()+i) &= *(x.data()+i); + } + return *this; + } + + auto& operator|=(const IntervalMatrix_& x) + { + if(!x.is_empty()) + { + if(this->is_empty()) + *this = x; + else + for(size_t i = 0 ; i < this->size() ; i++) + *(this->data()+i) |= *(x.data()+i); + } + return *this; + } + + auto& operator+=(const IntervalMatrix_& x) + { + (*this).noalias() += x;//.template cast(); + return *this; + } + + auto& operator-=(const IntervalMatrix_& x) + { + (*this).noalias() -= x;//.template cast(); + return *this; + } + }; + + class IntervalMatrix : public IntervalMatrix_<> { public: + explicit IntervalMatrix(size_t nb_rows, size_t nb_cols) + : IntervalMatrix_<>(nb_rows, nb_cols) + { } + + IntervalMatrix(const IntervalMatrix_<>& x) + : IntervalMatrix_<>(x) + { } + + template + explicit IntervalMatrix(const Matrix_& v) + : IntervalMatrix_<>(v) + { + + } }; } // namespace codac diff --git a/src/core/2/domains/interval/codac2_IntervalVector.h b/src/core/2/domains/interval/codac2_IntervalVector.h index 65bb8167f..28488b8c9 100644 --- a/src/core/2/domains/interval/codac2_IntervalVector.h +++ b/src/core/2/domains/interval/codac2_IntervalVector.h @@ -28,144 +28,64 @@ namespace codac2 using Eigen::Dynamic; template - class IntervalVector_ : public Eigen::Matrix + class IntervalVector_ : public IntervalMatrix_ { public: IntervalVector_() - : Eigen::Matrix() + : IntervalMatrix_() { } IntervalVector_(size_t n) - : Eigen::Matrix(n) + : IntervalMatrix_(n,1) { assert(N == Dynamic || N == n); } template - explicit IntervalVector_(const Vector_& v) - : Eigen::Matrix(v.size()) + IntervalVector_(const Matrix_& v) + : IntervalMatrix_(v.size(),1) { static_assert(N == M || N == -1 || M == -1); - for(size_t i = 0 ; i < size() ; i++) + for(size_t i = 0 ; i < IntervalMatrix_::size() ; i++) (*this)[i] = v[i]; } //explicit IntervalVector_(const Interval& xi) - // : Eigen::Matrix() + // : IntervalMatrix_() //{ // (*this)(0,0) = xi; //} IntervalVector_(std::initializer_list l) - : Eigen::Matrix() + : IntervalMatrix_(l.size(),1) { - assert(l.size() == size()); + assert(l.size() == this->size()); size_t i = 0; for(const auto& li : l) (*this)[i++] = li; } + template + IntervalVector_(const IntervalMatrix_& x) + : IntervalMatrix_(x) + { + assert(M == Dynamic || M == N); + } + template IntervalVector_(const Eigen::MatrixBase& other) - : Eigen::Matrix(other) + : IntervalMatrix_(other) { } // This method allows you to assign Eigen expressions to IntervalVector_ template IntervalVector_& operator=(const Eigen::MatrixBase& other) { - this->Eigen::Matrix::operator=(other); + this->IntervalMatrix_::operator=(other); return *this; } - constexpr size_t size() const - { - return this->Eigen::Matrix::size(); - } - - bool is_empty() const - { - for(size_t i = 0 ; i < size() ; i++) - if((*this)[i].is_empty()) - return true; - return false; - } - - bool intersects(const IntervalVector_& x) const - { - // todo: case of Dynamic vs Fixed - for(size_t i = 0 ; i < size() ; i++) - if(!(*this)[i].intersects(x[i])) - return false; - return true; - } - - std::pair,IntervalVector_> bisect(float ratio = 0.49) const - { - assert(Interval(0,1).interior_contains(ratio)); - size_t i = largest_dim(); - assert((*this)[i].is_bisectable()); - - auto p = std::make_pair(*this,*this); - auto pi = (*this)[i].bisect(ratio); - p.first[i] = pi.first; - p.second[i] = pi.second; - return p; - } - - size_t largest_dim() const - { - return to_codac1(*this).extr_diam_index(false); - } - - double volume() const - { - if(is_empty()) - return 0.; - double v = 0.; - for(size_t i = 0 ; i < size() ; i++) - { - if((*this)[i].is_unbounded()) return POS_INFINITY; - if((*this)[i].is_degenerated()) return 0.; - v += log((*this)[i].diam()); - } - return exp(v); - } - - Vector_ lb() const - { - assert(!is_empty()); // todo: use nan instead of assert? - Vector_ lb(size()); - for(size_t i = 0 ; i < size() ; i++) - lb[i] = (*this)[i].lb(); - return lb; - } - - Vector_ ub() const - { - assert(!is_empty()); // todo: use nan instead of assert? - Vector_ ub(size()); - for(size_t i = 0 ; i < size() ; i++) - ub[i] = (*this)[i].ub(); - return ub; - } - - Vector_ mid() const - { - assert(!is_empty()); // todo: use nan instead of assert? - Vector_ m(size()); - for(size_t i = 0 ; i < size() ; i++) - m[i] = (*this)[i].mid(); - return m; - } - - void set_empty() - { - for(size_t i = 0 ; i < size() ; i++) - (*this)[i].set_empty(); - } - static IntervalVector_ empty_set(size_t n = N) { IntervalVector_ x(n); @@ -173,54 +93,6 @@ namespace codac2 return x; } - IntervalVector_& inflate(double r) - { - assert(r >= 0.); - for(size_t i = 0 ; i < size() ; i++) - (*this)[i].inflate(r); - return *this; - } - - IntervalVector_& operator&=(const IntervalVector_& x) - { - if(!is_empty()) - { - if(x.is_empty()) - set_empty(); - - else - for(size_t i = 0 ; i < size() ; i++) - (*this)[i] &= x[i]; - } - return *this; - } - - IntervalVector_& operator|=(const IntervalVector_& x) - { - if(!x.is_empty()) - { - if(is_empty()) - *this = x; - - else - for(size_t i = 0 ; i < size() ; i++) - (*this)[i] |= x[i]; - } - return *this; - } - - IntervalVector_& operator+=(const Vector_& x) - { - (*this).noalias() += x.template cast(); - return *this; - } - - IntervalVector_& operator-=(const Vector_& x) - { - (*this).noalias() -= x.template cast(); - return *this; - } - template IntervalVector_ subvector() const { @@ -231,7 +103,7 @@ namespace codac2 IntervalVector_<> subvector(size_t start_index, size_t end_index) const { assert(end_index >= 0 && start_index >= 0); - assert(end_index < size() && start_index <= end_index); + assert(end_index < this->size() && start_index <= end_index); IntervalVector_<> s(end_index-start_index+1); for(size_t i = 0 ; i < s.size() ; i++) @@ -245,6 +117,18 @@ namespace codac2 assert(I >= 0 && I < N && M+I <= N); this->template block(I,0) << x; } + + auto& operator+=(const IntervalVector_& x) + { + (*this).noalias() += x;//.template cast(); + return *this; + } + + auto& operator-=(const IntervalVector_& x) + { + (*this).noalias() -= x;//.template cast(); + return *this; + } }; template @@ -304,7 +188,7 @@ namespace codac2 { public: - explicit IntervalVector(int n) + explicit IntervalVector(size_t n) : IntervalVector_<>(n) { } diff --git a/src/core/2/domains/paving/codac2_Paving.h b/src/core/2/domains/paving/codac2_Paving.h index f1450fe35..027295a09 100644 --- a/src/core/2/domains/paving/codac2_Paving.h +++ b/src/core/2/domains/paving/codac2_Paving.h @@ -12,6 +12,7 @@ #ifndef __CODAC2_PAVING_H__ #define __CODAC2_PAVING_H__ +#include #include #include #include "codac2_Interval.h" @@ -41,6 +42,15 @@ namespace codac2 return _x; } + bool is_empty() const + { + if(is_leaf()) + return _x.is_empty(); + + else + return (_left && _left->is_empty()) && (_right && _right->is_empty()); + } + bool is_leaf() const { return not _left && not _right; diff --git a/src/core/2/variables/codac2_Matrix.h b/src/core/2/variables/codac2_Matrix.h index cd4cad8f3..c18cdbb82 100644 --- a/src/core/2/variables/codac2_Matrix.h +++ b/src/core/2/variables/codac2_Matrix.h @@ -18,30 +18,39 @@ namespace codac2 { + using Eigen::Dynamic; + template - class Matrix : public Eigen::Matrix + class Matrix_ : public Eigen::Matrix { public: - Matrix() + Matrix_() { } + + Matrix_(size_t nb_rows, size_t nb_cols) + : Eigen::Matrix(nb_rows, nb_cols) + { + assert(R == Dynamic || R == nb_rows); + assert(C == Dynamic || C == nb_cols); + } template - Matrix(const Eigen::MatrixBase& other) + Matrix_(const Eigen::MatrixBase& other) : Eigen::Matrix(other) { } - // This method allows you to assign Eigen expressions to MyVectorType + // This method allows you to assign Eigen expressions to Matrix_ template - Matrix& operator=(const Eigen::MatrixBase& other) + Matrix_& operator=(const Eigen::MatrixBase& other) { this->Eigen::Matrix::operator=(other); return *this; } - static Matrix eye() + static Matrix_ eye() { return Eigen::Matrix::Identity(); } diff --git a/src/core/2/variables/codac2_Vector.h b/src/core/2/variables/codac2_Vector.h index 212741fe3..61b273cd0 100644 --- a/src/core/2/variables/codac2_Vector.h +++ b/src/core/2/variables/codac2_Vector.h @@ -12,9 +12,6 @@ #ifndef __CODAC2_VECTOR_H__ #define __CODAC2_VECTOR_H__ -#include -#include -#include #include #include @@ -23,7 +20,7 @@ namespace codac2 using Eigen::Dynamic; template - class Vector_ : public Eigen::Matrix + class Vector_ : public Matrix_ { public: @@ -33,25 +30,32 @@ namespace codac2 } Vector_(size_t n) - : Eigen::Matrix(n) + : Matrix_(n,1) { assert(N == Dynamic || N == n); } - Vector_(std::initializer_list l) : Eigen::Matrix(l.size()) + Vector_(std::initializer_list l) : Matrix_(l.size(),1) { assert(N == l.size() || N == -1); size_t i = 0; for(double li : l) (*this)(i++,0) = li; } + + template + Vector_(const Matrix_& x) + : Matrix_(x) + { + assert(M == Dynamic || M == N); + } template Vector_(const Eigen::MatrixBase& other) - : Eigen::Matrix(other) + : Matrix_(other) { } - // This method allows you to assign Eigen expressions to MyVector_Type + // This method allows you to assign Eigen expressions to Vector_ template Vector_& operator=(const Eigen::MatrixBase& other) { @@ -59,9 +63,9 @@ namespace codac2 return *this; } - Matrix as_diag() const + Matrix_ as_diag() const { - return Matrix(Eigen::Matrix(this->asDiagonal())); + return Matrix_(Eigen::Matrix(this->asDiagonal())); } static Vector_ zeros() @@ -72,7 +76,7 @@ namespace codac2 }; template - Matrix diag(const Vector_ v) + Matrix_ diag(const Vector_ v) { return v.as_diag(); } From b24b033c635962a95369ee489066d31b5c1d9dc4 Mon Sep 17 00:00:00 2001 From: SimonRohou Date: Fri, 8 Sep 2023 19:53:18 +0200 Subject: [PATCH 2/8] [core] templated Interval/Vector/Matrix --- .../domains/interval/codac2_IntervalMatrix.h | 184 ++++++++++++++++-- .../domains/interval/codac2_IntervalVector.h | 28 ++- src/core/2/variables/codac2_Matrix.h | 1 + 3 files changed, 198 insertions(+), 15 deletions(-) diff --git a/src/core/2/domains/interval/codac2_IntervalMatrix.h b/src/core/2/domains/interval/codac2_IntervalMatrix.h index dbcba4c66..1c52316dd 100644 --- a/src/core/2/domains/interval/codac2_IntervalMatrix.h +++ b/src/core/2/domains/interval/codac2_IntervalMatrix.h @@ -28,7 +28,7 @@ namespace codac2 public: IntervalMatrix_() - : Eigen::Matrix() + : Eigen::Matrix(R,C) { } @@ -40,6 +40,13 @@ namespace codac2 assert(C == Dynamic || C == nb_cols); } + IntervalMatrix_(size_t nb_rows, size_t nb_cols, const Interval& x) + : IntervalMatrix_(nb_rows, nb_cols) + { + for(size_t i = 0 ; i < size() ; i++) + *(this->data()+i) = x; + } + template IntervalMatrix_(const Eigen::MatrixBase& other) : Eigen::Matrix(other) @@ -66,14 +73,127 @@ namespace codac2 return false; } + bool is_flat() const + { + if(is_empty()) return true; + for(size_t i = 0 ; i < size() ; i++) + if((this->data()+i)->is_degenerated()) // don't use diam() because of roundoff + return true; + return false; + } + + bool is_unbounded() const + { + if(is_empty()) return false; + for(size_t i = 0 ; i < size() ; i++) + if((this->data()+i)->is_unbounded()) + return true; + return false; + } + + bool is_subset(const IntervalMatrix_& x) const + { + if(is_empty()) return true; + if(x.is_empty()) return false; + for(size_t i = 0 ; i < size() ; i++) + if(!(this->data()+i)->is_subset(*(x.data()+i))) + return false; + return true; + } + + bool is_strict_subset(const IntervalMatrix_& x) const + { + if(x.is_empty()) return false; + if(is_empty()) return true; + bool one_dim_strict_subset = false; + for(size_t i = 0 ; i < size() ; i++) + { + if((this->data()+i)->is_strict_subset(*(x.data()+i))) + one_dim_strict_subset = true; + if(!(this->data()+i)->is_subset(*(x.data()+i))) + return false; + } + return one_dim_strict_subset; + } + + bool is_superset(const IntervalMatrix_& x) const + { + return x.is_subset(*this); + } + + bool is_strict_superset(const IntervalMatrix_& x) const + { + return x.is_strict_subset(*this); + } + + bool contains(const Matrix_& x) const + { + if(is_empty()) + return false; + for(size_t i = 0 ; i < size() ; i++) + if(!(this->data()+i)->contains(*(x.data()+i))) + return false; + return true; + } + + bool interior_contains(const IntervalMatrix_& x) const + { + if(is_empty()) + return false; + for(size_t i = 0 ; i < size() ; i++) + if(!(this->data()+i)->interior_contains(*(x.data()+i))) + return false; + return true; + } + bool intersects(const IntervalMatrix_& x) const { + if(is_empty() || x.is_empty()) + return false; for(size_t i = 0 ; i < size() ; i++) if(!(this->data()+i)->intersects(*(x.data()+i))) return false; return true; } + bool overlaps(const IntervalMatrix_& x) const + { + if(is_empty() || x.is_empty()) + return false; + for(size_t i = 0 ; i < size() ; i++) + if(!(this->data()+i)->overlaps(*(x.data()+i))) + return false; + return true; + } + + bool is_disjoint(const IntervalMatrix_& x) const + { + if(is_empty() || x.is_empty()) + return true; + for(size_t i = 0 ; i < size() ; i++) + if((this->data()+i)->is_disjoint(*(x.data()+i))) + return true; + return false; + } + + bool is_bisectable() const + { + for(size_t i = 0 ; i < size() ; i++) + if((this->data()+i)->is_bisectable()) + return true; + return false; + } + + double min_diam() const + { + return (this->data()+extr_diam_index(true))->diam(); + } + + double max_diam() const + { + return (this->data()+extr_diam_index(false))->diam(); + } + size_t thinnest_diam_index() const { return extr_diam_index(true); @@ -160,11 +280,10 @@ namespace codac2 return selected_index; } - auto bisect(float ratio = 0.49) const + auto bisect(size_t i, float ratio = 0.49) const { - assert(Interval(0,1).interior_contains(ratio)); - size_t i = largest_diam_index(); assert((this->data()+i)->is_bisectable()); + assert(Interval(0,1).interior_contains(ratio)); auto p = std::make_pair(*this,*this); auto pi = (this->data()+i)->bisect(ratio); @@ -173,6 +292,11 @@ namespace codac2 return p; } + auto bisect_largest_dim(float ratio = 0.49) const + { + return bisect(largest_diam_index(), ratio); + } + double volume() const { if(this->is_empty()) @@ -190,7 +314,7 @@ namespace codac2 Matrix_ lb() const { assert(!this->is_empty()); // todo: use nan instead of assert? - Matrix_ lb; + Matrix_ lb(this->rows(), this->cols()); for(size_t i = 0 ; i < this->size() ; i++) *(lb.data()+i) = (this->data()+i)->lb(); return lb; @@ -199,7 +323,7 @@ namespace codac2 Matrix_ ub() const { assert(!this->is_empty()); // todo: use nan instead of assert? - Matrix_ ub; + Matrix_ ub(this->rows(), this->cols()); for(size_t i = 0 ; i < this->size() ; i++) *(ub.data()+i) = (this->data()+i)->ub(); return ub; @@ -208,22 +332,44 @@ namespace codac2 Matrix_ mid() const { assert(!this->is_empty()); // todo: use nan instead of assert? - Matrix_ m; + Matrix_ mid(this->rows(), this->cols()); for(size_t i = 0 ; i < this->size() ; i++) - *(m.data()+i) = (this->data()+i)->mid(); - return m; + *(mid.data()+i) = (this->data()+i)->mid(); + return mid; } - void set_empty() + Matrix_ rad() const + { + assert(!this->is_empty()); // todo: use nan instead of assert? + Matrix_ rad(this->rows(), this->cols()); + for(size_t i = 0 ; i < this->size() ; i++) + *(rad.data()+i) = (this->data()+i)->rad(); + return rad; + } + + Matrix_ diam() const + { + assert(!this->is_empty()); // todo: use nan instead of assert? + Matrix_ diam(this->rows(), this->cols()); + for(size_t i = 0 ; i < this->size() ; i++) + *(diam.data()+i) = (this->data()+i)->diam(); + return diam; + } + + void init(const Interval& x) { for(size_t i = 0 ; i < this->size() ; i++) - (this->data()+i)->set_empty(); + *(this->data()+i) = x; + } + + void set_empty() + { + init(Interval::empty_set()); } static IntervalMatrix_ empty_set(size_t nb_rows = R, size_t nb_cols = C) { - IntervalMatrix_ x(nb_rows, nb_cols); - x.set_empty(); + IntervalMatrix_ x(nb_rows, nb_cols, Interval::empty_set()); return x; } @@ -261,6 +407,18 @@ namespace codac2 return *this; } + auto operator&(const IntervalMatrix_& x) + { + auto y = *this; + return y &= x; + } + + auto operator|(const IntervalMatrix_& x) + { + auto y = *this; + return y |= x; + } + auto& operator+=(const IntervalMatrix_& x) { (*this).noalias() += x;//.template cast(); diff --git a/src/core/2/domains/interval/codac2_IntervalVector.h b/src/core/2/domains/interval/codac2_IntervalVector.h index 28488b8c9..804d99811 100644 --- a/src/core/2/domains/interval/codac2_IntervalVector.h +++ b/src/core/2/domains/interval/codac2_IntervalVector.h @@ -41,6 +41,12 @@ namespace codac2 { assert(N == Dynamic || N == n); } + + IntervalVector_(size_t n, const Interval& x) + : IntervalMatrix_(n,1,x) + { + assert(N == Dynamic || N == n); + } template IntervalVector_(const Matrix_& v) @@ -57,6 +63,24 @@ namespace codac2 // (*this)(0,0) = xi; //} + IntervalVector_(size_t n, double bounds[][2]) + : IntervalMatrix_(n,1) + { + for(size_t i = 0 ; i < n ; i++) + { + if(bounds == 0) // in case the user called IntervalVector(n,0) and 0 is interpreted as NULL! + (*this)[i] = Interval::zero(); + else + (*this)[i] = Interval(bounds[i][0],bounds[i][1]); + } + } + + IntervalVector_(double bounds[][2]) + : IntervalVector_(this->size(), bounds) + { + + } + IntervalVector_(std::initializer_list l) : IntervalMatrix_(l.size(),1) { @@ -64,6 +88,7 @@ namespace codac2 size_t i = 0; for(const auto& li : l) (*this)[i++] = li; + // todo: use thias as faster? std::copy(l.begin(), l.end(), vec); } template @@ -88,8 +113,7 @@ namespace codac2 static IntervalVector_ empty_set(size_t n = N) { - IntervalVector_ x(n); - x.set_empty(); + IntervalVector_ x(n, Interval::empty_set()); return x; } diff --git a/src/core/2/variables/codac2_Matrix.h b/src/core/2/variables/codac2_Matrix.h index c18cdbb82..c2a2ccb18 100644 --- a/src/core/2/variables/codac2_Matrix.h +++ b/src/core/2/variables/codac2_Matrix.h @@ -26,6 +26,7 @@ namespace codac2 public: Matrix_() + : Eigen::Matrix() { } From 06d8a99732ca8ec2a436efb608ebc2a540ecf5e0 Mon Sep 17 00:00:00 2001 From: SimonRohou Date: Fri, 8 Sep 2023 23:39:57 +0200 Subject: [PATCH 3/8] [core] templated Interval/Vector/Matrix with tests --- src/core/2/domains/interval/codac2_Interval.h | 4 +- .../domains/interval/codac2_IntervalMatrix.h | 54 +- .../domains/interval/codac2_IntervalVector.h | 115 +- .../2/domains/interval/codac2_cart_prod.h | 97 ++ src/core/2/variables/codac2_Matrix.h | 4 +- src/core/2/variables/codac2_Vector.h | 18 +- src/core/CMakeLists.txt | 1 + tests/core/CMakeLists.txt | 57 +- tests/core/tests_codac2_intervalvector.cpp | 1166 +++++++++++++++++ 9 files changed, 1417 insertions(+), 99 deletions(-) create mode 100644 src/core/2/domains/interval/codac2_cart_prod.h create mode 100644 tests/core/tests_codac2_intervalvector.cpp diff --git a/src/core/2/domains/interval/codac2_Interval.h b/src/core/2/domains/interval/codac2_Interval.h index c829cb992..2b417a05a 100644 --- a/src/core/2/domains/interval/codac2_Interval.h +++ b/src/core/2/domains/interval/codac2_Interval.h @@ -18,6 +18,8 @@ namespace codac2 { + const double oo = POS_INFINITY; + using codac::Interval; } // namespace codac @@ -35,7 +37,7 @@ namespace Eigen IsComplex = 0, IsInteger = 0, IsSigned = 1, - RequireInitialization = 0, + RequireInitialization = 1, ReadCost = 1, AddCost = 3, MulCost = 3 diff --git a/src/core/2/domains/interval/codac2_IntervalMatrix.h b/src/core/2/domains/interval/codac2_IntervalMatrix.h index 1c52316dd..76ecb92f7 100644 --- a/src/core/2/domains/interval/codac2_IntervalMatrix.h +++ b/src/core/2/domains/interval/codac2_IntervalMatrix.h @@ -28,30 +28,31 @@ namespace codac2 public: IntervalMatrix_() - : Eigen::Matrix(R,C) + : Eigen::Matrix() { } - IntervalMatrix_(size_t nb_rows, size_t nb_cols) + explicit IntervalMatrix_(size_t nb_rows, size_t nb_cols) : Eigen::Matrix(nb_rows, nb_cols) { - assert(R == Dynamic || R == nb_rows); - assert(C == Dynamic || C == nb_cols); + assert(R == Dynamic || R == (int)nb_rows); + assert(C == Dynamic || C == (int)nb_cols); } - IntervalMatrix_(size_t nb_rows, size_t nb_cols, const Interval& x) + explicit IntervalMatrix_(size_t nb_rows, size_t nb_cols, const Interval& x) : IntervalMatrix_(nb_rows, nb_cols) { for(size_t i = 0 ; i < size() ; i++) *(this->data()+i) = x; } + // This constructor allows you to construct IntervalMatrix_ from Eigen expressions template IntervalMatrix_(const Eigen::MatrixBase& other) : Eigen::Matrix(other) { } - + // This method allows you to assign Eigen expressions to IntervalMatrix_ template IntervalMatrix_& operator=(const Eigen::MatrixBase& other) @@ -211,7 +212,7 @@ namespace codac2 // Author: Gilles Chabert double d = min ? POS_INFINITY : -1; // -1 to be sure that even a 0-diameter interval can be selected - size_t selected_index = -1; + int selected_index = -1; bool unbounded = false; assert(!this->is_empty() && "Diameter of an empty IntervalVector is undefined"); @@ -244,7 +245,7 @@ namespace codac2 // The unbounded intervals are not considered if we look for the minimal diameter // and some bounded intervals have been found (selected_index!=-1) - if(unbounded && (!min || selected_index==-1)) + if(unbounded && (!min || selected_index == -1)) { double pt = min ? NEG_INFINITY : POS_INFINITY; // keep the point less/most distant from +oo (we normalize if the lower bound is -oo) selected_index = i; @@ -381,6 +382,23 @@ namespace codac2 return *this; } + bool operator==(const IntervalMatrix_& x) const + { + if(x.size() != this->size()) + return false; + if(is_empty() || x.is_empty()) + return is_empty() && x.is_empty(); + for(size_t i = 0 ; i < this->size() ; i++) + if(*(this->data()+i) != *(x.data()+i)) + return false; + return true; + } + + bool operator!=(const IntervalMatrix_& x) const + { + return !(*this == x); + } + auto& operator&=(const IntervalMatrix_& x) { if(!this->is_empty()) @@ -407,6 +425,18 @@ namespace codac2 return *this; } + auto operator+(const IntervalMatrix_& x) + { + auto y = *this; + return y += x; + } + + auto operator-(const IntervalMatrix_& x) + { + auto y = *this; + return y -= x; + } + auto operator&(const IntervalMatrix_& x) { auto y = *this; @@ -432,6 +462,14 @@ namespace codac2 } }; + template + auto operator-(const IntervalMatrix_& x) + { + auto y = x; + y.init(0.); + return y -= x; + } + class IntervalMatrix : public IntervalMatrix_<> { public: diff --git a/src/core/2/domains/interval/codac2_IntervalVector.h b/src/core/2/domains/interval/codac2_IntervalVector.h index 804d99811..6b4a8b2f3 100644 --- a/src/core/2/domains/interval/codac2_IntervalVector.h +++ b/src/core/2/domains/interval/codac2_IntervalVector.h @@ -12,6 +12,7 @@ #ifndef __CODAC2_INTERVALVECTOR_H__ #define __CODAC2_INTERVALVECTOR_H__ +#include #include #include #include @@ -36,20 +37,24 @@ namespace codac2 : IntervalMatrix_() { } - IntervalVector_(size_t n) + explicit IntervalVector_(size_t n) : IntervalMatrix_(n,1) { - assert(N == Dynamic || N == n); + assert(N == Dynamic || N == (int)n); } - IntervalVector_(size_t n, const Interval& x) + explicit IntervalVector_(size_t n, const Interval& x) : IntervalMatrix_(n,1,x) { - assert(N == Dynamic || N == n); + assert(N == Dynamic || N == (int)n); } + + explicit IntervalVector_(const Interval& x) + : IntervalMatrix_(1,1,x) + { } template - IntervalVector_(const Matrix_& v) + explicit IntervalVector_(const Matrix_& v) : IntervalMatrix_(v.size(),1) { static_assert(N == M || N == -1 || M == -1); @@ -57,13 +62,7 @@ namespace codac2 (*this)[i] = v[i]; } - //explicit IntervalVector_(const Interval& xi) - // : IntervalMatrix_() - //{ - // (*this)(0,0) = xi; - //} - - IntervalVector_(size_t n, double bounds[][2]) + explicit IntervalVector_(size_t n, double bounds[][2]) : IntervalMatrix_(n,1) { for(size_t i = 0 ; i < n ; i++) @@ -75,7 +74,7 @@ namespace codac2 } } - IntervalVector_(double bounds[][2]) + explicit IntervalVector_(double bounds[][2]) : IntervalVector_(this->size(), bounds) { @@ -84,7 +83,7 @@ namespace codac2 IntervalVector_(std::initializer_list l) : IntervalMatrix_(l.size(),1) { - assert(l.size() == this->size()); + assert(N == Dynamic || (int)l.size() == N); size_t i = 0; for(const auto& li : l) (*this)[i++] = li; @@ -92,17 +91,18 @@ namespace codac2 } template - IntervalVector_(const IntervalMatrix_& x) + explicit IntervalVector_(const IntervalMatrix_& x) : IntervalMatrix_(x) { assert(M == Dynamic || M == N); } + // This constructor allows you to construct IntervalVector_ from Eigen expressions template IntervalVector_(const Eigen::MatrixBase& other) : IntervalMatrix_(other) { } - + // This method allows you to assign Eigen expressions to IntervalVector_ template IntervalVector_& operator=(const Eigen::MatrixBase& other) @@ -111,9 +111,10 @@ namespace codac2 return *this; } - static IntervalVector_ empty_set(size_t n = N) + static IntervalVector_ empty_set() { - IntervalVector_ x(n, Interval::empty_set()); + IntervalVector_ x; + x.set_empty(); return x; } @@ -174,40 +175,6 @@ namespace codac2 return x_; } - template - IntervalVector_ cart_prod(const IntervalVector_& x1, const Interval& x2) - { - IntervalVector_ x; - x << x1,x2; - return x; - } - - template - IntervalVector_ cart_prod(const Interval& x1, const IntervalVector_& x2) - { - IntervalVector_ x; - x << x1,x2; - return x; - } - - template - IntervalVector_ cart_prod(const IntervalVector_& x1, const IntervalVector_& x2) - { - IntervalVector_ x; - x << x1,x2; - return x; - } - - template - auto cart_prod(const T1& x1, const T2& x2, const Args&... xi) // recursive variadic function - { - auto x_ = cart_prod(x1, x2); - if constexpr(sizeof...(xi) > 0) - return cart_prod(x_, xi...); - else - return x_; - } - class IntervalVector : public IntervalVector_<> { public: @@ -216,15 +183,57 @@ namespace codac2 : IntervalVector_<>(n) { } - IntervalVector(const IntervalVector_<>& x) + explicit IntervalVector(size_t n, const Interval& x) + : IntervalVector_<>(n, x) + { } + + explicit IntervalVector(const Interval& x) + : IntervalVector_<>({x}) + { } + + explicit IntervalVector(const IntervalVector_<>& x) : IntervalVector_<>(x) { } template explicit IntervalVector(const Vector_& v) : IntervalVector_<>(v) + { } + + explicit IntervalVector(size_t n, double bounds[][2]) + : IntervalVector_<>(n, bounds) + { } + + IntervalVector(std::initializer_list l) + : IntervalVector_<>(l) + { } + + // This constructor allows you to construct IntervalVector from Eigen expressions + template + IntervalVector(const Eigen::MatrixBase& other) + : IntervalVector_<>(other) + { } + + // This method allows you to assign Eigen expressions to IntervalVector + template + IntervalVector& operator=(const Eigen::MatrixBase& other) { + this->IntervalVector_<>::operator=(other); + return *this; + } + + void resize(size_t n) + { + // With resize of Eigen, the data is reallocated and all previous values are lost. + auto save = *this; + IntervalVector_<>::resize(n); + for(size_t i = 0 ; i < min(save.size(),n) ; i++) + (*this)[i] = save[i]; + } + static IntervalVector empty_set(size_t n) + { + return IntervalVector(n, Interval::empty_set()); } }; diff --git a/src/core/2/domains/interval/codac2_cart_prod.h b/src/core/2/domains/interval/codac2_cart_prod.h new file mode 100644 index 000000000..2f499b2fe --- /dev/null +++ b/src/core/2/domains/interval/codac2_cart_prod.h @@ -0,0 +1,97 @@ +/** + * \file + * + * ---------------------------------------------------------------------------- + * \date 2023 + * \author Simon Rohou + * \copyright Copyright 2023 Codac Team + * \license This program is distributed under the terms of + * the GNU Lesser General Public License (LGPL). + */ + +#ifndef __CODAC2_CARTPROD_H__ +#define __CODAC2_CARTPROD_H__ + +#include "codac2_Interval.h" +#include "codac2_IntervalVector.h" + +namespace codac2 +{ + IntervalVector cart_prod_dyn(const Interval& x1, const Interval& x2) + { + return IntervalVector({x1,x2}); + } + + auto cart_prod_static(const Interval& x1, const Interval& x2) + { + return IntervalVector_<2>({x1,x2}); + } + + auto cart_prod_dyn(const IntervalVector& x1, const Interval& x2) + { + IntervalVector x(x1.size()+1); + x << x1,x2; + return x; + } + + template + auto cart_prod_static(const IntervalVector_& x1, const Interval& x2) + { + IntervalVector_ x; + x << x1,x2; + return x; + } + + auto cart_prod_dyn(const Interval& x1, const IntervalVector& x2) + { + IntervalVector x(x2.size()+1); + x << x1,x2; + return x; + } + + template + auto cart_prod_static(const Interval& x1, const IntervalVector_& x2) + { + IntervalVector_ x; + x << x1,x2; + return x; + } + + auto cart_prod_dyn(const IntervalVector& x1, const IntervalVector& x2) + { + IntervalVector x(x1.size()+x2.size()); + x << x1,x2; + return x; + } + + template + auto cart_prod_static(const IntervalVector_& x1, const IntervalVector_& x2) + { + IntervalVector_ x; + x << x1,x2; + return x; + } + + template + IntervalVector_ cart_prod(const T1& x1, const T2& x2, const Args&... xi) // recursive variadic function + { + auto x_ = cart_prod_static(x1, x2); + if constexpr(sizeof...(xi) > 0) + return cart_prod(x_, xi...); + else + return x_; + } + + template + IntervalVector cart_prod(const T1& x1, const T2& x2, const Args&... xi) // recursive variadic function + { + IntervalVector x_ = cart_prod_dyn(IntervalVector(x1), IntervalVector(x2)); + if constexpr(sizeof...(xi) > 0) + return cart_prod(x_, xi...); + else + return x_; + } + +} // namespace codac + +#endif \ No newline at end of file diff --git a/src/core/2/variables/codac2_Matrix.h b/src/core/2/variables/codac2_Matrix.h index c2a2ccb18..65bc70114 100644 --- a/src/core/2/variables/codac2_Matrix.h +++ b/src/core/2/variables/codac2_Matrix.h @@ -34,8 +34,8 @@ namespace codac2 Matrix_(size_t nb_rows, size_t nb_cols) : Eigen::Matrix(nb_rows, nb_cols) { - assert(R == Dynamic || R == nb_rows); - assert(C == Dynamic || C == nb_cols); + assert(R == Dynamic || R == (int)nb_rows); + assert(C == Dynamic || C == (int)nb_cols); } template diff --git a/src/core/2/variables/codac2_Vector.h b/src/core/2/variables/codac2_Vector.h index 61b273cd0..271efc37f 100644 --- a/src/core/2/variables/codac2_Vector.h +++ b/src/core/2/variables/codac2_Vector.h @@ -32,12 +32,12 @@ namespace codac2 Vector_(size_t n) : Matrix_(n,1) { - assert(N == Dynamic || N == n); + assert(N == Dynamic || N == (int)n); } Vector_(std::initializer_list l) : Matrix_(l.size(),1) { - assert(N == l.size() || N == -1); + assert(N == (int)l.size() || N == -1); size_t i = 0; for(double li : l) (*this)(i++,0) = li; @@ -99,7 +99,7 @@ namespace codac2 : Vector_<>(n) { } - Vector(const Vector_<>& x) + Vector(const Vector& x) : Vector_<>(x) { } @@ -108,11 +108,15 @@ namespace codac2 { } template - explicit Vector(const Vector_& v) + Vector(const Vector_& v) : Vector_<>(v) - { - - } + { } + + template + Vector(const Matrix_& v) + : Vector_<>(v) + { } + }; } // namespace codac diff --git a/src/core/CMakeLists.txt b/src/core/CMakeLists.txt index fdf90f992..33c832456 100644 --- a/src/core/CMakeLists.txt +++ b/src/core/CMakeLists.txt @@ -183,6 +183,7 @@ ${CMAKE_CURRENT_SOURCE_DIR}/2/domains/interval/codac2_Interval.h ${CMAKE_CURRENT_SOURCE_DIR}/2/domains/interval/codac2_IntervalMatrix.h ${CMAKE_CURRENT_SOURCE_DIR}/2/domains/interval/codac2_IntervalVector.h + ${CMAKE_CURRENT_SOURCE_DIR}/2/domains/interval/codac2_cart_prod.h ${CMAKE_CURRENT_SOURCE_DIR}/2/domains/tube/codac2_AbstractConstTube.h ${CMAKE_CURRENT_SOURCE_DIR}/2/domains/tube/codac2_AbstractSlice.cpp ${CMAKE_CURRENT_SOURCE_DIR}/2/domains/tube/codac2_AbstractSlice.h diff --git a/tests/core/CMakeLists.txt b/tests/core/CMakeLists.txt index 868ad7da5..2bedd5190 100644 --- a/tests/core/CMakeLists.txt +++ b/tests/core/CMakeLists.txt @@ -6,35 +6,36 @@ set(TESTS_NAME codac-tests-core) list(APPEND SRC_TESTS ${CMAKE_CURRENT_SOURCE_DIR}/main.cpp - ${CMAKE_CURRENT_SOURCE_DIR}/tests_codac2_tubes.cpp + #${CMAKE_CURRENT_SOURCE_DIR}/tests_codac2_tubes.cpp + ${CMAKE_CURRENT_SOURCE_DIR}/tests_codac2_intervalvector.cpp - ${CMAKE_CURRENT_SOURCE_DIR}/tests_predefined_tubes.cpp - ${CMAKE_CURRENT_SOURCE_DIR}/tests_predefined_tubes.h - ${CMAKE_CURRENT_SOURCE_DIR}/tests_arithmetic.cpp - ${CMAKE_CURRENT_SOURCE_DIR}/tests_cn.cpp - ${CMAKE_CURRENT_SOURCE_DIR}/tests_ctc_box.cpp - ${CMAKE_CURRENT_SOURCE_DIR}/tests_ctc_cart_prod.cpp - ${CMAKE_CURRENT_SOURCE_DIR}/tests_ctc_delay.cpp - ${CMAKE_CURRENT_SOURCE_DIR}/tests_ctc_deriv.cpp - ${CMAKE_CURRENT_SOURCE_DIR}/tests_ctc_chain.cpp - ${CMAKE_CURRENT_SOURCE_DIR}/tests_ctc_eval.cpp - ${CMAKE_CURRENT_SOURCE_DIR}/tests_ctc_picard.cpp - ${CMAKE_CURRENT_SOURCE_DIR}/tests_ctc_lohner.cpp - ${CMAKE_CURRENT_SOURCE_DIR}/tests_ctc_static.cpp - ${CMAKE_CURRENT_SOURCE_DIR}/tests_definition.cpp - ${CMAKE_CURRENT_SOURCE_DIR}/tests_functions.cpp - ${CMAKE_CURRENT_SOURCE_DIR}/tests_integration.cpp - ${CMAKE_CURRENT_SOURCE_DIR}/tests_operators.cpp - ${CMAKE_CURRENT_SOURCE_DIR}/tests_geometry.cpp - ${CMAKE_CURRENT_SOURCE_DIR}/tests_polygons.cpp - ${CMAKE_CURRENT_SOURCE_DIR}/tests_serialization.cpp - ${CMAKE_CURRENT_SOURCE_DIR}/tests_slices_structure.cpp - ${CMAKE_CURRENT_SOURCE_DIR}/tests_trajectory.cpp - ${CMAKE_CURRENT_SOURCE_DIR}/tests_values.cpp - ${CMAKE_CURRENT_SOURCE_DIR}/tests_sep_polygon.cpp - ${CMAKE_CURRENT_SOURCE_DIR}/tests_sep_qinterprojf.cpp - ${CMAKE_CURRENT_SOURCE_DIR}/tests_sep_fixpoint_proj.cpp - ${CMAKE_CURRENT_SOURCE_DIR}/tests_sep_polar.cpp + #${CMAKE_CURRENT_SOURCE_DIR}/tests_predefined_tubes.cpp + #${CMAKE_CURRENT_SOURCE_DIR}/tests_predefined_tubes.h + #${CMAKE_CURRENT_SOURCE_DIR}/tests_arithmetic.cpp + #${CMAKE_CURRENT_SOURCE_DIR}/tests_cn.cpp + #${CMAKE_CURRENT_SOURCE_DIR}/tests_ctc_box.cpp + #${CMAKE_CURRENT_SOURCE_DIR}/tests_ctc_cart_prod.cpp + #${CMAKE_CURRENT_SOURCE_DIR}/tests_ctc_delay.cpp + #${CMAKE_CURRENT_SOURCE_DIR}/tests_ctc_deriv.cpp + #${CMAKE_CURRENT_SOURCE_DIR}/tests_ctc_chain.cpp + #${CMAKE_CURRENT_SOURCE_DIR}/tests_ctc_eval.cpp + #${CMAKE_CURRENT_SOURCE_DIR}/tests_ctc_picard.cpp + #${CMAKE_CURRENT_SOURCE_DIR}/tests_ctc_lohner.cpp + #${CMAKE_CURRENT_SOURCE_DIR}/tests_ctc_static.cpp + #${CMAKE_CURRENT_SOURCE_DIR}/tests_definition.cpp + #${CMAKE_CURRENT_SOURCE_DIR}/tests_functions.cpp + #${CMAKE_CURRENT_SOURCE_DIR}/tests_integration.cpp + #${CMAKE_CURRENT_SOURCE_DIR}/tests_operators.cpp + #${CMAKE_CURRENT_SOURCE_DIR}/tests_geometry.cpp + #${CMAKE_CURRENT_SOURCE_DIR}/tests_polygons.cpp + #${CMAKE_CURRENT_SOURCE_DIR}/tests_serialization.cpp + #${CMAKE_CURRENT_SOURCE_DIR}/tests_slices_structure.cpp + #${CMAKE_CURRENT_SOURCE_DIR}/tests_trajectory.cpp + #${CMAKE_CURRENT_SOURCE_DIR}/tests_values.cpp + #${CMAKE_CURRENT_SOURCE_DIR}/tests_sep_polygon.cpp + #${CMAKE_CURRENT_SOURCE_DIR}/tests_sep_qinterprojf.cpp + #${CMAKE_CURRENT_SOURCE_DIR}/tests_sep_fixpoint_proj.cpp + #${CMAKE_CURRENT_SOURCE_DIR}/tests_sep_polar.cpp ) add_executable(${TESTS_NAME} ${SRC_TESTS}) diff --git a/tests/core/tests_codac2_intervalvector.cpp b/tests/core/tests_codac2_intervalvector.cpp new file mode 100644 index 000000000..7b90cefe3 --- /dev/null +++ b/tests/core/tests_codac2_intervalvector.cpp @@ -0,0 +1,1166 @@ +#include "catch_interval.hpp" +#include "vibes.h" + +#include "codac2_Vector.h" +#include "codac2_IntervalVector.h" +#include "codac2_cart_prod.h" + +using namespace Catch; +using namespace Detail; +using namespace std; +using namespace codac2; + + +TEST_CASE("Test codac2::IntervalVector") +{ + // These tests come from the IBEX library (G. Chabert) + + SECTION("cons01") + { + IntervalVector x(2); + x[0]=Interval::all_reals(); + x[1]=Interval::all_reals(); + CHECK(x==IntervalVector(2)); + CHECK(x==IntervalVector(x)); + IntervalVector y(2); + y = x; + CHECK(x==y); + } + + SECTION("cons02") + { + IntervalVector x(2); + x[0]=Interval(0,1); + x[1]=Interval(0,1); + CHECK(x==IntervalVector(2,Interval(0,1))); + CHECK(x==IntervalVector(x)); + IntervalVector y(2); + y = x; + CHECK(x==y); + } + + SECTION("cons03") + { + IntervalVector x(2); + x[0]=Interval(0,1); + x[1]=Interval(2,3); + CHECK(x==IntervalVector(x)); + IntervalVector y(2); + y = x; + CHECK(x==y); + } + + SECTION("cons04") + { + double bounds[][2] = {{0,1},{2,3}}; + IntervalVector x(2); + x[0]=Interval(0,1); + x[1]=Interval(2,3); + CHECK(x==IntervalVector(2,bounds)); + IntervalVector y(2,bounds); + y = x; + CHECK(x==y); + } + + SECTION("cons05") + { + IntervalVector x(2); + x[0].set_empty(); + x[1].set_empty(); + CHECK(x==IntervalVector::empty_set(2)); + CHECK(x.is_empty()); + } + + SECTION("consInitList") + { + IntervalVector x{ + {1.0, 2.0}, + {2.0, 3.0}, + {4} + }; + CHECK(x.size() == 3); + CHECK(x[0] == Interval(1.0, 2.0)); + CHECK(x[1] == Interval(2.0, 3.0)); + CHECK(x[2] == Interval(4.0, 4.0)); + } + + SECTION("set_empty01") + { + IntervalVector x(2); + CHECK(!x.is_empty()); + x.set_empty(); + CHECK(x.is_empty()); + } + + SECTION("is_empty01") + { + CHECK(IntervalVector::empty_set(2).is_empty()); + } + + SECTION("is_empty02") + { + CHECK(!IntervalVector(2).is_empty()); + } + + SECTION("resize01") + { + IntervalVector x(1); + x[0]=Interval(1,2); + x.resize(3); + CHECK(x.size()==3); + CHECK(x[0]==Interval(1,2)); + CHECK(x[1]==Interval::all_reals()); + CHECK(x[2]==Interval::all_reals()); + } + + SECTION("resize02") + { + IntervalVector x(1); + x[0]=Interval(1,2); + x.resize(1); + CHECK(x.size()==1); + CHECK(x[0]==Interval(1,2)); + } + + SECTION("resize03") + { + IntervalVector x(2); + x[0]=Interval(1,2); + x.set_empty(); + x.resize(3); + CHECK(x.size()==3); + CHECK(x.is_empty()); + CHECK(x[2]==Interval::all_reals()); + } + + SECTION("resize04") + { + IntervalVector x(5); + x[0]=Interval(1,2); + x[1]=Interval(3,4); + x.resize(2); + CHECK(x.size()==2); + CHECK(x[0]==Interval(1,2)); + CHECK(x[1]==Interval(3,4)); + } + + static double _x[][2]={{0,1},{2,3},{4,5}}; + + SECTION("subvector01") + { + double _x01[][2]={{0,1},{2,3}}; + CHECK(IntervalVector(3,_x).subvector(0,1)==IntervalVector(2,_x01)); + } + + SECTION("subvector02") + { + double _x12[][2]={{2,3},{4,5}}; + CHECK(IntervalVector(3,_x).subvector(1,2)==IntervalVector(2,_x12)); + } + + SECTION("subvector03") + { + double _x11[][2]={{2,3}}; + CHECK(IntervalVector(3,_x).subvector(1,1)==IntervalVector(1,_x11)); + } + + SECTION("subvector04") + { + double _x22[][2]={{4,5}}; + CHECK(IntervalVector(3,_x).subvector(2,2)==IntervalVector(1,_x22)); + } + + SECTION("subvector05") + { + CHECK(IntervalVector(3,_x).subvector(0,2)==IntervalVector(3,_x)); + } + + SECTION("subvector06") + { + CHECK(IntervalVector::empty_set(3).subvector(1,2).is_empty()); + } + + SECTION("cart_prod01") + { + CHECK(codac2::cart_prod(IntervalVector(3,_x),IntervalVector::empty_set(3)).is_empty()); + CHECK(codac2::cart_prod(IntervalVector::empty_set(3),IntervalVector(3,_x)).is_empty()); + CHECK(codac2::cart_prod(Interval(1,2),Interval(2,3),IntervalVector({{5,6},{8,9}}),Interval(5,6)) == IntervalVector({{1,2},{2,3},{5,6},{8,9},{5,6}})); + } + + SECTION("inter01") + { + double _x1[][2]={{0,2},{4,6}}; + double _x2[][2]={{1,3},{5,7}}; + double _res[][2]={{1,2},{5,6}}; + CHECK(((IntervalVector(2,_x1)) &=IntervalVector(2,_x2))==IntervalVector(2,_res)); + CHECK(((IntervalVector(2,_x1)) & IntervalVector(2,_x2))==IntervalVector(2,_res)); + } + + SECTION("staticcartprod01") + { + IntervalVector_<2> x{{1,2},{3,5}}; + IntervalVector_<3> y{{1,2},{3,5},{-oo,oo}}; + Interval z{6,7}; + CHECK(cart_prod<6>(x,y,z)==IntervalVector({{1,2},{3,5},{1,2},{3,5},{-oo,oo},{6,7}})); + CHECK(cart_prod<2>(z,z)==IntervalVector({{6,7},{6,7}})); + } + + SECTION("inter02") + { + double _x1[][2]={{0,2},{4,6}}; + double _x2[][2]={{1,3},{7,8}}; + CHECK(((IntervalVector(2,_x1)) &=IntervalVector(2,_x2)).is_empty()); + CHECK(((IntervalVector(2,_x1)) & IntervalVector(2,_x2)).is_empty()); + } + + SECTION("inter03") + { + double _x1[][2]={{0,2},{4,6}}; + CHECK(((IntervalVector(2,_x1)) &=IntervalVector::empty_set(2)).is_empty()); + CHECK(((IntervalVector(2,_x1)) & IntervalVector::empty_set(2)).is_empty()); + } + + SECTION("inter04") + { + double _x1[][2]={{0,2},{4,6}}; + CHECK(((IntervalVector::empty_set(2)) &=IntervalVector(2,_x1)).is_empty()); + CHECK(((IntervalVector::empty_set(2)) & IntervalVector(2,_x1)).is_empty()); + } + + SECTION("hull01") + { + double _x1[][2]={{0,1},{4,5}}; + double _x2[][2]={{2,3},{6,7}}; + double _res[][2]={{0,3},{4,7}}; + CHECK(((IntervalVector(2,_x1)) |=IntervalVector(2,_x2))==IntervalVector(2,_res)); + CHECK(((IntervalVector(2,_x1)) | IntervalVector(2,_x2))==IntervalVector(2,_res)); + } + + SECTION("hull02") + { + double _x1[][2]={{0,1},{4,5}}; + IntervalVector x1(2,_x1); + CHECK((x1 |= x1)==x1); + CHECK((x1 | x1)==x1); + } + + SECTION("hull03") + { + double _x1[][2]={{0,2},{4,6}}; + IntervalVector x1(2,_x1); + CHECK((x1 |=IntervalVector::empty_set(2))==x1); + CHECK((x1 | IntervalVector::empty_set(2))==x1); + } + + SECTION("hull04") + { + double _x1[][2]={{0,2},{4,6}}; + IntervalVector x1(2,_x1); + CHECK(((IntervalVector::empty_set(2)) |= x1)==x1); + CHECK(((IntervalVector::empty_set(2)) | x1)==x1); + } + + SECTION("eq01") + { + IntervalVector x(3,_x); + CHECK(x==x); + CHECK(!(x!=x)); + } + + SECTION("eq02") + { + IntervalVector x(3,_x); + double _x01[][2]={{0,1},{2,3}}; + IntervalVector x1(2,_x01); + CHECK(!(x==x1)); + CHECK(x!=x1); + } + + SECTION("eq03") + { + double _x1[][2]={{0,1},{4,5}}; + double _x2[][2]={{2,3},{6,7}}; + IntervalVector x1(2,_x1); + IntervalVector x2(2,_x2); + x1.set_empty(); + x2.set_empty(); + CHECK(x1==x2); + CHECK(!(x1!=x2)); + } + + SECTION("eq04") + { + CHECK(IntervalVector::empty_set(2)==IntervalVector::empty_set(2)); + CHECK(IntervalVector::empty_set(2)!=IntervalVector::empty_set(3)); + IntervalVector x(2); + x.set_empty(); + CHECK(IntervalVector::empty_set(2)==x); + } + + SECTION("mid01") + { + IntervalVector x(3,_x); + Vector m=x.mid(); + CHECK(m[0]==0.5); + CHECK(m[1]==2.5); + CHECK(m[2]==4.5); + } + + SECTION("is_flat01") + { + CHECK(!IntervalVector(3,_x).is_flat()); + } + + SECTION("is_flat02") + { + CHECK(IntervalVector::empty_set(3).is_flat()); + } + + SECTION("is_flat03") + { + CHECK(IntervalVector(1,Interval(0,0)).is_flat()); + CHECK(!IntervalVector(1,Interval(0,1)).is_flat()); + } + + SECTION("is_flat04") + { + double _x1[][2]={{0,1},{2,2},{3,4}}; + CHECK(IntervalVector(3,_x1).is_flat()); + } + + SECTION("is_flat05") + { + double _x1[][2]={{0,1},{2,3},{4,4}}; + CHECK(IntervalVector(3,_x1).is_flat()); + } + + SECTION("is_unbounded01") + { + CHECK(!IntervalVector::empty_set(3).is_unbounded()); + } + + SECTION("is_unbounded02") + { + double _x1[][2]={{0,1},{0,2},{NEG_INFINITY,0}}; + CHECK(IntervalVector(3,_x1).is_unbounded()); + } + + SECTION("is_unbounded03") + { + double _x1[][2]={{0,1},{0,2}}; + CHECK(!IntervalVector(2,_x1).is_unbounded()); + } + + SECTION("is_unbounded04") + { + CHECK(IntervalVector(1).is_unbounded()); + } + + SECTION("is_subset01") + { + double _x1[][2]={{0,2},{2,4}}; + double _x2[][2]={{0,1},{3,4}}; + IntervalVector x1(2,_x1); + IntervalVector x2(2,_x2); + + CHECK(x1.is_superset(x2)); + CHECK(x2.is_subset(x1)); + CHECK(x1.is_strict_superset(x2)); + //CHECK(!x2.is_strict_interior_subset(x1)); + } + + SECTION("is_subset02") + { + double _x1[][2]={{0,2},{2,4}}; + double _x2[][2]={{1,1},{3,4}}; + IntervalVector x1(2,_x1); + IntervalVector x2(2,_x2); + + CHECK(x1.is_superset(x2)); + CHECK(x2.is_subset(x1)); + CHECK(x1.is_strict_superset(x2)); + //CHECK(!x2.is_strict_interior_subset(x1)); + } + + SECTION("is_subset03") + { + double _x1[][2]={{0,2},{2,4}}; + double _x2[][2]={{0,1},{3,3}}; + IntervalVector x1(2,_x1); + IntervalVector x2(2,_x2); + + CHECK(x1.is_superset(x2)); + CHECK(x2.is_subset(x1)); + CHECK(x1.is_strict_superset(x2)); + //CHECK(!x2.is_strict_interior_subset(x1)); + } + + SECTION("is_subset04") + { + double _x1[][2]={{0,2},{2,4}}; + double _x2[][2]={{1,1},{3,3}}; + IntervalVector x1(2,_x1); + IntervalVector x2(2,_x2); + + CHECK(x1.is_superset(x2)); + CHECK(x2.is_subset(x1)); + CHECK(x1.is_strict_superset(x2)); + //CHECK(x2.is_strict_interior_subset(x1)); + } + + SECTION("is_subset05") + { + double _x1[][2]={{0,2},{2,4}}; + IntervalVector x1(2,_x1); + IntervalVector x2(IntervalVector::empty_set(2)); + + CHECK(x1.is_superset(x2)); + CHECK(x2.is_subset(x1)); + CHECK(x1.is_strict_superset(x2)); + //CHECK(x2.is_strict_interior_subset(x1)); + } + + SECTION("is_subset06") + { + double _x2[][2]={{1,1},{3,3}}; + + IntervalVector x1(IntervalVector::empty_set(2)); + IntervalVector x2(2,_x2); + + CHECK(!x1.is_superset(x2)); + CHECK(!x2.is_subset(x1)); + CHECK(!x1.is_strict_superset(x2)); + //CHECK(!x2.is_strict_interior_subset(x1)); + } + + SECTION("is_subset07") + { + double _x1[][2]={{0,2},{2,4}}; + double _x2[][2]={{1,1},{3,5}}; + + IntervalVector x1(2,_x1); + IntervalVector x2(2,_x2); + + CHECK(!x1.is_superset(x2)); + CHECK(!x2.is_subset(x1)); + CHECK(!x1.is_strict_superset(x2)); + //CHECK(!x2.is_strict_interior_subset(x1)); + } + + SECTION("extr_diam_index01") + { + double _x1[][2]={{0,2},{0,1},{0,3}}; + IntervalVector x1(3,_x1); + CHECK(x1.extr_diam_index(true)==1); + CHECK(x1.extr_diam_index(false)==2); + CHECK(x1.min_diam()==1); + CHECK(x1.max_diam()==3); + } + + SECTION("extr_diam_index02") + { + double _x1[][2]={{0,1},{0,3},{0,2}}; + IntervalVector x1(3,_x1); + CHECK(x1.extr_diam_index(true)==0); + CHECK(x1.extr_diam_index(false)==1); + CHECK(x1.min_diam()==1); + CHECK(x1.max_diam()==3); + } + + SECTION("extr_diam_index03") + { + double _x1[][2]={{0,3},{0,2},{0,1}}; + IntervalVector x1(3,_x1); + CHECK(x1.extr_diam_index(true)==2); + CHECK(x1.extr_diam_index(false)==0); + CHECK(x1.min_diam()==1); + CHECK(x1.max_diam()==3); + } + + SECTION("extr_diam_index04") + { + double _x1[][2]={{0,1},{0,2},{NEG_INFINITY,0}}; + IntervalVector x1(3,_x1); + CHECK(x1.extr_diam_index(true)==0); + CHECK(x1.extr_diam_index(false)==2); + CHECK(x1.min_diam()==1); + CHECK(x1.max_diam()==POS_INFINITY); + } + + SECTION("extr_diam_index05") + { + double _x1[][2]={{NEG_INFINITY,0}}; + IntervalVector x1(1,_x1); + CHECK(x1.extr_diam_index(true)==0); + CHECK(x1.extr_diam_index(false)==0); + CHECK(x1.min_diam()==POS_INFINITY); + CHECK(x1.max_diam()==POS_INFINITY); + } + + SECTION("extr_diam_index06") + { + double _x1[][2]={{NEG_INFINITY,0},{0,1},{NEG_INFINITY,1},{1,3}}; + IntervalVector x1(4,_x1); + CHECK(x1.extr_diam_index(true)==1); + CHECK(x1.extr_diam_index(false)==2); + CHECK(x1.min_diam()==1); + CHECK(x1.max_diam()==POS_INFINITY); + } + + SECTION("extr_diam_index07") + { + double _x1[][2]={{NEG_INFINITY,0},{-2,POS_INFINITY},{NEG_INFINITY,1}}; + IntervalVector x1(3,_x1); + CHECK(x1.extr_diam_index(true)==0); + CHECK(x1.extr_diam_index(false)==1); + CHECK(x1.min_diam()==POS_INFINITY); + CHECK(x1.max_diam()==POS_INFINITY); + } + + SECTION("extr_diam_index08") + { + double _x1[][2]={{NEG_INFINITY,0},{NEG_INFINITY,1},{-2,POS_INFINITY}}; + IntervalVector x1(3,_x1); + CHECK(x1.extr_diam_index(true)==0); + CHECK(x1.extr_diam_index(false)==2); + CHECK(x1.min_diam()==POS_INFINITY); + CHECK(x1.max_diam()==POS_INFINITY); + } + + SECTION("extr_diam_index09") + { + double _x1[][2]={{-2,POS_INFINITY},{NEG_INFINITY,0},{NEG_INFINITY,1}}; + IntervalVector x1(3,_x1); + CHECK(x1.extr_diam_index(true)==1); + CHECK(x1.extr_diam_index(false)==0); + CHECK(x1.min_diam()==POS_INFINITY); + CHECK(x1.max_diam()==POS_INFINITY); + } + + SECTION("extr_diam_index10") + { + double _x1[][2]={{-2,POS_INFINITY},{NEG_INFINITY,1},{NEG_INFINITY,0}}; + IntervalVector x1(3,_x1); + CHECK(x1.extr_diam_index(true)==2); + CHECK(x1.extr_diam_index(false)==0); + CHECK(x1.min_diam()==POS_INFINITY); + CHECK(x1.max_diam()==POS_INFINITY); + } + + SECTION("volume01") + { + double _x1[][2]={{0,1},{0,POS_INFINITY}}; + CHECK(IntervalVector(2,_x1).volume()==POS_INFINITY); + } + + SECTION("volume02") + { + double _x1[][2]={{0,1},{1,1}}; + CHECK(IntervalVector(2,_x1).volume()==0); + } + + SECTION("volume03") + { + double _x1[][2]={{0,2},{2,5},{4,8}}; + CHECK(Approx(24.0)==IntervalVector(3,_x1).volume()); + } + + SECTION("minus01") + { + double _x1[][2]={{0,3},{0,2},{0,1}}; + double _x2[][2]={{-3,0},{-2,0},{-1,0}}; + CHECK((-IntervalVector(3,_x1))==IntervalVector(3,_x2)); + } + + SECTION("minus02") + { + double _x1[][2]={{0,1},{0,POS_INFINITY}}; + double _x2[][2]={{-1,0},{NEG_INFINITY,0}}; + CHECK(-IntervalVector(2,_x1)==IntervalVector(2,_x2)); + } + + SECTION("minus03") + { + CHECK(-IntervalVector::empty_set(2)==IntervalVector::empty_set(2)); + } + + SECTION("add01") + { + double _x1[][2]={{0,3},{0,2},{0,1}}; + double _x2[][2]={{0,1},{0,1},{0,1}}; + double _x3[][2]={{0,4},{0,3},{0,2}}; + + IntervalVector x1(3,_x1); + IntervalVector x2(3,_x2); + IntervalVector x3(3,_x3); + IntervalVector e(IntervalVector::empty_set(3)); + + CHECK(x1+x2==x3); + //CHECK(x1+e,e); + CHECK((x1+e).is_empty()); + //CHECK(IntervalVector(x1)+=e,e); + CHECK((IntervalVector(x1)+=e).is_empty()); + + //CHECK(e+x1,e); + CHECK((e+x1).is_empty()); + //CHECK(e+=x1,e); + CHECK((e+=x1).is_empty()); + //CHECK(e+e,e); + CHECK((e+e).is_empty()); + //CHECK(e+=e,e); + CHECK((e+=e).is_empty()); + + CHECK((IntervalVector(x1)+=x2)==x3); + //CHECK(IntervalVector(x1)+=e,e); + CHECK((IntervalVector(x1)+=e).is_empty()); + + CHECK((IntervalVector(x2)+=x1)==x3); + } + + SECTION("sub01") + { + double _x1[][2]={{0,3},{0,2},{0,1}}; + double _x2[][2]={{0,1},{0,1},{0,1}}; + double _x3[][2]={{-1,3},{-1,2},{-1,1}}; + IntervalVector x1(3,_x1); + IntervalVector x2(3,_x2); + IntervalVector x3(3,_x3); + IntervalVector e(IntervalVector::empty_set(3)); + + CHECK(x1-x2==x3); + CHECK(x2-x1==-x3); + //CHECK(x1-e,e); + CHECK((x1-e).is_empty()); + //CHECK(IntervalVector(x1)-=e,e); + CHECK((IntervalVector(x1)-=e).is_empty()); + + //CHECK(e-x1,e); + CHECK((e-x1).is_empty()); + //CHECK(e-=x1,e); + CHECK((e-=x1).is_empty()); + //CHECK(e-e,e); + CHECK((e-e).is_empty()); + //CHECK(e-=e,e); + CHECK((e-=e).is_empty()); + + CHECK((IntervalVector(x1)-=x2)==x3); + CHECK((IntervalVector(x2)-=x1)==-x3); + } +} + +#if 0 + + + +void TestIntervalVector::compl01() { + double _b[][2]={{0,1},{0,1}}; + IntervalVector b(2,_b); + IntervalVector* c; + int n=b.complementary(c); + + CHECK(n==4); + CHECK(c[0].size()==2); + + CHECK(c[0][0]==Interval::neg_reals()); + CHECK(c[0][1]==Interval::all_reals()); + + CHECK(c[1][0]==Interval(1,POS_INFINITY)); + CHECK(c[1][1]==Interval::all_reals()); + + CHECK(c[2][0]==Interval(0,1)); + CHECK(c[2][1]==Interval::neg_reals()); + + CHECK(c[3][0]==Interval(0,1)); + CHECK(c[3][1]==Interval(1,POS_INFINITY)); + + delete[] c; +} + +/** + * complementary of an empty box = (-oo,oo)x...(-oo,oo) + */ +void TestIntervalVector::compl02() { + + IntervalVector* c; + int n=IntervalVector::empty_set(2).complementary(c); + CHECK(n==1); + CHECK(c[0].size()==2); + + CHECK(c[0][0]==Interval::all_reals()); + CHECK(c[0][1]==Interval::all_reals()); + + delete[] c; +} + +bool TestIntervalVector::test_diff(int n, double _x[][2], double _y[][2], int m, double _z[][2], bool compactness, bool debug) { + IntervalVector x(n,_x); + IntervalVector y(n,_y); + IntervalMatrix mz(m,n,_z); + IntervalVector* c; + int nn=x.diff(y,c,compactness); + if (debug) { + cout << x << " diff " << y << " gives:" << endl; + for (int i=0; i& expected, bool compactness = true) +{ + auto c = x.diff(y, compactness); - CHECK(c[0][0]==Interval::all_reals()); - CHECK(c[0][1]==Interval::all_reals()); + CHECK(!c.empty()); + CHECK(c.size()==expected.size()); + CHECK(c.front().size()==x.size()); - delete[] c; -} + auto it = c.begin(); + while(it != c.end()) + { + bool is_same = false; + for(const auto& ri : expected) + if(ri == *it) + { + is_same = true; + break; + } -bool TestIntervalVector::test_diff(int n, double _x[][2], double _y[][2], int m, double _z[][2], bool compactness, bool debug) { - IntervalVector x(n,_x); - IntervalVector y(n,_y); - IntervalMatrix mz(m,n,_z); - IntervalVector* c; - int nn=x.diff(y,c,compactness); - if (debug) { - cout << x << " diff " << y << " gives:" << endl; - for (int i=0; i Date: Mon, 11 Sep 2023 15:57:14 +0200 Subject: [PATCH 5/8] [core] testing IntervalMatrix --- .../domains/interval/codac2_IntervalMatrix.h | 78 ++- .../domains/interval/codac2_IntervalVector.h | 31 +- tests/core/CMakeLists.txt | 1 + tests/core/tests_codac2_intervalmatrix.cpp | 483 ++++++++++++++++++ 4 files changed, 565 insertions(+), 28 deletions(-) create mode 100644 tests/core/tests_codac2_intervalmatrix.cpp diff --git a/src/core/2/domains/interval/codac2_IntervalMatrix.h b/src/core/2/domains/interval/codac2_IntervalMatrix.h index 559b64d79..e74721fb4 100644 --- a/src/core/2/domains/interval/codac2_IntervalMatrix.h +++ b/src/core/2/domains/interval/codac2_IntervalMatrix.h @@ -51,6 +51,42 @@ namespace codac2 for(size_t i = 0 ; i < size() ; i++) *(this->data()+i) = x; } + + explicit IntervalMatrix_(size_t nb_rows, size_t nb_cols, double bounds[][2]) + : IntervalMatrix_(nb_rows, nb_cols) + { + size_t k = 0; + for(size_t i = 0 ; i < nb_rows ; i++) + for(size_t j = 0 ; j < nb_cols ; j++) + { + if(bounds == 0) // in case the user called IntervalVector(n,0) and 0 is interpreted as NULL! + (*this)(i,j) = Interval::zero(); + else + (*this)(i,j) = Interval(bounds[k][0],bounds[k][1]); + k++; + } + } + + IntervalMatrix_(std::initializer_list> l) + : IntervalMatrix_() + { + assert(R == Dynamic || (int)l.size() == R); + int cols = -1; + for(const auto& ri : l) { + assert(cols == -1 || cols == (int)ri.size()); + cols = (int)ri.size(); + } + this->resize(l.size(),cols); + size_t i = 0; + for(const auto& ri : l) + { + size_t j = 0; + for(const auto& ci : ri) + (*this)(i,j++) = ci; + i++; + } + // todo: use thias as faster? std::copy(l.begin(), l.end(), vec); + } // This constructor allows you to construct IntervalMatrix_ from Eigen expressions template @@ -71,6 +107,16 @@ namespace codac2 return this->Eigen::Matrix::size(); } + void resize(size_t nb_rows, size_t nb_cols) + { + // With resize of Eigen, the data is reallocated and all previous values are lost. + auto copy = *this; + this->Eigen::Matrix::resize(nb_rows, nb_cols); + for(size_t i = 0 ; i < min(copy.rows(),nb_rows) ; i++) + for(size_t j = 0 ; j < min(copy.cols(),nb_cols) ; j++) + (*this)(i,j) = copy(i,j); + } + bool is_empty() const { for(size_t i = 0 ; i < size() ; i++) @@ -79,6 +125,12 @@ namespace codac2 return false; } + static IntervalMatrix_ empty_set(size_t nb_rows = R, size_t nb_cols = C) + { + IntervalMatrix_ x(nb_rows, nb_cols, Interval::empty_set()); + return x; + } + bool is_flat() const { if(is_empty()) return true; @@ -369,12 +421,6 @@ namespace codac2 init(Interval::empty_set()); } - static IntervalMatrix_ empty_set(size_t nb_rows = R, size_t nb_cols = C) - { - IntervalMatrix_ x(nb_rows, nb_cols, Interval::empty_set()); - return x; - } - auto& inflate(double r) { assert(r >= 0.); @@ -385,7 +431,7 @@ namespace codac2 bool operator==(const IntervalMatrix_& x) const { - if(x.size() != this->size()) + if(x.size() != this->size() || x.rows() != this->rows() || x.cols() != this->cols()) return false; if(is_empty() || x.is_empty()) return is_empty() && x.is_empty(); @@ -491,15 +537,31 @@ namespace codac2 : IntervalMatrix_<>(nb_rows, nb_cols) { } + + explicit IntervalMatrix(size_t nb_rows, size_t nb_cols, const Interval& x) + : IntervalMatrix_<>(nb_rows, nb_cols, x) + { } + + explicit IntervalMatrix(size_t nb_rows, size_t nb_cols, double bounds[][2]) + : IntervalMatrix_<>(nb_rows, nb_cols, bounds) + { } + IntervalMatrix(const IntervalMatrix_<>& x) : IntervalMatrix_<>(x) { } + IntervalMatrix(std::initializer_list> l) + : IntervalMatrix_<>(l) + { } + template explicit IntervalMatrix(const Matrix_& v) : IntervalMatrix_<>(v) - { + { } + static IntervalMatrix empty_set(size_t nb_rows, size_t nb_cols) + { + return IntervalMatrix_<>::empty_set(nb_rows,nb_cols); } }; diff --git a/src/core/2/domains/interval/codac2_IntervalVector.h b/src/core/2/domains/interval/codac2_IntervalVector.h index 072ac383e..40a1e8cce 100644 --- a/src/core/2/domains/interval/codac2_IntervalVector.h +++ b/src/core/2/domains/interval/codac2_IntervalVector.h @@ -68,17 +68,9 @@ namespace codac2 } explicit IntervalVector_(size_t n, double bounds[][2]) - : IntervalMatrix_(n,1) - { - for(size_t i = 0 ; i < n ; i++) - { - if(bounds == 0) // in case the user called IntervalVector(n,0) and 0 is interpreted as NULL! - (*this)[i] = Interval::zero(); - else - (*this)[i] = Interval(bounds[i][0],bounds[i][1]); - } - } - + : IntervalMatrix_(n,1,bounds) + { } + explicit IntervalVector_(double bounds[][2]) : IntervalVector_(this->size(), bounds) { @@ -118,9 +110,12 @@ namespace codac2 static IntervalVector_ empty_set(size_t n = N) { - IntervalVector_ x(n); - x.set_empty(); - return x; + return IntervalMatrix_::empty_set(n,1); + } + + void resize(size_t n) + { + this->IntervalMatrix_::resize(n,1); } template @@ -298,16 +293,12 @@ namespace codac2 void resize(size_t n) { - // With resize of Eigen, the data is reallocated and all previous values are lost. - auto save = *this; - IntervalVector_<>::resize(n); - for(size_t i = 0 ; i < min(save.size(),n) ; i++) - (*this)[i] = save[i]; + this->IntervalVector_<>::resize(n); } static IntervalVector empty_set(size_t n) { - return IntervalVector_<>::empty_set(n); + return IntervalMatrix_<>::empty_set(n,1); } }; diff --git a/tests/core/CMakeLists.txt b/tests/core/CMakeLists.txt index 2bedd5190..4cc10e560 100644 --- a/tests/core/CMakeLists.txt +++ b/tests/core/CMakeLists.txt @@ -8,6 +8,7 @@ list(APPEND SRC_TESTS ${CMAKE_CURRENT_SOURCE_DIR}/main.cpp #${CMAKE_CURRENT_SOURCE_DIR}/tests_codac2_tubes.cpp ${CMAKE_CURRENT_SOURCE_DIR}/tests_codac2_intervalvector.cpp + ${CMAKE_CURRENT_SOURCE_DIR}/tests_codac2_intervalmatrix.cpp #${CMAKE_CURRENT_SOURCE_DIR}/tests_predefined_tubes.cpp #${CMAKE_CURRENT_SOURCE_DIR}/tests_predefined_tubes.h diff --git a/tests/core/tests_codac2_intervalmatrix.cpp b/tests/core/tests_codac2_intervalmatrix.cpp new file mode 100644 index 000000000..8b8640a27 --- /dev/null +++ b/tests/core/tests_codac2_intervalmatrix.cpp @@ -0,0 +1,483 @@ +#include "catch_interval.hpp" +#include "vibes.h" + +#include "codac2_Matrix.h" +#include "codac2_IntervalVector.h" +#include "codac2_IntervalMatrix.h" + +using namespace Catch; +using namespace Detail; +using namespace std; +using namespace codac2; + +// Most of these tests come from the IBEX library (G. Chabert) +// They have been revised to fit the codac2::IntervalMatrix class + +IntervalMatrix M1() +{ + IntervalMatrix m(2,3); + double _r1[][2]={{0,1},{0,2},{0,3}}; + double _r2[][2]={{-1,0},{-2,0},{-3,0}}; + IntervalVector r1(3,_r1); + IntervalVector r2(3,_r2); + m.row(0)=r1; + m.row(1)=r2; + return m; +} + + +IntervalMatrix M2() // the transpose of M1 +{ + IntervalMatrix m(3,2); + double _c1[][2]={{0,1},{-1,0}}; + double _c2[][2]={{0,2},{-2,0}}; + double _c3[][2]={{0,3},{-3,0}}; + IntervalVector c1(2,_c1); + IntervalVector c2(2,_c2); + IntervalVector c3(2,_c3); + m.row(0)=c1; + m.row(1)=c2; + m.row(2)=c3; + return m; +} + +IntervalMatrix M3() // non-null intersection with M1 +{ + IntervalMatrix m(2,3); + double _r1[][2]={{1,2},{1,2},{2,4}}; + double _r2[][2]={{-2,-1},{-2,-1},{-4,-2}}; + IntervalVector r1(3,_r1); + IntervalVector r2(3,_r2); + m.row(0)=r1; + m.row(1)=r2; + return m; +} + + +TEST_CASE("Tests from IBEX IntervalMatrix") +{ + SECTION("eq01") + { + IntervalMatrix m(2,3); + IntervalMatrix m2(3,2); + CHECK(m!=m2); + CHECK(!(m==m2)); + } + + SECTION("eq02") + { + IntervalMatrix m(3,2); + IntervalMatrix m2(2,2); + CHECK(m!=m2); + CHECK(!(m==m2)); + } + + SECTION("eq03") + { + IntervalMatrix m(2,3); + IntervalMatrix m2(2,3); + + CHECK(m.rows()==2); + CHECK(m.cols()==3); + CHECK(m2.rows()==2); + CHECK(m2.cols()==3); + + m(0,0)=1; + m(0,1)=2; + m(0,2)=3; + m(1,0)=4; + m(1,1)=5; + m(1,2)=6; + m2(0,0)=1; + m2(0,1)=2; + m2(0,2)=3; + m2(1,0)=4; + m2(1,1)=5; + m2(1,2)=6; + + CHECK(m==m2); + CHECK(!(m!=m2)); + + m2(1,2)=7; + CHECK(m!=m2); + CHECK(!(m==m2)); + } + + SECTION("eq04") + { + IntervalMatrix m(2,3); + IntervalMatrix m2(2,3); + m(1,1)=-1; + m2(1,1)=-2; + CHECK(m!=m2); + CHECK(!(m==m2)); + m.set_empty(); + m2.set_empty(); + CHECK(m==m2); + CHECK(!(m!=m2)); + } + + SECTION("cons01") + { + IntervalMatrix m(2,3); + CHECK(m.rows()==2); + CHECK(m.cols()==3); + CHECK(m(0,0)==Interval::all_reals()); + CHECK(m(0,1)==Interval::all_reals()); + CHECK(m(0,2)==Interval::all_reals()); + CHECK(m(1,0)==Interval::all_reals()); + CHECK(m(1,1)==Interval::all_reals()); + CHECK(m(1,2)==Interval::all_reals()); + CHECK(m==IntervalMatrix(m)); + CHECK(m==(IntervalMatrix(2,3)=m)); + } + + SECTION("cons02") + { + IntervalMatrix m(2,3); + double _r1[][2]={{0,1},{0,2},{0,3}}; + double _r2[][2]={{-1,0},{-2,0},{-3,0}}; + IntervalVector r1(3,_r1); + IntervalVector r2(3,_r2); + m.row(0) = r1; + m.row(1) = r2; + + double _c1[][2]={{0,1},{-1,0}}; + double _c2[][2]={{0,2},{-2,0}}; + double _c3[][2]={{0,3},{-3,0}}; + IntervalVector c1(2,_c1); + IntervalVector c2(2,_c2); + IntervalVector c3(2,_c3); + + CHECK(m.rows()==2); + CHECK(m.cols()==3); + // not supported CHECK(m[0]==r1); + // not supported CHECK(m[1]==r2); +// not working CHECK(m.row(0)==r1); +// not working CHECK(m.row(1)==r2); + CHECK(m.col(0)==c1); + CHECK(m.col(1)==c2); + CHECK(m.col(2)==c3); + CHECK(m(0,0)==Interval(0,1)); + CHECK(m(0,1)==Interval(0,2)); + CHECK(m(0,2)==Interval(0,3)); + CHECK(m(1,0)==Interval(-1,0)); + CHECK(m(1,1)==Interval(-2,0)); + CHECK(m(1,2)==Interval(-3,0)); + CHECK(m==IntervalMatrix(m)); + CHECK(m==(IntervalMatrix(2,3)=m)); + } + + SECTION("cons03") + { + Interval x(-1,2); + IntervalMatrix m(2,3,x); + + CHECK(m.rows()==2); + CHECK(m.cols()==3); + for (int i=0; i<2; i++) { + for (int j=0; j<3; j++) + CHECK(m(i,j)==x); + } + + CHECK(m==IntervalMatrix(m)); + CHECK(m==(IntervalMatrix(2,3)=m)); + } + + SECTION("cons04") + { + double _m[][2]={ {0,1}, {0,2}, {0,3}, + {-1,0},{-2,0},{-3,0} }; + IntervalMatrix m(2,3,_m); + CHECK(m==M1()); + } + + SECTION("consInitList") + { + IntervalMatrix m{ + {{0,1}, {0,2}, {0,3}}, + {{-1,0},{-2,0},{-3,0}} + }; + CHECK(m == M1()); + } + + SECTION("empty01") + { + CHECK(IntervalMatrix::empty_set(2,3).rows()==2); + CHECK(IntervalMatrix::empty_set(2,3).cols()==3); + CHECK(IntervalMatrix(IntervalMatrix::empty_set(2,3))==IntervalMatrix::empty_set(2,3)); + CHECK((IntervalMatrix(2,3)=IntervalMatrix::empty_set(2,3))==IntervalMatrix::empty_set(2,3)); + } + + SECTION("is_empty01") + { + CHECK(!IntervalMatrix(2,3).is_empty()); + } + + SECTION("is_empty02") + { + CHECK(IntervalMatrix::empty_set(2,3).is_empty()); + } + + SECTION("set_empty01") + { + IntervalMatrix m(2,3); + m.set_empty(); + CHECK(m.is_empty()); + } + + // intersection of a matrix with itself + SECTION("inter01") + { + CHECK((M1()&=M1())==M1()); + } + + // intersection of two overlapping matrices + SECTION("inter02") + { + double _m[][2]={{1,1}, {1,2}, {2,3}, + {-1,-1},{-2,-1},{-3,-2}}; + + CHECK((M1()&=M3())==IntervalMatrix(2,3,_m)); + } + + // intersection of two non-overlapping matrices + SECTION("inter03") + { + IntervalMatrix m3(M3()); + m3(1,2)=Interval(-5,-4); + CHECK((M1()&=m3).is_empty()); + } + + SECTION("set_col01") + { + IntervalMatrix m(M1()); + + IntervalVector v(2); + v[0]=Interval(1,2); + v[1]=Interval(-2,-1); + + m.col(1)=v; + + double _m2[][2]={ {0,1}, {1,2}, {0,3}, + {-1,0},{-2,-1},{-3,0} }; + IntervalMatrix m2(2,3,_m2); + + CHECK(m==m2); + } + + SECTION("rows01") + { + CHECK(M1().block(0,0,2,3)==M1()); + } + + SECTION("rows02") + { + double _r0[][2]={ {0,1}, {0,2}, {0,3} }; + CHECK(M1().block(0,0,1,3)==IntervalMatrix(1,3,_r0)); + } + + SECTION("rows03") + { + double _r1[][2]={ {-1,0},{-2,0},{-3,0} }; + CHECK(M1().block(1,0,1,3)==IntervalMatrix(1,3,_r1)); + } + + SECTION("cols01") + { + CHECK(M1().block(0,0,2,3)==M1()); + } + + SECTION("cols02") + { + double _c0[][2]={ {0,1}, {-1,0} }; + CHECK(M1().block(0,0,2,1)==IntervalMatrix(2,1,_c0)); + } + + SECTION("cols03") + { + double _c1[][2]={ {0,2}, {-2,0} }; + CHECK(M1().block(0,1,2,1)==IntervalMatrix(2,1,_c1)); + } + + SECTION("cols04") + { + double _c2[][2]={ {0,3}, {-3,0} }; + CHECK(M1().block(0,2,2,1)==IntervalMatrix(2,1,_c2)); + } + + SECTION("cols04") + { + double _c12[][2]={ {0,2}, {0,3}, {-2,0}, {-3,0} }; + CHECK(M1().block(0,1,2,2)==IntervalMatrix(2,2,_c12)); + } + + SECTION("resize01") + { + IntervalMatrix m(2,2); + double _r1[][2]={{0,1},{0,2}}; + double _r2[][2]={{-1,0},{-2,0}}; + IntervalVector r1(2,_r1); + IntervalVector r2(2,_r2); + m.row(0)=r1; + m.row(1)=r2; + m.resize(2,3); + m(0,2)=Interval(0,3); + m(1,2)=Interval(-3,0); + + CHECK(m==M1()); + } + + SECTION("resize02") + { + IntervalMatrix m(1,3); + double _r1[][2]={{0,1},{0,2},{0,3}}; + IntervalVector r1(3,_r1); + m.row(0)=r1; + m.resize(2,3); + m(1,0)=Interval(-1,0); + m(1,1)=Interval(-2,0); + m(1,2)=Interval(-3,0); + + CHECK(m==M1()); + } + + SECTION("resize03") + { + IntervalMatrix e(IntervalMatrix::empty_set(1,1)); + e.resize(2,3); + CHECK(e.is_empty()); + } + + SECTION("minus01") + { + IntervalMatrix m(M1()); + IntervalMatrix m2(-m); + for (int i=0; i<2; i++) { + for (int j=0; j<3; j++) { + CHECK(m2(i,j)==-m(i,j)); + } + } + } + + SECTION("minus02") + { + CHECK(-IntervalMatrix::empty_set(2,3).is_empty()); + } + + SECTION("add01") + { + IntervalMatrix m(M1()); + IntervalMatrix m2(m+m); + + for (int i=0; i<2; i++) { + for (int j=0; j<3; j++) { + CHECK(m2(i,j)==m(i,j)+m(i,j)); + } + } + + CHECK(m2==(IntervalMatrix(m)+=m)); + } + + SECTION("add02") + { + IntervalMatrix m1(IntervalMatrix::empty_set(2,3)); + IntervalMatrix m2(2,3); + + CHECK((m1+m2).is_empty()); + CHECK((m1+=m2).is_empty()); + CHECK((m2+=m1).is_empty()); + } + + SECTION("sub01") + { + IntervalMatrix m(M1()); + IntervalMatrix m2(m-m); + for (int i=0; i<2; i++) { + for (int j=0; j<3; j++) { + CHECK(m2(i,j)==m(i,j)-m(i,j)); + } + } + + CHECK(m2==(IntervalMatrix(m)-=m)); + } + + SECTION("sub02") + { + IntervalMatrix m1(IntervalMatrix::empty_set(2,3)); + IntervalMatrix m2(2,3); + + CHECK((m1-m2).is_empty()); + CHECK((m1-=m2).is_empty()); + CHECK((m2-=m1).is_empty()); + } + + SECTION("mul01") + { + IntervalMatrix m(M1()); + IntervalMatrix m2(M2()); + IntervalMatrix m3(m*m2); + CHECK(m3.rows()==2); + CHECK(m3.cols()==2); + + for (int i=0; i<2; i++) { + for (int j=0; j<2; j++) + CHECK(m3(i,j)==m(i,0)*m2(0,j)+m(i,1)*m2(1,j)+m(i,2)*m2(2,j)); + } + + CHECK(m3==(IntervalMatrix(m)*=m2)); + } + + SECTION("mul02") + { + IntervalMatrix m1(IntervalMatrix::empty_set(2,3)); + IntervalMatrix m2(3,2); + + CHECK(IntervalMatrix(m1*m2).is_empty()); + CHECK(IntervalMatrix(m1*=m2).is_empty()); + CHECK(IntervalMatrix(m2*=m1).is_empty()); + } + + SECTION("put01") + { + // deprecated in codac (use Eigen instead) IntervalMatrix M1=2*Matrix::eye(3); + // deprecated in codac (use Eigen instead) IntervalVector V1(3); + // deprecated in codac (use Eigen instead) V1[0]=3; V1[1]=4; V1[2]=5; + // deprecated in codac (use Eigen instead) IntervalMatrix res(4,4); + // deprecated in codac (use Eigen instead) res.put(0,0,M1); + // deprecated in codac (use Eigen instead) res.put(0,3,V1,false); + // deprecated in codac (use Eigen instead) res.put(3,0,Vector::ones(3),true); + // deprecated in codac (use Eigen instead) res[3][3]=6; + // deprecated in codac (use Eigen instead) double _expected[16] = { 2,0,0,3, + // deprecated in codac (use Eigen instead) 0,2,0,4, + // deprecated in codac (use Eigen instead) 0,0,2,5, + // deprecated in codac (use Eigen instead) 1,1,1,6 }; + // deprecated in codac (use Eigen instead) CHECK(res==(Matrix(4,4,_expected))); + } +} + +#if 0 + +// Tests from IBEX that are not (yet) considered in Codac: + +void TestIntervalMatrix::rad01() { + RNG::srand(1); + IntervalMatrix M=Matrix::rand(2); + Matrix R=M.rad(); + CHECK(R[0][0]==M[0][0].rad()); + CHECK(R[0][1]==M[0][1].rad()); + CHECK(R[1][0]==M[1][0].rad()); + CHECK(R[1][1]==M[1][1].rad()); +} + +void TestIntervalMatrix::diam01() { + RNG::srand(1); + IntervalMatrix M=Matrix::rand(2); + Matrix R=M.diam(); + CHECK(R[0][0]==M[0][0].diam()); + CHECK(R[0][1]==M[0][1].diam()); + CHECK(R[1][0]==M[1][0].diam()); + CHECK(R[1][1]==M[1][1].diam()); +} + +#endif \ No newline at end of file From e17547af69f0142f172b310cc5d23a4dfd39f564 Mon Sep 17 00:00:00 2001 From: damien-masse Date: Tue, 12 Dec 2023 15:24:26 +0100 Subject: [PATCH 6/8] mise a jour mineure --- .../2/contractors/codac2_CtcDiffInclusion.cpp | 10 ++- src/core/2/integration/codac2_IParals.cpp | 70 ++++++++++++++++--- src/core/2/integration/codac2_IParals.h | 31 +++++++- src/core/2/integration/codac2_expIMat.cpp | 41 ++++++++++- src/core/2/integration/codac2_expIMat.h | 10 +++ 5 files changed, 148 insertions(+), 14 deletions(-) diff --git a/src/core/2/contractors/codac2_CtcDiffInclusion.cpp b/src/core/2/contractors/codac2_CtcDiffInclusion.cpp index 6fa331052..3a8992eb4 100644 --- a/src/core/2/contractors/codac2_CtcDiffInclusion.cpp +++ b/src/core/2/contractors/codac2_CtcDiffInclusion.cpp @@ -155,9 +155,14 @@ namespace codac2 double tsteps = timeslice.diam(); + global_exp(jac,tsteps,true,time_dependent, ExpM,invExpM,tauExpM,IExpM,tauIExpM, VExpM,tauVExpM,IntAbs); +// std::cout << "jacobian : " << jac << "\n"; +// std::cout << "ExpM : " << ExpM << "\n"; +// std::cout << "IExpM : " << IExpM << "\n"; +// std::cout << "IntAbs : " << IntAbs << "\n"; /* other variables which needs to be kept */ Vector cent_tdiff(dim); Vector fun_evalc(dim); @@ -250,7 +255,10 @@ namespace codac2 } for (int i=0;iis_gate()) { IParals end = x.next_slice_ptr()->codomain(); end &= g_t1; /* ensure contractance */ x.next_slice_ptr()->set(end); } + x.set(frame); /* first the gate, then the frame... */ } if((t_propa & TimePropag::BACKWARD) && (!x.output_gate().is_unbounded())) diff --git a/src/core/2/integration/codac2_IParals.cpp b/src/core/2/integration/codac2_IParals.cpp index f30b15ddb..1292c73f9 100644 --- a/src/core/2/integration/codac2_IParals.cpp +++ b/src/core/2/integration/codac2_IParals.cpp @@ -24,6 +24,7 @@ namespace codac2 { + bool IParals::display_inv=false; IParals::IParals(int dim) : dim(dim), empty(false), nbmat(0), mats(), Vrhs(1) @@ -106,6 +107,23 @@ namespace codac2 { } return true; } + bool IParals::intersects(const IntervalVector& iv) const { + assert (dim>0); + if (this->empty) return false; + if (!this->bbox().intersects(iv)) return false; + for (unsigned int i=0;inbmat;i++) { + if (!this->rhs(i).intersects(this->Imat(i)*iv)) return false; + } + return true; + } + bool IParals::intersects(const IParals& ip) const { + assert (dim>0); + if (this->empty) return false; + if (!this->intersects(ip.bbox())) return false; + if (!ip.intersects(this->bbox())) return false; + return true; + } + const IntervalMatrix& IParals::getMat(unsigned int i) const { assert (i>=0 && imat(i); @@ -352,6 +370,7 @@ namespace codac2 { if (this->bbox().is_empty()) { this->set_empty(); return *this; } if (this->nbmat==0) { for (unsigned int i=0;iborrow_base(iv,i,iv.rhs(i)); } this->simplify(); @@ -360,28 +379,44 @@ namespace codac2 { unsigned int nbcommun=0; for (unsigned int i=0;inbmat;i++) { this->rhs(i) &= this->Imat(i) * iv.bbox(); - for (unsigned int j=0;jrhs(i).is_empty()) { this->set_empty(); return *this; } + } + bool useful[iv.nbmat]; + int nbuseful=iv.nbmat; + int b=-1; /* last useful */ + for (unsigned int j=0;jbbox(); + for (unsigned int i=0;inbmat;i++) { if (this->mats[i]==iv.mats[j]) { this->rhs(i) &= iv.rhs(j); nbcommun++; + useful[j]=false; + nbuseful--; + break; + } + this->rhs(i) &= (this->Imat(i)*iv.mat(j)) * iv.rhs(j); + if (this->rhs(i).is_empty()) { this->set_empty(); return *this; } + checkuseful &= (iv.Imat(j)*this->mat(i))*this->rhs(i); + if (checkuseful.is_subset(iv.rhs(j))) { + useful[j]=false; + nbuseful--; + break; } - else - this->rhs(i) &= (this->Imat(i)*iv.mat(j)) * iv.rhs(j); } - if (this->rhs(i).is_empty()) { this->set_empty(); return *this; } + if (useful[j]) b=j; } - int b=iv.nbmat-1; - if (iv.nbmat==nbcommun || - (this->nbmat==2 && (ctcG || this->mats[1]==iv.mats[b]))) { + if (b==-1 || iv.nbmat==nbcommun || + (this->nbmat==2 && ctcG)) { this->simplify(); return (*this); } - if (this->nbmat==1 && this->mats[0]==iv.mats[b]) b=0; /* iv.nbmat=2 ! */ IntervalVector newV = iv.rhs(b); newV &= iv.Imat(b)*this->bbox(); for (unsigned int i=0;inbmat;i++) { newV &= (iv.Imat(b)*this->mat(i)) * this->rhs(i); } +// std::cout << "borrow base meet other=0\n"; this->borrow_base(iv,b,newV); if (newV.is_empty()) this->simplify(); @@ -744,6 +779,19 @@ namespace codac2 { } return this->bbox().is_subset(iv.bbox()); } + bool IParals::is_subset(const IParals& ip) const { + if (this->empty) return true; + if (ip.empty) return false; + if (!this->is_subset(ip.bbox())) return false; + for (unsigned int i=0;ibbox(); + for (unsigned int j=0;jnbmat;j++) { + r &= (ip.Imat(i)*this->mat(j))*this->rhs(j); + } + if (!r.is_subset(ip.rhs(i))) return false; + } + return true; + } IParals& IParals::toPointMatrices() { @@ -1035,7 +1083,11 @@ namespace codac2 { if (iv.empty || iv.nbmat==0) { return (str << iv.bbox()); } str << "IParals: box " << iv.bbox() << "\n"; for (unsigned int i=0;inbmat<2) { this->mats.push_back(iv.mats[n]); this->Vrhs.push_back(V); @@ -397,10 +414,12 @@ namespace codac2 { (unsigned int n,const IntervalMatrix &M, const IntervalMatrix &rM) { assert(n>=0 && nmats[n]->first; + IntervalMatrix newMi = this->mats[n]->second*rM; + punctualize_coupleMatrix(newM,newMi); std::shared_ptr> nPtr = std::make_shared> - (std::pair(M*this->mats[n]->first, - this->mats[n]->second*rM)); + (std::pair(newM,newMi)); this->mats[n] = nPtr; } inline void IParals::set_empty() { @@ -408,6 +427,12 @@ namespace codac2 { for (unsigned int i=0;i<=this->nbmat;i++) (this->Vrhs[i]).set_empty(); } + inline bool IParals::set_display_inv(bool newVal) { + bool old = IParals::display_inv; + IParals::display_inv=newVal; + return old; + } + bool set_display_inv(bool newVal); codac::TubeVector to_codac1(codac2::Tube& tube); } diff --git a/src/core/2/integration/codac2_expIMat.cpp b/src/core/2/integration/codac2_expIMat.cpp index c2935fd59..84f7d29a0 100644 --- a/src/core/2/integration/codac2_expIMat.cpp +++ b/src/core/2/integration/codac2_expIMat.cpp @@ -322,6 +322,7 @@ IntervalMatrix inv_IntervalMatrix(const IntervalMatrix& M, bool prec) { return Res; } + // 2) Res := A^(-1)Res ( empty if A contains a singular matrix) // (still crude Gaussian elimination using Rows) @@ -390,6 +391,44 @@ void inv_IntervalMatrix(const IntervalMatrix& A, } } +/* punctualize a couple "matrix" - "inverse of matrix" : + starting from ([A],[Ai]) such that [A].[Ai] ~ Id, compute + [A'] ... + technique used : Am = [A].mid Aim = [Ai].mid + then compute [A'] = (Am [Aim])^-1 Am + we use the algorithm of inv_IntervalMatrix with simplification + (assume Am [Aim] is almost Id */ +void punctualize_coupleMatrix(IntervalMatrix& A, + IntervalMatrix &Ai) { + Ai = Ai.mid(); + A = A.mid(); + IntervalMatrix M_copy(A*Ai); + const int order = A.nb_rows(); + assert(order == A.nb_cols()); + for (int col=0;col... + technique used : Am = [A].mid Aim = [Ai].mid + then compute [A'] = (Am [Aim])^-1 Am + we use the algorithm of inv_IntervalMatrix with simplification + (assume Am [Aim] is almost Id */ +void punctualize_coupleMatrix(IntervalMatrix& A, + IntervalMatrix &Ai); + // Determination of a rank a set of vectors A, with Af the "maximal" set // of free vectors, and E = pseudo-inverse (A Et = id) // if needed, we complete to get a square matrix with other vectors From f0de7ef7a565e397b752077fa520e40e8d850189 Mon Sep 17 00:00:00 2001 From: damien-masse Date: Fri, 12 Jan 2024 09:51:22 +0100 Subject: [PATCH 7/8] =?UTF-8?q?Premi=C3=A8re=20version=20avec=20intflpl,?= =?UTF-8?q?=20validation(=3F)=20avec=20IParals.?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- CMakeLists.txt | 12 + .../2/contractors/codac2_CtcDiffInclusion.cpp | 473 ---------------- .../2/contractors/codac2_CtcDiffInclusion.h | 515 +++++++++++++++++- src/core/2/integration/codac2_IParals.cpp | 28 +- src/core/2/integration/codac2_IParals.h | 4 +- src/core/2/integration/codac2_IntPoly.cpp | 72 +++ src/core/2/integration/codac2_IntPoly.h | 42 ++ src/core/CMakeLists.txt | 10 +- 8 files changed, 652 insertions(+), 504 deletions(-) delete mode 100644 src/core/2/contractors/codac2_CtcDiffInclusion.cpp create mode 100644 src/core/2/integration/codac2_IntPoly.cpp create mode 100644 src/core/2/integration/codac2_IntPoly.h diff --git a/CMakeLists.txt b/CMakeLists.txt index 3dbf9b2a2..4853d7663 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -108,6 +108,18 @@ add_definitions(${EIGEN3_DEFINITIONS}) include_directories(${EIGEN3_INCLUDE_DIRS}) +################################################################################ +# Looking for intflpl +################################################################################ + + option(WITH_INTFLPL "Using intfpl (floating-point polyhedra)" OFF) + +if(WITH_INTFLPL) + find_package(intflpl REQUIRED) + message(STATUS "Found intflpl version ${INTFLPL_VERSION}") + include_directories(${INTFLPL_INCLUDE_DIRS}) +endif() + ################################################################################ # Looking for CAPD (if needed) diff --git a/src/core/2/contractors/codac2_CtcDiffInclusion.cpp b/src/core/2/contractors/codac2_CtcDiffInclusion.cpp deleted file mode 100644 index 3a8992eb4..000000000 --- a/src/core/2/contractors/codac2_CtcDiffInclusion.cpp +++ /dev/null @@ -1,473 +0,0 @@ -/** - * CtcDiffInclusion class - * ---------------------------------------------------------------------------- - * \date 2022 - * \author - * \copyright Copyright 2022 Codac Team - * \license This program is distributed under the terms of - * the GNU Lesser General Public License (LGPL). - */ - -#include "codac2_CtcDiffInclusion.h" -#include "codac2_IParals.h" -#include "codac2_expIMat.h" - -using namespace std; -using namespace ibex; - -namespace codac2 -{ - CtcDiffInclusion::CtcDiffInclusion(const TFunction& f) - : _f(f) - { - - } - - const TFunction& CtcDiffInclusion::f() const - { - return _f; - } - - const IntervalVector CtcDiffInclusion::eval_function(const Interval &tim, - const IntervalVector& cdom, const IntervalVector *u) const { - if (u!=NULL) - return _f.eval_vector(tim,cdom,*u); else - return _f.eval_vector(tim,cdom); - - } - const IntervalVector CtcDiffInclusion::eval_function(const Interval &tim, - const IParals& cdom, const IntervalVector* u) const { - return this->eval_function(tim, cdom.box(),u); - } - const IntervalVector CtcDiffInclusion::eval_function(double tim, - const IntervalVector& cdom, const IntervalVector* u) const { - const Interval timI(tim); - return this->eval_function(timI,cdom,u); - } - const IntervalVector CtcDiffInclusion::eval_function(double tim, - const IParals& cdom, const IntervalVector* u) const { - const Interval timI(tim); - return this->eval_function(timI,cdom.box(),u); - } - const IntervalVector CtcDiffInclusion::eval_function(double tim, - const Vector& cdom, const IntervalVector* u) const { - const Interval timI(tim); - const IntervalVector cdomI(cdom); - return this->eval_function(timI,cdomI,u); - } - - - // basic enclosing for evolution - IParals CtcDiffInclusion::extend_box_basic(const IParals& frame, - const IParals& startIV, const IntervalVector* u, - const Interval& tim, - double inflation_factor, - TimePropag t_propa, - int nb_tries) const { - double tstep = (t_propa == TimePropag::FORWARD ? tim.diam() - : -tim.diam()); - Interval btime(0.0,tstep); - if (t_propa == TimePropag::BACKWARD) btime = -btime; - double starttime = (t_propa == TimePropag::FORWARD ? tim.lb() - : tim.ub()); - double endtime = (t_propa == TimePropag::FORWARD ? tim.ub() - : tim.lb()); - /* estimation des pentes */ - IntervalVector k1 = this->eval_function(starttime,startIV,u); - IParals XB1 = sum_tau(startIV,(tstep/2.0)*k1); - IntervalVector k2 = this->eval_function(tim.mid(), XB1,u); - XB1 = sum_tau(XB1,(tstep/2.0)* k2); - IntervalVector k3 = this->eval_function(endtime, XB1,u); - /* compute approximation of the result */ - IParals Res = sum_tau(startIV, - (tstep*0.505)* ((k1|k2) + (k2|k3))); - - Res&=frame; - - if (nb_tries==0) return Res; - /* inflation */ -// Res.inflate(inflation_factor); -// Res &= frame; -// if (nb_tries<3) { - double ifact = nb_tries==0 ? 1.0 : inflation_factor; - if (nb_tries==2) ifact *= inflation_factor; - Res = sum_tau(startIV, (ifact * btime) * - (this->eval_function(tim,Res,u))); - Res&=frame; - // } - - return Res; - } - - - IntervalMatrix CtcDiffInclusion::jacobian(const IParals& codom, - const IntervalVector* u, - const Interval& tdom, - IntervalVector& tvec) const { - const Function& func = _f.getFunction(); - const IntervalVector& b = codom.box(); - if (u!=NULL) { - IntervalVector box(b.size()+u->size()+1); - box[0]=tdom; - box.put(1,b); - box.put(1+b.size(),*u); - IntervalMatrix bjac = func.jacobian(box); - tvec = bjac.col(0); - return bjac.cols(1,b.size()); - } else { - IntervalVector box(b.size()+1); - box[0]=tdom; - box.put(1,b); - IntervalMatrix bjac = func.jacobian(box); - tvec = bjac.col(0); - return bjac.cols(1,b.size()); - } - } - - - bool CtcDiffInclusion::compute_step(const IParals& frame, - const IntervalVector* u, - const IParals& actState, - IParals& tauState, - IParals& finState, - const Interval& timeslice, - TimePropag t_propa) const { - - const int dim=frame.get_dim(); - IntervalVector tdiff(dim); - IntervalMatrix jacM = this->jacobian(tauState,u,timeslice,tdiff); - if (t_propa == TimePropag::BACKWARD) { - jacM = -jacM; - tdiff = -tdiff; - } - const Matrix jac = jacM.mid(); - jacM -= jac; - bool time_dependent = !tdiff.is_zero(); - - IntervalMatrix ExpM(dim,dim); - IntervalMatrix invExpM(dim,dim); - IntervalMatrix tauExpM(dim,dim); - IntervalMatrix IExpM(dim,dim); - IntervalMatrix tauIExpM(dim,dim); - IntervalMatrix VExpM(dim,dim); - IntervalMatrix tauVExpM(dim,dim); - Matrix IntAbs(dim,dim); - - double tsteps = timeslice.diam(); - - - global_exp(jac,tsteps,true,time_dependent, - ExpM,invExpM,tauExpM,IExpM,tauIExpM, - VExpM,tauVExpM,IntAbs); -// std::cout << "jacobian : " << jac << "\n"; -// std::cout << "ExpM : " << ExpM << "\n"; -// std::cout << "IExpM : " << IExpM << "\n"; -// std::cout << "IntAbs : " << IntAbs << "\n"; - /* other variables which needs to be kept */ - Vector cent_tdiff(dim); - Vector fun_evalc(dim); - Vector vuncert(dim); - Vector cent_tauState(dim); - - bool ok=false; - bool reducing=false; - int nb_red=0; - int nb_enl=0; - while(!ok) { - /* computing uncert and other values */ - cent_tauState = tauState.mid(); - IParals ctauState = tauState - cent_tauState; - IntervalVector uncert = jacM * ctauState; - if (time_dependent) { - cent_tdiff = tdiff.mid(); - IntervalVector ctdiff = tdiff-cent_tdiff; - uncert += timeslice.rad() * ctdiff; - } - IntervalVector fun_eval = - this->eval_function(timeslice.mid(),cent_tauState,u); - fun_evalc = fun_eval.mid(); - uncert += (fun_eval-fun_evalc); - if (t_propa == TimePropag::BACKWARD) { - fun_evalc = -fun_evalc; - } - -// cout << "uncert : " << uncert << "\n"; - /* now computing the new "tau" states */ - vuncert = IntAbs * uncert.ub(); - // ``derivation'' of the center ( f(cent) * int exp(M\tau) ) - IntervalVector tauCent = tauIExpM * fun_evalc; - // non-autonomous factor - if (time_dependent) { - tauCent += tauVExpM * cent_tdiff; - } - /* evolution of the center */ - for (int i=0;ijacobian(tauState,u,timeslice,tdiff); - if (t_propa == TimePropag::BACKWARD) { - jacM = -jacM; - tdiff = -tdiff; - } - jacM = jacM - jac; /* maybe not centered */ - } - } - /* computing ``arrival'' states */ - IntervalVector evolCenter = IExpM * fun_evalc; - // non-autonomous factor - if (time_dependent) { - evolCenter += VExpM * cent_tdiff; - } - for (int i=0;i& x, const Tube* u, TimePropag t_propa) - { - bool sameslicing=false; - // Verifying that x and u share exactly the same tdomain and slicing: - - if (u!=nullptr) { - if (x.tdomain() == u->tdomain()) sameslicing=true; - else { -// std::cout << x.tdomain()->t0_tf() << " " << u->tdomain()->t0_tf() << "\n"; - assert(x.tdomain()->t0_tf().is_subset(u->tdomain()->t0_tf())); - } - // Verifying that the provided tubes are consistent with the function -// std::cout << _f.nb_var() << " " << x.size() << " " << u->size() << "\n"; - assert((size_t)_f.nb_var() == x.size()+u->size()); - } else - assert((size_t)_f.nb_var() == x.size()); - assert((size_t)_f.image_dim() == x.size()); - - for(auto& sx : x) // sx is a Slice of the Tube x - { - if(sx.is_gate()) // the slice may be on a degenerated temporal domain, i.e. a gate - continue; - if (sx.tslice().t0_tf().is_unbounded()) continue; - - // su is a Slice of the Tube u: - const shared_ptr du = - (u == nullptr ? nullptr : - (sameslicing ? - std::make_shared( - std::static_pointer_cast>(sx.tslice().slices().at(u))->codomain()) : - std::make_shared(u->eval(sx.tslice().t0_tf())) - ) - ); -// std::cout << "uslice : " << *du << "\n"; - // const shared_ptr> su = (u==NULL ? NULL : &(sx.tslice().slices().at(u))); -// const double dt = sx.t0_tf().diam(); - -//nullptr - // sx.set(su.codomain()); - - // ... - contract(sx,du,t_propa); - -/* - if(t_propa & TimePropag::FORWARD) - { - contract(sx,su,t_progag); - // Computations related to forward propagation - // ... - } - - if(t_propa & TimePropag::BACKWARD) - { - contract(sx,su,t_progag); - // Computations related to backward propagation - // ... - } -*/ - } - } - - void CtcDiffInclusion::contract_from_slice(Tube& x, const Tube* u, std::shared_ptr>& gate, TimePropag t_propa) - { - bool sameslicing=false; - // Verifying that x and u share exactly the same tdomain and slicing: - - if (u!=nullptr) { - if (x.tdomain() == u->tdomain()) sameslicing=true; - else assert(x.tdomain()->t0_tf().is_subset(u->tdomain()->t0_tf())); - // Verifying that the provided tubes are consistent with the function - assert((size_t)_f.nb_var() == x.size()+u->size()); - } else - assert((size_t)_f.nb_var() == x.size()); - assert((size_t)_f.image_dim() == x.size()); - - if (t_propa & TimePropag::FORWARD) { - std::shared_ptr> sx = gate->next_slice_ptr(); - while (sx!=NULL) { - if (!sx->is_gate() && !sx->tslice().t0_tf().is_unbounded()) { - const shared_ptr du = - (u == nullptr ? nullptr : - (sameslicing ? - std::make_shared( - std::static_pointer_cast>(sx->tslice().slices().at(u))->codomain()) : - std::make_shared(u->eval(sx->tslice().t0_tf())) - ) - ); - contract(*sx,du,TimePropag::FORWARD); - } - sx = sx->next_slice_ptr(); - } - } - if (t_propa & TimePropag::BACKWARD) { - std::shared_ptr> sx = gate->prev_slice_ptr(); - while (sx!=NULL) { - if (!sx->is_gate() && !sx->tslice().t0_tf().is_unbounded()) { - const shared_ptr du = - (u == nullptr ? nullptr : - (sameslicing ? - std::make_shared( - std::static_pointer_cast>(sx->tslice().slices().at(u))->codomain()) : - std::make_shared(u->eval(sx->tslice().t0_tf())) - ) - ); - contract(*sx,du,TimePropag::BACKWARD); - } - sx = sx->prev_slice_ptr(); - } - } - } - - void CtcDiffInclusion::contract(Slice& x, - const Slice& u, TimePropag t_propa) { - const shared_ptr du = - std::make_shared(u.codomain()); - this->contract(x,du,t_propa); - } - - void CtcDiffInclusion::contract(Slice& x, - const std::shared_ptr& uDom, TimePropag t_propa) - { - if (uDom!=nullptr) { - // Verifying that the provided slices are consistent with the function - assert((size_t)_f.nb_var() == x.size()+uDom->size()); - } else - assert((size_t)_f.nb_var() == x.size()); - assert((size_t)_f.image_dim() == x.size()); - if (x.is_gate() || x.is_empty()) return; - - Interval timeslice=x.t0_tf(); -// double dt = timeslice.diam(); - - // ... - IParals frame = x.codomain(); - const IntervalVector* uVect = (uDom==nullptr ? nullptr : - &(uDom->box())); - if((t_propa & TimePropag::FORWARD) && !(x.input_gate().is_unbounded())) - { - IParals approxIV(frame); - bool okreduced=false; /* becomes true when we have a contraction */ - double inflation_f = default_inflation_factor; /* initial inflation factor */ - int nb_tries=0; - IParals g_t0 = x.input_gate(); - IParals g_t1(g_t0); /* start matrices */ - while (true) { - if (!okreduced) { - approxIV = this->extend_box_basic(frame, - g_t0,uVect,timeslice,inflation_f, - TimePropag::FORWARD,nb_tries); - } - bool safe = compute_step(frame,uVect, - g_t0,approxIV,g_t1,timeslice,TimePropag::FORWARD); - if (!safe) { - if (nb_tries>2) { - inflation_f *= 1.5; - } - } else { - okreduced=true; - if (nb_tries>=5) break; - } - nb_tries++; - } - frame &= approxIV; /* ensure contractance */ - if (x.next_slice_ptr() && x.next_slice_ptr()->is_gate()) { - IParals end = x.next_slice_ptr()->codomain(); - end &= g_t1; /* ensure contractance */ - x.next_slice_ptr()->set(end); - } - x.set(frame); /* first the gate, then the frame... */ - } - - if((t_propa & TimePropag::BACKWARD) && (!x.output_gate().is_unbounded())) - { - IParals approxIV(frame); - bool okreduced=false; /* becomes true when we have a contraction */ - double inflation_f = default_inflation_factor; /* initial inflation factor */ - int nb_tries=0; - IParals g_t1 = x.output_gate(); - IParals g_t0(g_t1); /* start matrices */ - while (true) { - if (!okreduced) { - approxIV = this->extend_box_basic(frame, - g_t1,uVect,timeslice,inflation_f, - TimePropag::BACKWARD,nb_tries); - } - bool safe = compute_step(frame,uVect, - g_t1,approxIV,g_t0,timeslice,TimePropag::BACKWARD); - if (!safe) { - if (nb_tries>2) { - inflation_f *= 1.5; - } - } else { - okreduced=true; - if (nb_tries>=5) break; - } - nb_tries++; - } - frame &= approxIV; /* ensure contractance */ - x.set(frame); - if (x.prev_slice_ptr() && x.prev_slice_ptr()->is_gate()) { - IParals beg = x.prev_slice_ptr()->codomain(); - beg &= g_t0; /* ensure contractance */ - x.prev_slice_ptr()->set(beg); - } - } - } -} diff --git a/src/core/2/contractors/codac2_CtcDiffInclusion.h b/src/core/2/contractors/codac2_CtcDiffInclusion.h index 0dedbe989..351495ada 100644 --- a/src/core/2/contractors/codac2_CtcDiffInclusion.h +++ b/src/core/2/contractors/codac2_CtcDiffInclusion.h @@ -27,6 +27,8 @@ namespace codac2 * \class CtcDiffInclusion * \brief ... */ + template + class CtcDiffInclusion { public: @@ -38,46 +40,46 @@ namespace codac2 const codac::IntervalVector eval_function(const Interval &tim, const codac::IntervalVector& cdom, const codac::IntervalVector* u) const; const codac::IntervalVector eval_function(const Interval &tim, - const IParals& cdom, const codac::IntervalVector* u) const; + const _TpP& cdom, const codac::IntervalVector* u) const; /** evaluate function on box + time */ const codac::IntervalVector eval_function(double tim, const codac::IntervalVector& cdom, const codac::IntervalVector* u) const; const codac::IntervalVector eval_function(double tim, - const IParals& cdom, const codac::IntervalVector* u) const; + const _TpP& cdom, const codac::IntervalVector* u) const; /** evaluate function on point + time */ const codac::IntervalVector eval_function(double tim, const Vector& cdom, const codac::IntervalVector* u) const; - IntervalMatrix jacobian(const IParals& codom, + IntervalMatrix jacobian(const _TpP& codom, const codac::IntervalVector* u, const Interval& tdom, codac::IntervalVector& tvec) const; - IParals extend_box_basic(const IParals& frame, - const IParals& startIV, + _TpP extend_box_basic(const _TpP& frame, + const _TpP& startIV, const codac::IntervalVector* u, const Interval& tim, double inflation_factor, TimePropag t_propa, int nb_tries) const; - bool compute_step(const IParals& frame, + bool compute_step(const _TpP& frame, const codac::IntervalVector* u, - const IParals& actState, - IParals& tauState, - IParals& finState, + const _TpP& actState, + _TpP& tauState, + _TpP& finState, const Interval& timeslice, TimePropag t_propa) const; - void contract(Tube& x, const Tube* u, TimePropag t_propa = TimePropag::FORWARD | TimePropag::BACKWARD); - void contract_from_slice(Tube& x, const Tube* u, std::shared_ptr>& gate, TimePropag t_propa = TimePropag::FORWARD | TimePropag::BACKWARD); - void contract(Slice& x, const Slice& u, TimePropag t_propa = TimePropag::FORWARD | TimePropag::BACKWARD); - void contract(Slice& x, const std::shared_ptr& uDom, TimePropag t_propa = TimePropag::FORWARD | TimePropag::BACKWARD); + void contract(Tube<_TpP>& x, const Tube<_TpP>* u, TimePropag t_propa = TimePropag::FORWARD | TimePropag::BACKWARD); + void contract_from_slice(Tube<_TpP>& x, const Tube<_TpP>* u, std::shared_ptr>& gate, TimePropag t_propa = TimePropag::FORWARD | TimePropag::BACKWARD); + void contract(Slice<_TpP>& x, const Slice<_TpP>& u, TimePropag t_propa = TimePropag::FORWARD | TimePropag::BACKWARD); + void contract(Slice<_TpP>& x, const std::shared_ptr& uDom, TimePropag t_propa = TimePropag::FORWARD | TimePropag::BACKWARD); const TFunction& f() const; @@ -87,11 +89,496 @@ namespace codac2 const TFunction _f; /** inflation factor in extend_box_basic */ - constexpr static double default_inflation_factor = 1.01; + constexpr static double default_inflation_factor = 1.1; constexpr static int narrow_factor = 3; }; } + +/** + * CtcDiffInclusion code + */ + +#include "codac2_CtcDiffInclusion.h" + +using namespace std; +using namespace ibex; + +namespace codac2 +{ + templace + CtcDiffInclusion<_TpP>::CtcDiffInclusion(const TFunction& f) + : _f(f) + { + + } + + templace + const TFunction& CtcDiffInclusion<_TpP>::f() const + { + return _f; + } + + templace + const IntervalVector CtcDiffInclusion<_TpP>::eval_function(const Interval &tim, + const IntervalVector& cdom, const IntervalVector *u) const { + if (u!=NULL) + return _f.eval_vector(tim,cdom,*u); else + return _f.eval_vector(tim,cdom); + + } + templace + const IntervalVector CtcDiffInclusion<_TpP>::eval_function(const Interval &tim, + const _TpP& cdom, const IntervalVector* u) const { + return this->eval_function(tim, cdom.box(),u); + } + templace + const IntervalVector CtcDiffInclusion<_TpP>::eval_function(double tim, + const IntervalVector& cdom, const IntervalVector* u) const { + const Interval timI(tim); + return this->eval_function(timI,cdom,u); + } + templace + const IntervalVector CtcDiffInclusion<_TpP>::eval_function(double tim, + const _TpP& cdom, const IntervalVector* u) const { + const Interval timI(tim); + return this->eval_function(timI,cdom.box(),u); + } + templace + const IntervalVector CtcDiffInclusion<_TpP>::eval_function(double tim, + const Vector& cdom, const IntervalVector* u) const { + const Interval timI(tim); + const IntervalVector cdomI(cdom); + return this->eval_function(timI,cdomI,u); + } + + + // basic enclosing for evolution + templace + _TpP CtcDiffInclusion<_TpP>::extend_box_basic(const _TpP& frame, + const _TpP& startIV, const IntervalVector* u, + const Interval& tim, + double inflation_factor, + TimePropag t_propa, + int nb_tries) const { + double tstep = (t_propa == TimePropag::FORWARD ? tim.diam() + : -tim.diam()); + Interval btime(0.0,tstep); + if (t_propa == TimePropag::BACKWARD) btime = -btime; + double starttime = (t_propa == TimePropag::FORWARD ? tim.lb() + : tim.ub()); + double endtime = (t_propa == TimePropag::FORWARD ? tim.ub() + : tim.lb()); + /* estimation des pentes */ + IntervalVector k1 = this->eval_function(starttime,startIV,u); + _TpP XB1 = sum_tau(startIV,(tstep/2.0)*k1); + IntervalVector k2 = this->eval_function(tim.mid(), XB1,u); + XB1 = sum_tau(XB1,(tstep/2.0)* k2); + IntervalVector k3 = this->eval_function(endtime, XB1,u); + /* compute approximation of the result */ + _TpP Res = sum_tau(startIV, + (tstep*0.505)* ((k1|k2) + (k2|k3))); + + Res&=frame; + + if (nb_tries==0) return Res; + /* inflation */ +// Res.inflate(inflation_factor); +// Res &= frame; +// if (nb_tries<3) { + double ifact = nb_tries==0 ? 1.0 : inflation_factor; + if (nb_tries==2) ifact *= inflation_factor; + Res = sum_tau(startIV, (ifact * btime) * + (this->eval_function(tim,Res,u))); + Res&=frame; + // } + + return Res; + } + + + templace + IntervalMatrix CtcDiffInclusion<_TpP>::jacobian(const _TpP& codom, + const IntervalVector* u, + const Interval& tdom, + IntervalVector& tvec) const { + const Function& func = _f.getFunction(); + const IntervalVector& b = codom.box(); + if (u!=NULL) { + IntervalVector box(b.size()+u->size()+1); + box[0]=tdom; + box.put(1,b); + box.put(1+b.size(),*u); + IntervalMatrix bjac = func.jacobian(box); + tvec = bjac.col(0); + return bjac.cols(1,b.size()); + } else { + IntervalVector box(b.size()+1); + box[0]=tdom; + box.put(1,b); + IntervalMatrix bjac = func.jacobian(box); + tvec = bjac.col(0); + return bjac.cols(1,b.size()); + } + } + + + templace + bool CtcDiffInclusion<_TpP>::compute_step(const _TpP& frame, + const IntervalVector* u, + const _TpP& actState, + _TpP& tauState, + _TpP& finState, + const Interval& timeslice, + TimePropag t_propa) const { + + const int dim=frame.get_dim(); + IntervalVector tdiff(dim); + IntervalMatrix jacM = this->jacobian(tauState,u,timeslice,tdiff); + if (t_propa == TimePropag::BACKWARD) { + jacM = -jacM; + tdiff = -tdiff; + } + const Matrix jac = jacM.mid(); + jacM -= jac; + bool time_dependent = !tdiff.is_zero(); + + IntervalMatrix ExpM(dim,dim); + IntervalMatrix invExpM(dim,dim); + IntervalMatrix tauExpM(dim,dim); + IntervalMatrix IExpM(dim,dim); + IntervalMatrix tauIExpM(dim,dim); + IntervalMatrix VExpM(dim,dim); + IntervalMatrix tauVExpM(dim,dim); + Matrix IntAbs(dim,dim); + + double tsteps = timeslice.diam(); + + + global_exp(jac,tsteps,true,time_dependent, + ExpM,invExpM,tauExpM,IExpM,tauIExpM, + VExpM,tauVExpM,IntAbs); +// std::cout << "jacobian : " << jac << "\n"; +// std::cout << "ExpM : " << ExpM << "\n"; +// std::cout << "IExpM : " << IExpM << "\n"; +// std::cout << "IntAbs : " << IntAbs << "\n"; + /* other variables which needs to be kept */ + Vector cent_tdiff(dim); + Vector fun_evalc(dim); + Vector vuncert(dim); + Vector cent_tauState(dim); + + bool ok=false; + bool reducing=false; + int nb_red=0; + int nb_enl=0; + while(!ok) { + /* computing uncert and other values */ + cent_tauState = tauState.mid(); + _TpP ctauState = tauState - cent_tauState; + IntervalVector uncert = jacM * ctauState; + if (time_dependent) { + cent_tdiff = tdiff.mid(); + IntervalVector ctdiff = tdiff-cent_tdiff; + uncert += timeslice.rad() * ctdiff; + } + IntervalVector fun_eval = + this->eval_function(timeslice.mid(),cent_tauState,u); + fun_evalc = fun_eval.mid(); + uncert += (fun_eval-fun_evalc); + if (t_propa == TimePropag::BACKWARD) { + fun_evalc = -fun_evalc; + } + +// cout << "uncert : " << uncert << "\n"; + /* now computing the new "tau" states */ + vuncert = IntAbs * uncert.ub(); + // ``derivation'' of the center ( f(cent) * int exp(M\tau) ) + IntervalVector tauCent = tauIExpM * fun_evalc; + // non-autonomous factor + if (time_dependent) { + tauCent += tauVExpM * cent_tdiff; + } + /* evolution of the center */ + for (int i=0;ijacobian(tauState,u,timeslice,tdiff); + if (t_propa == TimePropag::BACKWARD) { + jacM = -jacM; + tdiff = -tdiff; + } + jacM = jacM - jac; /* maybe not centered */ + } + } + /* computing ``arrival'' states */ + IntervalVector evolCenter = IExpM * fun_evalc; + // non-autonomous factor + if (time_dependent) { + evolCenter += VExpM * cent_tdiff; + } + for (int i=0;i + void CtcDiffInclusion<_TpP>::contract(Slice<_TpP>& x, + const std::shared_ptr& uDom, TimePropag t_propa) + { + if (uDom!=nullptr) { + // Verifying that the provided slices are consistent with the function + assert((size_t)_f.nb_var() == x.size()+uDom->size()); + } else + assert((size_t)_f.nb_var() == x.size()); + assert((size_t)_f.image_dim() == x.size()); + if (x.is_gate() || x.is_empty()) return; + + Interval timeslice=x.t0_tf(); +// double dt = timeslice.diam(); + + // ... + _TpP frame = x.codomain(); + const IntervalVector* uVect = (uDom==nullptr ? nullptr : + &(uDom->box())); + if((t_propa & TimePropag::FORWARD) && !(x.input_gate().is_unbounded())) + { + _TpP approxIV(frame); + bool okreduced=false; /* becomes true when we have a contraction */ + double inflation_f = default_inflation_factor; /* initial inflation factor */ + int nb_tries=0; + _TpP g_t0 = x.input_gate(); + _TpP g_t1(g_t0); /* start matrices */ + while (true) { + if (!okreduced) { + approxIV = this->extend_box_basic(frame, + g_t0,uVect,timeslice,inflation_f, + TimePropag::FORWARD,nb_tries); + } + bool safe = compute_step(frame,uVect, + g_t0,approxIV,g_t1,timeslice,TimePropag::FORWARD); + if (!safe) { + if (nb_tries>2) { + inflation_f *= 1.5; + } + } else { + okreduced=true; + if (nb_tries>=5) break; + } + nb_tries++; + } + frame &= approxIV; /* ensure contractance */ + if (x.next_slice_ptr() && x.next_slice_ptr()->is_gate()) { + _TpP end = x.next_slice_ptr()->codomain(); + end &= g_t1; /* ensure contractance */ + x.next_slice_ptr()->set(end); + } + x.set(frame); /* first the gate, then the frame... */ + } + + if((t_propa & TimePropag::BACKWARD) && (!x.output_gate().is_unbounded())) + { + _TpP approxIV(frame); + bool okreduced=false; /* becomes true when we have a contraction */ + double inflation_f = default_inflation_factor; /* initial inflation factor */ + int nb_tries=0; + _TpP g_t1 = x.output_gate(); + _TpP g_t0(g_t1); /* start matrices */ + while (true) { + if (!okreduced) { + approxIV = this->extend_box_basic(frame, + g_t1,uVect,timeslice,inflation_f, + TimePropag::BACKWARD,nb_tries); + } + bool safe = compute_step(frame,uVect, + g_t1,approxIV,g_t0,timeslice,TimePropag::BACKWARD); + if (!safe) { + if (nb_tries>2) { + inflation_f *= 1.5; + } + } else { + okreduced=true; + if (nb_tries>=5) break; + } + nb_tries++; + } + frame &= approxIV; /* ensure contractance */ + x.set(frame); + if (x.prev_slice_ptr() && x.prev_slice_ptr()->is_gate()) { + _TpP beg = x.prev_slice_ptr()->codomain(); + beg &= g_t0; /* ensure contractance */ + x.prev_slice_ptr()->set(beg); + } + } + } + + templace + void CtcDiffInclusion<_TpP>::contract_from_slice(Tube<_TpP>& x, const Tube<_TpP>* u, std::shared_ptr>& gate, TimePropag t_propa) + { + bool sameslicing=false; + // Verifying that x and u share exactly the same tdomain and slicing: + + if (u!=nullptr) { + if (x.tdomain() == u->tdomain()) sameslicing=true; + else assert(x.tdomain()->t0_tf().is_subset(u->tdomain()->t0_tf())); + // Verifying that the provided tubes are consistent with the function + assert((size_t)_f.nb_var() == x.size()+u->size()); + } else + assert((size_t)_f.nb_var() == x.size()); + assert((size_t)_f.image_dim() == x.size()); + + if (t_propa & TimePropag::FORWARD) { + std::shared_ptr> sx = gate->next_slice_ptr(); + while (sx!=NULL) { + if (!sx->is_gate() && !sx->tslice().t0_tf().is_unbounded()) { + const shared_ptr du = + (u == nullptr ? nullptr : + (sameslicing ? + std::make_shared( + std::static_pointer_cast>(sx->tslice().slices().at(u))->codomain()) : + std::make_shared(u->eval(sx->tslice().t0_tf())) + ) + ); + contract(*sx,du,TimePropag::FORWARD); + } + sx = sx->next_slice_ptr(); + } + } + if (t_propa & TimePropag::BACKWARD) { + std::shared_ptr> sx = gate->prev_slice_ptr(); + while (sx!=NULL) { + if (!sx->is_gate() && !sx->tslice().t0_tf().is_unbounded()) { + const shared_ptr du = + (u == nullptr ? nullptr : + (sameslicing ? + std::make_shared( + std::static_pointer_cast>(sx->tslice().slices().at(u))->codomain()) : + std::make_shared(u->eval(sx->tslice().t0_tf())) + ) + ); + contract(*sx,du,TimePropag::BACKWARD); + } + sx = sx->prev_slice_ptr(); + } + } + } + + templace + void CtcDiffInclusion<_TpP>::contract(Slice<_TpP>& x, + const Slice<_TpP>& u, TimePropag t_propa) { + const shared_ptr du = + std::make_shared(u.codomain()); + this->contract(x,du,t_propa); + } + + templace + void CtcDiffInclusion<_TpP>::contract(Tube<_TpP>& x, const Tube<_TpP>* u, TimePropag t_propa) + { + bool sameslicing=false; + // Verifying that x and u share exactly the same tdomain and slicing: + + if (u!=nullptr) { + if (x.tdomain() == u->tdomain()) sameslicing=true; + else { +// std::cout << x.tdomain()->t0_tf() << " " << u->tdomain()->t0_tf() << "\n"; + assert(x.tdomain()->t0_tf().is_subset(u->tdomain()->t0_tf())); + } + // Verifying that the provided tubes are consistent with the function +// std::cout << _f.nb_var() << " " << x.size() << " " << u->size() << "\n"; + assert((size_t)_f.nb_var() == x.size()+u->size()); + } else + assert((size_t)_f.nb_var() == x.size()); + assert((size_t)_f.image_dim() == x.size()); + + for(auto& sx : x) // sx is a Slice<_TpP> of the Tube<_TpP> x + { + if(sx.is_gate()) // the slice may be on a degenerated temporal domain, i.e. a gate + continue; + if (sx.tslice().t0_tf().is_unbounded()) continue; + + // su is a Slice<_TpP> of the Tube<_TpP> u: + const shared_ptr du = + (u == nullptr ? nullptr : + (sameslicing ? + std::make_shared( + std::static_pointer_cast>(sx.tslice().slices().at(u))->codomain()) : + std::make_shared(u->eval(sx.tslice().t0_tf())) + ) + ); +// std::cout << "uslice : " << *du << "\n"; + // const shared_ptr> su = (u==NULL ? NULL : &(sx.tslice().slices().at(u))); +// const double dt = sx.t0_tf().diam(); + +//nullptr + // sx.set(su.codomain()); + + // ... + contract(sx,du,t_propa); + +/* + if(t_propa & TimePropag::FORWARD) + { + contract(sx,su,t_progag); + // Computations related to forward propagation + // ... + } + + if(t_propa & TimePropag::BACKWARD) + { + contract(sx,su,t_progag); + // Computations related to backward propagation + // ... + } +*/ + } + } +} + #endif diff --git a/src/core/2/integration/codac2_IParals.cpp b/src/core/2/integration/codac2_IParals.cpp index 1292c73f9..cff320fdb 100644 --- a/src/core/2/integration/codac2_IParals.cpp +++ b/src/core/2/integration/codac2_IParals.cpp @@ -224,7 +224,7 @@ namespace codac2 { this->bbox() = x; return *this; } - assert (dim=x.size()); + assert (dim==x.size()); if (x.is_empty()) { this->set_empty(); return *this; } this->empty=false; this->bbox() = x; @@ -550,7 +550,7 @@ namespace codac2 { /** union with a box */ IParals& IParals::operator|= (const IntervalVector& x) { - std::cout << "|=" << *this << " " << x << "\n"; +// std::cout << "|=" << *this << " " << x << "\n"; assert(dim>0); if (x.is_empty()) return *this; if (this->empty) { @@ -1040,36 +1040,36 @@ namespace codac2 { /** generate a list of (2D) points, the convex hull of which is an * (over)-approximation of the projection of the polyhedron */ - ConvexPolygon IParals::over_polygon(const Matrix& M) const { + ConvexPolygon over_polygon(const IParals &ip, const Matrix& M) { /* first we generate a projection of the parallelotope */ - if (this->empty) return ConvexPolygon(); + if (ip.empty) return ConvexPolygon(); /* just the first polygon (not manage intersection) */ - Vector V1(this->dim); + Vector V1(ip.dim); ConvexPolygon res; /* compute the projection for large dimension is a bit complex (but interesting), will do it dirty */ - for (unsigned int k=0;k<=this->nbmat;k++) { - bool val[this->dim]; + for (unsigned int k=0;k<=ip.nbmat;k++) { + bool val[ip.dim]; std::vector lpoints; - for (int i=0;idim;i++) { + for (int i=0;irhs(k)[i].lb(); + V1[i] = ip.rhs(k)[i].lb(); } while (true) { - if (knbmat) { - lpoints.push_back(codac::ThickPoint(M*(this->mat(k)*V1))); + if (k=0 && val[j]==true) { - V1[j]=this->rhs(k)[j].lb(); + V1[j]=ip.rhs(k)[j].lb(); val[j]=false; j--; } if (j<0) break; val[j]=true; - V1[j] = this->rhs(k)[j].ub(); + V1[j] = ip.rhs(k)[j].ub(); } ConvexPolygon a(lpoints); if (k==0) res=a; else res= res & a; diff --git a/src/core/2/integration/codac2_IParals.h b/src/core/2/integration/codac2_IParals.h index 6795078a9..3d5298ec9 100644 --- a/src/core/2/integration/codac2_IParals.h +++ b/src/core/2/integration/codac2_IParals.h @@ -39,6 +39,7 @@ namespace codac2 IParals sum_tau(const IParals& iv, const IntervalVector& V, bool keep=false); + ConvexPolygon over_polygon(const IParals &ip, const Matrix& M); class IParals { @@ -303,7 +304,7 @@ namespace codac2 /** generate a ConvexPolygon which is an overapproximation of the * projection of the polyhedron (basic form) */ - ConvexPolygon over_polygon(const Matrix& M) const; + friend ConvexPolygon over_polygon(const IParals &ip, const Matrix& M); /** display */ @@ -435,6 +436,7 @@ namespace codac2 { bool set_display_inv(bool newVal); codac::TubeVector to_codac1(codac2::Tube& tube); + } #endif diff --git a/src/core/2/integration/codac2_IntPoly.cpp b/src/core/2/integration/codac2_IntPoly.cpp new file mode 100644 index 000000000..b589553ac --- /dev/null +++ b/src/core/2/integration/codac2_IntPoly.cpp @@ -0,0 +1,72 @@ +/** + * \file representation of polyhedron as inter M_i X_i + * ---------------------------------------------------------------------------- + * \date 2022 + * \author Damien Massé + * \copyright Copyright 2022 Codac Team + * \license This program is distributed under the terms of + * the GNU Lesser General Public License (LGPL). + */ + + +#include +#include +#include +#include +#include +#include +#include +#include "codac_Interval.h" +#include "codac_IntervalVector.h" +#include "codac_IntervalMatrix.h" +#include "codac2_expIMat.h" +#include "codac_polygon_arithmetic.h" +#include "codac_ConvexPolygon.h" +#include "codac2_IntPoly.h" + +using namespace intflpl; +using namespace codac; + +namespace codac2 { + + /** generate a ConvexPolygon which is an overapproximation of the + * projection of the polyhedron (basic form) + */ + codac::ConvexPolygon over_polygon(const IntPoly &ip, const Matrix& M) { + if (ip.get_dim()==2) { + std::vector lx; + std::vector ly; + ip.vertices2D(lx,ly); +// std::cout << lx.size() << " " << ly.size() << "\n"; + std::vector vert; + for (unsigned int i=0;i=3 */ + return codac::ConvexPolygon(M*ip.box()); + } + } + + + codac::TubeVector to_codac1(codac2::Tube& tube) + { + codac::TubeVector x(tube.t0_tf(), tube.size()); + for(const auto& s : tube) + if(!s.t0_tf().is_unbounded()) + x.set(s.codomain().box(), s.t0_tf()); + for(const auto& s : tube) // setting gate (were overwritten) + if(s.t0_tf().is_degenerated()) + x.set(s.codomain().box(), s.t0_tf()); + return x; + } + + +} + diff --git a/src/core/2/integration/codac2_IntPoly.h b/src/core/2/integration/codac2_IntPoly.h new file mode 100644 index 000000000..d4e736e8d --- /dev/null +++ b/src/core/2/integration/codac2_IntPoly.h @@ -0,0 +1,42 @@ +/** + * \file representation of polyhedron as inter M_i X_i + * + * ---------------------------------------------------------------------------- + * \date 2022 + * \author Damien Massé + * \copyright Copyright 2022 Codac Team + * \license This program is distributed under the terms of + * the GNU Lesser General Public License (LGPL). + */ + +#ifndef __CODAC2_INTPOLY_H__ +#define __CODAC2_INTPOLY_H__ + +#include +#include +#include +#include "codac_TubeVector.h" +#include "codac_Interval.h" +#include "codac_IntervalVector.h" +#include "codac_IntervalMatrix.h" +#include "codac_ConvexPolygon.h" +#include "codac2_Tube.h" +#include "codac2_expIMat.h" +#include + +using namespace intflpl; + + +namespace codac2 +{ + + /** generate a ConvexPolygon which is an overapproximation of the + * projection of the polyhedron (basic form) + */ + codac::ConvexPolygon over_polygon(const IntPoly &ip, const Matrix& M); + + codac::TubeVector to_codac1(codac2::Tube& tube); + +} + +#endif diff --git a/src/core/CMakeLists.txt b/src/core/CMakeLists.txt index 43d93f014..531202cb6 100644 --- a/src/core/CMakeLists.txt +++ b/src/core/CMakeLists.txt @@ -199,7 +199,6 @@ ${CMAKE_CURRENT_SOURCE_DIR}/2/domains/tube/codac2_TubeComponent.h ${CMAKE_CURRENT_SOURCE_DIR}/2/domains/tube/codac2_TubeEvaluation.h ${CMAKE_CURRENT_SOURCE_DIR}/2/domains/paving/codac2_Paving.h - ${CMAKE_CURRENT_SOURCE_DIR}/2/contractors/codac2_CtcDiffInclusion.cpp ${CMAKE_CURRENT_SOURCE_DIR}/2/contractors/codac2_CtcDiffInclusion.h ${CMAKE_CURRENT_SOURCE_DIR}/2/contractors/codac2_CtcLinobs.cpp ${CMAKE_CURRENT_SOURCE_DIR}/2/contractors/codac2_CtcLinobs.h @@ -218,6 +217,13 @@ ${CMAKE_CURRENT_SOURCE_DIR}/2/contractors/codac2_CtcAction.h ) +if(WITH_INTFLPL) + list(APPEND CODAC2_SRC # adding utility functions for intflpl::IntPoly + ${CMAKE_CURRENT_SOURCE_DIR}/2/integration/codac2_IntPoly.cpp + ${CMAKE_CURRENT_SOURCE_DIR}/2/integration/codac2_IntPoly.h + ) +endif() + set(SRC ${CODAC1_SRC} ${CODAC2_SRC}) @@ -311,4 +317,4 @@ install(FILES ${CODAC_HDR} DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/codac) install(FILES ${CODAC_MAIN_HEADER} DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/codac) install(FILES ${CODAC2_HDR} DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/codac) - install(FILES ${CODAC2_MAIN_HEADER} DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/codac) \ No newline at end of file + install(FILES ${CODAC2_MAIN_HEADER} DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/codac) From 5e02e36074aaf73a1aa24bcf95ce96ea44cba331 Mon Sep 17 00:00:00 2001 From: damien-masse Date: Fri, 12 Jan 2024 10:14:36 +0100 Subject: [PATCH 8/8] Bugfix (compile error) --- .../2/contractors/codac2_CtcDiffInclusion.h | 36 +++++++++---------- 1 file changed, 18 insertions(+), 18 deletions(-) diff --git a/src/core/2/contractors/codac2_CtcDiffInclusion.h b/src/core/2/contractors/codac2_CtcDiffInclusion.h index 351495ada..d99553fcd 100644 --- a/src/core/2/contractors/codac2_CtcDiffInclusion.h +++ b/src/core/2/contractors/codac2_CtcDiffInclusion.h @@ -108,20 +108,20 @@ using namespace ibex; namespace codac2 { - templace + template CtcDiffInclusion<_TpP>::CtcDiffInclusion(const TFunction& f) : _f(f) { } - templace + template const TFunction& CtcDiffInclusion<_TpP>::f() const { return _f; } - templace + template const IntervalVector CtcDiffInclusion<_TpP>::eval_function(const Interval &tim, const IntervalVector& cdom, const IntervalVector *u) const { if (u!=NULL) @@ -129,24 +129,24 @@ namespace codac2 return _f.eval_vector(tim,cdom); } - templace + template const IntervalVector CtcDiffInclusion<_TpP>::eval_function(const Interval &tim, const _TpP& cdom, const IntervalVector* u) const { return this->eval_function(tim, cdom.box(),u); } - templace + template const IntervalVector CtcDiffInclusion<_TpP>::eval_function(double tim, const IntervalVector& cdom, const IntervalVector* u) const { const Interval timI(tim); return this->eval_function(timI,cdom,u); } - templace + template const IntervalVector CtcDiffInclusion<_TpP>::eval_function(double tim, const _TpP& cdom, const IntervalVector* u) const { const Interval timI(tim); return this->eval_function(timI,cdom.box(),u); } - templace + template const IntervalVector CtcDiffInclusion<_TpP>::eval_function(double tim, const Vector& cdom, const IntervalVector* u) const { const Interval timI(tim); @@ -156,7 +156,7 @@ namespace codac2 // basic enclosing for evolution - templace + template _TpP CtcDiffInclusion<_TpP>::extend_box_basic(const _TpP& frame, const _TpP& startIV, const IntervalVector* u, const Interval& tim, @@ -199,7 +199,7 @@ namespace codac2 } - templace + template IntervalMatrix CtcDiffInclusion<_TpP>::jacobian(const _TpP& codom, const IntervalVector* u, const Interval& tdom, @@ -225,7 +225,7 @@ namespace codac2 } - templace + template bool CtcDiffInclusion<_TpP>::compute_step(const _TpP& frame, const IntervalVector* u, const _TpP& actState, @@ -306,12 +306,12 @@ namespace codac2 /* next tauState, do NOT change the underlying matrices */ _TpP ntauState(actState); - cout << "avant ctau" << ntauState << "\ncent " << cent_tauState << - "\ntauExpM " << tauExpM << "\ntauCent " << tauCent << "\n"; +// cout << "avant ctau" << ntauState << "\ncent " << cent_tauState << +// "\ntauExpM " << tauExpM << "\ntauCent " << tauCent << "\n"; ntauState.ctau_mult_and_add(cent_tauState,tauExpM,tauCent); - cout << "ntauState" << ntauState << "\ntauState " << tauState << - "\n"; +// cout << "ntauState" << ntauState << "\ntauState " << tauState << +// "\n"; ntauState &= frame; @@ -369,7 +369,7 @@ namespace codac2 - templace + template void CtcDiffInclusion<_TpP>::contract(Slice<_TpP>& x, const std::shared_ptr& uDom, TimePropag t_propa) { @@ -459,7 +459,7 @@ namespace codac2 } } - templace + template void CtcDiffInclusion<_TpP>::contract_from_slice(Tube<_TpP>& x, const Tube<_TpP>* u, std::shared_ptr>& gate, TimePropag t_propa) { bool sameslicing=false; @@ -510,7 +510,7 @@ namespace codac2 } } - templace + template void CtcDiffInclusion<_TpP>::contract(Slice<_TpP>& x, const Slice<_TpP>& u, TimePropag t_propa) { const shared_ptr du = @@ -518,7 +518,7 @@ namespace codac2 this->contract(x,du,t_propa); } - templace + template void CtcDiffInclusion<_TpP>::contract(Tube<_TpP>& x, const Tube<_TpP>* u, TimePropag t_propa) { bool sameslicing=false;