-
Notifications
You must be signed in to change notification settings - Fork 59
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
Matching/optimization of arbitrary matrix elements and their propagation #147
Comments
I suppose Sergey knows more about this than I do, but I will try to express my thoughts.
I guess you can write your own yes, but if you look in cell 14 of that notebook, you can for your lattice call this function repeatedly and try to minimise arbitrary elements of this matrix using whatever optimisation routine you would like. Actually this function is quite odd because if you look at its definition at Lines 75 to 94 in fc4938a
MagneticLattice instance. So this is how you can get the Ts from it. But if you just do lattice.transfer_maps(energy) then it neatly returns B, R and T. I would prefer the latter solution because its tidier.
Yes I agree. It is not difficult to implement something like this though. Look in the above function to see how to loop the lattice and get the maps and yield the individual maps rather than combining them all into a single net one. If it doesn't already exist then perhaps you can write this and a test and send a pull request.
The solution to this is to attach |
Hi, Thanks for your thoughts, @st-walker ! Sorry for the late reply, it took me some time to manage to figure things out. I tried the suggested route of implementing my own routine using (small note - the "cleaner" format you suggested for extracting the R and T matrices was a quite new addition, just had to update my local files for it to run. Might I ask what's in the B-matrix?) As for the matrix element evolution, I hope to get to that eventually... It seems a bit inefficient to me to loop through the lattice again just to get this information, and it might be more efficient to do it as part of the current/usual matrix multiplication - at least if the option to do is called. I don't really know much about proper Python programming or actually using Github, but I can give it a shot. Worst case, I guess I could send one of you developers my code? Similar opinion on using Have a nice weekend!
|
Sorry for late reply, I was at a conference and then on holiday. I'm not sure what the B matrix is actually. If you want to send any code then that is fine. I am not sure what you have in mind regarding the matrix multiplication or an alternative. OK I see what you mean, dumping the whole beam is too much, maybe you want just a single aggregate property of the beam. Yes you can do that, just define your own physics procwss by inheriting
and then add the physics process to one or more Marker instances just like with Copy/SaveBeam |
Hi, twiss parameters which is returned from track function has many parameters like first and second moments. I just checked and it seems it does not have only tws.pp :) I will implement it in next commit. Cheers, |
Hi @st-walker , So, I finally got back into this topic. Regarding the optimization, I ended up writing my own, generalized optimizer based on the beam-based optimization in the example notebook. It's not perfect but it can take any input I want. On the topic of calculating the propagation of the matrix elements, I looked into
It's definitely not a "final" implementation, but it does almost what I want. I'd imagine that one could keep only the last loop (the one I'm trying to add) and then just output for example The Can you tell what might be going wrong here? |
Although it might not be the cause of your problem this is definitely wrong:
If it's not what's causing your bug (though I think it probably is), then you will have to send me something to run, because it's difficult to debug purely by eye. |
Hi, The suggested fix produces something which is completely different whereas it was almost correct before, up until the dipole in the lattice. I'll send you a run file by email!
|
OK yes my mistake, I meant something like this:
but my brain failed and I introduced the enumerate into the inner loop for some reason. That will indeed be a bug for arbitrary lattices, but in your case it doesn't make a difference because you have no duplicated elements. Anyway it's tricky to debug, my advice would be to try to create the smallest/simplest lattice that recreates this problem and then try to debug, and I can try to help you. |
Hi, I tried making a very short lattice, it still seems something weird happens at the dipole but before then it's all ok with the code I posted above. The code I use is pasted below.
For this code, the "normal" calculation of R11 at each lattice position and my implementation yield, respectively
Turning dipole edge angles and fint off did not solve it, but it's clearly related to the dipole. The third entry for the dipole is in my case the same as the first one, rather than the previous one. Hope this help clarify the problem. |
Hi @st-walker, Did you have time to check the code above? Best regards, |
Yes I had a look at it but it was not immediately clear to me why it is happening. I was stepping through the code line by line with a debugger, I spent maybe about half an hour on this. Did you try stepping through it with a debugger? I think the example is still too complicated. Does the bug exist with one drift and one dipole for example? If you can make the bug go away then this will help, if not, then at least your test case is simpler so that when you step through it with a debugger then it will be much easier for you to fix it. |
Hi, My Spyder refuses to recognize the breakpoints, so I haven't managed to step through it the way I would like. I can maybe try again, or with a different method, once I get more time to return to these simulations. Regardless, the bug looks exactly the same for the above example as for my originally more complicated setup, i.e. the discrepancy starts inside the dipole.
|
Hi Jonas, Stuart, all
I was looking at the code you were discussing in this thread.
First of all, I think that it would make sense to extend the function
"lattice_transfer_map()" with an additional kwarg and not to implement a
new method in the "MagneticLattice" class.
I would propose the following, which I think is quite simple, since it
requires few lines of code:
def lattice_transfer_map(lattice, energy, output_at_each_step: bool =
False):
"""
Function calculates transfer maps, the first and second orders (R,
T), for the whole lattice.
Second order matrices are attached to lattice object:
lattice.T_sym - symmetric second order matrix
lattice.T - second order matrix
lattice.R - linear R matrix
:param lattice: MagneticLattice
:param energy: the initial electron beam energy [GeV]
:param output_at_each_step: option to calculate the transfer matrix
elements along the beamline
:return: R - matrix
"""
Ba = np.zeros((6, 1))
Ra = np.eye(6)
Ta = np.zeros((6, 6, 6))
Bs, Rs, Ts = [Ba], [Ra], [Ta]
E = energy
for elem in lattice.sequence:
for Rb, Bb, Tb, tm in zip(elem.R(E), elem.B(E), elem.T(E),
elem.tms):
Ba, Ra, Ta = transfer_maps_mult(Ba, Ra, Ta, Bb, Rb, Tb)
#Ba = np.dot(Rb, Ba) + Bb
E += tm.get_delta_e()
Bs.append(Ba)
Rs.append(Ra)
Ts.append(Ta)
# TODO: Adding Attributes at runtime should be avoided
lattice.E = E
lattice.T_sym = Ta
lattice.T = Ta #unsym_matrix(deepcopy(Ta))
lattice.R = Ra
lattice.B = Ba
if output_at_each_step:
return Bs, Rs, Ts
return Ra
It actually seems to work properly if I compare it with the dispersion
calculated with the "twiss()" function and also seems to agree with the
"normal" results of the test lattice that Jonas proposed.
It's a bit weird that the kind of output changes so much if
"output_at_each_step" is selected (i.e. three lists instead of only the
total Ra matrix). But maybe this is better than requiring additional
kwargs to select the function's output...
What do you think?
Cheers,
Pau
…On 10/26/22 20:08, jonasbjorklundsvensson wrote:
Hi,
My Spyder refuses to recognize the breakpoints, so I haven't managed
to step through it the way I would like. I can maybe try again, or
with a different method, once I get more time to return to these
simulations.
Regardless, the bug looks exactly the same for the above example as
for my originally more complicated setup, i.e. the discrepancy starts
inside the dipole.
* Jonas
—
Reply to this email directly, view it on GitHub
<#147 (comment)>,
or unsubscribe
<https://github.com/notifications/unsubscribe-auth/AIC5DCAOC4IYWBBZL66JLLLWFFXSHANCNFSM56ZBXJFA>.
You are receiving this because you are subscribed to this
thread.Message ID: ***@***.***>
--
Pau Gonzalez Caminal
FLASHForward - Deutsches Elektronen Synchrotron (DESY)
Notkestr. 85, 22607 Hamburg
Bldg. 1E / Room 02.514
Phone: +49 40 8998 4960
***@***.***
|
Hi again,
just one more comment. @stuart: in your first answer to Jonas, you sent
a link to a particular line of code. Might it be that the commit you
were referring to was detached from the "official" branches? If I look
at master or dev_2021 I don't see that the method
MagneticLattice.transfer_maps() exists...
Am I wrong? Or is there a pull request on its way that I am not aware of?
Cheers,
Pau
…On 11/8/22 14:40, Pau Gonzalez Caminal wrote:
Hi Jonas, Stuart, all
I was looking at the code you were discussing in this thread.
First of all, I think that it would make sense to extend the function
"lattice_transfer_map()" with an additional kwarg and not to implement
a new method in the "MagneticLattice" class.
I would propose the following, which I think is quite simple, since it
requires few lines of code:
def lattice_transfer_map(lattice, energy, output_at_each_step: bool =
False):
"""
Function calculates transfer maps, the first and second orders (R,
T), for the whole lattice.
Second order matrices are attached to lattice object:
lattice.T_sym - symmetric second order matrix
lattice.T - second order matrix
lattice.R - linear R matrix
:param lattice: MagneticLattice
:param energy: the initial electron beam energy [GeV]
:param output_at_each_step: option to calculate the transfer
matrix elements along the beamline
:return: R - matrix
"""
Ba = np.zeros((6, 1))
Ra = np.eye(6)
Ta = np.zeros((6, 6, 6))
Bs, Rs, Ts = [Ba], [Ra], [Ta]
E = energy
for elem in lattice.sequence:
for Rb, Bb, Tb, tm in zip(elem.R(E), elem.B(E), elem.T(E),
elem.tms):
Ba, Ra, Ta = transfer_maps_mult(Ba, Ra, Ta, Bb, Rb, Tb)
#Ba = np.dot(Rb, Ba) + Bb
E += tm.get_delta_e()
Bs.append(Ba)
Rs.append(Ra)
Ts.append(Ta)
# TODO: Adding Attributes at runtime should be avoided
lattice.E = E
lattice.T_sym = Ta
lattice.T = Ta #unsym_matrix(deepcopy(Ta))
lattice.R = Ra
lattice.B = Ba
if output_at_each_step:
return Bs, Rs, Ts
return Ra
It actually seems to work properly if I compare it with the dispersion
calculated with the "twiss()" function and also seems to agree with
the "normal" results of the test lattice that Jonas proposed.
It's a bit weird that the kind of output changes so much if
"output_at_each_step" is selected (i.e. three lists instead of only
the total Ra matrix). But maybe this is better than requiring
additional kwargs to select the function's output...
What do you think?
Cheers,
Pau
On 10/26/22 20:08, jonasbjorklundsvensson wrote:
>
> Hi,
>
> My Spyder refuses to recognize the breakpoints, so I haven't managed
> to step through it the way I would like. I can maybe try again, or
> with a different method, once I get more time to return to these
> simulations.
>
> Regardless, the bug looks exactly the same for the above example as
> for my originally more complicated setup, i.e. the discrepancy starts
> inside the dipole.
>
> * Jonas
>
> —
> Reply to this email directly, view it on GitHub
> <#147 (comment)>,
> or unsubscribe
> <https://github.com/notifications/unsubscribe-auth/AIC5DCAOC4IYWBBZL66JLLLWFFXSHANCNFSM56ZBXJFA>.
> You are receiving this because you are subscribed to this
> thread.Message ID:
> ***@***.***>
>
--
Pau Gonzalez Caminal
FLASHForward - Deutsches Elektronen Synchrotron (DESY)
Notkestr. 85, 22607 Hamburg
Bldg. 1E / Room 02.514
Phone: +49 40 8998 4960
***@***.***
--
Pau Gonzalez Caminal
FLASHForward - Deutsches Elektronen Synchrotron (DESY)
Notkestr. 85, 22607 Hamburg
Bldg. 1E / Room 02.514
Phone: +49 40 8998 4960
***@***.***
|
Hi Pau, it is in the "dev" branch. "dev_2021" is merged with "dev" and not in development any more. I will remove it soon. Cheers, |
Hi Sergey,
Ups, I see. Thanks for the clarification! I didn't check the dates on
the history and, for some reason, I mistakenly thought that "dev_2021"
was the most up-to-date branch (despite of its name...).
Cheers,
Pau
…On 11/8/22 15:44, sergey-tomin wrote:
Hi Pau,
it is in the "dev" branch. "dev_2021" is merged with "dev" and not in
development any more. I will remove it soon.
"dev" and "master" will be also merged in about 1-2 months.
As for your suggestion. It's really easy to implement, and it hasn't
been implemented because we never needed it. If it can be useful to
anyone, I will do it in the next few days.
Cheers,
Sergey.
—
Reply to this email directly, view it on GitHub
<#147 (comment)>,
or unsubscribe
<https://github.com/notifications/unsubscribe-auth/AIC5DCDJQFGM76IHKODF2KDWHJRMHANCNFSM56ZBXJFA>.
You are receiving this because you commented.Message ID:
***@***.***>
--
Pau Gonzalez Caminal
FLASHForward - Deutsches Elektronen Synchrotron (DESY)
Notkestr. 85, 22607 Hamburg
Bldg. 1E / Room 02.514
Phone: +49 40 8998 4960
***@***.***
|
Added output_at_each_step to MagneticLattice.transfer_maps(). see in dev branch a9edf7e |
Hi guys,
I needed to do something similar to what Jonas was asking for in the
beginning of this thread and I worked out a possible way to do it by
enabling the inclusion of matrix elements as a constraint in the "match"
function. The way to define the constraints would be, for instance:
constraints = {'global': {'beta_x': ["<", 300], 'beta_y': ["<", 300]},
{obj: {'R11': -10, 'R12': 0}}
I tested it with some use cases for our beamline and it works fine for me.
One of the problems is that the calculation of the transfer maps
requires to pass the energy of the beam, so I included it as an extra
argument in the "match" function. What do you think?
Here the link to the code:
https://github.com/PauGC/ocelot/blob/0ba27db1e106d219b0d949afe6038510555a7a0f/ocelot/cpbd/match.py#L153
Best,
Pau
…On 11/9/22 12:55, sergey-tomin wrote:
Added output_at_each_step to MagneticLattice.transfer_maps(). see in
dev branch a9edf7e
<a9edf7e>
—
Reply to this email directly, view it on GitHub
<#147 (comment)>,
or unsubscribe
<https://github.com/notifications/unsubscribe-auth/AIC5DCEYPLR7GE7OSATZO73WHOGM7ANCNFSM56ZBXJFA>.
You are receiving this because you commented.Message ID:
***@***.***>
--
Pau Gonzalez Caminal
FLASHForward - Deutsches Elektronen Synchrotron (DESY)
Notkestr. 85, 22607 Hamburg
Bldg. 1E / Room 02.514
Phone: +49 40 8998 4960
***@***.***
|
Sorry, I forgot to push after adding the "energy" as an argument... Here
the new link:
https://github.com/PauGC/ocelot/blob/f6dc2e9334369475ee18babaf26bb9cc6c46dac7/ocelot/cpbd/match.py#L153
…On 11/11/22 13:04, Pau Gonzalez Caminal wrote:
Hi guys,
I needed to do something similar to what Jonas was asking for in the
beginning of this thread and I worked out a possible way to do it by
enabling the inclusion of matrix elements as a constraint in the
"match" function. The way to define the constraints would be, for
instance:
constraints = {'global': {'beta_x': ["<", 300], 'beta_y': ["<", 300]},
{obj: {'R11': -10, 'R12': 0}}
I tested it with some use cases for our beamline and it works fine for me.
One of the problems is that the calculation of the transfer maps
requires to pass the energy of the beam, so I included it as an extra
argument in the "match" function. What do you think?
Here the link to the code:
https://github.com/PauGC/ocelot/blob/0ba27db1e106d219b0d949afe6038510555a7a0f/ocelot/cpbd/match.py#L153
Best,
Pau
On 11/9/22 12:55, sergey-tomin wrote:
>
> Added output_at_each_step to MagneticLattice.transfer_maps(). see in
> dev branch a9edf7e
> <a9edf7e>
>
> —
> Reply to this email directly, view it on GitHub
> <#147 (comment)>,
> or unsubscribe
> <https://github.com/notifications/unsubscribe-auth/AIC5DCEYPLR7GE7OSATZO73WHOGM7ANCNFSM56ZBXJFA>.
> You are receiving this because you commented.Message ID:
> ***@***.***>
>
--
Pau Gonzalez Caminal
FLASHForward - Deutsches Elektronen Synchrotron (DESY)
Notkestr. 85, 22607 Hamburg
Bldg. 1E / Room 02.514
Phone: +49 40 8998 4960
***@***.***
--
Pau Gonzalez Caminal
FLASHForward - Deutsches Elektronen Synchrotron (DESY)
Notkestr. 85, 22607 Hamburg
Bldg. 1E / Room 02.514
Phone: +49 40 8998 4960
***@***.***
|
Hi,
I'm trying to set up some lattice matching/optimization, and I have some questions (and, I guess, feature requests) after looking at this for a little while. Sorry for mixing topics a bit, but I feel that they are related.
One of the things I want to do is to make sure that a section of my lattice is second-order achromatic, i.e. have T166 = T266 = 0 at the section end. However, it seems to me that only "standard" Twiss parameters, such as betas, alphas, and first-order dispersions, are available for matching. Would I need to write my own optimization routine for this, like in the more general accelerator optimization tutorial (https://nbviewer.org/github/ocelot-collab/ocelot/blob/master/demos/ipython_tutorials/accelerator_optim.ipynb)?
As a follow-up question/request: it would also be great to have matrix elements available for plotting along with the first-order optics. It is quite informative to view the propagation of, for example, T166 and T566 along with the beta functions. Also, the first-order angle-to-space matrix elements R12 and R34 would be useful to have - we are often dealing with quite divergent beams, which can also jitter quite a bit in angle, so limiting the maximum of such parameters along certain lattice sections would also be very useful.
It would also be nice to have more beam moments available along the lattice, for example I can easily plot the mean beam energy along the accelerator but the energy spread is, to my knowledge, not there.
Best regards,
Jonas
The text was updated successfully, but these errors were encountered: