pyat ReadWrite Insertion Devices into mat files

[oscarxblanco]

Dear all,

this PR is a workaround to a limitation in read and write capabilities on pyat when adding an Insertion Device map, already mentioned in #558 . The existing InsertionDeviceKickMap class always needed the field map file (here called, IdFile) to be read, and that made it incompatible with the save_mat function.

Here, the class InsertionDeviceKickMap has been modified to a function, and a new class called IdTable is created compatible with the structure stored in .mat files from pyat and matlab. Due to the non-existence of an equivalent class in AT matlab, the insertion device element is saved by changing it into a multipole which exists in both AT versions.

I expect no impact on the users, as the InsertionDeviceKickMap function uses exactly the same call and parameters, and the IdTable is saved in a backwards compatible format through the multipole class (no new elements are created on the .mat output).

The output to .m files is still not compatible. The output from ATmatlab transforms the element into a drift which could change the beam dynamics results when read back. I have done the same on pyAT to keep similar results on both sides, but, I discourage to use .m file for the moment.

The following diagram shows the data flow.

Best regards, o

sequenceDiagram
    participant IdFile
    participant pyAT
    participant .mat
    participant .m
    participant ATmatlab
    IdFile->> pyAT: InsertionDeviceKickMap (creates IdTable)
    pyAT ->> .mat: at.save_mat(IdTable to multipole)
    .mat ->>pyAT: at.load_mat
    IdFile->>ATmatlab:atidtable_dat
    ATmatlab->>.mat: save (VoidClass)
    .mat->>ATmatlab: load
    ATmatlab->>.m: atwritem(VoidClass to Drift)
    pyAT->>.m: at.save_m(IdTable to Drift)

[oscarxblanco]

Dear @lfarv ,

what is your opinion on this modifications ? Would you consider them ok to review ? Otherwise I could leave them as a draft for a longer discussion on how to pass new classes from one AT implementation to the other.

Best regards, o

[swhite2401]

Hello @oscarxblanco , let's wait for @lfarv to be back he may have good suggestions, but this seems very complicated to me.

I have a couple questions: -would it be useful to add a load_file() function and make the file optional in the constructor, such that in case it is missing you can still build the element? -it seems the kickmap is saved in the element attributes, in case they already exist (as is the case when loading the element from lattice file I suppose) why do you request the kickmap file to construct the element? In fact the passmethod does not need the file, just the kickmap arrays -in matlab can't we just write an element creation function to avoid changing the class to multipole? -for the .m files I think the same comments as above apply and this should work, however the full kickmap will be displayed but that is ok I believe no?

Cheers

[oscarxblanco]

Dear @swhite2401

Hello @oscarxblanco , let's wait for @lfarv to be back he may have good suggestions, but this seems very complicated to me.

I have a couple questions: -would it be useful to add a load_file() function and make the file optional in the constructor, such that in case it is missing you can still build the element?

Sure it could be an option that I have not tried out yet.

-it seems the kickmap is saved in the element attributes, in case they already exist (as is the case when loading the element from lattice file I suppose) why do you request the kickmap file to construct the element?

When starting from python it should have some kind of class to be saved. And you are right, the attributes are saved but the mapping from pyat classes to AT matlab elements is still missing.

In fact the passmethod does not need the file, just the kickmap arrays -in matlab can't we just write an element creation function to avoid changing the class to multipole?

Yes, I would also prefer this solution. But, I wrote this work around to avoid changing both pyat at AT matlab at the same time. This of course is possible.

-for the .m files I think the same comments as above apply and this should work, however the full kickmap will be displayed but that is ok I believe no?

Unfortunately not yet, there is a problem when reading back the array from an m file. The memory allocation is not fortran-aligned. I believe this could be solved but I wanted to raise the point before doing more modifications.

Cheers

[swhite2401]

My feeling is that we should go for an option that does not require external files to load a lattice file. I you want I can help with the implementation (after IPAC).

[oscarxblanco]

My feeling is that we should go for an option that does not require external files to load a lattice file. I you want I can help with the implementation (after IPAC).

Perfect ! Thank you

[oscarxblanco]

Quick update: This PR is still WIP, but, there is a small update wrt to pyat.

After conversations with @swhite2401 and @simoneliuzzo , the initial proposal to store the ID data in a multipole class is abandoned because it was too misleading.

In tag https://github.com/atcollab/at/pull/597/commits/bd4803854923a616decf7534cef82a44f924db8c An element with the class name idtable is created, and I already tested writing to and reading from, mat and m files on pyat. There is still a problem with the Energy parameter that is automatically ignored when reading and therefore needs to be included manually every time the lattice is read.

So, I'm half way through. The pyat modifications work. The matlab implementation of a compatible atidtable is missing. I will solve the issue wrt to matlab, and at the end try to solve the conflict on pyat/at/load/utils.py shown in the automatic check.

o

[swhite2401]

Thanks for the update, I have not arrived to this issue yet but I will soon be able to help if needed!

[oscarxblanco]

