add check_6d in findorbit4 and atlinopt4

[oscarxblanco]

This pull request adds a warning when atlinopt4d is used on a 6D lattice.

[oscarxblanco]

Dear @swhite2401 , this is a proposal to solve issue #788 . Who would review it ?

[oscarxblanco]

I also added a check on findorbit4.

[lfarv]

Hello @oscarxblanco This is a very sensitive topic. The check was removed because some users have been wrongly using atlinopt/atlinop4 with 6d lattices for long, and in very large software. They complained that it broke existing software, even if they got wrong results… "backward compatibility" question!

In my opinion, the check is necessary and should give an error, not only a warning: the function gives wrong results. The same applies of course to findorbit4 and findsyncorbit. Note however that, since atlinopt4 calls findorbit4, this results in a double check, with a scan of the whole lattice each time. Not a negligible time.

We can raise the question again: should we protect these functions against wrong use? and wait for comments…

My proposal a while ago was to set a warning (as you did) for some time, and then turn it later into an error. This leaves some time to upgrade user software. But I think that even the warning was considered too disturbing at that time.

If we go that way, a cleaner solution is to add a keyword to check_6d, like strict=true by default, and if strict is set to false, the function issues a warning instead of an error.

[oscarxblanco]

Dear @lfarv , I have added the strict flag. Could you have a look again ?

[swhite2401]

I add @simoneliuzzo as reviewer for this one

[simoneliuzzo]

Dear @oscarxblanco, @lfarv, @swhite2401 ,

I have had very big troubles with the 4D/6D evolution of AT. It simply resulted in me abandoning matlab AT for new projects in favor of pyAT. I consider my old code deprecated. It has served its time and can now be discontinued. I doubt I will ever resume my matlab software, rewriting it would be easyer. This is what I have started doing in the last years in pyAT (I am at about 15%). Luckily the wrong results I produced for many years were not sufficiently wrong to drive wrong conclusions. With pyAT I have learnt the habit to always create new environments, this is a solution to assure myself that I will have a copy of the code I used for a specific purpose frozen forever (at the detriment of disk space).

So I have no objections to adding the strict check and I approve this pull request without testing it.

