Python Odedata.num_expect方法代码示例

# 需要导入模块: from qutip.odedata import Odedata [as 别名]
# 或者: from qutip.odedata.Odedata import num_expect [as 别名]
def _generic_ode_solve(r, rho0, tlist, e_ops, opt, progress_bar):
    """
    Internal function for solving ME.
    """

    #
    # prepare output array
    #
    n_tsteps = len(tlist)
    dt = tlist[1] - tlist[0]
    e_sops_data = []

    output = Odedata()
    output.solver = "mesolve"
    output.times = tlist

    if opt.store_states:
        output.states = []

    if isinstance(e_ops, types.FunctionType):
        n_expt_op = 0
        expt_callback = True

    elif isinstance(e_ops, list):

        n_expt_op = len(e_ops)
        expt_callback = False

        if n_expt_op == 0:
            # fall back on storing states
            output.states = []
            opt.store_states = True
        else:
            output.expect = []
            output.num_expect = n_expt_op
            for op in e_ops:
                e_sops_data.append(spre(op).data)
                if op.isherm and rho0.isherm:
                    output.expect.append(np.zeros(n_tsteps))
                else:
                    output.expect.append(np.zeros(n_tsteps, dtype=complex))

    else:
        raise TypeError("Expectation parameter must be a list or a function")

    #
    # start evolution
    #
    progress_bar.start(n_tsteps)

    rho = Qobj(rho0)

    for t_idx, t in enumerate(tlist):
        progress_bar.update(t_idx)

        if not r.successful():
            break

        if opt.store_states or expt_callback:
            rho.data = vec2mat(r.y)

            if opt.store_states:
                output.states.append(Qobj(rho))

            if expt_callback:
                # use callback method
                e_ops(t, rho)

        for m in range(n_expt_op):
            if output.expect[m].dtype == complex:
                output.expect[m][t_idx] = expect_rho_vec(e_sops_data[m], r.y)
            else:
                output.expect[m][t_idx] = np.real(
                    expect_rho_vec(e_sops_data[m], r.y))

        r.integrate(r.t + dt)

    progress_bar.finished()

    if not opt.rhs_reuse and odeconfig.tdname is not None:
        try:
            os.remove(odeconfig.tdname + ".pyx")
        except:
            pass

    if opt.store_final_state:
        rho.data = vec2mat(r.y)
        output.final_state = Qobj(rho)

    return output

# 需要导入模块: from qutip.odedata import Odedata [as 别名]
# 或者: from qutip.odedata.Odedata import num_expect [as 别名]
def _generic_ode_solve(r, psi0, tlist, e_ops, opt, progress_bar,
                       state_norm_func=None, dims=None):
    """
    Internal function for solving ODEs.
    """

    #
    # prepare output array
    #
    n_tsteps = len(tlist)
    output = Odedata()
    output.solver = "sesolve"
    output.times = tlist

    if opt.store_states:
        output.states = []

    if isinstance(e_ops, types.FunctionType):
        n_expt_op = 0
        expt_callback = True

    elif isinstance(e_ops, list):

        n_expt_op = len(e_ops)
        expt_callback = False

        if n_expt_op == 0:
            # fallback on storing states
            output.states = []
            opt.store_states = True
        else:
            output.expect = []
            output.num_expect = n_expt_op
            for op in e_ops:
                if op.isherm:
                    output.expect.append(np.zeros(n_tsteps))
                else:
                    output.expect.append(np.zeros(n_tsteps, dtype=complex))
    else:
        raise TypeError("Expectation parameter must be a list or a function")

    #
    # start evolution
    #
    progress_bar.start(n_tsteps)

    dt = np.diff(tlist)
    for t_idx, t in enumerate(tlist):
        progress_bar.update(t_idx)

        if not r.successful():
            break

        if state_norm_func:
            data = r.y / state_norm_func(r.y)
            r.set_initial_value(data, r.t)

        if opt.store_states:
            output.states.append(Qobj(r.y,dims=dims))

        if expt_callback:
            # use callback method
            e_ops(t, Qobj(r.y, dims=psi0.dims))

        for m in range(n_expt_op):
            output.expect[m][t_idx] = cy_expect_psi(e_ops[m].data, r.y, e_ops[m].isherm)

        if t_idx < n_tsteps - 1:
            r.integrate(r.t + dt[t_idx])

    progress_bar.finished()

    if not opt.rhs_reuse and odeconfig.tdname is not None:
        try:
            os.remove(odeconfig.tdname + ".pyx")
        except:
            pass

    if opt.store_final_state:
        output.final_state = Qobj(r.y)

    return output