Dear all, I finally manage to save and load the insertion devices data into .mat and .m file formats for both pyat and AT matlab, obtaining the same tune variation from each file and code.

A new class called KickMap has been created . It contains the info required by IdTablePass plus additional info to keep track of the filename and normalization energy for the insertion device table. This implementation only reads two tables from the insertion device file. Others are not implemented yet because I don't have a working example with more than two tables per file.

The functions to create the insertion device have not changed their name: InsertionDeviceKickMap on pyat and atidtable_dat on matlab. They both use the KickMap class.

Other functions in matlab like, idtable, and idtable_global have not been modified and seemed to support more than two tables per insertion device file. Also, the function atidtable seems not to work properly.

I can provide test examples for the review of this PR. If tests should be integrated to the PR I would need some guidance about how it is prefered to do so as this calls pyat and matlab functions.

@swhite2401 Do you know who could review this PR ? Or would you prefer/suggest to modify it ?

sequenceDiagram
    participant IdFile
    participant pyAT
    participant .mat
    participant .m
    participant ATmatlab
    IdFile->> pyAT: InsertionDeviceKickMap (creates KickMap class)
    pyAT ->> .mat: at.save_mat(Kickmap class)
    .mat ->>pyAT: at.load_mat (KickMap class)
    IdFile->>ATmatlab : atidtable_dat (creates KickMap class)
    ATmatlab->>.mat: save (KickMap class)
    .mat->>ATmatlab: load (KickMap class)
    ATmatlab->>.m: atwritem(KickMap class)
    .m->>ATmatlab: evaluate function from file (KickMap class)
    pyAT->>.m: at.save_m(KickMap class)
   .m->>pyAT: at.load_m(KickMap class)

[swhite2401]

@oscarxblanco this looks much better indeed. I have a first question: why in pyAT are you using a function to instantiate the element class? Can't it be all handled in the class?

[swhite2401]

I think @lfarv should take a look as well

[oscarxblanco]

@oscarxblanco this looks much better indeed. I have a first question: why in pyAT are you using a function to instantiate the element class? Can't it be all handled in the class?

@swhite2401 I thought at that time that it was better to make a function that creates the class from a file, and make it available to the user. Having a class that contains that functions might be also possible but have not done it that way. Do you think that would be better ?

[swhite2401]

Not necessarily, but pyAT users are used to select ring by class for example: ring.get_unit32_index(at.Dipole) would return all dipoles. at.Dipole is the element class and is used to instantiate the element.

In you implementation one would then naturally try: ring.get_unit32_index(at.InsertionDeviceKickMap) while the correct way would be ring.get_unit32_index(at.KickMap)

This could be confusing, but maybe not a huge problem...

[oscarxblanco]

Seems important to follow the pyat logic. I will implement the changes.

[oscarxblanco]

@swhite2401 while working on the file saving for this PR, I got stuck in an issue it took me some time to understand. Sextupoles in pyat are read and then stored in memory with an additional parameter 'H'. When saving an .mat file using save_mat in pyat and reading in matlab, this 'H' is preserved. If the lattice is saved to an .m file using atwritem in matlab, 'H' is written as one of the sextupole properties. atsextupole('S1',1e-08,150500000,'H',150500000);...

