Iterators for QuMCWFEnsemble

iterators_mcwf.jl

immutable QuMCWF <: QuPropagatorMethod
end

type QuMCWFEnsemble{QA<:QuBase.AbstractQuArray}
    ntraj
    rho::QA
    decomp
    QuMCWFEnsemble(ntraj, rho, decomp) = new(1:1:ntraj, rho, decomp)
end

QuMCWFEnsemble{QA<:QuBase.AbstractQuArray}(rho::QA, decomp) = QuMCWFEnsemble{QA}(500, rho, decomp)
QuMCWFEnsemble{QM<:QuBase.AbstractQuMatrix}(rho::QM) = QuMCWFEnsemble{QM}(500, rho, eigfact(reshape(full(coeffs(rho)), size(rho))))
QuMCWFEnsemble{QV<:QuBase.AbstractQuVector}(psi::QV) =  QuMCWFEnsemble{QV}(500, psi, nothing)

type QuMCWFEnsembleState
    rho
    traj
    traj_state
end

function Base.start(mcwfensemble::QuMCWFEnsemble)
    init_state = mcwfensemble.rho
    traj_state = start(mcwfensemble.ntraj)
    traj, traj_state = next(mcwfensemble.ntraj, traj_state)
    return QuMCWFEnsembleState(init_state, traj, traj_state)
end

function Base.next(mcwfensemble::QuMCWFEnsemble, state::QuMCWFEnsembleState)
    traj, traj_state = next(mcwfensemble.ntraj, state.traj_state)
    next_state = draw(mcwfensemble)
    return next_state, QuMCWFEnsembleState(next_state, traj, traj_state)
end

Base.done(mcwfensemble::QuMCWFEnsemble, state::QuMCWFEnsembleState) = done(mcwfensemble.ntraj, state.traj_state)

function draw{QM<:QuBase.AbstractQuMatrix}(mcwfensemble::QuMCWFEnsemble{QM})
    r = rand() # draw random number from [0,1)
    pacc = 0.
    for i=1:length(mcwfensemble.decomp[:values])
        pacc = pacc + mcwfensemble.decomp[:values][i]
        if pacc >= r
            return QuArray(complex(mcwfensemble.decomp[:vectors][:,i]))
        end
    end
end

function draw{QV<:QuBase.AbstractQuVector}(mcwfensemble::QuMCWFEnsemble{QV})
    return copy(mcwfensemble.rho)
end

function eff_hamiltonian(lme::QuLindbladMasterEq)
    # nlop = length(lme.collapse_ops)
    heff = lme.hamiltonian
    for c_ops in lme.collapse_ops
       heff = heff - im*0.5*c_ops'*c_ops
    end
    return heff
end

function mcwfpropagate(prob::QuMCWF, eq::QuLindbladMasterEq, ensemble::QuMCWFEnsembleState, t, current_t, current_qustate)
    jtol = 1.e-6  # jump tolerance
    # get information of QME
    heff = eff_hamiltonian(eq)
    c_ops = eq.collapse_ops
    # initial state -> ensemble
    # qse = QuMCWFEnsemble(current_qustate)
    CQST = QuBase.similar_type(current_qustate)
    next_state = CQST(zeros(size(current_qustate)), bases(current_qustate))
    # for traj = 1:qse.ntraj
        # initial state vec for trajectory
	      # psi = next(qse, start(qse)
        psi = ensemble.rho
        eps = rand() # draw random number from [0,1)
        dt = t - current_t
	      accdt = 0.
        tj = 0.
        while accdt < dt
		        # propagate one time-step
            psi1 = expm(-im*heff*dt)*psi
    		    if abs(norm(psi1) - eps) < jtol
    			      PnS = 0.
                for cind in c_ops
    				        PnS = PnS + real(expectation(psi1, cind'*cind))
    		        end
                Pn = 0.
                for cind in c_ops
    				        Pn = Pn + real(expectation(psi1, cind'*cind))/PnS
    				        if Pn >= eps
                        psi = psi1*cind
    					          break
    				        end
    			      end
    			      normalize!(psi)
    			      eps = rand() # draw new random number from [0,1)
    			      accdt = accdt + dt
    			      dt = dt - accdt
    		    elseif norm(psi1) < eps
    			      dt = dt/2		# jump in the current interval
    		    else
    			      psi = copy(psi1) # no jump
                normalize!(psi)
    			      accdt = accdt + dt
    			      dt = dt - accdt
    	      end
        end
        println(psi)
        next_state = next_state + CQST(reshape(coeffs(psi*psi'), size(current_qustate)), bases(current_qustate))/ensemble.traj
    # end
    return next_state
end

function propagate(prob::QuMCWF, eq::QuLindbladMasterEq, t, current_t, current_qustate)
    ens = QuMCWFEnsemble(current_qustate)
    return mcwfpropagate(prob, eq, ens, t, current_t, current_qustate)
end

export QuMCWFEnsemble,
      QuMCWF