本文整理汇总了Python中qutip.qobj.isket函数的典型用法代码示例。如果您正苦于以下问题:Python isket函数的具体用法?Python isket怎么用?Python isket使用的例子?那么恭喜您, 这里精选的函数代码示例或许可以为您提供帮助。
在下文中一共展示了isket函数的15个代码示例,这些例子默认根据受欢迎程度排序。您可以为喜欢或者感觉有用的代码点赞,您的评价将有助于系统推荐出更棒的Python代码示例。
示例1: update
def update(self, rho):
"""
Calculate the probability function for the given state of an harmonic
oscillator (as density matrix)
"""
if isket(rho):
rho = ket2dm(rho)
self.data = np.zeros(len(self.xvecs[0]), dtype=complex)
M, N = rho.shape
for m in range(M):
k_m = pow(self.omega / pi, 0.25) / \
sqrt(2 ** m * factorial(m)) * \
exp(-self.xvecs[0] ** 2 / 2.0) * \
np.polyval(hermite(m), self.xvecs[0])
for n in range(N):
k_n = pow(self.omega / pi, 0.25) / \
sqrt(2 ** n * factorial(n)) * \
exp(-self.xvecs[0] ** 2 / 2.0) * \
np.polyval(hermite(n), self.xvecs[0])
self.data += np.conjugate(k_n) * k_m * rho.data[m, n]
示例2: _correlation_es_2t
def _correlation_es_2t(H, state0, tlist, taulist, c_ops, a_op, b_op, c_op):
"""
Internal function for calculating the three-operator two-time
correlation function:
<A(t)B(t+tau)C(t)>
using an exponential series solver.
"""
# the solvers only work for positive time differences and the correlators
# require positive tau
if state0 is None:
rho0 = steadystate(H, c_ops)
tlist = [0]
elif isket(state0):
rho0 = ket2dm(state0)
else:
rho0 = state0
if debug:
print(inspect.stack()[0][3])
# contruct the Liouvillian
L = liouvillian(H, c_ops)
corr_mat = np.zeros([np.size(tlist), np.size(taulist)], dtype=complex)
solES_t = ode2es(L, rho0)
# evaluate the correlation function
for t_idx in range(len(tlist)):
rho_t = esval(solES_t, [tlist[t_idx]])
solES_tau = ode2es(L, c_op * rho_t * a_op)
corr_mat[t_idx, :] = esval(expect(b_op, solES_tau), taulist)
return corr_mat
示例3: _subsystem_apply_reference
def _subsystem_apply_reference(state, channel, mask):
if isket(state):
state = ket2dm(state)
if isoper(channel):
full_oper = tensor([channel if mask[j]
else qeye(state.dims[0][j])
for j in range(len(state.dims[0]))])
return full_oper * state * full_oper.dag()
else:
# Go to Choi, then Kraus
# chan_mat = array(channel.data.todense())
choi_matrix = super_to_choi(channel)
vals, vecs = eig(choi_matrix.full())
vecs = list(map(array, zip(*vecs)))
kraus_list = [sqrt(vals[j]) * vec2mat(vecs[j])
for j in range(len(vals))]
# Kraus operators to be padded with identities:
k_qubit_kraus_list = product(kraus_list, repeat=sum(mask))
rho_out = Qobj(inpt=zeros(state.shape), dims=state.dims)
for operator_iter in k_qubit_kraus_list:
operator_iter = iter(operator_iter)
op_iter_list = [next(operator_iter).conj().T if mask[j]
else qeye(state.dims[0][j])
for j in range(len(state.dims[0]))]
full_oper = tensor(list(map(Qobj, op_iter_list)))
rho_out = rho_out + full_oper * state * full_oper.dag()
return Qobj(rho_out)
示例4: _sesolve_const
def _sesolve_const(H, psi0, tlist, e_ops, args, opt, progress_bar):
"""!
Evolve the wave function using an ODE solver
"""
if debug:
print(inspect.stack()[0][3])
if not isket(psi0):
raise TypeError("psi0 must be a ket")
#
# setup integrator.
#
initial_vector = psi0.full().ravel()
r = scipy.integrate.ode(cy_ode_rhs)
L = -1.0j * H
r.set_f_params(L.data.data, L.data.indices, L.data.indptr) # cython RHS
r.set_integrator('zvode', method=opt.method, order=opt.order,
atol=opt.atol, rtol=opt.rtol, nsteps=opt.nsteps,
first_step=opt.first_step, min_step=opt.min_step,
max_step=opt.max_step)
r.set_initial_value(initial_vector, tlist[0])
#
# call generic ODE code
#
return _generic_ode_solve(r, psi0, tlist, e_ops, opt,
progress_bar, norm, dims=psi0.dims)
示例5: _sesolve_func_td
def _sesolve_func_td(H_func, psi0, tlist, e_ops, args, opt, progress_bar):
"""!
Evolve the wave function using an ODE solver with time-dependent
Hamiltonian.
"""
if debug:
print(inspect.stack()[0][3])
if not isket(psi0):
raise TypeError("psi0 must be a ket")
#
# setup integrator
#
new_args = None
if type(args) is dict:
new_args = {}
for key in args:
if isinstance(args[key], Qobj):
new_args[key] = args[key].data
else:
new_args[key] = args[key]
elif type(args) is list or type(args) is tuple:
new_args = []
for arg in args:
if isinstance(arg, Qobj):
new_args.append(arg.data)
else:
new_args.append(arg)
if type(args) is tuple:
new_args = tuple(new_args)
else:
if isinstance(args, Qobj):
new_args = args.data
else:
new_args = args
initial_vector = psi0.full().ravel()
if not opt.rhs_with_state:
r = scipy.integrate.ode(cy_ode_psi_func_td)
else:
r = scipy.integrate.ode(cy_ode_psi_func_td_with_state)
r.set_integrator('zvode', method=opt.method, order=opt.order,
atol=opt.atol, rtol=opt.rtol, nsteps=opt.nsteps,
first_step=opt.first_step, min_step=opt.min_step,
max_step=opt.max_step)
r.set_initial_value(initial_vector, tlist[0])
r.set_f_params(H_func, new_args)
#
# call generic ODE code
#
return _generic_ode_solve(r, psi0, tlist, e_ops, opt, progress_bar, norm,
dims=psi0.dims)
示例6: qfunc
def qfunc(state, xvec, yvec, g=sqrt(2)):
"""Q-function of a given state vector or density matrix
at points `xvec + i * yvec`.
Parameters
----------
state : qobj
A state vector or density matrix.
xvec : array_like
x-coordinates at which to calculate the Wigner function.
yvec : array_like
y-coordinates at which to calculate the Wigner function.
g : float
Scaling factor for `a = 0.5 * g * (x + iy)`, default `g = sqrt(2)`.
Returns
--------
Q : array
Values representing the Q-function calculated over the specified range
[xvec,yvec].
"""
X, Y = meshgrid(xvec, yvec)
amat = 0.5 * g * (X + Y * 1j)
if not (isoper(state) or isket(state)):
raise TypeError('Invalid state operand to qfunc.')
qmat = zeros(size(amat))
if isket(state):
qmat = _qfunc_pure(state, amat)
elif isoper(state):
d, v = la.eig(state.full())
# d[i] = eigenvalue i
# v[:,i] = eigenvector i
qmat = zeros(np.shape(amat))
for k in arange(0, len(d)):
qmat1 = _qfunc_pure(v[:, k], amat)
qmat += real(d[k] * qmat1)
qmat = 0.25 * qmat * g ** 2
return qmat
示例7: _mesolve_const
def _mesolve_const(H, rho0, tlist, c_op_list, e_ops, args, opt,
progress_bar):
"""
Evolve the density matrix using an ODE solver, for constant hamiltonian
and collapse operators.
"""
if debug:
print(inspect.stack()[0][3])
#
# check initial state
#
if isket(rho0):
# if initial state is a ket and no collapse operator where given,
# fall back on the unitary schrodinger equation solver
if len(c_op_list) == 0 and isoper(H):
return _sesolve_const(H, rho0, tlist, e_ops, args, opt,
progress_bar)
# Got a wave function as initial state: convert to density matrix.
rho0 = ket2dm(rho0)
#
# construct liouvillian
#
if opt.tidy:
H = H.tidyup(opt.atol)
L = liouvillian(H, c_op_list)
#
# setup integrator
#
initial_vector = mat2vec(rho0.full()).ravel('F')
if issuper(rho0):
r = scipy.integrate.ode(_ode_super_func)
r.set_f_params(L.data)
else:
if opt.use_openmp and L.data.nnz >= qset.openmp_thresh:
r = scipy.integrate.ode(cy_ode_rhs_openmp)
r.set_f_params(L.data.data, L.data.indices, L.data.indptr,
opt.openmp_threads)
else:
r = scipy.integrate.ode(cy_ode_rhs)
r.set_f_params(L.data.data, L.data.indices, L.data.indptr)
# r = scipy.integrate.ode(_ode_rho_test)
# r.set_f_params(L.data)
r.set_integrator('zvode', method=opt.method, order=opt.order,
atol=opt.atol, rtol=opt.rtol, nsteps=opt.nsteps,
first_step=opt.first_step, min_step=opt.min_step,
max_step=opt.max_step)
r.set_initial_value(initial_vector, tlist[0])
#
# call generic ODE code
#
return _generic_ode_solve(r, rho0, tlist, e_ops, opt, progress_bar)
示例8: _sesolve_list_func_td
def _sesolve_list_func_td(H_list, psi0, tlist, expt_ops, args, opt,
progress_bar):
"""
Internal function for solving the master equation. See mesolve for usage.
"""
if debug:
print(inspect.stack()[0][3])
#
# check initial state
#
if not isket(psi0):
raise TypeError("The unitary solver requires a ket as initial state")
#
# construct liouvillian in list-function format
#
L_list = []
constant_func = lambda x, y: 1.0
# add all hamitonian terms to the lagrangian list
for h_spec in H_list:
if isinstance(h_spec, Qobj):
h = h_spec
h_coeff = constant_func
elif isinstance(h_spec, list):
h = h_spec[0]
h_coeff = h_spec[1]
else:
raise TypeError("Incorrect specification of time-dependent " +
"Hamiltonian (expected callback function)")
L_list.append([-1j * h.data, h_coeff])
L_list_and_args = [L_list, args]
#
# setup integrator
#
initial_vector = psi0.full()
r = scipy.integrate.ode(psi_list_td)
r.set_integrator('zvode', method=opt.method, order=opt.order,
atol=opt.atol, rtol=opt.rtol, nsteps=opt.nsteps,
first_step=opt.first_step, min_step=opt.min_step,
max_step=opt.max_step)
r.set_initial_value(initial_vector, tlist[0])
r.set_f_params(L_list_and_args)
#
# call generic ODE code
#
return _generic_ode_solve(r, psi0, tlist, expt_ops, opt, progress_bar)
示例9: plot_wigner_fock_distribution
def plot_wigner_fock_distribution(rho, fig=None, axes=None, figsize=(8, 4),
cmap=None, alpha_max=7.5, colorbar=False,
method='iterative'):
"""
Plot the Fock distribution and the Wigner function for a density matrix
(or ket) that describes an oscillator mode.
Parameters
----------
rho : :class:`qutip.qobj.Qobj`
The density matrix (or ket) of the state to visualize.
fig : a matplotlib Figure instance
The Figure canvas in which the plot will be drawn.
axes : a list of two matplotlib axes instances
The axes context in which the plot will be drawn.
figsize : (width, height)
The size of the matplotlib figure (in inches) if it is to be created
(that is, if no 'fig' and 'ax' arguments are passed).
cmap : a matplotlib cmap instance
The colormap.
alpha_max : float
The span of the x and y coordinates (both [-alpha_max, alpha_max]).
colorbar : bool
Whether (True) or not (False) a colorbar should be attached to the
Wigner function graph.
method : string {'iterative', 'laguerre', 'fft'}
The method used for calculating the wigner function. See the
documentation for qutip.wigner for details.
Returns
-------
fig, ax : tuple
A tuple of the matplotlib figure and axes instances used to produce
the figure.
"""
if not fig and not axes:
fig, axes = plt.subplots(1, 2, figsize=figsize)
if isket(rho):
rho = ket2dm(rho)
plot_fock_distribution(rho, fig=fig, ax=axes[0])
plot_wigner(rho, fig=fig, ax=axes[1], figsize=figsize, cmap=cmap,
alpha_max=alpha_max, colorbar=colorbar, method=method)
return fig, axes
示例10: plot_fock_distribution
def plot_fock_distribution(rho, offset=0, fig=None, ax=None,
figsize=(8, 6), title=None, unit_y_range=True):
"""
Plot the Fock distribution for a density matrix (or ket) that describes
an oscillator mode.
Parameters
----------
rho : :class:`qutip.qobj.Qobj`
The density matrix (or ket) of the state to visualize.
fig : a matplotlib Figure instance
The Figure canvas in which the plot will be drawn.
ax : a matplotlib axes instance
The axes context in which the plot will be drawn.
title : string
An optional title for the figure.
figsize : (width, height)
The size of the matplotlib figure (in inches) if it is to be created
(that is, if no 'fig' and 'ax' arguments are passed).
Returns
-------
fig, ax : tuple
A tuple of the matplotlib figure and axes instances used to produce
the figure.
"""
if not fig and not ax:
fig, ax = plt.subplots(1, 1, figsize=figsize)
if isket(rho):
rho = ket2dm(rho)
N = rho.shape[0]
ax.bar(np.arange(offset, offset + N) - .4, np.real(rho.diag()),
color="green", alpha=0.6, width=0.8)
if unit_y_range:
ax.set_ylim(0, 1)
ax.set_xlim(-.5 + offset, N + offset)
ax.set_xlabel('Fock number', fontsize=12)
ax.set_ylabel('Occupation probability', fontsize=12)
if title:
ax.set_title(title)
return fig, ax
示例11: _correlation_me_2t
def _correlation_me_2t(H, state0, tlist, taulist, c_ops, a_op, b_op, c_op,
args={}, options=Options()):
"""
Internal function for calculating the three-operator two-time
correlation function:
<A(t)B(t+tau)C(t)>
using a master equation solver.
"""
# the solvers only work for positive time differences and the correlators
# require positive tau
if state0 is None:
rho0 = steadystate(H, c_ops)
tlist = [0]
elif isket(state0):
rho0 = ket2dm(state0)
else:
rho0 = state0
if debug:
print(inspect.stack()[0][3])
rho_t = mesolve(H, rho0, tlist, c_ops, [],
args=args, options=options).states
corr_mat = np.zeros([np.size(tlist), np.size(taulist)], dtype=complex)
H_shifted, c_ops_shifted, _args = _transform_L_t_shift(H, c_ops, args)
if config.tdname:
_cython_build_cleanup(config.tdname)
rhs_clear()
for t_idx, rho in enumerate(rho_t):
if not isinstance(H, Qobj):
_args["_t0"] = tlist[t_idx]
corr_mat[t_idx, :] = mesolve(
H_shifted, c_op * rho * a_op, taulist, c_ops_shifted,
[b_op], args=_args, options=options
).expect[0]
if t_idx == 1:
options.rhs_reuse = True
if config.tdname:
_cython_build_cleanup(config.tdname)
rhs_clear()
return corr_mat
示例12: plot_qubism
def plot_qubism(ket, theme='light', how='pairs',
grid_iteration=1, legend_iteration=0,
fig=None, ax=None, figsize=(6, 6)):
"""
Qubism plot for pure states of many qudits.
Works best for spin chains, especially with even number of particles
of the same dimension.
Allows to see entanglement between first 2*k particles and the rest.
More information:
J. Rodriguez-Laguna, P. Migdal,
M. Ibanez Berganza, M. Lewenstein, G. Sierra,
"Qubism: self-similar visualization of many-body wavefunctions",
New J. Phys. 14 053028 (2012), arXiv:1112.3560,
http://dx.doi.org/10.1088/1367-2630/14/5/053028 (open access)
Parameters
----------
ket : Qobj
Pure state for plotting.
theme : 'light' (default) or 'dark'
Set coloring theme for mapping complex values into colors.
See: complex_array_to_rgb.
how : 'pairs' (default), 'pairs_skewed' or 'before_after'
Type of Qubism plotting.
Options:
'pairs' - typical coordinates,
'pairs_skewed' - for ferromagnetic/antriferromagnetic plots,
'before_after' - related to Schmidt plot (see also: plot_schmidt).
grid_iteration : int (default 1)
Helper lines to be drawn on plot.
Show tiles for 2*grid_iteration particles vs all others.
legend_iteration : int (default 0) or 'grid_iteration' or 'all'
Show labels for first 2*legend_iteration particles.
Option 'grid_iteration' sets the same number of particles
as for grid_iteration.
Option 'all' makes label for all particles.
Typically it should be 0, 1, 2 or perhaps 3.
fig : a matplotlib figure instance
The figure canvas on which the plot will be drawn.
ax : a matplotlib axis instance
The axis context in which the plot will be drawn.
figsize : (width, height)
The size of the matplotlib figure (in inches) if it is to be created
(that is, if no 'fig' and 'ax' arguments are passed).
Returns
-------
fig, ax : tuple
A tuple of the matplotlib figure and axes instances used to produce
the figure.
"""
if not isket(ket):
raise Exception("Qubism works only for pure states, i.e. kets.")
# add for dm? (perhaps a separate function, plot_qubism_dm)
if not fig and not ax:
fig, ax = plt.subplots(1, 1, figsize=figsize)
dim_list = ket.dims[0]
n = len(dim_list)
# for odd number of particles - pixels are rectangular
if n % 2 == 1:
ket = tensor(ket, Qobj([1] * dim_list[-1]))
dim_list = ket.dims[0]
n += 1
ketdata = ket.full()
if how == 'pairs':
dim_list_y = dim_list[::2]
dim_list_x = dim_list[1::2]
elif how == 'pairs_skewed':
dim_list_y = dim_list[::2]
dim_list_x = dim_list[1::2]
if dim_list_x != dim_list_y:
raise Exception("For 'pairs_skewed' pairs " +
"of dimensions need to be the same.")
elif how == 'before_after':
dim_list_y = list(reversed(dim_list[:(n // 2)]))
dim_list_x = dim_list[(n // 2):]
else:
raise Exception("No such 'how'.")
size_x = np.prod(dim_list_x)
size_y = np.prod(dim_list_y)
qub = np.zeros([size_x, size_y], dtype=complex)
for i in range(ketdata.size):
qub[_to_qubism_index_pair(i, dim_list, how=how)] = ketdata[i, 0]
#.........这里部分代码省略.........
示例13: _correlation_mc_2t
def _correlation_mc_2t(H, state0, tlist, taulist, c_ops, a_op, b_op, c_op,
args=None, options=Options()):
"""
Internal function for calculating the three-operator two-time
correlation function:
<A(t)B(t+tau)C(t)>
using a Monte Carlo solver.
"""
# the solvers only work for positive time differences and the correlators
# require positive tau
if state0 is None:
raise NotImplementedError("steady state not implemented for " +
"mc solver, please use `es` or `me`")
elif not isket(state0):
raise TypeError("state0 must be a state vector.")
psi0 = state0
if debug:
print(inspect.stack()[0][3])
psi_t_mat = mcsolve(
H, psi0, tlist, c_ops, [],
args=args, ntraj=options.ntraj[0], options=options, progress_bar=None
).states
corr_mat = np.zeros([np.size(tlist), np.size(taulist)], dtype=complex)
H_shifted, _args = _transform_H_t_shift(H, args)
# calculation of <A(t)B(t+tau)C(t)> from only knowledge of psi0 requires
# averaging over both t and tau
for t_idx in range(np.size(tlist)):
if not isinstance(H, Qobj):
_args["_t0"] = tlist[t_idx]
for trial_idx in range(options.ntraj[0]):
if isinstance(a_op, Qobj) and isinstance(c_op, Qobj):
if a_op.dag() == c_op:
# A shortcut here, requires only 1/4 the trials
chi_0 = (options.mc_corr_eps + c_op) * \
psi_t_mat[trial_idx, t_idx]
# evolve these states and calculate expectation value of B
c_tau = chi_0.norm()**2 * mcsolve(
H_shifted, chi_0/chi_0.norm(), taulist, c_ops, [b_op],
args=_args, ntraj=options.ntraj[1], options=options,
progress_bar=None
).expect[0]
# final correlation vector computed by combining the
# averages
corr_mat[t_idx, :] += c_tau/options.ntraj[1]
else:
# otherwise, need four trial wavefunctions
# (Ad+C)*psi_t, (Ad+iC)*psi_t, (Ad-C)*psi_t, (Ad-iC)*psi_t
if isinstance(a_op, Qobj):
a_op_dag = a_op.dag()
else:
# assume this is a number, ex. i.e. a_op = 1
# if this is not correct, the over-loaded addition
# operation will raise errors
a_op_dag = a_op
chi_0 = [(options.mc_corr_eps + a_op_dag +
np.exp(1j*x*np.pi/2)*c_op) *
psi_t_mat[trial_idx, t_idx]
for x in range(4)]
# evolve these states and calculate expectation value of B
c_tau = [
chi.norm()**2 * mcsolve(
H_shifted, chi/chi.norm(), taulist, c_ops, [b_op],
args=_args, ntraj=options.ntraj[1], options=options,
progress_bar=None
).expect[0]
for chi in chi_0
]
# final correlation vector computed by combining the averages
corr_mat[t_idx, :] += \
1/(4*options.ntraj[0]) * (c_tau[0] - c_tau[2] -
1j*c_tau[1] + 1j*c_tau[3])
if t_idx == 1:
options.rhs_reuse = True
return corr_mat
示例14: plot_schmidt
def plot_schmidt(ket, splitting=None,
labels_iteration=(3, 2),
theme='light',
fig=None, ax=None, figsize=(6, 6)):
"""
Plotting scheme related to Schmidt decomposition.
Converts a state into a matrix (A_ij -> A_i^j),
where rows are first particles and columns - last.
See also: plot_qubism with how='before_after' for a similar plot.
Parameters
----------
ket : Qobj
Pure state for plotting.
splitting : int
Plot for a number of first particles versus the rest.
If not given, it is (number of particles + 1) // 2.
theme : 'light' (default) or 'dark'
Set coloring theme for mapping complex values into colors.
See: complex_array_to_rgb.
labels_iteration : int or pair of ints (default (3,2))
Number of particles to be shown as tick labels,
for first (vertical) and last (horizontal) particles, respectively.
fig : a matplotlib figure instance
The figure canvas on which the plot will be drawn.
ax : a matplotlib axis instance
The axis context in which the plot will be drawn.
figsize : (width, height)
The size of the matplotlib figure (in inches) if it is to be created
(that is, if no 'fig' and 'ax' arguments are passed).
Returns
-------
fig, ax : tuple
A tuple of the matplotlib figure and axes instances used to produce
the figure.
"""
if not isket(ket):
raise Exception("Schmidt plot works only for pure states, i.e. kets.")
if not fig and not ax:
fig, ax = plt.subplots(1, 1, figsize=figsize)
dim_list = ket.dims[0]
if splitting is None:
splitting = (len(dim_list) + 1) // 2
if isinstance(labels_iteration, int):
labels_iteration = labels_iteration, labels_iteration
ketdata = ket.full()
dim_list_y = dim_list[:splitting]
dim_list_x = dim_list[splitting:]
size_x = np.prod(dim_list_x)
size_y = np.prod(dim_list_y)
ketdata = ketdata.reshape((size_y, size_x))
dim_list_small_x = dim_list_x[:labels_iteration[1]]
dim_list_small_y = dim_list_y[:labels_iteration[0]]
quadrants_x = np.prod(dim_list_small_x)
quadrants_y = np.prod(dim_list_small_y)
ticks_x = [size_x / quadrants_x * (i + 0.5)
for i in range(quadrants_x)]
ticks_y = [size_y / quadrants_y * (quadrants_y - i - 0.5)
for i in range(quadrants_y)]
labels_x = [_sequence_to_latex(_index_to_sequence(i*size_x // quadrants_x,
dim_list=dim_list_x))
for i in range(quadrants_x)]
labels_y = [_sequence_to_latex(_index_to_sequence(i*size_y // quadrants_y,
dim_list=dim_list_y))
for i in range(quadrants_y)]
ax.set_xticks(ticks_x)
ax.set_xticklabels(labels_x)
ax.set_yticks(ticks_y)
ax.set_yticklabels(labels_y)
ax.set_xlabel("last particles")
ax.set_ylabel("first particles")
ax.imshow(complex_array_to_rgb(ketdata, theme=theme),
interpolation="none",
extent=(0, size_x, 0, size_y))
return fig, ax
示例15: floquet_markov_mesolve
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 = Options()
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 = Result()
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)
#
# 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, True))
else:
# calculate all the expectation values, or output rho if
# no operators
if n_expt_op == 0:
if floquet_basis:
#.........这里部分代码省略.........