Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

VariableMultipole compatibility between matlab AT and pyat #839

Open
wants to merge 73 commits into
base: master
Choose a base branch
from

Conversation

oscarxblanco
Copy link
Contributor

This pull request address issue #838 .
The pyat element class VariableMultipole is modified.
The attribute to_dict is checked in order to determine if the element is already created, if it has not been created the attributes are assigned as before.
This also fixes a bug in the determination of MaxOrder. Before it could be None, now it is limited to be at least zero.

@oscarxblanco oscarxblanco added Python For python AT code bug fix labels Oct 3, 2024
@oscarxblanco oscarxblanco added the WIP work in progress label Oct 3, 2024
@ZeusMarti
Copy link
Contributor

Not totally related, but I think it is also an small bug in In Matlab, I get this weird behavior:

KIDDK = atvariablemultipole('KIDDK','arbitrary','VariableThinMPolePass',...
'AmplitudeB',[0 0 0 0],...
'AmplitudeA',[0 0 0 0],...
'FuncA',[1 0],...
'FuncB',[1 0])

gives:

struct with fields:

   FamName: 'KIDDK'
PassMethod: 'VariableThinMPolePass'
    Length: 0
     Class: 'VariableMultipole'
      Mode: 2
  PolynomA: []
  PolynomB: []
AmplitudeB: [0 0 0 0]
AmplitudeA: [0 0 0 0]
     FuncA: [1 0]
     FuncB: [1 0]
  MaxOrder: -2
 NSamplesA: 2
 NSamplesB: 2

If MaxOrder is fixed to a value then it is ok.

@oscarxblanco
Copy link
Contributor Author

Dear @ZeusMarti ,
I confirm it. The same weird behaviour is in pyat, I'm working on a solution in pyat compatible with matlab. Once I have it I will fix the bug in matlab.

@oscarxblanco oscarxblanco removed the WIP work in progress label Oct 4, 2024
@oscarxblanco
Copy link
Contributor Author

Dear @swhite2401 ,
who could review this pull request ?

@oscarxblanco oscarxblanco changed the title only create attributes if self is not an element VariableMultipole compatibility between matlab AT and pyat Oct 4, 2024
@swhite2401
Copy link
Contributor

@oscarxblanco , I can do it. In fact there are some changes that you did that I don't understand.
Also I like that important input argument are not integrated in kwargs so that they appear in the API doc. I may propose changes to the __init__()

Copy link
Contributor

@swhite2401 swhite2401 left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Hi @oscarxblanco , here is a first round of comments

Periodic=bool,
)

def __init__(self, family_name: str, **kwargs: dict[str, any]):
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

You don't need to specify the type of kwargs this is standard python

Could you please restore the defaults for AmplitudeA/B and the mode? I think it is better if they appear in the API doc, in fact I would be tempted to add more of these keyword arguments rather than removing them

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Dear @swhite2401 , if I remove the annotation I get the following flake8 warning. This is the reason why I included, should I ignore it ? I may be using flake8 with a wrong configuration or I may use another tool ?

/variable_elements.py:42:44: ANN003 Missing type annotation for **kwargs

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Could you please restore the defaults for AmplitudeA/B and the mode? I think it is better if they appear in the API doc, in fact I would be tempted to add more of these keyword arguments rather than removing them

Dear @swhite2401 , the reason why I remove it is because of an incompatibility between AT matlab and this implementation of the element in pyat.
If I save a lattice in matlab and try to read it in pyat I get the error here below, I'll try to explain the reason in detail, but, in few words: one cannot have args and kwargs with the same name

When reading an element from matlab it is initially put in a dictionary elem_dict, and then the _BUILD_ATTRIBUTES are stripped to pass them to the class instantiation as elem_args. See the lines below

at/pyat/at/load/utils.py

Lines 368 to 369 in 268994f

elem_args = [elem_dict.pop(attr, None) for attr in cls._BUILD_ATTRIBUTES]
element = cls(*(arg for arg in elem_args if arg is not None), **elem_dict)

I could not put the AmplutudeA, AmplitudeB and Mode in the element BUILD_ATTRIBUTES because they are not garanteed to exist in matlab, see the following line

function elem=atvariablemultipole(fname,varargin)

