Skip to content

Commit

Permalink
Fix some style issue
Browse files Browse the repository at this point in the history
  • Loading branch information
Ericgig committed Mar 22, 2023
1 parent f328930 commit f0954e4
Show file tree
Hide file tree
Showing 3 changed files with 81 additions and 72 deletions.
34 changes: 21 additions & 13 deletions qutip/solver/sode/_noise.py
Original file line number Diff line number Diff line change
Expand Up @@ -2,10 +2,12 @@

__all__ = []


class _Noise:
"""
Weiner process generator used for tests.
"""

def __init__(self, T, dt, num=1):
N = int(np.round(T / dt))
self.T = T
Expand All @@ -17,19 +19,22 @@ def dw(self, dt):
"""
Ito integral I(i).
"""
N = int(np.round(dt /self.dt))
N = int(np.round(dt / self.dt))
return self.noise.reshape(-1, N, self.num).sum(axis=1)

def dz(self, dt):
"""
Ito integral I(0, i).
"""
N = int(np.round(dt /self.dt))
return np.einsum(
"ijk,j->ik",
self.noise.reshape(-1, N, self.num),
np.arange(N-0.5, 0, -1)
) * self.dt
N = int(np.round(dt / self.dt))
return (
np.einsum(
"ijk,j->ik",
self.noise.reshape(-1, N, self.num),
np.arange(N - 0.5, 0, -1),
)
* self.dt
)

def dW(self, dt):
"""
Expand All @@ -38,12 +43,15 @@ def dW(self, dt):
N = int(np.round(dt / self.dt))
noise = self.noise.copy()
if noise.shape[0] % N:
noise = noise[:-(noise.shape[0] % N)]
noise = noise[: -(noise.shape[0] % N)]
out = np.empty((noise.shape[0] // N, 2, self.num), dtype=float)
out[:, 0, :] = noise.reshape(-1, N, self.num).sum(axis=1)
out[:, 1, :] = np.einsum(
"ijk,j->ik",
self.noise.reshape(-1, N, self.num),
np.arange(N-0.5, 0, -1)
) * self.dt
out[:, 1, :] = (
np.einsum(
"ijk,j->ik",
self.noise.reshape(-1, N, self.num),
np.arange(N - 0.5, 0, -1),
)
* self.dt
)
return out
12 changes: 7 additions & 5 deletions qutip/solver/sode/sode.py
Original file line number Diff line number Diff line change
Expand Up @@ -44,6 +44,7 @@ class SIntegrator(Integrator):
keys, not the full options object passed to the solver. Options' keys
included here will be supported by the :cls:SolverOdeOptions.
"""

def set_state(self, t, state0, generator):
"""
Set the state of the SODE solver.
Expand Down Expand Up @@ -88,7 +89,6 @@ def integrate(self, t, copy=True):
"""
raise NotImplementedError


def mcstep(self, t, copy=True):
raise NotImplementedError

Expand All @@ -97,6 +97,7 @@ class _Explicit_Simple_Integrator(SIntegrator):
"""
Stochastic evolution solver
"""

