Variable multipole

[swhite2401]

This PR is a follow-up of #503.
It introduces the possibility to have time dependent multipoles using 3 different excitation modes:

• sine excitation
• gaussian white noise
• user defined turn-by-turn kick list (repeated)

A new passmethod was added together with element creation for both matlab and python.

[swhite2401]

Below an python example script

import numpy
import matplotlib.pyplot as plt
import at

path = '../lattice/'
filename = 'S28F.mat'
key = 'LOW_EMIT_RING_INJ'
latticef = path + filename
ring = at.load_lattice(latticef, key=key)
ring.radiation_off(cavity_pass='CavityPass')
tunes = ring.get_tune()
frev = ring.get_revolution_frequency()

print(tunes)
print(frev)

fring = at.fastring._fring(ring)

order = 2
ampb = 1e-5;
ampa = [0.0, 1e-5];
frequency = (0.01) * frev
#acmpole = at.VariableMultipole('ACMPOLE', AmplitudeB=amp, FrequencyB=frequency)
#acmpole = at.VariableMultipole('ACMPOLE', AmplitudeB=amp, mode=at.ACMode.WHITENOISE)
acmpole = at.VariableMultipole('ACMPOLE', AmplitudeB=ampb, FuncB=[0.0, 1.0, 0.0], mode=at.ACMode.ARBITRARY)
#acmpole = at.VariableMultipole('ACMPOLE', MaxOrder=order,
#                               AmplitudeB=ampb, FrequencyB=frequency,
#                               AmplitudeA=ampa, FrequencyA=frequency)

fring.append(acmpole)
print(acmpole)
xout = numpy.squeeze(at.tracking.lattice_pass(fring, numpy.zeros(6), nturns=998))
print(acmpole)
plt.plot(xout[0, :])
plt.show()

[swhite2401]

Same for matlab

clear;

v = load('./S28F.mat');
ring = v.LOW_EMIT_RING;
ring = atradoff(ring,'CAVIPASS','RFCavityPass');
frev = atGetRingProperties(ring, 'revolution_frequency');
tunes = tunechrom(ring);

freq = 0.01*frev;
amp = 1e-5;

acm = atvariablemultipole('ACM','AmplitudeB',amp,'FrequencyB',freq);
acm = atvariablemultipole('ACM','AmplitudeB',amp,'Mode','WHITENOISE');
acm = atvariablemultipole('ACM','AmplitudeB',amp,'Mode','ARBITRARY','FuncB',[0.0 1 0.0]);

ring = [ring;{acm}];

