Source code for qutip.fortran.mcsolve_f90

# This file is part of QuTiP: Quantum Toolbox in Python.
#
#    Copyright (c) 2011 and later, Paul D. Nation and Robert J. Johansson.
#    All rights reserved.
#
#    Redistribution and use in source and binary forms, with or without
#    modification, are permitted provided that the following conditions are
#    met:
#
#    1. Redistributions of source code must retain the above copyright notice,
#       this list of conditions and the following disclaimer.
#
#    2. Redistributions in binary form must reproduce the above copyright
#       notice, this list of conditions and the following disclaimer in the
#       documentation and/or other materials provided with the distribution.
#
#    3. Neither the name of the QuTiP: Quantum Toolbox in Python nor the names
#       of its contributors may be used to endorse or promote products derived
#       from this software without specific prior written permission.
#
#    THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
#    "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
#    LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A
#    PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
#    HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
#    SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
#    LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
#    DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
#    THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
#    (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
#    OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
#
#    Original Code by Arne Grimsmo (2012): github.com/arnelg/qutipf90mc
###############################################################################

import numpy as np

from qutip.fortran import qutraj_run as qtf90
from qutip.qobj import Qobj
from qutip.mcsolve import _mc_data_config
from qutip.solver import Options, Result, config
from qutip.settings import debug
import qutip.settings

if debug:
    import inspect
    import os

# Working precision
wpr = np.dtype(np.float64)
wpc = np.dtype(np.complex128)


