From 2c07462a82a4fc0e8fc7da171609a3144165a5c7 Mon Sep 17 00:00:00 2001 From: damien-masse Date: Thu, 22 Oct 2026 21:38:24 +0200 Subject: [PATCH] mix? --- examples/02_centered_form/main.cpp | 4 +- .../core/contractors/codac2_py_CtcInverse.h | 18 +-- src/core/actions/codac2_OctaSym_operator.h | 20 +-- src/core/contractors/codac2_CtcInverse.h | 49 ++++---- src/core/contractors/codac2_CtcInverseNotIn.h | 12 +- src/core/domains/tube/codac2_Tube_operator.h | 8 +- .../functions/analytic/codac2_AnalyticExpr.h | 17 +-- .../analytic/codac2_AnalyticFunction.h | 119 +++++++++++++++--- .../analytic/codac2_AnalyticFunction_impl.h | 50 ++++++-- .../functions/analytic/codac2_AnalyticType.h | 76 ++++++++++- .../analytic/codac2_analytic_constants.h | 8 +- .../analytic/codac2_analytic_variables.h | 2 +- src/core/operators/codac2_abs.h | 19 +-- src/core/operators/codac2_acos.h | 22 ++-- src/core/operators/codac2_arith_add.h | 54 ++++---- src/core/operators/codac2_arith_div.h | 65 ++++++---- src/core/operators/codac2_arith_mul.h | 101 +++++++++------ src/core/operators/codac2_arith_sub.h | 56 +++++---- src/core/operators/codac2_asin.h | 22 ++-- src/core/operators/codac2_atan.h | 14 ++- src/core/operators/codac2_atan2.h | 20 +-- src/core/operators/codac2_ceil.h | 8 +- src/core/operators/codac2_chi.h | 25 ++-- src/core/operators/codac2_component.h | 51 ++++---- src/core/operators/codac2_cos.h | 14 ++- src/core/operators/codac2_cosh.h | 14 ++- src/core/operators/codac2_cross_prod.h | 38 ++++-- src/core/operators/codac2_det.h | 89 +++++++------ src/core/operators/codac2_exp.h | 14 ++- src/core/operators/codac2_extend.h | 10 +- src/core/operators/codac2_flatten.h | 8 +- src/core/operators/codac2_floor.h | 8 +- src/core/operators/codac2_log.h | 22 ++-- src/core/operators/codac2_mat.h | 34 +++-- src/core/operators/codac2_max.h | 18 +-- src/core/operators/codac2_min.h | 17 +-- src/core/operators/codac2_mod.h | 8 +- src/core/operators/codac2_pow.h | 39 +++--- src/core/operators/codac2_sign.h | 12 +- src/core/operators/codac2_sin.h | 14 ++- src/core/operators/codac2_sinh.h | 14 ++- src/core/operators/codac2_sqr.h | 14 ++- src/core/operators/codac2_sqrt.h | 22 ++-- src/core/operators/codac2_subvector.h | 20 +-- src/core/operators/codac2_tan.h | 18 +-- src/core/operators/codac2_tanh.h | 14 ++- src/core/operators/codac2_transpose.h | 14 +-- src/core/operators/codac2_vec.h | 30 +++-- src/core/trajectory/codac2_Traj_operator.h | 14 ++- .../codac2_tests_AnalyticFunction.cpp | 36 ++++-- 50 files changed, 881 insertions(+), 514 deletions(-) diff --git a/examples/02_centered_form/main.cpp b/examples/02_centered_form/main.cpp index 8486b0079..58608f417 100644 --- a/examples/02_centered_form/main.cpp +++ b/examples/02_centered_form/main.cpp @@ -16,6 +16,6 @@ int main() } }; - CtcInverse ctc(f, {0,0}); + CtcInverse ctc(f, {0,0},true,true,false); DefaultFigure::pave({{0,2},{2,4},{0,10}}, ctc, 0.004); -} \ No newline at end of file +} diff --git a/python/src/core/contractors/codac2_py_CtcInverse.h b/python/src/core/contractors/codac2_py_CtcInverse.h index cdecfccae..1bbc4f55f 100644 --- a/python/src/core/contractors/codac2_py_CtcInverse.h +++ b/python/src/core/contractors/codac2_py_CtcInverse.h @@ -41,26 +41,26 @@ void export_CtcInverse(py::module& m, exported_ctc .def(py::init( - [](const AnalyticFunction& f, const D& y, bool with_centered_form) + [](const AnalyticFunction& f, const D& y, bool with_centered_form, bool with_mixed_diff) { - return std::make_unique(unaryize_function(f), y, with_centered_form); + return std::make_unique(unaryize_function(f), y, with_centered_form, with_mixed_diff); } ), - CTCINVERSE_Y_X0X_CTCINVERSE_CONST_ANALYTICFUNCTION_OUTPUTTYPE_REF_CONST_Y_REF_BOOL_BOOL, - "f"_a, "y"_a, "with_centered_form"_a = true) + CTCINVERSE_Y_X0X_CTCINVERSE_CONST_ANALYTICFUNCTION_OUTPUTTYPE_REF_CONST_Y_REF_BOOL_BOOL_BOOL, + "f"_a, "y"_a, "with_centered_form"_a = true, "with_mixed_diff"_a=false) ; if constexpr(std::is_same_v) // contractors are only associated with interval vectors { exported_ctc .def(py::init( - [](const AnalyticFunction& f, const CtcBase& c, bool with_centered_form) + [](const AnalyticFunction& f, const CtcBase& c, bool with_centered_form, bool with_mixed_diff) { - return std::make_unique(unaryize_function(f), c.copy(), with_centered_form); + return std::make_unique(unaryize_function(f), c.copy(), with_centered_form, with_mixed_diff); } ), - CTCINVERSE_Y_X0X_CTCINVERSE_CONST_ANALYTICFUNCTION_OUTPUTTYPE_REF_CONST_C_REF_BOOL_BOOL, - "f"_a, "c"_a, "with_centered_form"_a = true); + CTCINVERSE_Y_X0X_CTCINVERSE_CONST_ANALYTICFUNCTION_OUTPUTTYPE_REF_CONST_C_REF_BOOL_BOOL_BOOL, + "f"_a, "c"_a, "with_centered_form"_a = true, "with_mixed_diff"_a=false); } CONTRACT_METHODS(exported_ctc, C, @@ -286,4 +286,4 @@ void export_CtcInverse(py::module& m, exported_ctc.def("fnc", &C::fnc, CONST_ANALYTICFUNCTION_OUTPUTTYPE_REF_CTCINVERSE_Y_X0X_FNC_CONST); } -} \ No newline at end of file +} diff --git a/src/core/actions/codac2_OctaSym_operator.h b/src/core/actions/codac2_OctaSym_operator.h index 18040dd18..f0509df71 100644 --- a/src/core/actions/codac2_OctaSym_operator.h +++ b/src/core/actions/codac2_OctaSym_operator.h @@ -37,22 +37,24 @@ namespace codac2 { assert((Index)s.size() == x1.m.size()); return { - fwd(s, x1.a), + fwd(s, x1.get_a()), x1.def_domain }; } - static inline VectorType fwd_centered(const OctaSym& s, const VectorType& x1) + static inline VectorType fwd_centered([[maybe_unused]] bool mixed_diff, const OctaSym& s, const VectorType& x1) { assert((Index)s.size() == x1.m.size()); auto da = x1.da; - for(size_t i = 0 ; i < s.size() ; i++) + for(size_t i = 0 ; i < s.size() ; i++) { da.row(i) = sign(s[i])*x1.da.row(std::abs(s[i])-1); + } return { fwd(s, x1.m), - fwd(s, x1.a), + vector_map(x1.a, + [=](const IntervalVector& v1)->IntervalVector {return fwd(s,v1);}), da, x1.def_domain }; @@ -86,19 +88,19 @@ namespace codac2 return OperationExprBase>::replace_arg(old_arg_id, new_expr); } - VectorType fwd_eval(ValuesMap& v, Index total_input_size, bool natural_eval) const + VectorType fwd_eval(ValuesMap& v, Index total_input_size, bool natural_eval,bool mixed_diff) const { if(natural_eval) return AnalyticExpr::init_value( - v, OctaSymOp::fwd_natural(_s, std::get<0>(this->_x)->fwd_eval(v, total_input_size, natural_eval))); + v, OctaSymOp::fwd_natural(_s, std::get<0>(this->_x)->fwd_eval(v, total_input_size, natural_eval, mixed_diff))); else return AnalyticExpr::init_value( - v, OctaSymOp::fwd_centered(_s, std::get<0>(this->_x)->fwd_eval(v, total_input_size, natural_eval))); + v, OctaSymOp::fwd_centered(mixed_diff, _s, std::get<0>(this->_x)->fwd_eval(v, total_input_size, natural_eval, mixed_diff))); } void bwd_eval(ValuesMap& v) const { - OctaSymOp::bwd(_s, AnalyticExpr::value(v).a, std::get<0>(this->_x)->value(v).a); + OctaSymOp::bwd(_s, AnalyticExpr::value(v).get_a(), std::get<0>(this->_x)->value(v).get_a()); std::get<0>(this->_x)->bwd_eval(v); } @@ -126,4 +128,4 @@ namespace codac2 const OctaSym _s; }; -} \ No newline at end of file +} diff --git a/src/core/contractors/codac2_CtcInverse.h b/src/core/contractors/codac2_CtcInverse.h index 1d3c7c6d3..50c779997 100644 --- a/src/core/contractors/codac2_CtcInverse.h +++ b/src/core/contractors/codac2_CtcInverse.h @@ -32,15 +32,15 @@ namespace codac2 template requires IsCtcBaseOrPtr - CtcInverse(const AnalyticFunction& f, const C& ctc_y, bool with_centered_form = true, bool is_not_in = false) - : Ctc,X0,X...>(f.input_size()), _f(f), _ctc_y(ctc_y), _with_centered_form(with_centered_form), _is_not_in(is_not_in) + CtcInverse(const AnalyticFunction& f, const C& ctc_y, bool with_centered_form = true, bool with_mixed_diff = false, bool is_not_in = false) + : Ctc,X0,X...>(f.input_size()), _f(f), _ctc_y(ctc_y), _with_centered_form(with_centered_form), _with_mixed_diff(with_mixed_diff), _is_not_in(is_not_in) { assert_release([&]() { return f.output_size() == size_of(ctc_y); }() && "CtcInverse: invalid dimension of image argument ('y' or 'ctc_y')"); } - CtcInverse(const AnalyticFunction& f, const Y& y, bool with_centered_form = true, bool is_not_in = false) - : CtcInverse(f, CtcWrapper(y), with_centered_form, is_not_in) + CtcInverse(const AnalyticFunction& f, const Y& y, bool with_centered_form = true, bool with_mixed_diff = false, bool is_not_in = false) + : CtcInverse(f, CtcWrapper(y), with_centered_form, with_mixed_diff, is_not_in) { } void contract(X0& x0, X&... x) const @@ -62,7 +62,7 @@ namespace codac2 // Forward/backward algorithm: // [1/4] Forward evaluation - _f.expr()->fwd_eval(v, _f.input_size(), !_with_centered_form); + _f.expr()->fwd_eval(v, _f.input_size(), !_with_centered_form, _with_mixed_diff); auto& val_expr = _f.expr()->value(v); if(_is_not_in && !val_expr.def_domain) @@ -73,7 +73,7 @@ namespace codac2 // [2/4] Performing top contraction (from the root of the tree) with // the contractor expressing the output domain Y. // The underlying constraint is: f(x) \in [y] - ctc_y.contract(val_expr.a); + ctc_y.contract(val_expr.get_a()); // [3/4 - optional] // The contraction can be significantly improved using a centered form @@ -89,8 +89,8 @@ namespace codac2 { X0 x_mid = X0(x0.mid()); - assert(val_expr.a.size() == val_expr.m.size()); - IntervalVector fm { val_expr.a - val_expr.m }; + assert(val_expr.get_a().size() == val_expr.m.size()); + IntervalVector fm { val_expr.get_a() - val_expr.m }; if constexpr(std::is_same_v) { @@ -126,6 +126,7 @@ namespace codac2 const AnalyticFunction _f; const Collection> _ctc_y; bool _with_centered_form; + bool _with_mixed_diff; bool _is_not_in = false; }; @@ -134,63 +135,63 @@ namespace codac2 // ScalarType - CtcInverse(const AnalyticFunction&, std::initializer_list, bool = true, bool = false) -> + CtcInverse(const AnalyticFunction&, std::initializer_list, bool = true, bool = false, bool = false) -> CtcInverse; template - CtcInverse(const AnalyticFunction&, std::initializer_list, bool = true, bool = false) -> + CtcInverse(const AnalyticFunction&, std::initializer_list, bool = true, bool = false, bool = false) -> CtcInverse; - CtcInverse(const AnalyticFunction&, const Interval&, bool = true, bool = false) -> + CtcInverse(const AnalyticFunction&, const Interval&, bool = true, bool = false, bool = false) -> CtcInverse; - CtcInverse(const AnalyticFunction&, double, bool = true, bool = false) -> + CtcInverse(const AnalyticFunction&, double, bool = true, bool = false, bool = false) -> CtcInverse; template requires IsCtcBaseOrPtr - CtcInverse(const AnalyticFunction&, const C&, bool = true, bool = false) -> + CtcInverse(const AnalyticFunction&, const C&, bool = true, bool = false, bool = false) -> CtcInverse; // VectorType - CtcInverse(const AnalyticFunction&, std::initializer_list, bool = true, bool = false) -> + CtcInverse(const AnalyticFunction&, std::initializer_list, bool = true, bool = false, bool = false) -> CtcInverse; - CtcInverse(const AnalyticFunction&, std::initializer_list>, bool = true, bool = false) -> + CtcInverse(const AnalyticFunction&, std::initializer_list>, bool = true, bool = false, bool = false) -> CtcInverse; - CtcInverse(const AnalyticFunction&, const Vector&, bool = true, bool = false) -> + CtcInverse(const AnalyticFunction&, const Vector&, bool = true, bool = false, bool = false) -> CtcInverse; - CtcInverse(const AnalyticFunction&, const IntervalVector&, bool = true, bool = false) -> + CtcInverse(const AnalyticFunction&, const IntervalVector&, bool = true, bool = false, bool = false) -> CtcInverse; template requires (OtherDerived::RowsAtCompileTime == -1 && OtherDerived::ColsAtCompileTime == 1) - CtcInverse(const AnalyticFunction&, const Eigen::MatrixBase&, bool = true, bool = false) -> + CtcInverse(const AnalyticFunction&, const Eigen::MatrixBase&, bool = true, bool = false, bool = false) -> CtcInverse; template requires IsCtcBaseOrPtr - CtcInverse(const AnalyticFunction&, const C&, bool = true, bool = false) -> + CtcInverse(const AnalyticFunction&, const C&, bool = true, bool = false, bool = false) -> CtcInverse; // MatrixType - CtcInverse(const AnalyticFunction&, std::initializer_list>, bool = true, bool = false) -> + CtcInverse(const AnalyticFunction&, std::initializer_list>, bool = true, bool = false, bool = false) -> CtcInverse; - CtcInverse(const AnalyticFunction&, std::initializer_list>>, bool = true, bool = false) -> + CtcInverse(const AnalyticFunction&, std::initializer_list>>, bool = true, bool = false, bool = false) -> CtcInverse; template requires (OtherDerived::RowsAtCompileTime == -1 && OtherDerived::ColsAtCompileTime == -1) - CtcInverse(const AnalyticFunction&, const Eigen::MatrixBase&, bool = true, bool = false) -> + CtcInverse(const AnalyticFunction&, const Eigen::MatrixBase&, bool = true, bool = false, bool = false) -> CtcInverse; template requires IsCtcBaseOrPtr - CtcInverse(const AnalyticFunction&, const C&, bool = true, bool = false) -> + CtcInverse(const AnalyticFunction&, const C&, bool = true, bool = false, bool = false) -> CtcInverse; -} \ No newline at end of file +} diff --git a/src/core/contractors/codac2_CtcInverseNotIn.h b/src/core/contractors/codac2_CtcInverseNotIn.h index 1abb5eda1..05bf0cb45 100644 --- a/src/core/contractors/codac2_CtcInverseNotIn.h +++ b/src/core/contractors/codac2_CtcInverseNotIn.h @@ -23,7 +23,7 @@ namespace codac2 using Y = typename Wrapper::Domain; // Y_ is a possible Eigen expr. type using OutputType = typename ExprType::Type; - CtcInverseNotIn(const AnalyticFunction& f, const Y& y, bool with_centered_form = true) + CtcInverseNotIn(const AnalyticFunction& f, const Y& y, bool with_centered_form = true, bool with_mixed_diff = false) : CtcUnion(f.input_size()) { const bool is_not_in = true; @@ -32,17 +32,17 @@ namespace codac2 { Y w = complem_y & f_codomain; if(!w.is_empty() && !w.is_subset(y)) - *this |= CtcInverse(f, complem_y, with_centered_form, is_not_in); + *this |= CtcInverse(f, complem_y, with_centered_form, with_mixed_diff, is_not_in); } } template requires IsCtcBaseOrPtr - CtcInverseNotIn(const AnalyticFunction& f, const C& ctc_compl, bool with_centered_form = true) + CtcInverseNotIn(const AnalyticFunction& f, const C& ctc_compl, bool with_centered_form = true, bool with_mixed_diff = false) : CtcUnion(f.input_size()) { const bool is_not_in = true; - *this |= CtcInverse(f, ctc_compl, with_centered_form, is_not_in); + *this |= CtcInverse(f, ctc_compl, with_centered_form, with_mixed_diff, is_not_in); } void contract(X0& x0, X&... x) const @@ -52,7 +52,7 @@ namespace codac2 }; // Template deduction guides - // Same as CtcInverse + // Same as CtcInverse ( but without one boolean? how did this work? [DM] ) // ScalarType @@ -115,4 +115,4 @@ namespace codac2 requires IsCtcBaseOrPtr CtcInverseNotIn(const AnalyticFunction&, const C&, bool = true, bool = false) -> CtcInverseNotIn; -} \ No newline at end of file +} diff --git a/src/core/domains/tube/codac2_Tube_operator.h b/src/core/domains/tube/codac2_Tube_operator.h index 4666475cf..9ecf07613 100644 --- a/src/core/domains/tube/codac2_Tube_operator.h +++ b/src/core/domains/tube/codac2_Tube_operator.h @@ -31,7 +31,7 @@ namespace codac2 return { fwd(x1,x2.m), - fwd(x1,x2.a), + fwd(x1,x2.get_a()), d, x2.def_domain && x1.tdomain()->t0_tf().is_superset(x2.m) }; @@ -67,15 +67,15 @@ namespace codac2 return OperationExprBase>::replace_arg(old_arg_id, new_expr); } - T fwd_eval(ValuesMap& v, Index total_input_size, bool natural_eval) const + T fwd_eval(ValuesMap& v, Index total_input_size, bool natural_eval, bool mixed_diff) const { return AnalyticExpr::init_value( - v, TubeOp::fwd(_x1, std::get<0>(this->_x)->fwd_eval(v, total_input_size, natural_eval))); + v, TubeOp::fwd(_x1, std::get<0>(this->_x)->fwd_eval(v, total_input_size, natural_eval,mixed_diff))); } void bwd_eval(ValuesMap& v) const { - TubeOp::bwd(_x1, AnalyticExpr::value(v).a, std::get<0>(this->_x)->value(v).a); + TubeOp::bwd(_x1, AnalyticExpr::value(v).get_a(), std::get<0>(this->_x)->value(v).get_a()); std::get<0>(this->_x)->bwd_eval(v); } diff --git a/src/core/functions/analytic/codac2_AnalyticExpr.h b/src/core/functions/analytic/codac2_AnalyticExpr.h index 3c8ebeb8b..bcbe4a70c 100644 --- a/src/core/functions/analytic/codac2_AnalyticExpr.h +++ b/src/core/functions/analytic/codac2_AnalyticExpr.h @@ -26,7 +26,7 @@ namespace codac2 { public: - virtual T fwd_eval(ValuesMap& v, Index total_input_size, bool natural_eval) const = 0; + virtual T fwd_eval(ValuesMap& v, Index total_input_size, bool natural_eval, bool mixed_diff) const = 0; virtual void bwd_eval(ValuesMap& v) const = 0; virtual std::pair output_shape() const = 0; @@ -86,18 +86,21 @@ namespace codac2 return OperationExprBase...>::replace_arg(old_arg_id, new_expr); } - Y fwd_eval(ValuesMap& v, Index total_input_size, bool natural_eval) const override + Y fwd_eval(ValuesMap& v, Index total_input_size, bool natural_eval, bool mixed_diff) const override { return std::apply( - [this,&v,total_input_size,natural_eval](auto &&... x) + [this,&v,total_input_size,natural_eval,mixed_diff](auto &&... x) { if(natural_eval) return AnalyticExpr::init_value(v, - C::fwd_natural(x->fwd_eval(v, total_input_size, natural_eval)...)); + C::fwd_natural(x->fwd_eval(v, total_input_size, natural_eval,mixed_diff)...)); - else + else if (!mixed_diff) return AnalyticExpr::init_value(v, - C::fwd_centered(x->fwd_eval(v, total_input_size, natural_eval)...)); + C::fwd_centered(false,x->fwd_eval(v, total_input_size, natural_eval,mixed_diff)...)); + else + return AnalyticExpr::init_value(v, + C::fwd_centered(true,x->fwd_eval(v, total_input_size, natural_eval,mixed_diff)...)); }, this->_x); } @@ -108,7 +111,7 @@ namespace codac2 std::apply([&v,y](auto &&... x) { - C::bwd(y.a, x->value(v).a...); + C::bwd(y.a.back(), x->value(v).a.back()...); }, this->_x); std::apply([&v](auto &&... x) diff --git a/src/core/functions/analytic/codac2_AnalyticFunction.h b/src/core/functions/analytic/codac2_AnalyticFunction.h index 334be86e7..0e651975b 100644 --- a/src/core/functions/analytic/codac2_AnalyticFunction.h +++ b/src/core/functions/analytic/codac2_AnalyticFunction.h @@ -29,7 +29,9 @@ namespace codac2 { NATURAL = 0x01, CENTERED = 0x02, - DEFAULT = 0x03 // corresponds to (NATURAL|CENTERED) + DEFAULT = 0x03, // corresponds to (NATURAL|CENTERED) + MIXED = 0x04, // note: not compatible with CENTERED + MIXED_FULL = 0x05 // corresponds to (NATURAL|MIXED) }; inline EvalMode operator&(EvalMode a, EvalMode b) @@ -38,6 +40,36 @@ namespace codac2 inline EvalMode operator|(EvalMode a, EvalMode b) { return static_cast(static_cast(a) | static_cast(b)); } + /* partial mid to build the a field */ + inline Interval partial_mid(const Interval &a,[[maybe_unused]] Index i) { + return a; + } + inline IntervalVector partial_mid(const IntervalVector &a,Index i) { + IntervalVector m(a); + Index sz = m.size(); + for (Index j=i+1;j { // Mainly used to take advantage of initializer lists in AnalyticFunction constructors. @@ -102,12 +134,12 @@ namespace codac2 { case EvalMode::NATURAL: { - return eval_(x...).a; + return eval_(x...).get_a(); } case EvalMode::CENTERED: { - auto x_ = eval_(x...); + auto x_ = eval_(x...); auto flatten_x = IntervalVector(cart_prod(x...)); assert(x_.da.rows() == x_.a.size() && x_.da.cols() == flatten_x.size()); @@ -125,34 +157,65 @@ namespace codac2 } } + case EvalMode::MIXED: + case EvalMode::MIXED_FULL: + { + auto x_ = eval_(x...); + // If the centered form is not available for this expression... + if(m==EvalMode::MIXED_FULL && (x_.da.size() == 0 // .. because some parts have not yet been implemented, + || !x_.def_domain)) // .. or due to restrictions in the derivative definition domain + return x_.get_a(); // natural evaluation + + else { + auto flatten_x = IntervalVector(cart_prod(x...)); + assert(x_.da.rows() == x_.a.size() && x_.da.cols() == flatten_x.size()); + + if constexpr(std::is_same_v) { + Interval ret = x_.m + (x_.da*(flatten_x-flatten_x.mid()))[0]; + return (m==EvalMode::MIXED_FULL ? x_.get_a() & ret : ret); + } + else if constexpr(std::is_same_v) { + IntervalVector ret = x_.m + (x_.da*(flatten_x-flatten_x.mid())).col(0); + return (m==EvalMode::MIXED_FULL ? x_.get_a() & ret : ret); + } + else + { + static_assert(std::is_same_v); + IntervalMatrix ret = x_.m + (x_.da*(flatten_x-flatten_x.mid())) + .reshaped(x_.m.rows(), x_.m.cols()); + return (m==EvalMode::MIXED_FULL ? x_.get_a() & ret : ret); + } + } + } + case EvalMode::DEFAULT: default: { - auto x_ = eval_(x...); + auto x_ = eval_(x...); // If the centered form is not available for this expression... if(x_.da.size() == 0 // .. because some parts have not yet been implemented, || !x_.def_domain) // .. or due to restrictions in the derivative definition domain - return x_.a; // natural evaluation + return x_.get_a(); // natural evaluation else { auto flatten_x = IntervalVector(cart_prod(x...)); if constexpr(std::is_same_v) - return x_.a & (x_.m + (x_.da*(flatten_x-flatten_x.mid()))[0]); + return x_.get_a() & (x_.m + (x_.da*(flatten_x-flatten_x.mid()))[0]); else if constexpr(std::is_same_v) { assert(x_.da.rows() == x_.a.size() && x_.da.cols() == flatten_x.size()); - return x_.a & (x_.m + (x_.da*(flatten_x-flatten_x.mid())).col(0)); + return x_.get_a() & (x_.m + (x_.da*(flatten_x-flatten_x.mid())).col(0)); } else { static_assert(std::is_same_v); assert(x_.da.rows() == x_.a.size() && x_.da.cols() == flatten_x.size()); - return x_.a & (x_.m +(x_.da*(flatten_x-flatten_x.mid())) + return x_.get_a() & (x_.m +(x_.da*(flatten_x-flatten_x.mid())) .reshaped(x_.m.rows(),x_.m.cols())); } } @@ -172,11 +235,11 @@ namespace codac2 ); } - template + template auto diff(const Args&... x) const { check_valid_inputs(x...); - return eval_(x...).da; + return eval_(x...).da; } template @@ -209,7 +272,7 @@ namespace codac2 return y; } - template + template requires std::is_same_v && ((!std::is_same_v::Type>) && ...) Parallelepiped parallelepiped_eval(const Args&... x) const { @@ -222,10 +285,11 @@ namespace codac2 IntervalVector Y = this->eval(((typename Wrapper::Domain)(x)).mid()...); Vector c = Y.mid(); - Matrix A = this->diff(((typename Wrapper::Domain)(x)).mid()...).mid(); + Matrix A = this->diff(((typename Wrapper::Domain)(x)).mid()...).mid(); // Maximum error computation - double rho = error_peibos(Y, c, this->diff(x...), A, cart_prod(x...)); + double rho = error_peibos(Y, c, this->diff(x...), A, cart_prod(x...)); + // TODO : check the correction of this // Inflation of the parallelepiped Matrix A_inf = inflate_flat_parallelepiped(A, (cart_prod(x...).template cast()).rad(), rho); @@ -300,30 +364,47 @@ namespace codac2 for(Index k = p ; k < p+size_of(x) ; k++) d(k-p,k) = 1.; + + std::vector vecx; + vecx.reserve(this->args().total_size()); + + typename D_TYPE::Domain midx = typename D_TYPE::Domain(x).mid(); + Index l; + for(l=0;l(partial_mid))(x,l-p)); + } + for(;largs().total_size();l++) { + vecx.push_back(x); + } v[this->args()[i]->unique_id()] = - std::make_shared(typename D_TYPE::Domain(x).mid(), x, d, true); + std::make_shared(midx, + vecx, d, true); } template void intersect_value_from_arg_map(const ValuesMap& v, D& x, Index i) const { assert(v.find(this->args()[i]->unique_id()) != v.end() && "argument cannot be found"); - x &= std::dynamic_pointer_cast::Type>(v.at(this->args()[i]->unique_id()))->a; + x &= std::dynamic_pointer_cast::Type>(v.at(this->args()[i]->unique_id()))->get_a(); } - template + template auto eval_(const Args&... x) const { ValuesMap v; if constexpr(sizeof...(Args) == 0) - return this->expr()->fwd_eval(v, 0, NATURAL_EVAL); + return this->expr()->fwd_eval(v, 0, NATURAL_EVAL, MIXED_DIFF); else { fill_from_args(v, x...); - return this->expr()->fwd_eval(v, this->input_size(), NATURAL_EVAL); + return this->expr()->fwd_eval(v, this->input_size(), NATURAL_EVAL, MIXED_DIFF); } } @@ -372,4 +453,4 @@ namespace codac2 AnalyticFunction(const FunctionArgsList&, const T&) -> AnalyticFunction::Type>; -} \ No newline at end of file +} diff --git a/src/core/functions/analytic/codac2_AnalyticFunction_impl.h b/src/core/functions/analytic/codac2_AnalyticFunction_impl.h index 3ec22ddd3..f84b9a9d6 100644 --- a/src/core/functions/analytic/codac2_AnalyticFunction_impl.h +++ b/src/core/functions/analytic/codac2_AnalyticFunction_impl.h @@ -45,12 +45,24 @@ namespace codac2 { case EvalMode::NATURAL: { - return eval_(x...).a; + return eval_(x...).get_a(); } case EvalMode::CENTERED: { - auto x_ = eval_(x...); + auto x_ = eval_(x...); + auto flatten_x = cart_prod(x...); + assert(x_.da.rows() == x_.a.size() && x_.da.cols() == flatten_x.size()); + + if constexpr(std::is_same_v) + return x_.m + (x_.da*(flatten_x-flatten_x.mid()))[0]; + else + return x_.m + (x_.da*(flatten_x-flatten_x.mid())).col(0); + } + + case EvalMode::MIXED: + { + auto x_ = eval_(x...); auto flatten_x = cart_prod(x...); assert(x_.da.rows() == x_.a.size() && x_.da.cols() == flatten_x.size()); @@ -60,6 +72,26 @@ namespace codac2 return x_.m + (x_.da*(flatten_x-flatten_x.mid())).col(0); } + case EvalMode::FULL_MIXED: + { + auto x_ = eval_(x...); + + if(x_.da.size() == 0) // if the centered form is not available for this expression + return eval(EvalMode::NATURAL, x...); + + else + { + auto flatten_x = cart_prod(x...); + if constexpr(std::is_same_v) + return x_.get_a() & (x_.m + (x_.da*(flatten_x-flatten_x.mid()))[0]); + else + { + assert(x_.da.rows() == x_.a.size() && x_.da.cols() == flatten_x.size()); + return x_.get_a() & (x_.m + (x_.da*(flatten_x-flatten_x.mid())).col(0)); + } + } + } + case EvalMode::DEFAULT: default: { @@ -93,11 +125,11 @@ namespace codac2 template requires std::is_base_of_v - template + template inline auto AnalyticFunction::diff(const Args&... x) const { check_valid_inputs(x...); - return eval_(x...).da; + return eval_(x...).da; } template @@ -178,7 +210,7 @@ namespace codac2 inline void AnalyticFunction::intersect_value_from_arg_map(const ValuesMap& v, D& x, Index i) const { assert(v.find(this->args()[i]->unique_id()) != v.end() && "argument cannot be found"); - x &= std::dynamic_pointer_cast::Type>(v.at(this->args()[i]->unique_id()))->a; + x &= std::dynamic_pointer_cast::Type>(v.at(this->args()[i]->unique_id()))->get_a(); } template @@ -192,18 +224,18 @@ namespace codac2 template requires std::is_base_of_v - template + template inline auto AnalyticFunction::eval_(const Args&... x) const { ValuesMap v; if constexpr(sizeof...(Args) == 0) - return this->expr()->fwd_eval(v, 0, NATURAL_EVAL); + return this->expr()->fwd_eval(v, 0, NATURAL_EVAL, false); else { fill_from_args(v, x...); - return this->expr()->fwd_eval(v, cart_prod(x...).size(), NATURAL_EVAL); // todo: improve size computation + return this->expr()->fwd_eval(v, cart_prod(x...).size(), NATURAL_EVAL, MIXED_DIFF); // todo: improve size computation } } @@ -218,4 +250,4 @@ namespace codac2 assert_release(this->_args.total_size() == n && "Invalid arguments: wrong number of input arguments"); } -} \ No newline at end of file +} diff --git a/src/core/functions/analytic/codac2_AnalyticType.h b/src/core/functions/analytic/codac2_AnalyticType.h index 2df5f0338..df0b7790c 100644 --- a/src/core/functions/analytic/codac2_AnalyticType.h +++ b/src/core/functions/analytic/codac2_AnalyticType.h @@ -9,6 +9,8 @@ #pragma once +#include +#include #include "codac2_Interval.h" #include "codac2_Vector.h" #include "codac2_Matrix.h" @@ -30,27 +32,55 @@ namespace codac2 using Domain = D; D m; - D a; + std::vector a; IntervalMatrix da; bool def_domain; AnalyticType() = delete; AnalyticType(const D& a_, bool def_domain_) - : a(a_), def_domain(def_domain_) + : a(1,a_), def_domain(def_domain_) { } AnalyticType(const D& m_, const D& a_, const IntervalMatrix& da_, bool def_domain_) + : m(m_), a(1,a_), da(da_), def_domain(def_domain_) + { } + + AnalyticType(const D& m_, const std::vector& a_, const IntervalMatrix& da_, bool def_domain_) : m(m_), a(a_), da(da_), def_domain(def_domain_) { } AnalyticType& operator&=(const AnalyticType& x) { - a &= x.a; + for (Index i=0;iget_a(true,i) &= x.get_a(true,i); + } // restore this? da &= x.da; def_domain &= x.def_domain; return *this; } + + inline const D& get_a() const { + return a.back(); + } + + inline const D& get_a(bool mixed_diff, size_t i) const { + if (!mixed_diff || a.size()==1) return a.back(); + assert_release(i; @@ -62,4 +92,42 @@ namespace codac2 { return ((a.da.size() == 0) || ...); } -} \ No newline at end of file + + /* functional map of vector for operator evaluation of .a field */ + template + inline std::vector vector_map(const std::vector &t, + std::functionf) { + std::vector ret(t.size()); + auto r1=t.cbegin(); + auto rt=ret.begin(); + while (rt!=ret.end()) { + (*rt) = f(*r1); + rt++; r1++; + } + return ret; + } + template + inline std::vector vector_map(const std::vector &t1, + const std::vector &t2, + std::functionf) { + std::vector ret(std::max(t1.size(),t2.size())); + for (size_t i=0;i + inline std::vector vector_map(const std::vector &t1, + const std::vector &t2, + const std::vector &t3, + std::functionf) { + std::vector ret(std::max({ t1.size(),t2.size(),t3.size() })); + for (size_t i=0;i>(*this); } - T fwd_eval(ValuesMap& v, Index total_input_size, bool natural_eval) const override + T fwd_eval(ValuesMap& v, Index total_input_size, bool natural_eval,[[maybe_unused]] bool mixed_diff) const override { if(natural_eval) return AnalyticExpr::init_value(v, T( @@ -43,7 +43,7 @@ namespace codac2 true )); - else + else return AnalyticExpr::init_value(v, T( // the mid is not considered for const values in centered form expression: _x, @@ -57,7 +57,7 @@ namespace codac2 void bwd_eval(ValuesMap& v) const override { - AnalyticExpr::value(v).a &= _x; + AnalyticExpr::value(v).get_a() &= _x; } std::pair output_shape() const override @@ -102,4 +102,4 @@ namespace codac2 } } -#include "codac2_analytic_constants_impl.h" \ No newline at end of file +#include "codac2_analytic_constants_impl.h" diff --git a/src/core/functions/analytic/codac2_analytic_variables.h b/src/core/functions/analytic/codac2_analytic_variables.h index 3d942ccb2..83b66e4a5 100644 --- a/src/core/functions/analytic/codac2_analytic_variables.h +++ b/src/core/functions/analytic/codac2_analytic_variables.h @@ -33,7 +33,7 @@ namespace codac2 return AnalyticExpr::unique_id(); } - T fwd_eval(ValuesMap& v, [[maybe_unused]] Index total_input_size, [[maybe_unused]] bool natural_eval) const + T fwd_eval(ValuesMap& v, [[maybe_unused]] Index total_input_size, [[maybe_unused]] bool natural_eval, [[maybe_unused]] bool mixed_diff) const { return AnalyticExpr::value(v); } diff --git a/src/core/operators/codac2_abs.h b/src/core/operators/codac2_abs.h index 452342b67..9c3eacda4 100644 --- a/src/core/operators/codac2_abs.h +++ b/src/core/operators/codac2_abs.h @@ -31,7 +31,7 @@ namespace codac2 static Interval fwd(const Interval& x1); static ScalarType fwd_natural(const ScalarType& x1); - static ScalarType fwd_centered(const ScalarType& x1); + static ScalarType fwd_centered(bool mixed_diff,const ScalarType& x1); static void bwd(const Interval& y, Interval& x1); }; @@ -54,25 +54,28 @@ namespace codac2 inline ScalarType AbsOp::fwd_natural(const ScalarType& x1) { return { - fwd(x1.a), - x1.def_domain && !x1.a.contains(0.) // def domain of the derivative of abs + fwd(x1.get_a()), + x1.def_domain && !x1.get_a().contains(0.) + // def domain of the derivative of abs }; } - inline ScalarType AbsOp::fwd_centered(const ScalarType& x1) + inline ScalarType AbsOp::fwd_centered(bool mixed_diff,const ScalarType& x1) { if(centered_form_not_available_for_args(x1)) return fwd_natural(x1); IntervalMatrix d(1,x1.da.size()); - for(Index i = 0 ; i < d.size() ; i++) - d(0,i) = (x1.a/abs(x1.a))*x1.da(0,i); + for(Index i = 0 ; i < d.size() ; i++) { + const Interval& x1a = x1.get_a(mixed_diff,i); + d(0,i) = (x1a/abs(x1a))*x1.da(0,i); + } return { fwd(x1.m), - fwd(x1.a), + vector_map(x1.a,&fwd), d, - x1.def_domain && !x1.a.contains(0.) // def domain of the derivative of abs + x1.def_domain && !x1.get_a().contains(0.) // def domain of the derivative of abs }; } diff --git a/src/core/operators/codac2_acos.h b/src/core/operators/codac2_acos.h index e4f57dc05..a57f0b040 100644 --- a/src/core/operators/codac2_acos.h +++ b/src/core/operators/codac2_acos.h @@ -31,7 +31,7 @@ namespace codac2 static Interval fwd(const Interval& x1); static ScalarType fwd_natural(const ScalarType& x1); - static ScalarType fwd_centered(const ScalarType& x1); + static ScalarType fwd_centered(bool mixed_diff,const ScalarType& x1); static void bwd(const Interval& y, Interval& x1); }; @@ -54,28 +54,30 @@ namespace codac2 inline ScalarType AcosOp::fwd_natural(const ScalarType& x1) { return { - fwd(x1.a), - x1.a.is_subset({-1,1}) // def domain of acos - && !x1.a.contains(1.) // def domain of the derivative of acos + fwd(x1.get_a()), + x1.get_a().is_subset({-1,1}) // def domain of acos + && !x1.get_a().contains(1.) // def domain of the derivative of acos && x1.def_domain }; } - inline ScalarType AcosOp::fwd_centered(const ScalarType& x1) + inline ScalarType AcosOp::fwd_centered(bool mixed_diff, const ScalarType& x1) { if(centered_form_not_available_for_args(x1)) return fwd_natural(x1); IntervalMatrix d(1,x1.da.size()); - for(Index i = 0 ; i < d.size() ; i++) - d(0,i) = -x1.da(0,i)/sqrt(1.-sqr(x1.a)); + for(Index i = 0 ; i < d.size() ; i++) { + const Interval& x1a = x1.get_a(mixed_diff,i); + d(0,i) = -x1.da(0,i)/sqrt(1.-sqr(x1a)); + } return { fwd(x1.m), - fwd(x1.a), + vector_map(x1.a,&fwd), d, - x1.a.is_subset({-1,1}) // def domain of acos - && !x1.a.contains(1.) // def domain of the derivative of acos + x1.get_a().is_subset({-1,1}) // def domain of acos + && !x1.get_a().contains(1.) // def domain of the derivative of acos && x1.def_domain }; } diff --git a/src/core/operators/codac2_arith_add.h b/src/core/operators/codac2_arith_add.h index 41acede06..d10ce769a 100644 --- a/src/core/operators/codac2_arith_add.h +++ b/src/core/operators/codac2_arith_add.h @@ -34,17 +34,17 @@ namespace codac2 static Interval fwd(const Interval& x1); static ScalarType fwd_natural(const ScalarType& x1); - static ScalarType fwd_centered(const ScalarType& x1); + static ScalarType fwd_centered(bool mixed_diff, const ScalarType& x1); static void bwd(const Interval& y, Interval& x1); static IntervalVector fwd(const IntervalVector& x1); static VectorType fwd_natural(const VectorType& x1); - static VectorType fwd_centered(const VectorType& x1); + static VectorType fwd_centered(bool mixed_diff, const VectorType& x1); static void bwd(const IntervalVector& y, IntervalVector& x1); static IntervalMatrix fwd(const IntervalMatrix& x1); static MatrixType fwd_natural(const MatrixType& x1); - static MatrixType fwd_centered(const MatrixType& x1); + static MatrixType fwd_centered(bool mixed_diff, const MatrixType& x1); static void bwd(const IntervalMatrix& y, IntervalMatrix& x1); static std::pair output_shape(const std::pair &s1); @@ -67,17 +67,17 @@ namespace codac2 static Interval fwd(const Interval& x1, const Interval& x2); static ScalarType fwd_natural(const ScalarType& x1, const ScalarType& x2); - static ScalarType fwd_centered(const ScalarType& x1, const ScalarType& x2); + static ScalarType fwd_centered(bool mixed_diff, const ScalarType& x1, const ScalarType& x2); static void bwd(const Interval& y, Interval& x1, Interval& x2); static IntervalVector fwd(const IntervalVector& x1, const IntervalVector& x2); static VectorType fwd_natural(const VectorType& x1, const VectorType& x2); - static VectorType fwd_centered(const VectorType& x1, const VectorType& x2); + static VectorType fwd_centered(bool mixed_diff, const VectorType& x1, const VectorType& x2); static void bwd(const IntervalVector& y, IntervalVector& x1, IntervalVector& x2); static IntervalMatrix fwd(const IntervalMatrix& x1, const IntervalMatrix& x2); static MatrixType fwd_natural(const MatrixType& x1, const MatrixType& x2); - static MatrixType fwd_centered(const MatrixType& x1, const MatrixType& x2); + static MatrixType fwd_centered(bool mixed_diff, const MatrixType& x1, const MatrixType& x2); static void bwd(const IntervalMatrix& y, IntervalMatrix& x1, IntervalMatrix& x2); }; @@ -134,16 +134,17 @@ namespace codac2 inline ScalarType AddOp::fwd_natural(const ScalarType& x1) { return { - fwd(x1.a), + fwd(x1.get_a()), x1.def_domain }; } - inline ScalarType AddOp::fwd_centered(const ScalarType& x1) + inline ScalarType AddOp::fwd_centered([[maybe_unused]] bool mixed_diff, const ScalarType& x1) { return { fwd(x1.m), - fwd(x1.a), + vector_map(x1.a, + static_cast(&fwd)), x1.da, x1.def_domain }; @@ -165,16 +166,17 @@ namespace codac2 inline VectorType AddOp::fwd_natural(const VectorType& x1) { return { - fwd(x1.a), + fwd(x1.get_a()), x1.def_domain }; } - inline VectorType AddOp::fwd_centered(const VectorType& x1) + inline VectorType AddOp::fwd_centered([[maybe_unused]] bool mixed_diff, const VectorType& x1) { return { fwd(x1.m), - fwd(x1.a), + vector_map(x1.a, + static_cast(&fwd)), x1.da, x1.def_domain }; @@ -193,19 +195,20 @@ namespace codac2 inline MatrixType AddOp::fwd_natural(const MatrixType& x1) { return { - fwd(x1.a), + fwd(x1.get_a()), x1.def_domain }; } - inline MatrixType AddOp::fwd_centered(const MatrixType& x1) + inline MatrixType AddOp::fwd_centered([[maybe_unused]] bool mixed_diff, const MatrixType& x1) { if(centered_form_not_available_for_args(x1)) return fwd_natural(x1); return { fwd(x1.m), - fwd(x1.a), + vector_map(x1.a, + static_cast(&fwd)), x1.da, // unary plus x1.def_domain }; @@ -224,12 +227,12 @@ namespace codac2 inline ScalarType AddOp::fwd_natural(const ScalarType& x1, const ScalarType& x2) { return { - fwd(x1.a, x2.a), + fwd(x1.get_a(), x2.get_a()), x1.def_domain && x2.def_domain }; } - inline ScalarType AddOp::fwd_centered(const ScalarType& x1, const ScalarType& x2) + inline ScalarType AddOp::fwd_centered([[maybe_unused]] bool mixed_diff,const ScalarType& x1, const ScalarType& x2) { if(centered_form_not_available_for_args(x1,x2)) return fwd_natural(x1,x2); @@ -237,7 +240,8 @@ namespace codac2 assert(x1.da.size() == x2.da.size()); return { fwd(x1.m, x2.m), - fwd(x1.a, x2.a), + vector_map(x1.a,x2.a, + static_cast(&fwd)), x1.da + x2.da, x1.def_domain && x2.def_domain }; @@ -261,12 +265,12 @@ namespace codac2 inline VectorType AddOp::fwd_natural(const VectorType& x1, const VectorType& x2) { return { - fwd(x1.a, x2.a), + fwd(x1.get_a(), x2.get_a()), x1.def_domain && x2.def_domain }; } - inline VectorType AddOp::fwd_centered(const VectorType& x1, const VectorType& x2) + inline VectorType AddOp::fwd_centered([[maybe_unused]] bool mixed_diff,const VectorType& x1, const VectorType& x2) { if(centered_form_not_available_for_args(x1,x2)) return fwd_natural(x1,x2); @@ -274,7 +278,8 @@ namespace codac2 assert(x1.da.rows() == x2.da.rows() && x1.da.cols() == x2.da.cols()); return { fwd(x1.m, x2.m), - fwd(x1.a, x2.a), + vector_map(x1.a,x2.a, + static_cast(&fwd)), x1.da + x2.da, x1.def_domain && x2.def_domain }; @@ -296,19 +301,20 @@ namespace codac2 inline MatrixType AddOp::fwd_natural(const MatrixType& x1, const MatrixType& x2) { return { - fwd(x1.a, x2.a), + fwd(x1.get_a(), x2.get_a()), x1.def_domain && x2.def_domain }; } - inline MatrixType AddOp::fwd_centered(const MatrixType& x1, const MatrixType& x2) + inline MatrixType AddOp::fwd_centered([[maybe_unused]] bool mixed_diff,const MatrixType& x1, const MatrixType& x2) { if(centered_form_not_available_for_args(x1,x2)) return fwd_natural(x1,x2); return { fwd(x1.m, x2.m), - fwd(x1.a, x2.a), + vector_map(x1.a,x2.a, + static_cast(&fwd)), x1.da + x2.da, x1.def_domain && x2.def_domain }; diff --git a/src/core/operators/codac2_arith_div.h b/src/core/operators/codac2_arith_div.h index 5fa23b9ae..23603e573 100644 --- a/src/core/operators/codac2_arith_div.h +++ b/src/core/operators/codac2_arith_div.h @@ -36,17 +36,17 @@ namespace codac2 static Interval fwd(const Interval& x1, const Interval& x2); static ScalarType fwd_natural(const ScalarType& x1, const ScalarType& x2); - static ScalarType fwd_centered(const ScalarType& x1, const ScalarType& x2); + static ScalarType fwd_centered(bool mixed_diff, const ScalarType& x1, const ScalarType& x2); static void bwd(const Interval& y, Interval& x1, Interval& x2); static IntervalVector fwd(const IntervalVector& x1, const Interval& x2); static VectorType fwd_natural(const VectorType& x1, const ScalarType& x2); - static VectorType fwd_centered(const VectorType& x1, const ScalarType& x2); + static VectorType fwd_centered(bool mixed_diff, const VectorType& x1, const ScalarType& x2); static void bwd(const IntervalVector& y, IntervalVector& x1, Interval& x2); static IntervalMatrix fwd(const IntervalMatrix& x1, const Interval& x2); static MatrixType fwd_natural(const MatrixType& x1, const ScalarType& x2); - static MatrixType fwd_centered(const MatrixType& x1, const ScalarType& x2); + static MatrixType fwd_centered(bool mixed_diff, const MatrixType& x1, const ScalarType& x2); static void bwd(const IntervalMatrix& y, IntervalMatrix& x1, Interval& x2); }; @@ -81,12 +81,12 @@ namespace codac2 inline ScalarType DivOp::fwd_natural(const ScalarType& x1, const ScalarType& x2) { return { - fwd(x1.a, x2.a), - x1.def_domain && x2.def_domain && !x2.a.contains(0.) // def domain of the derivative of div + fwd(x1.get_a(), x2.get_a()), + x1.def_domain && x2.def_domain && !x2.get_a().contains(0.) // def domain of the derivative of div }; } - inline ScalarType DivOp::fwd_centered(const ScalarType& x1, const ScalarType& x2) + inline ScalarType DivOp::fwd_centered(bool mixed_diff, const ScalarType& x1, const ScalarType& x2) { if(centered_form_not_available_for_args(x1,x2)) return fwd_natural(x1,x2); @@ -94,14 +94,18 @@ namespace codac2 assert(x1.da.size() == x2.da.size()); IntervalMatrix d(1,x1.da.size()); - for(Index i = 0 ; i < d.size() ; i++) - d(0,i) = (x1.da(0,i)*x2.a-x1.a*x2.da(0,i))/sqr(x2.a); + for(Index i = 0 ; i < d.size() ; i++) { + const Interval &x1a = x1.get_a(mixed_diff,i); + const Interval &x2a = x2.get_a(mixed_diff,i); + d(0,i) = (x1.da(0,i)*x2a-x1a*x2.da(0,i))/sqr(x2a); + } return { fwd(x1.m, x2.m), - fwd(x1.a, x2.a), + vector_map(x1.a,x2.a, + static_cast(&DivOp::fwd)), d, - x1.def_domain && x2.def_domain && !x2.a.contains(0.) // def domain of the derivative of div + x1.def_domain && x2.def_domain && !x2.get_a().contains(0.) // def domain of the derivative of div }; } @@ -127,12 +131,12 @@ namespace codac2 inline VectorType DivOp::fwd_natural(const VectorType& x1, const ScalarType& x2) { return { - fwd(x1.a, x2.a), - x1.def_domain && x2.def_domain && !x2.a.contains(0.) // def domain of the derivative of div + fwd(x1.get_a(), x2.get_a()), + x1.def_domain && x2.def_domain && !x2.get_a().contains(0.) // def domain of the derivative of div }; } - inline VectorType DivOp::fwd_centered(const VectorType& x1, const ScalarType& x2) + inline VectorType DivOp::fwd_centered(bool mixed_diff,const VectorType& x1, const ScalarType& x2) { if(centered_form_not_available_for_args(x1,x2)) return fwd_natural(x1,x2); @@ -141,15 +145,20 @@ namespace codac2 assert(x1.da.cols() == x2.da.cols()); IntervalMatrix d(x1.da.rows(),x1.da.cols()); - for(Index j = 0 ; j < d.cols() ; j++) - for(Index i = 0 ; i < d.rows() ; i++) - d(i,j) = x1.da(i,j)/x2.a - x1.a[i]*x2.da(0,j)/sqr(x2.a); + for(Index j = 0 ; j < d.cols() ; j++) { + const IntervalVector& x1a = x1.get_a(mixed_diff,j); + const Interval& x2a = x2.get_a(mixed_diff,j); + for(Index i = 0 ; i < d.rows() ; i++) + d(i,j) = x1.da(i,j)/x2a - x1a[i]*x2.da(0,j)/sqr(x2a); + } + return { fwd(x1.m, x2.m), - fwd(x1.a, x2.a), + vector_map(x1.a,x2.a, + static_cast(&DivOp::fwd)), d, - x1.def_domain && x2.def_domain && !x2.a.contains(0.) // def domain of the derivative of div + x1.def_domain && x2.def_domain && !x2.get_a().contains(0.) // def domain of the derivative of div }; } @@ -168,12 +177,12 @@ namespace codac2 inline MatrixType DivOp::fwd_natural(const MatrixType& x1, const ScalarType& x2) { return { - fwd(x1.a, x2.a), - x1.def_domain && x2.def_domain && !x2.a.contains(0.) // def domain of the derivative of div + fwd(x1.get_a(), x2.get_a()), + x1.def_domain && x2.def_domain && !x2.get_a().contains(0.) // def domain of the derivative of div }; } - inline MatrixType DivOp::fwd_centered(const MatrixType& x1, const ScalarType& x2) + inline MatrixType DivOp::fwd_centered(bool mixed_diff, const MatrixType& x1, const ScalarType& x2) { if(centered_form_not_available_for_args(x1,x2)) return fwd_natural(x1,x2); @@ -182,16 +191,20 @@ namespace codac2 assert(x1.da.cols() == x2.da.cols()); IntervalMatrix d(x1.da.rows(),x1.da.cols()); - for(Index j = 0 ; j < d.cols() ; j++) + for(Index j = 0 ; j < d.cols() ; j++) { + const IntervalMatrix& x1a = x1.get_a(mixed_diff,j); + const Interval& x2a = x2.get_a(mixed_diff,j); for(Index i = 0 ; i < d.rows() ; i++) - d(i,j) = x1.da(i,j)/x2.a - - x1.a.reshaped()[i]*x2.da(0,j)/sqr(x2.a); + d(i,j) = x1.da(i,j)/x2a - + x1a.reshaped()[i]*x2.da(0,j)/sqr(x2a); + } return { fwd(x1.m, x2.m), - fwd(x1.a, x2.a), + vector_map(x1.a,x2.a, + static_cast(&DivOp::fwd)), d, - x1.def_domain && x2.def_domain && !x2.a.contains(0.) // def domain of the derivative of div + x1.def_domain && x2.def_domain && !x2.get_a().contains(0.) // def domain of the derivative of div }; } diff --git a/src/core/operators/codac2_arith_mul.h b/src/core/operators/codac2_arith_mul.h index 14fd49447..62cf68cf1 100644 --- a/src/core/operators/codac2_arith_mul.h +++ b/src/core/operators/codac2_arith_mul.h @@ -52,17 +52,17 @@ namespace codac2 static Interval fwd(const Interval& x1, const Interval& x2); static ScalarType fwd_natural(const ScalarType& x1, const ScalarType& x2); - static ScalarType fwd_centered(const ScalarType& x1, const ScalarType& x2); + static ScalarType fwd_centered(bool mixed_diff,const ScalarType& x1, const ScalarType& x2); static void bwd(const Interval& y, Interval& x1, Interval& x2); static IntervalVector fwd(const Interval& x1, const IntervalVector& x2); static VectorType fwd_natural(const ScalarType& x1, const VectorType& x2); - static VectorType fwd_centered(const ScalarType& x1, const VectorType& x2); + static VectorType fwd_centered(bool mixed_diff,const ScalarType& x1, const VectorType& x2); static void bwd(const IntervalVector& y, Interval& x1, IntervalVector& x2); static IntervalVector fwd(const IntervalVector& x1, const Interval& x2); static VectorType fwd_natural(const VectorType& x1, const ScalarType& x2); - static VectorType fwd_centered(const VectorType& x1, const ScalarType& x2); + static VectorType fwd_centered(bool mixed_diff,const VectorType& x1, const ScalarType& x2); static void bwd(const IntervalVector& y, IntervalVector& x1, Interval& x2); static Interval fwd(const IntervalRow& x1, const IntervalVector& x2); @@ -71,17 +71,17 @@ namespace codac2 static IntervalMatrix fwd(const Interval& x1, const IntervalMatrix& x2); static MatrixType fwd_natural(const ScalarType& x1, const MatrixType& x2); - static MatrixType fwd_centered(const ScalarType& x1, const MatrixType& x2); + static MatrixType fwd_centered(bool mixed_diff,const ScalarType& x1, const MatrixType& x2); static void bwd(const IntervalMatrix& y, Interval& x1, IntervalMatrix& x2); static IntervalVector fwd(const IntervalMatrix& x1, const IntervalVector& x2); static VectorType fwd_natural(const MatrixType& x1, const VectorType& x2); - static VectorType fwd_centered(const MatrixType& x1, const VectorType& x2); + static VectorType fwd_centered(bool mixed_diff,const MatrixType& x1, const VectorType& x2); static void bwd(const IntervalVector& y, IntervalMatrix& x1, IntervalVector& x2); static IntervalMatrix fwd(const IntervalMatrix& x1, const IntervalMatrix& x2); static MatrixType fwd_natural(const MatrixType& x1, const MatrixType& x2); - static MatrixType fwd_centered(const MatrixType& x1, const MatrixType& x2); + static MatrixType fwd_centered(bool mixed_diff,const MatrixType& x1, const MatrixType& x2); static void bwd(const IntervalMatrix& y, IntervalMatrix& x1, IntervalMatrix& x2); }; @@ -134,12 +134,12 @@ namespace codac2 inline ScalarType MulOp::fwd_natural(const ScalarType& x1, const ScalarType& x2) { return { - fwd(x1.a, x2.a), + fwd(x1.get_a(), x2.get_a()), x1.def_domain && x2.def_domain }; } - inline ScalarType MulOp::fwd_centered(const ScalarType& x1, const ScalarType& x2) + inline ScalarType MulOp::fwd_centered(bool mixed_diff, const ScalarType& x1, const ScalarType& x2) { if(centered_form_not_available_for_args(x1,x2)) return fwd_natural(x1,x2); @@ -149,11 +149,16 @@ namespace codac2 IntervalMatrix d(1,x1.da.cols()); for(Index i = 0 ; i < d.size() ; i++) - d(0,i) = x1.da(0,i)*x2.a + x1.a*x2.da(0,i); + { + const Interval &x1a = x1.get_a(mixed_diff,i); + const Interval &x2a = x2.get_a(mixed_diff,i); + d(0,i) = x1.da(0,i)*x2a + x1a*x2.da(0,i); + } return { fwd(x1.m, x2.m), - fwd(x1.a, x2.a), + vector_map(x1.a,x2.a, + static_cast(&fwd)), d, x1.def_domain && x2.def_domain }; @@ -173,28 +178,32 @@ namespace codac2 inline VectorType MulOp::fwd_natural(const ScalarType& x1, const VectorType& x2) { return { - fwd(x1.a, x2.a), + fwd(x1.get_a(), x2.get_a()), x1.def_domain && x2.def_domain }; } - inline VectorType MulOp::fwd_centered(const ScalarType& x1, const VectorType& x2) + inline VectorType MulOp::fwd_centered(bool mixed_diff,const ScalarType& x1, const VectorType& x2) { if(centered_form_not_available_for_args(x1,x2)) return fwd_natural(x1,x2); assert(x1.da.rows() == 1); assert(x1.da.cols() == x2.da.cols()); - assert(x2.a.size() == x2.da.rows()); + assert(x2.get_a().size() == x2.da.rows()); IntervalMatrix d(x2.da.rows(),x2.da.cols()); for(Index i = 0 ; i < d.rows() ; i++) - for(Index j = 0 ; j < d.cols() ; j++) - d(i,j) = x1.da(0,j)*x2.a[i]+x1.a*x2.da(i,j); + for(Index j = 0 ; j < d.cols() ; j++) { + const Interval& x1a = x1.get_a(mixed_diff,j); + const IntervalVector& x2a = x2.get_a(mixed_diff,j); + d(i,j) = x1.da(0,j)*x2a[i]+x1a*x2.da(i,j); + } return { fwd(x1.m, x2.m), - fwd(x1.a, x2.a), + vector_map(x1.a,x2.a, + static_cast(&fwd)), d, x1.def_domain && x2.def_domain }; @@ -217,9 +226,9 @@ namespace codac2 return MulOp::fwd_natural(x2,x1); } - inline VectorType MulOp::fwd_centered(const VectorType& x1, const ScalarType& x2) + inline VectorType MulOp::fwd_centered(bool mixed_diff,const VectorType& x1, const ScalarType& x2) { - return MulOp::fwd_centered(x2,x1); + return MulOp::fwd_centered(mixed_diff,x2,x1); } inline void MulOp::bwd(const IntervalVector& y, IntervalVector& x1, Interval& x2) @@ -249,26 +258,30 @@ namespace codac2 inline MatrixType MulOp::fwd_natural(const ScalarType& x1, const MatrixType& x2) { return { - fwd(x1.a, x2.a), + fwd(x1.get_a(), x2.get_a()), x1.def_domain && x2.def_domain }; } - inline MatrixType MulOp::fwd_centered(const ScalarType& x1, const MatrixType& x2) + inline MatrixType MulOp::fwd_centered(bool mixed_diff,const ScalarType& x1, const MatrixType& x2) { if(centered_form_not_available_for_args(x1,x2)) return fwd_natural(x1,x2); assert(x2.da.cols() == x1.da.cols()); IntervalMatrix d(x2.da.rows(),x2.da.cols()); - for (Index j=0; j()[i]+x1.a*x2.da(i,j); + d(i,j) = x1.da(0,j)*x2a.reshaped()[i]+x1a*x2.da(i,j); + } } return { fwd(x1.m, x2.m), - fwd(x1.a, x2.a), + vector_map(x1.a,x2.a, + static_cast(&fwd)), d, x1.def_domain && x2.def_domain }; @@ -292,28 +305,32 @@ namespace codac2 inline VectorType MulOp::fwd_natural(const MatrixType& x1, const VectorType& x2) { return { - fwd(x1.a, x2.a), + fwd(x1.get_a(), x2.get_a()), x1.def_domain && x2.def_domain }; } - inline VectorType MulOp::fwd_centered(const MatrixType& x1, const VectorType& x2) + inline VectorType MulOp::fwd_centered(bool mixed_diff,const MatrixType& x1, const VectorType& x2) { if(centered_form_not_available_for_args(x1,x2)) return fwd_natural(x1,x2); assert(x2.da.cols() == x1.da.cols()); - IntervalMatrix d = IntervalMatrix::zero(x1.a.rows(),x1.da.cols()); - for (Index j=0; j(x1.a,x2.a, + static_cast(&fwd)), d, x1.def_domain && x2.def_domain }; @@ -328,32 +345,36 @@ namespace codac2 inline MatrixType MulOp::fwd_natural(const MatrixType& x1, const MatrixType& x2) { return { - fwd(x1.a, x2.a), + fwd(x1.get_a(), x2.get_a()), x1.def_domain && x2.def_domain }; } - inline MatrixType MulOp::fwd_centered(const MatrixType& x1, const MatrixType& x2) + inline MatrixType MulOp::fwd_centered(bool mixed_diff,const MatrixType& x1, const MatrixType& x2) { if(centered_form_not_available_for_args(x1,x2)) return fwd_natural(x1,x2); assert(x2.da.cols() == x1.da.cols()); - IntervalMatrix d = IntervalMatrix::zero(x1.a.rows()*x2.a.cols(),x1.da.cols()); - for (Index j=0; j(x1.a,x2.a, + static_cast(&fwd)), d, x1.def_domain && x2.def_domain }; diff --git a/src/core/operators/codac2_arith_sub.h b/src/core/operators/codac2_arith_sub.h index afe6e441b..523b4a65f 100644 --- a/src/core/operators/codac2_arith_sub.h +++ b/src/core/operators/codac2_arith_sub.h @@ -34,17 +34,17 @@ namespace codac2 static Interval fwd(const Interval& x1); static ScalarType fwd_natural(const ScalarType& x1); - static ScalarType fwd_centered(const ScalarType& x1); + static ScalarType fwd_centered(bool mixed_diff,const ScalarType& x1); static void bwd(const Interval& y, Interval& x1); static IntervalVector fwd(const IntervalVector& x1); static VectorType fwd_natural(const VectorType& x1); - static VectorType fwd_centered(const VectorType& x1); + static VectorType fwd_centered(bool mixed_diff,const VectorType& x1); static void bwd(const IntervalVector& y, IntervalVector& x1); static IntervalMatrix fwd(const IntervalMatrix& x1); static MatrixType fwd_natural(const MatrixType& x1); - static MatrixType fwd_centered(const MatrixType& x1); + static MatrixType fwd_centered(bool mixed_diff,const MatrixType& x1); static void bwd(const IntervalMatrix& y, IntervalMatrix& x1); // Binary operations @@ -65,17 +65,17 @@ namespace codac2 static Interval fwd(const Interval& x1, const Interval& x2); static ScalarType fwd_natural(const ScalarType& x1, const ScalarType& x2); - static ScalarType fwd_centered(const ScalarType& x1, const ScalarType& x2); + static ScalarType fwd_centered(bool mixed_diff,const ScalarType& x1, const ScalarType& x2); static void bwd(const Interval& y, Interval& x1, Interval& x2); static IntervalVector fwd(const IntervalVector& x1, const IntervalVector& x2); static VectorType fwd_natural(const VectorType& x1, const VectorType& x2); - static VectorType fwd_centered(const VectorType& x1, const VectorType& x2); + static VectorType fwd_centered(bool mixed_diff,const VectorType& x1, const VectorType& x2); static void bwd(const IntervalVector& y, IntervalVector& x1, IntervalVector& x2); static IntervalMatrix fwd(const IntervalMatrix& x1, const IntervalMatrix& x2); static MatrixType fwd_natural(const MatrixType& x1, const MatrixType& x2); - static MatrixType fwd_centered(const MatrixType& x1, const MatrixType& x2); + static MatrixType fwd_centered(bool mixed_diff,const MatrixType& x1, const MatrixType& x2); static void bwd(const IntervalMatrix& y, IntervalMatrix& x1, IntervalMatrix& x2); }; @@ -131,16 +131,17 @@ namespace codac2 inline ScalarType SubOp::fwd_natural(const ScalarType& x1) { return { - fwd(x1.a), + fwd(x1.get_a()), x1.def_domain }; } - inline ScalarType SubOp::fwd_centered(const ScalarType& x1) + inline ScalarType SubOp::fwd_centered([[maybe_unused]] bool mixed_diff,const ScalarType& x1) { return { fwd(x1.m), - fwd(x1.a), + vector_map(x1.a, + static_cast(&fwd)), -x1.da, x1.def_domain }; @@ -160,16 +161,17 @@ namespace codac2 inline VectorType SubOp::fwd_natural(const VectorType& x1) { return { - fwd(x1.a), + fwd(x1.get_a()), x1.def_domain }; } - inline VectorType SubOp::fwd_centered(const VectorType& x1) + inline VectorType SubOp::fwd_centered([[maybe_unused]] bool mixed_diff,const VectorType& x1) { return { fwd(x1.m), - fwd(x1.a), + vector_map(x1.a, + static_cast(&fwd)), -x1.da, x1.def_domain }; @@ -190,16 +192,17 @@ namespace codac2 inline MatrixType SubOp::fwd_natural(const MatrixType& x1) { return { - fwd(x1.a), + fwd(x1.get_a()), x1.def_domain }; } - inline MatrixType SubOp::fwd_centered(const MatrixType& x1) + inline MatrixType SubOp::fwd_centered([[maybe_unused]] bool mixed_diff, const MatrixType& x1) { return { fwd(x1.m), - fwd(x1.a), + vector_map(x1.a, + static_cast(&fwd)), -x1.da, x1.def_domain }; @@ -220,12 +223,12 @@ namespace codac2 inline ScalarType SubOp::fwd_natural(const ScalarType& x1, const ScalarType& x2) { return { - fwd(x1.a, x2.a), + fwd(x1.get_a(), x2.get_a()), x1.def_domain && x2.def_domain }; } - inline ScalarType SubOp::fwd_centered(const ScalarType& x1, const ScalarType& x2) + inline ScalarType SubOp::fwd_centered([[maybe_unused]] bool mixed_diff, const ScalarType& x1, const ScalarType& x2) { if(centered_form_not_available_for_args(x1,x2)) return fwd_natural(x1,x2); @@ -233,7 +236,8 @@ namespace codac2 assert(x1.da.rows() == x2.da.rows() && x1.da.cols() == x2.da.cols()); return { fwd(x1.m, x2.m), - fwd(x1.a, x2.a), + vector_map(x1.a,x2.a, + static_cast(&fwd)), x1.da - x2.da, x1.def_domain && x2.def_domain }; @@ -257,12 +261,12 @@ namespace codac2 inline VectorType SubOp::fwd_natural(const VectorType& x1, const VectorType& x2) { return { - fwd(x1.a, x2.a), + fwd(x1.get_a(), x2.get_a()), x1.def_domain && x2.def_domain }; } - inline VectorType SubOp::fwd_centered(const VectorType& x1, const VectorType& x2) + inline VectorType SubOp::fwd_centered([[maybe_unused]] bool mixed_diff,const VectorType& x1, const VectorType& x2) { if(centered_form_not_available_for_args(x1,x2)) return fwd_natural(x1,x2); @@ -270,7 +274,8 @@ namespace codac2 assert(x1.da.rows() == x2.da.rows() && x1.da.cols() == x2.da.cols()); return { fwd(x1.m, x2.m), - fwd(x1.a, x2.a), + vector_map(x1.a,x2.a, + static_cast(&fwd)), x1.da - x2.da, x1.def_domain && x2.def_domain }; @@ -293,17 +298,18 @@ namespace codac2 { assert(x1.a.cols() == x2.a.cols() && x1.a.rows() == x2.a.rows()); return { - fwd(x1.a, x2.a), + fwd(x1.get_a(), x2.get_a()), x1.def_domain && x2.def_domain }; } - inline MatrixType SubOp::fwd_centered(const MatrixType& x1, const MatrixType& x2) + inline MatrixType SubOp::fwd_centered([[maybe_unused]] bool mixed_diff,const MatrixType& x1, const MatrixType& x2) { - assert(x1.a.cols() == x2.a.cols() && x1.a.rows() == x2.a.rows()); + assert(x1.get_a().cols() == x2.get_a().cols() && x1.get_a().rows() == x2.get_a().rows()); return { fwd(x1.m, x2.m), - fwd(x1.a, x2.a), + vector_map(x1.a,x2.a, + static_cast(&fwd)), x1.da - x2.da, x1.def_domain && x2.def_domain }; diff --git a/src/core/operators/codac2_asin.h b/src/core/operators/codac2_asin.h index a3e4bdc68..6f58c8a25 100644 --- a/src/core/operators/codac2_asin.h +++ b/src/core/operators/codac2_asin.h @@ -31,7 +31,7 @@ namespace codac2 static Interval fwd(const Interval& x1); static ScalarType fwd_natural(const ScalarType& x1); - static ScalarType fwd_centered(const ScalarType& x1); + static ScalarType fwd_centered(bool mixed_diff, const ScalarType& x1); static void bwd(const Interval& y, Interval& x1); }; @@ -54,28 +54,30 @@ namespace codac2 inline ScalarType AsinOp::fwd_natural(const ScalarType& x1) { return { - fwd(x1.a), - x1.a.is_subset({-1,1}) // def domain of asin - && !x1.a.contains(1.) // def domain of the derivative of asin + fwd(x1.get_a()), + x1.get_a().is_subset({-1,1}) // def domain of asin + && !x1.get_a().contains(1.) // def domain of the derivative of asin && x1.def_domain }; } - inline ScalarType AsinOp::fwd_centered(const ScalarType& x1) + inline ScalarType AsinOp::fwd_centered(bool mixed_diff,const ScalarType& x1) { if(centered_form_not_available_for_args(x1)) return fwd_natural(x1); IntervalMatrix d(1,x1.da.size()); - for(Index i = 0 ; i < d.size() ; i++) - d(0,i) = x1.da(0,i)/sqrt(1.-sqr(x1.a)); + for(Index i = 0 ; i < d.size() ; i++) { + const Interval& x1a = x1.get_a(mixed_diff,i); + d(0,i) = x1.da(0,i)/sqrt(1.-sqr(x1a)); + } return { fwd(x1.m), - fwd(x1.a), + vector_map(x1.a,&fwd), d, - x1.a.is_subset({-1,1}) // def domain of asin - && x1.a != 1. // def domain of the derivative of asin + x1.get_a().is_subset({-1,1}) // def domain of asin + && !x1.get_a().contains(1.) // def domain of the derivative of asin && x1.def_domain }; } diff --git a/src/core/operators/codac2_atan.h b/src/core/operators/codac2_atan.h index 24280986d..f1ebaaecb 100644 --- a/src/core/operators/codac2_atan.h +++ b/src/core/operators/codac2_atan.h @@ -31,7 +31,7 @@ namespace codac2 static Interval fwd(const Interval& x1); static ScalarType fwd_natural(const ScalarType& x1); - static ScalarType fwd_centered(const ScalarType& x1); + static ScalarType fwd_centered(bool mixed_diff,const ScalarType& x1); static void bwd(const Interval& y, Interval& x1); }; @@ -54,23 +54,25 @@ namespace codac2 inline ScalarType AtanOp::fwd_natural(const ScalarType& x1) { return { - fwd(x1.a), + fwd(x1.get_a()), x1.def_domain }; } - inline ScalarType AtanOp::fwd_centered(const ScalarType& x1) + inline ScalarType AtanOp::fwd_centered(bool mixed_diff,const ScalarType& x1) { if(centered_form_not_available_for_args(x1)) return fwd_natural(x1); IntervalMatrix d(1,x1.da.size()); - for(Index i = 0 ; i < d.size() ; i++) - d(0,i) = x1.da(0,i)/(1.+sqr(x1.a)); + for(Index i = 0 ; i < d.size() ; i++) { + const Interval& x1a = x1.get_a(mixed_diff,i); + d(0,i) = x1.da(0,i)/(1.+sqr(x1a)); + } return { fwd(x1.m), - fwd(x1.a), + vector_map(x1.a,&fwd), d, x1.def_domain }; diff --git a/src/core/operators/codac2_atan2.h b/src/core/operators/codac2_atan2.h index 199d6c3f9..5517b55d7 100644 --- a/src/core/operators/codac2_atan2.h +++ b/src/core/operators/codac2_atan2.h @@ -33,7 +33,7 @@ namespace codac2 static Interval fwd(const Interval& x1, const Interval& x2); static ScalarType fwd_natural(const ScalarType& x1, const ScalarType& x2); - static ScalarType fwd_centered(const ScalarType& x1, const ScalarType& x2); + static ScalarType fwd_centered(bool mixed_diff, const ScalarType& x1, const ScalarType& x2); static void bwd(const Interval& y, Interval& x1, Interval& x2); }; @@ -56,13 +56,13 @@ namespace codac2 inline ScalarType Atan2Op::fwd_natural(const ScalarType& x1, const ScalarType& x2) { return { - fwd(x1.a, x2.a), + fwd(x1.get_a(), x2.get_a()), x1.def_domain && x2.def_domain - && !(x1.a.contains(0.) && x2.a.contains(0.)) // def domain of the derivative of atan2 + && !(x1.get_a().contains(0.) && x2.get_a().contains(0.)) // def domain of the derivative of atan2 }; } - inline ScalarType Atan2Op::fwd_centered(const ScalarType& x1, const ScalarType& x2) + inline ScalarType Atan2Op::fwd_centered(bool mixed_diff, const ScalarType& x1, const ScalarType& x2) { if(centered_form_not_available_for_args(x1,x2)) return fwd_natural(x1,x2); @@ -71,15 +71,17 @@ namespace codac2 assert(x1.da.rows() == x2.da.rows() && x1.da.cols() == x2.da.cols()); IntervalMatrix d(1,x1.da.cols()); - for(Index i = 0 ; i < d.size() ; i++) - d(0,i) = (-x1.a*x2.da(0,i)/(sqr(x2.a)+sqr(x1.a)))+(x2.a*x1.da(0,i)/(sqr(x2.a)+sqr(x1.a))); - + for(Index i = 0 ; i < d.size() ; i++) { + const Interval& x1a = x1.get_a(mixed_diff,i); + const Interval& x2a = x2.get_a(mixed_diff,i); + d(0,i) = (-x1a*x2.da(0,i)/(sqr(x2a)+sqr(x1a)))+(x2a*x1.da(0,i)/(sqr(x2a)+sqr(x1a))); + } return { fwd(x1.m, x2.m), - fwd(x1.a, x2.a), + vector_map(x1.a,x2.a,&fwd), d, x1.def_domain && x2.def_domain - && !(x1.a.contains(0.) && x2.a.contains(0.)) // def domain of the derivative of atan2 + && !(x1.get_a().contains(0.) && x2.get_a().contains(0.)) // def domain of the derivative of atan2 }; } diff --git a/src/core/operators/codac2_ceil.h b/src/core/operators/codac2_ceil.h index df5866ab3..f5bdde8f7 100644 --- a/src/core/operators/codac2_ceil.h +++ b/src/core/operators/codac2_ceil.h @@ -31,7 +31,7 @@ namespace codac2 static Interval fwd(const Interval& x1); static ScalarType fwd_natural(const ScalarType& x1); - static ScalarType fwd_centered(const ScalarType& x1); + static ScalarType fwd_centered(bool mixed_diff,const ScalarType& x1); static void bwd(const Interval& y, Interval& x1); }; @@ -54,16 +54,16 @@ namespace codac2 inline ScalarType CeilOp::fwd_natural(const ScalarType& x1) { return { - fwd(x1.a), + fwd(x1.get_a()), x1.def_domain }; } - inline ScalarType CeilOp::fwd_centered(const ScalarType& x1) + inline ScalarType CeilOp::fwd_centered([[maybe_unused]] bool mixed_diff, const ScalarType& x1) { return { fwd(x1.m), - fwd(x1.a), + vector_map(x1.a,&fwd), IntervalMatrix(0,0), // not supported yet for auto diff x1.def_domain }; diff --git a/src/core/operators/codac2_chi.h b/src/core/operators/codac2_chi.h index 066e4e792..1e54d221c 100644 --- a/src/core/operators/codac2_chi.h +++ b/src/core/operators/codac2_chi.h @@ -50,28 +50,37 @@ namespace codac2 static inline T fwd_natural(const ScalarType& x1, const T& x2, const T& x3) { return { - fwd(x1.a,x2.a,x3.a), + fwd(x1.get_a(),x2.get_a(),x3.get_a()), x1.def_domain - && ((x1.a.lb()>0.0) ? x2.def_domain : true) - && ((x1.a.ub()<=0.0) ? x3.def_domain : true) + && ((x1.get_a().lb()>0.0) ? x2.def_domain : true) + && ((x1.get_a().ub()<=0.0) ? x3.def_domain : true) }; } template - static inline T fwd_centered(const ScalarType& x1, const T& x2, const T& x3) + static inline T fwd_centered(bool mixed_diff,const ScalarType& x1, const T& x2, const T& x3) { if(centered_form_not_available_for_args(x1,x2,x3)) return fwd_natural(x1,x2,x3); assert(x2.da.rows() == x3.da.rows()); + assert(x2.da.cols() == x3.da.cols()); + + IntervalMatrix d(x3.da.rows(),x3.da.cols()); + for(Index j = 0 ; j < d.cols() ; j++) { + const Interval& x1a = x1.get_a(mixed_diff,j); + if (x1a.lb()>0) d.col(j) = x3.da.col(j); + else if (x1a.ub()<=0) d.col(j) = x2.da.col(j); + else d.col(j) = (IntervalVector(x2.da.col(j)) | IntervalVector(x3.da.col(j))); + } return { fwd(x1.m,x2.m,x3.m), - fwd(x1.a,x2.a,x3.a), - (x1.a.lb()>0 ? x3.da : ((x1.a.ub()<=0) ? x2.da : (x2.da | x3.da))), + vector_map(x1.a,x2.a,x3.a,&fwd), + d, x1.def_domain - && ((x1.a.lb()>0.0) ? x2.def_domain : true) - && ((x1.a.ub()<=0.0) ? x3.def_domain : true) + && ((x1.get_a().lb()>0.0) ? x2.def_domain : true) + && ((x1.get_a().ub()<=0.0) ? x3.def_domain : true) }; } diff --git a/src/core/operators/codac2_component.h b/src/core/operators/codac2_component.h index 5558209b0..8db9a80f6 100644 --- a/src/core/operators/codac2_component.h +++ b/src/core/operators/codac2_component.h @@ -47,12 +47,12 @@ namespace codac2 static Interval fwd(const IntervalVector& x1, Index i); static ScalarType fwd_natural(const VectorType& x1, Index i); - static ScalarType fwd_centered(const VectorType& x1, Index i); + static ScalarType fwd_centered(bool mixed_diff,const VectorType& x1, Index i); static void bwd(const Interval& y, IntervalVector& x1, Index i); static Interval fwd(const IntervalMatrix& x1, Index i, Index j); static ScalarType fwd_natural(const MatrixType& x1, Index i, Index j); - static ScalarType fwd_centered(const MatrixType& x1, Index i, Index j); + static ScalarType fwd_centered(bool mixed_diff,const MatrixType& x1, Index i, Index j); static void bwd(const Interval& y, IntervalMatrix& x1, Index i, Index j); }; @@ -81,19 +81,22 @@ namespace codac2 return OperationExprBase>::replace_arg(old_arg_id, new_expr); } - ScalarType fwd_eval(ValuesMap& v, Index total_input_size, bool natural_eval) const override + ScalarType fwd_eval(ValuesMap& v, Index total_input_size, bool natural_eval, bool mixed_diff) const override { if(natural_eval) return AnalyticExpr::init_value( - v, ComponentOp::fwd_natural(std::get<0>(this->_x)->fwd_eval(v, total_input_size, natural_eval), _i)); - else + v, ComponentOp::fwd_natural(std::get<0>(this->_x)->fwd_eval(v, total_input_size, natural_eval,mixed_diff), _i)); + else if (!mixed_diff) + return AnalyticExpr::init_value( + v, ComponentOp::fwd_centered(false,std::get<0>(this->_x)->fwd_eval(v, total_input_size, natural_eval,mixed_diff), _i)); + else return AnalyticExpr::init_value( - v, ComponentOp::fwd_centered(std::get<0>(this->_x)->fwd_eval(v, total_input_size, natural_eval), _i)); + v, ComponentOp::fwd_centered(true,std::get<0>(this->_x)->fwd_eval(v, total_input_size, natural_eval,mixed_diff), _i)); } void bwd_eval(ValuesMap& v) const override { - ComponentOp::bwd(AnalyticExpr::value(v).a, std::get<0>(this->_x)->value(v).a, _i); + ComponentOp::bwd(AnalyticExpr::value(v).get_a(), std::get<0>(this->_x)->value(v).get_a(), _i); std::get<0>(this->_x)->bwd_eval(v); } @@ -156,19 +159,19 @@ namespace codac2 return OperationExprBase>::replace_arg(old_arg_id, new_expr); } - ScalarType fwd_eval(ValuesMap& v, Index total_input_size, bool natural_eval) const override + ScalarType fwd_eval(ValuesMap& v, Index total_input_size, bool natural_eval,bool mixed_diff) const override { if(natural_eval) return AnalyticExpr::init_value( - v, ComponentOp::fwd_natural(std::get<0>(this->_x)->fwd_eval(v, total_input_size, natural_eval), _i, _j)); + v, ComponentOp::fwd_natural(std::get<0>(this->_x)->fwd_eval(v, total_input_size, natural_eval,mixed_diff), _i, _j)); else return AnalyticExpr::init_value( - v, ComponentOp::fwd_centered(std::get<0>(this->_x)->fwd_eval(v, total_input_size, natural_eval), _i, _j)); + v, ComponentOp::fwd_centered(mixed_diff,std::get<0>(this->_x)->fwd_eval(v, total_input_size, natural_eval,mixed_diff), _i, _j)); } void bwd_eval(ValuesMap& v) const override { - ComponentOp::bwd(AnalyticExpr::value(v).a, std::get<0>(this->_x)->value(v).a, _i, _j); + ComponentOp::bwd(AnalyticExpr::value(v).get_a(), std::get<0>(this->_x)->value(v).get_a(), _i, _j); std::get<0>(this->_x)->bwd_eval(v); } @@ -225,22 +228,21 @@ namespace codac2 inline ScalarType ComponentOp::fwd_natural(const VectorType& x1, Index i) { - assert(i >= 0 && i < x1.a.rows()); return { - fwd(x1.a,i), + fwd(x1.get_a(),i), x1.def_domain }; } - inline ScalarType ComponentOp::fwd_centered(const VectorType& x1, Index i) + inline ScalarType ComponentOp::fwd_centered([[maybe_unused]] bool mixed_diff,const VectorType& x1, Index i) { if(centered_form_not_available_for_args(x1)) return fwd_natural(x1,i); - assert(i >= 0 && i < x1.a.rows()); return { fwd(x1.m,i), - fwd(x1.a,i), + vector_map(x1.a, + [=](const IntervalVector& v1)->Interval {return fwd(v1,i);}), x1.da.row(i), x1.def_domain }; @@ -263,25 +265,26 @@ namespace codac2 inline ScalarType ComponentOp::fwd_natural(const MatrixType& x1, Index i, Index j) { - assert(i >= 0 && i < x1.a.rows()); - assert(j >= 0 && j < x1.a.cols()); + assert(i >= 0 && i < x1.get_a().rows()); + assert(j >= 0 && j < x1.get_a().cols()); return { - fwd(x1.a,i,j), + fwd(x1.get_a(),i,j), x1.def_domain }; } - inline ScalarType ComponentOp::fwd_centered(const MatrixType& x1, Index i, Index j) + inline ScalarType ComponentOp::fwd_centered([[maybe_unused]] bool mixed_diff,const MatrixType& x1, Index i, Index j) { if(centered_form_not_available_for_args(x1)) return fwd_natural(x1,i,j); - assert(i >= 0 && i < x1.a.rows()); - assert(j >= 0 && j < x1.a.cols()); + assert(i >= 0 && i < x1.get_a().rows()); + assert(j >= 0 && j < x1.get_a().cols()); return { fwd(x1.m,i,j), - fwd(x1.a,i,j), - x1.da.row(x1.a.rows()*i+j), // centered form, ColMajor + vector_map(x1.a, + [=](const IntervalMatrix& v1)->Interval {return fwd(v1,i,j);}), + x1.da.row(x1.get_a().rows()*i+j), // centered form, ColMajor x1.def_domain }; } diff --git a/src/core/operators/codac2_cos.h b/src/core/operators/codac2_cos.h index 4389bfeb0..fc52ff71a 100644 --- a/src/core/operators/codac2_cos.h +++ b/src/core/operators/codac2_cos.h @@ -31,7 +31,7 @@ namespace codac2 static Interval fwd(const Interval& x1); static ScalarType fwd_natural(const ScalarType& x1); - static ScalarType fwd_centered(const ScalarType& x1); + static ScalarType fwd_centered(bool mixed_diff, const ScalarType& x1); static void bwd(const Interval& y, Interval& x1); }; @@ -54,23 +54,25 @@ namespace codac2 inline ScalarType CosOp::fwd_natural(const ScalarType& x1) { return { - fwd(x1.a), + fwd(x1.get_a()), x1.def_domain }; } - inline ScalarType CosOp::fwd_centered(const ScalarType& x1) + inline ScalarType CosOp::fwd_centered(bool mixed_diff, const ScalarType& x1) { if(centered_form_not_available_for_args(x1)) return fwd_natural(x1); IntervalMatrix d(1,x1.da.size()); - for(Index i = 0 ; i < d.size() ; i++) - d(0,i) = -sin(x1.a)*x1.da(0,i); + for(Index i = 0 ; i < d.size() ; i++) { + const Interval& x1a = x1.get_a(mixed_diff,i); + d(0,i) = -sin(x1a)*x1.da(0,i); + } return { fwd(x1.m), - fwd(x1.a), + vector_map(x1.a,&fwd), d, x1.def_domain }; diff --git a/src/core/operators/codac2_cosh.h b/src/core/operators/codac2_cosh.h index cd128f0e2..0a67bc2b4 100644 --- a/src/core/operators/codac2_cosh.h +++ b/src/core/operators/codac2_cosh.h @@ -31,7 +31,7 @@ namespace codac2 static Interval fwd(const Interval& x1); static ScalarType fwd_natural(const ScalarType& x1); - static ScalarType fwd_centered(const ScalarType& x1); + static ScalarType fwd_centered(bool mixed_diff, const ScalarType& x1); static void bwd(const Interval& y, Interval& x1); }; @@ -54,23 +54,25 @@ namespace codac2 inline ScalarType CoshOp::fwd_natural(const ScalarType& x1) { return { - fwd(x1.a), + fwd(x1.get_a()), x1.def_domain }; } - inline ScalarType CoshOp::fwd_centered(const ScalarType& x1) + inline ScalarType CoshOp::fwd_centered(bool mixed_diff, const ScalarType& x1) { if(centered_form_not_available_for_args(x1)) return fwd_natural(x1); IntervalMatrix d(1,x1.da.size()); - for(Index i = 0 ; i < d.size() ; i++) - d(0,i) = sinh(x1.a)*x1.da(0,i); + for(Index i = 0 ; i < d.size() ; i++) { + const Interval& x1a = x1.get_a(mixed_diff,i); + d(0,i) = sinh(x1a)*x1.da(0,i); + } return { fwd(x1.m), - fwd(x1.a), + vector_map(x1.a,&fwd), d, x1.def_domain }; diff --git a/src/core/operators/codac2_cross_prod.h b/src/core/operators/codac2_cross_prod.h index 824dc2af0..11f5cfcd4 100644 --- a/src/core/operators/codac2_cross_prod.h +++ b/src/core/operators/codac2_cross_prod.h @@ -12,6 +12,7 @@ #include "codac2_Interval.h" #include "codac2_AnalyticType.h" #include "codac2_AnalyticExprWrapper.h" +#include "codac2_arith_mul.h" namespace codac2 { @@ -31,7 +32,7 @@ namespace codac2 static IntervalVector fwd(const IntervalVector& x1, const IntervalVector& x2); static VectorType fwd_natural(const VectorType& x1, const VectorType& x2); - static VectorType fwd_centered(const VectorType& x1, const VectorType& x2); + static VectorType fwd_centered(bool mixed_diff,const VectorType& x1, const VectorType& x2); static void bwd(const IntervalVector& y, IntervalVector& x1, IntervalVector& x2); }; @@ -59,12 +60,12 @@ namespace codac2 inline VectorType CrossProdOp::fwd_natural(const VectorType& x1, const VectorType& x2) { return { - fwd(x1.a,x2.a), + fwd(x1.get_a(),x2.get_a()), x1.def_domain && x2.def_domain }; } - inline VectorType CrossProdOp::fwd_centered(const VectorType& x1, const VectorType& x2) + inline VectorType CrossProdOp::fwd_centered(bool mixed_diff,const VectorType& x1, const VectorType& x2) { if(centered_form_not_available_for_args(x1,x2)) return fwd_natural(x1,x2); @@ -72,25 +73,38 @@ namespace codac2 assert(x1.da.cols() == x2.da.cols()); IntervalMatrix d(3,x1.da.cols()); for (Index i=0;i(x1.a,x2.a,&fwd), d, x1.def_domain && x2.def_domain }; } - inline void CrossProdOp::bwd([[maybe_unused]] const IntervalVector& y, [[maybe_unused]] IntervalVector& x1, [[maybe_unused]] IntervalVector& x2) + inline void CrossProdOp::bwd(const IntervalVector& y, IntervalVector& x1, IntervalVector& x2) { assert_release(y.size() == 3 && x1.size() == 3 && x2.size() == 3); - // todo + IntervalRow u { 0.0, x2[2], -x2[1] }; + MulOp::bwd(y[0],u,x1); + if (x1.is_empty()) { x2.set_empty(); return; } + x2[2] &= u[1]; x2[1] &= -u[2]; + IntervalRow v { -x2[2], 0.0, x2[0] }; + MulOp::bwd(y[1],v,x1); + if (x1.is_empty()) { x2.set_empty(); return; } + x2[2] &= -v[0]; x2[0] &= v[2]; + IntervalRow w { x2[1], -x2[0], 0.0 }; + MulOp::bwd(y[2],w,x1); + if (x1.is_empty()) { x2.set_empty(); return; } + x2[1] &= w[0]; x2[0] &= -w[1]; } } diff --git a/src/core/operators/codac2_det.h b/src/core/operators/codac2_det.h index c11fb630b..83b99c4f0 100644 --- a/src/core/operators/codac2_det.h +++ b/src/core/operators/codac2_det.h @@ -17,6 +17,7 @@ #include "codac2_arith_sub.h" #include "codac2_arith_mul.h" #include "codac2_inversion.h" +#include namespace codac2 { @@ -49,7 +50,7 @@ namespace codac2 static Interval fwd(const IntervalMatrix& x); static ScalarType fwd_natural(const MatrixType& x); - static ScalarType fwd_centered(const MatrixType& x); + static ScalarType fwd_centered(bool mixed_diff,const MatrixType& x); static void bwd(const Interval& y, IntervalMatrix& x); // For two vectors (merged into a 2×2 matrix) @@ -61,7 +62,7 @@ namespace codac2 static Interval fwd(const IntervalVector& x1, const IntervalVector& x2); static ScalarType fwd_natural(const VectorType& x1, const VectorType& x2); - static ScalarType fwd_centered(const VectorType& x1, const VectorType& x2); + static ScalarType fwd_centered(bool mixed_diff,const VectorType& x1, const VectorType& x2); static void bwd(const Interval& y, IntervalVector& x1, IntervalVector& x2); // For three vectors (merged into a 3×3 matrix) @@ -73,7 +74,7 @@ namespace codac2 static Interval fwd(const IntervalVector& x1, const IntervalVector& x2, const IntervalVector& x3); static ScalarType fwd_natural(const VectorType& x1, const VectorType& x2, const VectorType& x3); - static ScalarType fwd_centered(const VectorType& x1, const VectorType& x2, const VectorType& x3); + static ScalarType fwd_centered(bool mixed_diff, const VectorType& x1, const VectorType& x2, const VectorType& x3); static void bwd(const Interval& y, IntervalVector& x1, IntervalVector& x2, IntervalVector& x3); }; @@ -120,45 +121,54 @@ namespace codac2 inline ScalarType DetOp::fwd_natural(const MatrixType& x) { return { - fwd(x.a), + fwd(x.get_a()), x.def_domain }; } - inline ScalarType DetOp::fwd_centered(const MatrixType& x) + inline ScalarType DetOp::fwd_centered(bool mixed_diff, const MatrixType& x) { if(centered_form_not_available_for_args(x)) return fwd_natural(x); - if(x.a.rows()==1) { + if(x.get_a().rows()==1) { return { fwd(x.m), - fwd(x.a), + vector_map(x.a, + static_cast(&fwd)), x.da, // keep matrix type for diff x.def_domain }; } IntervalMatrix d(1, x.da.cols()); - if (x.a.rows()==2) { /* otherwise, will fail afterwards */ + if (x.get_a().rows()==2) { /* otherwise, will fail afterwards */ for (Index i=0; i < d.cols() ; i++) { - d(0,i) = x.da(0,i)*x.a(1,1) + x.da(3,i)*x.a(0,0) - - x.da(1,i)*x.a(0,1) - x.da(2,i)*x.a(1,0); + const IntervalMatrix& xa = x.get_a(mixed_diff,i); + d(0,i) = x.da(0,i)*xa(1,1) + x.da(3,i)*xa(0,0) + - x.da(1,i)*xa(0,1) - x.da(2,i)*xa(1,0); /* note: ColMajor is assumed here */ } return { fwd(x.m), - fwd(x.a), + vector_map(x.a, + static_cast(&fwd)), d, x.def_domain }; } else { - auto ret1 = cofactor_matrix_enclosure(x.a); - for (Index i=0; i < d.cols() ; i++) - d(0,i) = ret1.first.reshaped().dot(x.da.col(i)); + std::vector cofact(x.a.size()); + std::vector detval(x.a.size()); + for (size_t i=0;i().dot(x.da.col(i)); + } return { fwd(x.m), - ret1.second, + detval, d, x.def_domain }; @@ -201,7 +211,7 @@ namespace codac2 inline ScalarType DetOp::fwd_natural(const VectorType& x1, const VectorType& x2) { IntervalMatrix a(2,2); - a.col(0) = x1.a; a.col(1) = x2.a; + a.col(0) = x1.get_a(); a.col(1) = x2.get_a(); return { fwd(a), @@ -209,26 +219,25 @@ namespace codac2 }; } - inline ScalarType DetOp::fwd_centered(const VectorType& x1, const VectorType& x2) + inline ScalarType DetOp::fwd_centered(bool mixed_diff, const VectorType& x1, const VectorType& x2) { if(centered_form_not_available_for_args(x1,x2)) return fwd_natural(x1,x2); - IntervalMatrix m(2,2); - m.col(0) = x1.m; m.col(1) = x2.m; - IntervalMatrix a(2,2); - a.col(0) = x1.a; a.col(1) = x2.a; - assert(x1.da.cols() == x2.da.cols()); IntervalMatrix d(1, x1.da.cols()); for (Index i=0; i < d.cols() ; i++) { - d(0,i) = x1.da(0,i)*x2.a[1] + x1.a[0]*x2.da(1,i) - - x1.da(1,i)*x2.a[0] - x1.a[1]*x2.da(0,i); + const IntervalVector& x1a = mixed_diff ? x1.a[i] : x1.get_a(); + const IntervalVector& x2a = mixed_diff ? x2.a[i] : x2.get_a(); + d(0,i) = x1.da(0,i)*x2a[1] + x1a[0]*x2.da(1,i) + - x1.da(1,i)*x2a[0] - x1a[1]*x2.da(0,i); } return { - fwd(m), - fwd(a), + fwd(x1.m,x2.m), + vector_map(x1.a,x2.a, + static_cast + (&fwd)), d, x1.def_domain && x2.def_domain }; @@ -255,7 +264,7 @@ namespace codac2 inline ScalarType DetOp::fwd_natural(const VectorType& x1, const VectorType& x2, const VectorType& x3) { IntervalMatrix a(3,3); - a.col(0) = x1.a; a.col(1) = x2.a; a.col(2) = x3.a; + a.col(0) = x1.get_a(); a.col(1) = x2.get_a(); a.col(2) = x3.get_a(); return { fwd(a), @@ -263,28 +272,36 @@ namespace codac2 }; } - inline ScalarType DetOp::fwd_centered(const VectorType& x1, const VectorType& x2, const VectorType& x3) + inline ScalarType DetOp::fwd_centered(bool mixed_diff, const VectorType& x1, const VectorType& x2, const VectorType& x3) { if(centered_form_not_available_for_args(x1,x2,x3)) return fwd_natural(x1,x2,x3); IntervalMatrix m(3,3); m.col(0) = x1.m; m.col(1) = x2.m; m.col(2) = x3.m; - IntervalMatrix a(3,3); - a.col(0) = x1.a; a.col(1) = x2.a; a.col(2) = x3.a; assert(x1.da.cols() == x2.da.cols()); assert(x1.da.cols() == x3.da.cols()); IntervalMatrix d(1, x1.da.cols()); - auto ret1 = cofactor_matrix_enclosure(a); + Index size_a = (mixed_diff ? d.cols() : 1); + std::vector cofact(size_a); + std::vector detval(size_a); + for (size_t i=0;i<(size_t) size_a;i++) { + IntervalMatrix a(3,3); + a.col(0) = x1.get_a(mixed_diff,i); + a.col(1) = x2.get_a(mixed_diff,i); + a.col(2) = x3.get_a(mixed_diff,i); + std::tie(cofact[i],detval[i])=cofactor_matrix_enclosure(a); + } for (Index i=0; i < d.cols() ; i++) { - d(0,i) = ret1.first.col(0).dot(x1.da.col(i)); - d(0,i) += ret1.first.col(1).dot(x2.da.col(i)); - d(0,i) += ret1.first.col(2).dot(x3.da.col(i)); + const IntervalMatrix& cof = mixed_diff ? cofact[i] : cofact.back(); + d(0,i) = cof.col(0).dot(x1.da.col(i)); + d(0,i) += cof.col(1).dot(x2.da.col(i)); + d(0,i) += cof.col(2).dot(x3.da.col(i)); } return { - fwd(m), - ret1.second, + fwd(x1.m,x2.m,x3.m), + detval, d, x1.def_domain && x2.def_domain && x3.def_domain }; diff --git a/src/core/operators/codac2_exp.h b/src/core/operators/codac2_exp.h index 19fd4cb8d..a269156ff 100644 --- a/src/core/operators/codac2_exp.h +++ b/src/core/operators/codac2_exp.h @@ -31,7 +31,7 @@ namespace codac2 static Interval fwd(const Interval& x1); static ScalarType fwd_natural(const ScalarType& x1); - static ScalarType fwd_centered(const ScalarType& x1); + static ScalarType fwd_centered(bool mixed_diff, const ScalarType& x1); static void bwd(const Interval& y, Interval& x1); }; @@ -54,23 +54,25 @@ namespace codac2 inline ScalarType ExpOp::fwd_natural(const ScalarType& x1) { return { - fwd(x1.a), + fwd(x1.get_a()), x1.def_domain }; } - inline ScalarType ExpOp::fwd_centered(const ScalarType& x1) + inline ScalarType ExpOp::fwd_centered(bool mixed_diff, const ScalarType& x1) { if(centered_form_not_available_for_args(x1)) return fwd_natural(x1); IntervalMatrix d(1,x1.da.size()); - for(Index i = 0 ; i < d.size() ; i++) - d(0,i) = x1.da(0,i)*exp(x1.a); + for(Index i = 0 ; i < d.size() ; i++) { + const Interval& x1a = x1.get_a(mixed_diff,i); + d(0,i) = x1.da(0,i)*exp(x1a); + } return { fwd(x1.m), - fwd(x1.a), + vector_map(x1.a,&fwd), d, x1.def_domain }; diff --git a/src/core/operators/codac2_extend.h b/src/core/operators/codac2_extend.h index faef2f7a4..e2cc77330 100644 --- a/src/core/operators/codac2_extend.h +++ b/src/core/operators/codac2_extend.h @@ -44,17 +44,19 @@ namespace codac2 static inline T fwd_natural(const T& x1, const T& x2) { return { - (x1.def_domain ? x1.a : fwd(x1.a,x2.a)), + (x1.def_domain ? x1.get_a() : fwd(x1.get_a(),x2.get_a())), x1.def_domain || x2.def_domain }; } template - static inline T fwd_centered(const T& x1, const T& x2) + static inline T fwd_centered([[maybe_unused]] bool mixed_diff, const T& x1, const T& x2) { return { (x1.def_domain ? x1.m : fwd(x1.m,x2.m)), - (x1.def_domain ? x1.a : fwd(x1.a,x2.a)), + (x1.def_domain ? x1.a : + vector_map(x1.a,x2.a, + static_cast(&fwd))), (x1.def_domain ? x1.da : (x1.da | x2.da)), x1.def_domain || x2.def_domain }; @@ -83,4 +85,4 @@ namespace codac2 return { std::make_shared>( (AnalyticExprWrapper)x1,(AnalyticExprWrapper)x2) }; } -} \ No newline at end of file +} diff --git a/src/core/operators/codac2_flatten.h b/src/core/operators/codac2_flatten.h index 3c0314643..b70ffba11 100644 --- a/src/core/operators/codac2_flatten.h +++ b/src/core/operators/codac2_flatten.h @@ -33,7 +33,7 @@ namespace codac2 static IntervalVector fwd(const IntervalMatrix& x1); static VectorType fwd_natural(const MatrixType& x1); - static VectorType fwd_centered(const MatrixType& x1); + static VectorType fwd_centered(bool mixed_diff, const MatrixType& x1); static void bwd(const IntervalVector& y, IntervalMatrix& x1); }; @@ -56,19 +56,19 @@ namespace codac2 inline VectorType FlattenOp::fwd_natural(const MatrixType& x1) { return { - fwd(x1.a), + fwd(x1.get_a()), x1.def_domain }; } - inline VectorType FlattenOp::fwd_centered(const MatrixType& x1) + inline VectorType FlattenOp::fwd_centered([[maybe_unused]] bool mixed_diff,const MatrixType& x1) { if(centered_form_not_available_for_args(x1)) return fwd_natural(x1); return { fwd(x1.m), - fwd(x1.a), + vector_map(x1.a,&fwd), x1.da, x1.def_domain }; diff --git a/src/core/operators/codac2_floor.h b/src/core/operators/codac2_floor.h index 12d99679e..d8f4cc406 100644 --- a/src/core/operators/codac2_floor.h +++ b/src/core/operators/codac2_floor.h @@ -31,7 +31,7 @@ namespace codac2 static Interval fwd(const Interval& x1); static ScalarType fwd_natural(const ScalarType& x1); - static ScalarType fwd_centered(const ScalarType& x1); + static ScalarType fwd_centered(bool mixed_type,const ScalarType& x1); static void bwd(const Interval& y, Interval& x1); }; @@ -54,19 +54,19 @@ namespace codac2 inline ScalarType FloorOp::fwd_natural(const ScalarType& x1) { return { - fwd(x1.a), + fwd(x1.get_a()), x1.def_domain }; } - inline ScalarType FloorOp::fwd_centered(const ScalarType& x1) + inline ScalarType FloorOp::fwd_centered([[maybe_unused]] bool mixed_diff, const ScalarType& x1) { if(centered_form_not_available_for_args(x1)) return fwd_natural(x1); return { fwd(x1.m), - fwd(x1.a), + vector_map(x1.a,&fwd), IntervalMatrix(0,0), // not supported yet for auto diff x1.def_domain }; diff --git a/src/core/operators/codac2_log.h b/src/core/operators/codac2_log.h index de82d4f52..fc44ef369 100644 --- a/src/core/operators/codac2_log.h +++ b/src/core/operators/codac2_log.h @@ -31,7 +31,7 @@ namespace codac2 static Interval fwd(const Interval& x1); static ScalarType fwd_natural(const ScalarType& x1); - static ScalarType fwd_centered(const ScalarType& x1); + static ScalarType fwd_centered(bool mixed_diff, const ScalarType& x1); static void bwd(const Interval& y, Interval& x1); }; @@ -57,28 +57,30 @@ namespace codac2 return fwd_natural(x1); return { - fwd(x1.a), - x1.a.is_subset({0,oo}) // def domain of log - && !x1.a.contains(0.) // def domain of the derivative of log + fwd(x1.get_a()), + x1.get_a().is_subset({0,oo}) // def domain of log + && !x1.get_a().contains(0.) // def domain of the derivative of log && x1.def_domain }; } - inline ScalarType LogOp::fwd_centered(const ScalarType& x1) + inline ScalarType LogOp::fwd_centered(bool mixed_diff,const ScalarType& x1) { if(centered_form_not_available_for_args(x1)) return fwd_natural(x1); IntervalMatrix d(1,x1.da.size()); - for(Index i = 0 ; i < d.size() ; i++) - d(0,i) = x1.da(0,i)/x1.a; + for(Index i = 0 ; i < d.size() ; i++) { + const Interval& x1a = x1.get_a(mixed_diff,i); + d(0,i) = x1.da(0,i)/x1a; + } return { fwd(x1.m), - fwd(x1.a), + vector_map(x1.a,&fwd), d, - x1.a.is_subset({0,oo}) // def domain of log - && !x1.a.contains(0.) // def domain of the derivative of log + x1.get_a().is_subset({0,oo}) // def domain of log + && !x1.get_a().contains(0.) // def domain of the derivative of log && x1.def_domain }; } diff --git a/src/core/operators/codac2_mat.h b/src/core/operators/codac2_mat.h index c61e1cbb4..42ca95957 100644 --- a/src/core/operators/codac2_mat.h +++ b/src/core/operators/codac2_mat.h @@ -59,7 +59,7 @@ namespace codac2 static inline MatrixType fwd_natural(const X&... x) { return { - MatrixOp::fwd(x.a...), + MatrixOp::fwd(x.get_a()...), (x.def_domain && ...) }; } @@ -74,21 +74,28 @@ namespace codac2 template requires (std::is_base_of_v && (std::is_base_of_v && ...)) - static inline MatrixType fwd_centered(const X1& x1, const X&... x) + static inline MatrixType fwd_centered(bool mixed_diff, const X1& x1, const X&... x) { if (centered_form_not_available_for_args(x1,x...)) return fwd_natural(x1,x...); - IntervalMatrix d(x1.a.size()*(1+sizeof...(X)),x1.da.cols()); + IntervalMatrix d(x1.get_a().size()*(1+sizeof...(X)),x1.da.cols()); Index l=0; d.topRows(x1.da.rows()) = x1.da; l += x1.da.rows(); ( MatrixOp::fill_diff_matrix(d,x.da,l) , ...); assert (l==d.rows()); + size_t mns = (mixed_diff ? std::max({ x1.a.size(), x.a.size()... }) : 1); + std::vector ret_a(mns); + for (size_t i=0;iinit_value(v, natural_fwd(v, total_input_size)); @@ -166,7 +173,7 @@ namespace codac2 bool centered_available = true; for(const auto& x : _xs) { - vals.push_back(x->fwd_eval(v, total_input_size, false)); + vals.push_back(x->fwd_eval(v, total_input_size, false, mixed_diff)); centered_available &= (vals.back().da.size() != 0); } @@ -174,9 +181,12 @@ namespace codac2 return this->init_value(v, natural_fwd(v, total_input_size)); const Index cols = static_cast(_xs.size()); - const Index rows = vals.empty() ? 0 : vals.front().a.size(); + const Index rows = vals.empty() ? 0 : vals.front().get_a().size(); const Index input_cols = vals.empty() ? total_input_size : vals.front().da.cols(); - IntervalMatrix m(rows, cols), a(rows, cols), da(rows*cols, input_cols); + IntervalMatrix m(rows, cols), da(rows*cols, input_cols); + std::vector + a(vals.empty() ? 0 : (mixed_diff ? input_cols : 1), + IntervalMatrix(rows,cols)); bool def_domain = true; Index l = 0; @@ -184,7 +194,9 @@ namespace codac2 { const auto& xj = vals[static_cast(j)]; m.col(j) = xj.m; - a.col(j) = xj.a; + for (size_t i=0;i(j)]->fwd_eval(v, total_input_size, true); - a.col(j) = xj.a; + auto xj = _xs[static_cast(j)]->fwd_eval(v, total_input_size, true, false); + a.col(j) = xj.get_a(); def_domain &= xj.def_domain; } diff --git a/src/core/operators/codac2_max.h b/src/core/operators/codac2_max.h index 514974f9f..7239497d5 100644 --- a/src/core/operators/codac2_max.h +++ b/src/core/operators/codac2_max.h @@ -31,7 +31,7 @@ namespace codac2 static Interval fwd(const Interval& x1, const Interval& x2); static ScalarType fwd_natural(const ScalarType& x1, const ScalarType& x2); - static ScalarType fwd_centered(const ScalarType& x1, const ScalarType& x2); + static ScalarType fwd_centered(bool mixed_diff, const ScalarType& x1, const ScalarType& x2); static void bwd(const Interval& y, Interval& x1, Interval& x2); }; @@ -47,12 +47,12 @@ namespace codac2 inline ScalarType MaxOp::fwd_natural(const ScalarType& x1, const ScalarType& x2) { return { - fwd(x1.a, x2.a), + fwd(x1.get_a(), x2.get_a()), x1.def_domain && x2.def_domain }; } - inline ScalarType MaxOp::fwd_centered(const ScalarType& x1, const ScalarType& x2) + inline ScalarType MaxOp::fwd_centered(bool mixed_diff, const ScalarType& x1, const ScalarType& x2) { if(centered_form_not_available_for_args(x1,x2)) return fwd_natural(x1,x2); @@ -61,15 +61,19 @@ namespace codac2 assert(x1.da.rows() == x2.da.rows() && x1.da.cols() == x2.da.cols()); IntervalMatrix d(1,x1.da.cols()); - for(Index i = 0 ; i < d.size() ; i++) - d(0,i) = chi(x1.a-x2.a, x2.da(0,i), x1.da(0,i)); + for(Index i = 0 ; i < d.size() ; i++) { + const Interval& x1a = x1.get_a(mixed_diff,i); + const Interval& x2a = x2.get_a(mixed_diff,i); + d(0,i) = chi(x1a-x2a, x2.da(0,i), x1.da(0,i)); + } return { fwd(x1.m, x2.m), - fwd(x1.a, x2.a), + vector_map(x1.a,x2.a,&fwd), d, x1.def_domain && x2.def_domain - && !x1.a.intersects(x2.a) // def domain of the derivative of max + && !x1.get_a().intersects(x2.get_a()) + // def domain of the derivative of max }; } diff --git a/src/core/operators/codac2_min.h b/src/core/operators/codac2_min.h index a8fee3df0..cf8cd0305 100644 --- a/src/core/operators/codac2_min.h +++ b/src/core/operators/codac2_min.h @@ -31,7 +31,7 @@ namespace codac2 static Interval fwd(const Interval& x1, const Interval& x2); static ScalarType fwd_natural(const ScalarType& x1, const ScalarType& x2); - static ScalarType fwd_centered(const ScalarType& x1, const ScalarType& x2); + static ScalarType fwd_centered(bool mixed_diff, const ScalarType& x1, const ScalarType& x2); static void bwd(const Interval& y, Interval& x1, Interval& x2); }; @@ -54,12 +54,12 @@ namespace codac2 inline ScalarType MinOp::fwd_natural(const ScalarType& x1, const ScalarType& x2) { return { - fwd(x1.a, x2.a), + fwd(x1.get_a(), x2.get_a()), x1.def_domain && x2.def_domain }; } - inline ScalarType MinOp::fwd_centered(const ScalarType& x1, const ScalarType& x2) + inline ScalarType MinOp::fwd_centered(bool mixed_diff, const ScalarType& x1, const ScalarType& x2) { if(centered_form_not_available_for_args(x1,x2)) return fwd_natural(x1,x2); @@ -68,15 +68,18 @@ namespace codac2 assert(x1.da.rows() == x2.da.rows() && x1.da.cols() == x2.da.cols()); IntervalMatrix d(1,x1.da.cols()); - for(Index i = 0 ; i < d.size() ; i++) - d(0,i) = chi(x1.a-x2.a, x1.da(0,i), x2.da(0,i)); + for(Index i = 0 ; i < d.size() ; i++) { + const Interval& x1a = x1.get_a(mixed_diff,i); + const Interval& x2a = x2.get_a(mixed_diff,i); + d(0,i) = chi(x1a-x2a, x1.da(0,i), x2.da(0,i)); + } return { fwd(x1.m, x2.m), - fwd(x1.a, x2.a), + vector_map(x1.a,x2.a,&fwd), d, x1.def_domain && x2.def_domain - && (x1.a != x2.a) // def domain of the derivative of min + && (x1.get_a() != x2.get_a()) // def domain of the derivative of min }; } diff --git a/src/core/operators/codac2_mod.h b/src/core/operators/codac2_mod.h index 02aefb661..1fb47b2c8 100644 --- a/src/core/operators/codac2_mod.h +++ b/src/core/operators/codac2_mod.h @@ -39,7 +39,7 @@ namespace codac2 static Interval fwd(const Interval& x1, const Interval& p); static ScalarType fwd_natural(const ScalarType& x1, const ScalarType& p); - static ScalarType fwd_centered(const ScalarType& x1, const ScalarType& p); + static ScalarType fwd_centered(bool mixed_diff, const ScalarType& x1, const ScalarType& p); static void bwd(const Interval& y, Interval& x1, Interval& p); static void fwd_bwd(Interval& y, Interval& x1, double p); @@ -99,19 +99,19 @@ namespace codac2 inline ScalarType ModOp::fwd_natural(const ScalarType& x1, const ScalarType& p) { return { - fwd(x1.a, p.a), + fwd(x1.get_a(), p.get_a()), x1.def_domain && p.def_domain }; } - inline ScalarType ModOp::fwd_centered(const ScalarType& x1, const ScalarType& p) + inline ScalarType ModOp::fwd_centered([[maybe_unused]] bool mixed_diff, const ScalarType& x1, const ScalarType& p) { if(centered_form_not_available_for_args(x1,p)) return fwd_natural(x1,p); return { fwd(x1.m,p.m), - fwd(x1.a,p.a), + vector_map(x1.a,p.a,&fwd), IntervalMatrix(0,0), // not supported yet for auto diff x1.def_domain && p.def_domain }; diff --git a/src/core/operators/codac2_pow.h b/src/core/operators/codac2_pow.h index 01c44aa0e..8ae8afa91 100644 --- a/src/core/operators/codac2_pow.h +++ b/src/core/operators/codac2_pow.h @@ -32,7 +32,7 @@ namespace codac2 static Interval fwd(const Interval& x1, const Interval& x2); static ScalarType fwd_natural(const ScalarType& x1, const ScalarType& x2); - static ScalarType fwd_centered(const ScalarType& x1, const ScalarType& x2); + static ScalarType fwd_centered(bool mixed_diff, const ScalarType& x1, const ScalarType& x2); static void bwd(const Interval& y, Interval& x1, int x2); static void bwd(const Interval& y, Interval& x1, Interval& x2); }; @@ -71,49 +71,54 @@ namespace codac2 inline ScalarType PowOp::fwd_natural(const ScalarType& x1, const ScalarType& x2) { - bool x2isint = x2.a.is_integer(); - bool x2positive = x2.a.lb()>=0.0; + bool x2isint = x2.get_a().is_integer(); + bool x2positive = x2.get_a().lb()>=0.0; return { - fwd(x1.a, x2.a), + fwd(x1.get_a(), x2.get_a()), x1.def_domain && x2.def_domain - && (x2isint || x1.a.lb()>=0.0) - && (x2positive || !x1.a.contains(0.0)) + && (x2isint || x1.get_a().lb()>=0.0) + && (x2positive || !x1.get_a().contains(0.0)) }; } - inline ScalarType PowOp::fwd_centered(const ScalarType& x1, const ScalarType& x2) + inline ScalarType PowOp::fwd_centered(bool mixed_diff, const ScalarType& x1, const ScalarType& x2) { if(centered_form_not_available_for_args(x1,x2)) return fwd_natural(x1,x2); - bool x2isint = x2.a.is_integer(); - bool x2positive = x2.a.lb()>=0.0; + bool x2isint = x2.get_a().is_integer(); + bool x2positive = x2.get_a().lb()>=0.0; IntervalMatrix d(1,x1.da.size()); - if (x2.a.is_degenerated()) + if (x2.get_a().is_degenerated()) { // To avoid calling log(x1.a), as it would return emptyset if x1.a<=0 - for(Index i = 0 ; i < d.size() ; i++) - d(0,i) = x2.a*x1.da(0,i)*pow(x1.a,x2.a-1.); + for(Index i = 0 ; i < d.size() ; i++) { + const Interval& x1a = x1.get_a(mixed_diff,i); + const Interval& x2a = x2.get_a(mixed_diff,i); + d(0,i) = x2a*x1.da(0,i)*pow(x1a,x2a-1.); + } } else { for(Index i = 0 ; i < d.size() ; i++) { - d(0,i) = x2.a*x1.da(0,i)*pow(x1.a,x2.a-1.) - + x2.da(0,i)*log(x1.a)*pow(x1.a,x2.a); + const Interval& x1a = x1.get_a(mixed_diff,i); + const Interval& x2a = x2.get_a(mixed_diff,i); + d(0,i) = x2a*x1.da(0,i)*pow(x1a,x2a-1.) + + x2.da(0,i)*log(x1a)*pow(x1a,x2a); // Not good when x1 close to 0 } } return { fwd(x1.m, x2.m), - fwd(x1.a, x2.a), + vector_map(x1.a,x2.a,&fwd), d, x1.def_domain && x2.def_domain - && (x2isint || x1.a.lb()>=0.0) - && (x2positive || !x1.a.contains(0.0)) + && (x2isint || x1.get_a().lb()>=0.0) + && (x2positive || !x1.get_a().contains(0.0)) }; } diff --git a/src/core/operators/codac2_sign.h b/src/core/operators/codac2_sign.h index 2d048bbe8..d3bdab27b 100644 --- a/src/core/operators/codac2_sign.h +++ b/src/core/operators/codac2_sign.h @@ -31,7 +31,7 @@ namespace codac2 static Interval fwd(const Interval& x1); static ScalarType fwd_natural(const ScalarType& x1); - static ScalarType fwd_centered(const ScalarType& x1); + static ScalarType fwd_centered(bool mixed_diff,const ScalarType& x1); static void bwd(const Interval& y, Interval& x1); }; @@ -54,21 +54,21 @@ namespace codac2 inline ScalarType SignOp::fwd_natural(const ScalarType& x1) { return { - fwd(x1.a), - x1.def_domain && !x1.a.contains(0.) // def domain of the derivative of sign + fwd(x1.get_a()), + x1.def_domain && !x1.get_a().contains(0.) // def domain of the derivative of sign }; } - inline ScalarType SignOp::fwd_centered(const ScalarType& x1) + inline ScalarType SignOp::fwd_centered([[maybe_unused]] bool mixed_diff, const ScalarType& x1) { if(centered_form_not_available_for_args(x1)) return fwd_natural(x1); return { fwd(x1.m), - fwd(x1.a), + vector_map(x1.a,&fwd), IntervalMatrix::zero(1,x1.da.size()), - x1.def_domain && !x1.a.contains(0.) // def domain of the derivative of sign + x1.def_domain && !x1.get_a().contains(0.) // def domain of the derivative of sign }; } diff --git a/src/core/operators/codac2_sin.h b/src/core/operators/codac2_sin.h index 749cbbd95..e9d9f40b9 100644 --- a/src/core/operators/codac2_sin.h +++ b/src/core/operators/codac2_sin.h @@ -31,7 +31,7 @@ namespace codac2 static Interval fwd(const Interval& x1); static ScalarType fwd_natural(const ScalarType& x1); - static ScalarType fwd_centered(const ScalarType& x1); + static ScalarType fwd_centered(bool mixed_diff, const ScalarType& x1); static void bwd(const Interval& y, Interval& x1); }; @@ -54,23 +54,25 @@ namespace codac2 inline ScalarType SinOp::fwd_natural(const ScalarType& x1) { return { - fwd(x1.a), + fwd(x1.get_a()), x1.def_domain }; } - inline ScalarType SinOp::fwd_centered(const ScalarType& x1) + inline ScalarType SinOp::fwd_centered(bool mixed_diff, const ScalarType& x1) { if(centered_form_not_available_for_args(x1)) return fwd_natural(x1); IntervalMatrix d(1,x1.da.size()); - for(Index i = 0 ; i < d.size() ; i++) - d(0,i) = cos(x1.a)*x1.da(0,i); + for(Index i = 0 ; i < d.size() ; i++) { + const Interval& x1a = x1.get_a(mixed_diff,i); + d(0,i) = cos(x1a)*x1.da(0,i); + } return { fwd(x1.m), - fwd(x1.a), + vector_map(x1.a,&fwd), d, x1.def_domain }; diff --git a/src/core/operators/codac2_sinh.h b/src/core/operators/codac2_sinh.h index 789f28f42..90a32e3cb 100644 --- a/src/core/operators/codac2_sinh.h +++ b/src/core/operators/codac2_sinh.h @@ -31,7 +31,7 @@ namespace codac2 static Interval fwd(const Interval& x1); static ScalarType fwd_natural(const ScalarType& x1); - static ScalarType fwd_centered(const ScalarType& x1); + static ScalarType fwd_centered(bool mixed_diff, const ScalarType& x1); static void bwd(const Interval& y, Interval& x1); }; @@ -54,23 +54,25 @@ namespace codac2 inline ScalarType SinhOp::fwd_natural(const ScalarType& x1) { return { - fwd(x1.a), + fwd(x1.get_a()), x1.def_domain }; } - inline ScalarType SinhOp::fwd_centered(const ScalarType& x1) + inline ScalarType SinhOp::fwd_centered(bool mixed_diff, const ScalarType& x1) { if(centered_form_not_available_for_args(x1)) return fwd_natural(x1); IntervalMatrix d(1,x1.da.size()); - for(Index i = 0 ; i < d.size() ; i++) - d(0,i) = cosh(x1.a)*x1.da(0,i); + for(Index i = 0 ; i < d.size() ; i++) { + const Interval& x1a = x1.get_a(mixed_diff,i); + d(0,i) = cosh(x1a)*x1.da(0,i); + } return { fwd(x1.m), - fwd(x1.a), + vector_map(x1.a,&fwd), d, x1.def_domain }; diff --git a/src/core/operators/codac2_sqr.h b/src/core/operators/codac2_sqr.h index 15fdf0770..95211df24 100644 --- a/src/core/operators/codac2_sqr.h +++ b/src/core/operators/codac2_sqr.h @@ -31,7 +31,7 @@ namespace codac2 static Interval fwd(const Interval& x1); static ScalarType fwd_natural(const ScalarType& x1); - static ScalarType fwd_centered(const ScalarType& x1); + static ScalarType fwd_centered(bool mixed_diff, const ScalarType& x1); static void bwd(const Interval& y, Interval& x1); }; @@ -54,12 +54,12 @@ namespace codac2 inline ScalarType SqrOp::fwd_natural(const ScalarType& x1) { return { - fwd(x1.a), + fwd(x1.get_a()), x1.def_domain }; } - inline ScalarType SqrOp::fwd_centered(const ScalarType& x1) + inline ScalarType SqrOp::fwd_centered(bool mixed_diff, const ScalarType& x1) { if(centered_form_not_available_for_args(x1)) return fwd_natural(x1); @@ -67,12 +67,14 @@ namespace codac2 assert(x1.da.rows() == 1); IntervalMatrix d(1,x1.da.cols()); - for(Index i = 0 ; i < d.size() ; i++) - d(0,i) = 2.*x1.a*x1.da(0,i); + for(Index i = 0 ; i < d.size() ; i++) { + const Interval& x1a = x1.get_a(mixed_diff,i); + d(0,i) = 2.*x1a*x1.da(0,i); + } return { fwd(x1.m), - fwd(x1.a), + vector_map(x1.a,&fwd), d, x1.def_domain }; diff --git a/src/core/operators/codac2_sqrt.h b/src/core/operators/codac2_sqrt.h index 8e79d7375..fcaf1a9cf 100644 --- a/src/core/operators/codac2_sqrt.h +++ b/src/core/operators/codac2_sqrt.h @@ -31,7 +31,7 @@ namespace codac2 static Interval fwd(const Interval& x1); static ScalarType fwd_natural(const ScalarType& x1); - static ScalarType fwd_centered(const ScalarType& x1); + static ScalarType fwd_centered(bool mixed_diff, const ScalarType& x1); static void bwd(const Interval& y, Interval& x1); }; @@ -54,28 +54,30 @@ namespace codac2 inline ScalarType SqrtOp::fwd_natural(const ScalarType& x1) { return { - fwd(x1.a), - x1.a.is_subset({0,oo}) // def domain of sqrt - && !x1.a.contains(0.) // def domain of the derivative of sqrt + fwd(x1.get_a()), + x1.get_a().is_subset({0,oo}) // def domain of sqrt + && !x1.get_a().contains(0.) // def domain of the derivative of sqrt && x1.def_domain }; } - inline ScalarType SqrtOp::fwd_centered(const ScalarType& x1) + inline ScalarType SqrtOp::fwd_centered(bool mixed_diff, const ScalarType& x1) { if(centered_form_not_available_for_args(x1)) return fwd_natural(x1); IntervalMatrix d(1,x1.da.size()); - for(Index i = 0 ; i < d.size() ; i++) - d(0,i) = x1.da(0,i)/(2.*sqrt(x1.a)); + for(Index i = 0 ; i < d.size() ; i++) { + const Interval& x1a = x1.get_a(mixed_diff,i); + d(0,i) = x1.da(0,i)/(2.*sqrt(x1a)); + } return { fwd(x1.m), - fwd(x1.a), + vector_map(x1.a,&fwd), d, - x1.a.is_subset({0,oo}) // def domain of sqrt - && !x1.a.contains(0.) // def domain of the derivative of sqrt + x1.get_a().is_subset({0,oo}) // def domain of sqrt + && !x1.get_a().contains(0.) // def domain of the derivative of sqrt && x1.def_domain }; } diff --git a/src/core/operators/codac2_subvector.h b/src/core/operators/codac2_subvector.h index c18eb7226..47cb2168e 100644 --- a/src/core/operators/codac2_subvector.h +++ b/src/core/operators/codac2_subvector.h @@ -35,7 +35,7 @@ namespace codac2 static IntervalVector fwd(const IntervalVector& x1, Index i, Index j); static VectorType fwd_natural(const VectorType& x1, Index i, Index j); - static VectorType fwd_centered(const VectorType& x1, Index i, Index j); + static VectorType fwd_centered(bool mixed_diff, const VectorType& x1, Index i, Index j); static void bwd(const IntervalVector& y, IntervalVector& x1, Index i, Index j); }; @@ -64,19 +64,19 @@ namespace codac2 return OperationExprBase>::replace_arg(old_arg_id, new_expr); } - VectorType fwd_eval(ValuesMap& v, Index total_input_size, bool natural_eval) const + VectorType fwd_eval(ValuesMap& v, Index total_input_size, bool natural_eval, bool mixed_diff) const { if(natural_eval) return AnalyticExpr::init_value( - v, SubvectorOp::fwd_natural(std::get<0>(this->_x)->fwd_eval(v, total_input_size, natural_eval), _i, _j)); + v, SubvectorOp::fwd_natural(std::get<0>(this->_x)->fwd_eval(v, total_input_size, natural_eval,false), _i, _j)); else return AnalyticExpr::init_value( - v, SubvectorOp::fwd_centered(std::get<0>(this->_x)->fwd_eval(v, total_input_size, natural_eval), _i, _j)); + v, SubvectorOp::fwd_centered(mixed_diff,std::get<0>(this->_x)->fwd_eval(v, total_input_size, natural_eval,mixed_diff), _i, _j)); } void bwd_eval(ValuesMap& v) const { - SubvectorOp::bwd(AnalyticExpr::value(v).a, std::get<0>(this->_x)->value(v).a, _i, _j); + SubvectorOp::bwd(AnalyticExpr::value(v).get_a(), std::get<0>(this->_x)->value(v).get_a(), _i, _j); std::get<0>(this->_x)->bwd_eval(v); } @@ -118,20 +118,22 @@ namespace codac2 { assert(i >= 0 && i < x1.a.rows() && j >= i && j < x1.a.rows()); return { - fwd(x1.a,i,j), + fwd(x1.get_a(),i,j), x1.def_domain }; } - inline VectorType SubvectorOp::fwd_centered(const VectorType& x1, Index i, Index j) + inline VectorType SubvectorOp::fwd_centered([[maybe_unused]] bool mixed_diff,const VectorType& x1, Index i, Index j) { if(centered_form_not_available_for_args(x1)) return fwd_natural(x1,i,j); - assert(i >= 0 && i < x1.a.rows() && j >= i && j < x1.a.rows()); + assert(i >= 0 && i < x1.get_a().rows() && + j >= i && j < x1.get_a().rows()); return { fwd(x1.m,i,j), - fwd(x1.a,i,j), + vector_map(x1.a, + [=](const IntervalVector& v1)->IntervalVector {return fwd(v1,i,j);}), x1.da.block(i,0,j-i+1,x1.da.cols()), x1.def_domain }; diff --git a/src/core/operators/codac2_tan.h b/src/core/operators/codac2_tan.h index 87795874a..968a1833c 100644 --- a/src/core/operators/codac2_tan.h +++ b/src/core/operators/codac2_tan.h @@ -31,7 +31,7 @@ namespace codac2 static Interval fwd(const Interval& x1); static ScalarType fwd_natural(const ScalarType& x1); - static ScalarType fwd_centered(const ScalarType& x1); + static ScalarType fwd_centered(bool mixed_diff, const ScalarType& x1); static void bwd(const Interval& y, Interval& x1); }; @@ -54,25 +54,27 @@ namespace codac2 inline ScalarType TanOp::fwd_natural(const ScalarType& x1) { return { - fwd(x1.a), - x1.def_domain && !cos(x1.a).contains(0.) // def domain of the derivative of tan + fwd(x1.get_a()), + x1.def_domain && !cos(x1.get_a()).contains(0.) // def domain of the derivative of tan }; } - inline ScalarType TanOp::fwd_centered(const ScalarType& x1) + inline ScalarType TanOp::fwd_centered(bool mixed_diff, const ScalarType& x1) { if(centered_form_not_available_for_args(x1)) return fwd_natural(x1); IntervalMatrix d(1,x1.da.size()); - for(Index i = 0 ; i < d.size() ; i++) - d(0,i) = x1.da(0,i)/sqr(cos(x1.a)); + for(Index i = 0 ; i < d.size() ; i++) { + const Interval& x1a = x1.get_a(mixed_diff,i); + d(0,i) = x1.da(0,i)/sqr(cos(x1a)); + } return { fwd(x1.m), - fwd(x1.a), + vector_map(x1.a,&fwd), d, - x1.def_domain && !cos(x1.a).contains(0.) // def domain of the derivative of tan + x1.def_domain && !cos(x1.get_a()).contains(0.) // def domain of the derivative of tan }; } diff --git a/src/core/operators/codac2_tanh.h b/src/core/operators/codac2_tanh.h index 35a977807..82afe8ff8 100644 --- a/src/core/operators/codac2_tanh.h +++ b/src/core/operators/codac2_tanh.h @@ -31,7 +31,7 @@ namespace codac2 static Interval fwd(const Interval& x1); static ScalarType fwd_natural(const ScalarType& x1); - static ScalarType fwd_centered(const ScalarType& x1); + static ScalarType fwd_centered(bool mixed_diff, const ScalarType& x1); static void bwd(const Interval& y, Interval& x1); }; @@ -54,23 +54,25 @@ namespace codac2 inline ScalarType TanhOp::fwd_natural(const ScalarType& x1) { return { - fwd(x1.a), + fwd(x1.get_a()), x1.def_domain }; } - inline ScalarType TanhOp::fwd_centered(const ScalarType& x1) + inline ScalarType TanhOp::fwd_centered(bool mixed_diff, const ScalarType& x1) { if(centered_form_not_available_for_args(x1)) return fwd_natural(x1); IntervalMatrix d(1,x1.da.size()); - for(Index i = 0 ; i < d.size() ; i++) - d(0,i) = x1.da(0,i)/sqr(cosh(x1.a)); + for(Index i = 0 ; i < d.size() ; i++) { + const Interval& x1a = x1.get_a(mixed_diff,i); + d(0,i) = x1.da(0,i)/sqr(cosh(x1a)); + } return { fwd(x1.m), - fwd(x1.a), + vector_map(x1.a,&fwd), d, x1.def_domain }; diff --git a/src/core/operators/codac2_transpose.h b/src/core/operators/codac2_transpose.h index e3b1ca130..3ce2cc870 100644 --- a/src/core/operators/codac2_transpose.h +++ b/src/core/operators/codac2_transpose.h @@ -33,7 +33,7 @@ namespace codac2 static IntervalMatrix fwd(const IntervalMatrix& x1); static MatrixType fwd_natural(const MatrixType& x1); - static MatrixType fwd_centered(const MatrixType& x1); + static MatrixType fwd_centered(bool mixed_diff, const MatrixType& x1); static void bwd(const IntervalMatrix& y, IntervalMatrix& x1); }; @@ -56,26 +56,26 @@ namespace codac2 inline MatrixType TransposeOp::fwd_natural(const MatrixType& x1) { return { - fwd(x1.a), + fwd(x1.get_a()), x1.def_domain }; } - inline MatrixType TransposeOp::fwd_centered(const MatrixType& x1) + inline MatrixType TransposeOp::fwd_centered([[maybe_unused]] bool mixed_diff, const MatrixType& x1) { if(centered_form_not_available_for_args(x1)) return fwd_natural(x1); IntervalMatrix d(x1.da.rows(),x1.da.cols()); - for (Index i=0; i(x1.a,&fwd), d, x1.def_domain }; diff --git a/src/core/operators/codac2_vec.h b/src/core/operators/codac2_vec.h index 8e862e5c2..6f72f3afc 100644 --- a/src/core/operators/codac2_vec.h +++ b/src/core/operators/codac2_vec.h @@ -50,14 +50,14 @@ namespace codac2 ((def_domain &= x.def_domain), ...); return { - fwd(x.a...), + fwd(x.get_a()...), def_domain }; } template requires (std::is_base_of_v && ...) - static inline VectorType fwd_centered(const X&... x) + static inline VectorType fwd_centered(bool mixed_diff, const X&... x) { if(centered_form_not_available_for_args(x...)) return fwd_natural(x...); @@ -69,9 +69,16 @@ namespace codac2 bool def_domain = true; ((def_domain &= x.def_domain), ...); + size_t sm = (mixed_diff ? std::max({ x.a.size()... }) : 1); + std::vector va(sm); + for(size_t i=0;iinit_value(v, natural_fwd(v, total_input_size)); @@ -149,7 +156,7 @@ namespace codac2 bool centered_available = true; for(const auto& x : _xs) { - vals.push_back(x->fwd_eval(v, total_input_size, false)); + vals.push_back(x->fwd_eval(v, total_input_size, false,mixed_diff)); centered_available &= (vals.back().da.size() != 0); } @@ -158,7 +165,9 @@ namespace codac2 const Index n = static_cast(_xs.size()); const Index input_cols = vals.empty() ? total_input_size : vals.front().da.cols(); - IntervalVector m(n), a(n); + IntervalVector m(n); + std::vector a(vals.empty() ? 0 : + (mixed_diff ? input_cols : 1)); IntervalMatrix da(n, input_cols); bool def_domain = true; @@ -166,11 +175,12 @@ namespace codac2 { const auto& xi = vals[static_cast(i)]; m[i] = xi.m; - a[i] = xi.a; + for (size_t j=0;jinit_value(v, VectorType(m, a, da, def_domain)); } @@ -232,8 +242,8 @@ namespace codac2 bool def_domain = true; for(Index i = 0 ; i < n ; ++i) { - auto xi = _xs[static_cast(i)]->fwd_eval(v, total_input_size, true); - a[i] = xi.a; + auto xi = _xs[static_cast(i)]->fwd_eval(v, total_input_size, true,false); + a[i] = xi.get_a(); def_domain &= xi.def_domain; } diff --git a/src/core/trajectory/codac2_Traj_operator.h b/src/core/trajectory/codac2_Traj_operator.h index 80e8fa625..5d2bb7817 100644 --- a/src/core/trajectory/codac2_Traj_operator.h +++ b/src/core/trajectory/codac2_Traj_operator.h @@ -45,7 +45,7 @@ namespace codac2 }; } - static typename TR::Type fwd_centered(const TR& x1, [[maybe_unused]] const SampledTraj& x1_deriv, const ScalarType& x2) + static typename TR::Type fwd_centered([[maybe_unused]] bool mixed_diff, const TR& x1, [[maybe_unused]] const SampledTraj& x1_deriv, const ScalarType& x2) { IntervalMatrix d(x1.size(),x2.da.cols()); @@ -66,9 +66,11 @@ namespace codac2 // d(i,j) = dx1[i]; // } + // Note: I don't even try to understand there [-- Damien] + return { fwd(x1,x2.m), - fwd(x1,x2.a), + fwd(x1,x2.get_a()), d, x2.def_domain && x1.tdomain().is_superset(x2.m) }; @@ -111,19 +113,19 @@ namespace codac2 return OperationExprBase>::replace_arg(old_arg_id, new_expr); } - T fwd_eval(ValuesMap& v, Index total_input_size, bool natural_eval) const + T fwd_eval(ValuesMap& v, Index total_input_size, bool natural_eval, bool mixed_diff) const { if(natural_eval) return AnalyticExpr::init_value( - v, TrajectoryOp::fwd_natural(_x1, std::get<0>(this->_x)->fwd_eval(v, total_input_size, natural_eval))); + v, TrajectoryOp::fwd_natural(_x1, std::get<0>(this->_x)->fwd_eval(v, total_input_size, natural_eval, mixed_diff))); else return AnalyticExpr::init_value( - v, TrajectoryOp::fwd_centered(_x1, _x1_deriv, std::get<0>(this->_x)->fwd_eval(v, total_input_size, natural_eval))); + v, TrajectoryOp::fwd_centered(mixed_diff, _x1, _x1_deriv, std::get<0>(this->_x)->fwd_eval(v, total_input_size, natural_eval, mixed_diff))); } void bwd_eval(ValuesMap& v) const { - TrajectoryOp::bwd(_x1, AnalyticExpr::value(v).a, std::get<0>(this->_x)->value(v).a); + TrajectoryOp::bwd(_x1, AnalyticExpr::value(v).get_a(), std::get<0>(this->_x)->value(v).get_a()); std::get<0>(this->_x)->bwd_eval(v); } diff --git a/tests/core/functions/analytic/codac2_tests_AnalyticFunction.cpp b/tests/core/functions/analytic/codac2_tests_AnalyticFunction.cpp index acdf9011e..e3b370298 100644 --- a/tests/core/functions/analytic/codac2_tests_AnalyticFunction.cpp +++ b/tests/core/functions/analytic/codac2_tests_AnalyticFunction.cpp @@ -23,10 +23,12 @@ auto create_f() TEST_CASE("AnalyticFunction") { - std::array modes { + std::array modes { EvalMode::NATURAL, EvalMode::CENTERED, - EvalMode::NATURAL | EvalMode::NATURAL + EvalMode::NATURAL | EvalMode::CENTERED, + EvalMode::MIXED, + EvalMode::MIXED | EvalMode::NATURAL }; for(const auto& m : modes) // try all available evaluation modes @@ -86,6 +88,7 @@ TEST_CASE("AnalyticFunction") { AnalyticFunction f({x}, vec(x,x)); CHECK(f.output_shape()==std::pair(2,1)); + std::cout << (long int)m << " " << f.eval(m,1) << "\n"; CHECK(Approx(f.eval(m, 1)) == IntervalVector({1,1})); } } @@ -211,19 +214,19 @@ TEST_CASE("AnalyticFunction") //.def("__pos__", {}(const VectorExpr& e1); CHECK(Approx(AnalyticFunction({v1}, +(v1+v1)).eval(m, IntervalVector({{0.},{-oo,5}}))) == - (m == EvalMode::CENTERED ? IntervalVector({{0.},{-oo,oo}}) : IntervalVector({{0.},{-oo,10}}))); + (m == EvalMode::CENTERED || m == EvalMode::MIXED ? IntervalVector({{0.},{-oo,oo}}) : IntervalVector({{0.},{-oo,10}}))); //.def(py::self + py::self); CHECK(Approx(AnalyticFunction({v1}, v1+v1).eval(m, IntervalVector({{0.},{-oo,5}}))) == - (m == EvalMode::CENTERED ? IntervalVector({{0.},{-oo,oo}}) : IntervalVector({{0.},{-oo,10}}))); + (m == EvalMode::CENTERED || m == EvalMode::MIXED ? IntervalVector({{0.},{-oo,oo}}) : IntervalVector({{0.},{-oo,10}}))); //.def("__radd__", {}(const VectorExpr& e1, const VectorVar& e2); CHECK(Approx(AnalyticFunction({v1}, (v1+v1)+v1).eval(m, IntervalVector({{0.},{-oo,5}}))) == - (m == EvalMode::CENTERED ? IntervalVector({{0.},{-oo,oo}}) : IntervalVector({{0.},{-oo,15}}))); + (m == EvalMode::CENTERED || m == EvalMode::MIXED ? IntervalVector({{0.},{-oo,oo}}) : IntervalVector({{0.},{-oo,15}}))); //.def("__radd__", {}(const VectorExpr& e1, const IntervalVector& e2); CHECK(Approx(AnalyticFunction({v1}, v1+(v1+v1)).eval(m, IntervalVector({{0.},{-oo,5}}))) == - (m == EvalMode::CENTERED ? IntervalVector({{0.},{-oo,oo}}) : IntervalVector({{0.},{-oo,15}}))); + (m == EvalMode::CENTERED || m == EvalMode::MIXED ? IntervalVector({{0.},{-oo,oo}}) : IntervalVector({{0.},{-oo,15}}))); //.def(- py::self); CHECK(Approx(AnalyticFunction({v1}, -(v1+v1)).eval(m, IntervalVector({{0.},{-oo,5}}))) == - (m == EvalMode::CENTERED ? IntervalVector({{0.},{-oo,oo}}) : IntervalVector({{0.},{-10,oo}}))); + (m == EvalMode::CENTERED || m == EvalMode::MIXED ? IntervalVector({{0.},{-oo,oo}}) : IntervalVector({{0.},{-10,oo}}))); //.def(py::self - py::self); CHECK(Approx(AnalyticFunction({v1,v2}, (v1-v2)).eval(m, IntervalVector({2,3}),Vector({1,5}))) == IntervalVector({1,-2})); //.def("__sub__", {}(const VectorExpr& e1, const VectorVar& e2); @@ -236,11 +239,27 @@ TEST_CASE("AnalyticFunction") CHECK(Approx(AnalyticFunction({v1,v2}, Vector({1,5})-(v1-v2)).eval(m, IntervalVector({2,3}),Vector({1,5}))) == IntervalVector({0,7})); //.def("__rmul__", {}(const VectorExpr& e1, const ScalarExpr& e2); CHECK(Approx(AnalyticFunction({v1,v2}, cos(v1[0])*(v2+v2)).eval(m, Vector({PI,-1}),Vector({2,3})),1e-9) == IntervalVector({-4,-6})); - //.def("__rmul__", {}(const VectorExpr& e1, const ScalarVar& e2); + //.def("__rmul__", {}(const VectorExpr& e1, const DomainVar& e2); CHECK(Approx(AnalyticFunction({x1}, x1*vec(3*x1,2*x1)).eval(m, 3),1e-9) == IntervalVector({27,18})); } } + { VectorVar x(2); + ScalarVar a,b; + AnalyticFunction f0({x}, x[1]); + std::cout << "Natural : " << f0.eval(EvalMode::NATURAL,IntervalVector({{0,1},{0,1}})) << "\n"; + std::cout << "Centered : " << f0.eval(EvalMode::CENTERED,IntervalVector({{0,1},{0,1}})) << "\n"; + std::cout << "Mixed : " << f0.eval(EvalMode::MIXED,IntervalVector({{0,1},{0,1}})) << "\n"; + AnalyticFunction f0p({a,b}, (a+b)*(a-b)); + std::cout << "Natural : " << f0p.eval(EvalMode::NATURAL,Interval({0,1}),Interval({0,1})) << "\n"; + std::cout << "Centered: " << f0p.eval(EvalMode::CENTERED,Interval({0,1}),Interval({0,1})) << "\n"; + std::cout << "Mixed : " << f0p.eval(EvalMode::MIXED,Interval({0,1}),Interval({0,1})) << "\n"; + AnalyticFunction f({x}, (x[0]+x[1])*(x[0]-x[1])); + std::cout << "Natural : " << f.eval(EvalMode::NATURAL,IntervalVector({{0,1},{0,1}})) << "\n"; + std::cout << "Centered : " << f.eval(EvalMode::CENTERED,IntervalVector({{0,1},{0,1}})) << "\n"; + std::cout << "Mixed : " << f.eval(EvalMode::MIXED,IntervalVector({{0,1},{0,1}})) << "\n"; + } + { VectorVar x(2); AnalyticFunction f({x}, x[0]*(x[0]+x[1])+sqr(x[1])); @@ -274,6 +293,7 @@ TEST_CASE("AnalyticFunction") AnalyticFunction f({x}, x-x); CHECK(f.eval(EvalMode::NATURAL,Interval(-1,1)) == Interval(-2,2)); CHECK(f.eval(EvalMode::CENTERED,Interval(-1,1)) == Interval(0)); + CHECK(f.eval(EvalMode::MIXED,Interval(-1,1)) == Interval(0)); CHECK(f.eval(Interval(-1,1)) == Interval(0)); }