diff --git a/.github/workflows/macosmatrix.yml b/.github/workflows/macosmatrix.yml index 90a45049d..da66789bb 100644 --- a/.github/workflows/macosmatrix.yml +++ b/.github/workflows/macosmatrix.yml @@ -54,7 +54,7 @@ jobs: # if: runner.os=='macOS' - run: brew install catch2 # Issues with binary packages when cross-compiling... if: (runner.os=='macOS')&&(matrix.cfg.cross!=true) - - run: brew install graphviz ; brew install --formula doxygen ; python -m pip install --upgrade pip ; pip install --upgrade wheel setuptools sphinx breathe sphinx_rtd_theme sphinx-tabs sphinx-issues sphinx-reredirects furo sphinx-math-dollar + - run: brew install graphviz ; brew install --formula doxygen ; python -m pip install --upgrade pip ; pip install --upgrade wheel setuptools sphinx breathe sphinx_rtd_theme sphinx-tabs sphinx-issues sphinx-reredirects furo sphinx-math-dollar sympy if: runner.os=='macOS' - run: | wget https://github.com/lebarsfa/ibex-lib/releases/download/ibex-2.8.9.20250626/ibex_${{ matrix.cfg.arch }}_${{ matrix.cfg.runtime }}.zip --no-check-certificate -nv @@ -80,7 +80,7 @@ jobs: - run: | pip install --no-deps --no-index *.whl python -c "import sys; print(sys.version)" ; python examples/02_centered_form/main.py - pip install numpy --prefer-binary + pip install numpy sympy --prefer-binary python -m unittest discover codac.tests cd build && ctest -C Release -V --output-on-failure cd .. diff --git a/.github/workflows/tests.yml b/.github/workflows/tests.yml index 8f644ec80..52fe8db4a 100644 --- a/.github/workflows/tests.yml +++ b/.github/workflows/tests.yml @@ -45,10 +45,10 @@ jobs: sudo sh -c 'echo "deb [trusted=yes] https://webperso.ensta.fr/packages/$(if [ -z "$(. /etc/os-release && echo $UBUNTU_CODENAME)" ]; then echo debian/$(. /etc/os-release && echo $VERSION_CODENAME); else echo ubuntu/$(. /etc/os-release && echo $UBUNTU_CODENAME); fi) ./" > /etc/apt/sources.list.d/ensta-bretagne.list' sudo apt update - sudo apt-get -y install flex bison catch2 # libeigen3-dev + sudo apt-get -y install flex bison catch2 pybind11-dev # libeigen3-dev # For documentation - pip install sphinx breathe sphinx-issues sphinx-tabs sphinx_rtd_theme + pip install sphinx breathe sphinx-issues sphinx-tabs sphinx_rtd_theme sympy sudo apt-get -y install doxygen graphviz # Doxygen might need to be upgraded @@ -61,7 +61,7 @@ jobs: #pip install --upgrade pyibex==1.8.0 # Fort Python tests - pip install numpy + pip install numpy sympy pwd ls @@ -121,7 +121,7 @@ jobs: cd ../examples cd 01_batman/ - mkdir build ; cd build ; cmake -DCMAKE_PREFIX_PATH="~/ibex-lib/build_install;~/codac/build_install" -DCMAKE_BUILD_TYPE=Debug .. ; make ; ./codac_example + mkdir build ; cd build ; cmake -DCMAKE_PREFIX_PATH="$HOME/ibex-lib/build_install;$HOME/codac/build_install" -DCMAKE_BUILD_TYPE=Debug .. ; make ; ./codac_example cd ../../02_centered_form/ - mkdir build ; cd build ; cmake -DCMAKE_PREFIX_PATH="~/ibex-lib/build_install;~/codac/build_install" -DCMAKE_BUILD_TYPE=Debug .. ; make ; ./codac_example \ No newline at end of file + mkdir build ; cd build ; cmake -DCMAKE_PREFIX_PATH="$HOME/ibex-lib/build_install;$HOME/codac/build_install" -DCMAKE_BUILD_TYPE=Debug .. ; make ; ./codac_example \ No newline at end of file diff --git a/.github/workflows/vcmatrix.yml b/.github/workflows/vcmatrix.yml index 14aa67dbc..62b58767f 100644 --- a/.github/workflows/vcmatrix.yml +++ b/.github/workflows/vcmatrix.yml @@ -56,7 +56,7 @@ jobs: if: runner.os=='Windows' #- run: choco install -y -r --no-progress eigen --version=3.4.0.20240224 ${{ matrix.cfg.choco_flags }} # if: runner.os=='Windows' - - run: choco install -y -r --no-progress graphviz doxygen.install & python -m pip install --upgrade pip & pip install --upgrade wheel setuptools sphinx breathe sphinx-issues sphinx-tabs sphinx_rtd_theme sphinx-reredirects furo sphinx-math-dollar + - run: choco install -y -r --no-progress graphviz doxygen.install & python -m pip install --upgrade pip & pip install --upgrade wheel setuptools sphinx breathe sphinx-issues sphinx-tabs sphinx_rtd_theme sphinx-reredirects furo sphinx-math-dollar sympy if: runner.os=='Windows' - run: | wget https://github.com/lebarsfa/ibex-lib/releases/download/ibex-2.8.9.20250626/ibex.2.8.9.20250626.nupkg --no-check-certificate -nv @@ -80,7 +80,7 @@ jobs: - run: | pip install --no-deps --no-index *.whl python -c "import sys; print(sys.version)" ; python examples/02_centered_form/main.py - pip install numpy --prefer-binary + pip install numpy sympy --prefer-binary python -m unittest discover codac.tests cd build && ctest -C Release -V --output-on-failure cd .. diff --git a/doc/manual/manual/extensions/index.rst b/doc/manual/manual/extensions/index.rst index 8652b2eaf..6d5c3fb18 100644 --- a/doc/manual/manual/extensions/index.rst +++ b/doc/manual/manual/extensions/index.rst @@ -5,7 +5,8 @@ Codac extensions .. toctree:: + sympy/index.rst capd/index.rst + .. ibex/index.rst -.. sympy/index.rst .. unsupported/index.rst \ No newline at end of file diff --git a/doc/manual/manual/extensions/sympy/index.rst b/doc/manual/manual/extensions/sympy/index.rst new file mode 100644 index 000000000..3ec521785 --- /dev/null +++ b/doc/manual/manual/extensions/sympy/index.rst @@ -0,0 +1,487 @@ +.. _sec-extensions-sympy: + +SymPy (for symbolic mathematics) +================================ + + Main authors: `Simon Rohou `_ + +The ``AnalyticFunction`` class already provides interval automatic differentiation through its ``.diff()`` method. +This built-in differentiation computes enclosures of Jacobian values over interval inputs and is the natural tool for guaranteed interval analysis. + +In some situations however, one may need an explicit *symbolic* expression obtained from an existing Codac analytic expression: + +* a simplified expression; +* a symbolic derivative; +* a symbolic gradient, Hessian, or Jacobian; +* a truncated Taylor expansion; +* a polynomial rewritten in Horner form; +* or a symbolic equality test between two ``AnalyticFunction`` objects. + +This is the purpose of the SymPy bridge presented on this page. +It converts a Codac analytic expression into a SymPy expression, applies symbolic transformations in SymPy, and converts the result back into a Codac ``AnalyticFunction``. + +The result remains a first-class Codac object: it can be evaluated, composed, differentiated by intervals, or reused in contractors exactly like any other ``AnalyticFunction``. + +Installing the ``codac-sympy`` extension +---------------------------------------- + +This feature is provided as a Codac extension. + +In C++, using it requires: + +* a Codac build configured with the same Python support as the Python binding, see in particular :ref:`sec-dev-info-binding`; +* and the Python package ``sympy`` available in the Python environment used by the extension. + +In Python, the extension is included by default in Codac. +The remaining requirement is therefore simply to install the Python package ``sympy``: + +.. code-block:: bash + + pip install sympy + +The same package is also required by the C++ extension, since the symbolic backend relies on the Python SymPy runtime. +The Python interpreter is managed internally by the extension in C++, and is of course already available when Codac is used from Python. +No explicit interpreter initialization is needed in user code. + +Finally, in order to include the features of the extension: + +.. tabs:: + + .. group-tab:: Python + + .. literalinclude:: src.py + :language: py + :start-after: [sympy-1-beg] + :end-before: [sympy-1-end] + :dedent: 0 + + .. group-tab:: C++ + + .. literalinclude:: src.cpp + :language: c++ + :start-after: [sympy-1-beg] + :end-before: [sympy-1-end] + :dedent: 0 + + .. code-tab:: matlab + + % todo + + +Simplification +-------------- + +Simplification is often the best first symbolic operation to apply. +It can remove obvious redundancies, expose simpler formulas, and produce expressions that are easier to read, compare, or differentiate. + +.. doxygenfunction:: codac2::sympy_simplify(const AnalyticFunction&) + :project: codac + +A classical trigonometric identity is simplified as expected: + +.. tabs:: + + .. group-tab:: Python + + .. literalinclude:: src.py + :language: py + :start-after: [sympy-2-beg] + :end-before: [sympy-2-end] + :dedent: 4 + + .. group-tab:: C++ + + .. literalinclude:: src.cpp + :language: c++ + :start-after: [sympy-2-beg] + :end-before: [sympy-2-end] + :dedent: 4 + + .. code-tab:: matlab + + % todo + +Simplification is also useful on purely algebraic expressions: + +.. tabs:: + + .. group-tab:: Python + + .. literalinclude:: src.py + :language: py + :start-after: [sympy-3-beg] + :end-before: [sympy-3-end] + :dedent: 4 + + .. group-tab:: C++ + + .. literalinclude:: src.cpp + :language: c++ + :start-after: [sympy-3-beg] + :end-before: [sympy-3-end] + :dedent: 4 + + .. code-tab:: matlab + + % todo + +Partial derivatives and symbolic derivatives +-------------------------------------------- + +The primitive symbolic differentiation operator is ``sympy_partial_diff``. +It returns the symbolic partial derivative of a scalar-valued function with respect to one input variable. + +.. doxygenfunction:: codac2::sympy_partial_diff(const AnalyticFunction&, const ScalarVar&) + :project: codac + +For scalar arguments: + +.. tabs:: + + .. group-tab:: Python + + .. literalinclude:: src.py + :language: py + :start-after: [sympy-4-beg] + :end-before: [sympy-4-end] + :dedent: 4 + + .. group-tab:: C++ + + .. literalinclude:: src.cpp + :language: c++ + :start-after: [sympy-4-beg] + :end-before: [sympy-4-end] + :dedent: 4 + + .. code-tab:: matlab + + % todo + +When a vector argument is involved, a direct component can be used as the differentiation variable. +Note that this notation is usually clearer than referring to a flattened input index, but this type of use is also possible. + +.. tabs:: + + .. group-tab:: Python + + .. literalinclude:: src.py + :language: py + :start-after: [sympy-5-beg] + :end-before: [sympy-5-end] + :dedent: 4 + + .. group-tab:: C++ + + .. literalinclude:: src.cpp + :language: c++ + :start-after: [sympy-5-beg] + :end-before: [sympy-5-end] + :dedent: 4 + + .. code-tab:: matlab + + % todo + +The function ``sympy_diff`` provides symbolic derivatives for scalar functions. +For a univariate scalar function, ``sympy_diff(f)`` returns the first derivative. +Higher-order derivatives are also available. + +.. doxygenfunction:: codac2::sympy_diff(const AnalyticFunction&,const ScalarVar&,Index) + :project: codac + +.. tabs:: + + .. group-tab:: Python + + .. literalinclude:: src.py + :language: py + :start-after: [sympy-6-beg] + :end-before: [sympy-6-end] + :dedent: 4 + + .. group-tab:: C++ + + .. literalinclude:: src.cpp + :language: c++ + :start-after: [sympy-6-beg] + :end-before: [sympy-6-end] + :dedent: 4 + + .. code-tab:: matlab + + % todo + +Gradient, Hessian and Jacobian +------------------------------ + +The symbolic gradient of a scalar function is obtained with ``sympy_gradient``. + +.. doxygenfunction:: codac2::sympy_gradient(const AnalyticFunction&) + :project: codac + +.. tabs:: + + .. group-tab:: Python + + .. literalinclude:: src.py + :language: py + :start-after: [sympy-7-beg] + :end-before: [sympy-7-end] + :dedent: 4 + + .. group-tab:: C++ + + .. literalinclude:: src.cpp + :language: c++ + :start-after: [sympy-7-beg] + :end-before: [sympy-7-end] + :dedent: 4 + + .. code-tab:: matlab + + % todo + +The symbolic Hessian matrix of a scalar function is obtained with ``sympy_hessian``. + +.. doxygenfunction:: codac2::sympy_hessian(const AnalyticFunction&) + :project: codac + +For the scalar function + +.. math:: + + f(x,y) = x y + \sin(x) + y^2, + +its Hessian is + +.. math:: + + H_f(x,y) = + \begin{pmatrix} + -\sin(x) & 1 \\ + 1 & 2 + \end{pmatrix}. + +In Codac: + +.. tabs:: + + .. group-tab:: Python + + .. literalinclude:: src.py + :language: py + :start-after: [sympy-8-beg] + :end-before: [sympy-8-end] + :dedent: 4 + + .. group-tab:: C++ + + .. literalinclude:: src.cpp + :language: c++ + :start-after: [sympy-8-beg] + :end-before: [sympy-8-end] + :dedent: 4 + + .. code-tab:: matlab + + % todo + +For a vector-valued function :math:`\mathbf{f} : \mathbb{R}^n \to \mathbb{R}^m`, symbolic differentiation returns a Jacobian matrix as another ``AnalyticFunction``: + +.. doxygenfunction:: codac2::sympy_diff(const AnalyticFunction&) + :project: codac + +.. math:: + + J_{\mathbf{f}}(\mathbf{x}) = + \begin{pmatrix} + \dfrac{\partial f_1}{\partial x_1} & \cdots & \dfrac{\partial f_1}{\partial x_n} \\ + \vdots & \ddots & \vdots \\ + \dfrac{\partial f_m}{\partial x_1} & \cdots & \dfrac{\partial f_m}{\partial x_n} + \end{pmatrix}. + +.. tabs:: + + .. group-tab:: Python + + .. literalinclude:: src.py + :language: py + :start-after: [sympy-9-beg] + :end-before: [sympy-9-end] + :dedent: 4 + + .. group-tab:: C++ + + .. literalinclude:: src.cpp + :language: c++ + :start-after: [sympy-9-beg] + :end-before: [sympy-9-end] + :dedent: 4 + + .. code-tab:: matlab + + % todo + +Symbolic equality +----------------- + +Symbolic equality is especially useful to validate a transformation, a differentiation result, or a rewriting. + +.. doxygenfunction:: codac2::sympy_equal(const AnalyticFunction&, const AnalyticFunction&) + :project: codac + +For example, the following two scalar functions are symbolically equal although they are expressed with different arguments: + +.. tabs:: + + .. group-tab:: Python + + .. literalinclude:: src.py + :language: py + :start-after: [sympy-10-beg] + :end-before: [sympy-10-end] + :dedent: 4 + + .. group-tab:: C++ + + .. literalinclude:: src.cpp + :language: c++ + :start-after: [sympy-10-beg] + :end-before: [sympy-10-end] + :dedent: 4 + + .. code-tab:: matlab + + % todo + +Series expansions +----------------- + +The bridge can also compute truncated Taylor series for scalar functions. + +.. doxygenfunction:: codac2::sympy_series(const AnalyticFunction&, const ScalarVar&, double, Index) + :project: codac + +For a univariate function: + +.. tabs:: + + .. group-tab:: Python + + .. literalinclude:: src.py + :language: py + :start-after: [sympy-11-beg] + :end-before: [sympy-11-end] + :dedent: 4 + + .. group-tab:: C++ + + .. literalinclude:: src.cpp + :language: c++ + :start-after: [sympy-11-beg] + :end-before: [sympy-11-end] + :dedent: 4 + + .. code-tab:: matlab + + % todo + +A multivariate function can also be expanded with respect to one scalar variable while the other arguments are treated as constants: + +.. tabs:: + + .. group-tab:: Python + + .. literalinclude:: src.py + :language: py + :start-after: [sympy-12-beg] + :end-before: [sympy-12-end] + :dedent: 4 + + .. group-tab:: C++ + + .. literalinclude:: src.cpp + :language: c++ + :start-after: [sympy-12-beg] + :end-before: [sympy-12-end] + :dedent: 4 + + .. code-tab:: matlab + + % todo + +Horner form +----------- + +For polynomial expressions, Horner form is a classical rewriting in interval methods. +It is usually much more favorable for interval evaluations because it reduces repeated occurrences of the variables and therefore helps reduce pessimism. +In practice, this may significantly improve interval polynomial evaluations. + +.. doxygenfunction:: codac2::sympy_horner(const AnalyticFunction&) + :project: codac + +For the polynomial + +.. math:: + + p(x) = 2x^5 + x^3 - 3x^2, + +SymPy can rewrite the expression in Horner form as + +.. math:: + + p(x) = x^2\bigl(-3+x(1+2x^2)\bigr). + +In Codac: + +.. tabs:: + + .. group-tab:: Python + + .. literalinclude:: src.py + :language: py + :start-after: [sympy-13-beg] + :end-before: [sympy-13-end] + :dedent: 4 + + .. group-tab:: C++ + + .. literalinclude:: src.cpp + :language: c++ + :start-after: [sympy-13-beg] + :end-before: [sympy-13-end] + :dedent: 4 + + .. code-tab:: matlab + + % todo + +Relationship with interval automatic differentiation +---------------------------------------------------- + +Symbolic differentiation with SymPy and interval automatic differentiation address two related but different needs. + +Interval automatic differentiation: + +* works directly on interval inputs; +* returns an ``IntervalMatrix`` enclosure of Jacobian values; +* is guaranteed and specifically designed for interval computations. + +Symbolic differentiation with SymPy: + +* works on the symbolic structure of the expression; +* returns another ``AnalyticFunction``; +* is useful when one needs an explicit derivative, Hessian, Jacobian, or rewritten expression that can later be evaluated, composed, or manipulated as a Codac function. + +Symbolic differentiation does not replace interval differentiation. +It complements it. + +Extending the bridge +-------------------- + +The SymPy bridge already covers a practical subset of Codac analytic expressions and symbolic transformations. +It is however intended to grow with user needs. + +If a missing operator, rewriting, simplification, or symbolic transformation would be useful for your applications, please contact the Codac developers. +This extension can be enriched progressively. diff --git a/doc/manual/manual/extensions/sympy/src.cpp b/doc/manual/manual/extensions/sympy/src.cpp new file mode 100644 index 000000000..229b8145e --- /dev/null +++ b/doc/manual/manual/extensions/sympy/src.cpp @@ -0,0 +1,176 @@ +/** + * Codac tests + * ---------------------------------------------------------------------------- + * \date 2024 + * \author Simon Rohou + * \copyright Copyright 2024 Codac Team + * \license GNU Lesser General Public License (LGPL) + */ + +#include +// [sympy-1-beg] +#include +// [sympy-1-end] + +using namespace std; +using namespace codac2; + +TEST_CASE("Sympy extension - manual") +{ + { + // [sympy-2-beg] + ScalarVar x; + AnalyticFunction f({x}, sqr(sin(x)) + sqr(cos(x))); + + auto fs = sympy_simplify(f); + assert(sympy_equal(fs, AnalyticFunction({x}, 1.))); + // [sympy-2-end] + } + + { + // [sympy-3-beg] + ScalarVar x; + AnalyticFunction f({x}, x*(x + 1.) - sqr(x)); + + auto fs = sympy_simplify(f); + assert(sympy_equal(fs, AnalyticFunction({x}, x))); + // [sympy-3-end] + } + + { + // [sympy-4-beg] + ScalarVar x,y; + AnalyticFunction f({x,y}, x*y + sin(x)); + + auto dfdx = sympy_partial_diff(f,x); + auto dfdy = sympy_partial_diff(f,y); + + assert(sympy_equal(dfdx, AnalyticFunction({x,y}, y + cos(x)))); + assert(sympy_equal(dfdy, AnalyticFunction({x,y}, x))); + // [sympy-4-end] + } + + { + // [sympy-5-beg] + VectorVar v(2); + AnalyticFunction g({v}, v[0]*v[1] + sin(v[0])); + + auto dg_dv0 = sympy_partial_diff(g,v[0]); + auto dg_dv1 = sympy_partial_diff(g,v[1]); + + assert(sympy_equal(dg_dv0, AnalyticFunction({v}, v[1] + cos(v[0])))); + assert(sympy_equal(dg_dv1, AnalyticFunction({v}, v[0]))); + // [sympy-5-end] + } + + { + // [sympy-6-beg] + ScalarVar x; + AnalyticFunction f({x}, cos(x)*x + sin(x)); + + auto df = sympy_diff(f); + auto d3f = sympy_diff(f, x, 3); + + assert(sympy_equal(df, AnalyticFunction({x}, 2.*cos(x) - x*sin(x)))); + assert(sympy_equal(d3f, AnalyticFunction({x}, x*sin(x) - 4.*cos(x)))); + // [sympy-6-end] + } + + { + // [sympy-7-beg] + VectorVar v(2); + AnalyticFunction g({v}, v[0]*v[1] + sin(v[0])); + + auto grad_g = sympy_gradient(g); + + assert(sympy_equal( + grad_g, + AnalyticFunction({v}, vec(v[1] + cos(v[0]), v[0])))); + // [sympy-7-end] + } + + { + // [sympy-8-beg] + ScalarVar x,y; + AnalyticFunction f({x,y}, x*y + sin(x) + sqr(y)); + + auto H = sympy_hessian(f); + + assert(sympy_equal( + H, + AnalyticFunction( + {x,y}, + mat( + vec(-sin(x), 1.), + vec(1., 2.) + )))); + // [sympy-8-end] + } + + { + // [sympy-9-beg] + VectorVar v(2); + + AnalyticFunction f({v}, { + v[0]*v[1] + sin(v[0]), + sqr(v[0]) + cos(v[1]) + }); + + auto J = sympy_diff(f); + + assert(sympy_equal( + J, + AnalyticFunction( + {v}, + mat( + vec(v[1] + cos(v[0]), 2.*v[0]), + vec(v[0], -sin(v[1])) + )))); + // [sympy-9-end] + } + + { + // [sympy-10-beg] + ScalarVar x,y; + VectorVar v(2); + + auto f = AnalyticFunction({x,y}, x + 2.*y); + auto g = AnalyticFunction({v}, v[0] + 2.*v[1]); + + assert(sympy_equal(f, g)); + // [sympy-10-end] + } + + { + // [sympy-11-beg] + ScalarVar x; + AnalyticFunction f({x}, 1. / (1. - x)); + + auto p = sympy_series(f, x, 0.0, 3); + + assert(sympy_equal(p, AnalyticFunction({x}, 1. + x + sqr(x) + x*sqr(x)))); + // [sympy-11-end] + } + + { + // [sympy-12-beg] + ScalarVar x,y; + AnalyticFunction f({x,y}, y + 1. / (1. - x)); + + auto p = sympy_series(f, x, 0.0, 2); + + assert(sympy_equal(p, AnalyticFunction({x,y}, y + 1. + x + sqr(x)))); + // [sympy-12-end] + } + + { + // [sympy-13-beg] + ScalarVar x; + AnalyticFunction p({x}, 2.*pow(x,5) + pow(x,3) - 3.*sqr(x)); + + auto h = sympy_horner(p); + + assert(sympy_equal(h, AnalyticFunction({x}, sqr(x)*(-3+x*(1+2*sqr(x)))))); + // [sympy-13-end] + } +} \ No newline at end of file diff --git a/doc/manual/manual/extensions/sympy/src.py b/doc/manual/manual/extensions/sympy/src.py new file mode 100644 index 000000000..3c75bfc05 --- /dev/null +++ b/doc/manual/manual/extensions/sympy/src.py @@ -0,0 +1,159 @@ +#!/usr/bin/env python + +# Codac tests +# ---------------------------------------------------------------------------- +# \date 2026 +# \author Simon Rohou +# \copyright Copyright 2026 Codac Team +# \license GNU Lesser General Public License (LGPL) + +import unittest +# [sympy-1-beg] +from codac import * +# sympy extension is already embedded in Codac +# [sympy-1-end] + + +class TestSympyManual(unittest.TestCase): + + def test_sympy_manual(self): + # [sympy-2-beg] + x = ScalarVar() + f = AnalyticFunction([x], sqr(sin(x)) + sqr(cos(x))) + + fs = sympy_simplify(f) + assert sympy_equal(fs, AnalyticFunction([x], 1.)) + # [sympy-2-end] + + # [sympy-3-beg] + x = ScalarVar() + f = AnalyticFunction([x], x*(x + 1.) - sqr(x)) + + fs = sympy_simplify(f) + assert sympy_equal(fs, AnalyticFunction([x], x)) + # [sympy-3-end] + + # [sympy-4-beg] + x = ScalarVar() + y = ScalarVar() + f = AnalyticFunction([x,y], x*y + sin(x)) + + dfdx = sympy_partial_diff(f, x) + dfdy = sympy_partial_diff(f, y) + + assert sympy_equal(dfdx, AnalyticFunction([x,y], y + cos(x))) + assert sympy_equal(dfdy, AnalyticFunction([x,y], x)) + # [sympy-4-end] + + # [sympy-5-beg] + v = VectorVar(2) + g = AnalyticFunction([v], v[0]*v[1] + sin(v[0])) + + dg_dv0 = sympy_partial_diff(g, v[0]) + dg_dv1 = sympy_partial_diff(g, v[1]) + + assert sympy_equal(dg_dv0, AnalyticFunction([v], v[1] + cos(v[0]))) + assert sympy_equal(dg_dv1, AnalyticFunction([v], v[0])) + # [sympy-5-end] + + # [sympy-6-beg] + x = ScalarVar() + f = AnalyticFunction([x], cos(x)*x + sin(x)) + + df = sympy_diff(f) + d3f = sympy_diff(f, x, 3) + + assert sympy_equal(df, AnalyticFunction([x], 2.*cos(x) - x*sin(x))) + assert sympy_equal(d3f, AnalyticFunction([x], x*sin(x) - 4.*cos(x))) + # [sympy-6-end] + + # [sympy-7-beg] + v = VectorVar(2) + g = AnalyticFunction([v], v[0]*v[1] + sin(v[0])) + + grad_g = sympy_gradient(g) + + assert sympy_equal( + grad_g, + AnalyticFunction([v], vec(v[1] + cos(v[0]), v[0]))) + # [sympy-7-end] + + # [sympy-8-beg] + x = ScalarVar() + y = ScalarVar() + f = AnalyticFunction([x,y], x*y + sin(x) + sqr(y)) + + H = sympy_hessian(f) + + assert sympy_equal( + H, + AnalyticFunction( + [x,y], + mat( + vec(-sin(x), 1.), + vec(1., 2.) + ))) + # [sympy-8-end] + + # [sympy-9-beg] + v = VectorVar(2) + + f = AnalyticFunction([v], [ + v[0]*v[1] + sin(v[0]), + sqr(v[0]) + cos(v[1]) + ]) + + J = sympy_diff(f) + + assert sympy_equal( + J, + AnalyticFunction( + [v], + mat( + vec(v[1] + cos(v[0]), 2.*v[0]), + vec(v[0], -sin(v[1])) + ))) + # [sympy-9-end] + + # [sympy-10-beg] + x = ScalarVar() + y = ScalarVar() + v = VectorVar(2) + + f = AnalyticFunction([x,y], x + 2.*y) + g = AnalyticFunction([v], v[0] + 2.*v[1]) + + assert sympy_equal(f, g) + # [sympy-10-end] + + # [sympy-11-beg] + x = ScalarVar() + f = AnalyticFunction([x], 1. / (1. - x)) + + p = sympy_series(f, x, 0.0, 3) + + assert sympy_equal(p, AnalyticFunction([x], 1. + x + sqr(x) + x*sqr(x))) + # [sympy-11-end] + + # [sympy-12-beg] + x = ScalarVar() + y = ScalarVar() + f = AnalyticFunction([x,y], y + 1. / (1. - x)) + + p = sympy_series(f, x, 0.0, 2) + + assert sympy_equal(p, AnalyticFunction([x,y], y + 1. + x + sqr(x))) + # [sympy-12-end] + + # [sympy-13-beg] + x = ScalarVar() + p = AnalyticFunction([x], 2.*pow(x,5) + pow(x,3) - 3.*sqr(x)) + + h = sympy_horner(p) + + assert sympy_equal(h, AnalyticFunction([x], sqr(x)*(-3 + x*(1 + 2*sqr(x))))) + # [sympy-13-end] + + +if __name__ == '__main__': + unittest.main() diff --git a/doc/manual/manual/functions/analytic/analytic_functions.rst b/doc/manual/manual/functions/analytic/analytic_functions.rst index 1752148f7..60fb121fd 100644 --- a/doc/manual/manual/functions/analytic/analytic_functions.rst +++ b/doc/manual/manual/functions/analytic/analytic_functions.rst @@ -317,6 +317,18 @@ The ``.diff()`` method can be used in the same way as the ``.eval()`` method. :dedent: 0 +Symbolic manipulations with SymPy +--------------------------------- + +Besides interval automatic differentiation, Codac also provides an optional SymPy-based extension for symbolic manipulations of ``AnalyticFunction`` objects. + +This extension makes it possible to simplify expressions symbolically, compute symbolic partial derivatives, gradients, Hessians and Jacobians, derive truncated Taylor series, rewrite polynomials in Horner form, and test symbolic equality between analytic expressions. The result of these operations is still represented as a Codac ``AnalyticFunction``, which means that the transformed expressions can then be evaluated, composed, or used in interval computations exactly like any other analytic function. + +This symbolic extension complements interval automatic differentiation rather than replacing it. Interval differentiation directly returns guaranteed enclosures of derivative values over boxes, whereas the SymPy extension returns explicit symbolic expressions that may then be manipulated or evaluated afterwards. + +See the dedicated page :ref:`sec-extensions-sympy`. + + Other properties ---------------- @@ -350,5 +362,4 @@ Let us consider a function :math:`[\mathbf{f}]:\mathbb{IR}^n\to\mathbb{IR}^m`, t In the case of multivariate functions, ``.input_size()`` returns the sum of the dimensions of the arguments. - -The ``AnalyticFunction`` class supports many mathematical operations, and the full set of operators that can be used is described in the next page. \ No newline at end of file +The ``AnalyticFunction`` class supports many mathematical operations. The full set of analytic operators is described in the next page, while symbolic manipulations based on SymPy are presented in :ref:`sec-extensions-sympy`. \ No newline at end of file diff --git a/doc/manual/manual/functions/analytic/src.cpp b/doc/manual/manual/functions/analytic/src.cpp index 5f48e424d..4eaad4549 100644 --- a/doc/manual/manual/functions/analytic/src.cpp +++ b/doc/manual/manual/functions/analytic/src.cpp @@ -19,7 +19,7 @@ TEST_CASE("AnalyticFunction - manual") { { // [1-beg] - ScalarVar x1, x2; + ScalarVar x1, x2("x2"); // last variable has a specific a name VectorVar v(3); // Example of scalar function: from R to R diff --git a/doc/manual/manual/functions/analytic/src.m b/doc/manual/manual/functions/analytic/src.m index 39740a987..d61657e08 100644 --- a/doc/manual/manual/functions/analytic/src.m +++ b/doc/manual/manual/functions/analytic/src.m @@ -2,7 +2,7 @@ % [1-beg] x1 = ScalarVar(); -x2 = ScalarVar(); +x2 = ScalarVar("x2"); % a variable with a name v = VectorVar(3); % Example of scalar function: from R to R diff --git a/doc/manual/manual/functions/analytic/src.py b/doc/manual/manual/functions/analytic/src.py index eb2df3ec5..886d7f59f 100644 --- a/doc/manual/manual/functions/analytic/src.py +++ b/doc/manual/manual/functions/analytic/src.py @@ -18,7 +18,7 @@ def tests_AnalyticFunction_manual(test): # [1-beg] x1 = ScalarVar() - x2 = ScalarVar() + x2 = ScalarVar("x2") # a variable with a name v = VectorVar(3) # Example of scalar function: from R to R diff --git a/doc/manual/manual/intervals/src.cpp b/doc/manual/manual/intervals/src.cpp index ca3508748..8b9ebca36 100644 --- a/doc/manual/manual/intervals/src.cpp +++ b/doc/manual/manual/intervals/src.cpp @@ -16,21 +16,18 @@ using namespace codac2; TEST_CASE("BoolInterval class - manual") { - { + /* // [boolinterval-class-1-beg] BoolInterval::FALSE; // certainly false BoolInterval::TRUE; // certainly true BoolInterval::UNKNOWN; // undetermined BoolInterval::EMPTY; // inconsistent / impossible // [boolinterval-class-1-end] - } - /* - { + // [boolinterval-class-2-beg] BoolInterval::UNKNOWN == BoolInterval::TRUE | BoolInterval::FALSE BoolInterval::EMPTY == BoolInterval::TRUE & BoolInterval::FALSE // [boolinterval-class-2-end] - } */ } diff --git a/python/codac/__init__.py b/python/codac/__init__.py index 853f37ac0..66a8d1735 100644 --- a/python/codac/__init__.py +++ b/python/codac/__init__.py @@ -1,4 +1,5 @@ from codac.core import * from codac.graphics import * +from codac.sympy import * from codac.unsupported import * from .version import __version__ \ No newline at end of file diff --git a/python/codac/core/__init__.py b/python/codac/core/__init__.py index 441128035..c3369ec99 100644 --- a/python/codac/core/__init__.py +++ b/python/codac/core/__init__.py @@ -17,7 +17,7 @@ def codac_error(message): class AnalyticFunction: def __init__(self, args, e=None): - if e: + if e is not None: if isinstance(e, (int,float,Interval,ScalarVar,ScalarExpr)): self.f = AnalyticFunction_Scalar(args,ScalarExpr(e)) elif isinstance(e, (Vector,IntervalVector,VectorVar,VectorExpr)): diff --git a/python/codac/sympy/__init__.py b/python/codac/sympy/__init__.py new file mode 100644 index 000000000..915a51847 --- /dev/null +++ b/python/codac/sympy/__init__.py @@ -0,0 +1,70 @@ +from codac._core._sympy import * +from codac import AnalyticFunction, AnalyticFunction_Scalar, AnalyticFunction_Vector, AnalyticFunction_Matrix + +def sympy_simplify(f): + if isinstance(f.f, (AnalyticFunction_Scalar)): + return AnalyticFunction(sympy_simplify_scalar(f.f)) + elif isinstance(f.f, (AnalyticFunction_Vector)): + return AnalyticFunction(sympy_simplify_vector(f.f)) + elif isinstance(f.f, (AnalyticFunction_Matrix)): + return AnalyticFunction(sympy_simplify_matrix(f.f)) + else: + codac_error("sympy_simplify: invalid function input") + +def sympy_horner(f): + if isinstance(f.f, (AnalyticFunction_Scalar)): + return AnalyticFunction(sympy_horner_scalar(f.f)) + elif isinstance(f.f, (AnalyticFunction_Vector)): + return AnalyticFunction(sympy_horner_vector(f.f)) + elif isinstance(f.f, (AnalyticFunction_Matrix)): + return AnalyticFunction(sympy_horner_matrix(f.f)) + else: + codac_error("sympy_horner: invalid function input") + +def sympy_partial_diff(f,x): + if isinstance(f.f, (AnalyticFunction_Scalar)): + return AnalyticFunction(sympy_partial_diff_(f.f,x)) + else: + codac_error("sympy_partial_diff: invalid function input") + +def sympy_diff(f,*args): + if isinstance(f.f, (AnalyticFunction_Scalar)): + return AnalyticFunction(sympy_diff_scalar(f.f,*args)) + elif isinstance(f.f, (AnalyticFunction_Vector)): + return AnalyticFunction(sympy_diff_vector(f.f,*args)) + else: + codac_error("sympy_diff: invalid function input") + +def sympy_gradient(f): + if isinstance(f.f, (AnalyticFunction_Scalar)): + return AnalyticFunction(sympy_gradient_(f.f)) + else: + codac_error("sympy_gradient: invalid function input") + +def sympy_hessian(f): + if isinstance(f.f, (AnalyticFunction_Scalar)): + return AnalyticFunction(sympy_hessian_(f.f)) + else: + codac_error("sympy_hessian: invalid function input") + +def sympy_series(f,*args): + if isinstance(f.f, (AnalyticFunction_Scalar)): + return AnalyticFunction(sympy_series_(f.f,*args)) + else: + codac_error("sympy_series: invalid function input") + +def sympy_equal(f,g): + if isinstance(f.f, (AnalyticFunction_Scalar)): + if not isinstance(g.f, (AnalyticFunction_Scalar)): + return False + return sympy_equal_scalar(f.f,g.f) + elif isinstance(f.f, (AnalyticFunction_Vector)): + if not isinstance(g.f, (AnalyticFunction_Vector)): + return False + return sympy_equal_vector(f.f,g.f) + elif isinstance(f.f, (AnalyticFunction_Matrix)): + if not isinstance(g.f, (AnalyticFunction_Matrix)): + return False + return sympy_equal_matrix(f.f,g.f) + else: + codac_error("sympy_equal: invalid function inputs") \ No newline at end of file diff --git a/python/setup.py.in b/python/setup.py.in index 6262e4b14..633f927e3 100644 --- a/python/setup.py.in +++ b/python/setup.py.in @@ -24,6 +24,7 @@ setup( '${PYTHON_PACKAGE_NAME}', '${PYTHON_PACKAGE_NAME}.core', '${PYTHON_PACKAGE_NAME}.graphics', + '${PYTHON_PACKAGE_NAME}.sympy', '${PYTHON_PACKAGE_NAME}.unsupported', '${PYTHON_PACKAGE_NAME}.tests' ], @@ -31,6 +32,7 @@ setup( '${PYTHON_PACKAGE_NAME}': [ '_core${PYTHON_MODULE_EXTENSION}', '_graphics${PYTHON_MODULE_EXTENSION}', + '_sympy${PYTHON_MODULE_EXTENSION}', '_unsupported${PYTHON_MODULE_EXTENSION}', ], }, diff --git a/python/src/core/CMakeLists.txt b/python/src/core/CMakeLists.txt index f29275756..c3386deaa 100644 --- a/python/src/core/CMakeLists.txt +++ b/python/src/core/CMakeLists.txt @@ -62,6 +62,8 @@ domains/tube/codac2_py_tube_cart_prod.cpp domains/codac2_py_BoolInterval.cpp + extensions/sympy/codac2_py_sympy.cpp + functions/analytic/codac2_py_analytic_variables.cpp functions/analytic/codac2_py_AnalyticFunction.h functions/analytic/codac2_py_AnalyticExprWrapper.h @@ -136,6 +138,7 @@ PRIVATE ${CMAKE_CURRENT_SOURCE_DIR}/domains/paving/ PRIVATE ${CMAKE_CURRENT_SOURCE_DIR}/domains/zonotope/ PRIVATE ${CMAKE_CURRENT_SOURCE_DIR}/domains/tube/ + PRIVATE ${CMAKE_CURRENT_SOURCE_DIR}/extensions/sympy/ PRIVATE ${CMAKE_CURRENT_SOURCE_DIR}/functions/ PRIVATE ${CMAKE_CURRENT_SOURCE_DIR}/functions/analytic/ PRIVATE ${CMAKE_CURRENT_SOURCE_DIR}/graphics/ @@ -148,7 +151,7 @@ ) target_link_libraries(_core - PRIVATE ${PROJECT_NAME}-core ${LIBS} Ibex::ibex + PRIVATE ${PROJECT_NAME}-core ${PROJECT_NAME}-sympy ${LIBS} Ibex::ibex ) # Copy the generated library in the package folder diff --git a/python/src/core/codac2_py_core.cpp b/python/src/core/codac2_py_core.cpp index 79c1adb82..582ff27e2 100644 --- a/python/src/core/codac2_py_core.cpp +++ b/python/src/core/codac2_py_core.cpp @@ -159,6 +159,9 @@ void export_trunc(py::module& m); void export_AnalyticTraj(py::module& m); void export_SampledTraj(py::module& m); +// Extension > sympy +void export_sympy(py::module& m); + PYBIND11_MODULE(_core, m) { @@ -325,7 +328,6 @@ PYBIND11_MODULE(_core, m) export_AnalyticTraj(m); export_SampledTraj(m); - m.def("srand", []() { srand(time(NULL)); @@ -338,4 +340,8 @@ PYBIND11_MODULE(_core, m) }, DOC_TO_BE_DEFINED, "seed"_a); -} + + // Extension > sympy + auto ms = m.def_submodule("_sympy"); // to keep a dedicated namespace + export_sympy(ms); +} \ No newline at end of file diff --git a/python/src/core/extensions/sympy/codac2_py_sympy.cpp b/python/src/core/extensions/sympy/codac2_py_sympy.cpp new file mode 100644 index 000000000..b86f2102c --- /dev/null +++ b/python/src/core/extensions/sympy/codac2_py_sympy.cpp @@ -0,0 +1,131 @@ +/** + * \file + * Codac binding (sympy) + * ---------------------------------------------------------------------------- + * \date 2026 + * \author Simon Rohou + * \copyright Copyright 2026 Codac Team + * \license GNU Lesser General Public License (LGPL) + */ + +#include +#include +#include +#include +#include "codac2_py_sympy_docs.h" // Generated file from Doxygen XML (doxygen2docstring.py): +#include "codac2_py_cast.h" + +using namespace codac2; +namespace py = pybind11; + + +void export_sympy(py::module& m) +{ + // sympy_simplify + + m.def("sympy_simplify_scalar", (AnalyticFunction (*)(const AnalyticFunction&))&codac2::sympy_simplify, + ANALYTICFUNCTION_SCALARTYPE_SYMPY_SIMPLIFY_CONST_ANALYTICFUNCTION_SCALARTYPE_REF, + "f"_a); + + m.def("sympy_simplify_vector", (AnalyticFunction (*)(const AnalyticFunction&))&codac2::sympy_simplify, + ANALYTICFUNCTION_VECTORTYPE_SYMPY_SIMPLIFY_CONST_ANALYTICFUNCTION_VECTORTYPE_REF, + "f"_a); + + m.def("sympy_simplify_matrix", (AnalyticFunction (*)(const AnalyticFunction&))&codac2::sympy_simplify, + ANALYTICFUNCTION_MATRIXTYPE_SYMPY_SIMPLIFY_CONST_ANALYTICFUNCTION_MATRIXTYPE_REF, + "f"_a); + + // sympy_horner + + m.def("sympy_horner_scalar", (AnalyticFunction (*)(const AnalyticFunction&))&codac2::sympy_horner, + ANALYTICFUNCTION_SCALARTYPE_SYMPY_HORNER_CONST_ANALYTICFUNCTION_SCALARTYPE_REF, + "f"_a); + + m.def("sympy_horner_vector", (AnalyticFunction (*)(const AnalyticFunction&))&codac2::sympy_horner, + ANALYTICFUNCTION_VECTORTYPE_SYMPY_HORNER_CONST_ANALYTICFUNCTION_VECTORTYPE_REF, + "f"_a); + + m.def("sympy_horner_matrix", (AnalyticFunction (*)(const AnalyticFunction&))&codac2::sympy_horner, + ANALYTICFUNCTION_MATRIXTYPE_SYMPY_HORNER_CONST_ANALYTICFUNCTION_MATRIXTYPE_REF, + "f"_a); + + // sympy_partial_diff + + m.def("sympy_partial_diff_", (AnalyticFunction (*)(const AnalyticFunction&,const ScalarVar&))&codac2::sympy_partial_diff, + ANALYTICFUNCTION_SCALARTYPE_SYMPY_PARTIAL_DIFF_CONST_ANALYTICFUNCTION_SCALARTYPE_REF_CONST_SCALARVAR_REF, + "f"_a, "x"_a); + + m.def("sympy_partial_diff_", (AnalyticFunction (*)(const AnalyticFunction&,const ScalarExpr&))&codac2::sympy_partial_diff, + ANALYTICFUNCTION_SCALARTYPE_SYMPY_PARTIAL_DIFF_CONST_ANALYTICFUNCTION_SCALARTYPE_REF_CONST_SCALARVAR_REF, + "f"_a, "x"_a); + + // sympy_diff + + m.def("sympy_diff_scalar", (AnalyticFunction (*)(const AnalyticFunction&))&codac2::sympy_diff, + ANALYTICFUNCTION_SCALARTYPE_SYMPY_DIFF_CONST_ANALYTICFUNCTION_SCALARTYPE_REF, + "f"_a); + + m.def("sympy_diff_scalar", (AnalyticFunction (*)(const AnalyticFunction&,const ScalarVar&))&codac2::sympy_diff, + ANALYTICFUNCTION_SCALARTYPE_SYMPY_DIFF_CONST_ANALYTICFUNCTION_SCALARTYPE_REF_CONST_SCALARVAR_REF, + "f"_a, "x"_a); + + m.def("sympy_diff_scalar", (AnalyticFunction (*)(const AnalyticFunction&,const ScalarExpr&))&codac2::sympy_diff, + ANALYTICFUNCTION_SCALARTYPE_SYMPY_DIFF_CONST_ANALYTICFUNCTION_SCALARTYPE_REF_CONST_SCALARVAR_REF, + "f"_a, "x"_a); + + m.def("sympy_diff_scalar", (AnalyticFunction (*)(const AnalyticFunction&,Index))&codac2::sympy_diff, + ANALYTICFUNCTION_SCALARTYPE_SYMPY_DIFF_CONST_ANALYTICFUNCTION_SCALARTYPE_REF_INDEX, + "f"_a, "order"_a); + + m.def("sympy_diff_scalar", (AnalyticFunction (*)(const AnalyticFunction&,const ScalarVar&,Index))&codac2::sympy_diff, + ANALYTICFUNCTION_SCALARTYPE_SYMPY_DIFF_CONST_ANALYTICFUNCTION_SCALARTYPE_REF_CONST_SCALARVAR_REF_INDEX, + "f"_a, "x"_a, "order"_a); + + m.def("sympy_diff_scalar", (AnalyticFunction (*)(const AnalyticFunction&,const ScalarExpr&,Index))&codac2::sympy_diff, + ANALYTICFUNCTION_SCALARTYPE_SYMPY_DIFF_CONST_ANALYTICFUNCTION_SCALARTYPE_REF_CONST_SCALARVAR_REF_INDEX, + "f"_a, "x"_a, "order"_a); + + m.def("sympy_diff_vector", (AnalyticFunction (*)(const AnalyticFunction&))&codac2::sympy_diff, + ANALYTICFUNCTION_MATRIXTYPE_SYMPY_DIFF_CONST_ANALYTICFUNCTION_VECTORTYPE_REF, + "f"_a); + + // sympy_gradient + + m.def("sympy_gradient_", (AnalyticFunction (*)(const AnalyticFunction&))&codac2::sympy_gradient, + ANALYTICFUNCTION_VECTORTYPE_SYMPY_GRADIENT_CONST_ANALYTICFUNCTION_SCALARTYPE_REF, + "f"_a); + + // sympy_hessian + + m.def("sympy_hessian_", (AnalyticFunction (*)(const AnalyticFunction&))&codac2::sympy_hessian, + ANALYTICFUNCTION_MATRIXTYPE_SYMPY_HESSIAN_CONST_ANALYTICFUNCTION_SCALARTYPE_REF, + "f"_a); + + // sympy_series + + m.def("sympy_series_", (AnalyticFunction (*)(const AnalyticFunction&,double,Index))&codac2::sympy_series, + ANALYTICFUNCTION_SCALARTYPE_SYMPY_SERIES_CONST_ANALYTICFUNCTION_SCALARTYPE_REF_DOUBLE_INDEX, + "f"_a, "center"_a, "order"_a); + + m.def("sympy_series_", (AnalyticFunction (*)(const AnalyticFunction&,const ScalarVar&,double,Index))&codac2::sympy_series, + ANALYTICFUNCTION_SCALARTYPE_SYMPY_SERIES_CONST_ANALYTICFUNCTION_SCALARTYPE_REF_CONST_SCALARVAR_REF_DOUBLE_INDEX, + "f"_a, "x"_a, "center"_a, "order"_a); + + m.def("sympy_series_", (AnalyticFunction (*)(const AnalyticFunction&,const ScalarExpr&,double,Index))&codac2::sympy_series, + ANALYTICFUNCTION_SCALARTYPE_SYMPY_SERIES_CONST_ANALYTICFUNCTION_SCALARTYPE_REF_CONST_SCALARVAR_REF_DOUBLE_INDEX, + "f"_a, "x"_a, "center"_a, "order"_a); + + // sympy_equal + + m.def("sympy_equal_scalar", (bool (*)(const AnalyticFunction&,const AnalyticFunction&))&codac2::sympy_equal, + BOOL_SYMPY_EQUAL_CONST_ANALYTICFUNCTION_SCALARTYPE_REF_CONST_ANALYTICFUNCTION_SCALARTYPE_REF, + "f"_a, "g"_a); + + m.def("sympy_equal_vector", (bool (*)(const AnalyticFunction&,const AnalyticFunction&))&codac2::sympy_equal, + BOOL_SYMPY_EQUAL_CONST_ANALYTICFUNCTION_VECTORTYPE_REF_CONST_ANALYTICFUNCTION_VECTORTYPE_REF, + "f"_a, "g"_a); + + m.def("sympy_equal_matrix", (bool (*)(const AnalyticFunction&,const AnalyticFunction&))&codac2::sympy_equal, + BOOL_SYMPY_EQUAL_CONST_ANALYTICFUNCTION_MATRIXTYPE_REF_CONST_ANALYTICFUNCTION_MATRIXTYPE_REF, + "f"_a, "g"_a); +} \ No newline at end of file diff --git a/scripts/docker/build_pybinding.sh b/scripts/docker/build_pybinding.sh index 7950be548..048640945 100755 --- a/scripts/docker/build_pybinding.sh +++ b/scripts/docker/build_pybinding.sh @@ -30,6 +30,7 @@ for PYBIN in /opt/python/cp3*/bin; do "${PYBIN}/python" -m pip install codac --no-deps --no-index -f /io/wheelhouse "${PYBIN}/python" ../examples/02_centered_form/main.py "${PYBIN}/python" -m pip install numpy --prefer-binary + "${PYBIN}/python" -m pip install sympy "${PYBIN}/python" -m unittest discover codac.tests make test ARGS="-V --output-on-failure" diff --git a/scripts/docker/build_pybinding_codac4matlab.sh b/scripts/docker/build_pybinding_codac4matlab.sh index 289feebb3..e9a688e88 100644 --- a/scripts/docker/build_pybinding_codac4matlab.sh +++ b/scripts/docker/build_pybinding_codac4matlab.sh @@ -31,6 +31,7 @@ for PYBIN in /opt/python/cp3*/bin; do "${PYBIN}/python" -c "import sys; print(sys.version); import codac4matlab; print(codac4matlab.__version__); from codac4matlab import *; print(IntervalVector([[-0.1],[0],[0.2]]))" #"${PYBIN}/python" ../examples/02_centered_form/main.py #"${PYBIN}/python" -m pip install numpy --prefer-binary + #"${PYBIN}/python" -m pip install sympy #"${PYBIN}/python" -m unittest discover codac.tests #make test ARGS="-V --output-on-failure" diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 34511403e..9843f864e 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -14,9 +14,9 @@ add_subdirectory(extensions/capd) endif() - #if(WITH_PYTHON) - # add_subdirectory(sympy) - #endif() + if(WITH_PYTHON) + add_subdirectory(extensions/sympy) + endif() # Generating PKG file @@ -95,6 +95,24 @@ # #endif() + if(WITH_PYTHON) + + file(APPEND ${CODAC_CMAKE_CONFIG_FILE} " + + # Optional 3rd party: + + find_path(CODAC_SYMPY_INCLUDE_DIR ${PROJECT_NAME}-sympy.h + PATH_SUFFIXES include/${PROJECT_NAME}-sympy) + set(CODAC_INCLUDE_DIRS \${CODAC_INCLUDE_DIRS} \${CODAC_SYMPY_INCLUDE_DIR}) + + find_library(CODAC_SYMPY_LIBRARY NAMES ${PROJECT_NAME}-sympy + PATH_SUFFIXES lib) + + set(CODAC_LIBRARIES \${CODAC_LIBRARIES} \${CODAC_SYMPY_LIBRARY}) + ") + + endif() + if(WITH_CAPD) file(APPEND ${CODAC_CMAKE_CONFIG_FILE} " diff --git a/src/core/CMakeLists.txt b/src/core/CMakeLists.txt index ffec45e23..f7fb12b23 100644 --- a/src/core/CMakeLists.txt +++ b/src/core/CMakeLists.txt @@ -99,8 +99,12 @@ ${CMAKE_CURRENT_SOURCE_DIR}/functions/codac2_FunctionArgsList.h ${CMAKE_CURRENT_SOURCE_DIR}/functions/codac2_FunctionBase.h ${CMAKE_CURRENT_SOURCE_DIR}/functions/codac2_VarBase.h + ${CMAKE_CURRENT_SOURCE_DIR}/functions/analytic/codac2_analytic_componentwise.cpp + ${CMAKE_CURRENT_SOURCE_DIR}/functions/analytic/codac2_analytic_componentwise.h ${CMAKE_CURRENT_SOURCE_DIR}/functions/analytic/codac2_analytic_constants.h ${CMAKE_CURRENT_SOURCE_DIR}/functions/analytic/codac2_analytic_constants_impl.h + ${CMAKE_CURRENT_SOURCE_DIR}/functions/analytic/codac2_analytic_flat_input_layout.cpp + ${CMAKE_CURRENT_SOURCE_DIR}/functions/analytic/codac2_analytic_flat_input_layout.h ${CMAKE_CURRENT_SOURCE_DIR}/functions/analytic/codac2_analytic_variables.cpp ${CMAKE_CURRENT_SOURCE_DIR}/functions/analytic/codac2_analytic_variables.h ${CMAKE_CURRENT_SOURCE_DIR}/functions/analytic/codac2_AnalyticExpr.h diff --git a/src/core/domains/tube/codac2_SlicedTube.h b/src/core/domains/tube/codac2_SlicedTube.h index 64bb6f2e9..4c771bc26 100644 --- a/src/core/domains/tube/codac2_SlicedTube.h +++ b/src/core/domains/tube/codac2_SlicedTube.h @@ -392,7 +392,7 @@ namespace codac2 friend inline std::ostream& operator<<(std::ostream& os, const SlicedTube& x) { os << x.t0_tf() - << "↦" << x.codomain() + << "↦" << (x.is_empty() ? x.empty_value() : x.codomain()) << ", " << x.nb_slices() << " slice" << (x.nb_slices() > 1 ? "s" : "") << std::flush; @@ -591,6 +591,19 @@ namespace codac2 return x; } + template + // requires std::is_same_v::Domain,T> + inline auto mid() const + { + SampledTraj m; + double t0 = _tdomain->t0_tf().lb(); + m.set((*this)(t0).mid(), t0); + for(const auto& s : *this) + if(!s.is_gate()) + m.set(s.output_gate().mid(),s.t0_tf().ub()); + return m; + } + public: diff --git a/src/core/functions/analytic/codac2_AnalyticExpr.h b/src/core/functions/analytic/codac2_AnalyticExpr.h index cc27110b2..fba88db74 100644 --- a/src/core/functions/analytic/codac2_AnalyticExpr.h +++ b/src/core/functions/analytic/codac2_AnalyticExpr.h @@ -151,5 +151,10 @@ namespace codac2 return b; } + + std::vector> children_expr_base() const override + { + return OperationExprBase...>::children_expr_base(); + } }; } diff --git a/src/core/functions/analytic/codac2_AnalyticExprWrapper.h b/src/core/functions/analytic/codac2_AnalyticExprWrapper.h index dce877121..451dca0a3 100644 --- a/src/core/functions/analytic/codac2_AnalyticExprWrapper.h +++ b/src/core/functions/analytic/codac2_AnalyticExprWrapper.h @@ -31,6 +31,9 @@ namespace codac2 AnalyticExprWrapper(const std::shared_ptr>& e) : std::shared_ptr>(e) { } + + AnalyticExprWrapper& operator=(const AnalyticExprWrapper&) = default; + AnalyticExprWrapper& operator=(AnalyticExprWrapper&&) noexcept = default; template requires std::is_base_of_v,V> diff --git a/src/core/functions/analytic/codac2_analytic_componentwise.cpp b/src/core/functions/analytic/codac2_analytic_componentwise.cpp new file mode 100644 index 000000000..77295d408 --- /dev/null +++ b/src/core/functions/analytic/codac2_analytic_componentwise.cpp @@ -0,0 +1,62 @@ +/** + * \file codac2_analytic_componentwise.cpp + * ---------------------------------------------------------------------------- + * \date 2026 + * \author Simon Rohou + * \copyright Copyright 2026 Codac Team + * \license GNU Lesser General Public License (LGPL) + */ + +#include +#include "codac2_assert.h" +#include "codac2_analytic_componentwise.h" + +namespace codac2 +{ + AnalyticFunction + map_scalar_entries(const AnalyticFunction& f, const ScalarExprTransform& transform) + { + return AnalyticFunction(f.args(), transform(ScalarExpr(f.expr()))); + } + + AnalyticFunction + map_scalar_entries(const AnalyticFunction& f, const ScalarExprTransform& transform) + { + const auto shape = f.output_shape(); + + assert(shape.second == 1 + && "map_scalar_entries(VectorType): only column-vector outputs are supported"); + + VectorExpr y(f.expr()); + std::vector entries; + entries.reserve(static_cast(shape.first)); + + for(Index i = 0 ; i < shape.first ; ++i) + entries.push_back(transform(y[i])); + + return AnalyticFunction(f.args(), vec(entries)); + } + + AnalyticFunction + map_scalar_entries(const AnalyticFunction& f, const ScalarExprTransform& transform) + { + const auto shape = f.output_shape(); + + MatrixExpr y(f.expr()); + std::vector cols; + cols.reserve(static_cast(shape.second)); + + for(Index j = 0 ; j < shape.second ; ++j) + { + std::vector col_entries; + col_entries.reserve(static_cast(shape.first)); + + for(Index i = 0 ; i < shape.first ; ++i) + col_entries.push_back(transform(y(i,j))); + + cols.push_back(vec(col_entries)); + } + + return AnalyticFunction(f.args(), mat(cols)); + } +} \ No newline at end of file diff --git a/src/core/functions/analytic/codac2_analytic_componentwise.h b/src/core/functions/analytic/codac2_analytic_componentwise.h new file mode 100644 index 000000000..0f785c6c1 --- /dev/null +++ b/src/core/functions/analytic/codac2_analytic_componentwise.h @@ -0,0 +1,59 @@ +/** + * \file codac2_analytic_componentwise.h + * ---------------------------------------------------------------------------- + * \date 2026 + * \author Simon Rohou + * \copyright Copyright 2026 Codac Team + * \license GNU Lesser General Public License (LGPL) + */ + +#pragma once + +#include + +#include "codac2_AnalyticFunction.h" + +namespace codac2 +{ + /** + * \brief Scalar transformation applied to the entries of an analytic function output. + * + * This callback takes one scalar analytic expression and returns the + * transformed scalar expression. + */ + using ScalarExprTransform = std::function; + + /** + * \brief Applies a scalar transformation to a scalar analytic function. + * + * \param f Scalar function. + * \param transform Scalar transformation. + * \return The transformed scalar function. + */ + AnalyticFunction + map_scalar_entries(const AnalyticFunction& f, const ScalarExprTransform& transform); + + /** + * \brief Applies a scalar transformation componentwise to a vector analytic function. + * + * The output shape of \p f is preserved. + * + * \param f Vector function. + * \param transform Scalar transformation. + * \return The transformed vector function. + */ + AnalyticFunction + map_scalar_entries(const AnalyticFunction& f, const ScalarExprTransform& transform); + + /** + * \brief Applies a scalar transformation entrywise to a matrix analytic function. + * + * The output shape of \p f is preserved. + * + * \param f Matrix function. + * \param transform Scalar transformation. + * \return The transformed matrix function. + */ + AnalyticFunction + map_scalar_entries(const AnalyticFunction& f, const ScalarExprTransform& transform); +} \ No newline at end of file diff --git a/src/core/functions/analytic/codac2_analytic_flat_input_layout.cpp b/src/core/functions/analytic/codac2_analytic_flat_input_layout.cpp new file mode 100644 index 000000000..31365483d --- /dev/null +++ b/src/core/functions/analytic/codac2_analytic_flat_input_layout.cpp @@ -0,0 +1,172 @@ +/** + * \file codac2_analytic_flat_input_layout.cpp + * ---------------------------------------------------------------------------- + * \date 2026 + * \author Simon Rohou + * \copyright Copyright 2026 Codac Team + * \license GNU Lesser General Public License (LGPL) + */ + +#include "codac2_analytic_flat_input_layout.h" + +#include "codac2_assert.h" +#include "codac2_analytic_variables.h" +#include "codac2_component.h" + +namespace codac2 +{ + bool FlatInputBinding::is_scalar() const + { + return type == std::type_index(typeid(ScalarType)); + } + + bool FlatInputBinding::is_vector() const + { + return type == std::type_index(typeid(VectorType)); + } + + bool FlatInputBinding::is_matrix() const + { + return type == std::type_index(typeid(MatrixType)); + } + + FlatInputLayout::FlatInputLayout(const FunctionArgsList& args) + { + Index flat = 0; + + for(const auto& arg : args) + { + if(std::dynamic_pointer_cast(arg)) + { + _bindings.emplace(arg->unique_id().id(), FlatInputBinding{typeid(ScalarType), flat, 1, 1}); + ++flat; + } + + else if(auto v = std::dynamic_pointer_cast(arg)) + { + _bindings.emplace(arg->unique_id().id(), FlatInputBinding{typeid(VectorType), flat, v->size(), 1}); + flat += v->size(); + } + + else if(auto m = std::dynamic_pointer_cast(arg)) + { + _bindings.emplace(arg->unique_id().id(), FlatInputBinding{typeid(MatrixType), flat, m->rows(), m->cols()}); + flat += m->rows() * m->cols(); + } + + else + assert_release(false && "FlatInputLayout: unsupported variable type in function argument list"); + } + + _size = flat; + } + + Index FlatInputLayout::size() const + { + return _size; + } + + bool FlatInputLayout::same_domain_as(const FlatInputLayout& other) const + { + return size() == other.size(); + } + + bool FlatInputLayout::same_domain_as(const FunctionArgsList& other_args) const + { + return same_domain_as(FlatInputLayout(other_args)); + } + + Index FlatInputLayout::flat_index_of(const ScalarVar& x) const + { + const auto& b = binding_of(x.unique_id()); + assert_release(b.is_scalar() && "FlatInputLayout::flat_index_of: expected a scalar input variable"); + return b.offset; + } + + Index FlatInputLayout::flat_index_of(const VectorVar& x, Index i) const + { + const auto& b = binding_of(x.unique_id()); + assert_release(b.is_vector() && "FlatInputLayout::flat_index_of: expected a vector input variable"); + assert_release(i >= 0 && i < b.rows && "FlatInputLayout::flat_index_of: vector component out of bounds"); + return b.offset + i; + } + + Index FlatInputLayout::flat_index_of(const MatrixVar& x, Index i, Index j) const + { + const auto& b = binding_of(x.unique_id()); + assert_release(b.is_matrix() && "FlatInputLayout::flat_index_of: expected a matrix input variable"); + assert_release(i >= 0 && i < b.rows && j >= 0 && j < b.cols + && "FlatInputLayout::flat_index_of: matrix component out of bounds"); + return b.offset + b.cols*i + j; + } + + bool FlatInputLayout::flat_index_of(const ScalarExpr& x, Index& flat_index) const + { + if(auto s = std::dynamic_pointer_cast(x)) + { + const auto* b = find_binding(s->unique_id()); + if(!b || !b->is_scalar()) + return false; + + flat_index = b->offset; + return true; + } + + if(auto cv = std::dynamic_pointer_cast>(x)) + { + const auto children = cv->children_expr_base(); + if(children.size() != 1) + return false; + + auto v = std::dynamic_pointer_cast(children[0]); + if(!v) + return false; + + const auto* b = find_binding(v->unique_id()); + if(!b || !b->is_vector()) + return false; + + if(cv->i() < 0 || cv->i() >= b->rows) + return false; + + flat_index = b->offset + cv->i(); + return true; + } + + if(auto cm = std::dynamic_pointer_cast>(x)) + { + const auto children = cm->children_expr_base(); + if(children.size() != 1) + return false; + + auto m = std::dynamic_pointer_cast(children[0]); + if(!m) + return false; + + const auto* b = find_binding(m->unique_id()); + if(!b || !b->is_matrix()) + return false; + + if(cm->i() < 0 || cm->i() >= b->rows || cm->j() < 0 || cm->j() >= b->cols) + return false; + + flat_index = b->offset + b->cols*cm->i() + cm->j(); + return true; + } + + return false; + } + + const FlatInputBinding& FlatInputLayout::binding_of(const ExprID& id) const + { + const auto* b = find_binding(id); + assert_release(b != nullptr && "FlatInputLayout::binding_of: unknown input expression identifier"); + return *b; + } + + const FlatInputBinding* FlatInputLayout::find_binding(const ExprID& id) const + { + auto it = _bindings.find(id.id()); + return it == _bindings.end() ? nullptr : &it->second; + } +} diff --git a/src/core/functions/analytic/codac2_analytic_flat_input_layout.h b/src/core/functions/analytic/codac2_analytic_flat_input_layout.h new file mode 100644 index 000000000..191c8f5f9 --- /dev/null +++ b/src/core/functions/analytic/codac2_analytic_flat_input_layout.h @@ -0,0 +1,175 @@ +/** + * \file codac2_analytic_flat_input_layout.h + * ---------------------------------------------------------------------------- + * \date 2026 + * \author Simon Rohou + * \copyright Copyright 2026 Codac Team + * \license GNU Lesser General Public License (LGPL) + */ + +#pragma once + +#include +#include + +#include "codac2_AnalyticExprWrapper.h" +#include "codac2_AnalyticType.h" +#include "codac2_ExprBase.h" +#include "codac2_FunctionArgsList.h" + +namespace codac2 +{ + class ScalarVar; + class VectorVar; + class MatrixVar; + + /** + * \brief Binding information associated with one input argument in a flattened input domain. + * + * An analytic function input argument may be scalar, vectorial or matricial. + * In a flattened input domain, each argument is represented by a contiguous + * block of scalar inputs. + * + * The \p type field stores the analytic category of the argument and is one of: + * - typeid(ScalarType) + * - typeid(VectorType) + * - typeid(MatrixType) + */ + struct FlatInputBinding + { + std::type_index type = typeid(ScalarType); //!< Analytic type of the bound input argument. + Index offset = 0; //!< First scalar index of the argument in the flattened domain. + Index rows = 1; //!< Number of rows of the argument block. + Index cols = 1; //!< Number of columns of the argument block. + + /** + * \brief Tests whether the binding corresponds to a scalar input argument. + * + * \return True if the binding corresponds to a scalar input argument. + */ + bool is_scalar() const; + + /** + * \brief Tests whether the binding corresponds to a vector input argument. + * + * \return True if the binding corresponds to a vector input argument. + */ + bool is_vector() const; + + /** + * \brief Tests whether the binding corresponds to a matrix input argument. + * + * \return True if the binding corresponds to a matrix input argument. + */ + bool is_matrix() const; + }; + + /** + * \brief Flattened layout associated with an analytic function input domain. + * + * This class provides a canonical flattened representation of a function + * argument list. Each scalar, vector or matrix input argument is assigned + * a contiguous block of scalar indices. + */ + class FlatInputLayout + { + public: + + /** + * \brief Builds the flattened layout associated with a function argument list. + * + * \param args Function argument list. + */ + explicit FlatInputLayout(const FunctionArgsList& args); + + /** + * \brief Returns the total number of scalar inputs in the flattened domain. + * + * \return Flattened input domain size. + */ + Index size() const; + + /** + * \brief Tests whether two layouts describe the same flattened input domain. + * + * Here, only the total number of flattened scalar inputs matters. + * + * \param other Other flattened layout. + * \return True if both layouts describe the same flattened input domain. + */ + bool same_domain_as(const FlatInputLayout& other) const; + + /** + * \brief Tests whether this layout matches the flattened domain induced by a function argument list. + * + * \param other_args Function argument list. + * \return True if both domains are identical once flattened. + */ + bool same_domain_as(const FunctionArgsList& other_args) const; + + /** + * \brief Returns the flat index associated with a scalar input variable. + * + * This method is restricted to scalar variables appearing directly in the + * function argument list. + * + * \param x Scalar input variable. + * \return Flat index associated with \p x. + */ + Index flat_index_of(const ScalarVar& x) const; + + /** + * \brief Returns the flat index associated with a direct component of a vector input variable. + * + * \param x Vector input variable. + * \param i Component index. + * \return Flat index associated with \p x[i]. + */ + Index flat_index_of(const VectorVar& x, Index i) const; + + /** + * \brief Returns the flat index associated with a direct component of a matrix input variable. + * + * \param x Matrix input variable. + * \param i Row index. + * \param j Column index. + * \return Flat index associated with \p x(i,j). + */ + Index flat_index_of(const MatrixVar& x, Index i, Index j) const; + + /** + * \brief Tries to resolve a scalar input expression into a flat input index. + * + * Supported expressions are: + * - a scalar input variable; + * - a direct component of a vector input variable; + * - a direct component of a matrix input variable. + * + * \param x Scalar input expression. + * \param flat_index Output flat input index. + * \return True if \p x could be resolved into a flat input index. + */ + bool flat_index_of(const ScalarExpr& x, Index& flat_index) const; + + /** + * \brief Returns the binding associated with an input expression identifier. + * + * \param id Input expression identifier. + * \return Binding associated with \p id. + */ + const FlatInputBinding& binding_of(const ExprID& id) const; + + /** + * \brief Returns the binding associated with an input expression identifier, if any. + * + * \param id Input expression identifier. + * \return Pointer to the associated binding, or nullptr if \p id does not belong to the layout. + */ + const FlatInputBinding* find_binding(const ExprID& id) const; + + private: + + std::unordered_map _bindings; //!< Bindings indexed by input expression identifier. + Index _size = 0; //!< Total number of scalar inputs in the flattened domain. + }; +} diff --git a/src/core/functions/codac2_ExprBase.h b/src/core/functions/codac2_ExprBase.h index eb3931929..a7b6fa08d 100644 --- a/src/core/functions/codac2_ExprBase.h +++ b/src/core/functions/codac2_ExprBase.h @@ -141,6 +141,19 @@ namespace codac2 */ virtual ~ExprBase() = default; + /** + * Returns the direct children of this expression node. + * + * The children are returned in the same order as the operands of the + * underlying operator. + * + * \return A vector of child expressions. + */ + virtual std::vector> children_expr_base() const + { + return {}; + } + protected: const ExprID _unique_id; //!< unique identifier for this expression @@ -216,6 +229,26 @@ namespace codac2 }, _x); } + /** + * Returns the direct children of this expression node. + * + * The children are returned in the same order as the operands of the + * underlying operator. + * + * \return A vector of child expressions. + */ + std::vector> children_expr_base() const + { + std::vector> children; + children.reserve(sizeof...(X)); + std::apply( + [&children](const auto&... x) + { + (children.push_back(std::static_pointer_cast(x)), ...); + }, _x); + return children; + } + protected: /** diff --git a/src/core/operators/codac2_component.h b/src/core/operators/codac2_component.h index 8c6db15a3..16ce50dcc 100644 --- a/src/core/operators/codac2_component.h +++ b/src/core/operators/codac2_component.h @@ -117,6 +117,16 @@ namespace codac2 { return true; } + + Index i() const + { + return _i; + } + + std::vector> children_expr_base() const override + { + return OperationExprBase>::children_expr_base(); + } protected: @@ -183,6 +193,21 @@ namespace codac2 return true; } + Index i() const + { + return _i; + } + + Index j() const + { + return _j; + } + + std::vector> children_expr_base() const override + { + return OperationExprBase>::children_expr_base(); + } + protected: const Index _i, _j; diff --git a/src/core/operators/codac2_mat.h b/src/core/operators/codac2_mat.h index 0be257079..c61e1cbb4 100644 --- a/src/core/operators/codac2_mat.h +++ b/src/core/operators/codac2_mat.h @@ -9,6 +9,8 @@ #pragma once +#include + #include "codac2_IntervalVector.h" #include "codac2_IntervalMatrix.h" #include "codac2_AnalyticType.h" @@ -63,7 +65,7 @@ namespace codac2 } static inline void fill_diff_matrix(IntervalMatrix &d, - const IntervalMatrix &dax, Index &l) { + const IntervalMatrix &dax, Index &l) { d.middleRows(l,dax.rows())=dax; l += dax.rows(); } @@ -71,7 +73,7 @@ namespace codac2 template requires (std::is_base_of_v - && (std::is_base_of_v && ...)) + && (std::is_base_of_v && ...)) static inline MatrixType fwd_centered(const X1& x1, const X&... x) { if (centered_form_not_available_for_args(x1,x...)) @@ -83,7 +85,7 @@ namespace codac2 l += x1.da.rows(); ( MatrixOp::fill_diff_matrix(d,x.da,l) , ...); assert (l==d.rows()); - + return { MatrixOp::fwd(x1.m,x.m...), MatrixOp::fwd(x1.a,x.a...), @@ -101,10 +103,177 @@ namespace codac2 } }; + namespace detail + { + inline void replace_vector_child(std::shared_ptr>& x, + const ExprID& old_arg_id, const std::shared_ptr& new_expr) + { + if(x->unique_id() == old_arg_id) + { + auto new_x = std::dynamic_pointer_cast>(new_expr); + assert_release(new_x); + x = new_x; + } + + else + x->replace_arg(old_arg_id, new_expr); + } + + class DynamicMatrixExpr final : public AnalyticExpr + { + public: + + explicit DynamicMatrixExpr(const std::vector& xs) + { + _xs.reserve(xs.size()); + for(const auto& x : xs) + { + auto vx = std::dynamic_pointer_cast>(x); + assert_release(vx); + _xs.push_back(vx); + } + } + + DynamicMatrixExpr(const DynamicMatrixExpr& e) + { + _xs.reserve(e._xs.size()); + for(const auto& x : e._xs) + { + auto vx = std::dynamic_pointer_cast>(x->copy()); + assert_release(vx); + _xs.push_back(vx); + } + } + + std::shared_ptr copy() const override + { + return std::make_shared(*this); + } + + void replace_arg(const ExprID& old_arg_id, const std::shared_ptr& new_expr) override + { + for(auto& x : _xs) + replace_vector_child(x, old_arg_id, new_expr); + } + + MatrixType fwd_eval(ValuesMap& v, Index total_input_size, bool natural_eval) const override + { + if(natural_eval) + return this->init_value(v, natural_fwd(v, total_input_size)); + + std::vector vals; + vals.reserve(_xs.size()); + bool centered_available = true; + for(const auto& x : _xs) + { + vals.push_back(x->fwd_eval(v, total_input_size, false)); + centered_available &= (vals.back().da.size() != 0); + } + + if(!centered_available) + 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 input_cols = vals.empty() ? total_input_size : vals.front().da.cols(); + IntervalMatrix m(rows, cols), a(rows, cols), da(rows*cols, input_cols); + bool def_domain = true; + Index l = 0; + + for(Index j = 0 ; j < cols ; ++j) + { + const auto& xj = vals[static_cast(j)]; + m.col(j) = xj.m; + a.col(j) = xj.a; + da.middleRows(l, xj.da.rows()) = xj.da; + l += xj.da.rows(); + def_domain &= xj.def_domain; + } + + return this->init_value(v, MatrixType(m, a, da, def_domain)); + } + + void bwd_eval(ValuesMap& v) const override + { + for(const auto& x : _xs) + x->bwd_eval(v); + } + + std::pair output_shape() const override + { + if(_xs.empty()) + return { 0, 0 }; + + const auto shape = _xs.front()->output_shape(); + assert(shape.second == 1); + return { shape.first, static_cast(_xs.size()) }; + } + + bool belongs_to_args_list(const FunctionArgsList& args) const override + { + bool ok = true; + for(const auto& x : _xs) + ok &= x->belongs_to_args_list(args); + return ok; + } + + std::string str(bool in_parentheses = false) const override + { + if(_xs.empty()) + return in_parentheses ? "([])" : "[]"; + + std::string s; + for(const auto& x : _xs) + s += "\t" + x->str() + ",\n"; + + s.pop_back(); + s.pop_back(); + + s = "[\n" + s + "\n]"; + return in_parentheses ? "(" + s + ")" : s; + } + + bool is_str_leaf() const override + { + return false; + } + + std::vector> children_expr_base() const override + { + std::vector> children; + children.reserve(_xs.size()); + for(const auto& x : _xs) + children.push_back(std::dynamic_pointer_cast(x)); + return children; + } + + private: + + MatrixType natural_fwd(ValuesMap& v, Index total_input_size) const + { + const Index cols = static_cast(_xs.size()); + const Index rows = _xs.empty() ? 0 : _xs.front()->output_shape().first; + IntervalMatrix a(rows, cols); + bool def_domain = true; + + for(Index j = 0 ; j < cols ; ++j) + { + auto xj = _xs[static_cast(j)]->fwd_eval(v, total_input_size, true); + a.col(j) = xj.a; + def_domain &= xj.def_domain; + } + + return { a, def_domain }; + } + + std::vector>> _xs; + }; + } + // Analytic operator // The following function can be used to build analytic expressions. - // Generic variadic case, cannot handle const values (int, double) for now + // Variadic (cannot handle const values (int, double) for now) template inline MatrixExpr @@ -113,4 +282,12 @@ namespace codac2 return { std::make_shared>( AnalyticOperationExpr(x...)) }; } + + // Dynamic: + + inline MatrixExpr + mat(const std::vector& x) + { + return { std::make_shared(x) }; + } } diff --git a/src/core/operators/codac2_vec.h b/src/core/operators/codac2_vec.h index 973a1ce4a..8e862e5c2 100644 --- a/src/core/operators/codac2_vec.h +++ b/src/core/operators/codac2_vec.h @@ -9,8 +9,11 @@ #pragma once +#include + #include "codac2_Interval.h" #include "codac2_IntervalVector.h" +#include "codac2_IntervalMatrix.h" #include "codac2_AnalyticType.h" #include "codac2_AnalyticExprWrapper.h" @@ -45,7 +48,7 @@ namespace codac2 { bool def_domain = true; ((def_domain &= x.def_domain), ...); - + return { fwd(x.a...), def_domain @@ -65,7 +68,7 @@ namespace codac2 bool def_domain = true; ((def_domain &= x.def_domain), ...); - + return { fwd(x.m...), fwd(x.a...), @@ -83,6 +86,164 @@ namespace codac2 } }; + namespace detail + { + inline void replace_scalar_child(std::shared_ptr>& x, + const ExprID& old_arg_id, const std::shared_ptr& new_expr) + { + if(x->unique_id() == old_arg_id) + { + auto new_x = std::dynamic_pointer_cast>(new_expr); + assert_release(new_x); + x = new_x; + } + + else + x->replace_arg(old_arg_id, new_expr); + } + + class DynamicVectorExpr final : public AnalyticExpr + { + public: + + explicit DynamicVectorExpr(const std::vector& xs) + { + _xs.reserve(xs.size()); + for(const auto& x : xs) + { + auto sx = std::dynamic_pointer_cast>(x); + assert_release(sx); + _xs.push_back(sx); + } + } + + DynamicVectorExpr(const DynamicVectorExpr& e) + { + _xs.reserve(e._xs.size()); + for(const auto& x : e._xs) + { + auto sx = std::dynamic_pointer_cast>(x->copy()); + assert_release(sx); + _xs.push_back(sx); + } + } + + std::shared_ptr copy() const override + { + return std::make_shared(*this); + } + + void replace_arg(const ExprID& old_arg_id, const std::shared_ptr& new_expr) override + { + for(auto& x : _xs) + replace_scalar_child(x, old_arg_id, new_expr); + } + + VectorType fwd_eval(ValuesMap& v, Index total_input_size, bool natural_eval) const override + { + if(natural_eval) + return this->init_value(v, natural_fwd(v, total_input_size)); + + std::vector vals; + vals.reserve(_xs.size()); + bool centered_available = true; + for(const auto& x : _xs) + { + vals.push_back(x->fwd_eval(v, total_input_size, false)); + centered_available &= (vals.back().da.size() != 0); + } + + if(!centered_available) + return this->init_value(v, natural_fwd(v, total_input_size)); + + 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); + IntervalMatrix da(n, input_cols); + bool def_domain = true; + + for(Index i = 0 ; i < n ; ++i) + { + const auto& xi = vals[static_cast(i)]; + m[i] = xi.m; + a[i] = xi.a; + da.row(i) = xi.da; + def_domain &= xi.def_domain; + } + + return this->init_value(v, VectorType(m, a, da, def_domain)); + } + + void bwd_eval(ValuesMap& v) const override + { + for(const auto& x : _xs) + x->bwd_eval(v); + } + + std::pair output_shape() const override + { + return { static_cast(_xs.size()), 1 }; + } + + bool belongs_to_args_list(const FunctionArgsList& args) const override + { + bool ok = true; + for(const auto& x : _xs) + ok &= x->belongs_to_args_list(args); + return ok; + } + + std::string str(bool in_parentheses = false) const override + { + if(_xs.empty()) + return in_parentheses ? "([])" : "[]"; + + std::string s; + for(const auto& x : _xs) + s += "\t" + x->str() + ",\n"; + + s.pop_back(); + s.pop_back(); + + s = "[\n" + s + "\n]"; + return in_parentheses ? "(" + s + ")" : s; + } + + bool is_str_leaf() const override + { + return false; + } + + std::vector> children_expr_base() const override + { + std::vector> children; + children.reserve(_xs.size()); + for(const auto& x : _xs) + children.push_back(std::dynamic_pointer_cast(x)); + return children; + } + + private: + + VectorType natural_fwd(ValuesMap& v, Index total_input_size) const + { + const Index n = static_cast(_xs.size()); + IntervalVector a(n); + 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; + def_domain &= xi.def_domain; + } + + return { a, def_domain }; + } + + std::vector>> _xs; + }; + } + // Analytic operator // The following functions can be used to build analytic expressions. @@ -96,6 +257,8 @@ namespace codac2 return const_value(x); } + // Variadic: + template requires ((std::is_same_v::Type,ScalarType>) && ...) inline VectorExpr @@ -103,4 +266,12 @@ namespace codac2 { return { std::make_shared::Type...>>(_add_to_vec(x)...) }; } + + // Dynamic: + + inline VectorExpr + vec(const std::vector& x) + { + return { std::make_shared(x) }; + } } diff --git a/src/extensions/sympy/CMakeLists.txt b/src/extensions/sympy/CMakeLists.txt index 08a8f67fe..d92bae572 100644 --- a/src/extensions/sympy/CMakeLists.txt +++ b/src/extensions/sympy/CMakeLists.txt @@ -2,74 +2,62 @@ # Codac - cmake configuration file # ================================================================== + set(CODAC_SYMPY_PUBLIC_HDR + ${CMAKE_CURRENT_SOURCE_DIR}/codac2_sympy.h + ) + + set(CODAC_SYMPY_PRIVATE_HDR + ${CMAKE_CURRENT_SOURCE_DIR}/codac2_sympy_bridge.h + ) + list(APPEND CODAC_SYMPY_SRC - ${CMAKE_CURRENT_SOURCE_DIR}/codac2_sympy_empty.cpp - ${CMAKE_CURRENT_SOURCE_DIR}/codac2_sympy_empty.h + ${CMAKE_CURRENT_SOURCE_DIR}/codac2_sympy_bridge.cpp + ${CODAC_SYMPY_PRIVATE_HDR} + ${CMAKE_CURRENT_SOURCE_DIR}/codac2_sympy.cpp + ${CODAC_SYMPY_PUBLIC_HDR} ) ################################################################################ # Looking for pybind11 ################################################################################ - add_subdirectory( - ${CMAKE_CURRENT_SOURCE_DIR}/../../python/pybind11 - ${CMAKE_CURRENT_BINARY_DIR}/pybind11 - ) + # Already added with WITH_PYTHON ################################################################################ # Create the target for libcodac-sympy ################################################################################ - #if(NOT CMAKE_CXX_STANDARD) - # set(CMAKE_CXX_STANDARD 20) - # set(CMAKE_CXX_STANDARD_REQUIRED ON) - #endif() + set(CMAKE_CXX_STANDARD 20) + set(CMAKE_CXX_STANDARD_REQUIRED ON) - #add_library(${PROJECT_NAME}-sympy - - pybind11_add_module(${PROJECT_NAME}-sympy MODULE - ${CODAC_SYMPY_SRC}) - target_include_directories(${PROJECT_NAME}-sympy PUBLIC - ${CMAKE_CURRENT_SOURCE_DIR}/ - ) - target_link_libraries(${PROJECT_NAME}-sympy PUBLIC) # pybind11::embed) - + add_library(${PROJECT_NAME}-sympy ${CODAC_SYMPY_SRC}) + target_link_libraries(${PROJECT_NAME}-sympy PUBLIC ${PROJECT_NAME}-core Ibex::ibex Eigen3::Eigen) + target_link_libraries(${PROJECT_NAME}-sympy PRIVATE pybind11::pybind11) + + ################################################################################ + # For the generation of the PKG file + ################################################################################ -################################################################################ -# For the generation of the PKG file -################################################################################ - set(CODAC_PKG_CONFIG_CFLAGS "${CODAC_PKG_CONFIG_CFLAGS} -I\${includedir}/${PROJECT_NAME}-sympy" PARENT_SCOPE) set(CODAC_PKG_CONFIG_LIBS "${CODAC_PKG_CONFIG_LIBS} -l${PROJECT_NAME}-sympy" PARENT_SCOPE) - -################################################################################ -# Installation of libcodac-sympy files -################################################################################ - -# Getting header files from sources + ################################################################################ + # Installation of libcodac-sympy files + ################################################################################ - foreach(srcfile ${CODAC_SYMPY_SRC}) - if(srcfile MATCHES "\\.h$" OR srcfile MATCHES "\\.hpp$") - list(APPEND CODAC_SYMPY_HDR ${srcfile}) - file(COPY ${srcfile} DESTINATION ${CMAKE_CURRENT_BINARY_DIR}/../../include) - endif() + foreach(srcfile ${CODAC_SYMPY_PUBLIC_HDR}) + file(COPY ${srcfile} DESTINATION ${CMAKE_CURRENT_BINARY_DIR}/../../../include) endforeach() -# Generating the file codac-sympy.h - set(CODAC_SYMPY_MAIN_HEADER ${CMAKE_CURRENT_BINARY_DIR}/codac-sympy.h) - set(CODAC_MAIN_SUBHEADERS ${CODAC_MAIN_SUBHEADERS} "codac-sympy.h" PARENT_SCOPE) file(WRITE ${CODAC_SYMPY_MAIN_HEADER} "/* This file is generated by CMake */\n\n") file(APPEND ${CODAC_SYMPY_MAIN_HEADER} "#pragma once\n\n") - foreach(header_path ${CODAC_SYMPY_HDR}) + foreach(header_path ${CODAC_SYMPY_PUBLIC_HDR}) get_filename_component(header_name ${header_path} NAME) file(APPEND ${CODAC_SYMPY_MAIN_HEADER} "#include <${header_name}>\n") endforeach() - file(COPY ${CODAC_SYMPY_MAIN_HEADER} DESTINATION ${CMAKE_CURRENT_BINARY_DIR}/../../include) - -# Install files in system directories + file(COPY ${CODAC_SYMPY_MAIN_HEADER} DESTINATION ${CMAKE_CURRENT_BINARY_DIR}/../../../include) install(TARGETS ${PROJECT_NAME}-sympy DESTINATION ${CMAKE_INSTALL_LIBDIR}) - install(FILES ${CODAC_SYMPY_HDR} DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/${PROJECT_NAME}-sympy) - install(FILES ${CODAC_SYMPY_MAIN_HEADER} DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/${PROJECT_NAME}-sympy) \ No newline at end of file + install(FILES ${CODAC_SYMPY_PUBLIC_HDR} DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/${PROJECT_NAME}-sympy) + install(FILES ${CODAC_SYMPY_MAIN_HEADER} DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/${PROJECT_NAME}-sympy) diff --git a/src/extensions/sympy/codac2_sympy.cpp b/src/extensions/sympy/codac2_sympy.cpp new file mode 100644 index 000000000..5226ca02b --- /dev/null +++ b/src/extensions/sympy/codac2_sympy.cpp @@ -0,0 +1,472 @@ +/** + * codac2_sympy.cpp + * ---------------------------------------------------------------------------- + * \date 2026 + * \author Simon Rohou, Maël Godard + * \copyright Copyright 2026 Codac Team + * \license GNU Lesser General Public License (LGPL) + */ + +#include "codac2_sympy.h" +#include "codac2_sympy_bridge.h" + +#include "codac2_analytic_componentwise.h" +#include "codac2_analytic_flat_input_layout.h" + +namespace codac2 +{ + namespace + { + using symbolic::detail::FlatSymbolTable; + using symbolic::detail::ScalarBridgeContext; + using symbolic::detail::SympyTransform; + + template + AnalyticFunction transform_componentwise( + const AnalyticFunction& f, + const SympyTransform& transform, + bool do_expand = true) + { + ScalarBridgeContext ctx(f.args()); + return map_scalar_entries(f, [&](const ScalarExpr& y) + { + return ctx.transform_scalar_expr(y, transform, do_expand); + }); + } + + Index flat_input_index_of(const FunctionArgsList& args, const ScalarExpr& x) + { + const FlatInputLayout layout(args); + Index flat_input_index = -1; + if(layout.flat_index_of(x, flat_input_index)) + return flat_input_index; + assert(false && "flat_input_index_of: expected a scalar input variable or a direct component of a vector/matrix input variable"); + } + + pybind11::object remap_to_reference_symbols( + pybind11::object expr, + const FlatSymbolTable& source_symbols, + const FlatSymbolTable& reference_symbols) + { + if(source_symbols.size() != reference_symbols.size()) + { + assert(false && "sympy_equal: inconsistent flattened symbol tables"); + } + + pybind11::dict subs; + for(Index k = 0 ; k < source_symbols.size() ; ++k) + subs[source_symbols.by_flat_index(k)] = reference_symbols.by_flat_index(k); + + return expr.attr("subs")(subs); + } + + bool sympy_zero(const pybind11::object& expr) + { + const pybind11::object& sympy = symbolic::detail::import_sympy(); + pybind11::object simplified = sympy.attr("simplify")(expr); + + try + { + simplified = sympy.attr("nsimplify")(simplified, pybind11::arg("rational")=true); + } + catch(const pybind11::error_already_set&) + { + // Keep the simplified expression if SymPy cannot rationalize it. + } + + pybind11::object is_zero = simplified.attr("is_zero"); + if(!is_zero.is_none()) + return pybind11::cast(is_zero); + + pybind11::object equals_zero = simplified.attr("equals")(pybind11::int_(0)); + if(!equals_zero.is_none()) + return pybind11::cast(equals_zero); + + return false; + } + + bool sympy_equal_scalar_expr( + const FunctionArgsList& args_f, + const ScalarExpr& yf, + const FunctionArgsList& args_g, + const ScalarExpr& yg) + { + const FlatInputLayout layout_f(args_f); + if(!layout_f.same_domain_as(args_g)) + return false; + + symbolic::detail::ensure_python_runtime(); + pybind11::gil_scoped_acquire gil; + + ScalarBridgeContext reference_ctx(args_f); + ScalarBridgeContext source_ctx(args_g); + + pybind11::object ef = reference_ctx.export_scalar(yf); + pybind11::object eg = source_ctx.export_scalar(yg); + eg = remap_to_reference_symbols(eg, source_ctx.symbols(), reference_ctx.symbols()); + + return sympy_zero(ef - eg); + } + + ScalarExpr sympy_partial_diff_expr( + const ScalarBridgeContext& ctx, + const ScalarExpr& y, + Index flat_input_index) + { + if(flat_input_index < 0 || flat_input_index >= ctx.symbols().size()) + { + assert(false && "sympy_partial_diff_expr: flat_input_index out of bounds"); + } + + return ctx.transform_scalar_expr( + y, + [flat_input_index](const pybind11::object& sympy, + const pybind11::object& ys, + const FlatSymbolTable& symbols) + { + return sympy.attr("diff")(ys, symbols.by_flat_index(flat_input_index)); + }); + } + } + + AnalyticFunction + sympy_simplify(const AnalyticFunction& f) + { + return transform_componentwise(f, + []([[maybe_unused]] const pybind11::object& sympy, const pybind11::object& ys, const FlatSymbolTable&) + { + return sympy.attr("simplify")(ys); + }, + true); + } + + AnalyticFunction + sympy_simplify(const AnalyticFunction& f) + { + return transform_componentwise(f, + []([[maybe_unused]] const pybind11::object& sympy, const pybind11::object& ys, const FlatSymbolTable&) + { + return sympy.attr("simplify")(ys); + }, + true); + } + + AnalyticFunction + sympy_simplify(const AnalyticFunction& f) + { + return transform_componentwise(f, + []([[maybe_unused]] const pybind11::object& sympy, const pybind11::object& ys, const FlatSymbolTable&) + { + return sympy.attr("simplify")(ys); + }, + true); + } + + AnalyticFunction + sympy_horner(const AnalyticFunction& f) + { + return transform_componentwise(f, + []([[maybe_unused]] const pybind11::object& sympy, const pybind11::object& ys, const FlatSymbolTable&) + { + return symbolic::detail::import_polyfuncs().attr("horner")(ys); + }, + false); + } + + AnalyticFunction + sympy_horner(const AnalyticFunction& f) + { + return transform_componentwise(f, + []([[maybe_unused]] const pybind11::object& sympy, const pybind11::object& ys, const FlatSymbolTable&) + { + return symbolic::detail::import_polyfuncs().attr("horner")(ys); + }, + false); + } + + AnalyticFunction + sympy_horner(const AnalyticFunction& f) + { + return transform_componentwise(f, + []([[maybe_unused]] const pybind11::object& sympy, const pybind11::object& ys, const FlatSymbolTable&) + { + return symbolic::detail::import_polyfuncs().attr("horner")(ys); + }, + false); + } + + AnalyticFunction + sympy_partial_diff(const AnalyticFunction& f, Index flat_input_index) + { + if(flat_input_index < 0 || flat_input_index >= f.input_size()) + { + assert(false && "sympy_partial_diff: flat_input_index out of bounds"); + } + + ScalarBridgeContext ctx(f.args()); + return AnalyticFunction( + f.args(), + sympy_partial_diff_expr(ctx, ScalarExpr(f.expr()), flat_input_index)); + } + + AnalyticFunction + sympy_partial_diff(const AnalyticFunction& f, const ScalarVar& x) + { + return sympy_partial_diff(f, ScalarExpr(x)); + } + + AnalyticFunction + sympy_partial_diff(const AnalyticFunction& f, const ScalarExpr& x) + { + return sympy_partial_diff(f, flat_input_index_of(f.args(),x)); + } + + AnalyticFunction + sympy_diff(const AnalyticFunction& f) + { + if(f.input_size() != 1) + { + assert(false && "sympy_diff(f): this overload supports only scalar functions with one flattened input. \ + Use sympy_partial_diff(f,j), sympy_diff(f,order), or sympy_gradient(f) for multivariate scalar functions."); + } + + return sympy_partial_diff(f, 0); + } + + AnalyticFunction + sympy_diff(const AnalyticFunction& f, const ScalarVar& x) + { + return sympy_diff(f, ScalarExpr(x)); + } + + AnalyticFunction + sympy_diff(const AnalyticFunction& f, const ScalarExpr& x) + { + return sympy_partial_diff(f, x); + } + + AnalyticFunction + sympy_diff(const AnalyticFunction& f, Index order) + { + if(order < 0) + { + assert(false && "sympy_diff(f,order): order must be nonnegative"); + } + + if(order == 0) + return f; + + if(f.input_size() != 1) + { + assert(false && "sympy_diff(f,order): this overload supports only scalar functions with one flattened input"); + } + + ScalarBridgeContext ctx(f.args()); + ScalarExpr y(f.expr()); + for(Index k = 0 ; k < order ; ++k) + y = sympy_partial_diff_expr(ctx, y, 0); + + return AnalyticFunction(f.args(), y); + } + + AnalyticFunction + sympy_diff(const AnalyticFunction& f, const ScalarVar& x, Index order) + { + return sympy_diff(f, ScalarExpr(x), order); + } + + AnalyticFunction + sympy_diff(const AnalyticFunction& f, const ScalarExpr& x, Index order) + { + if(order < 0) + { + assert(false && "sympy_diff(f,x,order): order must be nonnegative"); + } + + if(order == 0) + return f; + + const Index flat_input_index = flat_input_index_of(f.args(),x); + + ScalarBridgeContext ctx(f.args()); + ScalarExpr y(f.expr()); + for(Index k = 0 ; k < order ; ++k) + y = sympy_partial_diff_expr(ctx, y, flat_input_index); + + return AnalyticFunction(f.args(), y); + } + + AnalyticFunction + sympy_gradient(const AnalyticFunction& f) + { + ScalarBridgeContext ctx(f.args()); + std::vector entries; + entries.reserve(static_cast(f.input_size())); + const ScalarExpr y(f.expr()); + + for(Index j = 0 ; j < f.input_size() ; ++j) + entries.push_back(sympy_partial_diff_expr(ctx, y, j)); + + return AnalyticFunction(f.args(), vec(entries)); + } + + AnalyticFunction + sympy_hessian(const AnalyticFunction& f) + { + ScalarBridgeContext ctx(f.args()); + std::vector cols; + cols.reserve(static_cast(f.input_size())); + const ScalarExpr y(f.expr()); + + for(Index j = 0 ; j < f.input_size() ; ++j) + { + const ScalarExpr dj = sympy_partial_diff_expr(ctx, y, j); + std::vector col_entries; + col_entries.reserve(static_cast(f.input_size())); + for(Index i = 0 ; i < f.input_size() ; ++i) + col_entries.push_back(sympy_partial_diff_expr(ctx, dj, i)); + cols.push_back(vec(col_entries)); + } + + return AnalyticFunction(f.args(), mat(cols)); + } + + AnalyticFunction + sympy_diff(const AnalyticFunction& f) + { + const auto shape = f.output_shape(); + if(shape.second != 1) + { + assert(false && "sympy_diff(VectorType): only column-vector outputs are supported"); + } + + ScalarBridgeContext ctx(f.args()); + std::vector cols; + cols.reserve(static_cast(f.input_size())); + VectorExpr y(f.expr()); + + for(Index j = 0 ; j < f.input_size() ; ++j) + { + std::vector col_entries; + col_entries.reserve(static_cast(shape.first)); + for(Index i = 0 ; i < shape.first ; ++i) + col_entries.push_back(sympy_partial_diff_expr(ctx, y[i], j)); + cols.push_back(vec(col_entries)); + } + + return AnalyticFunction(f.args(), mat(cols)); + } + + AnalyticFunction + sympy_series(const AnalyticFunction& f, double center, Index order) + { + if(order < 0) + { + assert(false && "sympy_series: order must be nonnegative"); + } + + if(f.input_size() != 1) + { + assert(false && "sympy_series: this overload supports only scalar functions with one flattened input"); + } + + return AnalyticFunction( + f.args(), + symbolic::detail::transform_scalar_expr( + f.args(), ScalarExpr(f.expr()), + [center, order](const pybind11::object& sympy, + const pybind11::object& ys, + const FlatSymbolTable& symbols) + { + pybind11::object x0 = symbols.by_flat_index(0); + pybind11::object s = sympy.attr("series")(ys, x0, pybind11::float_(center), pybind11::int_(order+1)); + pybind11::object poly = s.attr("removeO")(); + return sympy.attr("expand")(poly); + })); + } + + AnalyticFunction + sympy_series(const AnalyticFunction& f, const ScalarVar& x, double center, Index order) + { + return sympy_series(f, ScalarExpr(x), center, order); + } + + AnalyticFunction + sympy_series(const AnalyticFunction& f, const ScalarExpr& x, double center, Index order) + { + if(order < 0) + { + assert(false && "sympy_series: order must be nonnegative"); + } + + const Index flat_input_index = flat_input_index_of(f.args(),x); + + return AnalyticFunction( + f.args(), + symbolic::detail::transform_scalar_expr( + f.args(), ScalarExpr(f.expr()), + [center, order, flat_input_index](const pybind11::object& sympy, + const pybind11::object& ys, + const FlatSymbolTable& symbols) + { + pybind11::object xj = symbols.by_flat_index(flat_input_index); + pybind11::object s = sympy.attr("series")(ys, xj, pybind11::float_(center), pybind11::int_(order+1)); + pybind11::object poly = s.attr("removeO")(); + return sympy.attr("expand")(poly); + })); + } + + bool sympy_equal(const AnalyticFunction& f, const AnalyticFunction& g) + { + return sympy_equal_scalar_expr(f.args(), ScalarExpr(f.expr()), g.args(), ScalarExpr(g.expr())); + } + + bool sympy_equal(const AnalyticFunction& f, const AnalyticFunction& g) + { + const FlatInputLayout layout_f(f.args()); + if(!layout_f.same_domain_as(g.args())) + return false; + + const auto shape_f = f.output_shape(); + const auto shape_g = g.output_shape(); + if(shape_f != shape_g) + return false; + + VectorExpr yf(f.expr()); + VectorExpr yg(g.expr()); + + for(Index i = 0 ; i < shape_f.first ; ++i) + { + if(!sympy_equal_scalar_expr(f.args(), yf[i], g.args(), yg[i])) + return false; + } + + return true; + } + + bool sympy_equal(const AnalyticFunction& f, const AnalyticFunction& g) + { + const FlatInputLayout layout_f(f.args()); + if(!layout_f.same_domain_as(g.args())) + return false; + + const auto shape_f = f.output_shape(); + const auto shape_g = g.output_shape(); + if(shape_f != shape_g) + return false; + + MatrixExpr yf(f.expr()); + MatrixExpr yg(g.expr()); + + for(Index j = 0 ; j < shape_f.second ; ++j) + { + for(Index i = 0 ; i < shape_f.first ; ++i) + { + if(!sympy_equal_scalar_expr(f.args(), yf(i,j), g.args(), yg(i,j))) + return false; + } + } + + return true; + } +} diff --git a/src/extensions/sympy/codac2_sympy.h b/src/extensions/sympy/codac2_sympy.h new file mode 100644 index 000000000..d0661c478 --- /dev/null +++ b/src/extensions/sympy/codac2_sympy.h @@ -0,0 +1,287 @@ +/** + * \file codac2_sympy.h + * ---------------------------------------------------------------------------- + * \date 2026 + * \author Simon Rohou, Maël Godard + * \copyright Copyright 2026 Codac Team + * \license GNU Lesser General Public License (LGPL) + */ + +#pragma once + +#include +#include "codac2_AnalyticFunction.h" + +namespace codac2 +{ + /** + * \brief Symbolically simplifies a scalar analytic function. + * + * \param f Function to simplify. + * \return A simplified analytic function. + */ + AnalyticFunction + sympy_simplify(const AnalyticFunction& f); + + /** + * \brief Symbolically simplifies a vector analytic function componentwise. + * + * \param f Function to simplify. + * \return A simplified analytic function. + */ + AnalyticFunction + sympy_simplify(const AnalyticFunction& f); + + /** + * \brief Symbolically simplifies a matrix analytic function componentwise. + * + * \param f Function to simplify. + * \return A simplified analytic function. + */ + AnalyticFunction + sympy_simplify(const AnalyticFunction& f); + + /** + * \brief Rewrites a scalar analytic function in Horner form when possible. + * + * \param f Function to rewrite. + * \return A rewritten analytic function. + */ + AnalyticFunction + sympy_horner(const AnalyticFunction& f); + + /** + * \brief Rewrites a vector analytic function in Horner form componentwise. + * + * \param f Function to rewrite. + * \return A rewritten analytic function. + */ + AnalyticFunction + sympy_horner(const AnalyticFunction& f); + + /** + * \brief Rewrites a matrix analytic function in Horner form componentwise. + * + * \param f Function to rewrite. + * \return A rewritten analytic function. + */ + AnalyticFunction + sympy_horner(const AnalyticFunction& f); + + /** + * \brief Returns the symbolic partial derivative of a scalar function. + * + * The derivative is taken with respect to the flattened scalar input + * of index \p flat_input_index. + * + * \param f Scalar function. + * \param flat_input_index Flattened input index. + * \return The symbolic partial derivative. + */ + AnalyticFunction + sympy_partial_diff(const AnalyticFunction& f, Index flat_input_index); + + /** + * \brief Returns the symbolic partial derivative of a scalar function + * with respect to a scalar input variable. + * + * This overload is restricted to scalar variables that appear directly + * in the function argument list. + * + * \param f Scalar function. + * \param x Scalar input variable. + * \return The symbolic partial derivative. + */ + AnalyticFunction + sympy_partial_diff(const AnalyticFunction& f, const ScalarVar& x); + + /** + * \brief Returns the symbolic partial derivative of a scalar function + * with respect to a scalar input expression. + * + * This overload supports either a scalar input variable, or a direct scalar + * component of a vectorial or matricial input variable. + * + * \param f Scalar function. + * \param x Scalar input expression. + * \return The symbolic partial derivative. + */ + AnalyticFunction + sympy_partial_diff(const AnalyticFunction& f, const ScalarExpr& x); + + /** + * \brief Returns the symbolic derivative of a scalar univariate function. + * + * \param f Scalar function. + * \return The symbolic derivative. + */ + AnalyticFunction + sympy_diff(const AnalyticFunction& f); + + /** + * \brief Returns the symbolic first derivative of a scalar function with + * respect to a scalar input variable. + * + * \param f Scalar function. + * \param x Scalar input variable. + * \return The symbolic derivative with respect to \p x. + */ + AnalyticFunction + sympy_diff(const AnalyticFunction& f, const ScalarVar& x); + + /** + * \brief Returns the symbolic first derivative of a scalar function with + * respect to a scalar input expression. + * + * This overload supports either a scalar input variable, or a direct scalar + * component of a vectorial or matricial input variable. + * + * \param f Scalar function. + * \param x Scalar input expression. + * \return The symbolic derivative with respect to \p x. + */ + AnalyticFunction + sympy_diff(const AnalyticFunction& f, const ScalarExpr& x); + + /** + * \brief Returns the symbolic derivative of given order for a scalar univariate function. + * + * \param f Scalar function. + * \param order Derivation order. + * \return The symbolic derivative of order \p order. + */ + AnalyticFunction + sympy_diff(const AnalyticFunction& f, Index order); + + /** + * \brief Returns the symbolic derivative of given order for a scalar + * function with respect to a scalar input variable. + * + * \param f Scalar function. + * \param x Scalar input variable. + * \param order Derivation order. + * \return The symbolic derivative of order \p order with respect to \p x. + */ + AnalyticFunction + sympy_diff(const AnalyticFunction& f, const ScalarVar& x, Index order); + + /** + * \brief Returns the symbolic derivative of given order for a scalar + * function with respect to a scalar input expression. + * + * This overload supports either a scalar input variable, or a direct scalar + * component of a vectorial or matricial input variable. + * + * \param f Scalar function. + * \param x Scalar input expression. + * \param order Derivation order. + * \return The symbolic derivative of order \p order with respect to \p x. + */ + AnalyticFunction + sympy_diff(const AnalyticFunction& f, const ScalarExpr& x, Index order); + + /** + * \brief Returns the symbolic gradient of a scalar function. + * + * \param f Scalar function. + * \return The symbolic gradient. + */ + AnalyticFunction + sympy_gradient(const AnalyticFunction& f); + + /** + * \brief Returns the symbolic Hessian matrix of a scalar function. + * + * \param f Scalar function. + * \return The symbolic Hessian matrix. + */ + AnalyticFunction + sympy_hessian(const AnalyticFunction& f); + + /** + * \brief Returns the symbolic Jacobian matrix of a vector function. + * + * \param f Vector function. + * \return The symbolic Jacobian matrix. + */ + AnalyticFunction + sympy_diff(const AnalyticFunction& f); + + /** + * \brief Returns a truncated Taylor series of a scalar univariate function. + * + * \param f Scalar function. + * \param center Expansion center. + * \param order Truncation order. + * \return The truncated Taylor series. + */ + AnalyticFunction + sympy_series(const AnalyticFunction& f, double center, Index order); + + /** + * \brief Returns a truncated Taylor series of a scalar function with + * respect to a scalar input variable. + * + * Other input variables are treated as constants. + * + * \param f Scalar function. + * \param x Scalar input variable. + * \param center Expansion center. + * \param order Truncation order. + * \return The truncated Taylor series. + */ + AnalyticFunction + sympy_series(const AnalyticFunction& f, const ScalarVar& x, double center, Index order); + + /** + * \brief Returns a truncated Taylor series of a scalar function with + * respect to a scalar input expression. + * + * Other input variables are treated as constants. This overload supports + * either a scalar input variable, or a direct scalar component of a + * vectorial or matricial input variable. + * + * \param f Scalar function. + * \param x Scalar input expression. + * \param center Expansion center. + * \param order Truncation order. + * \return The truncated Taylor series. + */ + AnalyticFunction + sympy_series(const AnalyticFunction& f, const ScalarExpr& x, double center, Index order); + + /** + * \brief Tests symbolic equality of two scalar analytic functions through SymPy. + * + * The two functions are compared on the same flattened domain: only the + * order and total number of flattened scalar inputs matter, not the original + * grouping into scalar, vector or matrix arguments. + * + * \param f First scalar function. + * \param g Second scalar function. + * \return True if the two functions are symbolically equal. + */ + bool sympy_equal(const AnalyticFunction& f, const AnalyticFunction& g); + + /** + * \brief Tests symbolic equality of two vector analytic functions through SymPy. + * + * Equality is checked componentwise on the same flattened domain. + * + * \param f First vector function. + * \param g Second vector function. + * \return True if the two functions are symbolically equal. + */ + bool sympy_equal(const AnalyticFunction& f, const AnalyticFunction& g); + + /** + * \brief Tests symbolic equality of two matrix analytic functions through SymPy. + * + * Equality is checked entrywise on the same flattened domain. + * + * \param f First matrix function. + * \param g Second matrix function. + * \return True if the two functions are symbolically equal. + */ + bool sympy_equal(const AnalyticFunction& f, const AnalyticFunction& g); +} diff --git a/src/extensions/sympy/codac2_sympy_bridge.cpp b/src/extensions/sympy/codac2_sympy_bridge.cpp new file mode 100644 index 000000000..38799f90e --- /dev/null +++ b/src/extensions/sympy/codac2_sympy_bridge.cpp @@ -0,0 +1,688 @@ +/** + * codac2_sympy_bridge.cpp + * ---------------------------------------------------------------------------- + * \date 2026 + * \author Simon Rohou, Maël Godard + * \copyright Copyright 2026 Codac Team + * \license GNU Lesser General Public License (LGPL) + */ + +#include "codac2_sympy_bridge.h" + +#include +#include +#include +#include +#include + +#include "codac2_assert.h" + +namespace codac2 +{ + namespace + { + using PosNode = AnalyticOperationExpr; + using NegNode = AnalyticOperationExpr; + using AddNode = AnalyticOperationExpr; + using SubNode = AnalyticOperationExpr; + using MulNode = AnalyticOperationExpr; + using DivNode = AnalyticOperationExpr; + using PowNode = AnalyticOperationExpr; + using SinNode = AnalyticOperationExpr; + using CosNode = AnalyticOperationExpr; + using ExpNode = AnalyticOperationExpr; + using LogNode = AnalyticOperationExpr; + using SqrtNode = AnalyticOperationExpr; + using SqrNode = AnalyticOperationExpr; + using TanNode = AnalyticOperationExpr; + using AsinNode = AnalyticOperationExpr; + using AcosNode = AnalyticOperationExpr; + using AtanNode = AnalyticOperationExpr; + using Atan2Node = AnalyticOperationExpr; + using SinhNode = AnalyticOperationExpr; + using CoshNode = AnalyticOperationExpr; + using TanhNode = AnalyticOperationExpr; + using AbsNode = AnalyticOperationExpr; + using SignNode = AnalyticOperationExpr; + using FloorNode = AnalyticOperationExpr; + using CeilNode = AnalyticOperationExpr; + using VecCompNode = AnalyticOperationExpr; + using MatCompNode = AnalyticOperationExpr; + + pybind11::tuple py_args(const pybind11::handle& obj) + { + return pybind11::reinterpret_borrow(obj.attr("args")); + } + + bool py_truth(const pybind11::handle& obj, const char* attr) + { + return pybind11::cast(obj.attr(attr)); + } + + double sympy_number_to_double(const pybind11::handle& obj) + { + pybind11::object borrowed = pybind11::reinterpret_borrow(obj); + pybind11::object as_float = symbolic::detail::import_builtins().attr("float")(borrowed); + return pybind11::cast(as_float); + } + + std::vector> maybe_children_expr_base(const std::shared_ptr& e) + { + if(!e) + return {}; + return e->children_expr_base(); + } + } + + namespace symbolic::detail + { + void ensure_python_runtime() + { + static std::once_flag once; + + if(Py_IsInitialized() != 0) + return; + + std::call_once(once, []() + { + if(Py_IsInitialized() == 0) + pybind11::initialize_interpreter(); + }); + } + + namespace + { + std::unordered_map& imported_modules_cache() + { + // Intentionally leaked: cached Python module objects must not be destroyed + // during C++ static teardown, because their destructors may run after the + // Python runtime has already started shutting down. + static auto* modules = new std::unordered_map(); + return *modules; + } + } + + const pybind11::object& import_module(const char* module_name) + { + ensure_python_runtime(); + pybind11::gil_scoped_acquire gil; + + auto& modules = imported_modules_cache(); + const std::string key(module_name); + auto it = modules.find(key); + if(it == modules.end()) + it = modules.emplace(key, pybind11::module_::import(module_name)).first; + + return it->second; + } + + const pybind11::object& import_sympy() + { + return import_module("sympy"); + } + + const pybind11::object& import_polyfuncs() + { + return import_module("sympy.polys.polyfuncs"); + } + + const pybind11::object& import_builtins() + { + return import_module("builtins"); + } + + pybind11::object normalize_sympy_expr(const pybind11::object& sympy, pybind11::object expr, bool do_expand) + { + expr = expr.attr("rewrite")(sympy.attr("sin")); + expr = expr.attr("rewrite")(sympy.attr("cos")); + expr = expr.attr("rewrite")(sympy.attr("sinh")); + expr = expr.attr("rewrite")(sympy.attr("cosh")); + if(do_expand) + expr = sympy.attr("expand")(expr); + return expr; + } + + ScalarBridgeContext::ScalarBridgeContext(const FunctionArgsList& args) + : _symbols(args), + _exporter(_symbols), + _importer(_symbols) + { + } + + const pybind11::object& ScalarBridgeContext::sympy() const + { + return import_sympy(); + } + + const FlatSymbolTable& ScalarBridgeContext::symbols() const + { + return _symbols; + } + + pybind11::object ScalarBridgeContext::export_scalar(const ScalarExpr& y) const + { + return _exporter.export_scalar(y); + } + + ScalarExpr ScalarBridgeContext::import_scalar(const pybind11::handle& obj) const + { + return _importer.import_scalar(obj); + } + + ScalarExpr ScalarBridgeContext::transform_scalar_expr( + const ScalarExpr& y, + const SympyTransform& transform, + bool do_expand) const + { + ensure_python_runtime(); + pybind11::gil_scoped_acquire gil; + const pybind11::object& sympy = import_sympy(); + pybind11::object ys = _exporter.export_scalar(y); + pybind11::object out = transform(sympy, ys, _symbols); + return _importer.import_scalar(normalize_sympy_expr(sympy, out, do_expand)); + } + + ScalarExpr transform_scalar_expr( + const FunctionArgsList& args, + const ScalarExpr& y, + const SympyTransform& transform, + bool do_expand) + { + ensure_python_runtime(); + pybind11::gil_scoped_acquire gil; + ScalarBridgeContext ctx(args); + return ctx.transform_scalar_expr(y, transform, do_expand); + } + + FlatSymbolTable::FlatSymbolTable(const FunctionArgsList& args) + : _layout(args) + { + _names.reserve(static_cast(_layout.size())); + + for(const auto& arg : args) + { + const auto& b = _layout.binding_of(arg->unique_id()); + + if(auto s = std::dynamic_pointer_cast(arg)) + { + const auto name = make_symbol_name(b.offset); + _names.push_back(name); + _codac_scalars.emplace(name, ScalarExpr(*s)); + continue; + } + + if(auto v = std::dynamic_pointer_cast(arg)) + { + VectorExpr vv(*v); + for(Index i = 0 ; i < v->size() ; ++i) + { + const auto name = make_symbol_name(b.offset + i); + _names.push_back(name); + _codac_scalars.emplace(name, vv[i]); + } + continue; + } + + if(auto m = std::dynamic_pointer_cast(arg)) + { + MatrixExpr mm(*m); + for(Index i = 0 ; i < m->rows() ; ++i) + { + for(Index j = 0 ; j < m->cols() ; ++j) + { + const auto name = make_symbol_name(b.offset + b.cols*i + j); + _names.push_back(name); + _codac_scalars.emplace(name, mm(i,j)); + } + } + continue; + } + + assert_release(false && "Unsupported variable type in FlatSymbolTable"); + } + } + + pybind11::object FlatSymbolTable::by_flat_index(Index k) const + { + assert_release(k >= 0 && k < static_cast(_names.size()) + && "Flat input index out of bounds"); + + ensure_python_runtime(); + pybind11::gil_scoped_acquire gil; + return import_sympy().attr("Symbol")(_names[static_cast(k)], pybind11::arg("real")=true); + } + + ScalarExpr FlatSymbolTable::codac_expr_by_name(const std::string& name) const + { + auto it = _codac_scalars.find(name); + assert_release(it != _codac_scalars.end() && "Unknown SymPy symbol in importer"); + return it->second; + } + + Index FlatSymbolTable::size() const + { + return _layout.size(); + } + + pybind11::object FlatSymbolTable::for_scalar_var(const ScalarVar& x) const + { + const auto& b = binding_of(x.unique_id()); + assert_release(b.is_scalar() && "Internal binding mismatch for scalar variable"); + return by_flat_index(b.offset); + } + + pybind11::object FlatSymbolTable::for_vector_component(const VectorVar& x, Index i) const + { + const auto& b = binding_of(x.unique_id()); + assert_release(b.is_vector() && "Internal binding mismatch for vector variable"); + assert_release(i >= 0 && i < b.rows && "Vector component out of bounds"); + return by_flat_index(b.offset + i); + } + + pybind11::object FlatSymbolTable::for_matrix_component(const MatrixVar& x, Index i, Index j) const + { + const auto& b = binding_of(x.unique_id()); + assert_release(b.is_matrix() && "Internal binding mismatch for matrix variable"); + assert_release(i >= 0 && i < b.rows && j >= 0 && j < b.cols + && "Matrix component out of bounds"); + + return by_flat_index(b.offset + b.cols*i + j); + } + + std::string FlatSymbolTable::make_symbol_name(Index flat_index) + { + std::ostringstream oss; + oss << "_codac_sym_" << flat_index; + return oss.str(); + } + + const FlatInputBinding& FlatSymbolTable::binding_of(const ExprID& id) const + { + return _layout.binding_of(id); + } + + SympyExporter::SympyExporter(const FlatSymbolTable& symbols) + : _symbols(symbols) + { + } + + pybind11::object SympyExporter::export_scalar(const ScalarExpr& e) const + { + ensure_python_runtime(); + pybind11::gil_scoped_acquire gil; + return export_node(std::static_pointer_cast(e)); + } + + pybind11::object SympyExporter::export_node(const std::shared_ptr& e) const + { + const pybind11::object& sympy = import_sympy(); + + if(auto c = std::dynamic_pointer_cast>(e)) + { + const double v = scalar_const_value(*c); + const double r = std::round(v); + if(std::abs(v - r) <= 0.) + return sympy.attr("Integer")(static_cast(r)); + return pybind11::float_(v); + } + + if(auto x = std::dynamic_pointer_cast(e)) + return _symbols.for_scalar_var(*x); + + if(auto op = std::dynamic_pointer_cast(e)) + { + auto ch = op->children_expr_base(); + return export_node(child_at(ch,0)); + } + + if(auto op = std::dynamic_pointer_cast(e)) + { + auto ch = op->children_expr_base(); + return -export_node(child_at(ch,0)); + } + + if(auto op = std::dynamic_pointer_cast(e)) + { + auto ch = op->children_expr_base(); + return export_node(child_at(ch,0)) + export_node(child_at(ch,1)); + } + + if(auto op = std::dynamic_pointer_cast(e)) + { + auto ch = op->children_expr_base(); + return export_node(child_at(ch,0)) - export_node(child_at(ch,1)); + } + + if(auto op = std::dynamic_pointer_cast(e)) + { + auto ch = op->children_expr_base(); + return export_node(child_at(ch,0)) * export_node(child_at(ch,1)); + } + + if(auto op = std::dynamic_pointer_cast(e)) + { + auto ch = op->children_expr_base(); + return export_node(child_at(ch,0)) / export_node(child_at(ch,1)); + } + + if(auto op = std::dynamic_pointer_cast(e)) + { + auto ch = op->children_expr_base(); + return sympy.attr("Pow")(export_node(child_at(ch,0)), export_node(child_at(ch,1))); + } + + if(auto op = std::dynamic_pointer_cast(e)) + { + auto ch = op->children_expr_base(); + return sympy.attr("sin")(export_node(child_at(ch,0))); + } + + if(auto op = std::dynamic_pointer_cast(e)) + { + auto ch = op->children_expr_base(); + return sympy.attr("cos")(export_node(child_at(ch,0))); + } + + if(auto op = std::dynamic_pointer_cast(e)) + { + auto ch = op->children_expr_base(); + return sympy.attr("exp")(export_node(child_at(ch,0))); + } + + if(auto op = std::dynamic_pointer_cast(e)) + { + auto ch = op->children_expr_base(); + return sympy.attr("log")(export_node(child_at(ch,0))); + } + + if(auto op = std::dynamic_pointer_cast(e)) + { + auto ch = op->children_expr_base(); + return sympy.attr("sqrt")(export_node(child_at(ch,0))); + } + + if(auto op = std::dynamic_pointer_cast(e)) + { + auto ch = op->children_expr_base(); + return sympy.attr("Pow")(export_node(child_at(ch,0)), 2); + } + + if(auto op = std::dynamic_pointer_cast(e)) + { + auto ch = op->children_expr_base(); + return sympy.attr("tan")(export_node(child_at(ch,0))); + } + + if(auto op = std::dynamic_pointer_cast(e)) + { + auto ch = op->children_expr_base(); + return sympy.attr("asin")(export_node(child_at(ch,0))); + } + + if(auto op = std::dynamic_pointer_cast(e)) + { + auto ch = op->children_expr_base(); + return sympy.attr("acos")(export_node(child_at(ch,0))); + } + + if(auto op = std::dynamic_pointer_cast(e)) + { + auto ch = op->children_expr_base(); + return sympy.attr("atan")(export_node(child_at(ch,0))); + } + + if(auto op = std::dynamic_pointer_cast(e)) + { + auto ch = op->children_expr_base(); + return sympy.attr("atan2")( + export_node(child_at(ch,0)), export_node(child_at(ch,1))); + } + + if(auto op = std::dynamic_pointer_cast(e)) + { + auto ch = op->children_expr_base(); + return sympy.attr("sinh")(export_node(child_at(ch,0))); + } + + if(auto op = std::dynamic_pointer_cast(e)) + { + auto ch = op->children_expr_base(); + return sympy.attr("cosh")(export_node(child_at(ch,0))); + } + + if(auto op = std::dynamic_pointer_cast(e)) + { + auto ch = op->children_expr_base(); + return sympy.attr("tanh")(export_node(child_at(ch,0))); + } + + if(auto op = std::dynamic_pointer_cast(e)) + { + auto ch = op->children_expr_base(); + return sympy.attr("Abs")(export_node(child_at(ch,0))); + } + + if(auto op = std::dynamic_pointer_cast(e)) + { + auto ch = op->children_expr_base(); + return sympy.attr("sign")(export_node(child_at(ch,0))); + } + + if(auto op = std::dynamic_pointer_cast(e)) + { + auto ch = op->children_expr_base(); + return sympy.attr("floor")(export_node(child_at(ch,0))); + } + + if(auto op = std::dynamic_pointer_cast(e)) + { + auto ch = op->children_expr_base(); + return sympy.attr("ceiling")(export_node(child_at(ch,0))); + } + + if(auto op = std::dynamic_pointer_cast(e)) + { + auto ch = op->children_expr_base(); + return export_vector_component(child_at(ch,0), op->i()); + } + + if(auto op = std::dynamic_pointer_cast(e)) + { + auto ch = op->children_expr_base(); + return export_matrix_component(child_at(ch,0), op->i(), op->j()); + } + + assert_release(false && "Unsupported Codac scalar node in SympyExporter"); + return pybind11::none(); + } + + pybind11::object SympyExporter::export_vector_component(const std::shared_ptr& e, Index i) const + { + if(auto v = std::dynamic_pointer_cast(e)) + return _symbols.for_vector_component(*v, i); + + if(auto ve = std::dynamic_pointer_cast>(e)) + { + auto shape = ve->output_shape(); + auto children = maybe_children_expr_base(e); + if(shape.second == 1 && i >= 0 && i < shape.first && children.size() == static_cast(shape.first)) + return export_node(child_at(children,i)); + } + + assert_release(false + && "Unsupported vector-component expression. Supported cases: direct VectorVar components, or vector expressions exposing scalar children via children_expr_base()."); + return pybind11::none(); + } + + pybind11::object SympyExporter::export_matrix_component(const std::shared_ptr& e, Index i, Index j) const + { + if(auto m = std::dynamic_pointer_cast(e)) + return _symbols.for_matrix_component(*m, i, j); + + if(auto me = std::dynamic_pointer_cast>(e)) + { + auto shape = me->output_shape(); + auto children = maybe_children_expr_base(e); + if(i >= 0 && i < shape.first && j >= 0 && j < shape.second + && children.size() == static_cast(shape.second)) + return export_vector_component(child_at(children,j), i); + } + + assert_release(false + && "Unsupported matrix-component expression. Supported cases: direct MatrixVar components, or matrix expressions exposing column children via children_expr_base()."); + return pybind11::none(); + } + + double SympyExporter::scalar_const_value(const ConstValueExpr& c) + { + const Interval& v = c.value(); + assert_release(v.is_degenerated() && "Only degenerate scalar constants can be exported to SymPy"); + return v.mid(); + } + + std::shared_ptr SympyExporter::child_at( + const std::vector>& children, Index i) + { + assert_release(i >= 0 && i < static_cast(children.size()) + && "Arity mismatch while exporting SymPy expression"); + return children[static_cast(i)]; + } + + SympyImporter::SympyImporter(const FlatSymbolTable& symbols) + : _symbols(symbols) + { + } + + ScalarExpr SympyImporter::import_scalar(const pybind11::handle& obj) const + { + ensure_python_runtime(); + pybind11::gil_scoped_acquire gil; + + if(py_truth(obj, "is_number")) + return const_value(sympy_number_to_double(obj)); + + if(py_truth(obj, "is_Symbol")) + return _symbols.codac_expr_by_name(pybind11::cast(pybind11::str(obj))); + + if(py_truth(obj, "is_Add")) + return import_add(obj); + + if(py_truth(obj, "is_Mul")) + return import_mul(obj); + + if(py_truth(obj, "is_Pow")) + return import_pow(obj); + + return import_function(obj); + } + + ScalarExpr SympyImporter::import_add(const pybind11::handle& obj) const + { + pybind11::tuple args = py_args(obj); + assert_release(!args.empty() && "Unexpected empty Add in SymPy importer"); + + ScalarExpr acc = import_scalar(args[0]); + for(pybind11::size_t i = 1 ; i < args.size() ; ++i) + acc = acc + import_scalar(args[i]); + return acc; + } + + ScalarExpr SympyImporter::import_mul(const pybind11::handle& obj) const + { + pybind11::tuple args = py_args(obj); + assert_release(!args.empty() && "Unexpected empty Mul in SymPy importer"); + + ScalarExpr acc = import_scalar(args[0]); + for(pybind11::size_t i = 1 ; i < args.size() ; ++i) + acc = acc * import_scalar(args[i]); + return acc; + } + + ScalarExpr SympyImporter::import_pow(const pybind11::handle& obj) const + { + pybind11::tuple args = py_args(obj); + assert_release(args.size() == 2 && "SymPy Pow node should have arity 2"); + + const pybind11::object& sympy = import_sympy(); + const auto base = import_scalar(args[0]); + pybind11::handle expo = args[1]; + + if(pybind11::isinstance(expo, sympy.attr("Integer"))) + { + const long n = pybind11::cast(expo); + + if(n == 0) return const_value(1.); + if(n == 1) return base; + if(n == 2) return sqr(base); + if(n == -1) return const_value(1.) / base; + if(n == -2) return const_value(1.) / sqr(base); + + ScalarExpr acc = const_value(1.); + if(n > 0) + { + for(long k = 0 ; k < n ; ++k) + acc = acc * base; + return acc; + } + else + { + for(long k = 0 ; k < -n ; ++k) + acc = acc * base; + return const_value(1.) / acc; + } + } + + if(pybind11::isinstance(expo, sympy.attr("Rational"))) + { + const long num = pybind11::cast(expo.attr("p")); + const long den = pybind11::cast(expo.attr("q")); + + if(num == 1 && den == 2) + return sqrt(base); + if(num == -1 && den == 2) + return const_value(1.) / sqrt(base); + } + + return pow(base, import_scalar(expo)); + } + + ScalarExpr SympyImporter::import_function(const pybind11::handle& obj) const + { + const pybind11::object& sympy = import_sympy(); + pybind11::object func = obj.attr("func"); + pybind11::tuple args = py_args(obj); + + if(args.size() == 1) + { + auto x = import_scalar(args[0]); + + if(func.is(sympy.attr("sin"))) return sin(x); + if(func.is(sympy.attr("cos"))) return cos(x); + if(func.is(sympy.attr("tan"))) return tan(x); + if(func.is(sympy.attr("asin"))) return asin(x); + if(func.is(sympy.attr("acos"))) return acos(x); + if(func.is(sympy.attr("atan"))) return atan(x); + if(func.is(sympy.attr("sinh"))) return sinh(x); + if(func.is(sympy.attr("cosh"))) return cosh(x); + if(func.is(sympy.attr("tanh"))) return tanh(x); + if(func.is(sympy.attr("exp"))) return exp(x); + if(func.is(sympy.attr("log"))) return log(x); + if(func.is(sympy.attr("sqrt"))) return sqrt(x); + if(func.is(sympy.attr("Abs"))) return abs(x); + if(func.is(sympy.attr("sign"))) return sign(x); + if(func.is(sympy.attr("floor"))) return floor(x); + if(func.is(sympy.attr("ceiling"))) return ceil(x); + } + else if(args.size() == 2) + { + auto x1 = import_scalar(args[0]); + auto x2 = import_scalar(args[1]); + + if(func.is(sympy.attr("atan2"))) return atan2(x1, x2); + } + + assert_release(false && "Unsupported SymPy node in importer"); + return const_value(0.); + } + } +} diff --git a/src/extensions/sympy/codac2_sympy_bridge.h b/src/extensions/sympy/codac2_sympy_bridge.h new file mode 100644 index 000000000..841061fd5 --- /dev/null +++ b/src/extensions/sympy/codac2_sympy_bridge.h @@ -0,0 +1,303 @@ +/** + * \file codac2_sympy_bridge.h + * ---------------------------------------------------------------------------- + * \date 2026 + * \author Simon Rohou, Maël Godard + * \copyright Copyright 2026 Codac Team + * \license GNU Lesser General Public License (LGPL) + */ + +#pragma once + +#include +#include + +#include "codac2_sympy.h" +#include "codac2_analytic_flat_input_layout.h" + +namespace codac2::symbolic::detail +{ + /** + * \brief Symbol table associated with a flattened analytic input domain. + * + * This class maps each scalar entry of a function input domain to a dedicated + * SymPy symbol, and provides the reverse mapping from such symbols to Codac + * scalar expressions. + */ + class FlatSymbolTable + { + public: + + /** + * \brief Builds the symbol table associated with a function argument list. + * + * \param args Function argument list. + */ + explicit FlatSymbolTable(const FunctionArgsList& args); + + /** + * \brief Returns the SymPy symbol associated with a flat input index. + * + * \param k Flat input index. + * \return SymPy symbol associated with \p k. + */ + pybind11::object by_flat_index(Index k) const; + + /** + * \brief Returns the Codac scalar expression associated with a symbol name. + * + * \param name SymPy symbol name. + * \return Corresponding Codac scalar expression. + */ + ScalarExpr codac_expr_by_name(const std::string& name) const; + + /** + * \brief Returns the number of scalar entries in the flattened input domain. + * + * \return Flattened input domain size. + */ + Index size() const; + + /** + * \brief Returns the SymPy symbol associated with a scalar variable. + * + * \param x Scalar variable. + * \return Corresponding SymPy symbol. + */ + pybind11::object for_scalar_var(const ScalarVar& x) const; + + /** + * \brief Returns the SymPy symbol associated with a vector component. + * + * \param x Vector variable. + * \param i Component index. + * \return Corresponding SymPy symbol. + */ + pybind11::object for_vector_component(const VectorVar& x, Index i) const; + + /** + * \brief Returns the SymPy symbol associated with a matrix component. + * + * \param x Matrix variable. + * \param i Row index. + * \param j Column index. + * \return Corresponding SymPy symbol. + */ + pybind11::object for_matrix_component(const MatrixVar& x, Index i, Index j) const; + + private: + + /** + * \brief Builds the symbol name associated with a flat input index. + * + * \param flat_index Flat input index. + * \return Symbol name associated with \p flat_index. + */ + static std::string make_symbol_name(Index flat_index); + + /** + * \brief Returns the flattened binding associated with an input identifier. + * + * \param id Input expression identifier. + * \return Binding associated with \p id. + */ + const FlatInputBinding& binding_of(const ExprID& id) const; + + FlatInputLayout _layout; //!< Flattened input layout. + std::vector _names; //!< Symbol names indexed by flat input index. + std::unordered_map _codac_scalars; //!< Reverse mapping from symbol names to Codac scalar expressions. + }; + + /** + * \brief Exports Codac scalar expressions to SymPy expressions. + */ + class SympyExporter + { + public: + + /** + * \brief Builds an exporter associated with a symbol table. + * + * \param symbols Flat symbol table. + */ + explicit SympyExporter(const FlatSymbolTable& symbols); + + /** + * \brief Exports a Codac scalar expression to SymPy. + * + * \param e Codac scalar expression. + * \return Corresponding SymPy expression. + */ + pybind11::object export_scalar(const ScalarExpr& e) const; + + private: + + pybind11::object export_node(const std::shared_ptr& e) const; + pybind11::object export_vector_component(const std::shared_ptr& e, Index i) const; + pybind11::object export_matrix_component(const std::shared_ptr& e, Index i, Index j) const; + static double scalar_const_value(const ConstValueExpr& c); + static std::shared_ptr child_at(const std::vector>& children, Index i); + + const FlatSymbolTable& _symbols; //!< Associated symbol table. + }; + + /** + * \brief Imports SymPy scalar expressions into Codac scalar expressions. + */ + class SympyImporter + { + public: + + /** + * \brief Builds an importer associated with a symbol table. + * + * \param symbols Flat symbol table. + */ + explicit SympyImporter(const FlatSymbolTable& symbols); + + /** + * \brief Imports a SymPy scalar expression. + * + * \param obj SymPy expression handle. + * \return Corresponding Codac scalar expression. + */ + ScalarExpr import_scalar(const pybind11::handle& obj) const; + + private: + + ScalarExpr import_add(const pybind11::handle& obj) const; + ScalarExpr import_mul(const pybind11::handle& obj) const; + ScalarExpr import_pow(const pybind11::handle& obj) const; + ScalarExpr import_function(const pybind11::handle& obj) const; + + const FlatSymbolTable& _symbols; //!< Associated symbol table. + }; + + /** + * \brief SymPy transformation applied to an exported scalar expression. + */ + using SympyTransform = std::function; + + /** + * \brief Context gathering the bridge objects needed for one scalar transformation. + */ + class ScalarBridgeContext + { + public: + + /** + * \brief Builds a scalar bridge context for a given function argument list. + * + * \param args Function argument list. + */ + explicit ScalarBridgeContext(const FunctionArgsList& args); + + /** + * \brief Returns the imported SymPy module. + * + * \return SymPy module. + */ + const pybind11::object& sympy() const; + + /** + * \brief Returns the associated flat symbol table. + * + * \return Flat symbol table. + */ + const FlatSymbolTable& symbols() const; + + /** + * \brief Exports a Codac scalar expression to SymPy. + * + * \param y Codac scalar expression. + * \return Corresponding SymPy expression. + */ + pybind11::object export_scalar(const ScalarExpr& y) const; + + /** + * \brief Imports a SymPy scalar expression into Codac. + * + * \param obj SymPy expression handle. + * \return Corresponding Codac scalar expression. + */ + ScalarExpr import_scalar(const pybind11::handle& obj) const; + + /** + * \brief Applies a SymPy transformation to a Codac scalar expression. + * + * \param y Codac scalar expression. + * \param transform SymPy transformation. + * \param do_expand Whether the transformed expression should be expanded. + * \return Transformed Codac scalar expression. + */ + ScalarExpr transform_scalar_expr(const ScalarExpr& y, const SympyTransform& transform, bool do_expand = true) const; + + private: + + FlatSymbolTable _symbols; //!< Flat symbol table. + SympyExporter _exporter; //!< Exporter. + SympyImporter _importer; //!< Importer. + }; + + /** + * \brief Ensures that the embedded Python runtime is initialized. + */ + void ensure_python_runtime(); + + /** + * \brief Imports a Python module. + * + * Imported modules are cached after the first call. + * + * \param module_name Python module name. + * \return Imported Python module. + */ + const pybind11::object& import_module(const char* module_name); + + /** + * \brief Imports the SymPy module. + * + * \return SymPy module. + */ + const pybind11::object& import_sympy(); + + /** + * \brief Imports the SymPy polyfuncs module. + * + * \return SymPy polyfuncs module. + */ + const pybind11::object& import_polyfuncs(); + + /** + * \brief Imports the Python builtins module. + * + * \return Python builtins module. + */ + const pybind11::object& import_builtins(); + + /** + * \brief Normalizes a SymPy expression before import into Codac. + * + * \param sympy SymPy module. + * \param expr SymPy expression. + * \param do_expand Whether the expression should be expanded. + * \return Normalized SymPy expression. + */ + pybind11::object normalize_sympy_expr(const pybind11::object& sympy, pybind11::object expr, bool do_expand = true); + + /** + * \brief Applies a SymPy transformation to a Codac scalar expression. + * + * \param args Function argument list. + * \param y Codac scalar expression. + * \param transform SymPy transformation. + * \param do_expand Whether the transformed expression should be expanded. + * \return Transformed Codac scalar expression. + */ + ScalarExpr transform_scalar_expr( + const FunctionArgsList& args, + const ScalarExpr& y, + const SympyTransform& transform, + bool do_expand = true); +} diff --git a/src/extensions/sympy/codac2_sympy_empty.cpp b/src/extensions/sympy/codac2_sympy_empty.cpp deleted file mode 100644 index 621b6026f..000000000 --- a/src/extensions/sympy/codac2_sympy_empty.cpp +++ /dev/null @@ -1,18 +0,0 @@ -/** - * \file codac2_sympy_empty.cpp - * ---------------------------------------------------------------------------- - * \date 2024 - * \author Simon Rohou - * \copyright Copyright 2024 Codac Team - * \license GNU Lesser General Public License (LGPL) - */ - -#include "codac2_sympy_empty.h" - -namespace codac2 -{ - namespace sympy - { - - } -} \ No newline at end of file diff --git a/src/extensions/sympy/codac2_sympy_empty.h b/src/extensions/sympy/codac2_sympy_empty.h deleted file mode 100644 index df1ebd334..000000000 --- a/src/extensions/sympy/codac2_sympy_empty.h +++ /dev/null @@ -1,18 +0,0 @@ -/** - * \file codac2_sympy_empty.h - * ---------------------------------------------------------------------------- - * \date 2024 - * \author Simon Rohou - * \copyright Copyright 2024 Codac Team - * \license GNU Lesser General Public License (LGPL) - */ - -#pragma once - -namespace codac2 -{ - namespace sympy - { - - } -} \ No newline at end of file diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index 38b8f5210..b3eded7af 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -116,6 +116,20 @@ file(COPY set(CODAC_LIBRARIES ${PROJECT_NAME}-core ${PROJECT_NAME}-graphics) +option(BUILD_SYMPY_EMBED_TESTS "Build C++ tests that require embedded Python for Sympy" ON) +if(WITH_PYTHON AND DEFINED PYBIND11_FINDPYTHON AND NOT PYBIND11_FINDPYTHON) + message(STATUS "Disabling Sympy C++ embed tests because PYBIND11_FINDPYTHON=OFF") + set(BUILD_SYMPY_EMBED_TESTS OFF) +endif() + +# Sympy test +if (WITH_PYTHON AND BUILD_SYMPY_EMBED_TESTS) + list(APPEND SRC_TESTS + extensions/sympy/codac2_tests_sympy + ../doc/manual/manual/extensions/sympy/src + ) +endif() + # CAPD test if (WITH_CAPD) list(APPEND SRC_TESTS @@ -142,6 +156,10 @@ foreach(SRC_TEST ${SRC_TESTS}) set(CODAC_HEADERS_DIR ${CMAKE_CURRENT_BINARY_DIR}/../include) target_include_directories(${TEST_NAME} SYSTEM PUBLIC ${CODAC_HEADERS_DIR}) target_link_libraries(${TEST_NAME} PUBLIC Ibex::ibex ${CODAC_LIBRARIES} PRIVATE Catch2::Catch2WithMain) + + if(WITH_PYTHON AND SRC_TEST MATCHES "extensions/sympy/") + target_link_libraries(${TEST_NAME} PRIVATE ${PROJECT_NAME}-sympy pybind11::embed) + endif() add_dependencies(check ${TEST_NAME}) add_test(NAME ${TEST_NAME}_cpp COMMAND ${TEST_NAME}) diff --git a/tests/extensions/sympy/codac2_tests_sympy.cpp b/tests/extensions/sympy/codac2_tests_sympy.cpp new file mode 100644 index 000000000..d7231a327 --- /dev/null +++ b/tests/extensions/sympy/codac2_tests_sympy.cpp @@ -0,0 +1,271 @@ +/** + * Codac tests + * ---------------------------------------------------------------------------- + * \date 2026 + * \author Simon Rohou + * \copyright Copyright 2026 Codac Team + * \license GNU Lesser General Public License (LGPL) + */ + +#include +#include + +using namespace codac2; + +TEST_CASE("Sympy") +{ + { + ScalarVar x; + AnalyticFunction f({x}, cos(x)*x + sin(x)); + auto df = sympy_diff(f); + auto dfx = sympy_diff(f, x); + + CHECK(sympy_equal(df, AnalyticFunction({x}, 2.*cos(x) - x*sin(x)))); + CHECK(sympy_equal(dfx, AnalyticFunction({x}, 2.*cos(x) - x*sin(x)))); + } + + { + ScalarVar x; + ScalarVar y; + AnalyticFunction f({x,y}, x*y + sin(x)); + + auto dfdx = sympy_partial_diff(f, x); + auto dfdy = sympy_partial_diff(f, y); + + CHECK(sympy_equal(dfdx, AnalyticFunction({x,y}, y + cos(x)))); + CHECK(sympy_equal(dfdy, AnalyticFunction({x,y}, x))); + } + + { + VectorVar v(2); + AnalyticFunction g({v}, v[0]*v[1] + sin(v[0])); + + auto dg_dv0 = sympy_partial_diff(g, v[0]); + auto dg_dv1 = sympy_partial_diff(g, v[1]); + auto grad_g = sympy_gradient(g); + + CHECK(sympy_equal(dg_dv0, AnalyticFunction({v}, v[1] + cos(v[0])))); + CHECK(sympy_equal(dg_dv1, AnalyticFunction({v}, v[0]))); + CHECK(sympy_equal(grad_g, AnalyticFunction({v}, vec(v[1] + cos(v[0]), v[0])))); + } + + { + VectorVar v(2); + AnalyticFunction f({v}, { + v[0]*v[1] + sin(v[0]), + sqr(v[0]) + cos(v[1]) + }); + + auto J = sympy_diff(f); + + CHECK(sympy_equal( + J, + AnalyticFunction( + {v}, + mat( + vec(v[1] + cos(v[0]), 2.*v[0]), + vec(v[0], -sin(v[1])) + )))); + } + + { + ScalarVar x; + AnalyticFunction f({x}, sqr(sin(x)) + sqr(cos(x))); + auto fs = sympy_simplify(f); + auto dfs = sympy_diff(fs); + + CHECK(sympy_equal(fs, AnalyticFunction({x}, 1.))); + CHECK(sympy_equal(dfs, AnalyticFunction({x}, 0.))); + } + + { + VectorVar v(2); + AnalyticFunction f({v}, { + v[0] + v[0], + sqr(sin(v[1])) + sqr(cos(v[1])) + }); + + auto fs = sympy_simplify(f); + + CHECK(sympy_equal(fs, AnalyticFunction({v}, { + 2.*v[0], + 1. + }))); + } + + { + ScalarVar x; + AnalyticFunction f({x}, sin(x)); + + auto d3f = sympy_diff(f, x, 3); + auto d0f = sympy_diff(f, x, 0); + + CHECK(sympy_equal(d3f, AnalyticFunction({x}, -cos(x)))); + CHECK(sympy_equal(d0f, f)); + } + + { + ScalarVar x; + ScalarVar y; + AnalyticFunction f({x,y}, x*y + sin(x) + sqr(y)); + + auto H = sympy_hessian(f); + + CHECK(sympy_equal( + H, + AnalyticFunction( + {x,y}, + mat( + vec(-sin(x), 1.), + vec(1., 2.) + )))); + } + + { + ScalarVar x; + AnalyticFunction f({x}, 1. / (1. - x)); + auto p = sympy_series(f, x, 0.0, 3); + + CHECK(sympy_equal(p, AnalyticFunction({x}, 1. + x + sqr(x) + x*sqr(x)))); + } + + { + ScalarVar x; + ScalarVar y; + AnalyticFunction f({x,y}, y + 1. / (1. - x)); + auto p = sympy_series(f, x, 0.0, 2); + + CHECK(sympy_equal(p, AnalyticFunction({x,y}, y + 1. + x + sqr(x)))); + } + + { + ScalarVar x; + AnalyticFunction f({x}, 9.*pow(x,4) + 8.*pow(x,3) + 7.*sqr(x) + 6.*x + 5.); + auto h = sympy_horner(f); + + CHECK(sympy_equal(h, f)); + } + + { + VectorVar v(2); + AnalyticFunction f({v}, { + 3.*pow(v[0],3) + 2.*sqr(v[0]) + v[0] + 1., + v[1]*(v[1] + 1.) + 2. + }); + + auto h = sympy_horner(f); + + CHECK(sympy_equal(h, f)); + } + + { + MatrixVar A(2,2); + AnalyticFunction f({A}, mat( + vec(A(0,0)*A(0,0) + 2.*A(0,0)*A(0,1) + 1., A(1,0) + 3.*A(1,1)), + vec(A(0,0) - A(0,1), A(1,1)*A(1,1) + A(1,1) + 2.) + )); + auto h = sympy_horner(f); + + CHECK(sympy_equal(h, f)); + } + + { + ScalarVar x; + + auto dtan = sympy_diff(AnalyticFunction({x}, tan(x))); + auto dasin = sympy_diff(AnalyticFunction({x}, asin(x))); + auto dacos = sympy_diff(AnalyticFunction({x}, acos(x))); + auto datan = sympy_diff(AnalyticFunction({x}, atan(x))); + auto dsinh = sympy_diff(AnalyticFunction({x}, sinh(x))); + auto dcosh = sympy_diff(AnalyticFunction({x}, cosh(x))); + auto dtanh = sympy_diff(AnalyticFunction({x}, tanh(x))); + + CHECK(sympy_equal(dtan, AnalyticFunction({x}, 1. / sqr(cos(x))))); + CHECK(sympy_equal(dasin, AnalyticFunction({x}, 1. / sqrt(1. - sqr(x))))); + CHECK(sympy_equal(dacos, AnalyticFunction({x}, -1. / sqrt(1. - sqr(x))))); + CHECK(sympy_equal(datan, AnalyticFunction({x}, 1. / (1. + sqr(x))))); + CHECK(sympy_equal(dsinh, AnalyticFunction({x}, cosh(x)))); + CHECK(sympy_equal(dcosh, AnalyticFunction({x}, sinh(x)))); + CHECK(sympy_equal(dtanh, AnalyticFunction({x}, 1. - sqr(tanh(x))))); + } + + { + ScalarVar x; + AnalyticFunction f({x}, abs(x)); + auto df = sympy_diff(f); + + CHECK(sympy_equal(df, AnalyticFunction({x}, sign(x)))); + } + + { + ScalarVar y; + ScalarVar x; + AnalyticFunction f({y,x}, atan2(y, x)); + + auto dfdy = sympy_partial_diff(f, y); + auto dfdx = sympy_partial_diff(f, x); + + CHECK(sympy_equal(dfdy, AnalyticFunction({y,x}, x / (sqr(x) + sqr(y))))); + CHECK(sympy_equal(dfdx, AnalyticFunction({y,x}, -y / (sqr(x) + sqr(y))))); + } + + + { + MatrixVar A(2,2); + AnalyticFunction f({A}, A(0,0)*A(1,1) + sin(A(0,1)) + A(1,0)); + + auto dfdA00 = sympy_partial_diff(f, A(0,0)); + auto dfdA01 = sympy_partial_diff(f, A(0,1)); + auto dfdA10 = sympy_partial_diff(f, A(1,0)); + auto dfdA11 = sympy_partial_diff(f, A(1,1)); + auto grad_f = sympy_gradient(f); + + CHECK(sympy_equal(dfdA00, AnalyticFunction({A}, A(1,1)))); + CHECK(sympy_equal(dfdA01, AnalyticFunction({A}, cos(A(0,1))))); + CHECK(sympy_equal(dfdA10, AnalyticFunction({A}, 1.))); + CHECK(sympy_equal(dfdA11, AnalyticFunction({A}, A(0,0)))); + CHECK(sympy_equal( + grad_f, + AnalyticFunction({A}, vec(A(1,1), cos(A(0,1)), 1., A(0,0))))); + } + + { + MatrixVar A(2,2); + AnalyticFunction f({A}, { + A(0,0) + A(1,1), + A(0,1)*A(1,0) + }); + + auto J = sympy_diff(f); + + CHECK(sympy_equal( + J, + AnalyticFunction( + {A}, + mat( + vec(1., 0.), + vec(0., A(1,0)), + vec(0., A(0,1)), + vec(1., 0.))))); + } + + { + ScalarVar x; + ScalarVar y; + VectorVar v(2); + MatrixVar A(2,1); + + CHECK(sympy_equal( + AnalyticFunction({x,y}, x + 2.*y), + AnalyticFunction({v}, v[0] + 2.*v[1]))); + + CHECK(sympy_equal( + AnalyticFunction({x,y}, vec(x+y, x-y)), + AnalyticFunction({v}, vec(v[0] + v[1], v[0] - v[1])))); + + CHECK(sympy_equal( + AnalyticFunction({x,y}, mat(vec(x+y, x-y))), + AnalyticFunction({A}, mat(vec(A(0,0) + A(1,0), A(0,0) - A(1,0)))))); + } + +} diff --git a/tests/extensions/sympy/codac2_tests_sympy.py b/tests/extensions/sympy/codac2_tests_sympy.py new file mode 100644 index 000000000..930839674 --- /dev/null +++ b/tests/extensions/sympy/codac2_tests_sympy.py @@ -0,0 +1,234 @@ +#!/usr/bin/env python + +# Codac tests +# ---------------------------------------------------------------------------- +# \date 2026 +# \author Simon Rohou +# \copyright Copyright 2026 Codac Team +# \license GNU Lesser General Public License (LGPL) + +import unittest +from codac import * + + +class TestSympy(unittest.TestCase): + + def test_sympy(self): + x = ScalarVar() + f = AnalyticFunction([x], cos(x)*x + sin(x)) + df = sympy_diff(f) + dfx = sympy_diff(f, x) + + self.assertTrue(sympy_equal(df, AnalyticFunction([x], 2.*cos(x) - x*sin(x)))) + self.assertTrue(sympy_equal(dfx, AnalyticFunction([x], 2.*cos(x) - x*sin(x)))) + + x = ScalarVar() + y = ScalarVar() + f = AnalyticFunction([x,y], x*y + sin(x)) + + dfdx = sympy_partial_diff(f, x) + dfdy = sympy_partial_diff(f, y) + + self.assertTrue(sympy_equal(dfdx, AnalyticFunction([x,y], y + cos(x)))) + self.assertTrue(sympy_equal(dfdy, AnalyticFunction([x,y], x))) + + v = VectorVar(2) + g = AnalyticFunction([v], v[0]*v[1] + sin(v[0])) + + dg_dv0 = sympy_partial_diff(g, v[0]) + dg_dv1 = sympy_partial_diff(g, v[1]) + grad_g = sympy_gradient(g) + + self.assertTrue(sympy_equal(dg_dv0, AnalyticFunction([v], v[1] + cos(v[0])))) + self.assertTrue(sympy_equal(dg_dv1, AnalyticFunction([v], v[0]))) + self.assertTrue(sympy_equal(grad_g, AnalyticFunction([v], vec(v[1] + cos(v[0]), v[0])))) + + v = VectorVar(2) + f = AnalyticFunction([v], [ + v[0]*v[1] + sin(v[0]), + sqr(v[0]) + cos(v[1]) + ]) + + J = sympy_diff(f) + + self.assertTrue(sympy_equal( + J, + AnalyticFunction( + [v], + mat( + vec(v[1] + cos(v[0]), 2.*v[0]), + vec(v[0], -sin(v[1])) + )))) + + x = ScalarVar() + f = AnalyticFunction([x], sqr(sin(x)) + sqr(cos(x))) + fs = sympy_simplify(f) + dfs = sympy_diff(fs) + + self.assertTrue(sympy_equal(fs, AnalyticFunction([x], 1.))) + self.assertTrue(sympy_equal(dfs, AnalyticFunction([x], 0.*x))) + + v = VectorVar(2) + f = AnalyticFunction([v], [ + v[0] + v[0], + sqr(sin(v[1])) + sqr(cos(v[1])) + ]) + + fs = sympy_simplify(f) + + self.assertTrue(sympy_equal(fs, AnalyticFunction([v], [ + 2.*v[0], + 1. + ]))) + + x = ScalarVar() + f = AnalyticFunction([x], sin(x)) + + d3f = sympy_diff(f, x, 3) + d0f = sympy_diff(f, x, 0) + + self.assertTrue(sympy_equal(d3f, AnalyticFunction([x], -cos(x)))) + self.assertTrue(sympy_equal(d0f, f)) + + x = ScalarVar() + y = ScalarVar() + f = AnalyticFunction([x,y], x*y + sin(x) + sqr(y)) + + H = sympy_hessian(f) + + self.assertTrue(sympy_equal( + H, + AnalyticFunction( + [x,y], + mat( + vec(-sin(x), 1.), + vec(1., 2.) + )))) + + x = ScalarVar() + f = AnalyticFunction([x], 1. / (1. - x)) + p = sympy_series(f, x, 0.0, 3) + + self.assertTrue(sympy_equal(p, AnalyticFunction([x], 1. + x + sqr(x) + x*sqr(x)))) + + x = ScalarVar() + y = ScalarVar() + f = AnalyticFunction([x,y], y + 1. / (1. - x)) + p = sympy_series(f, x, 0.0, 2) + + self.assertTrue(sympy_equal(p, AnalyticFunction([x,y], y + 1. + x + sqr(x)))) + + x = ScalarVar() + f = AnalyticFunction([x], 9.*pow(x,4) + 8.*pow(x,3) + 7.*sqr(x) + 6.*x + 5.) + h = sympy_horner(f) + + self.assertTrue(sympy_equal(h, f)) + + v = VectorVar(2) + f = AnalyticFunction([v], [ + 3.*pow(v[0],3) + 2.*sqr(v[0]) + v[0] + 1., + v[1]*(v[1] + 1.) + 2. + ]) + + h = sympy_horner(f) + + self.assertTrue(sympy_equal(h, f)) + + A = MatrixVar(2,2) + f = AnalyticFunction([A], mat( + vec(A(0,0)*A(0,0) + 2.*A(0,0)*A(0,1) + 1., A(1,0) + 3.*A(1,1)), + vec(A(0,0) - A(0,1), A(1,1)*A(1,1) + A(1,1) + 2.) + )) + h = sympy_horner(f) + + self.assertTrue(sympy_equal(h, f)) + + x = ScalarVar() + + dtan = sympy_diff(AnalyticFunction([x], tan(x))) + dasin = sympy_diff(AnalyticFunction([x], asin(x))) + dacos = sympy_diff(AnalyticFunction([x], acos(x))) + datan = sympy_diff(AnalyticFunction([x], atan(x))) + dsinh = sympy_diff(AnalyticFunction([x], sinh(x))) + dcosh = sympy_diff(AnalyticFunction([x], cosh(x))) + dtanh = sympy_diff(AnalyticFunction([x], tanh(x))) + + self.assertTrue(sympy_equal(dtan, AnalyticFunction([x], 1. / sqr(cos(x))))) + self.assertTrue(sympy_equal(dasin, AnalyticFunction([x], 1. / sqrt(1. - sqr(x))))) + self.assertTrue(sympy_equal(dacos, AnalyticFunction([x], -1. / sqrt(1. - sqr(x))))) + self.assertTrue(sympy_equal(datan, AnalyticFunction([x], 1. / (1. + sqr(x))))) + self.assertTrue(sympy_equal(dsinh, AnalyticFunction([x], cosh(x)))) + self.assertTrue(sympy_equal(dcosh, AnalyticFunction([x], sinh(x)))) + self.assertTrue(sympy_equal(dtanh, AnalyticFunction([x], 1. - sqr(tanh(x))))) + + x = ScalarVar() + f = AnalyticFunction([x], abs(x)) + df = sympy_diff(f) + + self.assertTrue(sympy_equal(df, AnalyticFunction([x], sign(x)))) + + y = ScalarVar() + x = ScalarVar() + f = AnalyticFunction([y,x], atan2(y, x)) + + dfdy = sympy_partial_diff(f, y) + dfdx = sympy_partial_diff(f, x) + + self.assertTrue(sympy_equal(dfdy, AnalyticFunction([y,x], x / (sqr(x) + sqr(y))))) + self.assertTrue(sympy_equal(dfdx, AnalyticFunction([y,x], -y / (sqr(x) + sqr(y))))) + + A = MatrixVar(2,2) + f = AnalyticFunction([A], A(0,0)*A(1,1) + sin(A(0,1)) + A(1,0)) + + dfdA00 = sympy_partial_diff(f, A(0,0)) + dfdA01 = sympy_partial_diff(f, A(0,1)) + dfdA10 = sympy_partial_diff(f, A(1,0)) + dfdA11 = sympy_partial_diff(f, A(1,1)) + grad_f = sympy_gradient(f) + + self.assertTrue(sympy_equal(dfdA00, AnalyticFunction([A], A(1,1)))) + self.assertTrue(sympy_equal(dfdA01, AnalyticFunction([A], cos(A(0,1))))) + self.assertTrue(sympy_equal(dfdA10, AnalyticFunction([A], 1.))) + self.assertTrue(sympy_equal(dfdA11, AnalyticFunction([A], A(0,0)))) + self.assertTrue(sympy_equal( + grad_f, + AnalyticFunction([A], vec(A(1,1), cos(A(0,1)), 1., A(0,0))))) + + A = MatrixVar(2,2) + f = AnalyticFunction([A], [ + A(0,0) + A(1,1), + A(0,1)*A(1,0) + ]) + + J = sympy_diff(f) + + self.assertTrue(sympy_equal( + J, + AnalyticFunction( + [A], + mat( + vec(1., 0.), + vec(0., A(1,0)), + vec(0., A(0,1)), + vec(1., 0.) + )))) + + x = ScalarVar() + y = ScalarVar() + v = VectorVar(2) + A = MatrixVar(2,1) + + self.assertTrue(sympy_equal( + AnalyticFunction([x,y], x + 2.*y), + AnalyticFunction([v], v[0] + 2.*v[1]))) + + self.assertTrue(sympy_equal( + AnalyticFunction([x,y], vec(x+y, x-y)), + AnalyticFunction([v], vec(v[0] + v[1], v[0] - v[1])))) + + self.assertTrue(sympy_equal( + AnalyticFunction([x,y], mat(vec(x+y, x-y))), + AnalyticFunction([A], mat(vec(A(0,0) + A(1,0), A(0,0) - A(1,0)))))) + +if __name__ == '__main__': + unittest.main()