# 需要导入模块: from qutip.odedata import Odedata [as 别名]
# 或者: from qutip.odedata.Odedata import num_expect [as 别名]
def floquet_markov_mesolve(R, ekets, rho0, tlist, e_ops, f_modes_table=None,
                           options=None, floquet_basis=True):
    """
    Solve the dynamics for the system using the Floquet-Markov master equation.
    """

    if options is None:
        opt = Odeoptions()
    else:
        opt = options

    if opt.tidy:
        R.tidyup()

    #
    # check initial state
    #
    if isket(rho0):
        # Got a wave function as initial state: convert to density matrix.
        rho0 = ket2dm(rho0)

    #
    # prepare output array
    #
    n_tsteps = len(tlist)
    dt = tlist[1] - tlist[0]

    output = Odedata()
    output.solver = "fmmesolve"
    output.times = tlist

    if isinstance(e_ops, FunctionType):
        n_expt_op = 0
        expt_callback = True

    elif isinstance(e_ops, list):

        n_expt_op = len(e_ops)
        expt_callback = False

        if n_expt_op == 0:
            output.states = []
        else:
            if not f_modes_table:
                raise TypeError("The Floquet mode table has to be provided " +
                                "when requesting expectation values.")

            output.expect = []
            output.num_expect = n_expt_op
            for op in e_ops:
                if op.isherm:
                    output.expect.append(np.zeros(n_tsteps))
                else:
                    output.expect.append(np.zeros(n_tsteps, dtype=complex))

    else:
        raise TypeError("Expectation parameter must be a list or a function")

    #
    # transform the initial density matrix to the eigenbasis: from
    # computational basis to the floquet basis
    #
    if ekets is not None:
        rho0 = rho0.transform(ekets, True)

    #
    # setup integrator
    #
    initial_vector = mat2vec(rho0.full())
    r = scipy.integrate.ode(cy_ode_rhs)
    r.set_f_params(R.data.data, R.data.indices, R.data.indptr)
    r.set_integrator('zvode', method=opt.method, order=opt.order,
                     atol=opt.atol, rtol=opt.rtol, max_step=opt.max_step)
    r.set_initial_value(initial_vector, tlist[0])

    #
    # start evolution
    #
    rho = Qobj(rho0)

    t_idx = 0
    for t in tlist:
        if not r.successful():
            break

        rho.data = vec2mat(r.y)

        if expt_callback:
            # use callback method
            if floquet_basis:
                e_ops(t, Qobj(rho))
            else:
                f_modes_table_t, T = f_modes_table
                f_modes_t = floquet_modes_t_lookup(f_modes_table_t, t, T)
                e_ops(t, Qobj(rho).transform(f_modes_t, False))
        else:
            # calculate all the expectation values, or output rho if
            # no operators
            if n_expt_op == 0:
                if floquet_basis:
#.........这里部分代码省略.........

