atcollab / at

Accelerator Toolbox
Apache License 2.0
48 stars 31 forks source link

ringparam always returns an emittance #583

Closed swhite2401 closed 1 year ago

swhite2401 commented 1 year ago

After PR #511 many users scripts broke because NaNs were returned for lattices considered as off-energy. This prevented using this function on lattice with errors following an orbit correction, this is a source of complaints and major disturbance to users. This PR proposes instead to always return the on-energy emittances and to modify the warning explicitly mentioned that the emittances returned are for the on-energy lattice. This restores backward compatibility while warning the user that he might be doing something wrong.

lfarv commented 1 year ago

April's fool ? Are you seriously considering deliberately returning wrong results ? I fear what the consequences of generalising this idea could be…

Seriously, what would you expect from such results, even for "compatibility" reasons ? @simoneliuzzo showed the error may be very large!

And back to compatibility, before #511 there was NO 'dp' input argument in ringpara or atsummary. So nobody can have used ringpara or atsummary off-momentum.

If "some user" has a real and urgent need for such a trick, you may consider making a fork, and urge him to carefully keep his results confidential…

swhite2401 commented 1 year ago

This function accepts off-momentum lattices: 6D with non nominal RF. So Let's keep this discussion correct please, this is not a joke or a trick.

Now for the function there are obvious problems and they need to be fixed, this proposal is one solution, that is for me perfectly acceptable, but there may be others. I just try to help people that suffer from these 'bugfixes', so please if you have better ideas propose them and let's try to improve this with an open mind. If we don't you are in fact perfectly right: new forks risk to appear and users disappear...

One question: If I understand well the code, you decided, to set the max. deviation from nominal RF for off-momentum at 2Hz ? Why is that? And why not 10Hz? Did you run simulation to define this number? the impact on the final result is machine dependent so does it even make sense to a non zero value?

simoneliuzzo commented 1 year ago

I have nothing against a warning, but changing the output is a clearly non-backward compatible action.

In my specific case the pull request #511 killed my code with an error.

What I need to do is to compute the bunch length for non-zero current on a lattice with errors and frequency not to the nominal setting. Could you please advise on how to proceed with these computations in the correct way?

In my opinion the problem is that the computation of the natural frequency of the SR is incorrect when using a lattice with errors. I use to set the correct frequency within my functions, but now this action returns an error concerning off-momentum that is clearly off scope in my case. A lattice with errors has a different nominal frequency than one without errors, and is not off-momentum.

thank you for your help best regards Simone

simoneliuzzo commented 1 year ago

By the way, PR #511 has been merged without the approval of any of the 3 requested reviewers.

swhite2401 commented 1 year ago

Alternatively we may call ohmi_envelop for 6D lattices, this would give the correct results also with errors, no approximation

simoneliuzzo commented 1 year ago

@swhite2401 , Yes, I think I will do like you say, and modify my code to use ATX instead of ringpara.

swhite2401 commented 1 year ago

Well what I was suggesting was to integrate it in ringparam, I give it a try and let you know

simoneliuzzo commented 1 year ago

I would also like to add that non-backward compatible bug fixes such as the one in PR #511 would deserve in my opinion a TAG before the change is merged. This would allow user as me to freeze their code to this given tag.

carmignani commented 1 year ago

Hello, I don't think this is an acceptable modification. I think here the problem is always the same: ringpara was a 4D function, as atlinopt, but it was accepting 6D lattices giving slightly wrong results. I was using in that way a lot, so sometimes I have problems running old scripts. But the fact that now bugs are not allowed can't be considered as breaking backward compatibility. Now we can give a 6D lattice and have good 6D results, except the emittance because the quadrupole component of the radiation is missing. In my opinion, until we don't implement properly the off-energy emittance computation including the quadrupoles in the radiation integrals we can't give an emittance result.

lfarv commented 1 year ago

What I need to do is to compute the bunch length for non-zero current on a lattice with errors and frequency not to the nominal setting. Could you please advise on how to proceed with these computations in the correct way?

ohmienvelope /atx: it's the only way. Even if improving the computation of radiation integrals, you'll never take coupling into account. And errors may introduce such coupling.

In my opinion the problem is that the computation of the natural frequency of the SR is incorrect when using a lattice with errors. I use to set the correct frequency within my functions, but now this action returns an error concerning off-momentum that is clearly off scope in my case. A lattice with errors has a different nominal frequency than one without errors, and is not off-momentum.

The nominal frequency is the frequency corresponding to the nominal path length. I does not depend on errors. If you have another definition, can you explicit it ?

