minv not accounting for parameters?

[mattsignorelli]

Say I have a GTPSA with 1 variable x and 1 parameter k, all to 2nd order, and I make a map:

m = (1+k)*x

I would expect mad_tpsa_minv(1, m, m_inv) to give m_inv = x/(1+k) = x - kx , however it is just returning x. Unless I am missing something?

[ldeniau]

Here is the code (in MAD-NG) showing the result you expected, but this is not a damap inversion. Parameters in damap are immutable 1st-order TPSA shared by all damaps of the same structure. This is why below I create a knob 'k' (external to the damap) from the parameters 'k', and a TPSA 'a' to show the inverse. Concerning mad_tpsa_inv, I will come back to you with another example, but inverting a damap should not clear the derivatives in the variables, which is the case in your example above for 'x' and k=1.

> m=MAD.damap{nv=1,mo=2,np=1,po=2,vn={'x','k'}}
> m.x=3
> k=2+m.k
> a=(1+k)*m.x
> a:print()

 -UNNAMED-:  R, NV =   1, MO =  2, NP =   1, PO =  2
 *******************************************************
     I   COEFFICIENT             ORDER   EXPONENTS
     1   9.0000000000000000E+00    0     0
     2   3.0000000000000000E+00    1     1
     3   3.0000000000000000E+00    1     0  2^1
     4   1.0000000000000000E+00    2     1  2^1
> a:inv():print()

 -UNNAMED-:  R, NV =   1, MO =  2, NP =   1, PO =  2
 *******************************************************
     I   COEFFICIENT             ORDER   EXPONENTS
     1   1.1111111111111110E-01    0     0
     2  -3.7037037037037035E-02    1     1
     3  -3.7037037037037035E-02    1     0  2^1
     4   1.2345679012345678E-02    2     2
     5   1.2345679012345678E-02    2     1  2^1
     6   1.2345679012345678E-02    2     0  2^2

[mattsignorelli]

I attached two programs. test_var.c inverts the map:

M(\vec{x}) =\begin{pmatrix} (1+x_2)x_1 \\\ x_2 \end{pmatrix}

treating $x_2$ as a variable.

test_param.c inverts the map:

M(\vec{x}) =\begin{pmatrix} (1+x_2)x_1 \end{pmatrix}

treating $x_2$ as a parameter . The results are different. Archive 2.zip

[mattsignorelli]

Here is the output for test_var.c:

Map1:
 (1+x2)*x1:  R, NV =   2, MO =  2
 *******************************************************
     I   COEFFICIENT             ORDER   EXPONENTS
     1   1.0000000000000000E+00    1     1 0
     2   0.0000000000000000E+00    1     0 1
     3   0.0000000000000000E+00    2     2 0
     4   1.0000000000000000E+00    2     1 1
     5   0.0000000000000000E+00    2     0 2

 x2      :  R, NV =   2, MO =  2
 *******************************************************
     I   COEFFICIENT             ORDER   EXPONENTS
     1   0.0000000000000000E+00    1     1 0
     2   1.0000000000000000E+00    1     0 1

Map1^(-1)
 x1 - x2*x1:  R, NV =   2, MO =  2
 *******************************************************
     I   COEFFICIENT             ORDER   EXPONENTS
     1   1.0000000000000000E+00    1     1 0
     2   0.0000000000000000E+00    1     0 1
     3   0.0000000000000000E+00    2     2 0
     4  -1.0000000000000000E+00    2     1 1
     5   0.0000000000000000E+00    2     0 2

 x2      :  R, NV =   2, MO =  2
 *******************************************************
     I   COEFFICIENT             ORDER   EXPONENTS
     1   0.0000000000000000E+00    1     1 0
     2   1.0000000000000000E+00    1     0 1

Map1*Map1^(-1)
 x1      :  R, NV =   2, MO =  2
 *******************************************************
     I   COEFFICIENT             ORDER   EXPONENTS
     1   1.0000000000000000E+00    1     1 0
     2   0.0000000000000000E+00    1     0 1

 x2      :  R, NV =   2, MO =  2
 *******************************************************
     I   COEFFICIENT             ORDER   EXPONENTS
     1   0.0000000000000000E+00    1     1 0
     2   1.0000000000000000E+00    1     0 1

And test_param.c:

Map1:
 (1+x2)*x1:  R, NV =   1, MO =  2, NP =   1, PO =  2
 *******************************************************
     I   COEFFICIENT             ORDER   EXPONENTS
     1   1.0000000000000000E+00    1     1
     2   0.0000000000000000E+00    1     0  2^1
     3   0.0000000000000000E+00    2     2
     4   1.0000000000000000E+00    2     1  2^1
     5   0.0000000000000000E+00    2     0  2^2

 x2      :  R, NV =   1, MO =  2, NP =   1, PO =  2
 *******************************************************
     I   COEFFICIENT             ORDER   EXPONENTS
     1   0.0000000000000000E+00    1     1
     2   1.0000000000000000E+00    1     0  2^1