# 需要导入模块: from qutip.odedata import Odedata [as 别名]
# 或者: from qutip.odedata.Odedata import num_expect [as 别名]
def fsesolve(H, psi0, tlist, e_ops=[], T=None, args={}, Tsteps=100):
    """
    Solve the Schrodinger equation using the Floquet formalism.

    Parameters
    ----------

    H : :class:`qutip.qobj.Qobj`
        System Hamiltonian, time-dependent with period `T`.

    psi0 : :class:`qutip.qobj`
        Initial state vector (ket).

    tlist : *list* / *array*
        list of times for :math:`t`.

    e_ops : list of :class:`qutip.qobj` / callback function
        list of operators for which to evaluate expectation values. If this
        list is empty, the state vectors for each time in `tlist` will be
        returned instead of expectation values.

    T : float
        The period of the time-dependence of the hamiltonian.

    args : dictionary
        Dictionary with variables required to evaluate H.

    Tsteps : integer
        The number of time steps in one driving period for which to
        precalculate the Floquet modes. `Tsteps` should be an even number.

    Returns
    -------

    output : :class:`qutip.odedata.Odedata`

        An instance of the class :class:`qutip.odedata.Odedata`, which
        contains either an *array* of expectation values or an array of
        state vectors, for the times specified by `tlist`.
    """

    if not T:
        # assume that tlist span exactly one period of the driving
        T = tlist[-1]

    # find the floquet modes for the time-dependent hamiltonian
    f_modes_0, f_energies = floquet_modes(H, T, args)

    # calculate the wavefunctions using the from the floquet modes
    f_modes_table_t = floquet_modes_table(f_modes_0, f_energies,
                                          np.linspace(0, T, Tsteps + 1),
                                          H, T, args)

    # setup Odedata for storing the results
    output = Odedata()
    output.times = tlist
    output.solver = "fsesolve"

    if isinstance(e_ops, FunctionType):
        output.num_expect = 0
        expt_callback = True

    elif isinstance(e_ops, list):

        output.num_expect = len(e_ops)
        expt_callback = False

        if output.num_expect == 0:
            output.states = []
        else:
            output.expect = []
            for op in e_ops:
                if op.isherm:
                    output.expect.append(np.zeros(len(tlist)))
                else:
                    output.expect.append(np.zeros(len(tlist), dtype=complex))

    else:
        raise TypeError("e_ops must be a list Qobj or a callback function")

    psi0_fb = psi0.transform(f_modes_0, True)
    for t_idx, t in enumerate(tlist):
        f_modes_t = floquet_modes_t_lookup(f_modes_table_t, t, T)
        f_states_t = floquet_states(f_modes_t, f_energies, t)
        psi_t = psi0_fb.transform(f_states_t, False)

        if expt_callback:
            # use callback method
            e_ops(t, psi_t)
        else:
            # calculate all the expectation values, or output psi if
            # no expectation value operators where defined
            if output.num_expect == 0:
                output.states.append(Qobj(psi_t))
            else:
                for e_idx, e in enumerate(e_ops):
                    output.expect[e_idx][t_idx] = expect(e, psi_t)

    return output

# 需要导入模块: from qutip.odedata import Odedata [as 别名]
# 或者: from qutip.odedata.Odedata import num_expect [as 别名]