rout = squeeze(ringpass(ring,zeros(6,1),998));
plot(rout(1,:))
'''

[swhite2401]

I have tried to include all the feature that were discussed in #503 , but of course this is all up for this discussion.
Please take a careful look at the matlab part, I am really not comfortable with it... thanks!

[lfarv]

There is a problem in the Matlab element creation:

atvariablemultipole('ACM','AmplitudeB',1.e-4,'FrequencyB',1.e3);

works, but

atvariablemultipole('ACM','SINE','AmplitudeB',1.e-4,'FrequencyB',1.e3);

crashes. decodeatargs assumes that positional arguments (those before the PassMethod) are numeric. Here the "mode" is not.

I'll commit a modified version which decomposes decodeatargs in its 2 steps and accepts the mode. The Matlab help also mentions the SEED and RAMPS keywords which apparently are not used. I removed them temporarily from the help.

[lfarv]

This loop rewrites completely the PolynomA/B fields of the element. It is not recommended for a passmethod to modify its input. Can't we allocate locally pola and polb(atMalloc/atFree) rather than using the element data fields? Then PolynomA and PolynomB can be removed the element fields.

[swhite2401]

I did implement it as you suggest to start but them re-introduce PolynomA/B because it was useful to keep track of their values for debugging.
I thought that users may find it useful as well to have access to these values so I left them in. They can be remove of course.
Is there a reason why a passmethod cannot change its element fields? We have several such passmethods for collective effects.

[lfarv]

OK, no problem. The only danger, and that's why the element attributes should better be "read-only" in the integrator, is that if a user modifies (shortens, removes…) these attributes, the integrator may crash… Since PolynomA and PolynomB are not documented, there's a small risk. We can live with that !

[swhite2401]

How do you handle this in standard elements. if MaxOrder>length(PolynomA/B) is crashes too right?

[lfarv]

I do not see clearly how the Ramps field interacts with the "mode" (SINE,WHITENOISE,ARBITRARY). Is it a 4th mode?

Also, is it possible to define a single pulse (which does not repeat as ARBITRARY does) ?

[lfarv]

There are 3 different implementations of the gaussian random generator:

• 2 in QuantDiffPass.c (random_normal, not used, and GenerateGaussian)
• 1 in VariableThinMPolePass.c (randn). This one is very similar to GenerateGaussian, except that it does not keep the 2nd generated value for the next call. So it's about twice slower, but avoids a static variable.

It would be nice to have a single implementation, in an #included C file, available to any integrator.

[swhite2401]

I do not see clearly how the Ramps field interacts with the "mode" (SINE,WHITENOISE,ARBITRARY). Is it a 4th mode?

Also, is it possible to define a single pulse (which does not repeat as ARBITRARY does) ?

Ramp is a multiplicator that allows to set the start and end of the excitation as well as linear ramps up and down.
See double get_amp(double amp, double *ramps, double t) that is called in the polynom generation function

[swhite2401]

There are 3 different implementations of the gaussian random generator:

* 2 in `QuantDiffPass.c` (`random_normal`, not used, and `GenerateGaussian`)

* 1 in `VariableThinMPolePass.c` (`randn`). This one is very similar to `GenerateGaussian`, except that it does not keep the 2nd generated value for the next call. So it's about twice slower, but avoids a static variable.

It would be nice to have a single implementation, in an #included C file, available to any integrator.

I fully agree, and in fact there is a general problem with random number generation as you pointed out in #502 that we need to solve. I would suggest to have a full review in a separate PR

[swhite2401]

The Matlab help also mentions the SEED and RAMPS keywords which apparently are not used.

Why do you say they are not used?

[lfarv]

Ramp is a multiplicator that allows to set the start and end of the excitation as well as linear ramps up and down.
See double get_amp(double amp, double *ramps, double t) that is called in the polynom generation function

Ok, I realised later how it works. May be it is worth a little more documentation?

[lfarv]

The Matlab help also mentions the SEED and RAMPS keywords which apparently are not used.

Why do you say they are not used?

Because I didn't look carefully enough… Sorry

[lfarv]

Hi @swhite2401 ! Concerning the random generator;

I fully agree that the global question of random generation with parallel processing should be addressed in a dedicated PR. For now, my proposal is:

• Since a gaussian generator is part of this PR, let's move into a separate C file: The multiprocessing problem will be easier solved if there is a single implementation, and this still allows this branch to be merged rapidly. The other clients of gaussian generator may be converted to the "common" one later.
• The seeding done from element parameters is not fully satisfactory: if there are several elements, it's not obvious which seed will be taken (the 1st along the lattice, I guess). And since it's done once, if you want to "replay" the sequence, you have to restart python or Matlab (may be "clear mex" would work on Matlab). It would be convenient to have a direct python/Matlab interface to the srand C function. I can deal with that in a separate PR.
• The "common" generator should be as simple as possible: no argument, mean and stddev are introduced outside. The one from QuantDiffPass.c looks like a good candidate after some cleaning.

[swhite2401]

Hello @lfarv !
Agreed, the subtlety is that the noise generator needs the same seed in all threads while quantum diffusion is the opposite... so it is really a complex problem because combined the white noise excitation may not be compatible with quantum diffusion elements.

I will do what you propose and include an additional C file to host random generators (don't forget the poisson random number generator as well!) in this PR, then we can merge as is.

[ZeusMarti]

Hi there,

This is a very helpful implementation for us, I have been using the branch and all seems OK to me. I just have a couple of small comments:

1. I agree with @lfarv it would be nice to be able to select or not the preriodicity of the excitation, in the case of a multipole injection kicker for example, if no preriodicity is seected, once the excitation finishes you can just set the PolynomA/B to 0.
2. I like the fetature of the PolynomA/B modification of the element after a run, after all it is a variable element! It was strange at first but I kind of like it now, of course this is just my personal preference.

Cheers,

[swhite2401]

1. I agree with @lfarv it would be nice to be able to select or not the preriodicity of the excitation, in the case of a multipole injection kicker for example, if no preriodicity is seected, once the excitation finishes you can just set the PolynomA/B to 0.

Hi Zeus,
You can already do this by using the ramps input.
If you set ramps=[t1,t1,t2,t2] the kick list will be played only at turn t1 and then set back to 0 (I admit haven't tried but it should work)

[swhite2401]

With t2 = t1+1

[ZeusMarti]

Ah, OK that should work, although it is not completely obvious to me from the help section as it is, I think other users may find it helpful if this behavior is explicitly described.
In this case the kick calculation is done for every turn but after t2 the effect is multiplied by a 0 right? May not be the most computationally efficient thing to do but I guess the overhead is ridiculous, right?

[swhite2401]

Correct, then I will add another optional argument periodic=True, if set to False the kick list will be played only once from turn 0, this is indeed cleaner.

[lfarv]

Back to the random generator: there is a very simple solution for OpenMP. It consists in replacing rand() by erand48().

double erand48(unsigned short xsubi[3]);

This function stores its state not in an internal buffer like rand, but in an external one (xsubi). This buffer can be declared as a "shared" variable for OpenMP, so that all threads will use the same one (for quantum diffusion). For variable multipole, the number is generated outside the particle loop, so no problem.

erand48 has other benefits: it works on 48 bits, returns a double in the range [0, 1) which is convenient for our use and by using different buffers you create different and independent streams of random numbers. There is also the equivalent drand48() using an internal buffer (no argument).

These function are declared in <stdlib.h> and should be standard. We have to check that they are available on all platforms, but if it's ok, they would be the best source for the "common" random generator.

It does not solve the MPI problem, though the idea of one buffer per thread with a different initial value is an idea…

[swhite2401]

Hello, I have move the random number generator to a separate file atrandom.c and added an optional Periodic argument that allows to select whether the kick list is repeated or not.
Also I have update the help of both python and matlab.

Any other comments or suggestions?

[swhite2401]

Back to the random generator: there is a very simple solution for OpenMP. It consists in replacing rand() by erand48().

Sounds good, but don't you think that implementing the same method for MPI and OpenMP would be simpler? We can use the rank for both methods

[lfarv]

Looks ok for me

[lfarv]

Back to the random generator:

I think that a random generator with an external state variable can solve the problem for both OpenMP and MPI, but possibly not in the same way…

1st problem: erand48 is missing in Microsoft Visual C (MSVC). There is a workaround using a similar rand_s Microsoft function, or we can use our own generator applicable to all platforms. I found a good candidate (Apache license): PCG.

• for MPI, we need a different state variable for each process. It could be allocated at the startup of each process (but how to access it in the integrator?), or one could create a new state variable before starting each particle loop.This is ok for quantum radiation, but not for the variable multipole.
• for OpenMP, either we use a "shared" state variable, or similarly we create a new state variable per thread just before the loop/ May be the shared variable implies some lock mechanism in OpenMP which would be detrimental for performance…

[swhite2401]

I think that a random generator with an external state variable can solve the problem for both OpenMP and MPI, but possibly not in the same way…

Ok let's continue this on a dedicated PR, you initiate it?

[swhite2401]

Looks ok for me

Good, can I merge? @ZeusMarti any other suggestions on your side?

[ZeusMarti]

It is all good from my side, thanks!