OctaSym can now be involved in analytic expressions

doc/manual/manual/functions/analytic/analytic_operators.rst

@@ -139,6 +139,12 @@ When operators are available for operations 2--3, then an ``AnalyticFunction`` c
 
 .. |MatrixOp| replace:: :raw-html:`<a href="https://github.com/codac-team/codac/blob/codac2/src/core/operators/codac2_mat.h"><code class="docutils literal notranslate"><span class="pre">MatrixOp</span></code></a>`
 
+.. |OctaSymOp| replace:: :raw-html:`<a href="https://github.com/codac-team/codac/blob/codac2/src/core/actions/codac2_OctaSym_operator.h"><code class="docutils literal notranslate"><span class="pre">OctaSymOp</span></code></a>`
+
+.. |TrajOp| replace:: :raw-html:`<a href="https://github.com/codac-team/codac/blob/codac2/src/core/trajectory/codac2_Traj_operator.h"><code class="docutils literal notranslate"><span class="pre">TrajectoryOp</span></code></a>`
+
+.. |TubeOp| replace:: :raw-html:`<a href="https://github.com/codac-team/codac/blob/codac2/src/core/domains/tube/codac2_Tube_operator.h"><code class="docutils literal notranslate"><span class="pre">TubeOp</span></code></a>`
+
 Only the :bg-ok:`✓` operators are supported at the moment.
 If you notice any mathematical operators missing from the list below, feel free to contribute to the library. You can submit your suggestions or pull requests on the `GitHub repository of Codac <https://github.com/codac-team/codac>`_.
 
@@ -266,6 +272,14 @@ If you notice any mathematical operators missing from the list below, feel free
   +-----------------------------------------------------+----------------------+---------------+-------------------------------------+--------+--------+-------+------------+
   | :math:`\left(\mathbf{x}_1,\mathbf{x}_2,\dots\right)`| ``mat(x1,x2,...)``   | |MatrixOp|    | ``x1``, ``...``: vector             ||okk|   ||okk|   ||okk|  ||nok|       |
   +-----------------------------------------------------+----------------------+---------------+-------------------------------------+--------+--------+-------+------------+
+  | :math:`\sigma(\mathbf{x})`                          | ``s(x)``             | |OctaSymOp|   | ``s``: action, ``x``: vector        ||okk|   ||okk|   ||okk|  ||okk|       |
+  +-----------------------------------------------------+----------------------+---------------+-------------------------------------+--------+--------+-------+------------+
+  | :bg-title:`Temporal operations`                                                                                                                                         |
+  +-----------------------------------------------------+----------------------+---------------+-------------------------------------+--------+--------+-------+------------+
+  | :math:`\mathbf{x}(t)`                               | ``x(t)``             | |TrajOp|      | ``x``: trajectory, ``t``: scalar    ||okk|   ||nok|   ||nok|  ||nok|       |
+  +-----------------------------------------------------+----------------------+---------------+-------------------------------------+--------+--------+-------+------------+
+  | :math:`[\mathbf{x}](t)`                             | ``x(t)``             | |TubeOp|      | ``x``: tube, ``t``: scalar          ||okk|   ||nok|   ||nok|  ||nok|       |
+  +-----------------------------------------------------+----------------------+---------------+-------------------------------------+--------+--------+-------+------------+
 
 | Note: the operator :math:`\det` is only available for :math:`1\times 1` and :math:`2\times 2` matrices.
 | Note: the operator :math:`\bmod` is only available for real periods (double precision), interval periods are not yet supported.

python/src/core/actions/codac2_py_OctaSym.cpp

@@ -14,10 +14,12 @@
 #include "codac2_py_Ctc.h"
 #include "codac2_py_Sep.h"
 #include <codac2_OctaSym.h>
+#include <codac2_OctaSym_operator.h>
 #include <codac2_template_tools.h>
 #include <codac2_CtcAction.h>
 #include <codac2_SepAction.h>
 #include "codac2_py_OctaSym_docs.h" // Generated file from Doxygen XML (doxygen2docstring.py):
+#include "codac2_py_cast.h"
 
 using namespace std;
 using namespace codac2;
@@ -75,6 +77,10 @@ void export_OctaSym(py::module& m)
       SAMPLEDTRAJ_T_OCTASYM_OPERATORCALL_CONST_SAMPLEDTRAJ_T_REF_CONST,
       "x"_a)
 