I would also like to add that non-backward compatible bug fixes such as the one in PR https://github.com/atcollab/at/pull/511 would deserve in my opinion a TAG before the change is merged. This would allow user as me to freeze their code to this given tag.

You can always go back to any point and freeze you version by using the commit hash: a tag is just a name associated with a commit.

Well what I was suggesting was to integrate it in ringparam, I give it a try and let you know

Why not, but be careful, this will not be backward compatible, unless you introduce new fields for ohmienvelope results… But then, why not using directly ohmienvelope if what you need is bunch length ?

simoneliuzzo commented 1 year ago

@lfarv ,

I have already started my migration away from matlab AT about 1 year ago. I plan to abandon it completely by the end of 2023. The updates and bug fixes made in matlab AT imply modifications that are comparable to rewriting my code. So I started rewriting it in pyAT. If also pyAT will undergo as often as matlab AT non-backward compatible changes, I will soon stop using it, and start again from scratch using the python xsuite code, presently developed at CERN.

I need bunch length with current, I use the function BL = atBunchLength(RING, I_bunch, Zn) . https://github.com/atcollab/at/blob/master/atmat/atphysics/LongitudinalDynamics/atBunchLength.m

It is this function that calls ringpara, and should call atx instead.

The tag is a way to inform users that from there on, there is something really different. Is not a simple commit. It is a peculiar commit, that needs to be remembered.

best regards "some-user"

lfarv commented 1 year ago

But the fact that now bugs are not allowed can't be considered as breaking backward compatibility.

I fully agree with @carmignani: returning wrong results is not an option, and the only way to improve the present state is to add the quadrupole contribution to the integrals. Whether this is worth the effort, knowing that ohmienvelope is available, and that anyway the integrals won't cover the H/V coupling is another question. But if someone does it it's perfect.

carmignani commented 1 year ago

I need bunch length with current, I use the function BL = atBunchLength(RING, I_bunch, Zn) . https://github.com/atcollab/at/blob/master/atmat/atphysics/LongitudinalDynamics/atBunchLength.m

Why do you want to use this function with a 6D lattice? This was allowed in the past, because ringpara allowed you to pass a 6D lattice and it was a bug. Now if you pass a 6D lattice ringpara computes the parameters correctly, except the emittance and the bunch length because the quadrupoles are not included in the integrals.

For this specific problem, I would add an atradoff inside atBunchLength. So the lattice is converted to 4D and there are no problems. I don't think this would create any problem to existing codes.

simoneliuzzo commented 1 year ago

The nominal frequency for a lattice with errors could be the frequency related to the electron path length, not to the lattice circumference.

For example using this function (from @lfarv, June 2016)

indd=find(atgetcells(ring,'Class','Bend'));

b=atgetfieldvalues(ring,indd,'BendingAngle');

o=findorbit6(ring,indd);
o1=findorbit6(ring,indd+1);

xave=(o1(1,:)+o(1,:))/2; % average closed orbit at dipole

