Skip to content

Mixed differentiation #549

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
damien-masse wants to merge 2 commits into
base: codac2_dev
Choose a base branch
from
Open

Mixed differentiation #549

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
2c07462
mix?
damien-masse Oct 22, 2026
d3ce3e9
Merge branch 'codac-team:codac2_dev' into codac2_dev
damien-masse Oct 22, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
4 changes: 2 additions & 2 deletions examples/02_centered_form/main.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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);
}
}
18 changes: 9 additions & 9 deletions python/src/core/contractors/codac2_py_CtcInverse.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -41,26 +41,26 @@ void export_CtcInverse(py::module& m,
exported_ctc

.def(py::init(
[](const AnalyticFunction<T>& f, const D& y, bool with_centered_form)
[](const AnalyticFunction<T>& f, const D& y, bool with_centered_form, bool with_mixed_diff)
{
return std::make_unique<C>(unaryize_function(f), y, with_centered_form);
return std::make_unique<C>(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<T,VectorType>) // contractors are only associated with interval vectors
{
exported_ctc
.def(py::init(
[](const AnalyticFunction<T>& f, const CtcBase<IntervalVector>& c, bool with_centered_form)
[](const AnalyticFunction<T>& f, const CtcBase<IntervalVector>& c, bool with_centered_form, bool with_mixed_diff)
{
return std::make_unique<C>(unaryize_function(f), c.copy(), with_centered_form);
return std::make_unique<C>(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,
Expand Down Expand Up @@ -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);
}
}
}
20 changes: 11 additions & 9 deletions src/core/actions/codac2_OctaSym_operator.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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<IntervalVector,IntervalVector>(x1.a,
[=](const IntervalVector& v1)->IntervalVector {return fwd(s,v1);}),
da,
x1.def_domain
};
Expand Down Expand Up @@ -86,19 +88,19 @@ namespace codac2
return OperationExprBase<AnalyticExpr<VectorType>>::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<VectorType>::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<VectorType>::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<VectorType>::value(v).a, std::get<0>(this->_x)->value(v).a);
OctaSymOp::bwd(_s, AnalyticExpr<VectorType>::value(v).get_a(), std::get<0>(this->_x)->value(v).get_a());
std::get<0>(this->_x)->bwd_eval(v);
}

Expand Down Expand Up @@ -126,4 +128,4 @@ namespace codac2

const OctaSym _s;
};
}
}
49 changes: 25 additions & 24 deletions src/core/contractors/codac2_CtcInverse.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -32,15 +32,15 @@ namespace codac2

template<typename C>
requires IsCtcBaseOrPtr<C,Y>
CtcInverse(const AnalyticFunction<OutputType>& f, const C& ctc_y, bool with_centered_form = true, bool is_not_in = false)
: Ctc<CtcInverse<Y,X0,X...>,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<OutputType>& f, const C& ctc_y, bool with_centered_form = true, bool with_mixed_diff = false, bool is_not_in = false)
: Ctc<CtcInverse<Y,X0,X...>,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<OutputType>& f, const Y& y, bool with_centered_form = true, bool is_not_in = false)
: CtcInverse(f, CtcWrapper<Y,Y>(y), with_centered_form, is_not_in)
CtcInverse(const AnalyticFunction<OutputType>& f, const Y& y, bool with_centered_form = true, bool with_mixed_diff = false, bool is_not_in = false)
: CtcInverse(f, CtcWrapper<Y,Y>(y), with_centered_form, with_mixed_diff, is_not_in)
{ }

void contract(X0& x0, X&... x) const
Expand All @@ -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)
Expand All @@ -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
Expand All @@ -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<Y,IntervalMatrix>)
{
Expand Down Expand Up @@ -126,6 +126,7 @@ namespace codac2
const AnalyticFunction<OutputType> _f;
const Collection<CtcBase<Y>> _ctc_y;
bool _with_centered_form;
bool _with_mixed_diff;
bool _is_not_in = false;
};

Expand All @@ -134,63 +135,63 @@ namespace codac2

// ScalarType

CtcInverse(const AnalyticFunction<ScalarType>&, std::initializer_list<double>, bool = true, bool = false) ->
CtcInverse(const AnalyticFunction<ScalarType>&, std::initializer_list<double>, bool = true, bool = false, bool = false) ->
CtcInverse<Interval,IntervalVector>;