+    .def("__call__", [](const OctaSym& a, const VectorExpr& x) { return a.operator()(x); },
+      VECTOREXPR_OCTASYM_OPERATORCALL_CONST_V_REF_CONST,
+      "x"_a)
+
     .def("__repr__", [](const OctaSym& s) {
           std::ostringstream stream;
           stream << s;

src/core/CMakeLists.txt

@@ -6,6 +6,7 @@
 
     ${CMAKE_CURRENT_SOURCE_DIR}/actions/codac2_OctaSym.cpp
     ${CMAKE_CURRENT_SOURCE_DIR}/actions/codac2_OctaSym.h
+    ${CMAKE_CURRENT_SOURCE_DIR}/actions/codac2_OctaSym_operator.h
 
     ${CMAKE_CURRENT_SOURCE_DIR}/contractors/codac2_Ctc.h
     ${CMAKE_CURRENT_SOURCE_DIR}/contractors/codac2_CtcAction.cpp

src/core/actions/codac2_OctaSym.cpp

@@ -41,14 +41,14 @@ OctaSym OctaSym::operator*(const OctaSym& s) const
   assert_release(size() == s.size());
   OctaSym a(*this);
   for(size_t i = 0 ; i < a.size() ; i++)
-    a[i] = _sign(s[i])*(*this)[std::abs((int)s[i])-1];
+    a[i] = sign(s[i])*(*this)[std::abs((int)s[i])-1];
   return a;
 }
 
 Matrix OctaSym::permutation_matrix() const
 { 
   Matrix m = Matrix::zero(size(),size());
   for(size_t i = 0 ; i < size() ; i++)
-    m(i,std::abs((*this)[i])-1) = _sign((*this)[i]);
+    m(i,std::abs((*this)[i])-1) = sign((*this)[i]);
   return m;
 }

src/core/actions/codac2_OctaSym.h

@@ -23,6 +23,7 @@ namespace codac2
   class SepBase;
   class SepAction;
   class SetExpr;
+  class OctaSymOp;
 
   /**
    * \class Action
@@ -46,19 +47,14 @@ namespace codac2
 
       Matrix permutation_matrix() const;
 
-      int _sign(int a) const
-      {
-        return (a > 0) ? 1 : ((a < 0) ? -1 : 0);
-      }
-
       template<typename Derived>
         requires (Derived::ColsAtCompileTime == 1)
       Mat<typename Derived::Scalar,-1,1> operator()(const Eigen::MatrixBase<Derived>& x) const
       {
         assert_release(x.size() == (Index)size());
         Mat<typename Derived::Scalar,-1,1> x_(x);
         for(size_t i = 0 ; i < size() ; i++)
-          x_[i] = _sign((*this)[i])*x[std::abs((*this)[i])-1];
+          x_[i] = sign((*this)[i])*x[std::abs((*this)[i])-1];
         return x_;
       }
 
@@ -84,6 +80,18 @@ namespace codac2
         return y;
       }
 
+      template<typename V>
+        // To avoid ambiguity with operator()(const Eigen::MatrixBase<Derived>& x):
+        requires (std::is_same_v<V,VectorExpr> || std::is_same_v<V,VectorVar>)
+      inline VectorExpr operator()(const V& x1) const
+      {
+        if constexpr(std::is_same_v<V,VectorExpr>)
+          assert_release((Index)this->size() == x1->output_shape().first);
+        else
+          assert_release((Index)this->size() == x1.output_shape().first);
+        return { std::make_shared<AnalyticOperationExpr<OctaSymOp,VectorType,VectorType>>(*this, x1) };
+      }
+
       friend std::ostream& operator<<(std::ostream& str, const OctaSym& s)
       {
         str << "(";

src/core/actions/codac2_OctaSym_operator.h

@@ -0,0 +1,129 @@
+/** 
+ *  \file codac2_OctaSym_operator.h
+ * ----------------------------------------------------------------------------
+ *  \date       2025
+ *  \author     Simon Rohou
+ *  \copyright  Copyright 2025 Codac Team
+ *  \license    GNU Lesser General Public License (LGPL)
+ */
+
+#pragma once
+
+#include "codac2_OctaSym.h"
+
+namespace codac2
+{
+  struct OctaSymOp
+  {
+    template<typename X1>
+    static inline std::string str(const X1& x1)
+    {
+      return "sym(" + x1->str() + ")";
+    }
+
+    template<typename X1>
+    static inline std::pair<Index,Index> output_shape([[maybe_unused]] const X1& s1)
+    {
+      return s1->output_shape();
+    }
+
+    static inline IntervalVector fwd(const OctaSym& s, const IntervalVector& x1)
+    {
+      assert((Index)s.size() == x1.size());
+      return s(x1);
+    }
+
+    static inline VectorType fwd_natural(const OctaSym& s, const VectorType& x1)
+    {
+      assert((Index)s.size() == x1.m.size());
+      return {
+        fwd(s, x1.a),
+        x1.def_domain
+      };
+    }
+
+    static inline VectorType fwd_centered(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++)
+        da.row(i) = sign(s[i])*x1.da.row(std::abs(s[i])-1);
+
+      return {
+        fwd(s, x1.m),
+        fwd(s, x1.a),
+        da,
+        x1.def_domain
+      };
+    }
+
+    static inline void bwd(const OctaSym& s, const IntervalVector& y, IntervalVector& x1)
+    {
+      assert((Index)s.size() == y.size() && (Index)s.size() == x1.size());
+      x1 &= s.invert()(y);
+    }
+  };
+
+
+  template<>
+  class AnalyticOperationExpr<OctaSymOp,VectorType,VectorType>
+    : public AnalyticExpr<VectorType>, public OperationExprBase<AnalyticExpr<VectorType>>
+  {
+    public:
+
+      AnalyticOperationExpr(const OctaSym& s, const VectorExpr& x1)
+        : OperationExprBase<AnalyticExpr<VectorType>>(x1), _s(s)
+      { }
+
+      std::shared_ptr<ExprBase> copy() const
+      {
+        return std::make_shared<AnalyticOperationExpr<OctaSymOp,VectorType,VectorType>>(*this);
+      }
+
+      void replace_arg(const ExprID& old_arg_id, const std::shared_ptr<ExprBase>& new_expr)
+      {
+        return OperationExprBase<AnalyticExpr<VectorType>>::replace_arg(old_arg_id, new_expr);
+      }
+
+      VectorType fwd_eval(ValuesMap& v, Index total_input_size, bool natural_eval) 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)));
+        else
+          return AnalyticExpr<VectorType>::init_value(
+            v, OctaSymOp::fwd_centered(_s, std::get<0>(this->_x)->fwd_eval(v, total_input_size, natural_eval)));
+      }
+
+      void bwd_eval(ValuesMap& v) const
+      {
+        OctaSymOp::bwd(_s, AnalyticExpr<VectorType>::value(v).a, std::get<0>(this->_x)->value(v).a);
+        std::get<0>(this->_x)->bwd_eval(v);
+      }
+
+      std::pair<Index,Index> output_shape() const {
+        return { _s.size(), 1 };
+      }
+
+      virtual bool belongs_to_args_list(const FunctionArgsList& args) const
+      {
+        return std::get<0>(this->_x)->belongs_to_args_list(args);
+      }
+
+      std::string str(bool in_parentheses = false) const
+      {
+        std::string s = "S"; // user cannot (yet) specify a name for the symmetry
+        return in_parentheses ? "(" + s + ")" : s;
+      }
+
+      virtual bool is_str_leaf() const
+      {
+        return true;
+      }
+
+    protected:
+
+      const OctaSym _s;
+  };
+}