integrator_options = {
"dt": 0.001,
"tol": 1e-10,
Expand All @@ -111,7 +112,7 @@ def __init__(self, rhs, options):
self.step_func = self.stepper(self.system).run

def integrate(self, t, copy=True):
delta_t = (t - self.t)
delta_t = t - self.t
if delta_t < 0:
raise ValueError("Stochastic integration time")
elif delta_t == 0:
Expand All @@ -125,9 +126,7 @@ def integrate(self, t, copy=True):
N += 1
dt = delta_t / N
dW = self.generator.normal(
0,
np.sqrt(dt),
size=(N, self.N_dw, self.system.num_collapse)
0, np.sqrt(dt), size=(N, self.N_dw, self.system.num_collapse)
)

self.state = self.step_func(self.t, self.state, dt, dW, N)
Expand Down Expand Up @@ -157,6 +156,7 @@ class _Implicit_Simple_Integrator(_Explicit_Simple_Integrator):
"""
Stochastic evolution solver
"""

integrator_options = {
"dt": 0.001,
"tol": 1e-10,
Expand Down Expand Up @@ -213,6 +213,7 @@ class PlatenSODE(_Explicit_Simple_Integrator):
- Order: strong 1, weak 2
"""

stepper = _sode.Platen
N_dw = 1

Expand All @@ -232,6 +233,7 @@ class PredCorr_SODE(_Explicit_Simple_Integrator):
(:math:`\\alpha=1/2`, :math:`\\eta=1/2`): ``'pc-euler-imp'``,
``'pc-euler-2'`` or ``'pred-corr-2'``
"""

integrator_options = {
"dt": 0.001,
"tol": 1e-10,
Expand Down
107 changes: 53 additions & 54 deletions qutip/solver/stochastic.py
Original file line number Diff line number Diff line change
Expand Up @@ -3,7 +3,7 @@
from .sode.ssystem import *
from .result import MultiTrajResult, Result, ExpectOp
from .multitraj import MultiTrajSolver
from ..import Qobj, QobjEvo, liouvillian, lindblad_dissipator
from .. import Qobj, QobjEvo, liouvillian, lindblad_dissipator
import numpy as np
from collections.abc import Iterable
from functools import partial
Expand All @@ -19,17 +19,14 @@ def _post_init(self, m_ops=(), dw_factor=(), heterodyne=False):
self.heterodyne = heterodyne
for op in m_ops:
f = self._e_op_func(op)
self.W.append([0.])
self.W.append([0.0])
self.m_expect.append([])
self.m_ops.append(ExpectOp(op, f, self.m_expect[-1].append))
self.add_processor(self.m_ops[-1]._store)

def add(self, t, state, noise):
super().add(t, state)
if (
noise is not None
and self.options['store_measurement']
):
if noise is not None and self.options["store_measurement"]:
for i, dW in enumerate(noise):
self.W[i].append(self.W[i][-1] + dW)

Expand Down Expand Up @@ -78,7 +75,9 @@ def measurement(self):
"""
dts = np.diff(self.times)
m_expect = np.array(self.m_expect)[:, 1:]
noise = np.einsum("i,ij,j->ij", self.dW_factor, np.diff(self.W, axis=1), (1/dts))
noise = np.einsum(
"i,ij,j->ij", self.dW_factor, np.diff(self.W, axis=1), (1 / dts)
)
if self.heterodyne:
m_expect = m_expect.reshape(-1, 2, m_expect.shape[1])
noise = noise.reshape(-1, 2, noise.shape[1])
Expand All @@ -89,11 +88,13 @@ class StochasticResult(MultiTrajResult):
def _post_init(self):
super()._post_init()

store_measurement = self.options['store_measurement']
keep_runs = self.options['keep_runs_results']
store_measurement = self.options["store_measurement"]
keep_runs = self.options["keep_runs_results"]

if not keep_runs and store_measurement:
self.add_processor(partial(self._reduce_attr, attr="wiener_process"))
self.add_processor(
partial(self._reduce_attr, attr="wiener_process")
)
self._wiener_process = []
self.add_processor(partial(self._reduce_attr, attr="dW"))
self._dW = []
Expand All @@ -115,10 +116,8 @@ def _trajectories_attr(self, attr):
"""
if hasattr(self, "_" + attr):
return getattr(self, "_" + attr)
elif self.options['keep_runs_results']:
return np.array([
getattr(traj, attr) for traj in self.trajectories
])
elif self.options["keep_runs_results"]:
return np.array([getattr(traj, attr) for traj in self.trajectories])
return None

@property
Expand Down Expand Up @@ -174,6 +173,7 @@ class _StochasticRHS:
the rouchon integrator need the part but does not use the usual drift and
diffusion computation.
"""

def __init__(self, issuper, H, sc_ops, c_ops, heterodyne):

if not isinstance(H, (Qobj, QobjEvo)) or not H.isoper:
Expand Down Expand Up @@ -212,15 +212,17 @@ def __init__(self, issuper, H, sc_ops, c_ops, heterodyne):
def __call__(self, options):
if self.issuper:
return StochasticOpenSystem(
self.H , self.sc_ops, self.c_ops, options.get("dt", 1.-6)
self.H, self.sc_ops, self.c_ops, options.get("dt", 1.0 - 6)
)
else:
return StochasticClosedSystem(self.H , self.sc_ops)
return StochasticClosedSystem(self.H, self.sc_ops)


def smesolve(H, rho0, tlist, c_ops=(), sc_ops=(), e_ops=(), m_ops=(),
args={}, ntraj=500, options=None, heterodyne=False,
seeds=None, target_tol=None, timeout=None):
def smesolve(
H, rho0, tlist, c_ops=(), sc_ops=(), heterodyne=False, /,
e_ops=(), args={}, ntraj=500, options=None,
seeds=None, target_tol=None, timeout=None,
):
"""
Solve stochastic master equation.
Expand Down Expand Up @@ -250,11 +252,6 @@ def smesolve(H, rho0, tlist, c_ops=(), sc_ops=(), e_ops=(), m_ops=(),
Callable signature must be, `f(t: float, state: Qobj)`.
See :func:`expect` for more detail of operator expectation.
m_ops : list of :class:`qutip.Qobj`
List of measurements operators, when not provided, they are computed
from the sc_ops according to whether homodyne or heterodyne detection
is used.
args : None / *dictionary*
Dictionary of parameters for time-dependent Hamiltonians and
collapse operators.
Expand Down Expand Up @@ -337,17 +334,20 @@ def smesolve(H, rho0, tlist, c_ops=(), sc_ops=(), e_ops=(), m_ops=(),
H = QobjEvo(H, args=args, tlist=tlist)
c_ops = [QobjEvo(c_op, args=args, tlist=tlist) for c_op in c_ops]
sc_ops = [QobjEvo(c_op, args=args, tlist=tlist) for c_op in sc_ops]
sol = SMESolver(H, sc_ops, c_ops=c_ops,
options=options, heterodyne=heterodyne)
if m_ops:
sol.m_ops = m_ops
return sol.run(rho0, tlist, ntraj, e_ops=e_ops,
seed=seeds, target_tol=target_tol, timeout=timeout)


def ssesolve(H, psi0, tlist, sc_ops=(), e_ops=(), m_ops=(),
args={}, ntraj=500, options=None, heterodyne=False,
seeds=None, target_tol=None, timeout=None):
sol = SMESolver(
H, sc_ops, c_ops=c_ops, options=options, heterodyne=heterodyne
)
return sol.run(
rho0, tlist, ntraj, e_ops=e_ops,
seed=seeds,target_tol=target_tol, timeout=timeout,
)


def ssesolve(
H, psi0, tlist, sc_ops=(), heterodyne=False, /,
e_ops=(), args={}, ntraj=500, options=None,
seeds=None, target_tol=None, timeout=None,
):
"""
Solve stochastic Schrodinger equation.
Expand All @@ -373,11 +373,6 @@ def ssesolve(H, psi0, tlist, sc_ops=(), e_ops=(), m_ops=(),
Callable signature must be, `f(t: float, state: Qobj)`.
See :func:`expect` for more detail of operator expectation.
m_ops : list of :class:`qutip.Qobj`
List of measurements operators, when not provided, they are computed
from the sc_ops according to whether homodyne or heterodyne detection
is used.
args : None / *dictionary*
Dictionary of parameters for time-dependent Hamiltonians and
collapse operators.
Expand Down Expand Up @@ -458,16 +453,17 @@ def ssesolve(H, psi0, tlist, sc_ops=(), e_ops=(), m_ops=(),
H = QobjEvo(H, args=args, tlist=tlist)
sc_ops = [QobjEvo(c_op, args=args, tlist=tlist) for c_op in sc_ops]
sol = SSESolver(H, sc_ops, options=options, heterodyne=heterodyne)
if m_ops:
sol.m_ops = m_ops
return sol.run(psi0, tlist, ntraj, e_ops=e_ops,
seed=seeds, target_tol=target_tol, timeout=timeout)
return sol.run(
psi0, tlist, ntraj, e_ops=e_ops,
seed=seeds, target_tol=target_tol, timeout=timeout,
)


class StochasticSolver(MultiTrajSolver):
"""
Generic stochastic solver.
"""

name = "StochasticSolver"
resultclass = StochasticResult
_avail_integrators = {}
Expand Down Expand Up @@ -499,9 +495,7 @@ def __init__(self, H, sc_ops, heterodyne, *, c_ops=(), options=None):
if heterodyne:
self._m_ops = []
for op in sc_ops:
self._m_ops += [
op + op.dag(), -1j * (op - op.dag())
]
self._m_ops += [op + op.dag(), -1j * (op - op.dag())]
self._dW_factors = np.ones(len(sc_ops) * 2) * 0.5**0.5
else:
self._m_ops = [op + op.dag() for op in sc_ops]
Expand Down Expand Up @@ -556,8 +550,10 @@ def m_ops(self, new_m_ops):
isinstance(op, Qobj) and op.dims == self.rhs.sc_ops[0].dims
for op in new_m_ops
):
raise ValueError("m_ops must be Qobj with the same dimensions"
" as the Hamiltonian")
raise ValueError(
"m_ops must be Qobj with the same dimensions"
" as the Hamiltonian"
)
self._m_ops = new_m_ops

@property
Expand Down Expand Up @@ -588,8 +584,11 @@ def _run_one_traj(self, seed, state, tlist, e_ops):
Run one trajectory and return the result.
"""
result = StochasticTrajResult(
e_ops, self.options, m_ops=self.m_ops, dw_factor=self.dW_factors,
heterodyne=self.heterodyne
e_ops,
self.options,
m_ops=self.m_ops,
dw_factor=self.dW_factors,
heterodyne=self.heterodyne,
)
generator = self._get_generator(seed)
self._integrator.set_state(tlist[0], state, generator)
Expand Down Expand Up @@ -630,9 +629,9 @@ def options(self):
brownian noise for each trajectories.
progress_bar: str {'text', 'enhanced', 'tqdm', ''}, default="text"
How to present the solver progress.
'tqdm' uses the python module of the same name and raise an error if
not installed. Empty string or False will disable the bar.
How to present the solver progress. 'tqdm' uses the python module
of the same name and raise an error if not installed. Empty string
or False will disable the bar.
progress_kwargs: dict, default={"chunk_size":10}
Arguments to pass to the progress_bar. Qutip's bars use
Expand Down

0 comments on commit f0954e4

Please sign in to comment.