#.........这里部分代码省略.........
        ntraj = options.ntraj

    if psi0.type != 'ket':
        raise Exception("Initial state must be a state vector.")
    odeconfig.options = options
    # set num_cpus to the value given in qutip.settings
    # if none in Odeoptions
    if not odeconfig.options.num_cpus:
        odeconfig.options.num_cpus = qutip.settings.num_cpus
    # set initial value data
    if options.tidy:
        odeconfig.psi0 = psi0.tidyup(options.atol).full()
    else:
        odeconfig.psi0 = psi0.full()
    odeconfig.psi0_dims = psi0.dims
    odeconfig.psi0_shape = psi0.shape
    # set general items
    odeconfig.tlist = tlist
    if isinstance(ntraj, (list, np.ndarray)):
        raise Exception("ntraj as list argument is not supported.")
    else:
        odeconfig.ntraj = ntraj
        # ntraj_list = [ntraj]
    # set norm finding constants
    odeconfig.norm_tol = options.norm_tol
    odeconfig.norm_steps = options.norm_steps

    if not options.rhs_reuse:
        odeconfig.soft_reset()
        # no time dependence
        odeconfig.tflag = 0
        # check for collapse operators
        if len(c_ops) > 0:
            odeconfig.cflag = 1
        else:
            odeconfig.cflag = 0
        # Configure data
        _mc_data_config(H, psi0, [], c_ops, [], [], e_ops, options, odeconfig)

    # Load Monte Carlo class
    mc = _MC_class()
    # Set solver type
    if (options.method == 'adams'):
        mc.mf = 10
    elif (options.method == 'bdf'):
        mc.mf = 22
    else:
        if debug:
            print('Unrecognized method for ode solver, using "adams".')
        mc.mf = 10
    # store ket and density matrix dims and shape for convenience
    mc.psi0_dims = psi0.dims
    mc.psi0_shape = psi0.shape
    mc.dm_dims = (psi0 * psi0.dag()).dims
    mc.dm_shape = (psi0 * psi0.dag()).shape
    # use sparse density matrices during computation?
    mc.sparse_dms = sparse_dms
    # run in serial?
    mc.serial_run = serial or (ntraj == 1)
    # are we doing a partial trace for returned states?
    mc.ptrace_sel = ptrace_sel
    if (ptrace_sel != []):
        if debug:
            print("ptrace_sel set to " + str(ptrace_sel))
            print("We are using dense density matrices during computation " +
                  "when performing partial trace. Setting sparse_dms = False")
            print("This feature is experimental.")
        mc.sparse_dms = False
        mc.dm_dims = psi0.ptrace(ptrace_sel).dims
        mc.dm_shape = psi0.ptrace(ptrace_sel).shape
    if (calc_entropy):
        if (ptrace_sel == []):
            if debug:
                print("calc_entropy = True, but ptrace_sel = []. Please set " +
                     "a list of components to keep when calculating average " +
                     "entropy of reduced density matrix in ptrace_sel. " +
                     "Setting calc_entropy = False.")
            calc_entropy = False
        mc.calc_entropy = calc_entropy

    # construct output Odedata object
    output = Odedata()

    # Run
    mc.run()
    output.states = mc.sol.states
    output.expect = mc.sol.expect
    output.col_times = mc.sol.col_times
    output.col_which = mc.sol.col_which
    if (hasattr(mc.sol, 'entropy')):
        output.entropy = mc.sol.entropy

    output.solver = 'Fortran 90 Monte Carlo solver'
    # simulation parameters
    output.times = odeconfig.tlist
    output.num_expect = odeconfig.e_num
    output.num_collapse = odeconfig.c_num
    output.ntraj = odeconfig.ntraj

    return output

# 需要导入模块: from qutip.odedata import Odedata [as 别名]
# 或者: from qutip.odedata.Odedata import num_expect [as 别名]
def _generic_ode_solve(r, psi0, tlist, expt_ops, opt,
                       state_vectorize, state_norm_func=None):
    """
    Internal function for solving ODEs.
    """

    #
    # prepare output array
    #
    n_tsteps = len(tlist)
    dt = tlist[1] - tlist[0]

    output = Odedata()
    output.solver = "mesolve"
    output.times = tlist

    if isinstance(expt_ops, types.FunctionType):
        n_expt_op = 0
        expt_callback = True

    elif isinstance(expt_ops, list):

        n_expt_op = len(expt_ops)
        expt_callback = False

        if n_expt_op == 0:
            output.states = []
        else:
            output.expect = []
            output.num_expect = n_expt_op
            for op in expt_ops:
                if op.isherm and psi0.isherm:
                    output.expect.append(np.zeros(n_tsteps))
                else:
                    output.expect.append(np.zeros(n_tsteps, dtype=complex))

    else:
        raise TypeError("Expectation parameter must be a list or a function")

    #
    # start evolution
    #
    psi = Qobj(psi0)

    t_idx = 0
    for t in tlist:
        if not r.successful():
            break

        if state_norm_func:
            psi.data = state_vectorize(r.y)
            state_norm = state_norm_func(psi.data)
            psi.data = psi.data / state_norm
            r.set_initial_value(r.y / state_norm, r.t)
        else:
            psi.data = state_vectorize(r.y)

        if expt_callback:
            # use callback method
            expt_ops(t, Qobj(psi))
        else:
            # calculate all the expectation values,
            # or output rho if no operators
            if n_expt_op == 0:
                output.states.append(Qobj(psi))  # copy psi/rho
            else:
                for m in range(0, n_expt_op):
                    output.expect[m][t_idx] = expect(expt_ops[m], psi)

        r.integrate(r.t + dt)
        t_idx += 1

    if not opt.rhs_reuse and odeconfig.tdname is not None:
        try:
            os.remove(odeconfig.tdname + ".pyx")
        except:
            pass

    return output

