codac2_CtcInverse_.h

 
#pragma once
 
#include <map>
#include "codac2_AnalyticFunction.h"
#include "codac2_Ctc.h"
#include "codac2_CtcWrapper.h"
#include "codac2_Collection.h"
#include "codac2_arith_mul.h"
#include "codac2_Wrapper.h"
 
namespace codac2
{
  class CtcNot;
  
  template<typename Y, typename... X>
    requires (sizeof...(X) > 0)
  class CtcInverse : public Ctc<CtcInverse<Y,X...>,X...>
  {
    public:
 
      using Ctc<CtcInverse<Y,X...>,X...>::contract;
 
      template<typename C>
        requires IsCtcBaseOrPtr<C,Y>
      CtcInverse(const AnalyticFunction<typename ExprType<Y>::Type>& f, const C& ctc_y, bool with_centered_form = true, bool is_not_in = false, double freq = 0.)
        : Ctc<CtcInverse<Y,X...>,X...>(f.args()[0]->size()), _f(f), _ctc_y(ctc_y), _with_centered_form(with_centered_form), _is_not_in(is_not_in), _freq(freq)
      {
        assert_release([&]() { return f.output_size() == size_of(ctc_y); }()
          && "CtcInverse: invalid dimension of image argument ('y' or 'ctc_y')");
      }
 
      CtcInverse(const AnalyticFunction<typename ExprType<Y>::Type>& f, const Y& y, bool with_centered_form = true, bool is_not_in = false, double freq = 0.)
        : CtcInverse(f, CtcWrapper<Y,Y>(y), with_centered_form, is_not_in, freq)
      { }
 
      void contract(X&... x) const
      {
        return contract_(*_ctc_y.front(), x...);
      }
 
      void contract_(const Y& y, X&... x) const
      {
        return contract_(CtcWrapper<Y,Y>(y), x...);
      }
 
      void contract_(const CtcBase<Y>& ctc_y, X&... x) const
      {
        //using X0b = std::tuple_element_t<0, std::tuple<X...>>;
        //X0b x0_before = std::get<0>(std::tie(x...));
//double md = std::get<0>(std::tie(x...)).max_diam();
bool using_nat = true;//_both || md >= _freq;
bool using_centered = true;//_both || _both_then_one || md < _freq;
 
        ValuesMap v;
        // Setting user values into a map before the tree evaluation
        _f.fill_from_args(v, x...);
 
        // Forward/backward algorithm:
 
          // [1/4] Forward evaluation
          _f.expr()->fwd_eval(v, _f.args().total_size(), !using_centered);
          auto& val_expr = _f.expr()->value(v);
 
          if(_is_not_in && !val_expr.def_domain)
            return; // <-- iota: if the input x is outside the definition 
          // domain of one of the involved expressions, or their combinations,
          // then the inner contraction is disabled.
          
          // [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);
 
          // [3/4 - optional]
          // The contraction can be significantly improved using a centered form
          // expression (enabled by default). This step must be processed before the
          // backward part of the FwdBwd algorithm (the .m, .a values must not be
          // changed before the centered evaluation).
 
          if(using_nat)
          {
        _f.expr()->bwd_eval(v); // recursive backward from the root to the leaves
        _f.intersect_from_args(v, x...); // updating input values
          }
 
 
          //if(_both || _both_then_one || md < _freq)
          //  if(_with_centered_form && val_expr.def_domain && !val_expr.da.is_unbounded() && val_expr.da.size() != 0)
          if(using_centered)
          {
            // todo: the above condition !val_expr.da.is_unbounded() should not be necesary,
            // possible bug in MulOp in case of unbounded domain?
            using X0 = std::tuple_element_t<0, std::tuple<X...>>;
 
            if constexpr(sizeof...(X) == 1 && std::is_same_v<X0,IntervalVector>)
            {
              X0& x_ = std::get<0>(std::tie(x...));
              X0 x_mid = X0(x_.mid());
 
              assert(val_expr.a.size() == val_expr.m.size());
 
              if constexpr(std::is_same_v<Y,IntervalMatrix>)
              {
                std::cout << "CtcInverse: matrices expressions not (yet) supported with centered form" << std::endl;
              }
 
              else
              {
                IntervalVector p = x_ - x_mid;
                IntervalVector fm { val_expr.a - val_expr.m };
                MulOp::bwd(fm, val_expr.da, p);
                x_ &= p + x_mid;
              }
            }
 
            else
            {
              // Centered form not (yet) implemented for multi-nonvector-arguments
            }
          }
          
        // [4/4] Backward evaluation
 
        /*X0b x0_after = std::get<0>(std::tie(x...));
        if(x0_after.is_empty())
          _data.push_back({x0_before.max_diam(), 0.});
        else
          _data.push_back({x0_before.max_diam(), x0_before.volume()/x0_after.volume()});*/
 
      }
 
