Clarification/bug .plug() usage with inverted maps

[WK-astrodynamics]

Hi,

I was using the plug command for a DACE::AlgebricVector either incorrectly or there is a bug. In the first case, I believe documentation could be improved for clarity. For the later case it could be a feature improvment.

The plug command as of V2.0.1 has this documentation:
Partial evaluation to replace given independent DA variable by value val. DA only.

Say I have a function evaluation with DA types, this is a map taking DA variable (inputs) to their evaluations (δx, δy, δz) --> (xf, yf, yz)
If I create the inversion of the map via .invert() my inverted map will be (xf, yf, yz) --> (δx, δy, δz).
Therefore, the variables on the left hand side are my new idependent variables for this inverted map object.

However, .plug does not act on these independent variables but on the origional ordering

It's best visualised with the following code and comments

#define WITH_ALGEBRAICMATRIX
#include "dace/dace.h"
#include <cmath>
#include <iostream>

#include "plug_inversion.h"

template<typename T>
DACE::AlgebraicVector<T> func(T x, T y, T z) {
    // simple math function
    // not realistics but demonstraits the idea
    return DACE::AlgebraicVector<T>{
            2.0*x + y,
            y + z,
            z + x,
    };
}

int main2() {
    const int ORDER  = 2;
    const int N_VARS = 3;
    DACE::DA::init(ORDER, N_VARS);
    DACE::DA::setEps(1e-300);

    // scale factor for DA independent variables
    const double SCALE = 1e-3;

    // setting inital expansion points
    DACE::DA x = 0.1 + SCALE * DACE::DA(1);
    DACE::DA y = 0.2 + SCALE * DACE::DA(2);
    DACE::DA z = 0.0 + SCALE * DACE::DA(3);

    // create a mapping of the function, func: R^3 -- > R^3; (x,y,z) |--> (xf, yf, zf)
    auto func_map = func(x, y, z);
    std::cout << "function evaluation map:\n" << func_map << "\n";

    // rearange the map, R^3 -- > R^3; (x,y,z) |--> (yf, x, z)
    DACE::AlgebraicVector<DACE::DA> map = {
            func_map[1],
            DACE::DA(1), // identity map below
            DACE::DA(3)
    };

    // create a inversion of the map, func^-1: R^3 -- > R^3; (yf, x, z) |--> (x,y,z)
    auto inv_map = map.invert();
    std::cout << "map inversion:\n" << inv_map << "\n";


    // say we wish to evalute the map when yf = 0.
    // I would assume that plug would do the evaluation for the first independent variable of the inverted map
    // not global DA
    auto constraint = inv_map.plug(1, 0.0); // should be .eval
    std::cout << "constraint:\n" << constraint << "\n";
    
    auto solution = map.eval(constraint);
    std::cout << "solution:\n" << solution << "\n";
    // but this is not true, notice that output values are not (x, y) |-- > (xf, yf, zf)
    // instead we have (y, z) |-- > (xf, yf, zf)
    
    return 1;
}

The problem is that the output of solution should be removing the variable yf or constrain yf in terms of identites shown in map. Since yf is the first independent variable of the inversed map. Thus, the solution being a composition, should have variables δx and δz but instead it has δy and δz.

Suggestion:

1. Potentially allow plug to evaluate in this manner making partial inversions of a DA objects nice instead of using .eval{0.0, DACE::DA(1), DACE::DA(3)},
2. Or, state in the documentation that the plug command is supposed work with the global original independent varaibales and not independent variables of the current DA object, unless they are the same.

Sorry if this made no sense, I can try to clarify what I mean.