def fast_rabi_terms(Ep, epsilonp, rm, xi, theta, unfolding,
                    matrix_form=False, file_name=None, return_code=False):
    r"""Return a fast function that returns the Rabi frequency terms.

    We test a basic two-level system.

    >>> import numpy as np
    >>> from scipy.constants import physical_constants
    >>> from sympy import Matrix
    >>> from fast.electric_field import electric_field_amplitude_top
    >>> from fast.symbolic import define_laser_variables, polarization_vector

    >>> Ne = 2
    >>> Nl = 1
    >>> a0 = physical_constants["Bohr radius"][0]
    >>> rm = [np.array([[0, 0], [a0, 0]]),
    ...       np.array([[0, 0], [0, 0]]),
    ...       np.array([[0, 0], [0, 0]])]
    >>> xi = np.array([[[0, 1], [1, 0]]])
    >>> theta = phase_transformation(Ne, Nl, rm, xi)

    We define symbolic variables to be used as token arguments.
    >>> Eps = [electric_field_amplitude_top(1e-3, 1e-3, 1, "SI")]
    >>> Ep, omega_laser = define_laser_variables(Nl)
    >>> epsilonps = [polarization_vector(0, 0, 0, 0, 1)]

    An map to unfold the density matrix.
    >>> unfolding = Unfolding(Ne, True, True, True)

    We obtain a function to calculate Rabi frequency terms.
    >>> rabi_terms = fast_rabi_terms(Ep, epsilonps, rm, xi, theta, unfolding)

    Apply this to a density matrix.
    >>> rhos = np.array([[0.6, 3+2j],
    ...                  [3-2j, 0.4]])
    >>> rhosv = unfolding(rhos)
    >>> rhs_rabi = rabi_terms(rhosv, Eps)
    >>> print(rhs_rabi)
    [-55680831.474          0.       2784041.5737]

    """
    if not unfolding.lower_triangular:
        mes = "It is very inefficient to solve using all components of the "
        mes += "density matrix. Better set lower_triangular=True in Unfolding."
        raise NotImplementedError(mes)
    if matrix_form and (not unfolding.real) and (unfolding.lower_triangular):
        mes = "It is not possible to express the equations in matrix form "
        mes += "for complex lower triangular components only."
        raise ValueError(mes)
    Nl = len(Ep)
    Ne = unfolding.Ne
    # We determine which arguments are constants.
    if True:
        try:
            Ep = np.array([complex(Ep[l]) for l in range(Nl)])
            variable_Ep = False
        except:
            variable_Ep = True

        try:
            epsilonp = [np.array([complex(epsilonp[l][i]) for i in range(3)])
                        for l in range(Nl)]
            variable_epsilonp = False
        except:
            variable_epsilonp = True
    # We unpack variables.
    if True:
        Nrho = unfolding.Nrho
        # Mu = unfolding.Mu
        IJ = unfolding.IJ
        normalized = unfolding.normalized
        lower_triangular = unfolding.lower_triangular

        # The conjugate stuff.
        Em = [ii.conjugate() for ii in Ep]
        if variable_epsilonp:
            epsilonm = Vector3D(epsilonp.args[0].conjugate())
        else:
            epsilonm = [ii.conjugate() for ii in epsilonp]
        # We convert rm to a numpy array
        rm = np.array([[[complex(rm[k][i, j])
                       for j in range(Ne)] for i in range(Ne)]
                       for k in range(3)])
        rp = np.array([rm[ii].T.conjugate() for ii in range(3)])
        rm_aux = Vector3D(IndexedBase("rm", shape=(Ne, Ne)))
        rp_aux = Vector3D(IndexedBase("rp", shape=(Ne, Ne)))

        # We define needed matrices.
        rho = define_density_matrix(Ne, explicitly_hermitian=lower_triangular,
                                    normalized=normalized)

        if variable_epsilonp:
            Omega = [[sum([xi[l, i, j] *
                          (Ep[l]*dot(epsilonp[l], rm_aux[i, j]) *
                          delta_greater(i, j) +
                           Em[l]*dot(epsilonm[l], rp_aux[i, j]) *
                          delta_lesser(i, j))
                           for l in range(Nl)])
                     for j in range(Ne)] for i in range(Ne)]
        else:
            Omega = [[sum([xi[l, i, j] *
                          (Ep[l]*dot(epsilonp[l], rm[:, i, j]) +
                           Em[l]*dot(epsilonm[l], rp[:, i, j]))
                           for l in range(Nl)])
                     for j in range(Ne)] for i in range(Ne)]
    # We establish the arguments of the output function.
    if True:
        code = ""
        code += "def rabi_terms("
        if not matrix_form: code += "rho, "
        if variable_Ep: code += "Ep, "
        if variable_epsilonp: code += "epsilonp, "
        if code[-2:] == ", ":
            code = code[:-2]
        code += "):\n"
        code += '    r"""A fast calculation of the Rabi terms."""\n'
    # We initialize the output and auxiliaries.
    if True:
        # We introduce the factor that multiplies all terms.
        if unfolding.real:
            code += "    fact = "+str(e_num/hbar_num)+"\n\n"
        else:
            code += "    fact = "+str(1j*e_num/hbar_num)+"\n\n"
        if variable_epsilonp:
            # We put rm and rp into the code
            np.set_printoptions(threshold=np.nan)
            code += "    rm = np."+rm.__repr__()+"\n\n"
            code += "    rp = np."+rp.__repr__()+"\n\n"
            code += "    def dot(epsilon, rij):\n"
            code += "        return epsilon[0]*rij[0]"
            code += " + epsilon[1]*rij[1]"
            code += " + epsilon[2]*rij[2]\n\n"
            np.set_printoptions(threshold=1000)

        if matrix_form:
            code += "    A = np.zeros(("+str(Nrho)+", "+str(Nrho)
            if not unfolding.real:
                code += "), complex)\n\n"
            else:
                code += "))\n\n"
            if unfolding.normalized:
                code += "    b = np.zeros(("+str(Nrho)
                if not unfolding.real:
                    code += "), complex)\n\n"
                else:
                    code += "))\n\n"
        else:
            code += "    rhs = np.zeros(("+str(Nrho)
            if not unfolding.real:
                code += "), complex)\n\n"
            else:
                code += "))\n\n"

    # We write code for the linear terms.
    rho00_terms = []
    for mu in range(Nrho):
        s, i, j = IJ(mu)
        for l in range(Nl):
            for k in range(Ne):
                if xi[l, i, k] == 1:
                    # There is a I* Omega_l,i,k * rho_k,j term.
                    u = k; v = j
                    args = (0.5*Omega[i][k], rho[k, j], s, i, j, k, u, v,
                            unfolding, matrix_form)
                    term_list = linear_get_coefficients(*args)
                    for term in term_list:
                        code += term_code(*term)

                    # We keep note that there was a term with rho00.
                    if k == 0 and j == 0:
                        rho00_terms += [args]

                if xi[l, k, j] == 1:
                    # There is a -I * Omega_l,k,j * rho_i,k term.
                    u = i; v = k
                    args = (-0.5*Omega[k][j], rho[i, k], s, i, j, k, u, v,
                            unfolding, matrix_form)
                    term_list = linear_get_coefficients(*args)
                    for term in term_list:
                        code += term_code(*term)
                    # We keep note that there was a term with rho00.
                    if i == 0 and k == 0:
                        rho00_terms += [args]

    # We write code for the independent terms.
    if unfolding.normalized:
        code += "\n    # Independent terms:\n"
        for term in rho00_terms:
            coef_list = independent_get_coefficients(*term)
            for coef in coef_list:
                code += term_code(*coef, linear=False)
    # We finish the code.
    if True:
        if matrix_form:
            if unfolding.normalized:
                code += "    A *= fact\n"
                code += "    b *= fact\n"
                code += "    return A, b\n"
            else:
                code += "    A *= fact\n"
                code += "    return A\n"
        else:
            code += "    rhs *= fact\n"
            code += "    return rhs\n"
    # We write the code to file if provided, and execute it.
    if True:
        if file_name is not None:
            f = file(file_name+".py", "w")
            f.write(code)
            f.close()

        rabi_terms = code
        if not return_code:
            exec rabi_terms
    return rabi_terms
