using specific guess pulse for grape/crab

[johny maartin]

Hi, is there a way to set the initial (guess) control field to be specific fields? It seems like currently optimize_pulse function is only able to deal with few different types of such as RND, SINE, SQUARE, TRIANGLE, SAW, etc. I'm wondering if we can start from some specific functions of pulse other than these. 

Thanks.

[Simon Cross]

Hi Johny,

I *think* that for CRAB this is possible by passing alg="CRAB" and a
list of PulseGenCrabFourier instances as the pulse_generator.

I recently reviewed a PR for adding this to the GRAPE solver too, but
it needs more work before it can be merged --
https://github.com/qutip/qutip/pull/1369.

A new PR would be very welcome, but hopefully this earlier work will
be helpful to you regardless.

Regards,
Simon

[johny maartin]

Hi, thank you for your reply. 

It seems like you're right that we can use PulseGenCrabFourier instances to specify initial parameters, although I haven't tried it yet. Do you know any good example notebook that uses PulseGen to specify initial pulse? Also, do you know how to generate parameters after the optimization is done? 

Regarding GRAPE, so I assume that currently we can't customize the initial control fields? If so, I look forward to hearing the updates in the near future!

Thanks.

[Simon Cross]

Hi Johny,

For GRAPE you can write your own PulseGen class and use that to supply
a custom initial pulse, but there can only be one initial pulse object
for all controls.

You may have seen them already, but under the "Optimal control"
heading at https://qutip.org/tutorials.html there are a set of example
notebooks for the control optimization in QuTiP. The "QFT (CRAB)"
example demonstrates some simple setting of the initial pulse using
pulsegen:

https://nbviewer.ipython.org/github/qutip/qutip-notebooks/blob/master/examples/control-pulseoptim-CRAB-QFT.ipynb

What the example doesn't show clearly is that guess_pulse can be a
list of pulsegen instances with one item per control.

Regards,
Simon

[johny maartin]

Hi, thanks for the reply. 

I see, so for GRAPE I can customize the initial pulse, but that should be the same initial fields for all fields. I want to make sure because I'm planning to begin with different sets of specific initial pulses (i.e. different u1(0), u2(0), ... etc) and I'd like to make sure that GRAPE is not supported for such task at least for now.

I have seen lots of other notebooks, but first time seeing what you mentioned. I'll check it out. Thank you for your help.

[Simon Cross]

> I see, so for GRAPE I can customize the initial pulse, but that should be the same initial fields for all fields. I want to make sure because I'm planning to begin with different sets of specific initial pulses (i.e. different u1(0), u2(0), ... etc) and I'd like to make sure that GRAPE is not supported for such task at least for now.

Correct. The PR https://github.com/qutip/qutip/pull/1369 was adding
exactly the support you asked for, so you could use that branch if you
wanted to.

[johny maartin]

Hi Simon, I have a follow-up question regarding CRAB function or its example notebook you mentioned: how can I print out the final (optimized) parameters? As far as I know, CRAB optimizes cost function wrt parameters of functions, e.g. if function is A*sin(wt) + B*cos(wt), then one choice would be optimizing for A, B, and I'm wondering how can I print out such parameters after the training is done. 

I have another question. I want to maximize Tr(\rho(T)*O), where O is some Hermitian operator and rho(T) is fully time-evolved state with given Hamiltonian. However, I don't see how to implement the CRAB algorithm for density matrix. What I tried is as follows: from (https://nbviewer.ipython.org/github/qutip/qutip-notebooks/blob/master/examples/control-pulseoptim-CRAB-QFT.ipynb), I changed U_0 = |0><0|, U_targ = O, and fid_type = 'TRACEDIFF', which I believe should maximize Tr(\rho(T)*O), right? But I keep getting a very non-physical density matrix that is not a quantum state. I ran the code with varying different parameters such as setting dyn_type = 'UNIT', replacing Lindbladian instead of Hamiltonian, and etc, but nothing actually works. 