# 需要导入模块: from qutip.odedata import Odedata [as 别名]
# 或者: from qutip.odedata.Odedata import num_expect [as 别名]

#.........这里部分代码省略.........
        #check if running in iPython and using Cython compiling (then no GUI to work around error)
        if odeconfig.options.gui and odeconfig.tflag in array([1,10,11]):
            try:
                __IPYTHON__
            except:
                pass
            else:
                odeconfig.options.gui=False    
        if qutip.settings.qutip_gui=="NONE":
            odeconfig.options.gui=False

        #check for collapse operators
        if c_terms>0:
            odeconfig.cflag=1
        else:
            odeconfig.cflag=0
    
        #Configure data
        _mc_data_config(H,psi0,h_stuff,c_ops,c_stuff,args,e_ops,options)
        if odeconfig.tflag in array([1,10,11]): #compile time-depdendent RHS code
            os.environ['CFLAGS'] = '-O3 -w'
            import pyximport
            pyximport.install(setup_args={'include_dirs':[numpy.get_include()]})
            if odeconfig.tflag in array([1,11]):
                code = compile('from '+odeconfig.tdname+' import cyq_td_ode_rhs,col_spmv,col_expect', '<string>', 'exec')
                exec(code, globals())
                odeconfig.tdfunc=cyq_td_ode_rhs
                odeconfig.colspmv=col_spmv
                odeconfig.colexpect=col_expect
            else:
                code = compile('from '+odeconfig.tdname+' import cyq_td_ode_rhs', '<string>', 'exec')
                exec(code, globals())
                odeconfig.tdfunc=cyq_td_ode_rhs
            try:
                os.remove(odeconfig.tdname+".pyx")
            except:
                print("Error removing pyx file.  File not found.")
        elif odeconfig.tflag==0:
            odeconfig.tdfunc=cyq_ode_rhs
    else:#setup args for new parameters when rhs_reuse=True and tdfunc is given
        #string based
        if odeconfig.tflag in array([1,10,11]):
            if any(args):
                odeconfig.c_args=[]
                arg_items=args.items()
                for k in range(len(args)):
                    odeconfig.c_args.append(arg_items[k][1])
        #function based
        elif odeconfig.tflag in array([2,3,20,22]):
            odeconfig.h_func_args=args
    
    
    #load monte-carlo class
    mc=_MC_class()
    #RUN THE SIMULATION
    mc.run()
    
    
    #AFTER MCSOLVER IS DONE --------------------------------------
    
    
    
    #-------COLLECT AND RETURN OUTPUT DATA IN ODEDATA OBJECT --------------#
    output=Odedata()
    output.solver='mcsolve'
    #state vectors
    if mc.psi_out is not None and odeconfig.options.mc_avg and odeconfig.cflag:
        output.states=parfor(_mc_dm_avg,mc.psi_out.T)
    elif mc.psi_out is not None:
        output.states=mc.psi_out
    #expectation values
    elif mc.expect_out is not None and odeconfig.cflag and odeconfig.options.mc_avg:#averaging if multiple trajectories
        if isinstance(ntraj,int):
            output.expect=mean(mc.expect_out,axis=0)
        elif isinstance(ntraj,(list,ndarray)):
            output.expect=[]
            for num in ntraj:
                expt_data=mean(mc.expect_out[:num],axis=0)
                data_list=[]
                if any([op.isherm==False for op in e_ops]):
                    for k in range(len(e_ops)):
                        if e_ops[k].isherm:
                            data_list.append(real(expt_data[k]))
                        else:
                            data_list.append(expt_data[k])
                else:
                    data_list=[data for data in expt_data]
                output.expect.append(data_list)
    else:#no averaging for single trajectory or if mc_avg flag (Odeoptions) is off
        if mc.expect_out is not None:        
            output.expect=mc.expect_out

    #simulation parameters
    output.times=odeconfig.tlist
    output.num_expect=odeconfig.e_num
    output.num_collapse=odeconfig.c_num
    output.ntraj=odeconfig.ntraj
    output.col_times=mc.collapse_times_out
    output.col_which=mc.which_op_out
    return output