      const AnalyticFunction<typename ExprType<Y>::Type>& function() const
      {
        return _f;
      }
 
      //mutable std::vector<std::pair<double,double>> _data;
 
      double _freq = 0.;
      bool _both = true;
      bool _both_then_one = true;
 
    protected:
 
      const AnalyticFunction<typename ExprType<Y>::Type> _f;
      const Collection<CtcBase<Y>> _ctc_y;
      bool _with_centered_form;
      bool _is_not_in = false;
  };
  
 
  // Template deduction guides
 
  // ScalarType
 
    CtcInverse(const AnalyticFunction<ScalarType>&, std::initializer_list<double>, bool = true, bool = false) -> 
      CtcInverse<Interval,IntervalVector>;
 
    template<typename Y>
    CtcInverse(const AnalyticFunction<ScalarType>&, std::initializer_list<Y>, bool = true, bool = false) -> 
      CtcInverse<Interval,IntervalVector>;
 
    CtcInverse(const AnalyticFunction<ScalarType>&, const Interval&, bool = true, bool = false) -> 
      CtcInverse<Interval,IntervalVector>;
 
    CtcInverse(const AnalyticFunction<ScalarType>&, double, bool = true, bool = false) -> 
      CtcInverse<Interval,IntervalVector>;
 
    template<typename C>
      requires IsCtcBaseOrPtr<C,Interval>
    CtcInverse(const AnalyticFunction<ScalarType>&, const C&, bool = true, bool = false) -> 
      CtcInverse<Interval,IntervalVector>;
 
  // VectorType
 
    CtcInverse(const AnalyticFunction<VectorType>&, std::initializer_list<double>, bool = true, bool = false) -> 
      CtcInverse<IntervalVector,IntervalVector>;
 
    CtcInverse(const AnalyticFunction<VectorType>&, std::initializer_list<std::initializer_list<double>>, bool = true, bool = false) -> 
      CtcInverse<IntervalVector,IntervalVector>;
 
    CtcInverse(const AnalyticFunction<VectorType>&, const Vector&, bool = true, bool = false) -> 
      CtcInverse<IntervalVector,IntervalVector>;
 
    CtcInverse(const AnalyticFunction<VectorType>&, const IntervalVector&, bool = true, 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<IntervalVector,IntervalVector>;
 
    template<typename C>
      requires IsCtcBaseOrPtr<C,IntervalVector>
    CtcInverse(const AnalyticFunction<VectorType>&, const C&, bool = true, bool = false) -> 
      CtcInverse<IntervalVector,IntervalVector>;
 
  // MatrixType
        
    CtcInverse(const AnalyticFunction<MatrixType>&, std::initializer_list<std::initializer_list<double>>, bool = true, bool = false) -> 
      CtcInverse<IntervalMatrix,IntervalVector>;
 
    CtcInverse(const AnalyticFunction<MatrixType>&, std::initializer_list<std::initializer_list<std::initializer_list<double>>>, bool = true, 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<IntervalMatrix,IntervalVector>;
 
    template<typename C>
      requires IsCtcBaseOrPtr<C,IntervalMatrix>
    CtcInverse(const AnalyticFunction<MatrixType>&, const C&, bool = true, bool = false) -> 
      CtcInverse<IntervalMatrix,IntervalVector>;
}