Add support for time-dependent operators

[goerz]

QNET allows to drive the numerical simulation of a symbolic master equation through several backends (qutip and QSD at the moment). We may wish to allow for time-dependent operators in such a numerical simulation. For example, we consider the Hamiltonian for a driven two-level system, H = H_0 + eps(t) * H_1, where H_0 = -(Delta/2)*SigmaZ, H_1 = SigmaX, and eps(t) is a time-dependent pulse, e.g. a Gaussian. In general, eps(t) may be a real or a complex function.

It is important to note that for the derivation of the master equation (QNET's main task), eps(t) is purely symbolic, and can be treated the same way as a numerical constant. Only at the level of the numerical simulation, eps(t) needs to be associated with its actual time-dependent definition. There are several possibilities to address time-dependent operators. For example,

1. Support for time-dependent operators could be pushed into the QNET modules that serve as drivers for the numerical backends. At the level of the QNET data structures, there would be no support for time-dependent pulses. In any case, when setting up the numerical simulation, the user has to map any symbolic constant to a specific numerical value. Support for time-dependent pulses would be added at this level: instead of mapping a symbol to a numerical value, it would be mapped to a time-dependent function. This would be backend-specific. For QSD, it would have to be mapped to a C++ snippet implementing the desired function. For qutip, the symbol would be mapped to a Python function eps(t). The specific time-dependence (i.e, that Gaussian shape, in our example) would be decided by the user when setting up the simulation.
2. QNET could recognize the occurrence of a time parameter (t) in the names of symbols. In the driver modules for the numerical backends, this would be taken as an indication that the symbol has to be mapped to a function instead of a numerical value. The workflow would otherwise follow the same procedure as above (i.e., the specific shape of the pulse would be decided in a backend-specific way by the user when setting up the simulation)
3. For either of the above two cases, as an alternative to leaving the pulse shape definition up to the user when setting up the simulation, a symbolic formula for the time dependence associated with eps(t) could be stored in the QNET data structures (but separate from the symbol eps(t)). Besides having to decide where in the data structure this extra information should be stored, this would then pose the problem that we need to translate the formula into backend-specific code (a Python function for qutip and a C++ code snippet defining a function for QSD), a task that does not appear completely trivial, although doable.
4. We could recognize only t as a symbol representing time, and fully write out eps(t) at the level of the QNET symbolic expression (i.e., instead of eps(t), the formula for the Gaussian would appear). In the backend, the only time-dependent function required would be f(t) = t in that case, which is easy to hard-code. However, the backend would then have to algebraically re-evalate the function at every application, which adds significant overhead and would unnecessarily slow down the simulation.

[goerz]

Regarding option (3) of storing the formula for the time dependence in the QNET data structures: One option would be to give the symbols additional properties, e.g. a boolean flag time_dependent, a (sympy) symbol time_symbol, and a sympy expression time_dependence. Arguably, the flag time_dependent is superfluous, as time-dependence is also indicated by the presence of the other two properties.

Currently, most (all?) symbols appearing in a QNET expressions are instances of sympy.Symbol. Obviously, this means we cannot add new properties directly, but we could add a new type TimeDependentSymbol to QNET that is a subclass of sympy.Symbol and adds the above properties. This seems like it could a a very elegant solution to the problem of how to store the time-dependence in a backend-independent way. We still leave open the question on how to translate the data in the time_dependence property into an appropriate backend-specific format.

[goerz]

Assuming we take the route proposed my previous comment, for the qutip backend, Sympy's lambdify routine should do the trick of converting a sympy expression to a Python function eps(t). For QSD, the situation is more complicated, but fundamentally it is not any different from the QNET-symbolic-expression-to-C++-code-conversion that we're already doing when converting the operators. It would just add another layer.

[goerz]

Added branch feature/timedep to address this.

[goerz]

After some reflection, my proposal for implementation would be the following:

Add TimeDependentSymbol class

1. In qnet.algebra.abstract_algebra module, define a new class TimeDependentSymbol that is subclass of sympy.Symbol, and behaves as follows:

   • Instantiate as e.g.

     t, t0, E0, sigma = sympy.symbols('t, t_0, E_0, sigma', real=True)
     gaussian = E0 * sympy.exp(-(t-t0)**2/sigma**2)
     eps = TimeDependentSymbol(r'\eps(t)', real=True, time_symbol=t, time_dependence=gaussian)

     The time_symbol and time_dependence arguments are mandatory.

   • has properties time_symbol and time_dependence that return the parameters given during instantiation

Account for time dependence in QSD backend

1. In the C++ program template of the QSDCodeGen class, move the {PARAMETERS} section to outside of the main() function, so that all numerical parameters are global variables. This allows to reference them both in the main function, and in any function that implements a time dependence, see below.

2. In QSDCodeGen.__init__, the syms attribute is set to all the free symbols in the Hamiltonian. Note that this will include instances of TimeDependentSymbol, as they are a subclass of Symbol.

   Augment the syms with all the free symbols used in the time dependence formulas:

   • Iterate over existing symbols in syms
   • If symbol has attribute time_dependence: add all free symbols in time_dependence to syms, excluding the time_symbol,.

3. Implement a new method QSDCodeGen._function_lines that returns a list of lines of C++ functions that define any time dependencies.

   • Iterate over any symbol of QSDCodeGen.syms.

   • If symbols has attribute time_dependence and time_symbol, generate a C++ function:

     {type} {symbol}_function(double {time_symbol}
     {
      {formula}
     }

     where {type} is double if symbol.is_real or Complex otherwise, and {formula} is the result of passing the time_dependence attribute to the QSDCodeGen._scalar_str method

4. Add {FUNCTIONS} to the QSDCodeGen program template, after {PARAMETERS}, but before the start of the main function. Insert the lines returned by _function_lines method for this placeholder.

5. modify the QSDCodeGen._parameters_lines method so that any symbol that has a time_dependence attribute is not assigned a value from the num_vals dictionary, but instead is assigned the function defined above, either as

   RealFunction {symbol} = {symbol}_function;

   or

   ComplexFunction {symbol} = {symbol}_function;

   depending on symbol.is_real.

Account for time dependence in qutip backend

The implementation of time dependence for the qutip backend uses the fact that qutip only allows for a simple (quasi-)linear time dependence of the form

H = H_0 + f_1(t) * H_1 + f_2(t) * H_2 + ...

In terms of the QNET operators, this means we only have to worry about a time dependence of a Hamiltonian that is of type OperatorPlus, and check the coefficients of any of its summands of type ScalarTimesOperator for time dependence.

1. Add a new method split_time_dependence to the class qnet.algebra.operator_algebra.OperatorPlus, which returns a nested list of the same structure as used by qutip . For example,

      H = H_0  =>  [H_0, ]
      H = H_0 + f_1(t) * H_1 => [H_0, [f_1(t), H_1]]
      H = H_0 + f_1(t) * H_1 + f_2(t) * H_2 => [H_0, [f_1(t), H_1], [f_2(t), H_2]]
      H = H_0 + f_1(t) * H_1 + f_2(t) * H_2 + c * H_c => [H_0 + c * H_c, [f_1(t), H_1], [f_2(t), H_2]]

   Note that f_1(t) etc. are still instances of TimeDependentSymbol.

2. Modify the qnet.algebra.circuit_algabra.SLH.HL_to_qutip method to return a nested list instead of a single object for H or any of the L (only if there is a time dependence). When processing H, go according to the following steps:

   • If H is not instance of OperatorPlus, convert directly with H.to_qutip
   • Otherwise, run H_split = H.split_time_dependence()
   • If H_split only has a single element (i.e., H really was time-independent), discard H_split and convert H directly with H.to_qutip
   • Otherwise, process each element of H_split (modifying H_split in place):
     • If element is single object, convert it directly with the to_qubip method
     • If element is list [f(t), H] where f(t) is an instance of TimeDependentSymbol and H is an Operator, convert f(t) to a python function via the sympy lambdify routine, and convert H with H.to_qutip

   The same procedure is also used for every L

[goerz]

@ntezak Any comments or suggestions? If not, I'll start implementing this as outlined above.

[ntezak]

Hey, I would prefer not introducing any additional classes specifically for time-dependent symbols as I think it is not necessary. Instead of having eps(t) be a more abstract "gaussian" time dependent coefficient, why not just use symbolic sympy expressions directly such as exp(-(t-t0)**2/2/s**2) where t0 and s are additional params?

[ntezak]

Here is some very quick/n/dirty code I wrote to convert an expression to the qutip time dependence: I've put it into this gist.

[goerz]

Yes, I agree, we can avoid the TimeDependentSymbol, without any drawbacks that I can see. Thinking it through, the proposal for the implementation details would then be as follows:

Account for time dependence in QSD backend

1. In the C++ program template of the QSDCodeGen class, move the {PARAMETERS} section to outside of the main() function, as in previous proposal
2. Add an optional parameter time_symbol as an argument (and property) to QSDCodeGen. If given, time_symbol is excluded from the syms attribute.

3. Implement a new method QSDCodeGen._function_lines that takes a parameter op (instance of qnet.algebra.operator_algebra.Operator, i.e. the Hamiltonian) that returns a list of lines of C++ functions that define any time dependencies.

   • Recursively search for any instances of ScalarTimesOperator in the operands of op, where the coefficient (coeff) includes the time_symbol

   • generate a C++ function:

     {type} tfunc{counter}(double {time_symbol}
     {
      {formula}
     }

     where {type} is double if the coefficient is real (coeff.is_real), or Complex otherwise, and {formula} is the result of passing the coeff attribute to the QSDCodeGen._scalar_str method. The {counter} is an internal counter that numbers different function in situations like H = H_0 + tfunc1 * H_1 + tfunc2 * H_2

   • Choose a time-function-placeholder, defaulting to e.g. u1 for tfunc1, or u1_2, u1_3, ... if u1 already is in syms.

   • Store the mapping

     coeff => (time-function-name, time-function-placeholder, is_real)

     in an attribute _tfuncs (note that the previous two steps can be skipped if coeff is already a key in _tfuncs.

4. Add {FUNCTIONS} to the QSDCodeGen program template, see previous proposal
5. In the _QSDCodeGen.operator_str method, for any term that is an instance of ScalarTimesOperator, and if its coefficient coeff is a key in the _tfuncs attribute, replace the coefficient with the time-function-placeholder

6. modify the QSDCodeGen._parameters_lines method so that any combination of time-function-name and time-function-placeholder in the _tfuncs attribute generates a definition of the form

   RealFunction {time-function-placeholder} = {time-function-name};

   or

   ComplexFunction {time-function-placeholder} = {time-function-name};

Account for time dependence in qutip backend

1. Add function from gist to qnet.algebra.circuit_algebra as _time_dependent_to_qutip, with the following modifications:
   • extra option convert_as='pyfunc' that can take value of either 'pyfunc' or 'str'. For value pyfunc, any coefficient is converted with sympy lambdify routine. For value 'str', the coefficient is converted using str, as currently in the gist (allowing for fast, compiled formula evaluation for a limited set of functions)
   • if operator is not time-dependent, return it directly, instead of [op, 1]
2. Modify the qnet.algebra.circuit_algabra.SLH.HL_to_qutip method to take an optional parameter time_symbol. If given, use the _time_dependent_to_qutip routine for the conversion of H and any L.

[goerz]

Kurt Jacobs pointed out that there is another situation that we should keep in mind as a future extension: there may be time dependencies that cannot easily be expressed as a sympy expression. Most importantly, these would be piecewise constant or similar discontinuous functions. For QSD, this would imply giving a custom C++ function for a given symbol, and for qutip, a custom Python function.

I think we could easily extend the above proposal to account for this. For the QSD backend, we could add an optional dict parameter make_time_dependent that maps symbols to C++ code snippets. Then,

• In QSDCodeGen._function_lines, the code snippet would be inserted, and the function name extracted
• In QSDCodeGen._parameters_lines, for symbols that are keys in_make_time_dependentdict, don't look for a numeric value fromnum_vals`, but connect to the function name

For the qutip backend,

• In the SLH.HL_to_qubit method, also add dict parameter make_time_dependent that maps symbols to Python functions
• Pass the make_time_dependent dict to _time_dependent_to_qutip, and treat is analogously to time-dependent coefficients.

I don't think we have to implement this at the moment, but can keep it as a plan for the future, if needed.

[ntezak]

Yeah, that sounds good! I would guess that such "input vectors" are most often used for piecewise constant expressions (at least that is how I have used them), i.e., you let some parameter value be determined by what time bin it falls into. If we allowed for a primitive function/expression type like that then we might still be able to combine this in expressions with other types of coefficient expressions/functions.

[goerz]

It seems that sympy does support Piecewise functions, and can convert them automatically into C or Fortran code (using conditional), see http://docs.sympy.org/latest/modules/printing.html#sympy.printing.ccode.ccode I haven't tested it yet, but it should work pretty much out of the box for both the qutip backend and the QSD backend