therefore, the call to cls with args in necessarily incomplete and the AmplitudeA, AmplitudeB and Mode are passed in a dictionary together with args that have the same name. Python does not handle this mix well, it removes the elements from the kwargs and leaves unset the args.

There are ways to workaround this situation when args and kwargs have the same name, but they require changing the way the class is instantiated in utils.py. I decided not to do it because that will change how all elements are instantiated.

A second possibility was to rename the element args in pyat, for example amplitudeA instead of AmplitudeA i.e. lower case in the first letter to make one small change. However, this would break compatibility.

I decided to keep the element without any requirement as in AT matlab, which forced me to do checks on the input parameters before assigning them inside __init__. This wouldn't break compatibility with previous versions of pyat, and AT matlab files could be read because the element also has no requirements except for the FamName.

If you prefer, I could do as you asked. Keep the args and break backwards compatibility. Or force the user to provide AmplitudeA, AmplitudeB and Mode in matlab, which will make possible to include them in BUILD_ATTRIBUTES that will be removed in the input dictionary and set before the class is instantiated.

Or maybe there is another solution ?

>>> ring = at.load_mat('test1.mat')
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "/home/oblanco/Documents/public/progs/pyenv/pyenvat/lib/python3.11/site-packages/at/load/matfile.py", line 202, in load_mat
    return Lattice(
           ^^^^^^^^
  File "/home/oblanco/Documents/public/progs/pyenv/pyenvat/lib/python3.11/site-packages/at/lattice/lattice_object.py", line 191, in __init__
    super(Lattice, self).__init__(elems)
  File "/home/oblanco/Documents/public/progs/pyenv/pyenvat/lib/python3.11/site-packages/at/lattice/lattice_object.py", line 1478, in params_filter
    for elem in elem_filter(params, *args):
  File "/home/oblanco/Documents/public/progs/pyenv/pyenvat/lib/python3.11/site-packages/at/load/matfile.py", line 140, in ringparam_filter
    for elem in elem_iterator(params, *args):
  File "/home/oblanco/Documents/public/progs/pyenv/pyenvat/lib/python3.11/site-packages/at/load/matfile.py", line 105, in _matfile_generator
    yield Element.from_dict(kwargs, index=index, check=check, quiet=quiet)
          ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/home/oblanco/Documents/public/progs/pyenv/pyenvat/lib/python3.11/site-packages/at/load/utils.py", line 365, in element_from_dict
    element = cls(*(arg for arg in elem_args if arg is not None), **elem_dict)
              ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/home/oblanco/Documents/public/progs/pyenv/pyenvat/lib/python3.11/site-packages/at/lattice/variable_elements.py", line 111, in __init__
    raise AtError("Please provide at least one amplitude for A or B")
at.lattice.utils.AtError: Please provide at least one amplitude for A or B

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

A second possibility was to rename the element args in pyat, for example amplitudeA instead of AmplitudeA i.e. lower case in the first letter to make one small change. However, this would break compatibility.

Hi @oscarxblanco , thanks for the detailed answers! My personnal opinion is that this was poorly implemented to start with (which is my fault in fact, sorry...) and that in such situation we should be allowed to break backward compatibility.
In all other pyAT element, args seem to have in fact different name that kwargs and we should probably stick to that convention to keep things simple. Similarly we could modify MATLAB to preserve compatibility with the new pyAT version.

This feature is relatively recent and it could well be that @ZeusMarti is the only user... we have used it at ESRF but I don't think we have saved lattice files including it. Would it be a big effort on your side to switch to a new signature?

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Dear @swhite2401 , I asked @ZeusMarti and he was not able to use it effectively because he intends to have a custom function where the time of the particle is also considered, i.e. not only a kick every turn for all particles.
I have started modifying the C code today to make some test, is this feature something you would accept for this element ? Or should I keep the custom function option as it is ?

o

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I have started modifying the C code today to make some test, is this feature something you would accept for this element ?

Sure! In fact I developed this for @ZeusMarti initially, so any improvement is perfectly fine.
You can already do a kick every n turns by using the table by the way, just fill zeros where you do not want any kick.

Do you want help for this?

super(VariableMultipole, self).__init__(family_name, **kwargs)

def _setmaxorder(self, ampa, ampb):
if len(kwargs) > 0:
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Why is this needed?

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This was needed in order to separate the case when the element is called inside pyat, therefore the object already exists and there are no kwargs, from the case when the element is instantiated using a dictionary (the dictionary comming from reading a .mat, or from user settings).

Sometimes the object is read, for example when printing in a terminal, the element is passed again (I guess this happens because it inherits __str__ from other element) and no modifications in the object are needed, I just need to pass it to super and find its way to __str__.


def _setmaxorder(self, ampa, ampb):
if len(kwargs) > 0:
self.FamName = family_name
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

family_name = kwargs.pop('FamName', family_name)

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I explained before, but, when the element comes from AT matlab, a dictionary is created but all BUILD_ATTRIBUTES are stripped. 'FamName' is a a BUILD_ATTRIBUTE of the class element and therefore it is not in kwargs. One has to take it from args.

raise AtError("Please provide at least one amplitude for A or B")
# start setting up Amplitudes and modes
# fist modes are called differently
modepyatinput = kwargs.pop("mode", ACMode.SINE)
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

You could replace all these lines by:

Self.Mode = kwargs.pop('Mode', int(mode))

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I cannot.
The issue here is that the argument is called 'mode' in pyat and 'Mode' in MATLAB.
I need to check in advance which is correct.

def _setmaxorder(self, ampa, ampb):
if len(kwargs) > 0:
self.FamName = family_name
if "AmplitudeA" not in kwargs and "AmplitudeB" not in kwargs:
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

AmplitudeA/B présence in kwargs also needs to be checked

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Sorry @swhite2401 , I don't understand. What you say is implemented in line 103 above.
You mean individually ? I thought this element should accept to have at least one of the two.

self.MaxOrder = kwargs.get("MaxOrder", 0)
amplitudea = None
amplitudeb = None
if "AmplitudeA" in kwargs:
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This lol ne and the ones below can be removed

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

What do you mean by lol ?

# fill in super class
super().__init__(family_name, **kwargs)

def _setmaxorder(self, ampa: np.ndarray, ampb: np.ndarray):
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

In what cases did this end up to be none? When ampa/b was filled with only zeros?

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I have to repeat the test but suppose you have the current implementation available in the master. You will have this behaviour with zeros.

>>> numpy.max(numpy.nonzero([0,0,0,0]))
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "/home/oblanco/Documents/public/progs/pyenv/pyenvat/lib/python3.11/site-packages/numpy/core/fromnumeric.py", line 2810, in max
    return _wrapreduction(a, np.maximum, 'max', axis, None, out,
           ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/home/oblanco/Documents/public/progs/pyenv/pyenvat/lib/python3.11/site-packages/numpy/core/fromnumeric.py", line 88, in _wrapreduction
    return ufunc.reduce(obj, axis, dtype, out, **passkwargs)
           ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
ValueError: zero-size array to reduction operation maximum which has no identity

I changed to mxa = np.max(np.append(np.nonzero(ampa), 0)) which is the maximum between the nonzero in ampa and zero, and I get sure ampa exists before doing it.


def _set_params(self, amplitude, mode, ab, **kwargs):
def _set_params(
self, amplitude: int or str, mode, a_b: str, **kwargs: dict[str, any]
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

No need to define the type of kwargs

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

It is flake8 who suggested. Should I ignore it ? Do you work with some custom flake8 settings for pyat development ?

@oscarxblanco oscarxblanco added the WIP work in progress label Nov 18, 2024
@oscarxblanco oscarxblanco added Matlab For Matlab/Octave AT code C For C code / pass methods and removed WIP work in progress labels Dec 20, 2024
@oscarxblanco
Copy link
Contributor Author

Dear @swhite2401 , I think this is ready for a second review.
Who should review it ?

I force the user to define name and mode.
The use of one amplitude value or array for A or B is checked at the element setup.
The MaxOrder option has been removed as it seemed to be overwritten in most cases by the element setup.
The seed option has been removed. See this discussion #879.
The modes SIN and ARBITRARY now consider the particle time delay. Additionally a Sinlimit option has been added to create a half-sin when setting it to zero.
Mode WHITENOISE, ignores the particle time delay. MEAN and STD have been added as options with default values 0 and 1.
Ramping has not been modified.
In mode ARBTRARY an additional time delay has been added as an option, default value zero. Also, the array of turn by turn values was extended to the 4th Taylor expansion wrt $\tau$ (the time delay, or the sixth coordinate divided by the speed of light).
Periodic is now only used in ARBITRARY mode.

I have done test creating the elements in matlab, saving and running in python. Also, creating the elements in python, saving and runing in matlab. And also, creating and running the elements on the same environment they are created. They all work, but as the number of options is large because the pass method combines three elements in one, I have checked only kicks and skew quadrupole components with varying modes.

@swhite2401
Copy link
Contributor

I can review it but I will be away for 2 weeks starting tomorrow.

Discussing with colleague here at ESRF we thought it would be nice that this element takes into account the time delay in multi-bunch mode (for instance to simulate an injection kicker pulse), I think this is possible with the present implementation since the multi-bunch beam is unfolded by default (see tracking/track.py) but it would be nice to confirm it.

I provide a review next year... unless @lfarv would like to take a look earlier?

@oscarxblanco
Copy link
Contributor Author

I'm not sure how multi bunch simulations are done, but, I have heard something similar from Gabriele at ALBA. So, I'll have a look.

I am also leaving today, and I'll be back the 7th of Jan/2025.

Happy hollydays.
o

@swhite2401
Copy link
Contributor

I can send you an example, happy holidays!

Comment on lines +237 to +242
def _set_white_noise(self, a_b: str, **kwargs: dict[str, any]) -> dict[str, any]:
if "Mean" + a_b not in kwargs:
kwargs["Mean" + a_b] = kwargs.get("Mean" + a_b, 0)
if "Std" not in kwargs:
kwargs["Std" + a_b] = kwargs.get("Std" + a_b, 1)
return kwargs
Copy link
Contributor

@lfarv lfarv Dec 23, 2024

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think that Mean and Std do not bring anything : the static part "Amplitude*Mean" can be provided by a standard ThinMultipole element, and for the dynamic part, Amplitude and Std are redundant. I would simplify both python and C by removing these attributes: Amplitude is sufficient.

Comment on lines +203 to +204
if "Amplitude" + key in kwargs and "amplitude" + key in kwargs:
raise AtError("Duplicated amplitude " + key + "parameters.")
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Why are you looking for "amplitudeA/B" if the attribute name is AmplitudeA/B" ?

'Ramps has to be a vector with 4 elements'
self.Ramps = ramps
super(VariableMultipole, self).__init__(family_name, **kwargs)
kwargs["Mode"] = kwargs.get("Mode", mode)
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

better written as:

kwargs.setdefault("Mode", mode)

Periodic=bool,
)

def __init__(self, family_name: str, mode: int, **kwargs: dict[str, any]):
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

"any" is undefined.

If you really want to give a type hint for kwargs, which is unusual, it should be dict[str, Any], with:

from typing import Any

kwargs["Ramps"] = ramps
super().__init__(family_name, **kwargs)

def _check_amplitudes(self, **kwargs: dict[str, any]) -> dict[str, any]:
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

All the following methods can be defined as static methods, or nested in init

Comment on lines +248 to +249
kwargs["Phase" + a_b] = kwargs.get("Phase" + a_b, 0)
kwargs["Sinlimit" + a_b] = kwargs.get("Sinlimit" + a_b, -1)
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Same comment as above:

kwargs.setdefault("Phase", 0.0)
kwargs.setdefault("Sinlimit", -1)

raise AtError("Please provide at least one amplitude for A or B")
return amp_aux

def _set_amplitude(self, amplitude: float or _array or None) -> _array:
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

_array is a function:

def _array(value, shape=(-1,), dtype=numpy.float64): ...

You cannot use it as a type. Here, the type could be float | Sequence[float]

Comment on lines +213 to +216
def _set_amplitude(self, amplitude: float or _array or None) -> _array:
if np.isscalar(amplitude):
amplitude = [amplitude]
return np.asarray(amplitude)
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Anyway the whole function is useless: when setting the Amplitude attribute, the conversion function (_array) will convert the value into a Fortran-aligned array of floats with shape (-1,) which is exactly what you are doing here.

@lfarv
Copy link
Contributor

lfarv commented Dec 23, 2024

@oscarxblanco : the documentation of the VariableMultipole class is a bit strange. Sphinx is tricky with indents and formatting…

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
bug fix C For C code / pass methods Matlab For Matlab/Octave AT code Python For python AT code
Projects
None yet
Development

Successfully merging this pull request may close these issues.

4 participants