However, load_m in pyat does not support the 'H' argument and throws the following error. super(Element, self).__setattr__( File "/home/sources/physmach/blanco-garcia/codes/pyenv/venv_pyat_from_sources/lib/python3.10/site-packages/at/lattice/elements.py", line 711, in H self.PolynomB[2] = strength AttributeError: In element S1, parameter H: 'Sextupole' object has no attribute 'PolynomB'

Do you think that it would be better to never store 'H' ? or to modify load_m in pyat so that 'H' is supported ?

[lfarv]

@oscarxblanco The .mat and .m formats defined by the Matlab specificities, and pyAT must mimic this behaviour.

in Matlab, both H (or K) and PolynomB are standard attributes, even if it's redundant. So in .mat files, both are dumped. For .m files, it's different, the text output mimics a call to the constructor of the element. H (or K) are positional arguments of the constructor and appear in atwritem as the 3rd argument. PolynomB may appear as a keyword argument, only if it is different from what H (or K) generate. But H should not. So what you report looks like a bug in atwritem.

However I cannot reproduce this behaviour:

>> sext = atsextupole('S6',0.4,-3.6,'NumIntSteps',20)

sext = 

  struct with fields:

        FamName: 'S6'
     PassMethod: 'StrMPoleSymplectic4Pass'
         Length: 0.400000000000000
          Class: 'Sextupole'
       PolynomB: [0 0 -3.600000000000000]
    NumIntSteps: 20
       PolynomA: [0 0 0]
       MaxOrder: 2

>> atwritem({sext})
ring={...
    atsextupole('S6',0.4,-3.6,'NumIntSteps',20);...
};

This is correct, only the non-default parameters appear as keyword arguments. I Tried on several lattices and always got the correct output. Could you give more details on the way to reproduce the bug?

[lfarv]

@oscarxblanco : I see that you modified at2str. However, what you added should not have any impact, so it must be something else…

[oscarxblanco]

@lfarv It is the result of saving .mat file in pyat, reading it in matlab, writing it to .m, and trying to read the .m file in pyat. Here is the sequence that reproduces the error (you will need to edit the CODEDIR environment variable)

python -c "import at; s1 = at.Sextupole('S6',0.4,-3.6); ring = at.Lattice([s1],energy='1'); at.save_mat(ring,'pyatsaveSext.mat');"
export CODEDIR=<mypathtoATmatlab>
matlab -nosplash -nodisplay -r "clear all; setatversion;load('pyatsaveSext.mat'); atwritem(RING,'matlabwriteSext.m');exit;"
python -c "import at; ring=at.load_m('matlabwriteSext.m')"

I use a script setatversion.m to set te atversion where the environment variable CODEDIR should contain the path to at:

%% setatversion.m
if isempty(getenv('CODEDIR'))
   fprintf('\nEnvironment variable CODEDIR is not set \n');
   fprintf('... leaving matlab execution.\n');
   return;
else
   codeDir = getenv('CODEDIR');
end

addpath(fullfile(codeDir,'atgithub/at/atmat')); % Set AT paths
atpath();

[lfarv]

@oscarxblanco: Ok, I got it. I prepare the fix.

[oscarxblanco]

@swhite2401 There will be a change for the user when calling the function to construct the ID because the function is now part of the class KickMap.

u20 = at.KickMap.InsertionDeviceKickMap('u20',10,"u20_g45mm_kicks_2022nov10.txt",2.75)

Apart from it, this PR should have no impact on the user. The easiest way I found to make visible InsertionDeviceKickMap before instantiating the class was by declaring a static method https://github.com/atcollab/at/blob/eb5bc53da070950db18a609782fd0da4eaf35ded/pyat/at/lattice/idtable_element.py#L39C5-L40C1

I attach a zip file where several combinations of read/write .m and .mat files in pyat and matlab are tested. saveIDtest01.zip

[lfarv]

Hello @oscarblanco I had a look at the python part and have some concern about 2 points:

1. __init__. The call to the constructor of the superclass is not in __init__ but in set_params. This is potentially dangerous, even if not fully relevant here. On one side you can create a KickMap object without calling it (example: KickMap()). In general, this may leave the "base" part of the object in an unpredictable state. On the other side, any call to set_params will call the superclass constructor, which is not desirable. So the call to super().__init__ should be in__init__ itself. I would even push to apply whenever possible the Matlab rule: the call to the superclass constructor must appear only once, outside of any branch. So that you are sure the you cannot miss it. But that's not required!
2. You need to use a "factory function" to build a KickMap object, instead of simply calling its constructor.

I have a proposal which avoids these two problems:

1. You could slightly modify the static InsertionDeviceKickMap method so that instead of returning a KickMap object, it simply returns a dictionary of all the data, excluding family_name. 13 items, from PassMethod to ytable. It may be a static method like now, or a normal private function. Let me call it from_text_file(Nslice, Filename, Energy), choose the name you like.

2. __init__: as it is now, it expects exactly 14 arguments(family_name + 13 others). You can check this to define 2 alternative signatures for the constructor:

   def __init__(self, family_name: str, *args, **kwargs):
      if len(args) < 13:
          # Get data from text file
          elemargs = from_text_data(*args, **kwargs)
      else:
          # Get data from arguments
          elemargs = dict(zip(_argnames, args))
          elemargs.update(kwargs)
      super().__init__(family_name, **elemargs)
   

   With:

   _argnames = ['PassMethod', 'filename_in',..., 'ytable']

And that's all. set_params disappears, you can check/convert the arguments types using the Element's _conversions dictionary:

def _anyarray(value):
    # Ensure proper ordering(F) and alignment(A) for "C" access in integrators
    return numpy.require(value, requirements=['F', 'A'])

class KickMap(Element):
    _BUILD_ATTRIBUTES  = ...
    _conversions = dict(Element._conversions, Nslice=int, xkick=_anyarray, ...)
    ...

[oscarxblanco]

@lfarv Thank you for the suggestion. I'm implementing the modifications and returning to the InsertionDeviceKickMap name for the class.

[oscarxblanco]

Dear @lfarv I have implemented the modifications you suggested, and it looks much better.

I attach a zip file with the scripts I used for testing. saveIDtest02.zip

NOTE: This pull request addresses the limitation when saving files since #558 . However, you mentioned additional blocks https://github.com/atcollab/at/pull/558#issuecomment-1420829127

The InsertionDeviceKickMap still does not implement the additional blocks xkick1 and ykick1, and assume they are filled with zeros. I don't have an example of an insertion device that uses these two parameters.

[oscarxblanco]

@swhite2401 do you have any suggestions ?

[swhite2401]

Hi @oscarxblanco , sorry I am just back from vacations. No objection on my side please go ahead with the merge!

[oscarxblanco]

@swhite2401 No problem