Map1^(-1)
 SHOULD BE x1 - x2*x1 ?? :  R, NV =   1, MO =  2, NP =   1, PO =  2
 *******************************************************
     I   COEFFICIENT             ORDER   EXPONENTS
     1   1.0000000000000000E+00    1     1
     2   0.0000000000000000E+00    1     0  2^1

 x2      :  R, NV =   1, MO =  2, NP =   1, PO =  2
 *******************************************************
     I   COEFFICIENT             ORDER   EXPONENTS
     1   0.0000000000000000E+00    1     1
     2   1.0000000000000000E+00    1     0  2^1

Map1*Map1^(-1)
 SHOULD BE x1 ??:  R, NV =   1, MO =  2, NP =   1, PO =  2
 *******************************************************
     I   COEFFICIENT             ORDER   EXPONENTS
     1   1.0000000000000000E+00    1     1
     2   0.0000000000000000E+00    1     0  2^1
     3   0.0000000000000000E+00    2     2
     4   1.0000000000000000E+00    2     1  2^1
     5   0.0000000000000000E+00    2     0  2^2

 x2      :  R, NV =   1, MO =  2, NP =   1, PO =  2
 *******************************************************
     I   COEFFICIENT             ORDER   EXPONENTS
     1   0.0000000000000000E+00    1     1
     2   1.0000000000000000E+00    1     0  2^1

[ldeniau]

I updated the library to give a better meaning of parameters (added mad_tpsa_setprm and enforced proper use of mad_tpsa_setvar), as I think the confusion is coming from there. I have attached four files where you can see at the top the definition of the map and the way I initialized them and at the bottom the output I obtained, which seems consistent to me unless I missed something. Original versions use an intermediate variable x2p (x2 plus something) like yours, while versions 2 use directly m. I also introduce values with different inverses to avoid 1^(-1) = 1 which could confuse the understanding.

test_var_vs_prm.zip

[ldeniau]

Your confusion may come from how a map inverse is computed, the steps are basically inverse first order like a matrix then split the linear part and the nonlinear part, and compose iteratively order by order the inverses, The algorithm is similar to the normal form but works on maps instead of Lie maps. Now when you have parameters, the composition must exclude them, i.e. variables are not injected in parameters and parameters are immutables, hence 1 var + 1 prm means composing an array of one TPSA only.

If you want to have a look at the code, quite simple, it's here (2 functions and ~100 lines): https://github.com/MethodicalAcceleratorDesign/MAD-NG/blob/32a63c9a98948905ac5ee9c56022d3bc61695913/src/mad_tpsa_minv.c#L119

[mattsignorelli]

Ok, so inverting a map does not include inversion of the parameters, only the variables? What exactly are parameters used for then, and are there any speedups for using parameters instead of variables for everything in certain scenarios (is there any optimization over knowing that dk/dx =0 where k is a knob and x is a map variable)?

[DavidSagan]

Your confusion may come from how a map inverse is computed,

It is not a question or confusion over how the inverse is computed, it is a question of usability. And for normal form, Matt and I and Etienne need an inverse that behaves in the way that Matt explained in the original post.

[ldeniau]

Your confusion may come from how a map inverse is computed,

It is not a question of confusion over how the inverse is computed, it is a question of usability. For normal form, Matt and I and Etienne need an inverse that behaves like Matt explained in the original post.

To my knowledge, no other library supports parameters as immutable shared TPSAs of order 1 in DA maps, and they were introduced to optimize speed and memory while making knobs (i.e. deferred expressions of TPSAs) easy for users to use, and fulfill these goals very well. So far, my parametric normal forms work well for many studies, and in fact, inverses with parameters work as expected, mainly for performance reasons. We can very easily calculate derivatives of Tune, Chroma, RDTs, etc. vs user knobs, including extracting the Jacobian vs user knobs for optimizers.

For now, the intended behavior can be obtained by replacing parameters with variables as in other codes (e.g. PTC/FPP, Matt examples). But Berz's package does not allow for handling many variables, nor variables with different orders, so the number of "parameters" as variables remains very limited.

That being said, I am starting to work on extracting parametric Hamiltonian from DA maps, and I agree that it would be much simpler to get the intended behavior you are asking for, as for now I am converting parameters to variables to obtain the proper behavior in the high-level Lua code.

Thus, I will soon improve the library to support this behavior and probably take this opportunity to improve some other aspects (from my TODO list). This will change the signature of the few functions accepting an array of TPSA as arguments.

[ldeniau]

The intended behavior has been implemented in 7cf8066d The function minv and pminv need an extra argument See attached examples in C and in MAD-NG.

test_vp.zip

./test_var

*** nv=2, mo=2, np=0, po=0

Map1:
 x1      :  R, NV =   2, MO =  2
 *******************************************************
     I   COEFFICIENT             ORDER   EXPONENTS
     1   2.0000000000000000E+00    0     0 0
     2   1.0000000000000000E+00    1     1 0
     3   0.0000000000000000E+00    1     0 1

 x2      :  R, NV =   2, MO =  2
 *******************************************************
     I   COEFFICIENT             ORDER   EXPONENTS
     1   0.0000000000000000E+00    1     1 0
     2   1.0000000000000000E+00    1     0 1