dl=sum(xave.*b');

Simulations of 50 possible lattice surveys show consistently +2kHz change of frequency. The EBS operation frequency is in fact 2kHz larger (352174152) than the one computed for a perfect ring (352172152).

simoneliuzzo commented 1 year ago

@carmignani I use BL = atBunchLength(RING, I_bunch, Zn) because I need bunch lengthening with current. There is nothing that stops me to do this in BL = atBunchLength(RING, I_bunch, Zn). The help does not mention the lattice has to be 4D. The function does not explicitly reject 6D lattices. For me this means the function works for any input ring. At least this is what "some-user" expects.

I have nothing against forcing atradoff or replacing ringpara with atx.

lfarv commented 1 year ago

@simoneliuzzo : going to PyAT is indeed a good option, but I don't think it makes a difference concerning bug fixes: every software has bugs, and each bug fix is by definition not backward compatible, you don't get the same result as before. And don't forget that python is evolving fast. Each python release has its set of deprecated features. For example, look at the python 3.11 release notes:

There is one python version per year, here is the official schedule. AT is following the python release schedule. As much as possible, we try to keep as many old versions as possible, but by now, we already dropped support for python 2.7, 2.5 and 3.6. Python 3.7 is almost 5 years old, its end-of-life is scheduled on 27 June 2023. Not talking of numpy, which also has its own schedule

So you are not finished with keeping your code up-to-date !

carmignani commented 1 year ago

I think I wrote atBunchLength and at the time I was ignoring the fact that ringpara was giving wrong results in case of 6D lattice. We could put atx instead of ringpara, but then the function would be slower. If we put atradoff inside atBunchLength we wouldn't have the good bunch length for off-energy lattices, but that would be the same with Simon's proposal. Probably the best would be to keep ringpara in case of 4D lattice and use atx in case of 6D lattice. I could do this modification, if you all agree.

simoneliuzzo commented 1 year ago

@carmignani ok for me

swhite2401 commented 1 year ago

@carmignani ok thanks, you can use this branch, so you have the discussion history. Could you please also modify the warning to point the user to atx?

carmignani commented 1 year ago

I can use this branch, but I think we should revert the commits about ringpara.

swhite2401 commented 1 year ago

@carmignani in fact, no please leave this branch to me, I have an idea from a simple integration of the quadrupole contribution, I'll put it here

lfarv commented 1 year ago

@carmignani : very good solution. Then ringpara can go back to its initial state, until someone looks for quadrupole integration.

carmignani commented 1 year ago

Ok, then if the quadrupole contribution will work, there will be no need to change atBunchLength. So I wait for you modification and then we can see if I have to update atBunchLength or not.

swhite2401 commented 1 year ago

Please modify atBunchLength in any case, it will take a bit of time to implement and test and I don't think speed is really something critical for atBunchLength

btw if you need bunch length only instead of atx you may use ohmienvelope directly to save some calculations.

carmignani commented 1 year ago

Ok, then I do another branch just for atBunchLength. I'll keep ringpara in case of 4D input and atx or ohmi_envelope in case of 6D. I have to check exactly what is needed.

lfarv commented 1 year ago

Please modify atBunchLength in any case, it will take a bit of time to implement and test and I don't think speed is really something critical for atBunchLength

btw if you need bunch length only instead of atx you may use ohmienvelope directly to save some calculations.

I agree, this way atBunchLength will be correct in any case.

carmignani commented 1 year ago

@simoneliuzzo

The help does not mention the lattice has to be 4D. The function does not explicitly reject 6D lattices.

I just checked and actually the help says: ring is the ring without radiation

swhite2401 commented 1 year ago

I have done a quick implementation in pyAT using the constant rho approximation, it is derived from the closed orbit. Results on the figure below, the error on the emittance is ~2% at for an RF shift of 300Hz, I will double check if I got everything right but I think this as good as it gets with this method. What do you think?

image

swhite2401 commented 1 year ago

Some complementary information: this is without radiation in quadrupoles, the error comes from the dipole.

image

lfarv commented 1 year ago

@swhite2401 : This looks excellent. Do you intend to implement it in both Matlab and python in this branch?

swhite2401 commented 1 year ago

I have done in it python I can propose something for matlab as well but I will need your help to clean it up.

Any idea why we have such error from the dipole? Can it be improved in some way? I haven't looked at the details of the implementation but could it be due to the DQ contribution for non ideal particles?

lfarv commented 1 year ago

In python there is already a get_radiation_integrals function, I hope you can just improve it. At the moment, it accepts only dp. It could be the opportunity to also allow dct and df.

In Matlab, there are 2 separated functions DipoleRadiation and WigglerRadiation. Either we add a third one, or we group everything inside a new global function built on top.

Now for the difference between both methods, no idea. It looks like you did your tests with an hmba lattice.May be looking at a simpler one (only 1 dipole type) could give hints… But I don't think the difference justifies too much work right now.

swhite2401 commented 1 year ago

Matlab and python implemented, this is the method that produced the figures above

swhite2401 commented 1 year ago

The matlab test are still failing... I have to figure this out

lfarv commented 1 year ago

Does anybody knows where these smileys come from ? nice, but ok, I got it !

swhite2401 commented 1 year ago

otherwise we should check for (abs(theta) < epsilon)rather than (theta == 0).

Yes in fact I wasn't so sure about this... I left it like this because 0.0 is forbidden as it would return NaN for I5. Very small values on the other hand should still give the correct result. But I will check that this does not diverge you are right.

swhite2401 commented 1 year ago

@lfarv, I have run scan down to 0.001Hz frequency shifts and everything looks stable to me, however just to be on the safe side I will put a minimum angle at 1.0e-7.

I have a question, now that things are correctly handle should we restore the dp input for matlab?

lfarv commented 1 year ago

should we restore the dp input for matlab?

What do you mean ? dp input where ? Anyway, everything looks correct, so I would say yes !

swhite2401 commented 1 year ago

What do you mean ? dp input where ? Anyway, everything looks correct, so I would say yes !

Sorry forget this, I thought you had removed the dp input to atsummary and ringparam since offmomentum was wrong, but in fact not.

Ok so I merge.