For @oscarxblanco, in case you are running matlab simulations for tolerances/robustness/commissioningsimulations, please recall that there is a python version of the SC toolkit by Thorsten Hellert developed at DESY (https://github.com/lmalina/pySC) by Lukas Malina.

ciao simone

[swhite2401]

Hello all, if @simoneliuzzo is now ready for the switch I propose to turn the warning into an error at the next release, maybe this could be mentioned in the warning message? What do you think?

[simoneliuzzo]

@swhite2401 @lfarv @oscarxblanco ,

ok for me to switch to error.

ciao Simone

[lfarv]

@oscarxblanco: I would put the "strict" keyword in the check_6d function itself. It's probably not necessary set it accessible to the user: following @swhite2401's proposal:

• In a 1st step, strict would be hard-coded to False in the check_6d call in findorbit4, with a message stating that in the future this wrong behaviour will be forbidden,
• In a second step, strict will be changed to True, so that wrong manipulations are forbidden.

Another question is to chose between protecting atlinopt4, or findorbit4, or both. It's a matter of performance: atlinopt calls findorbit at least 3 times, possibly more depending on options: protecting findorbit is more expensive ! To decide, I'd like to see a measurement of the impact.

[lfarv]

The full solution would be to add to findorbit4 (and also findsyncorbit) a 3-valued argument:

• -1 meaning "no check", for a fast call from atlinopt4 or similar, where the check_6d would have already been done,
• 0: loose check,
• 1: strict check.

-1 would be for internal use only, skipping the call to check_6d, 0 and 1 would be sent as a 'strict' argument to check_6d, which interprets it as false or true and raises either a warning or an error.

Example:

function atlinop4(ring, ...)
...
check_6d(ring, false,  'strict', true);  % or possibly strict=false temporarily
...
[...] = findorbit4(..., 'strict', -1);
...
[...] = findorbit4(..., 'strict', -1);
...
end

function findorbit4(ring, ...)
...
[strict, varargs] = getoption(varargs, 'strict', true);  % or possibly strict=false temporarily
...
if strict >= 0
check_6d(ring, false, 'strict', strict);
end
...
end

function check_6d(ring,...)
[force, varargs]=getflag(varargin,'force');
[strict, varargs] = getoption(varargs, 'strict', true);  % or possibly strict=false temporarily
is_6d=atGetRingProperties(ring,'is_6d');

if nargin > 1
    if xor(is_6d,enable)
        if force
            if enable
                ring=atenable_6d(ring);
            else
                ring=atdisable_6d(ring);
            end
        elseif strict
            error(...)
        else
            warning(...);
        end
    end
end
...
end

This way, we can protect most functions while minimising the penalty in speed. But it's more complicated to set up. What is your advice?

[oscarxblanco]

Dear @lfarv , I have rewritten check_6d based on your comments but quite different to what you showed above. Let me know if it would be acceptable, I tried to separate the flags from the actions on the ring.

With respect to the speed test, I did a simple loop calling findorbit4 with and without the 6d check. It is almost two times slower to do the check.

tic;
for i =1:1000
    findorbit4(RING,'strict',-1); % skip check
end
toc;
%%
RING4D = atdisable_6d(RING);
tic;
for i =1:1000
    findorbit4(RING4D,'strict',0);% do the check
end
toc;
Elapsed time is 7.123552 seconds. % skip check
Elapsed time is 11.762007 seconds. % do the check

With respect to the findsyncorbit I have doubts because I don't use it, should it be set to 6d by default, right ?

[lfarv]

@oscarxblanco: Your test confirms what I feared: the check is very expensive. Since findorbit4 is used in many functions, and often repeatedly, I think that simply adding a systematic check in findorbit4 is too penalising. So we are left with 2 solutions:

1. Put a check in atlinopt4 (and similar), and leave findorbit4 unprotected. The relative time increase in atlinopt4 should be less and hopefully acceptable. To be checked,
2. Introduce the "skip test" option described above in findorbit4. For atlinopt4, the timing will be exactly the same as in option 1, but findorbit4 will be protected in case it's called directly. To avoid unexpected slowdowns, all functions calling findorbit4 should also be identified and updated. This is more difficult…

Concerning findsyncorbit, it's exactly like findorbit4, it must be used only in 4d. findorbit4 solves for a given dp/p, while findsyncorbit solves for a given path lengthening, like what you get with a RF frequency change.

[lfarv]

@oscarxblanco : I have 2 comments on check_6d:

• the way you implement it, when force is true, it will switch the ring to the desired state even if it is already initially as desired: this is unnecessary
• personally I prefer that the strict keyword is given as a logical (true or false): it makes more sense. If it is given a numerical value, it will be implicitly interpreted according to the standard rule: 0 => false, anything else => true. But that's a matter of taste, you may disagree…

[oscarxblanco]

Dear @lfarv , I have modified check_6d to enable/disable 6d only if it is different to the state of is_6d. The speed test on findorbit4 shows that it is better but still slower by 40% wrt to having no check.

Elapsed time is 7.856496 seconds. % no check
Elapsed time is 10.870588 seconds. % with check

I will try to reduce this difference.

[lfarv]

I will try to reduce this difference.

I don't think you can make it significantly better. The only way is to skip the check when it's not necessary, because has been done earlier (see bullet 2 above)

[lfarv]

@oscarxblanco : nice job !

I found 1 error, in pytests.m

[~,o0]=findorbit4(r4,dp=dpp,'strict',-1);

should be either:

[~,o0]=findorbit4(r4,dp=dpp,strict=-1);

or:

[~,o0]=findorbit4(r4',dp',dpp,'strict',-1);

The new key=value syntax may be used here because this function can only be used with recent Matlab versions implanting the test framework. Elsewhere, we try to maintain the compatibility with old Matlab versions which do not support the syntax.

[lfarv]

To conclude, you must now protect atlinopt4 with check_6d(ring, false, 'strict', 0) (Let's start with a warning!), replacing the commented check_radiation

Looking at your modified files, I found also mcf and findm44 which should be protected the same way. But that may show unexpected warnings in existing (wrong) software. We got @simoneliuzzo's approval for atlinopt4, not for the others, so let's restrict to atlinopt4.

Ideally, findsyncorbit should follow exactly findorbit4

[simoneliuzzo]

Dear @lfarv and @oscarxblanco,

I approve any further modification concerning matlab. I am not concerned anymore. The wound was deep but it is now sealed by python replacement.

ciao Simone

[lfarv]

@simoneliuzzo : too bad, you are leaving me alone… We need more Matlab users ! @oscarxblanco ,I rely on you !

[simoneliuzzo]

Dear @lfarv, now that you provided also load_madx and load_mad8, I really lost all need for matlab!

[oscarxblanco]

@simoneliuzzo : too bad, you are leaving me alone… We need more Matlab users ! @oscarxblanco ,I rely on you !

Dear @lfarv , I will be using both AT matlab and pyat for looong time.

[oscarxblanco]

Dear @lfarv , I think I have covered the modifications required. I did again the speed test and results show that using check_6d inside findorbit4 would take around 40% more time, I only managed to improve a little bit the speed wrt to the previous speed test.

Elapsed time is 10.106664 seconds. % findorbit4, 5000 times, strict = -1
Elapsed time is 13.780269 seconds. % findorbit4, 5000 times, strict = 0
>> 13.78/10.11
ans =
    1.3630

I attach a file I edited to do tests. test_multiple_calls.m.txt

With respect to mcf and findm44 I have modified the function documentation to point out that the calculation is wrong for 6-d rings.

[oscarxblanco]

@lfarv , is this check implemented in pyat ?

[lfarv]

@lfarv , is this check implemented in pyat ?

Yes pyat is protected since day 1 ! But to make it fast, the 6d status is stored in the lattice object, so that no scan of the lattice is necessary. This has other drawbacks: it is possible to modify an element without the lattice being aware, so bypassing the Lattice enable_6d/disable_6d results in wrong information…

[oscarxblanco]

@lfarv , is this check implemented in pyat ?

Yes pyat is protected since day 1 ! But to make it fast, the 6d status is stored in the lattice object, so that no scan of the lattice is necessary. This has other drawbacks: it is possible to modify an element without the lattice being aware, so bypassing the Lattice enable_6d/disable_6d results in wrong information…

ok, understood. Thank you @lfarv