text:
module mpc

using QuantumOpticsBase, LinearAlgebra
using ..interaction, ..CollectiveSpins, ..quantum

import ..integrate

try
    using Optim
    global optimize = Optim.optimize
catch e
    if typeof(e) == ArgumentError
        println("Optim package not available. (Needed for calculation of squeezing parameter)")
    else
        rethrow(e)
    end
end

import ..meanfield: densityoperator

export MPCState, densityoperator

const spinbasis = SpinBasis(1//2)
const sigmax_ = dense(sigmax(spinbasis))
const sigmay_ = dense(sigmay(spinbasis))
const sigmaz_ = dense(sigmaz(spinbasis))
const sigmap_ = dense(sigmap(spinbasis))
const sigmam_ = dense(sigmam(spinbasis))

"""
Class describing a MPC state (Product state + Correlations).

The data layout is vector that in matrix form looks like

Cxx Cxy
Cyy Cxz
Czz Cyz

where the Cij are the appropriate correlation matrices.
The expectation values sx, sy and sz are the diagonals of
the matrices Cxx, Cyy and Czz, respectively.

# Arguments
* `N`: Number of spins.
* `data`: Vector of length (3*N)*(2*N+1).
"""
mutable struct MPCState{T1<:Int,T2<:Real}
    N::T1
    data::Vector{T2}
end

"""
    mpc.MPCState(N)

MPC state with all Pauli expectation values and correlations equal to zero.
"""
MPCState(N::T) where T<:Int = MPCState(N, zeros(float(T), (3*N)*(2*N+1)))

"""
    mpc.MPCState(data)

MPC state created from real valued vector of length (3*N)*(2*N+1).
"""
MPCState(data::Vector{<:Real}) = MPCState(dim(data), data)

"""
    mpc.MPCState(rho)

Create MPC state from density operator.
"""
function MPCState(rho::AbstractOperator)
    basis = rho.basis_l
    N = length(basis.bases)
    state = MPCState(N)
    sx, sy, sz, Cxx, Cyy, Czz, Cxy, Cxz, Cyz = splitstate(state)
    f(ind, op) = real(expect(embed(basis, ind, op), rho))
    for k=1:N
        sx[k] = f(k, sigmax_)
        sy[k] = f(k, sigmay_)
        sz[k] = f(k, sigmaz_)
        for l=1:N
            if k==l
                continue
            end
            Cxx[k,l] = f([k,l], [sigmax_, sigmax_])
            Cyy[k,l] = f([k,l], [sigmay_, sigmay_])
            Czz[k,l] = f([k,l], [sigmaz_, sigmaz_])
            Cxy[k,l] = f([k,l], [sigmax_, sigmay_])
            Cxz[k,l] = f([k,l], [sigmax_, sigmaz_])
            Cyz[k,l] = f([k,l], [sigmay_, sigmaz_])
        end
    end
    return state
end

"""
    mpc.blochstate(phi, theta[, N=1])

Product state of `N` single spin Bloch states.

All spins have the same azimuthal angle `phi` and polar angle `theta`.
"""
function blochstate(phi::Vector{T1}, theta::Vector{T2}) where {T1<:Real,T2<:Real}
    N = length(phi)
    @assert length(theta)==N
    state = MPCState(N)
    sx, sy, sz, Cxx, Cyy, Czz, Cxy, Cxz, Cyz = splitstate(state)
    for k=1:N
        sx[k] = cos(phi[k])*sin(theta[k])
        sy[k] = sin(phi[k])*sin(theta[k])
        sz[k] = cos(theta[k])
    end
    for k=1:N, l=1:N
        if k==l
            continue
        end
        Cxx[k,l] = sx[k]*sx[l]
        Cyy[k,l] = sy[k]*sy[l]
        Czz[k,l] = sz[k]*sz[l]
        Cxy[k,l] = sx[k]*sy[l]
        Cxz[k,l] = sx[k]*sz[l]
        Cyz[k,l] = sy[k]*sz[l]
    end
    return state
end

function blochstate(phi::Real, theta::Real, N::Int=1)
    state = MPCState(N)
    sx, sy, sz, Cxx, Cyy, Czz, Cxy, Cxz, Cyz = splitstate(state)
    for k=1:N
        sx[k] = cos(phi)*sin(theta)
        sy[k] = sin(phi)*sin(theta)
        sz[k] = cos(theta)
    end
    for k=1:N, l=1:N
        if k==l
            continue
        end
        Cxx[k,l] = sx[k]*sx[l]
        Cyy[k,l] = sy[k]*sy[l]
        Czz[k,l] = sz[k]*sz[l]
        Cxy[k,l] = sx[k]*sy[l]
        Cxz[k,l] = sx[k]*sz[l]
        Cyz[k,l] = sy[k]*sz[l]
    end
    return state
end

function integersqrt(N::Int)
    n = sqrt(N)
    if abs(trunc(Int, n)-n)>10*eps(n)
        error("N is not a square of an integer.")
    end
    return trunc(Int, n)
end

"""
    mpc.dim(state)

Number of spins described by this state.
"""
function dim(state::Vector{T}) where T<:Real
    x, rem = divrem(length(state), 3)
    @assert rem==0
    N, rem = divrem(-1 + integersqrt(1+8*x), 4)
    @assert rem==0
    return N
end

"""
    mpc.splitstate(N, data)
    mpc.splitstate(state)

Returns sx, sy, sz, Cxx, Cyy, Czz, Cxy, Cxz, Cyz.
"""
function splitstate(N::Int, data::Vector{<:Real})
    data = reshape(data, 3*N, 2*N+1)
    sx = view(data, 1:1*N, 2*N+1)
    sy = view(data, 1*N+1:2*N, 2*N+1)
    sz = view(data, 2*N+1:3*N, 2*N+1)
    Cxx = view(data, 1:1*N, 1:1*N)
    Cyy = view(data, 1*N+1:2*N, 1:1*N)
    Czz = view(data, 2*N+1:3*N, 1:1*N)
    Cxy = view(data, 1:1*N, 1*N+1:2*N)
    Cxz = view(data, 1*N+1:2*N, 1*N+1:2*N)
    Cyz = view(data, 2*N+1:3*N, 1*N+1:2*N)
    return sx, sy, sz, Cxx, Cyy, Czz, Cxy, Cxz, Cyz
end
splitstate(state::MPCState) = splitstate(state.N, state.data)

function covarianceoperator(productstate::Vector{<:DenseOpType}, operators::Vector{<:DenseOpType}, indices::Vector{Int})
    x = DenseOpType[(i in indices ? operators[something(findfirst(isequal(i), indices), 0)] : productstate[i]) for i=1:length(productstate)]
    return tensor(x...)
end

"""
    mpc.correlation2covariance(corstate)

Convert a MPCState from correlation form into covariance form.

Basically it just calculates Cov_ab = <s_a s_b> - <s_a> <s_b>.
"""
function correlation2covariance(corstate::MPCState)
    N = corstate.N
    covstate = MPCState(N)
    sx, sy, sz, Cxx, Cyy, Czz, Cxy, Cxz, Cyz = splitstate(corstate)
    covsx, covsy, covsz, Covxx, Covyy, Covzz, Covxy, Covxz, Covyz = splitstate(covstate)
    for k=1:N
        covsx[k] = sx[k]
        covsy[k] = sy[k]
        covsz[k] = sz[k]
    end
    for k=1:N, l=1:N
        if k==l
            continue
        end
        Covxx[k,l] = Cxx[k,l] - sx[k]*sx[l]
        Covyy[k,l] = Cyy[k,l] - sy[k]*sy[l]
        Covzz[k,l] = Czz[k,l] - sz[k]*sz[l]
        Covxy[k,l] = Cxy[k,l] - sx[k]*sy[l]
        Covxz[k,l] = Cxz[k,l] - sx[k]*sz[l]
        Covyz[k,l] = Cyz[k,l] - sy[k]*sz[l]
    end
    return covstate
end

"""
    mpc.covariance2correlation(covstate)

Convert a MPCState from covariance form into correlation form.

Basically it just calculates <s_a s_b> = Cov_ab + <s_a> <s_b>.
"""
function covariance2correlation(covstate::MPCState)
    N = covstate.N
    corstate = MPCState(N)
    sx, sy, sz, Cxx, Cyy, Czz, Cxy, Cxz, Cyz = splitstate(corstate)
    covsx, covsy, covsz, Covxx, Covyy, Covzz, Covxy, Covxz, Covyz = splitstate(covstate)
    for k=1:N
        sx[k] = covsx[k]
        sy[k] = covsy[k]
        sz[k] = covsz[k]
    end
    for k=1:N, l=1:N
        if k==l
            continue
        end
        Cxx[k,l] = Covxx[k,l] + sx[k]*sx[l]
        Cyy[k,l] = Covyy[k,l] + sy[k]*sy[l]
        Czz[k,l] = Covzz[k,l] + sz[k]*sz[l]
        Cxy[k,l] = Covxy[k,l] + sx[k]*sy[l]
        Cxz[k,l] = Covxz[k,l] + sx[k]*sz[l]
        Cyz[k,l] = Covyz[k,l] + sy[k]*sz[l]
    end
    return corstate
end

"""
    mpc.densityoperator(state)

Create density operator from MPCState.
"""
function densityoperator(state::MPCState)
    N = state.N
    covstate = correlation2covariance(state)
    sx, sy, sz, Cxx, Cyy, Czz, Cxy, Cxz, Cyz = splitstate(covstate)
    productstate = DenseOpType[densityoperator(sx[k], sy[k], sz[k]) for k=1:N]
    C(op1,op2,index1,index2) = covarianceoperator(productstate, [op1,op2], [index1,index2])
    ρ = reduce(tensor, productstate)
    for k=1:N, l=k+1:N
        ρ += 0.25*(
              Cxx[k,l]*C(sigmax_,sigmax_,k,l) + Cxy[l,k]*C(sigmay_,sigmax_,k,l) + Cxz[l,k]*C(sigmaz_,sigmax_,k,l)
            + Cxy[k,l]*C(sigmax_,sigmay_,k,l) + Cyy[k,l]*C(sigmay_,sigmay_,k,l) + Cyz[l,k]*C(sigmaz_,sigmay_,k,l)
            + Cxz[k,l]*C(sigmax_,sigmaz_,k,l) + Cyz[k,l]*C(sigmay_,sigmaz_,k,l) + Czz[k,l]*C(sigmaz_,sigmaz_,k,l))
    end
    return ρ
end

"""
    mpc.sx(state)

Sigma x expectation values of state.
"""
sx(x::MPCState) = view(reshape(x.data, 3*x.N, 2*x.N+1), 0*x.N+1:1*x.N, 2*x.N+1)

"""
    mpc.sy(state)

Sigma y expectation values of state.
"""
sy(x::MPCState) = view(reshape(x.data, 3*x.N, 2*x.N+1), 1*x.N+1:2*x.N, 2*x.N+1)

"""
    mpc.sz(state)

Sigma z expectation values of state.
"""
sz(x::MPCState) = view(reshape(x.data, 3*x.N, 2*x.N+1), 2*x.N+1:3*x.N, 2*x.N+1)

"""
    mpc.Cxx(state)

Sigmax-Sigmax correlation values of MPCState.
"""
Cxx(x::MPCState) = view(reshape(x.data, 3*x.N, 2*x.N+1), 0*x.N+1:1*x.N, 0*x.N+1:1*x.N)
"""
    mpc.Cyy(state)

Sigmay-Sigmay correlation values of MPCState.
"""
Cyy(x::MPCState) = view(reshape(x.data, 3*x.N, 2*x.N+1), 1*x.N+1:2*x.N, 0*x.N+1:1*x.N)
"""
    mpc.Czz(state)

Sigmaz-Sigmaz correlation values of MPCState.
"""
Czz(x::MPCState) = view(reshape(x.data, 3*x.N, 2*x.N+1), 2*x.N+1:3*x.N, 0*x.N+1:1*x.N)
"""
    mpc.Cxy(state)

Sigmax-Sigmay correlation values of MPCState.
"""
Cxy(x::MPCState) = view(reshape(x.data, 3*x.N, 2*x.N+1), 0*x.N+1:1*x.N, 1*x.N+1:2*x.N)
"""
    mpc.Cxz(state)

Sigmax-Sigmaz correlation values of MPCState.
"""
Cxz(x::MPCState) = view(reshape(x.data, 3*x.N, 2*x.N+1), 1*x.N+1:2*x.N, 1*x.N+1:2*x.N)
"""
    mpc.Cyz(state)

Sigmay-Sigmaz correlation values of MPCState.
"""
Cyz(x::MPCState) = view(reshape(x.data, 3*x.N, 2*x.N+1), 2*x.N+1:3*x.N, 1*x.N+1:2*x.N)

"""
    mpc.correlation(s1, s2, s3, C12, C13, C23)

3-spin correlation value.
"""
function correlation(s1::T, s2::T, s3::T, C12::T, C13::T, C23::T) where T<:Real
    return -2.0*s1*s2*s3 + s1*C23 + s2*C13 + s3*C12
end

"""
    mpc.timeevolution(T, S::SpinCollection, state0[; fout])

MPC time evolution.

# Arguments
* `T`: Points of time for which output will be generated.
* `S`: SpinCollection describing the system.
* `state0`: Initial MPCState.
* `fout` (optional): Function with signature fout(t, state) that is called
    whenever output should be generated.
"""
function timeevolution(T, S::SpinCollection, state0::MPCState; fout=nothing, kwargs...)
    N = length(S.spins)
    @assert N==state0.N
    Ω = interaction.OmegaMatrix(S)
    Γ = interaction.GammaMatrix(S)

    function f(dy, y, p, t)
        sx, sy, sz, Cxx, Cyy, Czz, Cxy, Cxz, Cyz = splitstate(N, y)
        dsx, dsy, dsz, dCxx, dCyy, dCzz, dCxy, dCxz, dCyz = splitstate(N, dy)
        @inbounds for k=1:N
            dsx[k] = -0.5*S.gammas[k]*sx[k]
            dsy[k] = -0.5*S.gammas[k]*sy[k]
            dsz[k] = -S.gammas[k]*(1.0+sz[k])
            for j=1:N
                if j==k
                    continue
                end
                dsx[k] += Ω[k,j]*Cyz[j,k] + 0.5*Γ[k,j]*Cxz[j,k]
                dsy[k] += -Ω[k,j]*Cxz[j,k] + 0.5*Γ[k,j]*Cyz[j,k]
                dsz[k] += Ω[k,j]*(Cxy[j,k]-Cxy[k,j]) - 0.5*Γ[k,j]*(Cxx[j,k] + Cyy[j,k])
            end
        end
        @inbounds for k=1:N, l=1:N
            if k==l
                continue
            end
            dCxx[k,l] = -sqrt(S.gammas[k]*S.gammas[l])*Cxx[k,l] + Γ[k,l]*(Czz[k,l]+0.5*sz[k]+0.5*sz[l])
            dCyy[k,l] = -sqrt(S.gammas[k]*S.gammas[l])*Cyy[k,l] + Γ[k,l]*(Czz[k,l]+0.5*sz[k]+0.5*sz[l])
            dCzz[k,l] = -2*sqrt(S.gammas[k]*S.gammas[l])*Czz[k,l] - sqrt(S.gammas[k]*S.gammas[l])*(sz[k]+sz[l]) + Γ[k,l]*(Cyy[k,l]+Cxx[k,l])
            dCxy[k,l] = Ω[k,l]*(sz[k]-sz[l]) - sqrt(S.gammas[k]*S.gammas[l])*Cxy[k,l]
            dCxz[k,l] = Ω[k,l]*sy[l] - 1.5*sqrt(S.gammas[k]*S.gammas[l])*Cxz[k,l] - S.gammas[k]*sx[k] - Γ[k,l]*(Cxz[l,k]+0.5*sx[l])
            dCyz[k,l] = -Ω[k,l]*sx[l] - 1.5*sqrt(S.gammas[k]*S.gammas[l])*Cyz[k,l] - S.gammas[k]*sy[k] - Γ[k,l]*(Cyz[l,k]+0.5*sy[l])
            for j=1:N
                if j==l || j==k
                    continue
                end
                Cxxx = correlation(sx[k], sx[l], sx[j], Cxx[k,l], Cxx[k,j], Cxx[l,j])
                Cxxy = correlation(sx[k], sx[l], sy[j], Cxx[k,l], Cxy[k,j], Cxy[l,j])
                Cxyx = correlation(sx[k], sy[l], sx[j], Cxy[k,l], Cxx[k,j], Cxy[j,l])
                Cxyy = correlation(sx[k], sy[l], sy[j], Cxy[k,l], Cxy[k,j], Cyy[l,j])
                Cxzx = correlation(sx[k], sz[l], sx[j], Cxz[k,l], Cxx[k,j], Cxz[j,l])
                Cxzy = correlation(sx[k], sz[l], sy[j], Cxz[k,l], Cxy[k,j], Cyz[j,l])
                Cyxx = correlation(sy[k], sx[l], sx[j], Cxy[l,k], Cxy[j,k], Cxx[l,j])
                Cyxy = correlation(sy[k], sx[l], sy[j], Cxy[l,k], Cyy[k,j], Cxy[l,j])
                Cyyx = correlation(sy[k], sy[l], sx[j], Cyy[k,l], Cxy[j,k], Cxy[j,l])
                Cyyy = correlation(sy[k], sy[l], sy[j], Cyy[k,l], Cyy[k,j], Cyy[l,j])
                Cyzx = correlation(sy[k], sz[l], sx[j], Cyz[k,l], Cxy[j,k], Cxz[j,l])
                Cyzy = correlation(sy[k], sz[l], sy[j], Cyz[k,l], Cyy[k,j], Cyz[j,l])
                Czxx = correlation(sz[k], sx[l], sx[j], Cxz[l,k], Cxz[j,k], Cxx[l,j])
                Czxy = correlation(sz[k], sx[l], sy[j], Cxz[l,k], Cyz[j,k], Cxy[l,j])
                Czyx = correlation(sz[k], sy[l], sx[j], Cyz[l,k], Cxz[j,k], Cxy[j,l])
                Czyy = correlation(sz[k], sy[l], sy[j], Cyz[l,k], Cyz[j,k], Cyy[l,j])
                Czzx = correlation(sz[k], sz[l], sx[j], Czz[k,l], Cxz[j,k], Cxz[j,l])
                Czzy = correlation(sz[k], sz[l], sy[j], Czz[k,l], Cyz[j,k], Cyz[j,l])

                dCxx[k,l] += Ω[k,j]*Czxy + Ω[l,j]*Cxzy + 0.5*Γ[k,j]*Czxx + 0.5*Γ[l,j]*Cxzx
                dCyy[k,l] += -Ω[k,j]*Czyx - Ω[l,j]*Cyzx + 0.5*Γ[k,j]*Czyy + 0.5*Γ[l,j]*Cyzy
                dCzz[k,l] += (Ω[k,j]*(Cyzx-Cxzy) + Ω[l,j]*(Czyx-Czxy)
                                - 0.5*Γ[k,j]*(Cxzx+Cyzy) - 0.5*Γ[l,j]*(Czxx+Czyy))
                dCxy[k,l] += Ω[k,j]*Czyy - Ω[l,j]*Cxzx + 0.5*Γ[k,j]*Czyx + 0.5*Γ[l,j]*Cxzy
                dCxz[k,l] += (Ω[k,j]*Czzy + Ω[l,j]*(Cxyx-Cxxy)
                                + 0.5*Γ[k,j]*Czzx - 0.5*Γ[l,j]*(Cxxx+Cxyy))
                dCyz[k,l] += (-Ω[k,j]*Czzx + Ω[l,j]*(Cyyx-Cyxy)
                                + 0.5*Γ[k,j]*Czzy - 0.5*Γ[l,j]*(Cyxx+Cyyy))
            end
        end
    end

    if isa(fout, Nothing)
        fout_ = (t, state) -> deepcopy(state)
    else
        fout_ = fout
    end

    return integrate(T, f, state0, fout_; kwargs...)
end

function axisangle2rotmatrix(axis::Vector{T}, angle::Real) where T<:Real
    x, y, z = axis
    c = cos(angle)
    s = sin(angle)
    t = 1-c
    R = zeros(T, 3, 3)
    R[1,1] = t*x^2 + c
    R[1,2] = t*x*y - z*s
    R[1,3] = t*x*z + y*s
    R[2,1] = t*x*y + z*s
    R[2,2] = t*y^2 + c
    R[2,3] = t*y*z - x*s
    R[3,1] = t*x*z - y*s
    R[3,2] = t*y*z + x*s
    R[3,3] = t*z^2 + c
    return R
end

"""
    mpc.rotate(axis, angles, state)

Rotations on the Bloch sphere for the given [`MPCState`](@ref).

# Arguments
* `axis`: Rotation axis.
* `angles`: Rotation angle(s).
* `state`: [`MPCState`](@ref) that should be rotated.
"""
function rotate(axis::Vector{T1}, angles::Vector{T2}, state::MPCState) where {T1<:Real, T2<:Real}
    N = state.N
    @assert length(axis)==3
    @assert length(angles)==N
    w = axis/norm(axis)
    sx, sy, sz, Cxx, Cyy, Czz, Cxy, Cxz, Cyz = splitstate(state)
    rotstate = deepcopy(state)
    sx_rot, sy_rot, sz_rot, Cxx_rot, Cyy_rot, Czz_rot, Cxy_rot, Cxz_rot, Cyz_rot = splitstate(rotstate)
    S_rot = splitstate(rotstate)
    R = [axisangle2rotmatrix(w, angle) for angle=angles]
    for k=1:N
        sx_rot[k], sy_rot[k], sz_rot[k] = R[k]*[sx[k], sy[k], sz[k]]
    end
    for k=1:N,l=1:N
        if k==l
            continue
        end
        S2 = [Cxx[k,l],Cxy[k,l],Cxz[k,l],
              Cxy[l,k],Cyy[k,l],Cyz[k,l],
              Cxz[l,k],Cyz[l,k],Czz[k,l]]
        S2_rot = kron(R[k], R[l])*S2
        Cxx_rot[k,l] = S2_rot[1]
        Cyy_rot[k,l] = S2_rot[5]
        Czz_rot[k,l] = S2_rot[9]
        Cxy_rot[k,l] = S2_rot[2]
        Cxz_rot[k,l] = S2_rot[3]
        Cyz_rot[k,l] = S2_rot[6]
    end
    return rotstate
end

rotate(axis::Vector{T}, angle::Real, state::MPCState) where {T<:Real} = rotate(axis, ones(float(T), state.N)*angle, state)

"""
    mpc.var_Sx(state0)

Variance of the total Sx operator for the given MPCState.
"""
function var_Sx(state0::MPCState)
    N = state0.N
    sx, sy, sz, Cxx, Cyy, Czz, Cxy, Cxz, Cyz = splitstate(state0)
    exp_Sx_2 = 0.
    exp_Sx2 = 0.
    for k=1:N, l=1:N
        exp_Sx_2 += sx[k]*sx[l]
        if k==l
            continue
        end
        exp_Sx2 += Cxx[k,l]
    end
    return 1.0/N + 1.0/N^2*exp_Sx2 - 1.0/N^2*exp_Sx_2
end

"""
    mpc.var_Sy(state)

Variance of the total Sy operator for the given MPCState.
"""
function var_Sy(state0::MPCState)
    N = state0.N
    sx, sy, sz, Cxx, Cyy, Czz, Cxy, Cxz, Cyz = splitstate(state0)
    exp_Sy_2 = 0.
    exp_Sy2 = 0.
    for k=1:N, l=1:N
        exp_Sy_2 += sy[k]*sy[l]
        if k==l
            continue
        end
        exp_Sy2 += Cyy[k,l]
    end
    return 1.0/N + 1.0/N^2*exp_Sy2 - 1.0/N^2*exp_Sy_2
end

"""
    mpc.var_Sz(state)

Variance of the total Sz operator for the given MPCState.
"""
function var_Sz(state0::MPCState)
    N = state0.N
    sx, sy, sz, Cxx, Cyy, Czz, Cxy, Cxz, Cyz = splitstate(state0)
    exp_Sz_2 = 0.
    exp_Sz2 = 0.
    for k=1:N, l=1:N
        exp_Sz_2 += sz[k]*sz[l]
        if k==l
            continue
        end
        exp_Sz2 += Czz[k,l]
    end
    return 1.0/N + 1.0/N^2*exp_Sz2 - 1.0/N^2*exp_Sz_2
end

"""
    mpc.squeeze_sx(χT, state)

Spin squeezing along sx.

# Arguments
* `χT`: Squeezing strength.
* `state0`: MPCState that should be squeezed.
"""
function squeeze_sx(χT::Real, state0::MPCState)
    T = [0,1.]
    N = state0.N
    χeff = 4.0*χT/N^2
    function f(dy, y, p, t)
        sx, sy, sz, Cxx, Cyy, Czz, Cxy, Cxz, Cyz = splitstate(N, y)
        dsx, dsy, dsz, dCxx, dCyy, dCzz, dCxy, dCxz, dCyz = splitstate(N, dy)
        for k=1:N
            dsx[k] = 0.
            dsy[k] = 0.
            dsz[k] = 0.
            for j=1:N
                if j==k
                    continue
                end
                dsy[k] += -χeff*Cxz[j,k]
                dsz[k] += χeff*Cxy[j,k]
            end
        end
        for k=1:N, l=1:N
            if k==l
                continue
            end
            dCxx[k,l] = 0.
            dCyy[k,l] = 0.
            dCzz[k,l] = 0.
            dCxy[k,l] = -χeff*sz[l]
            dCxz[k,l] = χeff*sy[l]
            dCyz[k,l] = 0.
            for j=1:N
                if j==l || j==k
                    continue
                end
                Cxyx = correlation(sx[k], sy[l], sx[j], Cxy[k,l], Cxx[k,j], Cxy[j,l])
                Cxzx = correlation(sx[k], sz[l], sx[j], Cxz[k,l], Cxx[k,j], Cxz[j,l])
                Cyyx = correlation(sy[k], sy[l], sx[j], Cyy[k,l], Cxy[j,k], Cxy[j,l])
                Cyzx = correlation(sy[k], sz[l], sx[j], Cyz[k,l], Cxy[j,k], Cxz[j,l])
                Czyx = correlation(sz[k], sy[l], sx[j], Cyz[l,k], Cxz[j,k], Cxy[j,l])
                Czzx = correlation(sz[k], sz[l], sx[j], Czz[k,l], Cxz[j,k], Cxz[j,l])

                dCyy[k,l] += -χeff*(Czyx+Cyzx)
                dCzz[k,l] += χeff*(Cyzx+Czyx)
                dCxy[k,l] += -χeff*Cxzx
                dCxz[k,l] += χeff*Cxyx
                dCyz[k,l] += χeff*(Czzx-Cyyx)
            end
        end
    end

    fout_(t, state::MPCState) = deepcopy(state)
    time_out, state_out = integrate(T, f, state0, fout_)

    return state_out[end]
end

"""
    mpc.squeeze(axis, χT, state0)

Spin squeezing along an arbitrary axis.

# Arguments
* `axis`: Squeezing axis.
* `χT`: Squeezing strength.
* `state0`: MPCState that should be squeezed.
"""
function squeeze(axis::Vector{T}, χT::Real, state0::MPCState) where T<:Real
    @assert length(axis)==3
    axis = axis/norm(axis)
    # Rotation into sigma_x
    w = float(T)[0., axis[3], -axis[2]]
    if norm(w)<1e-5
        state_squeezed = squeeze_sx(χT, state0)
    else
        α = acos(axis[1])
        state_rot = rotate(w, α, state0)
        state_rot_squeezed = squeeze_sx(χT, state_rot)
        state_squeezed = rotate(w, -α, state_rot_squeezed)
    end
    return state_squeezed
end

"""
    mpc.orthogonal_vectors(n)

Create 3 orthonormal vectors where one is in the given direction `n`.
"""
function orthogonal_vectors(n::Vector{<:Real})
    n = n/norm(n)
    v = (n[1]<n[2] ? [1.,0.,0.] : [0.,1.,0.])
    e1 = v - dot(n,v)*n
    e1 = e1/norm(e1)
    e2 = cross(n, e1)
    e2 = e2/norm(e2)
    return e1, e2
end

"""
    mpc.squeezing_parameter(state)

Calculate squeezing parameter for the given `state`.
"""
function squeezingparameter(state::MPCState)
    N = state.N
    sx, sy, sz, Cxx, Cyy, Czz, Cxy, Cxz, Cyz = map(sum, splitstate(state))
    n = 1.0/N*[sx, sy, sz]
    e1, e2 = orthogonal_vectors(n)
    function f(phi)
        nphi = cos(phi)*e1 + sin(phi)*e2
        nx = nphi[1]
        ny = nphi[2]
        nz = nphi[3]
        Sphi2 = 1.0/N^2*(N+nx*nx*Cxx + ny*ny*Cyy + nz*nz*Czz +
                2*nx*ny*Cxy + 2*nx*nz*Cxz + 2*ny*nz*Cyz)
        return Sphi2
    end
    varSmin = Optim.optimize(f, 0., 2*pi).f_minimum
    return sqrt(N*varSmin)/norm(n)
end

end # module


text:
function unwrap!(wrapped::AbstractArray{T,3}; regions=zeros(UInt8, size(wrapped)), keyargs...) where T
    weights = calculateweights(wrapped; keyargs...)
    @assert sum(weights) != 0 "Unwrap-weights are all zero!"

    regions .= grow_region_unwrap!(wrapped, weights; keyargs...)

    if haskey(keyargs, :correctglobal) && keyargs[:correctglobal]
        mask = if haskey(keyargs, :mask)
            keyargs[:mask]
        else
            trues(size(wrapped))
        end
        wrapped .-= (2π * median(round.(filter(isfinite, wrapped[mask]) ./ 2π))) # TODO time (sample)
    end
    return wrapped
end

function unwrap!(wrapped::AbstractArray; keyargs...)
    sz = size(wrapped)
    @assert ndims(wrapped) <= 3 "This is 3D (4D) unwrapping! data is $(ndims(wrapped))D"
    if ndims(wrapped) <= 2 # algorithm requires 3D input
        wrapped = reshape(wrapped, size(wrapped)..., ones(Int, 3-ndims(wrapped))...)
    end
    unwrap!(wrapped; keyargs...)
    return reshape(wrapped, sz)
end

"""
    unwrap(wrapped::AbstractArray; keyargs...)

ROMEO unwrapping.

###  Optional keyword arguments:

- `weights`: Options are [`:romeo`] | `:romeo2` | `:romeo3` | `:bestpath`.
- `mag`: The magnitude is used to improve the unwrapping-path.
- `mask`: Unwrapping is only performed inside the mask.
- `phase2`: A second reference phase image (possibly with different echo time).
    It is used for calculating the phasecoherence weight. This is automatically
    done for 4D multi-echo input and therefore not required.
- `TEs`: The echo times for multi-echo data. In the case of single-echo
    application with phase and the phase2 as a tuple (eg. (5, 10) or [5, 10]).
- `correctglobal=false`: If `true` corrects global n2π offsets.
- `individual=false`: If `true` perform individual unwrapping of echos.
- `template=2`: echo that is spatially unwrapped (if `individual` is `false`)
- `maxseeds=1`: higher values allow more seperate regions
- `merge_regions=false`: spatially merge neighboring regions after unwrapping
- `correct_regions=false`: bring each regions median closest to 0 by adding n2π
- `wrap_addition=0`: [0;π], allows 'linear unwrapping', neighbors can have more
    (π+wrap_addition) phase difference
- `temporal_uncertain_unwrapping=false`: uses spatial unwrapping on voxels that
    have high uncertainty values after temporal unwrapping

"""
unwrap(wrapped; keyargs...) = unwrap!(copy(wrapped); keyargs...)

function unwrap!(wrapped::AbstractArray{T,4}; TEs, individual=false,
        template=2, p2ref=ifelse(template==1, 2, template-1),
        temporal_uncertain_unwrapping=false, keyargs...) where T
    if individual return unwrap_individual!(wrapped; TEs, keyargs...) end
    ## INIT
    args = Dict{Symbol, Any}(keyargs)
    args[:phase2] = wrapped[:,:,:,p2ref]
    args[:TEs] = TEs[[template, p2ref]]
    if haskey(args, :mag)
        args[:mag] = args[:mag][:,:,:,template]
    end
    ## Calculate
    weights = calculateweights(view(wrapped,:,:,:,template); args...)
    unwrap!(view(wrapped,:,:,:,template); args..., weights) # rightmost keyarg takes precedence
    quality = similar(wrapped)
    V = falses(size(wrapped))
    for ieco in [(template-1):-1:1; (template+1):length(TEs)]
        iref = if (ieco < template) ieco+1 else ieco-1 end
        refvalue = wrapped[:,:,:,iref] .* (TEs[ieco] / TEs[iref])
        w = view(wrapped,:,:,:,ieco)
        w .= unwrapvoxel.(w, refvalue) # temporal unwrapping

        if temporal_uncertain_unwrapping # TODO extract as function
            quality[:,:,:,ieco] .= getquality.(w, refvalue)
            visited = quality[:,:,:,ieco] .< π/2
            mask = if haskey(keyargs, :mask)
                keyargs[:mask]
            else
                dropdims(sum(weights; dims=1); dims=1) .< 100
            end
            visited[.!mask] .= true
            V[:,:,:,ieco] = visited
            if any(visited) && !all(visited)
                edges = getseededges(visited)
                edges = filter(e -> weights[e] != 0, edges)
                grow_region_unwrap!(w, weights, visited, initqueue(edges, weights))
            end
        end
    end
    return wrapped#, quality, weights, V
end

function getquality(vox, ref)
    return abs(vox - ref)
end

function getseededges(visited::BitArray)
    stridelist = getdimoffsets(visited)
    edges = Int64[]
    for dim in 1:3, I in LinearIndices(visited)
        J = I + stridelist[dim]
        if checkbounds(Bool, visited, J) # borders should be no problem due to invalid weights
            if visited[I] + visited[J] == 1 # one visited and one not visited
                push!(edges, getedgeindex(I, dim))
            end
        end
    end
    return edges
end

initqueue(seed::Int, weights) = initqueue([seed], weights)
function initqueue(seeds, weights)
    pq = PQueue{eltype(seeds)}(NBINS)
    for seed in seeds
        enqueue!(pq, seed, weights[seed])
    end
    return pq
end

unwrap_individual(wrapped; keyargs...) = unwrap_individual!(copy(wrapped); keyargs...)
function unwrap_individual!(wrapped::AbstractArray{T,4}; TEs, keyargs...) where T
    args = Dict{Symbol,Any}(keyargs)
    for i in 1:length(TEs)
        e2 = if (i == 1) 2 else i-1 end
        if haskey(keyargs, :mag) args[:mag] = keyargs[:mag][:,:,:,i] end
        unwrap!(view(wrapped,:,:,:,i); phase2=wrapped[:,:,:,e2], TEs=TEs[[i,e2]], args...)
    end
    return wrapped
end

ext: jl


text:
function fprop!(::RelaxInterval, vt::Variable, g::DAT, b::IntervalCache{T}, k) where T<:Real
    i = first_index(g, k)
    l = lbd(b, i)
    u = ubd(b, i)
    b[k] = (l == u) ? Interval(l) : Interval(l, u)
    nothing
end

function fprop!(t::RelaxInterval, ex::Subexpression, g::DAT, b::IntervalCache{T}, k) where T<:Real
    b[k] = subexpression_set(t, b, first_index(g, k))
end

for (F, f) in ((DIV, :/), (ARH, :arh), (POW, :^))
    @eval fprop!(t::RelaxInterval, v::Val{$F}, g::DAT, b::IntervalCache{T}, k) where T<:Real = (b[k] = ($f)(set(t, b, child(g, 1, k)), set(b, child(g, 2, k))); nothing)
end

function fprop!(t::RelaxInterval, v::Val{MINUS}, g::DAT, b::IntervalCache{T}, k) where T<:Real
    x = child(g, 1, k)
    b[k] = is_binary(g, k) ? (set(t, b, x) - set(t, b, child(g, 2, k))) : - set(t, b, x)
end

for (F, f) in ((PLUS, :sum), (MIN, :minimum), (MAX, :maximum), (MULT, :prod))
    @eval fprop!(t::Interval, v::Val{$F}, g::DAT, b::IntervalCache{T}, k) where T<:Real = (b[k] = ($f)(i -> set(t, b, i), children(g, k)); nothing)
end

function fprop!(t::RelaxInterval, v::Val{USER}, g::DAT, b::IntervalCache{T}, k) where T<:Real
    f = user_univariate_operator(g, first_index(g, k))
    b[k] = f(set(t, b, child(g, 1, k)))
end

function fprop!(t::RelaxInterval, v::Val{USERN}, g::DAT, b::IntervalCache{T}, k) where T<:Real
    mv = user_multivariate_operator(g, first_index(g, k))
    set_input = set_input(t, b, arity(g, k))
    for c in children(g, k)
        set_input[i] = set(t, b, c)
    end
    b[k] = MOI.eval_objective(mv, set_input)::Interval{Float64}
end

for ft in UNIVARIATE_ATOM_TYPES
    f = UNIVARIATE_ATOM_DICT[ft]
    (f == :user || f == :+ || f == :-) && continue
    @eval function fprop!(t::RelaxInterval, v::Val{$ft}, g::DAT, b::IntervalCache{T}, k) where T<:Real
        x = child(g, 1, k)
        b[k] = ($f)(set(t, b, x))
    end
end

for (F, f) in ((LOWER_BND, :lower_bnd), (UPPER_BND, :upper_bnd))
    @eval function fprop!(t::RelaxInterval, v::Val{$F}, g::DAT, b::IntervalCache{T}, k) where T<:Real
        y = child(g, 2, k)
        if is_num(t, b, y)
            z = set(t, b, child(g, 1, k))
            b[k] = ($f)(z, num(t, b, y))
        end
        nothing
    end
end

function fprop!(t::RelaxInterval, v::Val{BND}, g::DAT, b::IntervalCache{T}, k) where T<:Real
    y = child(g, 2, k)
    r = child(g, 3, k)
    if is_num(t, b, y) && is_num(t, b, r)
        z = set(t, b, child(g, 1, k))
        b[k] = bnd(z, num(t, b, y), num(t, b, r))
    end
    nothing
end

max_stars_repo_name: PSORLab/RSActivationFunctions

max_issues_repo_name: PSORLab/RSActivationFunctions

max_forks_repo_name: MatthewStuber/EAGO


text:
abstract type AbstractChargeDensity{T <: SSDFloat} end

@inline function ChargeDensity(T::DataType, dict::Union{Dict{String, Any}, Dict{Any, Any}}, inputunit_dict::Dict)
    return ChargeDensity(T, Val{Symbol(dict["name"])}(), dict, inputunit_dict)
end

"""
    struct LinearChargeDensity{T <: SSDFloat} <: AbstractChargeDensity{T}

Simple charge density model which assumes a linear gradient in charge density in each dimension of a Cartesian coordinate system.
`offsets::NTuple{3, T}` are the charge densities at 0 and `gradients::NTuple{3, T}` are the linear slopes in each dimension, `x`, `y` and `z`.
"""
struct LinearChargeDensity{T <: SSDFloat} <: AbstractChargeDensity{T}
    offsets::NTuple{3, T}
    gradients::NTuple{3, T}
end

function ChargeDensity(T::DataType, t::Val{:linear}, dict::Union{Dict{String, Any}, Dict{Any, Any}}, inputunit_dict::Dict)
    unit_factor::T = 1
    gradient_unit_factor::T = 1
    if haskey(inputunit_dict, "length")
        lunit = inputunit_dict["length"]
        unit_factor = inv(ustrip(uconvert( internal_length_unit^3, 1 * lunit^3 )))
        gradient_unit_factor = inv(ustrip(uconvert( internal_length_unit^4, 1 * lunit^4 )))
    end
    return LinearChargeDensity{T}( dict, unit_factor, gradient_unit_factor )
end

function LinearChargeDensity{T}(dict::Union{Dict{String, Any}, Dict{Any, Any}}, unit_factor::T, gradient_unit_factor::T)::LinearChargeDensity{T} where {T <: SSDFloat}
    offsets, gradients = zeros(T,3), zeros(T,3)
    if prod(map(k -> k in ["x","y","z"], collect(keys(dict)))) @warn "Only x, y and z are supported in the linear charge density model.\nChange the charge density model in the config file or remove all other entries." end
    if haskey(dict, "x")     offsets[1] = geom_round(unit_factor * T(dict["x"]["init"]));     gradients[1] = geom_round(gradient_unit_factor * T(dict["x"]["gradient"]))    end
    if haskey(dict, "y")     offsets[2] = geom_round(unit_factor * T(dict["y"]["init"]));     gradients[2] = geom_round(gradient_unit_factor * T(dict["y"]["gradient"]))    end
    if haskey(dict, "z")     offsets[3] = geom_round(unit_factor * T(dict["z"]["init"]));     gradients[3] = geom_round(gradient_unit_factor * T(dict["z"]["gradient"]))    end
    LinearChargeDensity{T}( NTuple{3, T}(offsets), NTuple{3, T}(gradients) )
end

function get_charge_density(lcdm::LinearChargeDensity{T}, pt::AbstractCoordinatePoint{T})::T where {T <: SSDFloat}
    pt::CartesianPoint{T} = CartesianPoint(pt)
    ρ::T = 0
    for i in eachindex(lcdm.offsets)
        ρ += (lcdm.offsets[i] + pt[i] * lcdm.gradients[i]) #* T(1e16) # * T(1e10) * T(1e6) -> 1/cm^3 -> 1/m^3
    end
    return ρ
end

"""
    struct CylindricalChargeDensity{T <: SSDFloat} <: AbstractChargeDensity{T}

Simple charge density model which assumes a linear gradient in charge density in each spatial dimension of a cylindrical coordinate system.
`offsets::NTuple{3, T}` are the charge densities at 0 and `gradients::NTuple{3, T}` are the linear slopes in `r` and `z` direction.
"""
struct CylindricalChargeDensity{T <: SSDFloat} <: AbstractChargeDensity{T}
    offsets::NTuple{3, T}
    gradients::NTuple{3, T}
end

function ChargeDensity(T::DataType, t::Val{:cylindrical}, dict::Union{Dict{String, Any}, Dict{Any, Any}}, inputunit_dict::Dict)
    unit_factor::T = 1
    gradient_unit_factor::T = 1
    if haskey(inputunit_dict, "length")
        lunit = inputunit_dict["length"]
        unit_factor = inv(ustrip(uconvert( internal_length_unit^3, 1 * lunit^3 )))
        gradient_unit_factor = inv(ustrip(uconvert( internal_length_unit^4, 1 * lunit^4 )))
    end
    return CylindricalChargeDensity{T}( dict, unit_factor, gradient_unit_factor )
end

function CylindricalChargeDensity{T}(dict::Union{Dict{String, Any}, Dict{Any, Any}}, unit_factor::T, gradient_unit_factor::T)::CylindricalChargeDensity{T} where {T <: SSDFloat}
    offsets, gradients = zeros(T,3), zeros(T,3)
    if prod(map(k -> k in ["r","z"], collect(keys(dict)))) @warn "Only r and z are supported in the cylindrical charge density model.\nChange the charge density model in the config file or remove all other entries." end
    if haskey(dict, "r")     offsets[1] = geom_round(unit_factor * T(dict["r"]["init"]));     gradients[1] = geom_round(gradient_unit_factor * T(dict["r"]["gradient"]))    end
    if haskey(dict, "z")     offsets[3] = geom_round(unit_factor * T(dict["z"]["init"]));     gradients[3] = geom_round(gradient_unit_factor * T(dict["z"]["gradient"]))    end
    CylindricalChargeDensity{T}( NTuple{3, T}(offsets), NTuple{3, T}(gradients) )
end

function get_charge_density(lcdm::CylindricalChargeDensity{T}, pt::AbstractCoordinatePoint{T})::T where {T <: SSDFloat}
    pt::CylindricalPoint{T} = CylindricalPoint(pt)
    ρ::T = 0
    for i in eachindex(lcdm.offsets)
        ρ += (lcdm.offsets[i] + pt[i] * lcdm.gradients[i]) #* T(1e16) # * T(1e10) * T(1e6) -> 1/cm^3 -> 1/m^3
    end
    return ρ
end

"""
    struct ZeroChargeDensity{T <: SSDFloat} <: AbstractChargeDensity{T}

Returns always 0.
"""
struct ZeroChargeDensity{T <: SSDFloat} <: AbstractChargeDensity{T} end

function get_charge_density(lcdm::ZeroChargeDensity{T}, pt::AbstractCoordinatePoint{T})::T where {T <: SSDFloat}
    ρ::T = 0
    return ρ
end

"""
    struct ConstantChargeDensity{T <: SSDFloat} <: AbstractChargeDensity{T}

Returns always a fixed charge density.
"""
struct ConstantChargeDensity{T <: SSDFloat} <: AbstractChargeDensity{T}
    ρ::T
end

function get_charge_density(cdm::ConstantChargeDensity{T}, pt::AbstractCoordinatePoint{T})::T where {T <: SSDFloat}
    return cdm.ρ
end

function ChargeDensity(T::DataType, t::Val{:constant}, dict::Union{Dict{String, Any}, Dict{Any, Any}}, inputunit_dict::Dict)
    unit_factor::T = 1
    if haskey(inputunit_dict, "length")
        lunit = inputunit_dict["length"]
        unit_factor = inv(ustrip(uconvert( internal_length_unit^3, 1 * lunit^3 )))
    end
    return ConstantChargeDensity{T}( dict, unit_factor )
end
function ConstantChargeDensity{T}(dict::Union{Dict{String, Any}, Dict{Any, Any}}, unit_factor::T)::ConstantChargeDensity{T} where {T <: SSDFloat}
    ρ::T = if haskey(dict, "charge_density")
        geom_round(unit_factor * T(dict["charge_density"]))
    else
        T(0)
    end
    ConstantChargeDensity{T}( ρ )
end


text:
function constructBMRBlist(dict)
    N = length(dict)
    entries = Vector{String}(undef, N)
    molecule_names = Vector{String}(undef, N)

    i = 0
    for (key, value) in dict
        #
        #println("typeof(key) = ", typeof(key))
        #println("typeof(value[1]) = ", typeof(value[1]))

        i += 1
        entries[i] = String(key)
        molecule_names[i] = value[1]
    end

    return entries, molecule_names
end

### fetch the current list of molecules in the BMRB metabolomics database.
function getchBMRBmetabolomicslist(folder_path::String = joinpath(homedir(),"MEGAsync/data/NMR/BMRB"))

    a = "http://api.bmrb.io/current/search/get_all_values_for_tag/Chem_comp.Name?database=metabolomics"
    HTTP_result = HTTP.request("GET", a)
    dict = JSON3.read(HTTP_result.body)
    BMRB_list, list_names = constructBMRBlist(dict)

    # store queried results.
    isdir(folder_path) || mkdir(folder_path)

    save_path = joinpath(folder_path, "list.jld")
    JLD.save(save_path, "BMRB_list", BMRB_list,
                        "list_names", list_names)

    return BMRB_list, list_names
end

function downloadBMRBdata(target_names::Vector{String},
                BMRB_list::Vector{String},
                list_names::Vector{String},
                GISSMO_dir::String,
                BMRB_zip_dir::String,
                BMRB_data_dir::String;
                required_experiment_substring::String = "1D 1H",
                delay_time = 0.2)

    #
    N_molecules = length(target_names)

    entry_set = Vector{Vector{String}}(undef, N_molecules)
    query_success_flags = falses(N_molecules)
    entry_success_flags_set = Vector{BitArray{1}}(undef, N_molecules)

    for i = 1:N_molecules

        display_string = "Processing $(i)/$(N_molecules), $(target_names[i])"
        println(display_string)

        entry_set[i], query_success_flags[i],
            entry_success_flags_set[i] = downloadBMRBdata(target_names[i],
                BMRB_list,
                list_names,
                GISSMO_dir,
                BMRB_zip_dir,
                BMRB_data_dir;
                required_experiment_substring = required_experiment_substring,
                delay_time = delay_time)

        println()
        println()
    end

    return entry_set, query_success_flags, entry_success_flags_set
end

# for one target molecule.
# delay_time in seconds. To respect REST API rate limit.
function downloadBMRBdata(target_name::String,
                BMRB_list::Vector{String},
                list_names::Vector{String},
                GISSMO_dir::String,
                BMRB_zip_dir::String,
                BMRB_data_dir::String;
                required_experiment_substring::String = "1D 1H",
                delay_time = 0.2)

    ### search for the entries that contain the target name.
    inds, target_exist_flag = searchmoleculelist(list_names, target_name)
    search_result_entries = BMRB_list[inds]

    N = length(search_result_entries)
    entry_success_flags = falses(N)

    for n = 1:N

        entry = search_result_entries[n]
        println("Processing entry ", entry)

        # get metadata, so we can figure out which experiments to keep.
        URLs, all_metadata_JSON, inds_1D1H = queryBMRBexperimentdataURL(entry; target_string = required_experiment_substring)
        metadata_JSON = all_metadata_JSON[inds_1D1H]

        # download, extract, keep only 1D 1H entries. unzip_dir not used.
        unzip_dir = downloadandextractBMRBentry(entry, BMRB_zip_dir, BMRB_data_dir; delay = delay_time)
        entry_dir, experiment_paths, data_move_success_flag = keeponly1D1Hexperiments(entry, BMRB_zip_dir, BMRB_data_dir, unzip_dir, URLs)

        if !data_move_success_flag

            # report bad entry, skip.
            entry_success_flags[n] = false

        else
            # save raw metadata to storage.
            save_path = joinpath(entry_dir, "metadata.json")
            open(save_path, "w") do f
                JSON3.pretty(f, JSON3.write(metadata_JSON))
                println(f)
            end

            # cache the 1D 1H experiment directories for this entry.
            save_path = joinpath(entry_dir, "experiment_dirs.jld")
            JLD.save(save_path, "experiment_dirs", experiment_paths,
                                "molecule_name", target_name)

            ### store the sample metadata for each experiment.
            N_experiments = length(metadata_JSON)
            @assert length(experiment_paths) == N_experiments

            for j = 1:N_experiments

                # concentration.
                X = metadata_JSON[j][:Sample_component]
                solute_mM, ref_mM, ref_molecule_name,
                    solubility_gL_solute, molar_mass_solute,
                    solubility_gL_ref, molar_mass_ref = convertBMRBexperimentconcentrationinfo(X, target_name, GISSMO_dir)

                # conditions.
                pH = metadata_JSON[j][:Sample_condition_variable][:ph]
                temperature = metadata_JSON[j][:Sample_condition_variable][:temperature]

                # store.
                # debug.

                @task begin;
                    save_path = joinpath(experiment_paths[j], "metadata.jld")
                    JLD.save(save_path, "solute_mM", solute_mM,
                                    "ref_mM", ref_mM,
                                    "ref_molecule_name", ref_molecule_name,
                                    "solubility_gL_solute", solubility_gL_solute,
                                    "molar_mass_solute", molar_mass_solute,
                                    "solubility_gL_ref", solubility_gL_ref,
                                    "molar_mass_ref", molar_mass_ref,
                                    "pH", pH,
                                    "temperature", temperature)
                    # debug.
                    println("save_path= ", save_path)

                    save_path = joinpath(experiment_paths[j], "metadata.json")
                    open(save_path, "w") do f
                        JSON3.pretty(f, JSON3.write(metadata_JSON[j]))
                        println(f)
                    end
                end
                sleep(0.1)
            end

            entry_success_flags[n] = true
        end

    end

    return search_result_entries, target_exist_flag, entry_success_flags
end

# base_dir is the parent folder where the BMRB entries can be found.
function loadBMRBdata(entries::Vector{String},
                        base_dir::String)

    #
    full_paths = Vector{String}(undef, 0)

    #BMRB_entries = readdir(base_path)

    # loop over each specified entry.
    for i = 1:length(entries)
        entry = entries[i]

        # see if the entry exists on our local storage.
        entry_path = joinpath(base_dir, entry)
        if isdir(entry_path)

            experiment_tags = readdir(entry_path)

            # loop over the experiment for the current entry.
            entry_paths = collect( joinpath(entry_path, experiment_tags[k]) for k = 1:length(experiment_tags) )

            # push.
            push!(full_paths, entry_paths...)
        end
    end

    # remove paths to files.
    keep_flags = collect( !isfile(full_paths[i]) for i = 1:length(full_paths) )
    full_paths = full_paths[keep_flags]

    return full_paths
end

# X = metadata_JSON[i][:Sample_component]
function convertBMRBexperimentconcentrationinfo(X,
            molecule_name::String, GISSMO_dir::String)

    ### concentration
    ref_mM = NaN
    solute_mM = NaN

    ref_ind = findfirst(xx->(occursin("Reference", xx[:Type]) || occursin("reference", xx[:Type])), X)
    if length(ref_ind) < 1
        println("Warning, cannot find reference entry in provided sample composition metadata. Assigning NaN as reference concentration.")
        ref_mM = NaN
    end

    solute_ind = findfirst(xx->(occursin("Solute", xx[:Type]) || occursin("solute", xx[:Type])), X)
    if length(solute_ind) < 1
        println("Warning, cannot find solute entry in provided sample composition metadata. Assigning NaN as solute concentration.")
        solute_mM = NaN
    end

    A = X[solute_ind]

    solubility_gL_solute = solubilityinwater(molecule_name)
    molar_mass_solute = getmolarmass(GISSMO_dir, molecule_name)

    solute_mM = convertstringtomM(A[:Concentration_val],
                A[:Concentration_val_units], solubility_gL_solute, molar_mass_solute)

    #
    B = X[ref_ind]
    ref_molecule_name = B[:Mol_common_name]

    solubility_gL_ref = solubilityinwater(ref_molecule_name)

    ref_entry = getGISSMOentry(ref_molecule_name)
    molar_mass_ref = getmolarmass(GISSMO_dir, ref_molecule_name)

    ref_mM = convertstringtomM(B[:Concentration_val],
                B[:Concentration_val_units], solubility_gL_ref, molar_mass_ref)

    return solute_mM, ref_mM, ref_molecule_name,
        solubility_gL_solute, molar_mass_solute,
        solubility_gL_ref, molar_mass_ref
end

# download to zip dir and unzip to data dir. delay in seconds.
function downloadandextractBMRBentry(entry, BMRB_zip_dir, BMRB_data_dir; delay = 0.2)

    # download entry to zip directory..
    URL = "https://bmrb.io/metabolomics/mol_summary/zip_entry_directory.php?id=$(entry)"
    dest = joinpath(BMRB_zip_dir, "$(entry).zip")

    t = @task begin; isfile(dest) || download(URL, dest); println("done download ", entry); end
    schedule(t); wait(t)
    sleep(delay) # respect REST API query rate limit.

    # unzip to dataset directory.
    unzip_dir = joinpath(BMRB_data_dir, "$(entry)_unzip")
    t = @task begin; isdir(unzip_dir) || mkdir(unzip_dir); end
    schedule(t); wait(t)

    t = @task begin; InfoZIP.unzip(dest, unzip_dir); println("done unzip ", entry); end
    schedule(t); wait(t)

    return unzip_dir
end

# keep only 1D1H experiments from the unzipped directory.
function keeponly1D1Hexperiments(entry,
            BMRB_zip_dir, BMRB_data_dir, unzip_dir, URLs)

    ### delete directories that aren't 1D 1H.
    # move each of the experiment paths.
    entry_dir = joinpath(BMRB_data_dir, "$(entry)")
    t = @task begin; ispath(entry_dir) || mkpath(entry_dir); end
    schedule(t); wait(t)

    experiment_tmp_paths = collect( getlocal1D1Hexperimentpath(URLs[k], entry , unzip_dir) for k = 1:length(URLs) )
    problem_flags = collect( !ispath(experiment_tmp_paths[i]) for i = 1:length(experiment_tmp_paths) )

    experiment_paths = collect( getlocal1D1Hexperimentpath(URLs[k], entry , entry_dir) for k = 1:length(URLs) )

    if any(problem_flags)

        # zip archive and metadata URL not the same.
        # report bad entry. Remove the unzip folder without moving.
        println("Path issue for ", entry, " does not agree with metadata URL. Skip.")

        println("experiment_tmp_paths ", experiment_tmp_paths)
        println("experiment_paths ", experiment_paths)
        println("URLs ", URLs)
        println("problem_flags ", problem_flags)

        #rm(unzip_dir; force = true, recursive = true)

        return entry_dir, experiment_paths, false
    end

    t = @task begin; mvfolders(experiment_tmp_paths, experiment_paths); println("done moving ", entry); end
    schedule(t); wait(t)

    rm(unzip_dir; force = true, recursive = true)

    return entry_dir, experiment_paths, true
end

######### work-around for experiment ID mismatch with download URL issue.

# if we've X[k][:attribute] has target_string as a sub-string, then return k.
# return a vector of all k's that match.
function searchJSON3arrayformatch(X, attribute::Symbol, target_string::String)
    N = length(X)

    inds = Vector{Int}(undef, 0)

    for i = 1:length(X)
        if occursin(target_string, X[i][attribute])
            push!(inds, i)
        end
    end

    return inds
end

# assumes mkpath(dest) will work.
function mvfolders(experiment_tmp_paths, experiment_paths)

    for k = 1:length(experiment_paths)

        src = experiment_tmp_paths[k]
        dest = experiment_paths[k]

        #println("src ", src)
        #println("dest ", dest)

        t = @task begin; ispath(dest) || mkpath(dest); end
        schedule(t); wait(t)

        mv(src, dest; force = true)
    end
end

# get path to the 1D 1H experiment reference in the URL.
function getlocal1D1Hexperimentpath(URL::String,
            entry::String, entry_dir::String)::String

    # keep only the folder names after the entry in the URL.
    q = split(URL, '/')
    l = findfirst(xx->xx==entry, q)

    # assemble local path.
    experiment_path = entry_dir
    for i = l+1:length(q)
        experiment_path = joinpath(experiment_path, q[i])

        # if !ispath(experiment_path)
        #     # error.
        #     return ""
        # end
    end

    return experiment_path
end

# search for 1D 1H for the current entry.
function queryBMRBexperimentdataURL(entry::String;
                        target_string::String = "1D 1H")
    a = "http://api.bmrb.io/v2/entry/$(entry)/experiments"

    HTTP_result = HTTP.request("GET", a)
    experiments = JSON3.read(HTTP_result.body)

    inds = searchJSON3arrayformatch(experiments, :Name, target_string)

    URLs = Vector{String}(undef, 0)
    for k in inds
        tmp = experiments[k][:Experiment_file]
        inds2 = searchJSON3arrayformatch(tmp, :description, "experiment directory")

        for j in inds2
            push!(URLs, tmp[j][:url])
        end
    end

    return URLs, experiments, inds
end


text:
using MagneticLocSuchowiak
using Test

@testset "MagneticLocSuchowiak.jl" begin
    # Write your tests here.
end

max_stars_repo_name: SebastianSuchowiak/MagneticLocSuchowiak

max_issues_repo_name: SebastianSuchowiak/MagneticLocSuchowiak

max_forks_repo_name: SebastianSuchowiak/MagneticLocSuchowiak


text:
# sample from a double well in 1D using HMC
using Plots
using Printf
using Random
using LinearAlgebra
using ForwardDiff
using QuadGK
using BasicMD

include("potentials.jl")

β = 5.0;
x₀ = [-1.0];
seed = 100;
Δt = 1e-2;
n_iters = 10^5; # number of samples

V = x->DoubleWell(x);
cfg = ForwardDiff.GradientConfig(V, x₀);
gradV! = (gradV, x)-> ForwardDiff.gradient!(gradV, V, x, cfg);

sampler = EM(gradV!, β, Δt);

# define the recycling functions
a = [-1.0];
b = 0.9;
function inB(state::BasicMD.EMState)
    return state.x[1]>b
end
function restartA!(state::BasicMD.EMState)
    @. state.x = a;
    gradV!(state.∇V, a);
    state
end

recycler = Recycler(restartA!, inB);

Random.seed!(100);
X₀ = copy(x₀);
sample_trajectory!(X₀, sampler,recycler,  options=MDOptions(n_iters=n_iters));
@printf("In Place X after %d iterations: %g\n",n_iters, X₀[1])

Random.seed!(100);
Xvals = sample_trajectory(x₀, sampler,recycler, options=MDOptions(n_iters=n_iters,n_save_iters=n_iters));
X = Xvals[end];
@printf("X after %d iterations: %g\n",n_iters, X[1])

Random.seed!(100);
X_vals = sample_trajectory(x₀, sampler,recycler, options=MDOptions(n_iters=n_iters));
histogram([X[1] for X in X_vals],label="Samples",normalize=true)
xlabel!("x")
ylabel!("Frequency")

max_stars_repo_name: gideonsimpson/JuBasicMD

max_issues_repo_name: gideonsimpson/JuBasicMD

max_forks_repo_name: gideonsimpson/JuBasicMD


text:
# =================================================================
# Lotka-Volterra with nonlinear guard crossing model
# See https://easychair.org/publications/paper/nrdD
# =================================================================

using ReachabilityAnalysis, Plots

const T_lv = 3.64
const r₀ = 0.161
const ε  = 0.012
const X0 = IntervalBox( 1.3 ± ε, 1 .. 1)
const SB2 = symmetric_box(2, Float64)

@taylorize function lotka_volterra!(du, u, p, t)
    x, y = u
    xy = x * y
    du[1] = 3.0 * (x - xy)
    du[2] = xy - y
    return du
end

circle(x,y) = (x-1)^2 + (y-1)^2

function conditionNotOutside_Up(x,y; r₀=r₀)
    r₀² = r₀^2
    if y > 1.1
        return !(circle(x,y) > r₀²)
    else
        return false
    end
end
function conditionNotOutside_Down(x,y; r₀=r₀)
    r₀² = r₀^2
    if y < 0.9
        return !(circle(x,y) > r₀²)
    else
        return false
    end
end

lotka_volterra(Xini) = @ivp(x'= lotka_volterra!(x), dim:2, x(0) ∈ Xini)

@inline function lv_property(solz, r₀)
    # Time, time domain, x- and y-interval of the (minced) solutions
    timeI = getindex.(getfield.(solz,:Δt),1);
    domTs = domain.(getindex.(getfield.(solz,:X),1));
    xs = evaluate.(evaluate.(getindex.(getfield.(solz, :X),1), domTs), (SB2,));
    ys = evaluate.(evaluate.(getindex.(getfield.(solz, :X),2), domTs), (SB2,));

    # Get the indexes of the sets that are not outside the circle, i.e.,
    # that intersect the nonlinear guard
    rng_Up   = findall(conditionNotOutside_Up.(xs, ys, r₀=r₀));
    rng_Down = findall(conditionNotOutside_Down.(xs, ys, r₀=r₀));

    # Compute time intervals where the boxes are not outside the non-linear guard,
    # and the area of the region not outside the circle
    area = 0.0
    hull1 = emptyinterval()
    hull2 = emptyinterval()
    ib1 = [IntervalBox(0..0, 2)]
    ib2 = [IntervalBox(0..0, 2)]
    if !isempty(rng_Up)
        hull1 = hull(timeI[rng_Up]...)
        area = sum(diam.(xs[rng_Up]) .* diam.(ys[rng_Up]))
        ib1 = IntervalBox.(xs[rng_Up[1]:end], ys[rng_Up[1]:end])
    end
    if !isempty(rng_Down)
        hull2 = hull(timeI[rng_Down]...)
        area += sum(diam.(xs[rng_Down]) .* diam.(ys[rng_Down]))
        ib2 = IntervalBox.(xs[rng_Down[1]:end], ys[rng_Down[1]:end])
    end

    return (hull1, hull2, area, ib1, ib2)
end

# function lotka_volterra(; nsplit=4,
#                           ε = 0.012,
#                           ε_ext=1e-4, # threshold for the outer approximation
#                           n_int=50,   # number of directions for the inner approximation
#                           r₀ = 0.161) # radius circle

#     # generate external / internal polytopic approximations of the guard
#     B = Ball2([1.0, 1.0], r₀) # "exact"
#     B_ext = overapproximate(B, ε_ext) # outer approximation
#     B_int = underapproximate(B, PolarDirections(n_int)) # inner approximation
#     B_int = tohrep(convert(VPolygon, B_int)) # cast to Hrep
#     B_intᶜ = complement(B_int)

#     # define modes
#     aut = LightAutomaton(3)
#     outside = @system(x' = lotka_volterra!(x), dim: 2, x ∈ B_intᶜ)
#     inside = @system(x' = lotka_volterra!(x), dim: 2, x ∈ B_ext)
#     outside_unconstrained = @system(x' = lotka_volterra!(x), dim: 2, x ∈ Universe(2))

#     # define the transition graph
#     add_transition!(aut, 1, 2, 1)
#     add_transition!(aut, 2, 3, 2)
#     T_out_in = @map(x -> x, dim:2, x ∈ B_ext)
#     T_in_out = @map(x -> x, dim:2, x ∈ B_intᶜ)

#     # initial-value problem
#     H = HybridSystem(automaton=aut, modes=[outside, inside, outside_unconstrained],
#                                            resetmaps=[T_out_in, T_in_out])

#     # initial states with splitting
#     X0 = Hyperrectangle(low=[1.3-ε, 1.], high=[1.3+ε, 1.])
#     X0s = split(X0, [nsplit, 1])
#     X0st = [(X0s_i, 1) for X0s_i in X0s]

#     return InitialValueProblem(H, X0st)
# end

# @inline function lv_property(solz, ε_ext, r₀)

#     # Sets intersecting the nonlinear guard
#     B = Ball2([1.0, 1.0], r₀) # "exact"
#     B_ext = overapproximate(B, ε_ext) # outer approximation
#     intersecting_reachsets = []
#     for (i, Fi) in enumerate(solz)
#         for (j, Xj) in enumerate(Fi)
# 	        Xjp = convert(VPolygon, Xj)
#             !is_intersection_empty(Xjp, B_ext) && push!(intersecting_reachsets, (i, j))
#         end
#     end

#     # Compute time spent inside non-linear guard
#     times = [tspan(solz[ind[1]][ind[2]]) for ind in intersecting_reachsets]
#     tmin = minimum(tstart, times)
#     tmax = maximum(tend, times)
#     @show(tmin, tmax, tmax-tmin)

#     indxs = Int[]
#     for (i, e) in enumerate(tspan.(solz))
#         T_lv ∈ tspan(e) && push!(indxs, i)
#     end
#     chsolz = [solz[i](T_lv) |> convexify for i in indxs]
#     chlasth = ConvexHullArray([set(Ri) for Ri in chsolz]) |> box_approximation
#     return LazySets.area(chlasth), tmax-tmin
# end


text:
using DocStringExtensions: SIGNATURES

function importdata_nav(frames::Array{Int,1};
                        datadir::String="/media/data1/data/kayaks/20_gps_pos")

    # "/media/data1/data/kayaks/20_parsed" this data is windowed 10 elements referenced to the first
    # get position data
    #  logdir = wavedatadir*"/array_positions$(sasframe).csv"  # 85

    dataOut = zeros(length(frames),2);
    count = 0
    for i in frames
        count+=1
        navfile = datadir*"/nav$(i).csv"
        dataOut[count,:] = readdlm(navfile,',',Float64,'\n')
    end
    dataOut
end

function sanitycheck_nav(navdata)
    check = true
    errorout = 0;
    for i in 1:size(navdata,1)-1
        dist = abs.(navdata[i,:]-navdata[i+1,:]);
        if dist[1] > 3
            println("NAV_X is off by $(dist[1])")
            check = false
            errorout = i
        end
        if dist[2] > 3
            println("NAV_Y is off by $(dist[2])")
            check = false
            errorout = i
        end
    end
    check, errorout
end

"""
    $(SIGNATURES)

Load hydrophone waveform data from hard drive, after it has been processed from MOOS [a/b]log.  Processing from MOOS format can be done according to scripts/parseData.py -- Mei says so, call on @mc2922 if you can't find it...

Input:
- `frames` is an index into the time of the mission, i.e. an array of poses that you are calling (meaningless index unless 1Hz happens to be time).
- `element` describes which hydrophone to use on the physcal array (initial development had 5 elements and commonly using 2).
- `datadir` is the directory location where all the waveforms%d.csv exist.

Output:
- an array of waveforms where the first dimension is the length of the waveform in time (i.e. number of samples), and the second dimension is the number of frames you ask for.

Example:

```julia
frm = [1,2,3,4] # corresponds to files waveforms[1/2/3/4].csv
data = importdata_waveforms(frm, 2)
```

Future:
There might be a more memory 'in-place' implementation, which would change the API.
"""
function importdata_waveforms(frames::Array{Int,1},
                              element::Int;
                              datadir::String="/media/data1/data/kayaks/20_gps_pos")
    wavefile = datadir*"/waveform$(frames[1]).csv"
    tempData = readdlm(wavefile,',',Float64,'\n');
    dataOut = zeros(size(tempData,2),length(frames));
    count = 0
    for i in frames
        count+=1;
        wavefile = datadir*"/waveform$(i).csv"
        tempData = readdlm(wavefile,',',Float64,'\n');
        dataOut[:,count] = tempData[element,:];
    end
    dataOut
end

function cumulativeDrift!(cdata::AbstractArray,cstart::Array,σ::Array)
    cdrift = cstart
    for i in 1:size(cdata,1)
        cdata[i,:] += cdrift
        cdrift += σ.*randn(2)
    end
    nothing
end

function addsascluster_velpt!(fg::G,
                              poses::Vector{Symbol},
                              beacon::Symbol,
                              sasframes::Array{Int},
                              tstart;
                              rtkCov::Array{Float64,2}= diagm([0.1;0.1].^2),
                              element::Int=2,
                              autoinit::Bool=false,
                              datadir::String="/media/data1/data/kayaks/20_gps_pos",
                              N::Int=100  ) where G <: AbstractDFG

    if length(sasframes) != length(poses)
        error("Frames and Poses are different lengths")
    end

  numelems = length(poses)

  lss = union(ls(fg)...)

  posData = importdata_nav(sasframes,datadir=datadir);
  dposData = deepcopy(posData)
  cumulativeDrift!(dposData,[0.0;0],[0.05,0.05])
  waveformData = importdata_waveforms(sasframes,element,datadir=datadir);
  tcurrent = tstart*1_000_000

  for sym in poses
    addVariable!(fg, sym, DynPoint2(ut=tcurrent))
    tcurrent += 1_000_000
  end

  #Fix the first point
  #Interpolate velocity from GPS ground truth
  xdotp = posData[2,1] - posData[1,1];
  ydotp = posData[2,2] - posData[1,2];
  dpμ = [posData[1,1];posData[1,2];xdotp;ydotp];
  dpσ = diagm([0.1;0.1;0.1;0.1].^2)

  pp = DynPoint2VelocityPrior(MvNormal(dpμ,dpσ))
  addFactor!(fg, [poses[1];], pp, autoinit=true)
  IncrementalInference.doautoinit!(fg,[getVariable(fg,poses[1])])
  setValKDE!(getVariable(fg,poses[1]),kde!(getSample(getfnctype(getFactor(fg, Symbol("$(poses[1])f1"))),N)[1]))

  #Add odo factors
  for i in 2:numelems
      @show xdotp = dposData[i,1] - dposData[i-1,1]
      @show ydotp = dposData[i,2] - dposData[i-1,2]
      dpμ = [xdotp;ydotp;0;0];
      dpσ = diagm([0.5;0.5;0.1;0.1].^2)
      vp = VelPoint2VelPoint2(MvNormal(dpμ,dpσ))
      addFactor!(fg, [poses[i-1];poses[i]], vp, autoinit=true)
      IncrementalInference.doautoinit!(fg,[getVariable(fg,poses[i])])
      #factorsym = Symbol("$(poses[i-1])$(poses[i])f1")
      #setValKDE!(getVariable(fg,poses[i]),kde!(getSample(getfnctype(getVariable(fg, factorsym,nt=:fnc)),N)[1]))
      @show hackinit = ones(4,N).*[dposData[i-1,:];xdotp;ydotp] + rand(MvNormal(dpμ,dpσ),N)
      setVal!(getVariable(fg,poses[i]),hackinit)
  end

  sas2d = prepareSAS2DFactor(numelems, waveformData, rangemodel=:Correlator)
  addFactor!(fg, [beacon;poses], sas2d, autoinit=autoinit)

  dposData
end

#

function addsascluster_only_gps!(fg::G,
                                 poses::Vector{Symbol},
                                 beacon::Symbol,
                                 sasframes::Array{Int};
                                 rtkCov::Array{Float64,2}= diagm([0.1;0.1].^2),
                                 datadir::String="/media/data1/data/kayaks/20_gps_pos",
                                 element::Int=2,
                                 autoinit::Bool=false,
                                 N::Int=100  ) where G <: AbstractDFG

  # "/media/data1/data/kayaks/20_parsed" this data is windowed 10 elements referenced to the first
  # get position data
#  logdir = wavedatadir*"/array_positions$(sasframe).csv"  # 85

    if length(sasframes) != length(poses)
        error("Frames and Poses are different lengths")
    end

  numelems = length(poses)

  lss = union(ls(fg)...)

  posData = importdata_nav(sasframes, datadir=datadir);
  waveformData = importdata_waveforms(sasframes,element,datadir=datadir);
  count = 0

  for sym in poses
    count += 1
    addVariable!(fg, sym, Point2)
    pp = PriorPoint2(MvNormal(posData[count,:], rtkCov) )
    addFactor!(fg, [sym;], pp, autoinit=true)
    IncrementalInference.doautoinit!(fg,sym)
    setValKDE!(getVariable(fg,sym),kde!(getSample(getfnctype(getFactor(fg, Symbol("$(sym)f1"))),N)[1]))
  end

  sas2d = prepareSAS2DFactor(numelems, waveformData, rangemodel=:Correlator)
  addFactor!(fg, [beacon;poses], sas2d, autoinit=autoinit)

  nothing
end

function approxConvFwdBFRaw(fg,
                            poses,
                            beacon,
                            origin,
                            scale,
                            snrfloor;
                            N::Int=100  )
  #
  fctsym = Symbol(string(beacon,poses...,:f1))
  sas2d = IncrementalInference.getfnctype(getFactor(fg, fctsym))
  sas2d.debugging = true
  for ti in 1:Threads.nthreads() reset!(sas2d.threadreuse[ti].dbg) end
  predL1 = IncrementalInference.approxConv(fg, fctsym, beacon, N=N)

  allbeams = []
  for ti in 1:Threads.nthreads()
      for itr in eachindex(sas2d.threadreuse[ti].dbg.beams)
          push!(allbeams,sas2d.threadreuse[ti].dbg.beams[itr])
      end
  end
  @show size(allbeams[1])
  avg_beam = zeros(length(allbeams[1]))
  for i in 1:N
    avg_beam[:] = avg_beam + allbeams[i]
  end

  #avg_beam ./= sum(avg_beam)
  #avg_beam .-= quantile(avg_beam, snrfloor)
  #avg_beam[avg_beam .< 0.0] = 0.0

  circ_avg_beam = zeros(2,length(avg_beam))
  count = 0
  for th in range(0,2pi,length=length(avg_beam))
    count += 1
    circ_avg_beam[:,count] = R(th)*scale*[avg_beam[count];0.0]+origin
  end

  allbeams
end

function approxConvFwdBFlayer(fg,
                              poses,
                              beacon,
                              origin,
                              scale,
                              snrfloor::Float64=0.0;
                              N::Int=100  )
  #
  circ_avg_beam = approxConvFwdBFRaw(fg, poses, beacon, origin, scale, snrfloor, N=N)

  plPica = Gadfly.layer(x=circ_avg_beam[1,:], y=circ_avg_beam[2,:], Geom.path(),   Theme(default_color=colorant"magenta", line_width=2pt))

  plPica
end

"""
    $SIGNATURES

Plot side by side pair of poses and beacon for a particular SAS factor.

Notes:
- Specify non-default `filepath<:String`
- can show pdf via evince if `show::Bool`
"""
function plotSASPair(fg::G,
                     fsym::Symbol;
                     show::Bool=true,
                     filepath::AS="/tmp/test.pdf") where {G <: AbstractDFG, AS <: AbstractString}
    #
    PL = [];

    @show fsym
    @show vars = lsf(fg, fsym)
    beac = vars[1]

    for var in vars[2:end]
        X1 = getVal(getVariable(fg, var))
        push!(PL, layer(x=X1[1,:],y=X1[2,:], Geom.histogram2d))
    end

    # push!(PL, approxConvFwdBFlayer(fg, poses, :l1, posData[1,:], 2000 ))

    pl11 = Gadfly.plot(PL...);
    # pl.coord = Coord.Cartesian(xmin=10,xmax=30,ymin=-60,ymax=-45)

    pl21 = plotKDE(fg, vars[2:end], levels=1, dims=[1;2])

    L1v = getVariable(fg, beac)
    L1 = getVal(L1v)
    pl12 = Gadfly.plot(x=L1[1,:],y=L1[2,:], Geom.histogram2d)

    pl22 = plotKDE(getKDE(L1v), levels=3)

    # stack all images together
    plT = Gadfly.hstack(pl11,pl12)
    plB = Gadfly.hstack(pl21,pl22)
    plA = Gadfly.vstack(plT, plB)

    # if the user wants to see a file rendering
    if show
      plA |> PDF(filepath);  @async run(`evince $filepath`)
    end

    plA
end


text:
using ZMQ

server = Socket(REP)
client = Socket(REQ)
addr = "tcp://127.0.0.1:5555"

bind(server, addr)
connect(client, addr)

send(client, "request")
println("<-", "request")
msg = recv(server, String)
println("->", msg)
send(server, "response")
println("<-", "response")
close(server)
close(client)

max_stars_repo_name: cdluminate/MyNotes

max_issues_repo_name: cdluminate/MyNotes

max_forks_repo_name: cdluminate/MyNotes


text:
using PkgBenchmark, MAVs, NLOptControl

@benchgroup "psMethods"  ["integrationScheme","Nck"] begin
    for scheme in (:lgrExplicit,:lgrImplicit)
        for Nck in ([10,8,6],[12,10,8,6])
          c=defineCase(;(:mode=>:autoBench));
          setMisc!(c;integrationScheme=scheme,Nck=Nck,max_cpu_time=20.0)
          n=initializeAutonomousControl(c);
          @bench string(scheme), Nck  optimize!($n)
        end
    end
end


text:
#!/usr/bin/julia --color=yes
using ProgressMeter

    noflags = filter(x->'-' ∉x, ARGS)
if length(noflags) == 1
    info("Will use $(noflags[1]) as L.")
    const L        = parse(Int64,noflags[1])
else
    const L = 8
    info("Defaulting to L=$L")
end

function run_annealer(iters; T=1)
    cmd = `./exe/annealer -l $L -i $iters -m 2 -n 2 1 /tmp/debuganneal`
    val, time = @timed run(cmd)
    return time
end

if "-f" ∉ ARGS
    tt = linspace(1e2,1e5,10) |> floor
    times = @showprogress "Benchmarking..." [run_annealer(t) for t in tt]
else
    info("Skipping measurement.")
    tt = linspace(1e2,1e5,10) |> floor
    times = 0*tt
end

function fit_to_line(x, y, yerr=ones(x))
    L = length(x)

    br(x,σ) = sum(  x./(σ.^2) / length(x) )
    brone = br(ones(L), yerr)
    brx   = br(x      , yerr)
    brxx  = br(x.*x   , yerr)
    bry   = br(y      , yerr)
    brxy  = br(x.*y   , yerr)

    m_mean = ( brone*brxy - brx*bry )/
            ( brone*brxx - brx*brx )
    n_mean = mean(y) - m_mean * mean(x)

    covmatrix = zeros(2,2)

    D = brxx*brone - brx*brx

    covmatrix[2,2] = +1/(L*D) * brxx
    covmatrix[1,2] = -1/(L*D) * brx
    covmatrix[2,1] = -1/(L*D) * brx
    covmatrix[1,1] = +1/(L*D) * brone

    σ = diag(covmatrix) |> sqrt

    return m_mean, σ[1], n_mean, σ[2]
end

μm, σm, μn, σn = fit_to_line(tt,times)

@printf("t ∼ %.2lf ± %.2lf × mcsteps × L³  ns\n",
        1e9*μm/L^3, 1e9*σm/L^3)

println()

# Based in empirical data:
μm_memento30 = 126.99*L^3 / 1e9
μm_memento100 = 281.37*L^3 / 1e9

println("Approximate iterations (L=$L):")

print_with_color(:bold,"                                           Memento  \n")
print_with_color(:bold,"        Time         This PC        30 jobs        100 jobs\n")
@printf("    In 1 minute      %.1e        %.1e        %.1e\n" , 60/μm        , 60/μm_memento30        , 60/μm_memento100        )
@printf("    In 1 hour        %.1e        %.1e        %.1e\n" , 3600/μm      , 3600/μm_memento30      , 3600/μm_memento100      )
@printf("    In 1 day         %.1e        %.1e        %.1e\n" , 24*3600/μm   , 24*3600/μm_memento30   , 24*3600/μm_memento100   )
@printf("    In 3 days        %.1e        %.1e        %.1e\n" , 3*24*3600/μm , 3*24*3600/μm_memento30 , 3*24*3600/μm_memento100 )
@printf("    In 1 week        %.1e        %.1e        %.1e\n" , 7*24*3600/μm , 7*24*3600/μm_memento30 , 7*24*3600/μm_memento100 )

max_stars_repo_name: RedPointyJackson/tfg

max_issues_repo_name: RedPointyJackson/tfg

max_forks_repo_name: RedPointyJackson/tfg