src/core/tools/codac2_math.h

@@ -30,9 +30,11 @@ namespace codac2
 
   constexpr double PI = std::numbers::pi; // Need C++20
 
-  inline int sign(double x)
+  template<typename T>
+  requires std::is_arithmetic_v<T>
+  inline constexpr int sign(T x)
   {
-    return (x > 0) ? 1 : ((x < 0) ? -1 : 0);
+    return (x > T(0)) - (x < T(0));
   }
 
   inline int integer(double x)

tests/core/actions/codac2_tests_OctaSym.cpp

@@ -9,6 +9,7 @@
 
 #include <catch2/catch_test_macros.hpp>
 #include <codac2_OctaSym.h>
+#include <codac2_OctaSym_operator.h>
 
 using namespace std;
 using namespace codac2;
@@ -28,4 +29,14 @@ TEST_CASE("OctaSym")
 
   OctaSym c({-2,1,3});
   CHECK(c.permutation_matrix() == Matrix({{0,-1,0},{1,0,0},{0,0,1}}));
+}
+
+TEST_CASE("OctaSym as operator")
+{
+  OctaSym a({3,1,-2});
+  VectorVar x(3);
+  AnalyticFunction f({x}, a(2*x));
+  CHECK(f.eval(IntervalVector({{1},{2},{3}})) == IntervalVector({{6},{2},{-4}}));
+  CHECK(f.eval(IntervalVector({{-oo,oo},{-oo,oo},{-oo,oo}})) == IntervalVector(3));
+  CHECK(f.eval(IntervalVector::empty(3)) == IntervalVector::empty(3));
 }

tests/core/actions/codac2_tests_OctaSym.py

@@ -28,5 +28,14 @@ def test_OctaSym(self):
     c = OctaSym([-2,1,3])
     self.assertTrue(c.permutation_matrix() == Matrix([[0,-1,0],[1,0,0],[0,0,1]]))
 
+  def test_OctaSym_as_operator(self):
+
+    a = OctaSym([3,1,-2])
+    x = VectorVar(3)
+    f = AnalyticFunction([x], a(2*x))
+    self.assertTrue(f.eval(IntervalVector([[1],[2],[3]])) == IntervalVector([[6],[2],[-4]]))
+    self.assertTrue(f.eval(IntervalVector([[-oo,oo],[-oo,oo],[-oo,oo]])) == IntervalVector(3))
+    self.assertTrue(f.eval(IntervalVector.empty(3)) == IntervalVector.empty(3))
+
 if __name__ ==  '__main__':
   unittest.main()