I'd appreciate it if you can help me out with the two issues I'm having with. 

Thank you.

[Simon Cross]

Hi Johny,

Regarding the first question about the CRAB optimization variables, it
looks like each CRAB pulse is described by an instance of
pulsegen.PulseGenCrabFourier and these are available after the
optimization as optim.pulse_generator (as seen in the notebook). The
coefficients of the fourier modes are accessible via `.coeff`. The
frequency of each mode is 2 * pi * N / pulse_time for N = 1 .. number
of coefficients.

Regarding the unphysical density matrix -- I guess this has nothing to
do with the pulse optimization, but rather to do with whether the
evolution of the system is unitary, so perhaps examining the system
evolution is the right way to debug it? If you can post a small
example, perhaps someone can spot what is wrong. Of course it is also
possible that there is a bug, but the first step is to get a failing
example and then for you and maybe us to figure out what is happening.

Regards,
Simon

[johny maartin]

Hi, Simon. Much appreciated for your response. I figured out the problem with unphysical density matrix part a few days ago.

Regarding the coefficients of CRAB pulses, however, I'm still struggling to understand what are exactly these ".coeffs" related to the guess pulse. Basically, I produced a simple sine pulse as a (fixed) guess pulse for optim.pulse_generator[0] as follows:

crab_pgen = optim.pulse_generator[0]
guess_pgen = pulsegen.PulseGenSine(dyn=dyn)
guess_pgen.freq = 10/(evo_time*2*np.pi)
guess_pgen.start_phase = 0
guess_pgen.scaling = 1.0
guess_pgen.offset = 0
crab_pgen.guess_pulse = guess_pgen.gen_pulse()
crab_pgen.scaling = 1.0
crab_pgen.offset = 0.0

But whenever I print out "optim.pulse_generator[0].coeffs", it gives me different sets of parameters, which means that these parameters are not related to the guess sine pulse at all. Then what exactly do these parameters mean related to the guess pulse? 

The main reason I'm confused is because as soon as we set the initial parameters that the optimization begins with (i.e. '.coffes' and 'freqs'), I thought the optimization starts with initial pulses generated by that parameters. However, it seems like that's not the case. As said above, it seems like these parameters are independent with guess pulses, which is very confusing to me. 

Thanks.

[Simon Cross]

Hi Johny,

It looks like the guess pulse is used in an unexpected way by QuTiP's
CRAB algorithm.

If `guess_pulse_action` is not set it is entirely ignored (perhaps
this is the case you are encountering currently).

The other allowed values for `guess_pulse_action` are "ADD" and
"MODULATE". "ADD" adds the guess pulse to every generated CRAB pulse.
"MODULATE" multiplies every generated pulse by the guess pulse. This
is rather unexpected. Instead I would have expected that the CRAB
basis coefficients were initialized to the best fit to the guess
pulse.

The options that affect the default CRAB basis coefficients are
"randomize_coeffs" (which sets each coefficient to a random value in
[-scaling, scaling]). By default each coefficient is set to the
scaling.

I think this explains what is happening, but it might not be that
helpful. You could also just try setting `.coeffs` directly to what
you want them to be after creating the crab_pgen object and before
using it in the optimizer.

Regards,
Simon

[johny maartin]

Thanks for your comment. 

As you said, specifying " coeffs " would be the most straightforward way to set the initial pulse. But, I'm actually facing a very minor problem. How can I get rid of the effect of "ramping pulse"? I've set all my " coeffs " and " freqs " correctly (for CRABFourier, i.e. https://qutip.org/docs/latest/modules/qutip/control/pulsegen.html#PulseGenCrabFourier.gen_pulse), but apparently the effects coming from "ramping pulse" somehow mess my pulse shape up. I tried with setting it to False, 0, 1, ..., etc, but all of them didn't get rid of the effect. 

p.s.) I think simply setting " dyn.initialize_controls(init_amps) ", where init_amps is which control amplitudes I want to start with, might be the solution of initial pulse problem.  For instance, related to the last code I posted, the following additional code might solve the problem and thus the optimization starts with init_amps:

dyn = optim.dynamics

crab_pgen = optim.pulse_generator[0]
guess_pgen = pulsegen.PulseGenSine(dyn=dyn)
guess_pgen.freq = 10/(evo_time*2*np.pi)
guess_pgen.start_phase = 0
guess_pgen.scaling = 1.0
guess_pgen.offset = 0
crab_pgen.guess_pulse = guess_pgen.gen_pulse()
crab_pgen.scaling = 1.0
crab_pgen.offset = 0.0

init_amps = np.zeros([n_ts, n_ctrls])
for j in range(dyn.num_ctrls):
    pgen = optim.pulse_generator[j]
    pgen.init_pulse()
    init_amps[:, j] = pgen.guess_pulse
   
dyn.initialize_controls(init_amps)

But I'm not sure how to actually check whether the optimization done via the code above actually began with the "init_amps" or not. Do you think this solves the problem and/or any way to check it? 

[Simon Cross]

Hi Johny,

It appears from the code that ramping pulse defaults to not being used, so I'm not sure what you are seeing happen.

Regarding the initial amplitudes, I'm not sure whether that is used by the CRAB optimizer -- the CRAB pulses are defined by their coeffs.

Do you perhaps have a complete example on somewhere like GitHub that you could link to? Then I or someone else could look at it / try run it and see what is happening.

Regards,

Simon

[johny maartin]

Hi Simon, 

I put lots of time to figure this out. This was why: 

If we set bounds for "create_pulse_optimizer" function, e.g. amp_lbound=[-1.0,-1.0,-1.0,-1.0], amp_ubound=[1.0,1.0,1.0,1.0], this automatically sets "_bound_scale_cond" to "_BSC_ALL", which apply "tanh" function to the pulse. I think the problem of this function is that it shrinks the input pulse even though the pulse doesn't go over/below 1/ -1. So for instance, if I want to initiate my guess control as cosine function so that every value is within 1 and -1, and try to set bounds for amps during my training (so that the optimized control won't go over/below 1/-1), I couldn't properly generate the cosine function as my guess control via "gen_pulse()" since "apply_bound" function automatically shrinks the pulse. This has been the major problem to what I want to do with CRAB package, and I'm not quite sure how to resolve this issue.

You're right that it's hard to see the problems without looking at the codes directly. I can make a few examples in the github to address such issues. 

Thank you.

[Simon Cross]

Hi Johny,

Thank you for digging into this and sticking with it. If you have an example which you think demonstrates a bug, you can include it in an issue report at https://github.com/qutip/qutip/issues/new/choose and we can try fix it together.

Regards,

Simon

[johny maartin]

Sure thing. I'll post them after some more examinations. 

Just one more minor question: is it ever possible to customize the cost function (e.g. fidelity) for CRAB? 

[Simon Cross]

> Just one more minor question: is it ever possible to customize the cost function (e.g. fidelity) for CRAB?

There isn't a built-in way. One would need to write a new OptimizerCrab sub-class.

[Leeseok Kim]

Hi all, 

I personally think that customizing initial pulse for GRAPE is possible, just by making np.array and initialize that array as the initial pulse by using dyn.initialize_controls. 

I uploaded an issue, which is not quite relevant to this but includes a way to customize the initial pulse, in here: https://github.com/qutip/qutip/issues/1900. Feel free to check it out and please let me know if it's wrong. But I personally think this does the job of initializing the pulse well.

Best,

Leeseok

[Johnny Lin]

I think setting initial pulse for GRAPE is different from setting initial pulse for CRAB, since for CRAB one has to set the Fourier Series coefficients of the intended initial pulse. However, I don't think it's built in and I had to modify the source code. If someone knows a better way would really appreciate to hear about it!