Skip to content

Commit 8daa7f2

Browse files
authored
Merge pull request #301 from SimonRohou/codac2_dev
OctaSym can now be involved in analytic expressions
2 parents 4052eef + 862bcb0 commit 8daa7f2

File tree

9 files changed

+190
-10
lines changed

9 files changed

+190
-10
lines changed

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

Lines changed: 14 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -139,6 +139,12 @@ When operators are available for operations 2--3, then an ``AnalyticFunction`` c
139139

140140
.. |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>`
141141

142+
.. |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>`
143+
144+
.. |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>`
145+
146+
.. |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>`
147+
142148
Only the :bg-ok:`` operators are supported at the moment.
143149
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>`_.
144150

@@ -266,6 +272,14 @@ If you notice any mathematical operators missing from the list below, feel free
266272
+-----------------------------------------------------+----------------------+---------------+-------------------------------------+--------+--------+-------+------------+
267273
| :math:`\left(\mathbf{x}_1,\mathbf{x}_2,\dots\right)`| ``mat(x1,x2,...)`` | |MatrixOp| | ``x1``, ``...``: vector ||okk| ||okk| ||okk| ||nok| |
268274
+-----------------------------------------------------+----------------------+---------------+-------------------------------------+--------+--------+-------+------------+
275+
| :math:`\sigma(\mathbf{x})` | ``s(x)`` | |OctaSymOp| | ``s``: action, ``x``: vector ||okk| ||okk| ||okk| ||okk| |
276+
+-----------------------------------------------------+----------------------+---------------+-------------------------------------+--------+--------+-------+------------+
277+
| :bg-title:`Temporal operations` |
278+
+-----------------------------------------------------+----------------------+---------------+-------------------------------------+--------+--------+-------+------------+
279+
| :math:`\mathbf{x}(t)` | ``x(t)`` | |TrajOp| | ``x``: trajectory, ``t``: scalar ||okk| ||nok| ||nok| ||nok| |
280+
+-----------------------------------------------------+----------------------+---------------+-------------------------------------+--------+--------+-------+------------+
281+
| :math:`[\mathbf{x}](t)` | ``x(t)`` | |TubeOp| | ``x``: tube, ``t``: scalar ||okk| ||nok| ||nok| ||nok| |
282+
+-----------------------------------------------------+----------------------+---------------+-------------------------------------+--------+--------+-------+------------+
269283

270284
| Note: the operator :math:`\det` is only available for :math:`1\times 1` and :math:`2\times 2` matrices.
271285
| 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

Lines changed: 6 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -14,10 +14,12 @@
1414
#include "codac2_py_Ctc.h"
1515
#include "codac2_py_Sep.h"
1616
#include <codac2_OctaSym.h>
17+
#include <codac2_OctaSym_operator.h>
1718
#include <codac2_template_tools.h>
1819
#include <codac2_CtcAction.h>
1920
#include <codac2_SepAction.h>
2021
#include "codac2_py_OctaSym_docs.h" // Generated file from Doxygen XML (doxygen2docstring.py):
22+
#include "codac2_py_cast.h"
2123

2224
using namespace std;
2325
using namespace codac2;
@@ -75,6 +77,10 @@ void export_OctaSym(py::module& m)
7577
SAMPLEDTRAJ_T_OCTASYM_OPERATORCALL_CONST_SAMPLEDTRAJ_T_REF_CONST,
7678
"x"_a)
7779

80+
.def("__call__", [](const OctaSym& a, const VectorExpr& x) { return a.operator()(x); },
81+
VECTOREXPR_OCTASYM_OPERATORCALL_CONST_V_REF_CONST,
82+
"x"_a)
83+
7884
.def("__repr__", [](const OctaSym& s) {
7985
std::ostringstream stream;
8086
stream << s;

src/core/CMakeLists.txt

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -6,6 +6,7 @@
66

77
${CMAKE_CURRENT_SOURCE_DIR}/actions/codac2_OctaSym.cpp
88
${CMAKE_CURRENT_SOURCE_DIR}/actions/codac2_OctaSym.h
9+
${CMAKE_CURRENT_SOURCE_DIR}/actions/codac2_OctaSym_operator.h
910

1011
${CMAKE_CURRENT_SOURCE_DIR}/contractors/codac2_Ctc.h
1112
${CMAKE_CURRENT_SOURCE_DIR}/contractors/codac2_CtcAction.cpp

src/core/actions/codac2_OctaSym.cpp

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -41,14 +41,14 @@ OctaSym OctaSym::operator*(const OctaSym& s) const
4141
assert_release(size() == s.size());
4242
OctaSym a(*this);
4343
for(size_t i = 0 ; i < a.size() ; i++)
44-
a[i] = _sign(s[i])*(*this)[std::abs((int)s[i])-1];
44+
a[i] = sign(s[i])*(*this)[std::abs((int)s[i])-1];
4545
return a;
4646
}
4747

4848
Matrix OctaSym::permutation_matrix() const
4949
{
5050
Matrix m = Matrix::zero(size(),size());
5151
for(size_t i = 0 ; i < size() ; i++)
52-
m(i,std::abs((*this)[i])-1) = _sign((*this)[i]);
52+
m(i,std::abs((*this)[i])-1) = sign((*this)[i]);
5353
return m;
5454
}

src/core/actions/codac2_OctaSym.h

Lines changed: 14 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -23,6 +23,7 @@ namespace codac2
2323
class SepBase;
2424
class SepAction;
2525
class SetExpr;
26+
class OctaSymOp;
2627

2728
/**
2829
* \class Action
@@ -46,19 +47,14 @@ namespace codac2
4647

4748
Matrix permutation_matrix() const;
4849

49-
int _sign(int a) const
50-
{
51-
return (a > 0) ? 1 : ((a < 0) ? -1 : 0);
52-
}
53-
5450
template<typename Derived>
5551
requires (Derived::ColsAtCompileTime == 1)
5652
Mat<typename Derived::Scalar,-1,1> operator()(const Eigen::MatrixBase<Derived>& x) const
5753
{
5854
assert_release(x.size() == (Index)size());
5955
Mat<typename Derived::Scalar,-1,1> x_(x);
6056
for(size_t i = 0 ; i < size() ; i++)
61-
x_[i] = _sign((*this)[i])*x[std::abs((*this)[i])-1];
57+
x_[i] = sign((*this)[i])*x[std::abs((*this)[i])-1];
6258
return x_;
6359
}
6460

@@ -84,6 +80,18 @@ namespace codac2
8480
return y;
8581
}
8682

83+
template<typename V>
84+
// To avoid ambiguity with operator()(const Eigen::MatrixBase<Derived>& x):
85+
requires (std::is_same_v<V,VectorExpr> || std::is_same_v<V,VectorVar>)
86+
inline VectorExpr operator()(const V& x1) const
87+
{
88+
if constexpr(std::is_same_v<V,VectorExpr>)
89+
assert_release((Index)this->size() == x1->output_shape().first);
90+
else
91+
assert_release((Index)this->size() == x1.output_shape().first);
92+
return { std::make_shared<AnalyticOperationExpr<OctaSymOp,VectorType,VectorType>>(*this, x1) };
93+
}
94+
8795
friend std::ostream& operator<<(std::ostream& str, const OctaSym& s)
8896
{
8997
str << "(";
Lines changed: 129 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,129 @@
1+
/**
2+
* \file codac2_OctaSym_operator.h
3+
* ----------------------------------------------------------------------------
4+
* \date 2025
5+
* \author Simon Rohou
6+
* \copyright Copyright 2025 Codac Team
7+
* \license GNU Lesser General Public License (LGPL)
8+
*/
9+
10+
#pragma once
11+
12+
#include "codac2_OctaSym.h"
13+
14+
namespace codac2
15+
{
16+
struct OctaSymOp
17+
{
18+
template<typename X1>
19+
static inline std::string str(const X1& x1)
20+
{
21+
return "sym(" + x1->str() + ")";
22+
}
23+
24+
template<typename X1>
25+
static inline std::pair<Index,Index> output_shape([[maybe_unused]] const X1& s1)
26+
{
27+
return s1->output_shape();
28+
}
29+
30+
static inline IntervalVector fwd(const OctaSym& s, const IntervalVector& x1)
31+
{
32+
assert((Index)s.size() == x1.size());
33+
return s(x1);
34+
}
35+
36+
static inline VectorType fwd_natural(const OctaSym& s, const VectorType& x1)
37+
{
38+
assert((Index)s.size() == x1.m.size());
39+
return {
40+
fwd(s, x1.a),
41+
x1.def_domain
42+
};
43+
}
44+
45+
static inline VectorType fwd_centered(const OctaSym& s, const VectorType& x1)
46+
{
47+
assert((Index)s.size() == x1.m.size());
48+
49+
auto da = x1.da;
50+
for(size_t i = 0 ; i < s.size() ; i++)
51+
da.row(i) = sign(s[i])*x1.da.row(std::abs(s[i])-1);
52+
53+
return {
54+
fwd(s, x1.m),
55+
fwd(s, x1.a),
56+
da,
57+
x1.def_domain
58+
};
59+
}
60+
61+
static inline void bwd(const OctaSym& s, const IntervalVector& y, IntervalVector& x1)
62+
{
63+
assert((Index)s.size() == y.size() && (Index)s.size() == x1.size());
64+
x1 &= s.invert()(y);
65+
}
66+
};
67+
68+
69+
template<>
70+
class AnalyticOperationExpr<OctaSymOp,VectorType,VectorType>
71+
: public AnalyticExpr<VectorType>, public OperationExprBase<AnalyticExpr<VectorType>>
72+
{
73+
public:
74+
75+
AnalyticOperationExpr(const OctaSym& s, const VectorExpr& x1)
76+
: OperationExprBase<AnalyticExpr<VectorType>>(x1), _s(s)
77+
{ }
78+
79+
std::shared_ptr<ExprBase> copy() const
80+
{
81+
return std::make_shared<AnalyticOperationExpr<OctaSymOp,VectorType,VectorType>>(*this);
82+
}
83+
84+
void replace_arg(const ExprID& old_arg_id, const std::shared_ptr<ExprBase>& new_expr)
85+
{
86+
return OperationExprBase<AnalyticExpr<VectorType>>::replace_arg(old_arg_id, new_expr);
87+
}
88+
89+
VectorType fwd_eval(ValuesMap& v, Index total_input_size, bool natural_eval) const
90+
{
91+
if(natural_eval)
92+
return AnalyticExpr<VectorType>::init_value(
93+
v, OctaSymOp::fwd_natural(_s, std::get<0>(this->_x)->fwd_eval(v, total_input_size, natural_eval)));
94+
else
95+
return AnalyticExpr<VectorType>::init_value(
96+
v, OctaSymOp::fwd_centered(_s, std::get<0>(this->_x)->fwd_eval(v, total_input_size, natural_eval)));
97+
}
98+
99+
void bwd_eval(ValuesMap& v) const
100+
{
101+
OctaSymOp::bwd(_s, AnalyticExpr<VectorType>::value(v).a, std::get<0>(this->_x)->value(v).a);
102+
std::get<0>(this->_x)->bwd_eval(v);
103+
}
104+
105+
std::pair<Index,Index> output_shape() const {
106+
return { _s.size(), 1 };
107+
}
108+
109+
virtual bool belongs_to_args_list(const FunctionArgsList& args) const
110+
{
111+
return std::get<0>(this->_x)->belongs_to_args_list(args);
112+
}
113+
114+
std::string str(bool in_parentheses = false) const
115+
{
116+
std::string s = "S"; // user cannot (yet) specify a name for the symmetry
117+
return in_parentheses ? "(" + s + ")" : s;
118+
}
119+
120+
virtual bool is_str_leaf() const
121+
{
122+
return true;
123+
}
124+
125+
protected:
126+
127+
const OctaSym _s;
128+
};
129+
}

src/core/tools/codac2_math.h

Lines changed: 4 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -30,9 +30,11 @@ namespace codac2
3030

3131
constexpr double PI = std::numbers::pi; // Need C++20
3232

33-
inline int sign(double x)
33+
template<typename T>
34+
requires std::is_arithmetic_v<T>
35+
inline constexpr int sign(T x)
3436
{
35-
return (x > 0) ? 1 : ((x < 0) ? -1 : 0);
37+
return (x > T(0)) - (x < T(0));
3638
}
3739

3840
inline int integer(double x)

tests/core/actions/codac2_tests_OctaSym.cpp

Lines changed: 11 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -9,6 +9,7 @@
99

1010
#include <catch2/catch_test_macros.hpp>
1111
#include <codac2_OctaSym.h>
12+
#include <codac2_OctaSym_operator.h>
1213

1314
using namespace std;
1415
using namespace codac2;
@@ -28,4 +29,14 @@ TEST_CASE("OctaSym")
2829

2930
OctaSym c({-2,1,3});
3031
CHECK(c.permutation_matrix() == Matrix({{0,-1,0},{1,0,0},{0,0,1}}));
32+
}
33+
34+
TEST_CASE("OctaSym as operator")
35+
{
36+
OctaSym a({3,1,-2});
37+
VectorVar x(3);
38+
AnalyticFunction f({x}, a(2*x));
39+
CHECK(f.eval(IntervalVector({{1},{2},{3}})) == IntervalVector({{6},{2},{-4}}));
40+
CHECK(f.eval(IntervalVector({{-oo,oo},{-oo,oo},{-oo,oo}})) == IntervalVector(3));
41+
CHECK(f.eval(IntervalVector::empty(3)) == IntervalVector::empty(3));
3142
}

tests/core/actions/codac2_tests_OctaSym.py

Lines changed: 9 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -28,5 +28,14 @@ def test_OctaSym(self):
2828
c = OctaSym([-2,1,3])
2929
self.assertTrue(c.permutation_matrix() == Matrix([[0,-1,0],[1,0,0],[0,0,1]]))
3030

31+
def test_OctaSym_as_operator(self):
32+
33+
a = OctaSym([3,1,-2])
34+
x = VectorVar(3)
35+
f = AnalyticFunction([x], a(2*x))
36+
self.assertTrue(f.eval(IntervalVector([[1],[2],[3]])) == IntervalVector([[6],[2],[-4]]))
37+
self.assertTrue(f.eval(IntervalVector([[-oo,oo],[-oo,oo],[-oo,oo]])) == IntervalVector(3))
38+
self.assertTrue(f.eval(IntervalVector.empty(3)) == IntervalVector.empty(3))
39+
3140
if __name__ == '__main__':
3241
unittest.main()

0 commit comments

Comments
 (0)