[Feature] Composition of compatible DerivativeStructures

[Serrof]

A univariate DerivativeStructure can be composed on the right by any DerivativeStructure of the same maximum order via the existing "compose" method acting on their respective attributes "data". It would be useful to have a generalization of this where the left-hand side is a multivariate DerivativeStructure and the right-hand side is an array-like collection of DerivativeStructures (at least with the same number of variables). This is handy when functions are only known via their partial derivatives (in particular, no algebraic expression).

This is related to issue #190 as one needs the composition operator in order to compute the composition inverse.

[Serrof]

Heads up for an upcoming merge request (in the coming days) on a tentative fix for this issue

[maisonobe]

I am currently working on it. Did you start it also on your side?

[maisonobe]

The intended API is the following one:

    /** Rebase instance with respect to low level parameter functions.
     * <p>
     * The instance is considered to be a function of {@link #getFreeParameters()
     * n free parameters} up to order {@link #getOrder() o} \(f(p_0, p_1, \ldots p_{n-1})\).
     * Its {@link #getPartialDerivative(int...) partial derivatives} are therefore
     * \(f, \frac{\partial f}{\partial p_0}, \frac{\partial f}{\partial p_1}, \ldots
     * \frac{\partial^2 f}{\partial p_0^2}, \frac{\partial^2 f}{\partial p_0 p_1},
     * \ldots \frac{\partial^o f}{\partial p_{n-1}^o\). The free parameters
     * \(p_0, p_1, \ldots p_{n-1}\) are considered to be functions of \(m\) lower
     * level other parameters \(q_0, q_1, \ldots q_{m-1}\).
     * </p>
     * \( \begin{align}
     * p_0 &amp; = p_0(q_0, q_1, \ldots q_{m-1})\\
     * p_1 &amp; = p_1(q_0, q_1, \ldots q_{m-1})\\
     * p_{n-1} &amp; = p_{n-1}(q_0, q_1, \ldots q_{m-1})
     * \end{align}\)
     * <p>
     * This method compute the composition of the partial derivatives of \(f\)
     * and the partial derivatives of \(p_0, p_1, \ldots p_{n-1}, i.e. the
     * {@link #getPartialDerivative(int...) partial derivatives} of the value
     * returned will be
     * \(f, \frac{\partial f}{\partial q_0}, \frac{\partial f}{\partial q_1}, \ldots
     * \frac{\partial^2 f}{\partial q_0^2}, \frac{\partial^2 f}{\partial q_0 q_1},
     * \ldots \frac{\partial^o f}{\partial q_{m-1}^o\).
     * </p>
     * <p>
     * The number of parameters must match {@link #getFreeParameters()} and the
     * derivation orders of the instance and parameters must also match.
     * </p>
     * @param p base parameters with respect to which partial derivatives
     * were computed in the instance
     * @return derivative structure with partial derivatives computed
     * with respect to the lower level parameters used in the \(p_i\)
     * @since 2.2
     */
    public DerivativeStructure rebase(final DerivativeStructure... p) {

        MathUtils.checkDimension(getFreeParameters(), p.length);
        final double[][] pData = new double[p.length][];
        for (int i = 0; i < pData.length; ++i) {
            MathUtils.checkDimension(getOrder(), p[i].getOrder());
            pData[i] = p[i].data;
        }

        final DerivativeStructure result = p[0].factory.build();
        factory.getCompiler().rebase(data, p[0].factory.getCompiler(), pData, result.data);
        return result;

    }

[Serrof]

I did start on my side as well, sorry I should have warned you. I now have what I believe is a working version here: https://github.com/Serrof/hipparchus/tree/issue-191 Note that I implemented it at the DerivativeStructure level and did not deal with the Field version. In a nutshell, it introduces a dedicated class for Taylor expansions (represented by their polynomial coefficients) to go back and forth between them and sets of partial derivatives, the point being that the formula to compose the former is rather straight-forward! Let me know how we can coordinate our efforts.

[maisonobe]

I will take a more thorough look at it this week. My approach was completely different, going back to the compiled rules in the same way compileLowerIndirection and compileMultiplicationIndirection are set up in DSCompiler, with an additional cache of Rebaser engines that convert to a given number of parameters. My approach is not completed yet. I did not push it on Github to avoid triggering the continuous integration. It was a mistake as we both worked separately, I'm sorry for that.

I think I will complete what I have started, and then compare the two implementations. Then I'll keep the most efficient one as the workhorse and use the second one in the tests for cross validation, in addition to some more generic tests already set up.

[Serrof]

No worries. Your approach is probably more computationally efficient. Mine is rather light in terms of lines of code (and can still be optimized a little bit in its present state). Could be a temporary fix to the issue in order to reach the milestone for version 2.2 faster (as issue #190 depends on this as well), before something more permanent later? It's up to you.

[maisonobe]

For information, I am making good progress on this issue. A first attempt failed right after development and before validation as I noticed I made a big mistake computing the derivation formulas… In fact, they are much simpler than expected and easy to set up, once you choose the proper recursion direction (it's much simpler when recursing on derivation order than on number of parameters). So I am doing it right now and expect to finish this in the upcoming days.

[maisonobe]

I finally completed this, for both DerivativeStructure and FieldDerivativeStructure. Could you check on your side if it give the correct results and if the API suits your needs? I will take some days of rest before tackling issue #190. Once done, I guess we could publish an new Hipparchus version.

[Serrof]

Dear Luc, Thanks a lot. I'll have a look this weekend. Btw it would be awesome to have Hipparchus 2.2 in the next minor release of Orekit. Cheers, Romain.

[Serrof]

@maisonobe The API looks fine to me. I ran the same couple of tests I originally had in my fork of Hipparchus and they work. Moreover, issue #190 will validate further the implementation as it heavily relies on it.