Blink issue #426

00_README.md

This is some demo Julia code with nontrivial Unicode, used to illustrate Blink issue #426

Bases.jl

__precompile__()
module Bases
export Basis, length, nfactors

import Base: length, getindex
using IterTools: product

struct Basis{N}
    e::Vector{NTuple{N, Int}}
    dims::NTuple{N, Int}
end

length(bas::Basis) = length(bas.e)
nfactors{N}(bas::Basis{N}) = N
getindex(bas::Basis, ind::Int) = bas.e[ind]

function Basis(dims::Vararg{Int})
    e = (collect ∘ product)([0:d - 1 for d in reverse(dims)]...)
    return Basis(reverse.(e), dims)
end

end #module

Cavity.jl

module Cavity

export cavity

using ..Bases: Basis, length
using ..QOperators
using ..Simulations

function cavity(bas::Basis, Δ::Real, κ::Real, αin::Function, η::Real)
    a = MatrixOperator(destroy, bas, 1)
    a⁺ = a'
    id = eye(length(bas))

    function Heff(t, θ)
        α = θ[1]
        return (Δ * (a⁺ + conj(α) * id) * (a + α * id)
                + im * sqrt(κ) * (αin(t) * a⁺ + conj(αin(t)) * a))
    end
    Leff = [sqrt(κ) * a]
    Fs = [im * a⁺, -im * a]
    function θeom(t, θ)
        α = θ[1]
        αdot = -κ * α / 2
        return [αdot, conj(αdot)]
    end
    Ys = [-a, -a⁺]  # THIS IS THE CORRECT ORDER!!
    return MovingBasis(Heff, Fs, Ys, Leff, θeom, η, 2)
end

end  # module

Optomechanics.jl

__precompile__()
module Optomechanics

export single_node, two_nodes, eq_cavity_amplitude, eq_mech_displacement

using ..Bases
using ..QOperators
using ..Simulations

using Roots: fzero


"""
symbol ("hexagon") we use to initialize coeffs in an OpArray, as a placeholder.
Coefficients initialized with this are intended to be updated dynamically.
"""
const ⎔ = 1.1im   # an arbitrary Complex128