Map1: x1 <- (3+x2)*x1
 x1      :  R, NV =   2, MO =  2
 *******************************************************
     I   COEFFICIENT             ORDER   EXPONENTS
     1   6.0000000000000000E+00    0     0 0
     2   3.0000000000000000E+00    1     1 0
     3   2.0000000000000000E+00    1     0 1
     4   0.0000000000000000E+00    2     2 0
     5   1.0000000000000000E+00    2     1 1
     6   0.0000000000000000E+00    2     0 2

 x2      :  R, NV =   2, MO =  2
 *******************************************************
     I   COEFFICIENT             ORDER   EXPONENTS
     1   0.0000000000000000E+00    1     1 0
     2   1.0000000000000000E+00    1     0 1

Map1^(-1):
 x1      :  R, NV =   2, MO =  2
 *******************************************************
     I   COEFFICIENT             ORDER   EXPONENTS
     1   3.3333333333333331E-01    1     1 0
     2  -6.6666666666666663E-01    1     0 1
     3   0.0000000000000000E+00    2     2 0
     4  -1.1111111111111110E-01    2     1 1
     5   2.2222222222222221E-01    2     0 2

 x2      :  R, NV =   2, MO =  2
 *******************************************************
     I   COEFFICIENT             ORDER   EXPONENTS
     1   0.0000000000000000E+00    1     1 0
     2   1.0000000000000000E+00    1     0 1

Map1*Map1^(-1):
 x1      :  R, NV =   2, MO =  2
 *******************************************************
     I   COEFFICIENT             ORDER   EXPONENTS
     1   6.0000000000000000E+00    0     0 0
     2   1.0000000000000000E+00    1     1 0
     3   0.0000000000000000E+00    1     0 1

 x2      :  R, NV =   2, MO =  2
 *******************************************************
     I   COEFFICIENT             ORDER   EXPONENTS
     1   0.0000000000000000E+00    1     1 0
     2   1.0000000000000000E+00    1     0 1

./test_param

*** nv=1, mo=2, np=1, po=2

Map1:
 x1      :  R, NV =   1, MO =  2, NP =   1, PO =  2
 *******************************************************
     I   COEFFICIENT             ORDER   EXPONENTS
     1   2.0000000000000000E+00    0     0
     2   1.0000000000000000E+00    1     1
     3   0.0000000000000000E+00    1     0  2^1

 x2      :  R, NV =   1, MO =  2, NP =   1, PO =  2
 *******************************************************
     I   COEFFICIENT             ORDER   EXPONENTS
     1   0.0000000000000000E+00    1     1
     2   1.0000000000000000E+00    1     0  2^1

Map1: x1 <- (3+x2)*x1
 x1      :  R, NV =   1, MO =  2, NP =   1, PO =  2
 *******************************************************
     I   COEFFICIENT             ORDER   EXPONENTS
     1   6.0000000000000000E+00    0     0
     2   3.0000000000000000E+00    1     1
     3   2.0000000000000000E+00    1     0  2^1
     4   0.0000000000000000E+00    2     2
     5   1.0000000000000000E+00    2     1  2^1
     6   0.0000000000000000E+00    2     0  2^2

 x2      :  R, NV =   1, MO =  2, NP =   1, PO =  2
 *******************************************************
     I   COEFFICIENT             ORDER   EXPONENTS
     1   0.0000000000000000E+00    1     1
     2   1.0000000000000000E+00    1     0  2^1

Map1^(-1):
 x1      :  R, NV =   1, MO =  2, NP =   1, PO =  2
 *******************************************************
     I   COEFFICIENT             ORDER   EXPONENTS
     1   3.3333333333333331E-01    1     1
     2  -6.6666666666666663E-01    1     0  2^1
     3   0.0000000000000000E+00    2     2
     4  -1.1111111111111110E-01    2     1  2^1
     5   2.2222222222222221E-01    2     0  2^2

 x2      :  R, NV =   1, MO =  2, NP =   1, PO =  2
 *******************************************************
     I   COEFFICIENT             ORDER   EXPONENTS
     1   0.0000000000000000E+00    1     1
     2   1.0000000000000000E+00    1     0  2^1

Map1*Map1^(-1):
 x1      :  R, NV =   1, MO =  2, NP =   1, PO =  2
 *******************************************************
     I   COEFFICIENT             ORDER   EXPONENTS
     1   6.0000000000000000E+00    0     0
     2   1.0000000000000000E+00    1     1
     3   0.0000000000000000E+00    1     0  2^1

 x2      :  R, NV =   1, MO =  2, NP =   1, PO =  2
 *******************************************************
     I   COEFFICIENT             ORDER   EXPONENTS
     1   0.0000000000000000E+00    1     1
     2   1.0000000000000000E+00    1     0  2^1