Check elem track needs energy

[oscarxblanco]

Dear all,
this PR addresses issue #1015 , where tracking a well defined wiggler element returns nans because the energy is not defined.

As explained by @lfarv in issue #1013 , the energy could be pass in the track method, therefore, this PR checks that the energy parameter is used when any of the element has any of the two wiggler pass methods, GWigSymplecticPass or GWigSymplecticPass.

Additionally I changed numpy to np, I run black, isort and solved some flake8 messages.

[swhite2401]

From the last exchange I understand this one is ready right? Shall it be merged?

[lfarv]

@oscarxblanco, @swhite2401: As I said above, for this test to be really useful, it should also concern all the *RadPass methods and RFCavityPass. Up to you to decide to merge it as it is or to enlarge the test…

[oscarxblanco]

@lfarv and @swhite2401, RFCavityPass needs the Energy, but, BndPoleSymplectic4RadPass and StrMPoleSymplectic4RadPass don't need it.

I am adding RFCavityPass to this test.

@lfarv , is there any way to avoid hard coded pass methods inside the function ?
I think in matlab I could create a class with constants (https://es.mathworks.com/matlabcentral/answers/358826-best-method-to-define-constants-used-by-several-functions) where I list the pass methods that require Energy. And for pyat there might be already some place to put it.

[oscarxblanco]

I have tested three more pass methods and they don't need the Energy property, nor the energy parameter.

no ExactMultipoleRadPass.c
no IdTablePass.c
no QuantDiffPass.c

There is a long list of other pass methods that I would not know how to test because is not clear to me which element they are modelling.
Is there anyone willing to test them ?
The VariableThinMPolePass method needs to be fixed, somehow it crashes when using elempass, I'll do it in the branch #839 .

Apart from that, this could be merged and other methods could be added later, unless there is a suggestion on how to avoid to hard code pass methods inside elempass.m and the element_pass method in track.py.

QuadLinearFPass.c
QuadLinearPass.c
LongtAperturePass.c
AperturePass.c
BeamLoadingCavityPass.c
BeamMomentsPass.c
BendLinearPass.c
BndMPoleSymplectic4E2Pass.c
BndMPoleSymplectic4E2RadPass.c
BndMPoleSymplectic4Pass.c
BndMPoleSymplectic4QuantPass.c
BndMPoleSymplectic4RadPass.c
BndOldRadPass.c
BndStrMPoleSymplectic4Pass.c
CavityPass.c
ChangePRefPass.c
CorrectorPass.c
CrabCavityPass.c
DeltaQPass.c
DriftPass.c
EAperturePass.c
EnergyLossRadPass.c
ExactDriftPass.c
ExactMultipolePass.c
ExactMultipoleRadPass.c
ExactRectBendPass.c
ExactRectBendRadPass.c
ExactRectangularBendPass.c
ExactRectangularBendRadPass.c
ExactSectorBendPass.c
ExactSectorBendRadPass.c
IdTablePass.c
IdentityPass.c
ImpedanceTablePass.c
Matrix66Pass.c
MatrixTijkPass.c
QuantDiffPass.c
RFCavityPass.c
SimpleQuantDiffPass.c
SimpleRadiationRadPass.c
SliceMomentsPass.c
SolenoidLinearPass.c
StrMPoleSymplectic4Pass.c
StrMPoleSymplectic4QuantPass.c
StrMPoleSymplectic4RadPass.c
StrOldRadPass.c
TestRandomPass.c
ThinMPolePass.c
WakeFieldPass.c
WiggLinearPass.c
GWigSymplecticRadPass.c
GWigSymplecticPass.c
VariableThinMPolePass.c

[lfarv]

@oscarxblanco : I confirm: all the *RadPass methods need energy. For instance, in StrMPoleSymplectic4RadPass:

    gamma = atGamma(Param->energy, Elem->Energy, Param->rest_energy);

That makes sense: the radiated power depends highly on the particle energy ($\gamma^3$ or $\gamma^5$). So we have

    double rad_const = RAD_CONST*pow(gamma, 3);
    double diff_const = DIF_CONST*pow(gamma, 5);

So a rule could be:

• as soon as radiation is involved (*RadPass, *QuantPass), you need energy,
• RF cavities need energy (RFCavityPass, CrabCavityPass, BeamLoadingCavityPass,…)
• otherwise no, with an exception for GWigSymplecticPass, because of the field in T, which must be normalised with the energy. All the other methods use a normalised field in $m^{-n}$.

I don't see how to avoid a hard-coded definition of the passmethods needing energy. Most elements potentially radiating should derive from the _Radiative class, but only the passmethod (a string) tells whether they are effectively radiating or not. Up to now there is no exhaustive list of radiating passmethods (anybody can add his own private ones).

[oscarxblanco]

Thank you @lfarv for the explanation, you are right. I think I was oversimplifying the tests I did yesterday.

@lfarv , another alternative is to check for well defined energy inside the pass method, using your PR #1021 . This will avoid the hard coded list and make the pass methods a little more save to use at the expense of an extra check. What do you think ?

[lfarv]

@oscarxblanco : I think that the test inside the pass method is safer. But it's more modifications (~15 passmethods), and it's tricky: to avoid a memory leak, all the tests must be done before the element structure is allocated. At first glance, the call to atEnergy/atGamma must be duplicated and moved into both branches of the if(!Elem) test.

[oscarxblanco]

Dear @lfarv ,
I removed the python and matlab modifications, and for the wiggler pass method in C I tried to check Param->energy before memory is allocated in ExportMode struct elem *trackFunction(...) inside GWigSymplecticPass.c.

    if (Param->energy == 0){
            atError("Energy needs to be defined. Check lattice parameters or pass method options.\n");
            check_error();
    }

It works using your last modification in #1021 , but, it is not the duplicate call to atGamma/atEnergy you suggested because, different to matlab, in pyat there is an inline definition without checks.

at/atintegrators/atelem.c

Lines 209 to 210 in a3a8783

	#define atEnergy(ringenergy,elemenergy) (ringenergy)
	#define atGamma(ringenergy,elemenergy,rest_energy) ((rest_energy) == 0.0 ? 1.0E-9*(ringenergy)/__E0 : (ringenergy)/(rest_energy))

As for AT matlab, the energy check is already there, atGamma returns an error with a different message because I have not modified it.

at/atintegrators/atelem.c

Lines 77 to 97 in a3a8783

	double atEnergy(double ringenergy, double elemenergy)
	{
	if (ringenergy!=0.0)
	return ringenergy;
	else
	if (elemenergy!=0.0)
	return elemenergy;
	else {
	atError("Energy not defined.");
	return 0.0; /* Never reached but makes the compiler happy */
	}
	}
	double atGamma(double ringenergy, double elemenergy, double rest_energy)
	{
	double energy = atEnergy(ringenergy, elemenergy);
	if (rest_energy == 0.0)
	return 1.0E-9 * energy / __E0;
	else
	return energy / rest_energy;
	}

What do you think ? I leave here below the test code I wrote for pyat and AT matlab.

import at
import numpy as np

zin = np.zeros((6,1))
zin[1] = 1e-5

und = at.Wiggler('und', length=1, wiggle_period=1/100, b_max=0,Nmeth=2,Nstep=10)
ring = at.Lattice([und],energy=1e9)
print(ring.track(zin)[0].T)
print(und.track(zin,energy=1e9).T)
print(und.track(zin)) # this should fail

clear
%%
% cd ~/atfortest/at
cd atmat
atpath
atmexall
% cd ~/test_passmethod_needs_energy
%%
und = atwiggler('und', 1, 1/100, 0, 1e9);
ringpara = {atringparam('UND',energy=1e9)};
ring = [ringpara {und}];
zin = zeros(6,1);
zin(2) = 1e-5;
%%
ringpass(ring,zin)'
linepass(ring,zin)'
elempass(und,zin,energy=1e9)'
elempass(und,zin)
und.Energy = [];
elempass(und,zin) % this should fail

[lfarv]

A comment before you replicate this in all the passmethods.

In Matlab, the energy may be given by an element attribute. So you have to check the result of atEnergy/atGamma, and it has to be checked before the memory allocation:

    if (!Elem) {
        double *R1, *R2, *T1, *T2;
        double *By, *Bx;
        double Ltot, Lw, Bmax, Energy;
        int Nstep, Nmeth;
        int NHharm, NVharm;

        Ltot = atGetDouble(ElemData, "Length"); check_error();
        ...
        /* Optional fields */
        Energy=atGetOptionalDouble(ElemData,"Energy",Param->energy); check_error();
        ...

        gamma = atGamma(Param->energy, Energy, Param->rest_energy);
        if (gamma == 0) {
            atError("Energy needs to be defined. Check lattice parameters or pass method options.\n");
            check_error();
        }

        Elem = (struct elem*)atMalloc(sizeof(struct elem));
        ...
    }
    else {
        gamma = atGamma(Param->energy, Elem->Energy, Param->rest_energy);
    }

Note that the 2nd gamma computation is done because Param->energy may be used for energy ramping, and may change on each call. But there I don't think it needs to be check again.

[oscarxblanco]

@lfarv, right, I didn't thought about ramping. And I think I miss the case in matlab when the wiggler defines the energy. Thank you.

I just uploaded a new version, should I modify other pass methods as we did with this one ?

[oscarxblanco]

@lfarv , could you check the latest modification of RFCavityPass.c ? Is it OK for you ?

[lfarv]

@oscarxblanco : this looks correct.

Another option is to move the test inside the python atEnergy (similarly to the Matlab version) and follow the call with a check_error. This would make Matlab and python versions more similar and avoid a double test in Matlab. But atEnergy/atGamma still have to be moved to where you just put them. Up to you !

[lfarv]

@oscarxblanco : after checking RFCavityPass, I notice that unlike all other elements, RFCavityPass, as it is in the master branch, does use the element Energy field if necessary. Do we keep a special case for it?

Yes, if necessary. I would prefer to have all methods with similar checks of the energy but I will have to do try it out.

After looking at that, I think that the special case of RFCavity is related to the fact that the energy attribute is mandatory for RF cavities: it's one of the arguments of the RF constructor of the cavity element, which was made similar the the Matlab one.
However, I think that we should keep the same rule for all elements, and so ignore this element energy attribute for the RF cavity also.

[oscarxblanco]

Dear @lfarv, I think this PR is ready to be reviewed and merged if you agree.

The energy parameter is checked in every trackFunction implementation in C requiring it. This has an effect on pyAT for lattice and element tracking, and only in lattice tracking in AT matlab.

The behaviour of AT matlab elempass, for element tracking, has not changed because the element tracking is done in a separate mexFunction implementation.

The mexFunctions could also be modified to match the behaviour of pyAT at the expense of backwards compatibility, for example, tracking a single wiggler will always require the energy as a parameter in the tracking call and the wiggler energy field would be always ignored. This is what I was actually aiming at when starting this PR but at this point I think it would be longer to review. Should I implement it now ?

[lfarv]

Ok !