function two_nodes(bas::Basis, ϕ::Real, Δ::Real, κ::Real, Ω::Real, γ::Real,
                   g0::Real, αin::Real, kick::Function, η::Real,
                   kick_signs::Tuple{<:Int,<:Int}=(1, 1))
    const a1 = sparse(MatrixOperator(destroy, bas, 1))
    const b1 = sparse(MatrixOperator(destroy, bas, 2))
    const a2 = sparse(MatrixOperator(destroy, bas, 3))
    const b2 = sparse(MatrixOperator(destroy, bas, 4))
    const a1⁺ = sparse(a1')
    const b1⁺ = sparse(b1')
    const a2⁺ = sparse(a2')
    const b2⁺ = sparse(b2')

    const sqrtκ = sqrt(κ)
    const sqrtγ = sqrt(γ)
    const κhalf = κ / 2
    const κfourth = κ / 4
    const sqrt2inv = 1 / sqrt(2)
    const im_sqrt2inv = im / sqrt(2)
    const cotϕ_2 = cot(ϕ / 2)
    const cotϕ_22 = cotϕ_2 * cotϕ_2

    const a1⁺a1 = sparse(a1⁺ * a1)
    const a2⁺a2 = sparse(a2⁺ * a2)
    const b1⁺b1 = sparse(b1⁺ * b1)
    const b2⁺b2 = sparse(b2⁺ * b2)
    const Ω_b⁺b = sparse(Ω * (b1⁺b1 + b2⁺b2))
    const i_γ_b⁺²_minus_b²_4 = sparse(
        im * (γ / 4) * (b1⁺^2 - b1^2 + b2⁺^2 - b2^2))
    const b1⁺_plus_b1 = sparse(b1⁺ + b1)
    const b1⁺_minus_b1 = sparse(b1⁺ - b1)
    const a1⁺_plus_a1 = sparse(a1⁺ + a1)
    const a1⁺_minus_a1 = sparse(a1⁺ - a1)
    const b2⁺_plus_b2 = sparse(b2⁺ + b2)
    const b2⁺_minus_b2 = sparse(b2⁺ - b2)
    const a2⁺_plus_a2 = sparse(a2⁺ + a2)
    const a2⁺_minus_a2 = sparse(a2⁺ - a2)
    const a1⁺b1⁺_plus_b1 = sparse(a1⁺ * b1⁺_plus_b1)
    const a1b1⁺_plus_b1 = sparse(a1 * b1⁺_plus_b1)
    const a2⁺b2⁺_plus_b2 = sparse(a2⁺ * b2⁺_plus_b2)
    const a2b2⁺_plus_b2 = sparse(a2 * b2⁺_plus_b2)
    const a⁺a_b⁺_plus_b = sparse((a1⁺a1 * b1⁺_plus_b1 + a2⁺a2 * b2⁺_plus_b2))
    const q1 = sparse(sqrt2inv * b1⁺_plus_b1)
    const p1 = sparse(im_sqrt2inv * b1⁺_minus_b1)
    const x1 = sparse(sqrt2inv * a1⁺_plus_a1)
    const y1 = sparse(im_sqrt2inv * a1⁺_minus_a1)
    const q1² = sparse(q1 * q1)
    const p1² = sparse(p1 * p1)
    const x1² = sparse(x1 * x1)
    const y1² = sparse(y1 * y1)
    const q2 = sparse(sqrt2inv * b2⁺_plus_b2)
    const p2 = sparse(im_sqrt2inv * b2⁺_minus_b2)
    const x2 = sparse(sqrt2inv * a2⁺_plus_a2)
    const y2 = sparse(im_sqrt2inv * a2⁺_minus_a2)
    const q2² = sparse(q2 * q2)
    const p2² = sparse(p2 * p2)
    const x2² = sparse(x2 * x2)
    const y2² = sparse(y2 * y2)

    α1 = ⎔
    β1 = ⎔
    α2 = ⎔
    β2 = ⎔


    op_array = [
        [-im,                             i_γ_b⁺²_minus_b²_4],  # 1
        [-im * -g0 * conj(α2),            a2b2⁺_plus_b2],       # 2
        [-im * -g0 * α2 ,                 a2⁺b2⁺_plus_b2],      # 3
        [-im * -g0 * conj(α1),            a1b1⁺_plus_b1],       # 4
        [-im * -g0 * α1,                  a1⁺b1⁺_plus_b1],      # 5
        [-im * -g0,                       a⁺a_b⁺_plus_b],       # 6
        [-im,                             Ω_b⁺b],               # 7
        [-im * (-Δ - 2 * g0 * real(β2)),  a2⁺a2],               # 8
        [-im * (-Δ - 2 * g0 * real(β1)),  a1⁺a1]]               # 9

    function Heff(t, θ)
        α1, β1, α2, β2 = θ
        op_array[2][1] = -im * -g0 * conj(α2)
        op_array[3][1] = -im * -g0 * α2
        op_array[4][1] = -im * -g0 * conj(α1)
        op_array[5][1] = -im * -g0 * α1
        op_array[8][1] = -im * (-Δ - 2 * g0 * real(β2))
        op_array[9][1] = -im * (-Δ - 2 * g0 * real(β1))
        return op_array
    end

    Leff = [-sqrtκ * sqrt2inv * cotϕ_2 * (a1 + a2),
            sqrtκ * sqrt2inv * (a1 - a2),
            sqrtγ * b1, sqrtγ * b2]
    Fs = [im * a1⁺, im * b1⁺, im * a2⁺, im * b2⁺]

    const ks1 = kick_signs[1]
    const ks2 = kick_signs[2]

    function θeom(t, θ)
        α1, β1, α2, β2 = θ
        α1dot = (im * α1 * (Δ + 2 * g0 * real(β1)) + sqrtκ * αin
                 - κfourth * (cotϕ_22 * (α1 + α2) + (α1 - α2)))
        β1dot = im * (-Ω * β1 - γ * imag(β1) + g0 * abs2(α1)
                      + ks1 * sqrt2inv * kick(t))
        α2dot = (im * α2 * (Δ + 2 * g0 * real(β2)) + sqrtκ * αin
                 - κfourth * (cotϕ_22 * (α1 + α2) - (α1 - α2)))
        β2dot = im * (-Ω * β2 - γ * imag(β2) + g0 * abs2(α2)
                      + ks2 * sqrt2inv * kick(t))
        return [α1dot, β1dot, α2dot, β2dot]
    end
    Ys = [-a1, -b1, -a2, -b2]  # THIS IS THE CORRECT ORDER!!

    #Build observer object

    const Es = [
        x1, x1², y1, y1², q1, q1², p1, p1², a1⁺a1, b1⁺b1,
        x2, x2², y2, y2², q2, q2², p2, p2², a2⁺a2, b2⁺b2,
    ]

    return MovingBasis(Heff, Fs, Ys, Leff, θeom, η, 4), Es
end


function single_node(bas::Basis, Δ::Real, κ::Real, Ω::Real, γ::Real, g0::Real,
                     αin::Real, kick::Function, η::Real)
    const a = sparse(MatrixOperator(destroy, bas, 1))
    const b = sparse(MatrixOperator(destroy, bas, 2))
    const a⁺ = sparse(a')
    const b⁺ = sparse(b')

    const sqrtκ = sqrt(κ)
    const sqrtγ = sqrt(γ)
    const κhalf = κ / 2
    const sqrt2inv = 1 / sqrt(2)
    const im_sqrt2inv = im / sqrt(2)

    const a⁺a = sparse(a⁺ * a)
    const b⁺b = sparse(b⁺ * b)
    const Ω_b⁺b = sparse(Ω * b⁺b)
    const b⁺_plus_b = sparse(b⁺ + b)
    const b⁺_minus_b = sparse(b⁺ - b)
    const a⁺_plus_a = sparse(a⁺ + a)
    const a⁺_minus_a = sparse(a⁺ - a)
    const ab⁺_plus_b = sparse(a * b⁺_plus_b)
    const a⁺b⁺_plus_b = sparse(a⁺ * b⁺_plus_b)
    const i_γ_b⁺²_minus_b²_4 = sparse(im * (γ / 4) * (b⁺^2 - b^2))
    const a⁺ab⁺_plus_b = sparse(a⁺a * b⁺_plus_b)
    const q = sparse(sqrt2inv * b⁺_plus_b)
    const p = sparse(im_sqrt2inv * b⁺_minus_b)
    const x = sparse(sqrt2inv * a⁺_plus_a)
    const y = sparse(im_sqrt2inv * a⁺_minus_a)
    const q² = sparse(q * q)
    const p² = sparse(p * p)
    const x² = sparse(x * x)
    const y² = sparse(y * y)

    α = ⎔
    β = ⎔

    op_array = [
        [-im,                            i_γ_b⁺²_minus_b²_4], # 1
        [-im * -g0 * conj(α),            ab⁺_plus_b],         # 2
        [-im * -g0 * α,                  a⁺b⁺_plus_b ],       # 3
        [-im * -g0,                      a⁺ab⁺_plus_b],       # 4
        [-im,                            Ω_b⁺b],              # 5
        [-im * (-Δ - 2 * g0 * real(β)),  a⁺a]]                # 6

    function Heff(t, θ)
        α, β = θ
        op_array[6][1] = -im * (-Δ - 2 * g0 * real(β))
        op_array[3][1] = -im * (-g0 * α)
        op_array[2][1] = -im * (-g0 * conj(α))
        return op_array
    end

    Lfactor = [sqrtκ, sqrtγ]
    Leff = Lfactor.*[a,  b]

    Fs = [im * a⁺, im * b⁺]
    function θeom(t, θ)
        α, β = θ
        αdot = im * α * (Δ + 2 * g0 * real(β)) + sqrtκ * αin - κhalf * α
        βdot = im * (-Ω * β - γ * imag(β) + g0 * abs2(α) + sqrt2inv * kick(t))
        return [αdot, βdot]
    end
    Ys = [-a, -b]  # THIS IS THE CORRECT ORDER!!

    #Build observer object

    const Es = [x, x², y, y², q, q², p, p², a⁺a, b⁺b]

    return MovingBasis(Heff, Fs, Ys, Leff, θeom, η, 2) , Es
end

# Parameter utilities
function eq_cavity_amplitude(q, Δ, κ, g0, αin)
    return sqrt(κ) * αin / ((κ / 2) - im * (Δ + sqrt(2) * g0 * q))
end

function eq_radiation_force(q, Δ, κ, g0, αin)
    return sqrt(2) * g0 * abs2(eq_cavity_amplitude(q, Δ, κ, g0, αin))
end

function restoring_force(q, Ω)
    return -Ω * q
end

function eq_total_force(q, Δ, κ, Ω, g0, αin)
    return restoring_force(q, Ω) + eq_radiation_force(q, Δ, κ, g0, αin)
end

function eq_mech_displacement(Δ, κ, Ω, g0, αin, high::Bool=false)
    function f(qq)
        return eq_total_force(qq, Δ, κ, Ω, g0, αin)
    end
    if high
        qstart = -Δ / (sqrt(2) * g0)
    else
        qstart = 0
    end

    return fzero(f, qstart)
end

end #module

QOperators.jl

__precompile__()
module QOperators
export MatrixOperator, ManifoldOperator, destroy, ptrace

using ..Bases: Basis, length, nfactors

struct ManifoldOperator
    fun::Function
end

const Operator = Union{AbstractMatrix{<:Number}, ManifoldOperator}

function MatrixOperator(op::Function, bas::Basis, ind::Int)
    A = op(bas.dims[ind])
    B = zeros(eltype(A), length(bas), length(bas))
    for (i,ei) in enumerate(bas.e)
        for (j,ej) in enumerate(bas.e)
            if all(n==ind || ei[n]==ej[n] for n in 1:nfactors(bas))
                B[i,j] = getindex(A, ei[ind]+1, ej[ind]+1)
            end
        end
    end
    return B
end

destroy(N::Int) = sparse(diagm(sqrt.(1.:N - 1), 1)) #use sparse?

function ptrace(A::AbstractMatrix, bas::Basis, keep::Vector{<:Int})
    bas_red = Basis(unique([e[keep] for e in bas.e]), bas.dims[keep])
    subinds = [find(e->e[keep]==f, bas.e) for f in bas_red.e]
    A_red = zeros(eltype(A), length(bas_red), length(bas_red))
    for j in 1:length(bas_red)
        for k in 1:length(bas_red)
            subblock = view(A, subinds[j], subinds[k])
            A_red[j,k] = trace(subblock)
        end
    end
    return A_red, bas_red
end

end #module

QuantumManifold.jl

__precompile__()
module QuantumManifold
    include("Bases.jl")
    include("States.jl")
    include("QOperators.jl")
    include("Simulations.jl")
    include("Optomechanics.jl")
end

REQUIRE

julia 0.6

run_single_node.jl

#=
Simulation of a single optomechanical node

Usage:

    julia run_single_node.jl RF

`RF` is the runfolder, which must contain a file `params.dat` that provides
parameters for the simulation, e.g.

Delta: -15          # Optical detuning
kappa: 10           # Optical decay rate
Omega: 1            # Mechanical frequency
gamma: 1.2          # Mechanical decay rate
g0: 0.3162          # Optomechanical coupling
alphain: 15.81      # Optical drive amplitude
phi: 1.482          # Feedback beamsplitter mixing angle
tkick: 2.0          # Time at which to kick the system
kick_amplitude: 25  # Strength of kick
kick_duration: 0.1  # Duration of kick impulse
eta: 120            # Decay rate of deviation from manifold state
Tmax: 15.0          # End time of the simulation
nt: 900000          # Number of time steps
nopt: 6             # Size of optical Hilbert space
nmech: 9            # Size of mechanical Hilbert space
plotstep: 5000      # Save results only every `plotstep` time steps
writestates: true   # Save the full state along with other results?
seed: 1             # Random seed
=#
println("START at " * Dates.format(now(), "yyyy-mm-dd HH:MM:SS"))

using States: ManifoldState, expect
using Bases: Basis, length
using Simulations: integrate, MovingBasis
using DifferentialEquations: DiscreteCallback
using Optomechanics


"""
Read parameters from `runfolder/params.dat`, return them as a tuple
"""
function read_params(runfolder::String)
    print("Reading params...")
    params = Dict{String,String}()
    open(joinpath(runfolder, "params.dat")) do f_params
        for line in eachline(f_params)
            line = replace(line, r"#.*$", "")  # strip out comments
            key, val = map(strip, split(line, ':'))
            params[lowercase(key)] = val
        end
    end
    Δ = float(params["delta"])
    κ = float(params["kappa"])
    Ω = float(params["omega"])
    γ = float(params["gamma"])
    g0 = float(params["g0"])
    αin = float(params["alphain"])
    tkick = float(params["tkick"])
    kick_amplitude = float(params["kick_amplitude"])
    kick_duration = float(params["kick_duration"])
    η = float(params["eta"])
    Tmax = float(params["tmax"])
    nt = parse(Int, params["nt"])
    dt = Tmax / nt
    nopt = parse(Int, params["nopt"])
    nmech = parse(Int, params["nmech"])
    plotstep = parse(Int, get(params, "plotstep", "1"))
    writestates = parse(Bool, lowercase(get(params, "writestates", "false")))
    seed = parse(Int, get(params, "seed", "1"))

    println("Done.")
    return (Δ, κ, Ω, γ, g0, αin, Tmax, nt, dt, tkick, kick_duration,
            kick_amplitude, η, nopt, nmech, seed, writestates, plotstep)
end


"""
Set up the simulation, consisting of `sim` (the equation of motion), the array
of operators for which to calculate expectation values, and the initial state
"""
function initialize_simulation(
        Δ::Number, κ::Number, Ω::Number, γ::Number, g0::Number, αin::Number,
        tkick::Number, kick_duration::Number,
        kick_amplitude::Number, η::Number, nopt::Int, nmech::Int)
    print("Initialize...")

    function gaussian(x)
        return exp(-.5 * x * x) / sqrt(2 * π)
    end

    function pulse(t, t0, dt)
        return gaussian((t - t0) / dt) / dt
    end

    function kick(t)
        return kick_amplitude * pulse(t * Ω, tkick, kick_duration)
    end

    bas = Basis(nopt, nmech)
    sim, Es = single_node(bas, Δ, κ, Ω, γ, g0, αin, kick, η)

    q = eq_mech_displacement(Δ, κ, Ω, g0, αin)
    α = eq_cavity_amplitude(q, Δ, κ, g0, αin)
    θ = Array{Complex128}(2)
    θ[1] = α
    θ[2] = q / sqrt(2)

    groundstate = zeros(Complex128, length(bas))
    groundstate[1] = 1

    Ψ0 = ManifoldState(θ, groundstate)
    println("Done.")
    return sim::MovingBasis, Es::Array, Ψ0::ManifoldState
end


"""
Factory for a `write_output` function that can be used as part of
a `DiscreteCallback` for DifferentialEquations
"""
function make_write_output(
        plotstep::Int, nt::Int, writestates::Bool, sim::MovingBasis, Es::Array,
        f_times::IOStream, f_theta::IOStream, f_expvals::IOStream,
        f_wiener::IOStream, f_states_real=nothing, f_states_imag=nothing)
    n = 0
    const pdim = sim.pdim
    const ndim = length(sim.dA.L) # Noise dimension

    function write_line(t, θ, Ψ, Ws)
        #=println("t = " * repr(t))=#
        write(f_times, "$t\n")
        writedlm(f_theta, transpose(vcat(real(θ), imag(θ))), ',')
        if writestates
            writedlm(f_states_real, transpose(real(Ψ)), ',')
            writedlm(f_states_imag, transpose(imag(Ψ)), ',')
        end
        es = [real(expect(E, Ψ)) for E in Es]
        writedlm(f_expvals, transpose(es), ',')
        writedlm(f_wiener, transpose(vcat(real(Ws), imag(Ws))), ',')
    end

    return function write_output(integrator)
        if n == 0
            # write out initial state
            t = integrator.sol.t[1]
            θ = view(integrator.sol.u[1], 1:pdim)
            Ψ = view(integrator.sol.u[1], pdim + 1:endof(integrator.sol.u[1]))
            Ws = zeros(ndim)
            write_line(t, θ, Ψ, Ws)
        end
        n += 1
        if (mod(n, plotstep) == 0) || (n == nt)
            t = integrator.t
            θ = view(integrator.u, 1:pdim)
            Ψ = view(integrator.u, pdim + 1:endof(integrator.u))
            Ws = conj.(integrator.sol.W[end])
            write_line(t, θ, Ψ, Ws)
        end
    end
end


"""
Run a simulation based on the data in `runfolder/params.dat, write results to
csv files in the `runfolder`. Returns `sol` from DifferentialEquations's `solve`.
"""
function run_simulation(
        runfolder::String, sim::MovingBasis, Es::Array, Ψ0::ManifoldState,
        tspan::Tuple, nt::Int, dt::Number, writestates::Bool, plotstep::Int;
        seed::Int=1)

    print("Initialize I/O...")

    f_times = open(joinpath(runfolder, "./times.csv"), "w")
    f_theta = open(joinpath(runfolder, "./theta.csv"), "w")
    write(f_theta, "x,q,y,p\n")
    f_expvals = open(joinpath(runfolder, "./expvals.csv"), "w")
    write(f_expvals, "<x>,<x^2>,<y>,<y^2>,<q>,<q^2>,<p>,<p^2>,<a+ a>,<b+ b>\n")
    f_wiener = open(joinpath(runfolder, "./wiener.csv"), "w")
    write(f_wiener, "Wa_re,Wb_re,Wa_im,Wb_im\n")

    if writestates
        f_states_real = open(joinpath(runfolder, "./states_real.csv"), "w")
        f_states_imag = open(joinpath(runfolder, "./states_imag.csv"), "w")
    else
        f_states_real = nothing
        f_states_imag = nothing
    end

    cb = DiscreteCallback(
        (u,t,integrator) -> true,
        make_write_output(plotstep, nt, writestates, sim, Es, f_times, f_theta,
                          f_expvals, f_wiener, f_states_real, f_states_imag),
        save_positions=(false,false),
        initialize=(c,u,t,integrator) -> nothing)
    println("Done.")

    println("Simulate...")
    @time sol = integrate(sim, Ψ0, tspan, dt, seed=seed, callback=cb,
                          save_everystep=false, save_noise=true)
    println("...Done.")

    print("Finalize I/O...")
    close(f_times)
    close(f_expvals)
    close(f_theta)
    close(f_wiener)

    if writestates
        close(f_states_real)
        close(f_states_imag)
    end

    println("Done.")
    return sol
end


"""
Pad the given `state` defined on `basis` with zeroes in order to fit on
`new_basis`

The `state` is assumed to be a concatenation of `pdim` classical paramters,
followed by `n_hilbert` entries giving the complex amplitude for each basis
state in `basis`. The output will contain the first `pdim` entries unchanged.
"""
function resize_state(state::Array, basis::Basis, new_basis::Basis):: Array
    n_hilbert = prod(basis.dims)
    new_n_hilbert = prod(new_basis.dims)
    pdim = length(state) - n_hilbert
    @assert length(new_basis.dims) == length(basis.dims)
    @assert new_basis.dims >= basis.dims
    new_state = zeros(Complex128, pdim + new_n_hilbert)
    new_state[1:pdim] = state[1:pdim]
    j = 0
    for (i, etuple) in enumerate(basis.e)
        j += 1
        new_etuple = new_basis.e[j]
        while new_etuple != etuple
            j += 1
            new_etuple = new_basis.e[j]
        end
        new_state[pdim + j] = state[pdim + i]
    end
    return new_state
end


function main(args)

    runfolder = args[1]
    (Δ, κ, Ω, γ, g0, αin, Tmax, nt, dt, tkick, kick_duration, kick_amplitude, η,
     nopt, nmech, seed, writestates, plotstep) = read_params(runfolder)

    tspan = (0.0, Tmax)
    pdim = 2

    iteration = 0
    θerr = 1.0
    Ψerr = 1.0
    const ensure_convergence = true

    if ensure_convergence
        f_convergence = open(joinpath(runfolder, "./convergence.csv"), "w")
        write(f_convergence, "nopt,nmech,θerr,Ψerr\n")
    end

    while θerr > 1e-3 || Ψerr > 1e-3

        iteration += 1
        basis = Basis(nopt, nmech)
        n_hilbert = length(basis.e)
        sim, Es, Ψ0 = initialize_simulation(
            Δ, κ, Ω, γ, g0, αin, tkick, kick_duration, kick_amplitude, η, nopt,
            nmech)
        pdim = sim.pdim
        println("\nIteration $iteration: nmech = $nmech, nopt = $nopt")
        sol = run_simulation(
            runfolder, sim, Es, Ψ0, tspan, nt, dt, writestates, plotstep, seed=seed)
        state = sol.u[end]  # concatenation of θ, Ψ
        if iteration == 1
            if !ensure_convergence
                break
            end
            prev_basis = basis
            prev_state = zeros(Complex128, pdim + n_hilbert)
            prev_state[1:pdim] = 1.0  # so as not divide by zero
        end
        diff = state - resize_state(prev_state, prev_basis, basis)
        θerr = norm(diff[1:pdim]) / norm(prev_state[1:pdim])
        println("norm(Ψ) = $(norm(state[pdim+1:end]))") # DEBUG
        Ψerr = norm(diff[pdim+1:end])
        println("deviation from previous iteration( θ, Ψ): $θerr, $Ψerr")
        if ensure_convergence
            writedlm(f_convergence,
                     transpose([nopt, nmech, θerr, Ψerr]), ',')
        end

        # update for next iteration

        # Out strategy is to increase both nmech and nopt in iteration 2
        # (hopefully, just to confirm that the initial choice was ok). If that
        # doesn't confirm convergence, we alternatingly increase nmech and nopt

        prev_state = state
        prev_basis = basis
        #=if iteration == 1=#
            #=nmech += 1=#
            #=nopt += 1=#
        #=else=#
            #=nmech += iteration % 2=#
            #=nopt += (iteration + 1) % 2=#
        #=end=#

    end

    if ensure_convergence
        close(f_convergence)
    end

end


main(ARGS)

println("FINISHED at " * Dates.format(now(), "yyyy-mm-dd HH:MM:SS"))

run_two_nodes.jl

#=
Simulation of two optomechanical nodes in a feedback network

Usage:

    julia run_two_nodes.jl RF

`RF` is the runfolder, which must contain a file `params.dat` that provides
parameters for the simulation, e.g.

Delta: -15          # Optical detuning
kappa: 10           # Optical decay rate
Omega: 1            # Mechanical frequency
gamma: 1.2          # Mechanical decay rate
g0: 0.3162          # Optomechanical coupling
alphain: 15.81      # Optical drive amplitude
phi: 1.482          # Feedback beamsplitter mixing angle
tkick: 2.0          # Time at which to kick the system
kick_amplitude: 25  # Strength of kick
kick_duration: 0.1  # Duration of kick impulse
eta: 120            # Decay rate of deviation from manifold state
Tmax: 15.0          # End time of the simulation
nt: 900000          # Number of time steps
nopt: 6             # Size of optical Hilbert space
nmech: 9            # Size of mechanical Hilbert space
plotstep: 5000      # Save results only every `plotstep` time steps
writestates: true   # Save the full state along with other results?
seed: 1             # Random seed
=#
println("START at " * Dates.format(now(), "yyyy-mm-dd HH:MM:SS"))

using States: ManifoldState, expect
using Bases: Basis, length
using Simulations: integrate, MovingBasis
using DifferentialEquations: DiscreteCallback
using Optomechanics


"""
Read parameters from `runfolder/params.dat`, return them as a tuple
"""
function read_params(runfolder::String)
    print("Reading params...")
    params = Dict{String,String}()
    open(joinpath(runfolder, "params.dat")) do f_params
        for line in eachline(f_params)
            line = replace(line, r"#.*$", "")  # strip out comments
            key, val = map(strip, split(line, ':'))
            params[lowercase(key)] = val
        end
    end
    Δ = float(params["delta"])
    κ = float(params["kappa"])
    Ω = float(params["omega"])
    γ = float(params["gamma"])
    g0 = float(params["g0"])
    αin = float(params["alphain"])
    ϕ = float(params["phi"])
    tkick = float(params["tkick"])
    kick_amplitude = float(params["kick_amplitude"])
    kick_duration = float(params["kick_duration"])
    η = float(params["eta"])
    Tmax = float(params["tmax"])
    nt = parse(Int, params["nt"])
    dt = Tmax / nt
    nopt = parse(Int, params["nopt"])
    nmech = parse(Int, params["nmech"])
    plotstep = parse(Int, get(params, "plotstep", "1"))
    writestates = parse(Bool, lowercase(get(params, "writestates", "false")))
    seed = parse(Int, get(params, "seed", "1"))

    println("Done.")
    return (Δ, κ, Ω, γ, g0, αin, ϕ, Tmax, nt, dt, tkick, kick_duration,
            kick_amplitude, η, nopt, nmech, seed, writestates, plotstep)
end


"""
Set up the simulation, consisting of `sim` (the equation of motion), the array
of operators for which to calculate expectation values, and the initial state
"""
function initialize_simulation(
        Δ::Number, κ::Number, Ω::Number, γ::Number, g0::Number, αin::Number,
        ϕ::Number, tkick::Number, kick_duration::Number,
        kick_amplitude::Number, η::Number, nopt::Int, nmech::Int)
    print("Initialize...")

    function gaussian(x)
        return exp(-.5 * x * x) / sqrt(2 * π)
    end

    function pulse(t, t0, dt)
        return gaussian((t - t0) / dt) / dt
    end

    function kick(t)
        return kick_amplitude * pulse(t * Ω, tkick, kick_duration)
    end

    bas = Basis(nopt, nmech, nopt, nmech)
    sim, Es = two_nodes(bas, ϕ, Δ, κ, Ω, γ, g0, αin, kick, η)

    q = eq_mech_displacement(Δ, κ, Ω, g0, αin)
    α = eq_cavity_amplitude(q, Δ, κ, g0, αin)
    θ = Array{Complex128}(4)
    θ[1] = α
    θ[2] = q / sqrt(2)
    θ[3] = α
    θ[4] = q / sqrt(2)

    groundstate = zeros(Complex128, length(bas))
    groundstate[1] = 1

    Ψ0 = ManifoldState(θ, groundstate)
    println("Done.")
    return sim::MovingBasis, Es::Array, Ψ0::ManifoldState
end


"""
Factory for a `write_output` function that can be used as part of
a `DiscreteCallback` for DifferentialEquations
"""
function make_write_output(
        plotstep::Int, nt::Int, writestates::Bool, sim::MovingBasis, Es::Array,
        f_times::IOStream, f_theta::IOStream, f_expvals::IOStream,
        f_wiener::IOStream, f_states_real=nothing, f_states_imag=nothing)
    n = 0
    const pdim = sim.pdim
    const ndim = length(sim.dA.L) # Noise dimension

    function write_line(t, θ, Ψ, Ws)
        println("t = " * repr(t))
        write(f_times, "$t\n")
        writedlm(f_theta, transpose(vcat(real(θ), imag(θ))), ',')
        if writestates
            writedlm(f_states_real, transpose(real(Ψ)), ',')
            writedlm(f_states_imag, transpose(imag(Ψ)), ',')
        end
        es = [real(expect(E, Ψ)) for E in Es]
        writedlm(f_expvals, transpose(es), ',')
        writedlm(f_wiener, transpose(vcat(real(Ws), imag(Ws))), ',')
    end

    return function write_output(integrator)
        if n == 0
            # write out initial state
            t = integrator.sol.t[1]
            θ = view(integrator.sol.u[1], 1:pdim)
            Ψ = view(integrator.sol.u[1], pdim + 1:endof(integrator.sol.u[1]))
            Ws = zeros(ndim)
            write_line(t, θ, Ψ, Ws)
        end
        n += 1
        if (mod(n, plotstep) == 0) || (n == nt)
            t = integrator.t
            θ = view(integrator.u, 1:pdim)
            Ψ = view(integrator.u, pdim + 1:endof(integrator.u))
            Ws = conj.(integrator.sol.W[end])
            write_line(t, θ, Ψ, Ws)
        end
    end
end


"""
Run a simulation based on the data in `runfolder/params.dat, write results to
csv files in the `runfolder`. Returns `sol` from DifferentialEquations's `solve`.
"""
function run_simulation(
        runfolder::String, sim::MovingBasis, Es::Array, Ψ0::ManifoldState,
        tspan::Tuple, nt::Int, dt::Number, writestates::Bool, plotstep::Int;
        seed::Int=1)

    print("Initialize I/O...")

    f_times = open(joinpath(runfolder, "./times.csv"), "w")
    f_theta = open(joinpath(runfolder, "./theta.csv"), "w")
    write(f_theta, "x1,q1,x2,q2,y1,p1,y2,p2\n")
    f_expvals = open(joinpath(runfolder, "./expvals.csv"), "w")
    write(f_expvals, "<x1>,<x1^2>,<y1>,<y1^2>,<q1>,<q1^2>,<p1>,<p1^2>,<a1+ a1>,<b1+ b1>,<x2>,<x2^2>,<y2>,<y2^2>,<q2>,<q2^2>,<p2>,<p2^2>,<a2+ a2>,<b2+ b2>\n")
    f_wiener = open(joinpath(runfolder, "./wiener.csv"), "w")
    write(f_wiener, "W1_re,W2_re,W3_re,W4_re,W1_im,W2_im,W3_im,W4_im\n")

    if writestates
        f_states_real = open(joinpath(runfolder, "./states_real.csv"), "w")
        f_states_imag = open(joinpath(runfolder, "./states_imag.csv"), "w")
    else
        f_states_real = nothing
        f_states_imag = nothing
    end

    cb = DiscreteCallback(
        (u,t,integrator) -> true,
        make_write_output(plotstep, nt, writestates, sim, Es, f_times, f_theta,
                          f_expvals, f_wiener, f_states_real, f_states_imag),
        save_positions=(false,false),
        initialize=(c,u,t,integrator) -> nothing)
    println("Done.")

    println("Simulate...")
    @time sol = integrate(sim, Ψ0, tspan, dt, seed=seed, callback=cb,
                          save_everystep=false, save_noise=true)
    println("...Done.")

    print("Finalize I/O...")
    close(f_times)
    close(f_expvals)
    close(f_theta)
    close(f_wiener)

    if writestates
        close(f_states_real)
        close(f_states_imag)
    end

    println("Done.")
    return sol
end


"""
Pad the given `state` defined on `basis` with zeroes in order to fit on
`new_basis`

The `state` is assumed to be a concatenation of `pdim` classical paramters,
followed by `n_hilbert` entries giving the complex amplitude for each basis
state in `basis`. The output will contain the first `pdim` entries unchanged.
"""
function resize_state(state::Array, basis::Basis, new_basis::Basis):: Array
    n_hilbert = prod(basis.dims)
    new_n_hilbert = prod(new_basis.dims)
    pdim = length(state) - n_hilbert
    @assert length(new_basis.dims) == length(basis.dims)
    @assert new_basis.dims >= basis.dims
    new_state = zeros(Complex128, pdim + new_n_hilbert)
    new_state[1:pdim] = state[1:pdim]
    j = 0
    for (i, etuple) in enumerate(basis.e)
        j += 1
        new_etuple = new_basis.e[j]
        while new_etuple != etuple
            j += 1
            new_etuple = new_basis.e[j]
        end
        new_state[pdim + j] = state[pdim + i]
    end
    return new_state
end


function main(args)

    runfolder = args[1]
    (Δ, κ, Ω, γ, g0, αin, ϕ, Tmax, nt, dt, tkick, kick_duration, kick_amplitude, η,
     nopt, nmech, seed, writestates, plotstep) = read_params(runfolder)

    tspan = (0.0, Tmax)
    pdim = 2

    iteration = 0
    θerr = 1.0
    Ψerr = 1.0
    const ensure_convergence = false

    if ensure_convergence
        f_convergence = open(joinpath(runfolder, "./convergence.csv"), "w")
        write(f_convergence, "nopt,nmech,θerr,Ψerr\n")
    end

    while θerr > 1e-3 || Ψerr > 1e-3

        iteration += 1
        basis = Basis(nopt, nmech, nopt, nmech)
        n_hilbert = length(basis.e)
        sim, Es, Ψ0 = initialize_simulation(
            Δ, κ, Ω, γ, g0, αin, ϕ, tkick, kick_duration, kick_amplitude, η, nopt,
            nmech)
        pdim = sim.pdim
        println("\nIteration $iteration: nmech = $nmech, nopt = $nopt")
        sol = run_simulation(
            runfolder, sim, Es, Ψ0, tspan, nt, dt, writestates, plotstep, seed=seed)
        state = sol.u[end]  # concatenation of θ, Ψ
        if iteration == 1
            if !ensure_convergence
                break
            end
            prev_basis = basis
            prev_state = zeros(Complex128, pdim + n_hilbert)
            prev_state[1:pdim] = 1.0  # so as not divide by zero
        end
        diff = state - resize_state(prev_state, prev_basis, basis)
        θerr = norm(diff[1:pdim]) / norm(prev_state[1:pdim])
        println("norm(Ψ) = $(norm(state[pdim+1:end]))") # DEBUG
        Ψerr = norm(diff[pdim+1:end])
        println("deviation from previous iteration( θ, Ψ): $θerr, $Ψerr")
        if ensure_convergence
            writedlm(f_convergence,
                     transpose([nopt, nmech, θerr, Ψerr]), ',')
        end

        # update for next iteration

        # Out strategy is to increase both nmech and nopt in iteration 2
        # (hopefully, just to confirm that the initial choice was ok). If that
        # doesn't confirm convergence, we alternatingly increase nmech and nopt

        prev_state = state
        prev_basis = basis
        if iteration == 1
            nmech += 1
            nopt += 1
        else
            nmech += iteration % 2
            nopt += (iteration + 1) % 2
        end

    end

    if ensure_convergence
        close(f_convergence)
    end

end


main(ARGS)

println("FINISHED at " * Dates.format(now(), "yyyy-mm-dd HH:MM:SS"))

Simulations.jl

__precompile__()
module Simulations

export MovingBasis, SSE, integrate

using DifferentialEquations: SDEProblem, ImplicitEM, solve

using ..States: KetState, ManifoldState, length, expect

#abstract type Measurement end
#struct Heterodyne <: Measurement end

abstract type Evolution end


"""
We generally want to store Hamiltonians/EOMs as a array `[ [coeff, op], ...]`
so that we can efficiently apply it to a state, in-memory, via `apply_op_array`
"""
const OpArray = Vector{Vector{Any}}  # type alias
# We don't use an array of tuples because we'll want to update the
# coeffs dynamically


"""
symbol ("hexagon") we use to initialize coeffs in an OpArray, as a placeholder.
Coefficients initialized with this are intended to be updated dynamically.
"""
const ⎔ = 1.1im   # an arbitrary Complex128


"""
Function to be used as `op` in `OpArray` that simply multipies the state with
the `coeff`. That is, the identity operator in the context of `appy_op_array`
"""
function apply_scalar!(Ψ, a, Φ, Φtemp, add)
    if add
        Φ .+= a * Ψ
    else
        Φ .= a * Ψ
    end
end


"""
Description for the stochastic part of the evolution

This should generally be instantiated as `da = SSE(Ls)` where `Ls` is a list of
Lindblad Operators (Matrices). The other fields of SSE are for internal use
only.
"""
struct SSE <: Evolution
    L::Vector{<:AbstractMatrix{<:Number}}
    mL⁺Lsum_half::AbstractMatrix{<:Number}
    op_array::OpArray

    "Internal constructor for `da = SSE(Ls)` (the recommended instantiation)"
    function SSE(L::AbstractVector{<:AbstractMatrix{<:Number}})
        L⁺ = ctranspose.(L)
        mL⁺Lsum_half = -sum(L⁺.*L) / 2
        n = length(L)
        op_array = OpArray(n+2)
        op_array[1] = [1, mL⁺Lsum_half]
        for i in 1:n
            op_array[i+1] = [⎔, L[i]]
        end
        op_array[n+2] = [⎔, apply_scalar!]
        return new(L, mL⁺Lsum_half, op_array)
    end

end


"""
Use an SSE object `dA` as a function `dA(Ψ, ϕ, Φtemp, add)`
that evaluates

    Φ = (∑⟨L⟩L - ½∑L⁺L - ½∑|⟨L⟩|² ) Ψ   or   Φ += ...  (add = true)

in-place, using Φtemp as a temporary variable
"""
function (dA::SSE)(
        Ψ::AbstractVector{<:Number},
        Φ::AbstractVector{<:Number},
        Φtemp::AbstractVector{<:Number}, add=false)
    scalar = 0.0
    n = length(dA.L)
    for i in 1:n
        Lex = expect(dA.L[i], Ψ, Φtemp)
        dA.op_array[i+1][1] = conj(Lex)
        scalar += abs2(Lex)
    end
    dA.op_array[n+2][1] =  - 0.5 * scalar
    apply_op_array(Ψ, dA.op_array, Φ, Φtemp, add)
end


abstract type Simulation end


struct MovingBasis{
        E<:Evolution,
        T<:AbstractMatrix{<:Number},
        S<:AbstractMatrix{<:Number}} <: Simulation
    H::Function
    Fs::Vector{T}
    Ys::Vector{S}
    dA::E  # Constant
    θeom::Function  # (t,θ)
    η::Real
    pdim::Int
end


noiselength(sim::MovingBasis{<:SSE}) = length(sim.dA.L)


function MovingBasis(H::Function, Fs::Vector{<:AbstractMatrix{<:Number}},
                     Ys::Vector{<:AbstractMatrix{<:Number}},
                     Ls::Vector{<:AbstractMatrix{<:Number}}, θeom::Function,
                     η::Real, pdim::Int)
    return MovingBasis(H, Fs, Ys, SSE(Ls), θeom, η, pdim)
end


function drift(sim::MovingBasis{SSE})
    H = sim.H
    Fs = sim.Fs
    Ys = sim.Ys
    dA = sim.dA
    θeom = sim.θeom
    η = sim.η
    pdim = sim.pdim

    op_array = OpArray(0)
    for F in Fs
        push!(op_array, [⎔, F])
        push!(op_array, [⎔, F'])
    end

    n_hilbert = size(Fs[1], 1)
    temp = Array{Complex128}(n_hilbert)

    return function(t::Real, u::AbstractVector, du::AbstractVector)
        θ = @view u[1:pdim]
        ϕ = @view u[pdim + 1:end]
        dϕ = @view du[pdim + 1:end]
        du[1:pdim] = θeom(t, θ)
        for i in 1:pdim
            θc = -η * expect(Ys[i], ϕ, temp)
            du[i] += θc
            θc *= im
            op_array[2*(i-1) + 1][1] = θc
            op_array[2*(i-1) + 2][1] = -conj(θc)
        end
        apply_op_array(ϕ, op_array, dϕ, temp)
        add = true
        dA(ϕ, dϕ, temp, add)
        apply_op_array(ϕ, H(t, θ), dϕ, temp, add)
    end

end


function diffusion(sim::MovingBasis{SSE})
    dA = sim.dA
    pdim = sim.pdim
    return function(t::Real, u::AbstractVector, g::AbstractMatrix)
        g[1:pdim, :] = 0
        ϕ = @view u[pdim + 1:end]
        for (j, L) in enumerate(dA.L)
            gϕ = @view g[pdim + 1:end, j]
            A_mul_B!(gϕ, L, ϕ)
            gϕ .-= (ϕ' * gϕ) * ϕ
        end
    end
end


function integrate(sim::MovingBasis{SSE}, ψ0::ManifoldState{KetState},
                   tspan::Tuple{<:Real,<:Real}, dt::Real; kws...)
    u0 = vec(ψ0)
    nrp = Matrix{eltype(u0)}(length(u0), noiselength(sim))
    problem = SDEProblem(drift(sim), diffusion(sim), u0, tspan;
                         noise_rate_prototype=nrp)
    return integrate(problem; dt=dt, kws...)
end


function integrate(problem::SDEProblem; kws...)
    #return solve(problem, ImplicitEM(symplectic=true); kws...)
    return solve(problem; kws...)
end


function integrate(sim::Simulation, ntraj::Int; kws...)
    return [integrate(sim; kws...) for i in 1:ntraj]
end


"""
    apply_op_array(Ψ, op_array, Φ, Φtemp, add=false)

Calculate (in-place)

    Φ = ∑ aᵢ Oᵢ Ψ   or   Φ += ∑ aᵢ Oᵢ Ψ   (add=true)

where (aᵢ, Oᵢ) is the i'th element of `op_array`. Φtemp can be used as an
in-memory scratch space if necessary. The cases aᵢ ∈ {0, 1, -1} are handled
efficiently.
"""
function apply_op_array(
        Ψ::AbstractVector, op_array::OpArray, Φ::AbstractVector,
        Φtemp::AbstractVector, add::Bool=false)
    for (a, op) in op_array
        #=@assert a != ⎔=#
        if (abs(a) > 1e-10)
            apply_op(Ψ, a, op, Φ, Φtemp, add)
        else  # skip zero coefficients
            if !add
                Φ .= 0.0
            end
        end
        add = true
    end
end


"""
    apply_op(Ψ, a, Op, Φ, Φtemp, add)

Calculate

    Φ += a⋅Op Ψ     (add=true) or

    Φ = a⋅Op Ψ      (add=false)

in-place. `Op` may be a (sparse) matrix, or a function of the from
`Op(Ψ, a, Φ, Φtemp, add)` that calculates the effect of applying the operator.
In necessary Φtemp may be used as an in-memory scratch space.
"""
function apply_op(
        Ψ::AbstractVector, a::Number, Op::SparseMatrixCSC,
        Φ::AbstractVector, Φtemp::AbstractVector, add::Bool)
    if a == 1  # works for Real/Complex, too
        if add
            A_mul_B!(Φtemp, Op, Ψ)
            Φ .+= Φtemp
        else
            A_mul_B!(Φ, Op, Ψ)
        end
    elseif a == -1
        if add
            A_mul_B!(Φtemp, Op, Ψ)
            Φ .-= Φtemp
        else
            A_mul_B!(Φ, Op, Ψ)
            Φ .*= -1
        end
    else
        if add
            A_mul_B!(Φtemp, Op, Ψ)
            Φ .+= a * Φtemp
        else
            A_mul_B!(Φ, Op, Ψ)
            Φ .*= a
        end
    end
end


function apply_op(
        Ψ::AbstractVector, a::Number, Op::Function,
        Φ::AbstractVector, Φtemp::AbstractVector, add::Bool)
    Op(Ψ, a, Φ, Φtemp, add)
end

end #module

States.jl

__precompile__()
module States

import Base: vec, length

export KetState, ManifoldState, vec, length, pdim, qdim, ket2dm, expect

abstract type State end

length(ψ::State) = length(vec(ψ))

struct KetState <: State
    ket::AbstractVector{<:Number}
end

vec(ψ::KetState) = ψ.ket
qdim(ψ::KetState) = length(ψ)

struct DMState <: State
    dm::AbstractMatrix{<:Number}
end

vec(ψ::DMState) = vec(ψ.dm)  # Returns flattened matrix, useful in (O|S)DEs
qdim(ψ::DMState) = size(ψ.dm, 1)

struct ManifoldState{Q<:State} <: State
    param::AbstractVector{<:Number}
    state::Q
end

function ManifoldState(param::AbstractVector{<:Number},
                       ket::AbstractVector{<:Number})
    return ManifoldState(param, KetState(ket))
end

function ManifoldState(param::AbstractVector{<:Number},
                       dm::AbstractMatrix{<:Number})
    return ManifoldState(param, DMState(dm))
end

function ManifoldState(v::AbstractVector{<:Number}, pdim::Int, dm::Bool=false)
    param = @view v[1:pdim]
    state = @view v[pdim + 1:end]
    if dm
        qdim = round(Int, √length(state))
        state = reshape(state, (qdim, qdim))
    end
    return ManifoldState(param, state)
end

vec(ψ::ManifoldState) = vcat(ψ.param, vec(ψ.state))
pdim(ψ::ManifoldState) = length(ψ.param)
qdim(ψ::ManifoldState) = qdim(ψ.state)

ket2dm(ψ::KetState) = DMState(ψ.ket * ψ.ket')
ket2dm(ψ::ManifoldState{KetState}) = ManifoldState(ψ.param, ket2dm(ψ.ket))

function expect(A::AbstractMatrix{<:Number}, v::AbstractMatrix{<:Number})
    return trace(A * v)
end

function expect(A::AbstractMatrix{<:Number}, v::AbstractVector{<:Number})
    return v' * A * v
end

function expect(A::AbstractMatrix{<:Number}, v::AbstractMatrix{<:Number}, vtemp::AbstractMatrix{<:Number})
    A_mul_B!(vtemp, A, v)
    return trace(vtemp)
end

function expect(A::AbstractMatrix{<:Number}, v::AbstractVector{<:Number}, vtemp::AbstractVector{<:Number})
    A_mul_B!(vtemp, A, v)
    return dot(v, vtemp)
end

# vec(ψ::KetState) = ψ.ket
# vec(ψ::ManifoldState) = vcat(vec(ψ.cstate), vec(ψ.qstate))
# dim(cs::ClassicalState) = length(cs)
# dim(ψ::KetState) = length(ψ.ket)
# mat(A::Operator) = A.op
#
# function *(A::ManifoldOperator, ψ::ManifoldState)
#     return ManifoldState(ψ.cstate, A.fun(ψ.cstate)*ψ.qstate)
# end
#
# *(A::Operator, ψ::KetState) = KetState(mat(A) * vec(ψ))
#
# dot(ψ::ManifoldState, ϕ::ManifoldState) = dot(ψ.qstate, ϕ.qstate)
#
# dot(ψ::KetState, ϕ::KetState) = dot(vec(ψ), vec(ϕ))
#
# const AnyOperator = Union{AbstractMatrix{<:Number}, AbstractOperator}
# const AnyQuantumState = Union{AbstractVector{<:Number}, QuantumState}
# expect(A::AnyOperator, ψ::AnyQuantumState) = dot(ψ, A * ψ)

end #module