Improvedsampling

docs/src/jump_unraveling.md

@@ -48,6 +48,17 @@ BackAction.calculatechannelweights!(P::Vector{Float64}, psi::Vector{ComplexF64},
 ```@docs
 BackAction.calculatechannelweights!(P::Vector{Float64}, psi::Matrix{ComplexF64}, sys::System)
 ```
+
+```@docs
+BackAction.sampletauindex!(W::Vector{Float64}, Qs::Array{ComplexF64}, psi::Vector{ComplexF64},
+                         params::SimulParameters)
+```
+
+```@docs
+BackAction.sampletauindex!(W::Vector{Float64}, Qs::Array{ComplexF64}, psi::Matrix{ComplexF64},
+                         params::SimulParameters)
+```
+
 ### State Updates
 ```@docs
 BackAction.prejumpupdate!(V::Matrix{ComplexF64}, psi::Vector{ComplexF64}; normalize=false)
@@ -91,6 +102,12 @@ BackAction.run_singletrajectory_renewal(sys::System, params::SimulParameters, W:
 BackAction.gillipsiestep_returntau!(sys::System, params::SimulParameters, W::Vector{Float64},P::Vector{Float64}, Vs::Array{ComplexF64}, ts::Vector{Float64},t::Float64, psi::VecOrMat{ComplexF64}, traj::Trajectory )
 ```
 
+```@docs
+BackAction.gillipsiestep_returntau!(sys::System, params::SimulParameters, W::Vector{Float64},
+                        P::Vector{Float64}, Vs::Array{ComplexF64}, ts::Vector{Float64},
+                        t::Float64, psi::VecOrMat{ComplexF64}, traj::Trajectory, Qs::Array{ComplexF64}  )
+```
+
 ```@docs
 BackAction.writestate!(states::Array{ComplexF64}, psi::Union{Vector{ComplexF64}, Matrix{ComplexF64}}, counter::Int64)
 ```

src/functions_jump.jl

@@ -203,6 +203,56 @@ function postjumpupdate!(L::Matrix{ComplexF64}, psi::Matrix{ComplexF64}; normali
         end
 end
 
+"""
+
+```
+samplejumptime!(W::Vector{Float64}, Qs::Array{ComplexF64}, psi::VecOrMat{ComplexF64})
+```
+
+Sample a jump time index from the state `psi` (pure or mixed), modfying `W` to write on it.
+The technique is inversion sampling
+"""
+function sampletauindex!(W::Vector{Float64}, Qs::Array{ComplexF64}, psi::Vector{ComplexF64},
+                         params::SimulParameters)
+    # First, sample a random number and divide by dt to avoid multiplying by dt the weights
+    alpha = rand() / params.dt
+    u = 0.0
+    # now sum until alpha is exceded
+    # println(psi)
+    tau_index = 1
+    while u < alpha && tau_index < params.nsamples
+        u = u + real(dot(psi, Qs[:, :, tau_index], psi))
+        tau_index = tau_index + 1
+    end
+    return tau_index
+end
+
+"""
+
+```
+samplejumptime!(W::Vector{Float64}, Qs::Array{ComplexF64}, psi::VecOrMat{ComplexF64})
+```
+
+Sample a jump time index from the state `psi` (pure or mixed), modfying `W` to write on it.
+The technique is inversion sampling
+"""
+function sampletauindex!(W::Vector{Float64}, Qs::Array{ComplexF64}, psi::Matrix{ComplexF64},
+                         params::SimulParameters)
+    # First, sample a random number and divide by dt to avoid multiplying by dt the weights
+    alpha = rand() / params.dt
+    u = 0.0
+    # now sum until alpha is exceded
+    tau_index = 1
+    while u < alpha && tau_index < params.nsamples
+        u = u + real(tr(Qs[:, :, tau_index]* psi))
+        tau_index = tau_index + 1
+    end
+    return tau_index
+end
+
+
+
+
 """
 
 ```
@@ -212,13 +262,12 @@ gillipsiestep_returntau!(sys::System, params::SimulParameters, W::Vector{Float64
 
 ```
 Do a step of the Gillipsie algorithm, updating the state and the weights, and returning the
-obtained jump time.
+obtained jump time. In this version the time jump sampling is done by calling `StatsBase`.
 """
 function gillipsiestep_returntau!(sys::System, params::SimulParameters, W::Vector{Float64},
                         P::Vector{Float64}, Vs::Array{ComplexF64}, ts::Vector{Float64},
                         t::Float64, psi::VecOrMat{ComplexF64}, traj::Trajectory )
     #Sample jump time and  move state to pre-jump state
-    # println(psi)
     tau_index = StatsBase.sample(1:params.nsamples, StatsBase.weights(W))
     prejumpupdate!(Vs[:, :, tau_index], psi)
     # Sample jump channel
@@ -231,6 +280,44 @@ function gillipsiestep_returntau!(sys::System, params::SimulParameters, W::Vecto
     return tau
 
 end
+
+
+"""
+
+```
+gillipsiestep_returntau!(sys::System, params::SimulParameters, W::Vector{Float64},
+                        P::Vector{Float64}, Qs::Array{ComplexF64}, Vs::Array{ComplexF64},
+ ts::Vector{Float64},
+                        t::Float64, psi::VecOrMat{ComplexF64}, traj::Trajectory )
+
+```
+
+Do a step of the Gillipsie algorithm, updating the state and the weights, and returning the
+obtained jump time. In this version the time is extracted using inversion sampling instead of
+calling `StatsBase`.
+"""
+function gillipsiestep_returntau!(sys::System, params::SimulParameters, W::Vector{Float64},
+                        P::Vector{Float64}, Vs::Array{ComplexF64}, ts::Vector{Float64},
+                        t::Float64, psi::VecOrMat{ComplexF64}, traj::Trajectory, Qs::Array{ComplexF64}  )
+    tau_index = sampletauindex!(W, Qs, psi, params)
+    # in case the last index was at the last index, return already to avoid errors with dark states
+    if tau_index == params.nsamples
+        # push!(traj, DetectionClick(ts[tau_index], channel))
+        return tau_index
+    end
+    prejumpupdate!(Vs[:, :, tau_index], psi)
+    # Sample jump channel
+    calculatechannelweights!(P, psi, sys)
+    channel = StatsBase.sample(1:sys.NCHANNELS, StatsBase.weights(P))
+    # State update
+    postjumpupdate!(sys.Ls[channel], psi)
+    tau = ts[tau_index]
+    push!(traj, DetectionClick(tau, channel))
+    return tau
+
+end
+
+
 ############# Single Trajectory Routine ######################
 """
 ```
@@ -267,14 +354,9 @@ function run_singletrajectory(sys::System, params::SimulParameters,
     t::Float64 = 0
     channel = 0
     # Run the trajectory
-    calculatewtdweights!(W, Qs, psi, params)
+    # calculatewtdweights!(W, Qs, psi, params)
     while t < params.tf
-        t = t + gillipsiestep_returntau!(sys, params, W, P, Vs, ts, t, psi, traj)
-        # reSample WTD
-        calculatewtdweights!(W, Qs, psi, params)
-        if sum(W) < params.eps
-            break
-        end
+        t = t + gillipsiestep_returntau!(sys, params, W, P, Vs, ts, t, psi, traj, Qs)
     end
     return traj
 end

src/monitoring.jl

@@ -275,7 +275,7 @@ function monitoringoperator(t_given::Vector{Float64},
                 break
             end
         end
-        return states
+        return xis
     end
     # Evaluations before first jump
     while (t_given[counter] < traj[counter_c].time) && (counter <= ntimes)

test/runtests.jl

@@ -7,27 +7,29 @@ using Statistics
 using Random
 using Distributions
 import HypothesisTests
-@testset "Radiative Damping" begin
-    include("test_radiative_damping.jl")
-end
+@testset "Complete Test" begin
+    @testset "Radiative Damping" begin
+        include("test_radiative_damping.jl")
+    end
 
-@testset "states_atjumps" begin
-    include("test_states_atjumps.jl")
-end
+    @testset "states_atjumps" begin
+        include("test_states_atjumps.jl")
+    end
 
-@testset "states_att" begin
-    include("test_states_att.jl")
-end
+    @testset "states_att" begin
+        include("test_states_att.jl")
+    end
 
 
-@testset "Resonance Fluorescene" begin
-    include("test_resonance_fluorescene.jl")
-end
+    @testset "Resonance Fluorescene" begin
+        include("test_resonance_fluorescene.jl")
+    end
 
-@testset "Driven Qubit" begin
-    include("test_drive_qubit.jl")
-end
+    @testset "Driven Qubit" begin
+        include("test_drive_qubit.jl")
+    end
 
-@testset "Monitoring Operator" begin
-    include("test_monitoring.jl")
+    @testset "Monitoring Operator" begin
+        include("test_monitoring.jl")
+    end
 end

test/test_drive_qubit.jl

@@ -34,40 +34,38 @@ end
 # Trajectory Sampling
 sampled_trajectories = run_trajectories(sys, params);
 
-# Lindblad Evolution of Observables
-sigma = [BackAction.sigma_x, BackAction.sigma_y, BackAction.sigma_z]
-r0 = zeros(Float64, 3)
-for k in 1:3
-    r0[k] = real(tr(sigma[k]*psi0))
-end
-tspan = (0.0, params.tf)
-ntimes = 1000
-t_given = collect(LinRange(0, params.tf, ntimes));
+@testset "Driven Qubit: Expectation Value Convergence" begin
+    # Lindblad Evolution of Observables
+    sigma = [BackAction.sigma_x, BackAction.sigma_y, BackAction.sigma_z]
+    r0 = zeros(Float64, 3)
+    for k in 1:3
+        r0[k] = real(tr(sigma[k]*psi0))
+    end
+    tspan = (0.0, params.tf)
+    ntimes = 1000
+    t_given = collect(LinRange(0, params.tf, ntimes));
 
-# Analytical Solution
+    # Analytical Solution
 
-# Obtain the states between jumps
-sample = zeros(ComplexF64, sys.NLEVELS, sys.NLEVELS, ntimes, params.ntraj)
-for n in 1:params.ntraj
-    sample[:, :, :, n]  = BackAction.states_att(t_given, sampled_trajectories[n], sys,  params.psi0)
-end
-# Evaluate the observables
-r_sample = zeros(Float64, ntimes, 3, params.ntraj)
+    # Obtain the states between jumps
+    sample = zeros(ComplexF64, sys.NLEVELS, sys.NLEVELS, ntimes, params.ntraj)
+    for n in 1:params.ntraj
+        sample[:, :, :, n]  = BackAction.states_att(t_given, sampled_trajectories[n], sys,  params.psi0)
+    end
+    # Evaluate the observables
+    r_sample = zeros(Float64, ntimes, 3, params.ntraj)
 
-@time begin
-for j in 1:params.ntraj
-    for k in 1:3
-        for tn in 1:ntimes
+    for j in 1:params.ntraj
+        for k in 1:3
+            for tn in 1:ntimes
                 r_sample[tn, k, j] = real( tr(sigma[k] * sample[:, :, tn, j]) )
+            end
         end
     end
-end
-end
-# Average
-r_avg = dropdims(mean(r_sample, dims=3), dims=3);
+    # Average
+    r_avg = dropdims(mean(r_sample, dims=3), dims=3);
 
-r_analytical = BackAction.rk4(f_drivenqubit, r0, tspan, ntimes)
-@testset "Driven Qubit: Expectation Value Convergence" begin
+    r_analytical = BackAction.rk4(f_drivenqubit, r0, tspan, ntimes)
     for k in 1:ntimes
         @test abs( r_analytical[1, k]- r_avg[k, 1]) < 0.05
         @test abs( r_analytical[2, k]- r_avg[k, 2]) < 0.05