template<typename Y>
CtcInverse(const AnalyticFunction<ScalarType>&, std::initializer_list<Y>, bool = true, bool = false) ->
CtcInverse(const AnalyticFunction<ScalarType>&, std::initializer_list<Y>, bool = true, bool = false, bool = false) ->
CtcInverse<Interval,IntervalVector>;

CtcInverse(const AnalyticFunction<ScalarType>&, const Interval&, bool = true, bool = false) ->
CtcInverse(const AnalyticFunction<ScalarType>&, const Interval&, bool = true, bool = false, bool = false) ->
CtcInverse<Interval,IntervalVector>;

CtcInverse(const AnalyticFunction<ScalarType>&, double, bool = true, bool = false) ->
CtcInverse(const AnalyticFunction<ScalarType>&, double, bool = true, bool = false, bool = false) ->
CtcInverse<Interval,IntervalVector>;

template<typename C>
requires IsCtcBaseOrPtr<C,Interval>
CtcInverse(const AnalyticFunction<ScalarType>&, const C&, bool = true, bool = false) ->
CtcInverse(const AnalyticFunction<ScalarType>&, const C&, bool = true, bool = false, bool = false) ->
CtcInverse<Interval,IntervalVector>;

// VectorType

CtcInverse(const AnalyticFunction<VectorType>&, std::initializer_list<double>, bool = true, bool = false) ->
CtcInverse(const AnalyticFunction<VectorType>&, std::initializer_list<double>, bool = true, bool = false, bool = false) ->
CtcInverse<IntervalVector,IntervalVector>;

CtcInverse(const AnalyticFunction<VectorType>&, std::initializer_list<std::initializer_list<double>>, bool = true, bool = false) ->
CtcInverse(const AnalyticFunction<VectorType>&, std::initializer_list<std::initializer_list<double>>, bool = true, bool = false, bool = false) ->
CtcInverse<IntervalVector,IntervalVector>;

CtcInverse(const AnalyticFunction<VectorType>&, const Vector&, bool = true, bool = false) ->
CtcInverse(const AnalyticFunction<VectorType>&, const Vector&, bool = true, bool = false, bool = false) ->
CtcInverse<IntervalVector,IntervalVector>;

CtcInverse(const AnalyticFunction<VectorType>&, const IntervalVector&, bool = true, bool = false) ->
CtcInverse(const AnalyticFunction<VectorType>&, const IntervalVector&, bool = true, bool = false, bool = false) ->
CtcInverse<IntervalVector,IntervalVector>;

template<typename OtherDerived>
requires (OtherDerived::RowsAtCompileTime == -1 && OtherDerived::ColsAtCompileTime == 1)
CtcInverse(const AnalyticFunction<VectorType>&, const Eigen::MatrixBase<OtherDerived>&, bool = true, bool = false) ->
CtcInverse(const AnalyticFunction<VectorType>&, const Eigen::MatrixBase<OtherDerived>&, bool = true, bool = false, bool = false) ->
CtcInverse<IntervalVector,IntervalVector>;

template<typename C>
requires IsCtcBaseOrPtr<C,IntervalVector>
CtcInverse(const AnalyticFunction<VectorType>&, const C&, bool = true, bool = false) ->
CtcInverse(const AnalyticFunction<VectorType>&, const C&, bool = true, bool = false, bool = false) ->
CtcInverse<IntervalVector,IntervalVector>;

// MatrixType

CtcInverse(const AnalyticFunction<MatrixType>&, std::initializer_list<std::initializer_list<double>>, bool = true, bool = false) ->
CtcInverse(const AnalyticFunction<MatrixType>&, std::initializer_list<std::initializer_list<double>>, bool = true, bool = false, bool = false) ->
CtcInverse<IntervalMatrix,IntervalVector>;

CtcInverse(const AnalyticFunction<MatrixType>&, std::initializer_list<std::initializer_list<std::initializer_list<double>>>, bool = true, bool = false) ->
CtcInverse(const AnalyticFunction<MatrixType>&, std::initializer_list<std::initializer_list<std::initializer_list<double>>>, bool = true, bool = false, bool = false) ->
CtcInverse<IntervalMatrix,IntervalVector>;

template<typename OtherDerived>
requires (OtherDerived::RowsAtCompileTime == -1 && OtherDerived::ColsAtCompileTime == -1)
CtcInverse(const AnalyticFunction<MatrixType>&, const Eigen::MatrixBase<OtherDerived>&, bool = true, bool = false) ->
CtcInverse(const AnalyticFunction<MatrixType>&, const Eigen::MatrixBase<OtherDerived>&, bool = true, bool = false, bool = false) ->
CtcInverse<IntervalMatrix,IntervalVector>;