[docs]def mcsolve_f90(H, psi0, tlist, c_ops, e_ops, ntraj=None, options=Options(), sparse_dms=True, serial=False, ptrace_sel=[], calc_entropy=False): """ Monte-Carlo wave function solver with fortran 90 backend. Usage is identical to qutip.mcsolve, for problems without explicit time-dependence, and with some optional input: Parameters ---------- H : qobj System Hamiltonian. psi0 : qobj Initial state vector tlist : array_like Times at which results are recorded. ntraj : int Number of trajectories to run. c_ops : array_like ``list`` or ``array`` of collapse operators. e_ops : array_like ``list`` or ``array`` of operators for calculating expectation values. options : Options Instance of solver options. sparse_dms : boolean If averaged density matrices are returned, they will be stored as sparse (Compressed Row Format) matrices during computation if sparse_dms = True (default), and dense matrices otherwise. Dense matrices might be preferable for smaller systems. serial : boolean If True (default is False) the solver will not make use of the multiprocessing module, and simply run in serial. ptrace_sel: list This optional argument specifies a list of components to keep when returning a partially traced density matrix. This can be convenient for large systems where memory becomes a problem, but you are only interested in parts of the density matrix. calc_entropy : boolean If ptrace_sel is specified, calc_entropy=True will have the solver return the averaged entropy over trajectories in results.entropy. This can be interpreted as a measure of entanglement. See Phys. Rev. Lett. 93, 120408 (2004), Phys. Rev. A 86, 022310 (2012). Returns ------- results : Result Object storing all results from simulation. """ if ntraj is None: ntraj = options.ntraj if psi0.type != 'ket': raise Exception("Initial state must be a state vector.") config.options = options # set num_cpus to the value given in qutip.settings # if none in Options if not config.options.num_cpus: config.options.num_cpus = qutip.settings.num_cpus # set initial value data if options.tidy: config.psi0 = psi0.tidyup(options.atol).full() else: config.psi0 = psi0.full() config.psi0_dims = psi0.dims config.psi0_shape = psi0.shape # set general items config.tlist = tlist if isinstance(ntraj, (list, np.ndarray)): raise Exception("ntraj as list argument is not supported.") else: config.ntraj = ntraj # ntraj_list = [ntraj] # set norm finding constants config.norm_tol = options.norm_tol config.norm_steps = options.norm_steps if not options.rhs_reuse: config.soft_reset() # no time dependence config.tflag = 0 # check for collapse operators if len(c_ops) > 0: config.cflag = 1 else: config.cflag = 0 # Configure data _mc_data_config(H, psi0, [], c_ops, [], [], e_ops, options, config) # 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 Result object output = Result() # 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 = config.tlist output.num_expect = config.e_num output.num_collapse = config.c_num output.ntraj = config.ntraj return output
class _MC_class(): def __init__(self): self.cpus = config.options.num_cpus self.nprocs = self.cpus self.sol = Result() self.mf = 10 # If returning density matrices, return as sparse or dense? self.sparse_dms = True # Run in serial? self.serial_run = False self.ntraj = config.ntraj self.ntrajs = [] self.seed = None self.psi0_dims = None self.psi0_shape = None self.dm_dims = None self.dm_shape = None self.unravel_type = 2 self.ptrace_sel = [] self.calc_entropy = False def parallel(self): from multiprocessing import Process, Queue, JoinableQueue if debug: print(inspect.stack()[0][3]) self.ntrajs = [] for i in range(self.cpus): self.ntrajs.append(min(int(np.floor(float(self.ntraj) / self.cpus)), self.ntraj - sum(self.ntrajs))) cnt = sum(self.ntrajs) while cnt < self.ntraj: for i in range(self.cpus): self.ntrajs[i] += 1 cnt += 1 if (cnt >= self.ntraj): break self.ntrajs = np.array(self.ntrajs) self.ntrajs = self.ntrajs[np.where(self.ntrajs > 0)] self.nprocs = len(self.ntrajs) sols = [] processes = [] resq = JoinableQueue() resq.join() if debug: print("Number of cpus: " + str(self.cpus)) print("Trying to start " + str(self.nprocs) + " process(es).") print("Number of trajectories for each process: " + str(self.ntrajs)) for i in range(self.nprocs): p = Process(target=self.evolve_serial, args=((resq, self.ntrajs[i], i, self.seed * (i + 1)),)) p.start() processes.append(p) cnt = 0 while True: try: sols.append(resq.get()) resq.task_done() cnt += 1 if (cnt >= self.nprocs): break except KeyboardInterrupt: break except: pass resq.join() for proc in processes: try: proc.join() except KeyboardInterrupt: if debug: print("Cancel thread on keyboard interrupt") proc.terminate() proc.join() resq.close() return sols def serial(self): if debug: print(inspect.stack()[0][3]) self.nprocs = 1 self.ntrajs = [self.ntraj] if debug: print("Running in serial.") print("Number of trajectories: " + str(self.ntraj)) sol = self.evolve_serial((0, self.ntraj, 0, self.seed)) return [sol] def run(self): if debug: print(inspect.stack()[0][3]) from numpy.random import random_integers if (config.c_num == 0): # force one trajectory if no collapse operators config.ntraj = 1 self.ntraj = 1 # Set unravel_type to 1 to integrate without collapses self.unravel_type = 1 if (config.e_num == 0): # If we are returning states, and there are no # collapse operators, set average_states to False to return # ket vectors instead of density matrices config.options.average_states = False # generate a random seed, useful if e.g. running with MPI self.seed = random_integers(1e8) if (self.serial_run): # run in serial sols = self.serial() else: # run in paralell sols = self.parallel() # gather data self.sol = _gather(sols) def evolve_serial(self, args): if debug: print(inspect.stack()[0][3] + ":" + str(os.getpid())) # run ntraj trajectories for one process via fortran # get args queue, ntraj, instanceno, rngseed = args # initialize the problem in fortran _init_tlist() _init_psi0() if (self.ptrace_sel != []): _init_ptrace_stuff(self.ptrace_sel) _init_hamilt() if (config.c_num != 0): _init_c_ops() if (config.e_num != 0): _init_e_ops() # set options qtf90.qutraj_run.n_c_ops = config.c_num qtf90.qutraj_run.n_e_ops = config.e_num qtf90.qutraj_run.ntraj = ntraj qtf90.qutraj_run.unravel_type = self.unravel_type qtf90.qutraj_run.average_states = config.options.average_states qtf90.qutraj_run.average_expect = config.options.average_expect qtf90.qutraj_run.init_result(config.psi0_shape[0], config.options.atol, config.options.rtol, mf=self.mf, norm_steps=config.norm_steps, norm_tol=config.norm_tol) # set optional arguments qtf90.qutraj_run.order = config.options.order qtf90.qutraj_run.nsteps = config.options.nsteps qtf90.qutraj_run.first_step = config.options.first_step qtf90.qutraj_run.min_step = config.options.min_step qtf90.qutraj_run.max_step = config.options.max_step qtf90.qutraj_run.norm_steps = config.options.norm_steps qtf90.qutraj_run.norm_tol = config.options.norm_tol # use sparse density matrices during computation? qtf90.qutraj_run.rho_return_sparse = self.sparse_dms # calculate entropy of reduced density matrice? qtf90.qutraj_run.calc_entropy = self.calc_entropy # run show_progress = 1 if debug else 0 qtf90.qutraj_run.evolve(instanceno, rngseed, show_progress) # construct Result instance sol = Result() sol.ntraj = ntraj # sol.col_times = qtf90.qutraj_run.col_times # sol.col_which = qtf90.qutraj_run.col_which-1 sol.col_times, sol.col_which = self.get_collapses(ntraj) if (config.e_num == 0): sol.states = self.get_states(len(config.tlist), ntraj) else: sol.expect = self.get_expect(len(config.tlist), ntraj) if (self.calc_entropy): sol.entropy = self.get_entropy(len(config.tlist)) if (not self.serial_run): # put to queue queue.put(sol) queue.join() # deallocate stuff # finalize() return sol # Routines for retrieving data data from fortran def get_collapses(self, ntraj): if debug: print(inspect.stack()[0][3]) col_times = np.zeros((ntraj), dtype=np.ndarray) col_which = np.zeros((ntraj), dtype=np.ndarray) if (config.c_num == 0): # no collapses return col_times, col_which for i in range(ntraj): qtf90.qutraj_run.get_collapses(i + 1) times = qtf90.qutraj_run.col_times which = qtf90.qutraj_run.col_which if times is None: times = np.array([]) if which is None: which = np.array([]) else: which = which - 1 col_times[i] = np.array(times, copy=True) col_which[i] = np.array(which, copy=True) return col_times, col_which def get_states(self, nstep, ntraj): if debug: print(inspect.stack()[0][3]) from scipy.sparse import csr_matrix if (config.options.average_states): states = np.array([Qobj()] * nstep) if (self.sparse_dms): # averaged sparse density matrices for i in range(nstep): qtf90.qutraj_run.get_rho_sparse(i + 1) val = qtf90.qutraj_run.csr_val col = qtf90.qutraj_run.csr_col - 1 ptr = qtf90.qutraj_run.csr_ptr - 1 m = qtf90.qutraj_run.csr_nrows k = qtf90.qutraj_run.csr_ncols states[i] = Qobj(csr_matrix((val, col, ptr), (m, k)).toarray(), dims=self.dm_dims, shape=self.dm_shape) else: # averaged dense density matrices for i in range(nstep): states[i] = Qobj(qtf90.qutraj_run.sol[0, i, :, :], dims=self.dm_dims, shape=self.dm_shape) else: # all trajectories as kets if (ntraj == 1): states = np.array([Qobj()] * nstep, dtype=object) for i in range(nstep): states[i] = Qobj(np.matrix( qtf90.qutraj_run.sol[0, 0, i, :]).transpose(), dims=self.psi0_dims, shape=self.psi0_shape) else: states = np.array([np.array([Qobj()] * nstep, dtype=object)] * ntraj) for traj in range(ntraj): for i in range(nstep): states[traj][i] = Qobj(np.matrix( qtf90.qutraj_run.sol[0, traj, i, :]).transpose(), dims=self.psi0_dims, shape=self.psi0_shape) return states def get_expect(self, nstep, ntraj): if debug: print(inspect.stack()[0][3]) if (config.options.average_expect): expect = [] for j in range(config.e_num): if config.e_ops_isherm[j]: expect += [np.real(qtf90.qutraj_run.sol[j, 0, :, 0])] else: expect += [qtf90.qutraj_run.sol[j, 0, :, 0]] else: expect = np.array([[np.array([0. + 0.j] * nstep)] * config.e_num] * ntraj) for j in range(config.e_num): expect[:, j, :] = qtf90.qutraj_run.sol[j, :, :, 0] return expect def get_entropy(self, nstep): if debug: print(inspect.stack()[0][3]) if (not self.calc_entropy): raise Exception('get_entropy: calc_entropy=False. Aborting.') entropy = np.array([0.] * nstep) entropy[:] = qtf90.qutraj_run.reduced_state_entropy[:] return entropy def finalize(): # not in use... if debug: print(inspect.stack()[0][3]) qtf90.qutraj_run.finalize_work() qtf90.qutraj_run.finalize_sol() def _gather(sols): # gather list of Result objects, sols, into one. sol = Result() # sol = sols[0] ntraj = sum([a.ntraj for a in sols]) sol.col_times = np.zeros((ntraj), dtype=np.ndarray) sol.col_which = np.zeros((ntraj), dtype=np.ndarray) sol.col_times[0:sols[0].ntraj] = sols[0].col_times sol.col_which[0:sols[0].ntraj] = sols[0].col_which sol.states = np.array(sols[0].states, dtype=object) sol.expect = np.array(sols[0].expect) if (hasattr(sols[0], 'entropy')): sol.entropy = np.array(sols[0].entropy) sofar = 0 for j in range(1, len(sols)): sofar = sofar + sols[j - 1].ntraj sol.col_times[sofar:sofar + sols[j].ntraj] = ( sols[j].col_times) sol.col_which[sofar:sofar + sols[j].ntraj] = ( sols[j].col_which) if (config.e_num == 0): if (config.options.average_states): # collect states, averaged over trajectories sol.states += np.array(sols[j].states) else: # collect states, all trajectories sol.states = np.vstack((sol.states, np.array(sols[j].states))) else: if (config.options.average_expect): # collect expectation values, averaged for i in range(config.e_num): sol.expect[i] += np.array(sols[j].expect[i]) else: # collect expectation values, all trajectories sol.expect = np.vstack((sol.expect, np.array(sols[j].expect))) if (hasattr(sols[j], 'entropy')): if (config.options.average_states or config.options.average_expect): # collect entropy values, averaged sol.entropy += np.array(sols[j].entropy) else: # collect entropy values, all trajectories sol.entropy = np.vstack((sol.entropy, np.array(sols[j].entropy))) if (config.options.average_states or config.options.average_expect): if (config.e_num == 0): sol.states = sol.states / len(sols) else: sol.expect = list(sol.expect / len(sols)) inds = np.where(config.e_ops_isherm)[0] for jj in inds: sol.expect[jj] = np.real(sol.expect[jj]) if (hasattr(sols[0], 'entropy')): sol.entropy = sol.entropy / len(sols) # convert sol.expect array to list and fix dtypes of arrays if (not config.options.average_expect) and config.e_num != 0: temp = [list(sol.expect[ii]) for ii in range(ntraj)] for ii in range(ntraj): for jj in np.where(config.e_ops_isherm)[0]: temp[ii][jj] = np.real(temp[ii][jj]) sol.expect = temp # convert to list/array to be consistent with qutip mcsolve sol.states = list(sol.states) return sol # # Functions to initialize the problem in fortran # def _init_tlist(): Of = _realarray_to_fortran(config.tlist) qtf90.qutraj_run.init_tlist(Of, np.size(Of)) def _init_psi0(): # Of = _qobj_to_fortranfull(config.psi0) Of = _complexarray_to_fortran(config.psi0) qtf90.qutraj_run.init_psi0(Of, np.size(Of)) def _init_ptrace_stuff(sel): psi0 = Qobj(config.psi0, dims=config.psi0_dims, shape=config.psi0_shape) qtf90.qutraj_run.init_ptrace_stuff(config.psi0_dims[0], np.array(sel) + 1, psi0.ptrace(sel).shape[0]) def _init_hamilt(): # construct effective non-Hermitian Hamiltonian # H_eff = H - 0.5j*sum([c_ops[i].dag()*c_ops[i] # for i in range(len(c_ops))]) # Of = _qobj_to_fortrancsr(H_eff) # qtf90.qutraj_run.init_hamiltonian(Of[0],Of[1], # Of[2],Of[3],Of[4]) d = np.size(config.psi0) qtf90.qutraj_run.init_hamiltonian( _complexarray_to_fortran(config.h_data), config.h_ind + 1, config.h_ptr + 1, d, d) def _init_c_ops(): d = np.size(config.psi0) n = config.c_num first = True for i in range(n): # Of = _qobj_to_fortrancsr(c_ops[i]) # qtf90.qutraj_run.init_c_ops(i+1,n,Of[0],Of[1], # Of[2],Of[3],Of[4],first) qtf90.qutraj_run.init_c_ops( i + 1, n, _complexarray_to_fortran(config.c_ops_data[i]), config.c_ops_ind[i] + 1, config.c_ops_ptr[i] + 1, d, d, first) first = False def _init_e_ops(): d = np.size(config.psi0) # n = config.e_num n = len(config.e_ops_data) first = True for i in range(n): # Of = _qobj_to_fortrancsr(e_ops[i]) # qtf90.qutraj_run.init_e_ops(i+1,n,Of[0],Of[1], # Of[2],Of[3],Of[4],first) qtf90.qutraj_run.init_e_ops( i + 1, n, _complexarray_to_fortran(config.e_ops_data[i]), config.e_ops_ind[i] + 1, config.e_ops_ptr[i] + 1, d, d, first) first = False # # Misc. converison functions # def _realarray_to_fortran(a): datad = np.array(a, dtype=wpr) return datad def _complexarray_to_fortran(a): datad = np.array(a, dtype=wpc) return datad def _qobj_to_fortranfull(A): datad = np.array(A.data.toarray(), dtype=wpc) return datad def _qobj_to_fortrancsr(A): data = np.array(A.data.data, dtype=wpc) indices = np.array(A.data.indices) indptr = np.array(A.data.indptr) m = A.data.shape[0] k = A.data.shape[1] return data, indices + 1, indptr + 1, m, k