def fast_rabi_terms ( Ep , epsilonp , rm , xi , theta , unfolding , matrix_form = False , file_name = None , return_code = False ) : if not unfolding . lower_triangular : mes = "It is very inefficient to solve using all components of the " mes += "density matrix. Better set lower_triangular=True in Unfolding." raise NotImplementedError ( mes ) if matrix_form and ( not unfolding . real ) and ( unfolding . lower_triangular ) : mes = "It is not possible to express the equations in matrix form " mes += "for complex lower triangular components only." raise ValueError ( mes ) Nl = len ( Ep ) Ne = unfolding . Ne if True : try : Ep = np . array ( [ complex ( Ep [ l ] ) for l in range ( Nl ) ] ) variable_Ep = False except : variable_Ep = True try : epsilonp = [ np . array ( [ complex ( epsilonp [ l ] [ i ] ) for i in range ( 3 ) ] ) for l in range ( Nl ) ] variable_epsilonp = False except : variable_epsilonp = True if True : Nrho = unfolding . Nrho IJ = unfolding . IJ normalized = unfolding . normalized lower_triangular = unfolding . lower_triangular Em = [ ii . conjugate ( ) for ii in Ep ] if variable_epsilonp : epsilonm = Vector3D ( epsilonp . args [ 0 ] . conjugate ( ) ) else : epsilonm = [ ii . conjugate ( ) for ii in epsilonp ] rm = np . array ( [ [ [ complex ( rm [ k ] [ i , j ] ) for j in range ( Ne ) ] for i in range ( Ne ) ] for k in range ( 3 ) ] ) rp = np . array ( [ rm [ ii ] . T . conjugate ( ) for ii in range ( 3 ) ] ) rm_aux = Vector3D ( IndexedBase ( "rm" , shape = ( Ne , Ne ) ) ) rp_aux = Vector3D ( IndexedBase ( "rp" , shape = ( Ne , Ne ) ) ) rho = define_density_matrix ( Ne , explicitly_hermitian = lower_triangular , normalized = normalized ) if variable_epsilonp : Omega = [ [ sum ( [ xi [ l , i , j ] * ( Ep [ l ] * dot ( epsilonp [ l ] , rm_aux [ i , j ] ) * delta_greater ( i , j ) + Em [ l ] * dot ( epsilonm [ l ] , rp_aux [ i , j ] ) * delta_lesser ( i , j ) ) for l in range ( Nl ) ] ) for j in range ( Ne ) ] for i in range ( Ne ) ] else : Omega = [ [ sum ( [ xi [ l , i , j ] * ( Ep [ l ] * dot ( epsilonp [ l ] , rm [ : , i , j ] ) + Em [ l ] * dot ( epsilonm [ l ] , rp [ : , i , j ] ) ) for l in range ( Nl ) ] ) for j in range ( Ne ) ] for i in range ( Ne ) ] if True : code = "" code += "def rabi_terms(" if not matrix_form : code += "rho, " if variable_Ep : code += "Ep, " if variable_epsilonp : code += "epsilonp, " if code [ - 2 : ] == ", " : code = code [ : - 2 ] code += "):\n" code += '    r"""A fast calculation of the Rabi terms."""\n' if True : if unfolding . real : code += "    fact = " + str ( e_num / hbar_num ) + "\n\n" else : code += "    fact = " + str ( 1j * e_num / hbar_num ) + "\n\n" if variable_epsilonp : np . set_printoptions ( threshold = np . nan ) code += "    rm = np." + rm . __repr__ ( ) + "\n\n" code += "    rp = np." + rp . __repr__ ( ) + "\n\n" code += "    def dot(epsilon, rij):\n" code += "        return epsilon[0]*rij[0]" code += " + epsilon[1]*rij[1]" code += " + epsilon[2]*rij[2]\n\n" np . set_printoptions ( threshold = 1000 ) if matrix_form : code += "    A = np.zeros((" + str ( Nrho ) + ", " + str ( Nrho ) if not unfolding . real : code += "), complex)\n\n" else : code += "))\n\n" if unfolding . normalized : code += "    b = np.zeros((" + str ( Nrho ) if not unfolding . real : code += "), complex)\n\n" else : code += "))\n\n" else : code += "    rhs = np.zeros((" + str ( Nrho ) if not unfolding . real : code += "), complex)\n\n" else : code += "))\n\n" rho00_terms = [ ] for mu in range ( Nrho ) : s , i , j = IJ ( mu ) for l in range ( Nl ) : for k in range ( Ne ) : if xi [ l , i , k ] == 1 : u = k ; v = j args = ( 0.5 * Omega [ i ] [ k ] , rho [ k , j ] , s , i , j , k , u , v , unfolding , matrix_form ) term_list = linear_get_coefficients ( * args ) for term in term_list : code += term_code ( * term ) if k == 0 and j == 0 : rho00_terms += [ args ] if xi [ l , k , j ] == 1 : u = i ; v = k args = ( - 0.5 * Omega [ k ] [ j ] , rho [ i , k ] , s , i , j , k , u , v , unfolding , matrix_form ) term_list = linear_get_coefficients ( * args ) for term in term_list : code += term_code ( * term ) if i == 0 and k == 0 : rho00_terms += [ args ] if unfolding . normalized : code += "\n    # Independent terms:\n" for term in rho00_terms : coef_list = independent_get_coefficients ( * term ) for coef in coef_list : code += term_code ( * coef , linear = False ) if True : if matrix_form : if unfolding . normalized : code += "    A *= fact\n" code += "    b *= fact\n" code += "    return A, b\n" else : code += "    A *= fact\n" code += "    return A\n" else : code += "    rhs *= fact\n" code += "    return rhs\n" if True : if file_name is not None : f = file ( file_name + ".py" , "w" ) f . write ( code ) f . close ( ) rabi_terms = code if not return_code : exec rabi_terms return rabi_terms
r"""Return a fast function that returns the Rabi frequency terms.

    We test a basic two-level system.

    >>> import numpy as np
    >>> from scipy.constants import physical_constants
    >>> from sympy import Matrix
    >>> from fast.electric_field import electric_field_amplitude_top
    >>> from fast.symbolic import define_laser_variables, polarization_vector

    >>> Ne = 2
    >>> Nl = 1
    >>> a0 = physical_constants["Bohr radius"][0]
    >>> rm = [np.array([[0, 0], [a0, 0]]),
    ...       np.array([[0, 0], [0, 0]]),
    ...       np.array([[0, 0], [0, 0]])]
    >>> xi = np.array([[[0, 1], [1, 0]]])
    >>> theta = phase_transformation(Ne, Nl, rm, xi)

    We define symbolic variables to be used as token arguments.
    >>> Eps = [electric_field_amplitude_top(1e-3, 1e-3, 1, "SI")]
    >>> Ep, omega_laser = define_laser_variables(Nl)
    >>> epsilonps = [polarization_vector(0, 0, 0, 0, 1)]

    An map to unfold the density matrix.
    >>> unfolding = Unfolding(Ne, True, True, True)

    We obtain a function to calculate Rabi frequency terms.
    >>> rabi_terms = fast_rabi_terms(Ep, epsilonps, rm, xi, theta, unfolding)

    Apply this to a density matrix.
    >>> rhos = np.array([[0.6, 3+2j],
    ...                  [3-2j, 0.4]])
    >>> rhosv = unfolding(rhos)
    >>> rhs_rabi = rabi_terms(rhosv, Eps)
    >>> print(rhs_rabi)
    [-55680831.474          0.       2784041.5737]