template<typename C>
requires IsCtcBaseOrPtr<C,IntervalMatrix>
CtcInverse(const AnalyticFunction<MatrixType>&, const C&, bool = true, bool = false) ->
CtcInverse(const AnalyticFunction<MatrixType>&, const C&, bool = true, bool = false, bool = false) ->
CtcInverse<IntervalMatrix,IntervalVector>;
}
}
12 changes: 6 additions & 6 deletions src/core/contractors/codac2_CtcInverseNotIn.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -23,7 +23,7 @@ namespace codac2
using Y = typename Wrapper<Y_>::Domain; // Y_ is a possible Eigen expr. type
using OutputType = typename ExprType<Y>::Type;

CtcInverseNotIn(const AnalyticFunction<OutputType>& f, const Y& y, bool with_centered_form = true)
CtcInverseNotIn(const AnalyticFunction<OutputType>& f, const Y& y, bool with_centered_form = true, bool with_mixed_diff = false)
: CtcUnion<X0,X...>(f.input_size())
{
const bool is_not_in = true;
Expand All @@ -32,17 +32,17 @@ namespace codac2
{
Y w = complem_y & f_codomain;
if(!w.is_empty() && !w.is_subset(y))
*this |= CtcInverse<Y,X0,X...>(f, complem_y, with_centered_form, is_not_in);
*this |= CtcInverse<Y,X0,X...>(f, complem_y, with_centered_form, with_mixed_diff, is_not_in);
}
}

template<typename C>
requires IsCtcBaseOrPtr<C,Y>
CtcInverseNotIn(const AnalyticFunction<OutputType>& f, const C& ctc_compl, bool with_centered_form = true)
CtcInverseNotIn(const AnalyticFunction<OutputType>& f, const C& ctc_compl, bool with_centered_form = true, bool with_mixed_diff = false)
: CtcUnion<X0,X...>(f.input_size())
{
const bool is_not_in = true;
*this |= CtcInverse<Y,X0,X...>(f, ctc_compl, with_centered_form, is_not_in);
*this |= CtcInverse<Y,X0,X...>(f, ctc_compl, with_centered_form, with_mixed_diff, is_not_in);
}

void contract(X0& x0, X&... x) const
Expand All @@ -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

Expand Down Expand Up @@ -115,4 +115,4 @@ namespace codac2
requires IsCtcBaseOrPtr<C,IntervalMatrix>
CtcInverseNotIn(const AnalyticFunction<MatrixType>&, const C&, bool = true, bool = false) ->
CtcInverseNotIn<IntervalMatrix,IntervalVector>;
}
}
8 changes: 4 additions & 4 deletions src/core/domains/tube/codac2_Tube_operator.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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)
};
Expand Down Expand Up @@ -67,15 +67,15 @@ namespace codac2
return OperationExprBase<AnalyticExpr<ScalarType>>::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<T>::init_value(
v, TubeOp<TU>::fwd(_x1, std::get<0>(this->_x)->fwd_eval(v, total_input_size, natural_eval)));
v, TubeOp<TU>::fwd(_x1, std::get<0>(this->_x)->fwd_eval(v, total_input_size, natural_eval,mixed_diff)));
}

void bwd_eval(ValuesMap& v) const
{
TubeOp<TU>::bwd(_x1, AnalyticExpr<T>::value(v).a, std::get<0>(this->_x)->value(v).a);
TubeOp<TU>::bwd(_x1, AnalyticExpr<T>::value(v).get_a(), std::get<0>(this->_x)->value(v).get_a());
std::get<0>(this->_x)->bwd_eval(v);
}

Expand Down
17 changes: 10 additions & 7 deletions src/core/functions/analytic/codac2_AnalyticExpr.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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<Index,Index> output_shape() const = 0;

Expand Down Expand Up @@ -86,18 +86,21 @@ namespace codac2
return OperationExprBase<AnalyticExpr<X>...>::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<Y>::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<Y>::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<Y>::init_value(v,
C::fwd_centered(true,x->fwd_eval(v, total_input_size, natural_eval,mixed_diff)...));
},
this->_x);
}
Expand All @@ -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)
Expand Down
Loading
Loading