text:
###############################################
#                                             #
# Four state model:                           #
# 1. Solve                                    #
# 2. Generate output for figures and tables   #
#                                             #
###############################################

"""
    withpolicy(policy; TT=52, ndims=20, thetai=0.38, thetas=0.0044)
    withpolicy(policy, nprocs; TT=52, ndims=20, thetai=0.38, thetas=0.0044)

Solve the model under a user-specified government policy.
In all cases, the government quarantines individuals that are known to be infected.
The output is a [`Dict`](https://docs.julialang.org/en/v1/base/collections/#Dictionaries) containing the necessary variables and time-series for generating the figures.

# Arguments:
- `policy::String`: Label for the government policy. Takes values
    - `notest`: no testing at all.
    - `untargettest`: random testing.
    - `targettest`: targeted testing under access to contact-tracing technology.
    - `isolate`: no testing, but isolation of certain fraction of population.
- `nprocs::Int64`: Optional. If specified, uses the module [Distributed](https://docs.julialang.org/en/v1/manual/distributed-computing/) for parallelizing the computations. *Requires prior setup of `nprocs` workers*.
- `TT::Integer`: Number of periods.
- `ndims::Integer`: Number of grid points in each dimension.
- `algorithm`: [NLopt](https://github.com/JuliaOpt/NLopt.jl) optimization algorithm.
- `thetai::Float`: Probability that infected person sends "infected" signal.
- `thetas::Float`: Probability that non-infected person sends "infected" signal.
"""
function withpolicy(policy::String; TT=52, ndims=20, thetai=0.38, thetas=0.0044)

    println("Solving model for policy = '$policy'. Counting iterations ($TT in total)...")

    # Par.ISO = 0, Par.TEST = 0: Model with no testing and isolation
    # Par.ISO = 0, Par.TEST = 1: Model with testing
    # Par.ISO = 1, Par.TEST = 0: Model with isolation and no testing

    ISO, TEST, theti, thets = [0, 0, 0, 0]; # may replace in next step
    algorithm = NLopt.LN_BOBYQA;            # may replace in next step

    if policy == "notest"
        ISO = 0;
        TEST = 0;
        theti = 1;
        thets = 1;
        algorithm = NLopt.LN_BOBYQA;
    elseif policy == "untargettest"
        ISO = 0;
        TEST = 1;
        theti = 1;
        thets = 1;
        algorithm = NLopt.LD_SLSQP;
    elseif policy =="targettest"
        ISO = 0;
        TEST = 1;
        theti = thetai;
        thets = thetas;
        algorithm = NLopt.LN_BOBYQA;
    elseif policy == "isolate"
        ISO = 1;
        TEST = 0;
        theti = thetai;
        thets = thetas;
        algorithm = NLopt.LD_SLSQP;
    else
        println("Invalid policy.")
    end

    # Parameters
    OPT   = 0;    # 0 for fmincon, 1 for global, 2 for interp
    beta  = .99;  # Discount factor
    gamma = 7/18; # Exit rate
    delta = .005; # Death rate
    sigma = 2;    # CES parameter on production function
    B1    = 1;    # Productivity of work sector
    B2    = .1;   # Productivity of home sector
    nu    = 1;    # Testing cost parameter
    c     = (1/.01)*50*52/63000 .* YY_noep(B1, B2, sigma, 1); # Testing cost parameter
    z     = 7500*52/63000 .* YY_noep(B1, B2, sigma, 1);       # Treatment cost
    d     = 15*3*52*(1/.96^(1/52)-1);                         # Death cost
    xi    = .8;           # Infected productivity
    tautilde = 13/18*1/2; # Prob. of becoming symptomatic
    pi2   = .01;          # Infection rate at home
    pi1   = 1.5;          # Infection rate at work
    ut    = 1;
    # Combine into (immutable) NamedTuple
    Par = (ISO=ISO, TEST=TEST, OPT=OPT, beta=beta, gamma=gamma, delta=delta, sigma=sigma, B1=B1, B2=B2,
    nu=nu, c=c, z=z, d=d, xi=xi, tautilde=tautilde, pi2=pi2, pi1=pi1, thets=thets, theti=theti, ut=ut)

    # Grids
    nSS = ndims;
    nII = ndims;
    nRR = ndims;
    Icurve = -2;

    GridS = range(0, 1, length=nSS); #Susceptible
    GridR = range(0, 1, length=nRR); #Recovered
    GridI = 10 .^(range(Icurve, 0, length=nII)); #Infected but not known to be so
    GridI = GridI .- GridI[1];
    GridI[nII] = 1;

    Grid = (nS=nSS, nI=nII, nR=nRR, Icurve=Icurve, S=GridS, R=GridR, I=GridI);

    ## Solve model via backward induction
    rets = zeros(Int8, (TT,nSS,nII,nRR));
    ret_value = Dict(:SUCCESS => 1, :STOPVAL_REACHED => 2, :FTOL_REACHED => 3, :XTOL_REACHED => 4, :MAXEVAL_REACHED => 5, :MAXTIME_REACHED => 6,
    :FAILURE => -1, :INVALID_ARGS => -2, :OUT_OF_MEMORY => -3, :ROUNDOFF_LIMITED => -4, :FORCED_STOP => -5);
    ret_message = Dict("1" => :SUCCESS, "2" => :STOPVAL_REACHED, "3" => :FTOL_REACHED, "4" => :XTOL_REACHED, "5" => :MAXEVAL_REACHED, "6" => :MAXTIME_REACHED,
    "-1" => :FAILURE, "-2" => :INVALID_ARGS, "-3" => :OUT_OF_MEMORY, "-4" => :ROUNDOFF_LIMITED, "-5" => :FORCED_STOP);

    # Empty matrices for results
    V0 = zeros(TT,nSS,nII,nRR); #value function
    qq = zeros(TT,nSS,nII,nRR); #mass tested/isolated
    ll = zeros(TT,nSS,nII,nRR); #mass of unknown types assigned to work sector
    mm = zeros(TT,nSS,nII,nRR); #mass of known types assigned to work sector
    V00 = zeros(nSS,nII,nRR);
    qqg = zeros(nSS,nII,nRR);
    llg = zeros(nSS,nII,nRR);
    mmg = zeros(nSS,nII,nRR);
    V00T = zeros(nSS,nII,nRR);
    x0 = [0.0, 0.0, 0.0];

    # Value in last period
    for i in 1:nSS
        for j in 1:nII
            for k in 1:nRR
                S = Grid.S[i];
                I = Grid.I[j];
                R = Grid.R[k];
                V00T[i,j,k] = VTerm(Par,S,I,R);
            end
        end
    end

    # Backward induction
    for i in 1:TT
        println(i)
        if i>1
            V00[:,:,:] = V0[i-1,:,:,:];
            qqg[:,:,:] = qq[i-1,:,:,:];
            llg[:,:,:] = ll[i-1,:,:,:];
            mmg[:,:,:] = qq[i-1,:,:,:];
        else
            V00[:,:,:] = V00T;
        end
        for j in 1:nSS
            for k in 1:nII
                for l in 1:nRR
                    S = Grid.S[j];
                    I = Grid.I[k];
                    R = Grid.R[l];
                    if S+I>1.2 || S+I+R>1.2
                        V0[i,j,k,l] = NaN;
                        qq[i,j,k,l] = NaN;
                        ll[i,j,k,l] = NaN;
                        mm[i,j,k,l] = NaN;
                    else
                        lb = [0,0,0]; # [tau, lambda, mu]
                        if Par.ISO == 0 && Par.TEST == 0
                            ub = [0,1,1]; # tau=0
                        else
                            ub = [1,1,1];
                        end
                        if i>1
                            x0 = [qqg[j,k,l],llg[j,k,l],mmg[j,k,l]];
                        else
                            x0 = [0,.9,.9];
                        end
                        opt = Opt(algorithm, 3);
                        # Local, derivative-based, finite differences
                        if algorithm == NLopt.LD_SLSQP
                            opt.min_objective = (x,g) -> Bellman_NLopt(Par, Grid, S, I, R, V00, x, g);
                        # Local, derivative-free
                        else
                            opt.min_objective = (x,g) -> Bellman(Par, Grid, S, I, R, V00, x);
                        end
                        opt.lower_bounds = lb;
                        opt.upper_bounds = ub;
                        opt.xtol_abs = 1.0e-8 .* ones(3);
                        opt.maxeval=2500;
                        (optf, optx, ret) = optimize(opt, x0)
                        # If local optimizer did not converge, perform global search
                        #if ret==:FORCED_STOP
                        #    # Local optimizer
                        #    local_opt = Opt(NLopt.LD_SLSQP, 3);
                        #    local_opt.xtol_abs = 1.0e-8 .* ones(3);
                        #    local_opt.maxeval=2500;
                        #    # Global optimizer
                        #    opt = Opt(NLopt.AUGLAG, 3);
                        #    opt.min_objective = (x,g) -> Bellman_NLopt(Par, Grid, S, I, R, V00, x, g);
                        #    #opt.min_objective = (x,g) -> Bellman(Par, Grid, S, I, R, V00, x);
                        #    opt.local_optimizer = local_opt;
                        #    opt.lower_bounds = lb;
                        #    opt.upper_bounds = ub;
                        #    opt.xtol_abs = 1.0e-8 .* ones(3);
                        #    opt.maxeval=2500;
                        #    #opt.population = 20;
                        #    (optf, optx, ret) = optimize(opt, x0)
                        #end
                        rets[i,j,k,l] = ret_value[ret];
                        qq[i,j,k,l] = optx[1];
                        ll[i,j,k,l] = optx[2];
                        mm[i,j,k,l] = optx[3];
                        V0[i,j,k,l] = -Bellman(Par,Grid,S,I,R,V00,optx);
                    end
                end
            end
        end
    end

    # Simulate optimal path
    I0 = 0.02;
    S0 = 1-I0/(1-Par.tautilde);
    R0  = 0;
    IT = zeros(TT);
    ST = zeros(TT);
    RT = zeros(TT);
    qT = zeros(TT);
    lamT = zeros(TT);
    muT = zeros(TT);
    Del = zeros(TT);
    VT = zeros(TT);
    YT = zeros(TT);
    L1 = zeros(TT);
    L2 = zeros(TT);
    C1T = zeros(TT);
    C2T = zeros(TT);
    Rknot = zeros(TT);
    #Test = zeros(TT);
    masstest = zeros(TT);
    mctest = zeros(TT);
    ctest = zeros(TT);
    IT[1] = I0;
    ST[1] = S0;
    RT[1] = R0;
    Del = 0;
    for i = 1:TT

        S = ST[i];
        I = IT[i];
        R = RT[i];
        IN = 1-S-R;

        if i<TT
            V00[:,:,:] = V0[TT-i,:,:,:];
        else
            V00[:,:,:] = V00T;
        end
        lb = [0,0,0];
        if Par.ISO == 0 && Par.TEST == 0
            ub = [0,1,1];
        else
            ub = [1,1,1];
        end
        if i>2
            x0 = [qT[i-1],lamT[i-1],muT[i-1]];
        else
            if Par.ISO == 0 && Par.TEST == 0
                x0 = [0,1,1];
            else
                x0 = [.3,.3,.3];
            end
        end
        opt = Opt(NLopt.LN_BOBYQA, 3);
        opt.min_objective = (x,g) -> Bellman(Par, Grid, S, I, R, V00, x);
        opt.lower_bounds = lb;
        opt.upper_bounds = ub;
        opt.xtol_abs = 1.0e-8 .* ones(3);
        opt.maxeval=2500;
        (optf, optx, ret) = optimize(opt, x0)
        qT[i] = optx[1];
        lamT[i] = optx[2];
        muT[i] = optx[3];
        VT[i] = -Bellman(Par,Grid,S,I,R,V00,optx);
        if Par.TEST == 1
            Del = qT[i].*Par.theti;
        else
            Del = 0;
        end
        if Par.ISO == 1
            LI1 = lamT[i].*((1-qT[i].*Par.thets).*S+(1-qT[i].*Par.theti).*I)+muT[i].*(1-Par.delta).*R;
            LI2 = (1-lamT[i]).*((1-qT[i].*Par.thets).*S+(1-qT[i].*Par.theti).*I)+(1-muT[i]).*(1-Par.delta).*R;
            if LI1>0
                Ipr1 = (lamT[i].^2 .*Par.pi1.*(1-qT[i].*Par.theti).*I)./LI1;
                RIpr1 = (lamT[i].^2 .*Par.pi1.*(1-qT[i].*Par.theti))./LI1;
            else
                Ipr1 = 0;
                RIpr1 = 0;
            end
            if LI2>0
                Ipr2 = ((1-lamT[i]).^2 .*Par.pi2.*(1-qT[i].*Par.theti).*I)./LI2;
                RIpr2 = ((1-lamT[i]).^2 .*Par.pi2.*(1-qT[i].*Par.theti))./LI2;
            else
                Ipr2 = 0;
            end
            Ipr = Ipr1 + Ipr2;
            Ipr = min(Ipr,1);
            if i<TT
                IT[i+1] = min(max((1-Par.tautilde).*(1-Par.gamma).*I+(1-qT[i]*Par.thets).*S.*Ipr,0),1);
                ST[i+1] = min(max((1-(1-qT[i].*Par.thets).*Ipr).*S,0),1);
                RT[i+1] = max(min(R + Par.gamma.*IN,1),0);
            end
        else
            LI1 = lamT[i].*(S+(1-Del).*I)+muT[i].*(1-Par.delta).*R;
            LI2 = (1-lamT[i]).*(S+(1-Del).*I)+(1-muT[i]).*(1-Par.delta).*R;
            if LI1>0
                Ipr1 = (lamT[i].^2 .*Par.pi1.*(1-Del).*I)./LI1;
                RIpr1 = (lamT[i].^2 .*Par.pi1)./LI1;
            else
                Ipr1 = 0;
                RIpr1 = 0;
            end
            if LI2>0
                Ipr2 = ((1-lamT[i]).^2 .*Par.pi2.*(1-Del).*I)./LI2;
                RIpr2 = ((1-lamT[i]).^2 .*Par.pi2)./LI2;
            else
                Ipr2 = 0;
                RIpr2 = 0;
            end
            Ipr = Ipr1 + Ipr2;
            Ipr = min(Ipr,1);
            if i<TT
                IT[i+1] = min(max((1-Del).*(1-Par.tautilde).*(1-Par.gamma).*I+S.*Ipr,0),1);
                ST[i+1] = min(max((1-Ipr).*S,0),1);
                RT[i+1] = max(min(R + Par.gamma.*IN,1),0);
            end
        end

        YT[i] = YY(Par,S,I,R,qT[i],lamT[i],muT[i]);
        C1T[i], C2T[i] = CC(Par,S,I,R,qT[i],lamT[i],muT[i]);
        if Par.ISO == 0
            L1[i] = (lamT[i].*(S+Par.xi.*(1-Del).*I)+muT[i].*(1-Par.delta).*R);
            L2[i] = ((1-lamT[i]).*(S+Par.xi.*(1-Del).*I)+(1-muT[i]).*(1-Par.delta).*R);
        else
            L1[i] = (lamT[i].*(S.*(1-qT[i].*Par.thets)+Par.xi.*(1-qT[i].*Par.theti).*I)+muT[i].*(1-Par.delta).*R);
            L2[i] = ((1-lamT[i]).*(S.*(1-qT[i].*Par.thets)+Par.xi.*(1-qT[i].*Par.theti).*I)+(1-muT[i]).*(1-Par.delta).*R);
        end
        if Par.TEST == 1
            mctest[i] = Par.c.*qT[i].*(Par.thets.*S+Par.theti.*I).*63000 ./52*1 ./YY_noep(Par,1);
            ctest[i] = Par.c.*(qT[i].*(Par.thets.*S+Par.theti.*I))^(1+Par.nu)./(1+Par.nu).*63000 ./52*1. /YY_noep(Par,1);
            masstest[i] = qT[i].*(Par.thets.*ST[i]+Par.theti.*IT[i]);
        else
            mctest[i] = 0;
            ctest[i] = 0;
            masstest[i] = qT[i].*(Par.thets.*ST[i]+Par.theti.*IT[i]);
        end
        if I>0
            Rknot[i]=(1 ./(Par.gamma)).*Ipr./(1-S-R);
        else
            Rknot[i] = 0;
        end

    end

    # Calculate output variables
    #cumdeaths = maximum(cumsum(Par.gamma.*Par.delta.*max.(1 .- ST .- RT, 0))); # Cum deaths
    #Ycumsum1 = YT[1]-YY_noep(Par,1);
    #for i = 2:TT
    #    Ycumsum1 = Ycumsum1 + Par.beta.^(i-1).*(YT[i]-YY_noep(Par,1));
    #end
    #Ycumsum1 = (1-Par.beta).*Ycumsum1; # Cum output loss rel to no ep. ss
    #compvar = 1/YY_noep(Par,1)*exp((1-Par.beta).*(log(YT[1]-C1T[1])-C2T[1])+Par.beta*(VT[2]));
    #compvar = 1-compvar; # Welfare loss rel to no ep. ss
    #S52 = ST[52]; # Susceptible at 52 weeks
    #I52 = 1-ST[52]-RT[52]; # Infected at 52 weeks
    #Imax = maximum(1 .-ST .-RT); #  Max number infected
    #R01 = Rknot[1]; # R0 at 1 week
    #R26 = Rknot[26]; # R0 at 26 weeks
    #R52 = Rknot[52]; # R0 at  52 weeks
    #println(ST[TT])
    #println(1-ST[TT]-RT[TT])
    #println(Rknot[TT])
    #println(Imax)
    #percentmaxeval = round(nmaxeval / (nsolve) * 100, digits=2);
    #println("MAXEVAL_REACHED $nmaxeval times ($percentmaxeval%).")
    #percentforcedstop = round(nforcedstop / (nsolve) * 100, digits=2);
    #println("FORCED_STOP $nforcedstop times ($percentforcedstop%).")

    nsolve = sum(rets .!= 0);
    for value in setdiff(unique(rets), [0])
        nvalue = sum(rets .== value);
        println("$(string(ret_message[string(value)])) $nvalue times ($(round(nvalue/nsolve*100, digits=2))%).")
    end

    return Dict("Par" => Par, "TT" => TT, "ST" => ST, "IT" => IT, "RT" => RT,
    "YT" => YT, "C1T" => C1T, "Rknot" => Rknot, "masstest" => masstest, "mctest" => mctest)

end


text:
# draw robustness evaluation pictures

using DataFrames, Gadfly

folderloc = "/home/dehann/Videos/"

NN = Int[25,50,75,100,125,150,175,200,250,300,350,400,500,600,700,800,900,1000]

NN = Int[25,50,150,350,700]

# N = NN[12]

alllikelihoods = DataFrame[]
allerrors = DataFrame[]

computetime = Float64[]

for N in NN

nLK = "slamedonutAllLK_$(N).txt"
fi = joinpath(folderloc, nLK)
df = readtable(fi, header=false)

push!(alllikelihoods, DataFrame(
variable = union(["L$(i)" for i in 3:4],["X$(i)" for i in 1:13]),
likelihood = df[:x1].data[3:end],
particles = N*ones(length(df[:x1].data[3:end])),
Particles="$(N)"
))

nPE = "slamedonutAllPE_$(N).txt"
fi = joinpath(folderloc, nPE)
df = readtable(fi, header=false)

push!(allerrors, DataFrame(
variable = union(["L$(i)" for i in 3:4],["X$(i)" for i in 1:13]),
distance = df[:x1].data[3:end],
particles = N*ones(length(df[:x1].data[3:end])),
Particles="$(N)"
))

nPE = "slamedonutBatchTime_$(N).txt"
fi = joinpath(folderloc, nPE)
fid = open(fi,"r")
vals = readline(fid)
close(fid)
val = parse(Float64,vals)
push!(computetime, val)

end

vDF = vcat(alllikelihoods[end:-1:1]...)

pl = plot(vDF,Geom.point,color=:Particles,
x=:variable,y=:likelihood,
Guide.xticks(orientation=:vertical))

draw(PDF("slamedonutLikelihoods.pdf",12cm,7cm),pl)

vDF = vcat(alllikelihoods...)

pl = plot(vDF, Geom.boxplot,
x=:particles, y=:likelihood,
Coord.cartesian(xmax=1025))

draw(PDF("slamedonutallLK.pdf",12cm,7cm),pl)

vDF = vcat(allerrors[end:-1:1]...)

pl = plot(vDF,Geom.point,color=:Particles,
x=:variable,y=:distance,
Guide.xticks(orientation=:vertical))

draw(PDF("slamedonutPoseErrs.pdf",12cm,7cm),pl)

vDF = vcat(allerrors...)

pl = plot(vDF, Geom.boxplot,
x=:particles, y=:distance,
Coord.cartesian(xmax=1025))

draw(PDF("slamedonutallPE.pdf",12cm,7cm),pl)

DFcomp = DataFrame(
particles=NN,
seconds=computetime,
time="measured"
)
DFcompE = DataFrame(
particles=NN,
seconds=computetime/5,
time="possible"
)

DFc = vcat(DFcomp, DFcompE)

pl = plot(DFc,x=:particles,y=:seconds, color=:time, Geom.line)

draw(PDF("slamedonutcomputetimes.pdf",12cm,7cm),pl)

#

ext: jl


text:
using BenchmarkTools
import Flux: Flux, gradient, gpu
using CUDA
using Zygote: pullback

fw(m, x) = m(x)
bw(back) = back(1f0)
fwbw(m, ps, x) = gradient(() -> sum(m(x)), ps)

fw_recur(m, X::AbstractVector{<:AbstractArray}) = [m(x_t) for x_t in X]
fw_recur(m, X::AbstractArray{<:Number}) = m(X)
bw_recur(back) = back(1f0)
fwbw_recur(m, ps, X) = gradient(() -> sum(sum(fw_recur(m, X))), ps)

function run_benchmark_recur(model, x; cuda=false, show=false)

    t_fw, t_bw, t_fwbw = if cuda
        model_gpu = model |> gpu
        x_gpu = x |> gpu

        ps = Flux.params(model_gpu)
        y, back = pullback(() -> sum(sum(fw_recur(model_gpu, x_gpu))), ps)

        CUDA.allowscalar(false)
        # CUDA.device!(3)
        fw_recur(model_gpu, x_gpu); GC.gc(); CUDA.reclaim(); #warmup
        t_fw = @benchmark CUDA.@sync(fw_recur($model_gpu, $x_gpu)) teardown=(GC.gc(); CUDA.reclaim())

        bw_recur(back); GC.gc(); CUDA.reclaim(); #warmup
        t_bw = @benchmark CUDA.@sync(bw_recur($back)) teardown=(GC.gc(); CUDA.reclaim())

        fwbw_recur(model_gpu, ps, x_gpu); GC.gc(); CUDA.reclaim(); #warmup
        t_fwbw = @benchmark CUDA.@sync(fwbw_recur($model_gpu, $ps, $x_gpu)) teardown=(GC.gc(); CUDA.reclaim())

        t_fw, t_bw, t_fwbw
    else
        ps = Flux.params(model)
        y, back = pullback(() -> sum(sum(fw_recur(model, x))), ps)

        Flux.reset!(model)
        fw_recur(model, x)  #warmup
        t_fw = @benchmark fw_recur($model, $x)

        Flux.reset!(model)
        bw_recur(back)  #warmup
        t_bw = @benchmark bw_recur($back)

        Flux.reset!(model)
        fwbw_recur(model, ps, x) # warmup
        t_fwbw = @benchmark fwbw_recur($model, $ps, $x)

        Flux.reset!(model)
        t_fw, t_bw, t_fwbw
    end

    if show
        println("\t\tforward: ")
        @show t_fw
        println("\t\tbackward: ")
        @show t_bw
        println("\t\tforw and back: ")
        @show t_fwbw
    end

    (fw=t_fw, bw=t_bw, fwbw=t_fwbw)
end


text:
using JuLIP, ForwardDiff, StaticArrays, BenchmarkTools, Test, LinearAlgebra

# simple EAM-like potential
f(R) = sqrt( 1.0 + sum( exp(-norm(r)) for r in R ) )

# hand-coded gradient
function f_d(R::AbstractVector{<:JVec})
   ∇f = zeros(JVecF, length(R))
   ρ̄ = sum( exp(-norm(r)) for r in R )
   dF = 0.5 * (1.0 + ρ̄)^(-0.5)
   return [ dF * (- exp(-norm(r))) * r/norm(r)  for r in R ]
end

# ForwardDiff gradient
f_fd(R) =  collect(ForwardDiff.gradient( T -> f(vecs(T)),  collect(mat(R)[:]) ) |> vecs)

# turn it into an AD potential
V = JuLIP.Potentials.ADPotential(f, Inf)

# a test configuration
nneigs = 30
R0 = [ @SVector rand(3) for n = 1:nneigs ]

# check that all of them are the same
println(@test V(R0) == f(R0))
println(@test (@D V(R0)) == f_fd(R0))
println(@test f_d(R0) ≈ f_fd(R0))

# timings
print("Timing for      f: ")
@btime f($R0);
print("Timing for      V: ")
@btime V($R0);
print("Timing for    f_d: ")
@btime f_d($R0);
print("Timing for   f_fd: ")
@btime f_fd($R0);
print("Timing for   @D V: ")
@btime (@D $V($R0));


text:
module JaynesCummings

global const h_bar = 1.0#1.05457173E-34 # in Joules*seconds

export a, a_dagger

include("gen_initialstate.jl")
include("rabi_hamiltonian.jl")
include("gen_timeevoarray.jl")
include("calc_fullevo.jl")
include("partialtrace.jl")
include("print_matrix.jl")
include("get_eig_hamiltonian.jl")
include("get_prob.jl")

function a(cutoffN) # Truncated matrix for the lowering operator
  out = zeros(Complex128,cutoffN,cutoffN);
  for i = 1:cutoffN-1
      out[i,i+1] = sqrt(i);
  end
  return out
end

function a_dagger(cutoffN) # Truncated matrix for the raising operator
  out = zeros(Complex128,cutoffN,cutoffN)
  for i = 1:cutoffN-1
    out[i+1,i] = sqrt(i);
  end
  return out
end

end


text:
# Gives the fidelity and success probability of several well known
# distillation protocols for 2 copies of a bipartite state consisting of two qubits

export toBellBasis
export BBPSSWParam
export DEJMPSParam
export EPLDParam

import Base.convert
convert(::Type{Float64}, x::Array{Float64,1}) = x[1]
convert(::Type{Float64}, x::Array{Float64,2}) = x[1]

toMatrix(x::VecOrMat) = x*x'

e0 = [1; 0]
e1 = [0; 1]

proj0 = toMatrix(e0)
proj1 = toMatrix(e1)

phi_plus = 1/√2 * (kron(e0,e0) + kron(e1,e1))
phi_min  = 1/√2 * (kron(e0,e0) - kron(e1,e1))
psi_plus = 1/√2 * (kron(e0,e1) + kron(e1,e0))
psi_min  = 1/√2 * (kron(e0,e1) - kron(e1,e0))
basis = Any[phi_plus, phi_min, psi_plus, psi_min]

function toBellBasis(rho::AbstractMatrix)
  output = zeros(size(rho))
  i = 1
  for x in basis
    j = 1
    for y in basis
      output[i, j] = (x' * rho * y)[1]
      j += 1
    end
    i += 1
  end
  return output
end

""" `(F, p_succ) = DEJMPSParam(rho)`

Determines the performance of DEJMPS for two copies of the input state *rho* (consisting of two qubits) and returns the fidelity and success probability achieved.
The procedure optimises over local single qubit rotations that permute the bell diagonal coefficients before the CNOTs to achieve the highest output fidelity.
returns (F_out, p_succ)
"""
function DEJMPSParam(rho::AbstractMatrix)
  M = toBellBasis(rho);
  m1 = M[1, 1] # phi_+
  m2 = M[2, 2] # phi_-
  m3 = M[3, 3] # psi_+
  m4 = M[4, 4] # psi_-

  p1, p2, p3, p4 = sort([m1; m2; m3; m4], rev = true)

  p = (p1 + p4)^2 + (p2 + p3)^2
  F = (p1^2 + p4^2)/p
  return (F, p)
end

""" `(F, p_succ) = BBPSSWParam(F)`

Determines the performance of BBPSSW for two copies of the input state (consisting of two qubits) with fidelity *F* each.
returns (F_out, p_succ)
"""
function BBPSSWParam(F::Real)
  p_succ = F^2 + 2F*(1-F)/3 + 5((1-F)/3)^2
  F_out = (F^2 + ((1-F)/3)^2)/p_succ
  return (F_out, p_succ)
end

""" `(F, p_succ) = BBPSSWParam(rho)`

Determines the performance of BBPSSW for two copies of the input state *rho* (consisting of two qubits).
returns (F_out, p_succ)
"""
function BBPSSWParam(rho::AbstractMatrix)
  F = max(diag(toBellBasis(rho))...)
  return BBPSSWParam(F)
end

""" `(F, p_succ) = EPLDParam(p, pd)`

Determines the performance of EPL-D for rStateCorrPhase.
The inputs are the *p* and *pd* parameters of rStateCorrPhase.
returns (F_out, p_succ)
"""
function EPLDParam(p::Number, pd::Number)
  return (pd, 0.5*p^2)
end

""" `(F, p_succ) = EPLDParam(rho)`

Determines the performance of EPL-D for an arbitrary bipartite state of dimensions nA=nB=4.
returns (F_out, p_succ)
"""
function EPLDParam(rho::AbstractMatrix)

  #Define the CNOT operation:
  CNOT = [1 0 0 0; 0 1 0 0; 0 0 0 1; 0 0 1 0]

  #Alice and Bob apply CNOT:
  rho2 = kron(CNOT, CNOT) * rho * kron(CNOT, CNOT)'

  #Alice and Bob measure the flags and postselect them to be 11.
  #Projector on the success space on one side
  projSucc = kron(eye(2), proj1)

  #Success probability:
  p_succ = trace(kron(projSucc,projSucc) * rho2)

  #State on 2 qubit space after postselection:
  rhoOut = partialtrace((kron(projSucc,projSucc) * rho2 * kron(projSucc,projSucc)')/p_succ, [2,4], [2,2,2,2])

  #Output fidelity
  M = toBellBasis(rhoOut)

  F = max(M[1, 1], M[2, 2], M[3, 3], M[4, 4])

  return (F, p_succ)
end


text:
import RootSolvers
const RS = RootSolvers
import Plots

FT = Float64;
import Thermodynamics as TD
const ICP = TD.InternalClimaParams
import CLIMAParameters as CP
struct EarthParameterSet <: CP.AbstractEarthParameterSet end;
const param_set = EarthParameterSet();

TD.print_warning() = false
profiles = TD.TestedProfiles.PhaseEquilProfiles(param_set, Array{FT});
(; ρ, q_tot) = profiles
T_true = profiles.T
prof_pts = (ρ, T_true, q_tot)

dims = (6, 6, 6);
ρ = range(min(ρ...), stop = max(ρ...), length = dims[1]);
T_true = range(min(T_true...), stop = max(T_true...), length = dims[2]);
q_tot = range(min(q_tot...), stop = max(q_tot...), length = dims[3]);

ρ_all = Array{FT}(undef, prod(dims));
T_true_all = Array{FT}(undef, prod(dims));
q_tot_all = Array{FT}(undef, prod(dims));

linear_indices = LinearIndices((1:dims[1], 1:dims[2], 1:dims[3]));

numerical_methods = (
    RS.SecantMethod,
    RS.NewtonsMethod,
    # RS.NewtonsMethodAD, # we need to relax diagonalization somewhere
    # RS.RegulaFalsiMethod, # hit assertion error, bounds need adjusted for 3D space
)

ts = Dict(
    NM => Array{Union{TD.ThermodynamicState, Nothing}}(undef, prod(dims))
    for NM in numerical_methods
)
ts_no_err = Dict(
    NM => Array{TD.ThermodynamicState}(undef, prod(dims)) for
    NM in numerical_methods
)

@inbounds for i in linear_indices.indices[1]
    @inbounds for j in linear_indices.indices[2]
        @inbounds for k in linear_indices.indices[3]
            n = linear_indices[i, j, k]
            ρ_all[n] = ρ[i]
            T_true_all[n] = T_true[j]
            q_tot_all[n] = q_tot[k]
            e_int = TD.internal_energy(
                param_set,
                T_true[j],
                TD.PhasePartition(q_tot[k]),
            )

            @inbounds for NM in numerical_methods
                TD.error_on_non_convergence() = false
                ts_no_err[NM][n] = TD.PhaseEquil_dev_only(
                    param_set,
                    ρ[i],
                    e_int,
                    q_tot[k];
                    sat_adjust_method = NM,
                    maxiter = 10,
                )
                TD.error_on_non_convergence() = true
                # @show n / prod(dims) * 100
                try
                    ts[NM][n] = TD.PhaseEquil_dev_only(
                        param_set,
                        ρ[i],
                        e_int,
                        q_tot[k];
                        sat_adjust_method = NM,
                        maxiter = 10,
                    )
                catch
                    ts[NM][n] = nothing
                end
            end
        end
    end
end

# folder = "sat_adjust_analysis"
folder = @__DIR__
mkpath(folder)

function save_binned_surface_plots(
    x,
    y,
    z,
    title,
    filename,
    n_plots = (3, 3),
    z_label_prefix = "z",
    n_digits = 5;
    xlims = (:auto, :auto),
    ylims = (:auto, :auto),
    label = label,
    ref_points = nothing,
)
    n_z_partitions = prod(n_plots)
    local z_min_global, z_max_global
    try
        z_min_global = min(z...)
        z_max_global = max(z...)
    catch
        z_min_global = 0
        z_max_global = 1
    end
    Δz = (z_max_global - z_min_global) / n_z_partitions
    z_min = ntuple(i -> z_min_global + (i - 1) * Δz, n_z_partitions)
    z_max = ntuple(i -> z_min_global + (i) * Δz, n_z_partitions)
    p = []
    for i in 1:n_z_partitions
        mask = z_min[i] .<= z .<= z_max[i]
        x_i = x[mask]
        y_i = y[mask]
        sz_min = string(z_min[i])[1:min(n_digits, length(string(z_min[i])))]
        sz_max = string(z_max[i])[1:min(n_digits, length(string(z_max[i])))]
        p_i = Plots.plot(
            x_i,
            y_i,
            title = "$(title), in ($sz_min, $sz_max)",
            seriestype = :scatter,
            markersize = 5,
            label = label,
            xlims = xlims,
            ylims = ylims,
        )
        if ref_points ≠ nothing
            ref_mask = z_min[i] .<= ref_points[3] .<= z_max[i]
            Plots.plot!(
                ref_points[1][ref_mask],
                ref_points[2][ref_mask],
                seriestype = :scatter,
                markersize = 5,
                label = "ref points",
            )
        end
        push!(p, p_i)
    end
    Plots.plot(p..., layout = n_plots, legend = false)
    Plots.savefig(filename)
end;
# Full 3D scatter plot
function plot3D(ts_no_err, ts, NM; converged)
    mask = converged ? ts .≠ nothing : ts .== nothing

    c_name = converged ? "converged" : "non_converged"
    label = converged ? "converged" : "non-converged"
    casename = converged ? "converged" : "non-converged"
    nm_name = nameof(NM)
    filename = "3DSpace_$(c_name)_$nm_name.svg"

    ρ_mask = ρ_all[mask]
    q_tot_mask = q_tot_all[mask]
    T_mask = T_true_all[mask]
    pts = (ρ_mask, T_mask, q_tot_mask)
    Plots.plot(
        pts...,
        color = "blue",
        seriestype = :scatter,
        markersize = 7,
        label = casename,
    )
    Plots.plot!(
        prof_pts...,
        color = "red",
        seriestype = :scatter,
        markersize = 7,
        label = "tested thermo profiles",
    )
    Plots.plot!(
        xlabel = "Density",
        ylabel = "Temperature",
        zlabel = "Total specific humidity",
        title = "3D input to PhaseEquil",
        xlims = (min(ρ_all...), max(ρ_all...)),
        ylims = (min(T_true_all...), max(T_true_all...)),
        zlims = (min(q_tot_all...), max(q_tot_all...)),
        camera = (50, 50),
        # camera = (50,70),
    )
    Plots.savefig(joinpath(folder, filename))
end

# 2D binned scatter plots
function plot2D_slices(ts_no_err, ts, NM; converged)
    mask = converged ? ts .≠ nothing : ts .== nothing
    ρ_mask = ρ_all[mask]
    T_mask = T_true_all[mask]
    q_tot_mask = q_tot_all[mask]
    c_name = converged ? "converged" : "non_converged"
    label = converged ? "converged" : "non-converged"
    short_name = converged ? "C" : "NC"
    nm_name = nameof(NM)
    filename = "2DSlice_$(c_name)_$nm_name.svg"
    filename = joinpath(folder, filename)
    save_binned_surface_plots(
        ρ_mask,
        T_mask,
        q_tot_mask,
        short_name,
        filename;
        xlims = (min(ρ_all...), max(ρ_all...)),
        ylims = (min(T_true_all...), max(T_true_all...)),
        label = label,
        ref_points = prof_pts,
    )
end

# https://github.com/jheinen/GR.jl/issues/278#issuecomment-587090846
ENV["GKSwstype"] = "nul"

for NM in numerical_methods
    plot3D(ts_no_err[NM], ts[NM], NM; converged = false)
    plot3D(ts_no_err[NM], ts[NM], NM; converged = true)
    plot2D_slices(ts_no_err[NM], ts[NM], NM; converged = true)
    plot2D_slices(ts_no_err[NM], ts[NM], NM; converged = false)
end

convergence_percent = Dict()
for NM in numerical_methods
    convergence_percent[NM] = count(ts[NM] .≠ nothing) / length(ts[NM])
end
println("Convergence percentages:")
for (k, v) in convergence_percent
    println("$k = $v")
end


text:
using Distributions

export GaussianKernel, BoxKernel, RotationalKernel, evaluate

"""
    abstract type Kernel

Kernels serve as the backbone of broadening and convolution operations. All kernels have a common interface:

    evaluate(::Kernel; maxwidth=8, length=100)

which provides the evaluated PDF of the kernel. Note that the user should not have to provide the evaluated points, so each kernel should have a way of automatically choosing its range. For instance, a Gaussian kernel may choose to automatically evaluate 1000 points from -8σ to 8σ.
"""
abstract type Kernel end

evaluate(k::Kernel; maxwidth = 8, length = 101) = k.(_bounds(k, maxwidth, length))

## Implementation Note
# For each new kernel, implement a _pdf function and _bounds function as below
# which serves as the functional form of your kernel this way, we can limit the
# amount of redefinitions of evaluate.

"""
    GaussianKernel(σ::Number) <: Kernel

This is a simple Gaussian kernel which is evaluated as a normal distribution with variance σ^2 and a pdf of

    ``p(x) = \\frac{1}{\\sigma\\sqrt{2\\pi}} exp\\left(-\\frac{x^2}{2\\sigma^2}\\right)``

The default evaluation width is ± 4σ.
"""
struct GaussianKernel <: Kernel
    σ::Number
end
(k::GaussianKernel)(x) = Normal(0, k.σ).pdf(x)
_bounds(k::GaussianKernel, width, n) = range(-width / 2 * k.σ, width / 2 * k.σ, length = n)

"""
    BoxKernel(width::Number) <: Kernel

This is a simple Box (top-hat) kernel with a pdf of

    ``p(x) = \\begin{cases}
        1/width & \abs{x} < width/2 \\
        0 & \\
    \\end{cases}``

The default evaluation width is ± 4*width
"""
struct BoxKernel <: Kernel
    width::Number
end
(k::BoxKernel)(x) = Uniform(-k.width / 2, k.width / 2).pdf(x)
_bounds(k::BoxKernel, width, n) = range(-width / 2 * k.width, width / 2 * k.width, length = n)

struct RotationalKernel{T <: Number} <: Kernel
    dv::T
    vsini::T
end

RotationalKernel(dv, vsini) = RotationalKernel(promote(dv, vsini)...)


text:
export extractatoms
export LowRankChol, ShiftChol, SVDChol
export MomentMatrixWeightSolver, MomentVectorWeightSolver

using RowEchelon
using SemialgebraicSets

# Solve the system
# y = [U 0] * y
# where y = x[end:-1:1]
# which is
# y = U * β
# The code for solving the system with reordered schur has been moved to SemialgebraicSets.jl
#function solve_system(U, x)
#    m, r = size(U)
#    @assert m == length(x)
#    n = nvariables(x)
#    v = variables(x)
#    pivots = [findfirst(j -> U[j, i] != 0, 1:m) for i in 1:r]
#    if any(pivots .== 0)
#        keep = pivots .> 0
#        pivots = pivots[keep]
#        r = length(pivots)
#        U = U[:, keep]
#    end
#    β = x[m+1-pivots]
#    function multisearch_check(y)
#        idxs = multisearch(x, y)
#        if any(idxs .== 0)
#            error("Missing monomials $(y[idxs .== 0]) in $(x)")
#        end
#        idxs
#    end
#    Ns = [U[m+1-reverse(multisearch_check(v[i] * β)), :] for i in 1:n]
#    λ = rand(n)
#    λ /= sum(λ)
#    N = sum(λ .* Ns)
#    Z = schurfact(N)[:Z]
#    vals = [Vector{Float64}(n) for j in 1:r]
#    for j in 1:r
#        qj = Z[:, j]
#        for i in 1:n
#            vals[j][i] = dot(qj, Ns[i] * qj)
#        end
#    end
#    r, vals
#end

function build_system(U::AbstractMatrix, basis::MB.MonomialBasis, ztol, args...)
    # System is
    # y = [U 0] * y
    # where y = x[end:-1:1]
    # which is
    # y = U * β
    m = length(basis)
    r = size(U, 1)
    pivots = [findfirst(j -> U[i, j] != 0, 1:m) for i in 1:r]
    if any(pivots .== 0)
        keep = pivots .> 0
        pivots = pivots[keep]
        r = length(pivots)
        U = U[keep, :]
    end
    monos = basis.monomials
    β = monovec(monos[m + 1 .- pivots]) # monovec makes sure it stays sorted, TypedPolynomials wouldn't let it sorted
    function equation(i)
        if iszero(r) # sum throws ArgumentError: reducing over an empty collection is not allowed, if r is zero
            z = zero(eltype(β)) * zero(eltype(U))
            s = z + z # For type stability
        else
            s = sum(j -> β[r+1-j] * U[j, i], 1:r)
        end
        s - monos[m+1-i]
    end
    system = filter(p -> maxdegree(p) > 0, map(equation, 1:length(monos)))
    # Type instability here :(
    if mindegree(monos) == maxdegree(monos) # Homogeneous
        projectivealgebraicset(system, Buchberger(ztol), args...)
    else
        algebraicset(system, Buchberger(ztol), args...)
    end
end

"""
    LowRankChol

Method for computing a ``r \\times n`` matrix `U` of a ``n \\times n`` rank ``r`` psd matrix `Q` such that `Q = U'U`.
"""
abstract type LowRankChol end

"""
    ShiftChol <: LowRankChol

Shift the matrix by `shift` times the identity matrix before cholesky.
"""
struct ShiftChol{T} <: LowRankChol
    shift::T
end

"""
    SVDChol <: LowRankChol

Use SVD decomposition.
"""
struct SVDChol <: LowRankChol end

"""
    MultivariateMoments.lowrankchol(Q::AbstractMatrix, dec::LowRankChol, ranktol)

Returns a ``r \\times n`` matrix ``U`` of a ``n \\times n`` rank ``r`` positive semidefinite matrix ``Q`` such that ``Q = U^\\top U``.
The rank of ``Q`` is the number of singular values larger than `ranktol```{} \\cdot \\sigma_1`` where ``\\sigma_1`` is the largest singular value.
"""
function lowrankchol end

function lowrankchol(M::AbstractMatrix, dec::ShiftChol, ranktol)
    m = LinearAlgebra.checksquare(M)
    U = cholesky(M + dec.shift * I).U
    σs = map(i -> (U[i, i])^2, 1:m)
    nM = maximum(σs)
    tol = nM * ranktol
    rm = findall(σs .<= tol)
    if isempty(rm)
        cM = ranktol
    else
        cM = maximum(σs[rm]) / nM
    end
    nM, cM, U[σs .> tol, :]
end
function lowrankchol(M::AbstractMatrix, dec::SVDChol, ranktol)
    F = svd(M)
    nM = F.S[1] # norm of M
    tol = nM * ranktol
    r = something(findfirst(σ2 -> σ2 <= tol, F.S), 0)
    if iszero(r)
        cM = ranktol
        r = length(F.S)
    else
        cM = F.S[r] / nM
        r -= 1
    end
    nM, cM, (F.U[:, 1:r] * Diagonal(sqrt.(F.S[1:r])))'
end

"""
    MultivariateMoments.computesupport!(ν::MomentMatrix, ranktol, [dec])

Computes the `support` field of `ν`.
The `ranktol` and `dec` parameters are passed as is to the [`lowrankchol`](@ref) function.
"""
function computesupport!(μ::MomentMatrix, ranktol::Real, dec::LowRankChol, args...)
    # We reverse the ordering so that the first columns corresponds to low order monomials
    # so that we have more chance that low order monomials are in β and then more chance
    # v[i] * β to be in μ.x
    M = getmat(μ)
    m = LinearAlgebra.checksquare(M)
    M = M[m:-1:1, m:-1:1]
    nM, cM, U = lowrankchol(M, dec, ranktol)
    W = Matrix(U)
    # If M is multiplied by λ, W is multiplied by √λ
    # so we take √||M|| = √nM
    rref!(W, √(nM) * cM / sqrt(m))
    #r, vals = solve_system(U', μ.x)
    μ.support = build_system(W, μ.basis, √cM, args...) # TODO determine what is better between ranktol and sqrt(ranktol) here
end

function computesupport!(μ::MomentMatrix, ranktol::Real, args...)
    return computesupport!(μ::MomentMatrix, ranktol::Real, SVDChol(), args...)
end

# Determines weight

"""
    struct MomentMatrixWeightSolver
        rtol::T
        atol::T
    end

Given a moment matrix `ν` and the atom centers,
determine the weights by solving a linear system over all the moments
of the moment matrix, keeping duplicates (e.g., entries corresponding to the same monomial).

If the moment values corresponding to the same monomials are known to be equal
prefer [`MomentVectorWeightSolver`](@ref) instead.
"""
struct MomentMatrixWeightSolver
end

function solve_weight(ν::MomentMatrix{T}, centers, solver::MomentMatrixWeightSolver) where {T}
    vars = variables(ν)
    A = Matrix{T}(undef, length(ν.Q.Q), length(centers))
    vbasis = vectorized_basis(ν)
    for i in eachindex(centers)
        η = dirac(vbasis.monomials, vars => centers[i])
        A[:, i] = moment_matrix(η, ν.basis.monomials).Q.Q
    end
    return A \ ν.Q.Q
end

"""
    struct MomentVectorWeightSolver{T}
        rtol::T
        atol::T
    end

Given a moment matrix `ν` and the atom centers, first convert the moment matrix
to a vector of moments, using [`measure(ν; rtol=rtol, atol=atol)`](@ref measure)
and then determine the weights by solving a linear system over the monomials obtained.

If the moment values corresponding to the same monomials can have small differences,
[`measure`](@ref) can throw an error if `rtol` and `atol` are not small enough.
Alternatively to tuning these tolerances [`MomentVectorWeightSolver`](@ref) can be used instead.
"""
struct MomentVectorWeightSolver{T}
    rtol::T
    atol::T
end
function MomentVectorWeightSolver{T}(; rtol=Base.rtoldefault(T), atol=zero(T)) where {T}
    return MomentVectorWeightSolver{T}(rtol, atol)
end
function MomentVectorWeightSolver(; rtol=nothing, atol=nothing)
    if rtol === nothing && atol === nothing
        return MomentVectorWeightSolver{Float64}()
    elseif rtol !== nothing
        if atol === nothing
            return MomentVectorWeightSolver{typeof(rtol)}(; rtol=rtol)
        else
            return MomentVectorWeightSolver{typeof(rtol)}(; rtol=rtol, atol=atol)
        end
    else
        return MomentVectorWeightSolver{typeof(atol)}(; atol=atol)
    end
end

function solve_weight(ν::MomentMatrix{T}, centers, solver::MomentVectorWeightSolver) where {T}
    μ = measure(ν; rtol=solver.rtol, atol=solver.atol)
    vars = variables(μ)
    A = Matrix{T}(undef, length(μ.x), length(centers))
    for i in eachindex(centers)
        A[:, i] = dirac(μ.x, vars => centers[i]).a
    end
    return A \ μ.a
end

"""
    extractatoms(ν::MomentMatrix, ranktol, [dec::LowRankChol], [solver::SemialgebraicSets.AbstractAlgebraicSolver])

Return an `AtomicMeasure` with the atoms of `ν` if it is atomic or `nothing` if
`ν` is not atomic. The `ranktol` and `dec` parameters are passed as is to the
[`lowrankchol`](@ref) function. By default, `dec` is an instance of
[`SVDChol`](@ref). The extraction relies on the solution of a system of
algebraic equations. using `solver`. For instance, given a
[`MomentMatrix`](@ref), `μ`, the following extract atoms using a `ranktol` of
`1e-4` for the low-rank decomposition and homotopy continuation to solve the
obtained system of algebraic equations.
```julia
using HomotopyContinuation
solver = SemialgebraicSetsHCSolver(; compile = false)
atoms = extractatoms(ν, 1e-4, solver)
```
If no solver is given, the default solver of SemialgebraicSets is used which
currently computes the Gröbner basis, then the multiplication matrices and
then the Schur decomposition of a random combination of these matrices.
For floating point arithmetics, homotopy continuation is recommended as it is
more numerically stable than Gröbner basis computation.
"""
function extractatoms(ν::MomentMatrix{T}, ranktol, args...; weight_solver = MomentMatrixWeightSolver()) where T
    computesupport!(ν, ranktol, args...)
    supp = ν.support
    if !iszerodimensional(supp)
        return nothing
    end
    centers = collect(supp)
    r = length(centers)
    weights = solve_weight(ν, centers, weight_solver)
    isf = isfinite.(weights)
    weights = weights[isf]
    centers = centers[isf]
    if isempty(centers)
        nothing
    else
        AtomicMeasure(variables(ν), WeightedDiracMeasure.(centers, weights))
    end
end


text:
# This doesn't run under "runtests.jl" but does under CI
using ImageMorphology
using Test

@testset "multithreaded" begin
    @test Threads.nthreads() > 1
    @testset "feature_transform" begin
        img = rand(100, 128) .> 0.9
        @test feature_transform(img; nthreads=Threads.nthreads()) ==
            feature_transform(img; nthreads=1)
        # Since the threaded implementation handles two dimensions specially, we should check 0d and 1d
        img = reshape([true])
        @test feature_transform(img; nthreads=Threads.nthreads()) ==
            feature_transform(img; nthreads=1) ==
            reshape([CartesianIndex()])
        img = rand(100) .> 0.9
        @test feature_transform(img; nthreads=Threads.nthreads()) ==
            feature_transform(img; nthreads=1)
    end
end


text:
####
#### Particle marginal Metropolis-Hastings sampler.
####

"""
    PMMH(n_iters::Int, smc_alg:::SMC, parameters_algs::Tuple{MH})

Particle independant Metropolis–Hastings and
Particle marginal Metropolis–Hastings samplers.

Note that this method is particle-based, and arrays of variables
must be stored in a [`TArray`](@ref) object.

Usage:

```julia
alg = PMMH(100, SMC(20, :v1), MH(1,:v2))
alg = PMMH(100, SMC(20, :v1), MH(1,(:v2, (x) -> Normal(x, 1))))
```

Arguments:

- `n_iters::Int` : Number of iterations to run.
- `smc_alg:::SMC` : An [`SMC`](@ref) algorithm to use.
- `parameters_algs::Tuple{MH}` : An [`MH`](@ref) algorithm, which includes a
sample space specification.
"""
mutable struct PMMH{space, A<:Tuple} <: InferenceAlgorithm
    n_iters::Int               # number of iterations
    algs::A                 # Proposals for state & parameters
end
function PMMH(n_iters::Int, algs::Tuple, space::Tuple)
    return PMMH{space, typeof(algs)}(n_iters, algs)
end
function PMMH(n_iters::Int, smc_alg::SMC, parameter_algs...)
    return PMMH(n_iters, tuple(parameter_algs..., smc_alg), ())
end

PIMH(n_iters::Int, smc_alg::SMC) = PMMH(n_iters, tuple(smc_alg), ())

function Sampler(alg::PMMH, model::Model, s::Selector)
    info = Dict{Symbol, Any}()
    spl = Sampler(alg, info, s)

    alg_str = "PMMH"
    n_samplers = length(alg.algs)
    samplers = Array{Sampler}(undef, n_samplers)

    space = Set{Symbol}()

    for i in 1:n_samplers
        sub_alg = alg.algs[i]
        if isa(sub_alg, Union{SMC, MH})
            samplers[i] = Sampler(sub_alg, model, Selector(Symbol(typeof(sub_alg))))
        else
            error("[$alg_str] unsupport base sampling algorithm $alg")
        end
        if typeof(sub_alg) == MH && sub_alg.n_iters != 1
            warn(
                "[$alg_str] number of iterations greater than 1" *
                "is useless for MH since it is only used for its proposal"
            )
        end
        space = union(space, sub_alg.space)
    end

    info[:old_likelihood_estimate] = -Inf # Force to accept first proposal
    info[:old_prior_prob] = 0.0
    info[:samplers] = samplers

    return spl
end

function step(model, spl::Sampler{<:PMMH}, vi::VarInfo, is_first::Bool)
    violating_support = false
    proposal_ratio = 0.0
    new_prior_prob = 0.0
    new_likelihood_estimate = 0.0
    old_θ = copy(vi[spl])

    Turing.DEBUG && @debug "Propose new parameters from proposals..."
    for local_spl in spl.info[:samplers][1:end-1]
        Turing.DEBUG && @debug "$(typeof(local_spl)) proposing $(local_spl.alg.space)..."
        propose(model, local_spl, vi)
        if local_spl.info[:violating_support] violating_support=true; break end
        new_prior_prob += local_spl.info[:prior_prob]
        proposal_ratio += local_spl.info[:proposal_ratio]
    end

    if violating_support
        # do not run SMC if going to refuse anyway
        accepted = false
    else
        Turing.DEBUG && @debug "Propose new state with SMC..."
        vi, _ = step(model, spl.info[:samplers][end], vi)
        new_likelihood_estimate = spl.info[:samplers][end].info[:logevidence][end]

        Turing.DEBUG && @debug "Decide whether to accept..."
        accepted = mh_accept(
          spl.info[:old_likelihood_estimate] + spl.info[:old_prior_prob],
          new_likelihood_estimate + new_prior_prob,
          proposal_ratio,
        )
    end

    if accepted
        spl.info[:old_likelihood_estimate] = new_likelihood_estimate
        spl.info[:old_prior_prob] = new_prior_prob
    else                      # rejected
        vi[spl] = old_θ
    end

    return vi, accepted
end

function sample(  model::Model,
                  alg::PMMH;
                  save_state=false,         # flag for state saving
                  resume_from=nothing,      # chain to continue
                  reuse_spl_n=0             # flag for spl re-using
                )

    spl = Sampler(alg, model)
    if resume_from !== nothing
        spl.selector = resume_from.info[:spl].selector
    end
    alg_str = "PMMH"

    # Number of samples to store
    sample_n = spl.alg.n_iters

    # Init samples
    time_total = zero(Float64)
    samples = Array{Sample}(undef, sample_n)
    weight = 1 / sample_n
    for i = 1:sample_n
        samples[i] = Sample(weight, Dict{Symbol, Any}())
    end

    # Init parameters
    vi = if resume_from === nothing
        vi_ = VarInfo(model)
    else
        resume_from.info[:vi]
    end
    n = spl.alg.n_iters

    # PMMH steps
    accept_his = Bool[]
    PROGRESS[] && (spl.info[:progress] = ProgressMeter.Progress(n, 1, "[$alg_str] Sampling...", 0))
    for i = 1:n
      Turing.DEBUG && @debug "$alg_str stepping..."
      time_elapsed = @elapsed vi, is_accept = step(model, spl, vi, i==1)

      if is_accept # accepted => store the new predcits
          samples[i].value = Sample(vi, spl).value
      else         # rejected => store the previous predcits
          samples[i] = samples[i - 1]
      end

      time_total += time_elapsed
      push!(accept_his, is_accept)
      if PROGRESS[]
        haskey(spl.info, :progress) && ProgressMeter.update!(spl.info[:progress], spl.info[:progress].counter + 1)
      end
    end

    println("[$alg_str] Finished with")
    println("  Running time    = $time_total;")
    accept_rate = sum(accept_his) / n  # calculate the accept rate
    println("  Accept rate         = $accept_rate;")

    if resume_from !== nothing   # concat samples
      pushfirst!(samples, resume_from.info[:samples]...)
    end
    c = Chain(-Inf, samples)       # wrap the result by Chain

    if save_state               # save state
      c = save(c, spl, model, vi, samples)
    end

    c
end

####
#### IMCMC Sampler.
####

"""
    IPMCMC(n_particles::Int, n_iters::Int, n_nodes::Int, n_csmc_nodes::Int)

Particle Gibbs sampler.

Note that this method is particle-based, and arrays of variables
must be stored in a [`TArray`](@ref) object.

Usage:

```julia
IPMCMC(100, 100, 4, 2)
```

Arguments:

- `n_particles::Int` : Number of particles to use.
- `n_iters::Int` : Number of iterations to employ.
- `n_nodes::Int` : The number of nodes running SMC and CSMC.
- `n_csmc_nodes::Int` : The number of CSMC nodes.
```

A paper on this can be found [here](https://arxiv.org/abs/1602.05128).
"""
mutable struct IPMCMC{T, F} <: InferenceAlgorithm
  n_particles::Int         # number of particles used
  n_iters::Int         # number of iterations
  n_nodes::Int         # number of nodes running SMC and CSMC
  n_csmc_nodes::Int         # number of nodes CSMC
  resampler::F           # function to resample
  space::Set{T}      # sampling space, emtpy means all
end
IPMCMC(n1::Int, n2::Int) = IPMCMC(n1, n2, 32, 16, resample_systematic, Set())
IPMCMC(n1::Int, n2::Int, n3::Int) = IPMCMC(n1, n2, n3, Int(ceil(n3/2)), resample_systematic, Set())
IPMCMC(n1::Int, n2::Int, n3::Int, n4::Int) = IPMCMC(n1, n2, n3, n4, resample_systematic, Set())
function IPMCMC(n1::Int, n2::Int, n3::Int, n4::Int, space...)
  _space = isa(space, Symbol) ? Set([space]) : Set(space)
  IPMCMC(n1, n2, n3, n4, resample_systematic, _space)
end

function Sampler(alg::IPMCMC, s::Selector)
  info = Dict{Symbol, Any}()
  spl = Sampler(alg, info, s)
  # Create SMC and CSMC nodes
  samplers = Array{Sampler}(undef, alg.n_nodes)
  # Use resampler_threshold=1.0 for SMC since adaptive resampling is invalid in this setting
  default_CSMC = CSMC(alg.n_particles, 1, alg.resampler, alg.space)
  default_SMC = SMC(alg.n_particles, alg.resampler, 1.0, false, alg.space)

  for i in 1:alg.n_csmc_nodes
    samplers[i] = Sampler(default_CSMC, Selector(Symbol(typeof(default_CSMC))))
  end
  for i in (alg.n_csmc_nodes+1):alg.n_nodes
    samplers[i] = Sampler(default_SMC, Selector(Symbol(typeof(default_CSMC))))
  end

  info[:samplers] = samplers

  return spl
end

function step(model, spl::Sampler{<:IPMCMC}, VarInfos::Array{VarInfo}, is_first::Bool)
    # Initialise array for marginal likelihood estimators
    log_zs = zeros(spl.alg.n_nodes)

    # Run SMC & CSMC nodes
    for j in 1:spl.alg.n_nodes
        reset_num_produce!(VarInfos[j])
        VarInfos[j] = step(model, spl.info[:samplers][j], VarInfos[j])[1]
        log_zs[j] = spl.info[:samplers][j].info[:logevidence][end]
    end

    # Resampling of CSMC nodes indices
    conditonal_nodes_indices = collect(1:spl.alg.n_csmc_nodes)
    unconditonal_nodes_indices = collect(spl.alg.n_csmc_nodes+1:spl.alg.n_nodes)
    for j in 1:spl.alg.n_csmc_nodes
        # Select a new conditional node by simulating cj
        log_ksi = vcat(log_zs[unconditonal_nodes_indices], log_zs[j])
        ksi = exp.(log_ksi .- maximum(log_ksi))
        c_j = wsample(ksi) # sample from Categorical with unormalized weights

        if c_j < length(log_ksi) # if CSMC node selects another index than itself
            conditonal_nodes_indices[j] = unconditonal_nodes_indices[c_j]
            unconditonal_nodes_indices[c_j] = j
        end
    end
    nodes_permutation = vcat(conditonal_nodes_indices, unconditonal_nodes_indices)

    VarInfos[nodes_permutation]
end

function sample(model::Model, alg::IPMCMC)

  spl = Sampler(alg)

  # Number of samples to store
  sample_n = alg.n_iters * alg.n_csmc_nodes

  # Init samples
  time_total = zero(Float64)
  samples = Array{Sample}(undef, sample_n)
  weight = 1 / sample_n
  for i = 1:sample_n
    samples[i] = Sample(weight, Dict{Symbol, Any}())
  end

  # Init parameters
  vi = empty!(VarInfo(model))
  VarInfos = Array{VarInfo}(undef, spl.alg.n_nodes)
  for j in 1:spl.alg.n_nodes
    VarInfos[j] = deepcopy(vi)
  end
  n = spl.alg.n_iters

  # IPMCMC steps
  if PROGRESS[] spl.info[:progress] = ProgressMeter.Progress(n, 1, "[IPMCMC] Sampling...", 0) end
  for i = 1:n
    Turing.DEBUG && @debug "IPMCMC stepping..."
    time_elapsed = @elapsed VarInfos = step(model, spl, VarInfos, i==1)

    # Save each CSMS retained path as a sample
    for j in 1:spl.alg.n_csmc_nodes
      samples[(i-1)*alg.n_csmc_nodes+j].value = Sample(VarInfos[j], spl).value
    end

    time_total += time_elapsed
    if PROGRESS[]
      haskey(spl.info, :progress) && ProgressMeter.update!(spl.info[:progress], spl.info[:progress].counter + 1)
    end
  end

  println("[IPMCMC] Finished with")
  println("  Running time    = $time_total;")

  Chain(0.0, samples) # wrap the result by Chain
end


text:
export boundary_point, tangent_cone

struct BoundaryPoint{S} <: SymbolicVariable
    set::S
    symbol::Symbol
end
struct TangentCone{S} <: SymbolicVariable
    set::S
    symbol::Symbol
end
boundary_point(S, s::Symbol) = BoundaryPoint(S, s)
tangent_cone(S, s::Symbol) = TangentCone(S, s)
function tangent_cone(S, b::BoundaryPoint)
    @assert S === b.set
    tangent_cone(S, b.symbol)
end
variablify(t::BoundaryPoint) = BoundaryPoint(variablify(t.set), t.symbol)
need_variablify(t::TangentCone) = need_variablify(t.set)
variablify(t::TangentCone) = TangentCone(variablify(t.set), t.symbol)
create_spaces(b::BoundaryPoint, spaces::Spaces) = create_spaces(b.set, spaces)
create_spaces(t::TangentCone, spaces::Spaces) = create_spaces(t.set, spaces)
function create_spaces(c::MembershipConstraint, spaces::Spaces)
    sub = create_spaces(c.member, spaces)
    sup = create_spaces(c.set, spaces)
    return merge_spaces(spaces, sub, sup)
end
clear_spaces(b::BoundaryPoint) = clear_spaces(b.set)
clear_spaces(t::TangentCone) = clear_spaces(t.set)
function clear_spaces(c::MembershipConstraint)
    clear_spaces(c.member)
    clear_spaces(c.set)
end
Sets.perspective_variable(b::BoundaryPoint) = Sets.perspective_variable(b.set)
Sets.perspective_variable(t::TangentCone) = Sets.perspective_variable(t.set)
function Sets.perspective_variable(c::MembershipConstraint)
    return synchronize_perspective(
        Sets.perspective_variable(c.member),
        Sets.perspective_variable(c.set)
    )
end
function set_space(
    space::Space,
    ::MembershipConstraint{
        <:Sets.LinearImage{<:BoundaryPoint},
        <:Sets.LinearImage{<:TangentCone}}
)
    return set_space(space, DualSpace)
end
using LinearAlgebra
function linear_algebraic_surface(set::Sets.PolarOf{<:Sets.Ellipsoid}, A, E)
    Q = set.set.Q
    return psd_constraint(Symmetric(-A * Q * E' - E * Q * A'))
end
function linear_algebraic_surface(set::Sets.PolarOf{<:Sets.ConvexPolySet}, A, E)
    v = MP.variables(set.set.p)
    DynamicPolynomials.@polyvar z[1:size(E, 1)]
    Ez = E' * z
    return JuMP.build_constraint(
        error,
        -dot(A' * z, [d(v => Ez) for d in MP.differentiate(set.set.p, v)]),
        SOSCone(),
    )
end
function JuMP.build_constraint(
    _error::Function,
    member::Sets.LinearImage{BoundaryPoint{S}},
    set::Sets.LinearImage{TangentCone{S}}
   ) where S <: Sets.PolarOf{<:Union{Sets.Ellipsoid, Sets.ConvexPolySet}}
    _set = member.set.set
    @assert set.set.set === _set
    return linear_algebraic_surface(_set, member.A, set.A)
end
function add_linear_algebraic_surface_domain(
    model, set::Sets.Ellipsoid, domain, A, E)
    Q = set.Q
    P = Symmetric(-A * Q * E' - E * Q * A')
    return _add_constraint_or_not(model, psd_in_domain(model, P, E' \ domain))
end
function add_linear_algebraic_surface(model, set::Sets.PolarOf{<:Sets.Piecewise}, A, E)
    for (i, si) in enumerate(set.set.sets)
        add_linear_algebraic_surface_domain(
            model, si, set.set.pieces[i], A, E
        )
    end
end
function JuMP.add_constraint(
    model::JuMP.Model,
    constraint::MembershipConstraint{
        <:Sets.LinearImage{<:BoundaryPoint},
        <:Sets.LinearImage{<:TangentCone}
    }
)
    _set = constraint.member.set.set
    @assert constraint.set.set.set === _set
    return add_linear_algebraic_surface(model, _set, constraint.member.A, constraint.set.A)
end


text:
using ..Serde: QDiffPayloadSerde, QDiffStaticPayloadSerde

"""
    QDiffConfig{T <: Unsigned}(
        base::BaseConfig,
        chunksize::Int,
        is_client::Bool
    ) <: Configuration

Apply the same quantization algorithm used by `QuantizedConfig`, but also
further reduces the size of the paylaod by adding client-side "diff-
deduplication", where unchanged chunks of data are not sent back to the server.
This configuration instanciates a [QDiffPayloadSerde](@ref).
"""
struct QDiffConfig{T <: Unsigned} <: Configuration
    base::BaseConfig
    payload_serde::QDiffPayloadSerde

    QDiffConfig{T}(
        base::BaseConfig,
        chunksize::Int,
        is_client::Bool
    ) where T <: Unsigned = new(
        base,
        QDiffPayloadSerde{T}(chunksize, is_client)
    )
end

"""
    QDiffStaticConfif(base, chunksize, is_client)

Same as QDiffConfig with a static quantization range of [-1.0, 1.0].
"""
struct QDiffStaticConfig{T <: Unsigned} <: Configuration
    base::BaseConfig
    payload_serde::QDiffStaticPayloadSerde

    QDiffStaticConfig{T}(
        base::BaseConfig,
        chunksize::Int,
        is_client::Bool
    ) where T <: Unsigned = new(
        base,
        QDiffStaticPayloadSerde{T}(chunksize, is_client)
    )
end


text:
# Autogenerated wrapper script for TinyXML_jll for i686-linux-musl-cxx11
export libtinyxml

JLLWrappers.@generate_wrapper_header("TinyXML")
JLLWrappers.@declare_library_product(libtinyxml, "libtinyxml.so")
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_library_product(
        libtinyxml,
        "lib/libtinyxml.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
using .Microclimate

"""
    apply_environment!(plant::Plant, env::AbstractMicroclimate, organs, u, t)

Apply a microclimate from Microclimates.jl to a `Plant`.
"""
apply_environment!(env::AbstractMicroclimate, plant::Plant, organs, u, t) = begin
    envind = calc_envind(plant, t)
    shoot, root = organs
    # Get the climate at a particular height and time
    # Use half the height for the median temperature/windspeed/humidity etc.
    # relating temperature to the distribution of biomass may be more accurate.
    shootenv = MicroclimInstant(env, envind, height(shoot)/2)
    rootenv = MicroclimInstant(env, envind, depth(root)/2)
    # Update any environment dependent variables
    apply_environment!(assimilation_pars(shoot), shoot, u, shootenv, rootenv)
    apply_environment!(assimilation_pars(root), root, u, rootenv, shootenv)
    nothing
end

# We need an environment data point for t = 0 even if we
# never use it - we may need to interpolate at t = 0.5
calc_envind(o, t) = ustrip(o.environment_start[] + t) + 1

"""
    apply_environment!(plant::Plant, env::MicroclimControl, organs, u, t)

Apply an environment to a Plant using live user input for variables,
instead of data.
"""
apply_environment!(env::MicroclimControl, plant::Plant, organs, u, t) = begin
    shoot, root = organs
    apply_environment!(assimilation_pars(shoot), shoot, u, env, env)
    apply_environment!(assimilation_pars(root), root, u, env, env)
    nothing
end

"""
    apply_environment!(a::AbstractKooijmanPhoto, o, u, shootenv, rootenv)

Apply environment to shoot with `AbstractKooijmanPhoto` assimilation model.
"""
apply_environment!(a::AbstractKooijmanPhoto, o, u, shootenv, rootenv) = begin
    v = a.vars
    update_temp!(o, airtemperature(shootenv))
    v.J_L_F = radiation(shootenv) * parconv
    # Use the mean soil water potential between zero and the full root depth
    # swp = mean_soilwaterpotential(rootenv.microclimate, depth(o), rootenv.t)
    swp = if typeof(rootenv) <: MicroclimControl
        soilwaterpotential(rootenv)
    else
        layermax(soilwaterpotential(rootenv.microclimate), rootenv)
    end
    v.soilwaterpotential = swp
    set_swp!(o, swp)
end
"""

Apply environment to shoot with any `AbstractAssim` assimilation model.
This simply adds the air temperature.
"""
apply_environment!(a::AbstractAssim, o, u, shootenv, rootenv) =
    update_temp!(o, airtemperature(shootenv))
"""
    apply_environment!(a::AbstractNAssim, o, u, rootenv, shootenv)

Apply environment to root with `AbstractNAssim` assimilation model.
This simply adds the soil temperature.
"""
apply_environment!(a::AbstractNAssim, o, u, rootenv, shootenv) =
    update_temp!(o, soiltemperature(rootenv))


text:
import RowEchelon
import LinearAlgebra

#=
Homework 8
Section 4.2
Problems 18
=#

A = [5 -2 3 ;
     -1 0 -1 ;
     0 -2 -2 ;
     -5 7 2]

sizeA = size(A)
println("Original Matrix, A:")
println(A)
print("NulA is a subspace of R")
println(sizeA[2])
print("ColA is a subspace of R")
println(sizeA[1])


text:
using ProgressMeter

include("./plot_ic_reward.jl")

folder = "M:/globus/InterestTB/Experiment1/"
ic = ItemCollection(folder)

# Experiment parameters
log_interval = 50
num_steps = 2000

p = plot()

configs = [
    Dict(
        :search => Dict("demon_learner" => "Q", "demon_update" => "ETB"),
        :label => "Q ETB"
    )
    Dict(
        :search => Dict("demon_learner" => "Q", "demon_update" => "EmphESARSA"),
        :label => "Q EmphESARSA"
    )
    Dict(
        :search => Dict("demon_learner" => "Q", "demon_update" => "InterestTB"),
        :label => "Q InterestTB"
    )
    Dict(
        :search => Dict("demon_learner" => "SR", "demon_update" => "ETB"),
        :label => "SR ETB"
    )
    Dict(
        :search => Dict("demon_learner" => "SR", "demon_update" => "EmphESARSA"),
        :label => "SR EmphESARSA"
    )
    Dict(
        :search => Dict("demon_learner" => "SR", "demon_update" => "InterestTB"),
        :label => "SR InterestTB"
    )
]

@showprogress for config in configs
    local search_ic = search(ic, config[:search])
    local best_ic = get_best(search_ic, ["demon_alpha_init", "demon_eta"], :ttmaze_direct_error)
    plot_ic_reward!(p, best_ic, :ttmaze_direct_error, log_interval; smooth_step=1, label=config[:label])
end

savefig("plotting/chunlok/generated_plots/Experiment1_rewards.svg")

max_stars_repo_name: MatthewMcLeod/curiosity

max_issues_repo_name: MatthewMcLeod/curiosity

max_forks_repo_name: MatthewMcLeod/curiosity


text:
using ..TimeScales: find_path
using EarthOrientation: getΔUT1, EOP_DATA, get
using LeapSeconds: offset_tai_utc, offset_utc_tai
using MuladdMacro

export
    NoOffsetError,
    NoPathError,
    getoffset,
    insideleap

struct NoPathError <: Base.Exception
    in_scale::String
    out_scale::String
end

function Base.showerror(io::IO, err::NoPathError)
    ins = err.in_scale
    out = err.out_scale
    print(io, "No conversion path between '$ins' and '$out' available.")
end

struct NoOffsetError <: Base.Exception
    in_scale::String
    out_scale::String
    in_type::String
    out_type::String
end

function Base.showerror(io::IO, err::NoOffsetError)
    ins = err.in_scale
    out = err.out_scale
    it = err.in_type
    ot = err.out_type
    print(io, "No conversion '$ins->$out' available. ",
          "If one of these is a custom time scale, you may need to define ",
          "`AstroTime.Epochs.getoffset(::$it, ::$ot, second, fraction, args...)`.")
end

function getoffset(s1::TimeScale, s2::TimeScale, _, _)
    err = NoOffsetError(string(s1), string(s2),
                        string(typeof(s1)), string(typeof(s2)))
    throw(err)
end

function j2000(second, fraction)
    (fraction + second) / SECONDS_PER_DAY * days
end

function getoffset(ep::Epoch{S}, scale::TimeScale) where S<:TimeScale
    path = find_path(S(), scale)
    isempty(path) && throw(NoPathError(string(S()), string(scale)))
    length(path) == 2 && return getoffset(S(), scale, ep.second, ep.fraction)
    total_offset = 0.0
    second = ep.second
    fraction = ep.fraction
    for i in 1:length(path) - 1
        offset::Float64 = getoffset(path[i], path[i+1], second, fraction)
        total_offset += offset
        second, fraction = apply_offset(second, fraction, offset)
    end
    return total_offset
end

"""
    getoffset(ep::Epoch, scale::TimeScale)

For a given epoch `ep` return the offset between its time scale and
another time `scale` in seconds.

# Example

```jldoctest
julia> tai = TAIEpoch(2000, 1, 1)
2000-01-01T00:00:00.000 TAI

julia> getoffset(tai, TT)
32.184
```
"""
function getoffset(ep::Epoch{S}, scale::TimeScale, args...) where S<:TimeScale
    return getoffset(S(), scale, ep.second, ep.fraction, args...)
end

@inline function apply_offset(second::Int64,
                              fraction::T,
                              from::S1,
                              to::S2)::Tuple{Int64, T} where {T, S1<:TimeScale, S2<:TimeScale}
    path = find_path(from, to)
    isempty(path) && throw(NoPathError(string(from), string(to)))
    length(path) == 2 && return _apply_offset(second, fraction, from, to)
    return _apply_offset((second, fraction), path...)
end

@inline function _apply_offset(second::Int64,
                               fraction::T,
                               from::S1,
                               to::S2)::Tuple{Int64, T} where {T, S1<:TimeScale, S2<:TimeScale}
    return apply_offset(second, fraction, getoffset(from, to, second, fraction))
end

@generated function _apply_offset(sf, path...)
    expr = :(sf)
    for i in 1:length(path) - 1
        s1 = path[i]
        s2 = path[i+1]
        expr = :(_apply_offset($expr..., $s1(), $s2()))
    end
    return quote
        Base.@_inline_meta
        $expr
    end
end

######
# TT #
######

const OFFSET_TAI_TT = 32.184

"""
    getoffset(TAI, TT, args...)

Return the fixed offset between [`TAI`](@ref) and [`TT`](@ref) in seconds.

# Example

```jldoctest
julia> getoffset(TAI, TT)
32.184
```
"""
getoffset(::InternationalAtomicTime, ::TerrestrialTime, args...) = OFFSET_TAI_TT

"""
    getoffset(TT, TAI, args...)

Return the fixed offset between [`TT`](@ref) and [`TAI`](@ref) in seconds.

# Example

```jldoctest
julia> getoffset(TT, TAI)
-32.184
```
"""
getoffset(::TerrestrialTime, ::InternationalAtomicTime, args...) = -OFFSET_TAI_TT

#######
# TCG #
#######

const JD77_SEC = -7.25803167816e8
const LG_RATE = 6.969290134e-10

"""
    getoffset(TCG, TT, second, fraction)

Return the linear offset between [`TCG`](@ref) and [`TT`](@ref) for the
current epoch (`second` after J2000 and `fraction`) in seconds.

# Example

```jldoctest
julia> getoffset(TCG, TT, 0, 0.0)
-0.5058332856685995
```
"""
function getoffset(::GeocentricCoordinateTime, ::TerrestrialTime, second, fraction)
    dt = second - JD77_SEC + fraction
    return -LG_RATE * dt
end

"""
    getoffset(TT, TCG, second, fraction)

Return the linear offset between [`TT`](@ref) and [`TCG`](@ref) for the
current epoch (`second` after J2000 and `fraction`) in seconds.

# Example

```jldoctest
julia> getoffset(TT, TCG, 0, 0.0)
0.5058332860211293
```
"""
function getoffset(::TerrestrialTime, ::GeocentricCoordinateTime, second, fraction)
    rate = LG_RATE / (1.0 - LG_RATE)
    dt = second - JD77_SEC + fraction
    return rate * dt
end

#######
# TCB #
#######

const LB_RATE = 1.550519768e-8

"""
    getoffset(TCB, TDB, second, fraction)

Return the linear offset between [`TCB`](@ref) and [`TDB`](@ref) for the
current epoch (`second` after J2000 and `fraction`) in seconds.

# Example

```jldoctest
julia> getoffset(TCB, TDB, 0, 0.0)
-11.253721593757295
```
"""
function getoffset(::BarycentricCoordinateTime, ::BarycentricDynamicalTime, second, fraction)
    dt = second - JD77_SEC + fraction
    return -LB_RATE * dt
end

"""
    getoffset(TDB, TCB, second, fraction)

Return the linear offset between [`TDB`](@ref) and [`TCB`](@ref) for the
current epoch (`second` after J2000 and `fraction`) in seconds.

# Example

```jldoctest
julia> getoffset(TDB, TCB, 0, 0.0)
11.253721768248475
```
"""
function getoffset(::BarycentricDynamicalTime, ::BarycentricCoordinateTime, second, fraction)
    rate = LB_RATE / (1.0 - LB_RATE)
    dt = second - JD77_SEC + fraction
    return rate * dt
end

#######
# UTC #
#######

function getleap(jd0)
    jd, frac = divrem(jd0, 1.0)
    jd += ifelse(frac >= 0.5, 1, -1) * 0.5
    drift0 = offset_tai_utc(jd)
    drift12 = offset_tai_utc(jd + 0.5)
    drift24 = offset_tai_utc(jd + 1.5)

    return drift24 - (2.0drift12 - drift0)
end

function getleap(::CoordinatedUniversalTime, date::Date)
    getleap(j2000(date) + value(J2000_TO_JULIAN))
end

getleap(::TimeScale, ::Date) = 0.0
getleap(ep::Epoch{CoordinatedUniversalTime}) = getleap(ep |> julian |> value)
getleap(::Epoch) = 0.0

function insideleap(jd0)
    jd1 = jd0 + 1 / SECONDS_PER_DAY
    o1 = offset_tai_utc(jd0)
    o2 = offset_tai_utc(jd1)
    return o1 != o2
end

"""
    insideleap(ep::UTCEpoch)

Check whether the UTC epoch `ep` is located during the introduction of a leap
second.

# Example

```jldoctest
julia> insideleap(UTCEpoch(2020, 3, 20))
false

julia> insideleap(UTCEpoch(2016, 12, 31, 23, 59, 60.5))
true
```
"""
insideleap(ep::Epoch{CoordinatedUniversalTime}) = insideleap(ep |> julian |> value)

insideleap(::Epoch) = false

"""
    getoffset(UTC, TAI, second, fraction)

Return the offset between [`UTC`](@ref) and [`TAI`](@ref) for the
current epoch (`second` after J2000 and `fraction`) in seconds.
For dates after `1972-01-01` this is the number of leap seconds.

# Example

```jldoctest
julia> getoffset(UTC, TAI, 0, 0.0)
32.0
```
"""
@inline function getoffset(::CoordinatedUniversalTime, ::InternationalAtomicTime,
                           second, fraction)
    jd = value(j2000(second, fraction) + J2000_TO_JULIAN)
    return -offset_utc_tai(jd)
end

"""
    getoffset(TAI, UTC, second, fraction)

Return the offset between [`TAI`](@ref) and [`UTC`](@ref) for the
current epoch (`second` after J2000 and `fraction`) in seconds.
For dates after `1972-01-01` this is the number of leap seconds.

# Example

```jldoctest
julia> getoffset(TAI, UTC, 0, 0.0)
-32.0
```
"""
@inline function getoffset(::InternationalAtomicTime, ::CoordinatedUniversalTime,
                           second, fraction)
    jd = value(j2000(second, fraction) + J2000_TO_JULIAN)
    return -offset_tai_utc(jd)
end

#######
# UT1 #
#######

"""
    getoffset(UTC, UT1, second, fraction[, eop])

Return the offset between [`UTC`](@ref) and [`UT1`](@ref) for the
current epoch (`second` after J2000 and `fraction`) in seconds.
Optionally, a custom Earth orientation data struct `eop` can be provided,
see [EarthOrientation.jl](https://github.com/JuliaAstro/EarthOrientation.jl).

# Example

```jldoctest; setup=:(AstroTime.load_test_eop())
julia> getoffset(UTC, UT1, 0, 0.0)
0.3550253556501879
```
"""
@inline function getoffset(::CoordinatedUniversalTime, ::UniversalTime,
                           second, fraction, eop=get(EOP_DATA))
    utc = value(j2000(second, fraction) + J2000_TO_JULIAN)
    return getΔUT1(eop, utc)
end

"""
    getoffset(UT1, UTC, second, fraction[, eop])

Return the offset between [`UT1`](@ref) and [`UTC`](@ref) for the
current epoch (`second` after J2000 and `fraction`) in seconds.
Optionally, a custom Earth orientation data struct `eop` can be provided,
see [EarthOrientation.jl](https://github.com/JuliaAstro/EarthOrientation.jl).

# Example

```jldoctest; setup=:(AstroTime.load_test_eop())
julia> getoffset(UT1, UTC, 0, 0.0)
-0.3550253592514352
```
"""
@inline function getoffset(::UniversalTime, ::CoordinatedUniversalTime,
                           second, fraction, eop=get(EOP_DATA))
    ut1 = value(j2000(second, fraction) + J2000_TO_JULIAN)
    utc = ut1 - getΔUT1(eop, ut1) / SECONDS_PER_DAY
    return -getΔUT1(eop, utc)
end

#######
# TDB #
#######

const k = 1.657e-3
const eb = 1.671e-2
const m₀ = 6.239996
const m₁ = 1.99096871e-7

"""
    getoffset(TT, TDB, second, fraction[, eop])

Return the offset between [`TT`](@ref) and [`TDB`](@ref) for the
current epoch (`second` after J2000 and `fraction`) in seconds.
This routine is accurate to ~40 microseconds over the interval 1900-2100.

!!! note
    An accurate transformation between TDB and TT depends on the
    trajectory of the observer. For two observers fixed on Earth's surface
    the quantity TDB-TT can differ by as much as ~4 microseconds. See
    [`here`](@ref getoffset(::TerrestrialTime, ::BarycentricDynamicalTime, second, fraction, elong, u, v)).

# Example

```jldoctest
julia> getoffset(TT, TDB, 0, 0.0)
-7.273677619130569e-5
```

### References ###

- [https://www.cv.nrao.edu/~rfisher/Ephemerides/times.html#TDB](https://www.cv.nrao.edu/~rfisher/Ephemerides/times.html#TDB)
- [Issue #26](https://github.com/JuliaAstro/AstroTime.jl/issues/26)
"""
@inline function getoffset(::TerrestrialTime, ::BarycentricDynamicalTime,
                           second, fraction)
    tt = fraction + second
    g = m₀ + m₁ * tt
    return k * sin(g + eb * sin(g))
end

"""
    getoffset(TDB, TT, second, fraction[, eop])

Return the offset between [`TDB`](@ref) and [`TT`](@ref) for the
current epoch (`second` after J2000 and `fraction`) in seconds.
This routine is accurate to ~40 microseconds over the interval 1900-2100.

!!! note
    An accurate transformation between TDB and TT depends on the
    trajectory of the observer. For two observers fixed on Earth's surface
    the quantity TDB-TT can differ by as much as ~4 microseconds. See
    [`here`](@ref getoffset(::BarycentricDynamicalTime, ::TerrestrialTime, second, fraction, elong, u, v)).

# Example

```jldoctest
julia> getoffset(TDB, TT, 0, 0.0)
7.273677616693264e-5
```

### References ###

- [https://www.cv.nrao.edu/~rfisher/Ephemerides/times.html#TDB](https://www.cv.nrao.edu/~rfisher/Ephemerides/times.html#TDB)
- [Issue #26](https://github.com/JuliaAstro/AstroTime.jl/issues/26)
"""
@inline function getoffset(::BarycentricDynamicalTime, ::TerrestrialTime,
                           second, fraction)
    tdb = fraction + second
    tt = tdb
    offset = 0.0
    for _ in 1:3
        g = m₀ + m₁ * tt
        offset = -k * sin(g + eb * sin(g))
        tt = tdb + offset
    end
    return offset
end

include(joinpath("constants", "tdb.jl"))

"""
    getoffset(TDB, TT, second, fraction[, eop])

Return the offset between [`TDB`](@ref) and [`TT`](@ref) for the
current epoch (`second` after J2000 and `fraction`) for an observer on earth
in seconds.

### Arguments ###

- `second`, `fraction`: Current epoch
- `elong`: Longitude (east positive, radians)
- `u`: Distance from Earth's spin axis (km)
- `v`: Distance north of equatorial plane (km)

# Example

```jldoctest
julia> getoffset(TDB, TT, 0, 0.0, π, 6371.0, 0.0)
9.928419814101153e-5
```

### References ###

- [ERFA](https://github.com/liberfa/erfa/blob/master/src/dtdb.c)
"""
function getoffset(::BarycentricDynamicalTime, ::TerrestrialTime, second, fraction, elong, u, v)
    tt = TTEpoch(second, fraction)
    ut = fractionofday(UT1Epoch(tt))
    t = (second + fraction) / (SECONDS_PER_CENTURY * 10.0)
    # Convert UT to local solar time in radians.
    tsol = mod(ut, 1.0) * 2π  + elong

    # FUNDAMENTAL ARGUMENTS:  Simon et al. 1994.
    # Combine time argument (millennia) with deg/arcsec factor.
    w = t / 3600.0
    # Sun Mean Longitude.
    elsun = deg2rad(mod(280.46645683 + 1296027711.03429 * w, 360.0))
    # Sun Mean Anomaly.
    emsun = deg2rad(mod(357.52910918 + 1295965810.481 * w, 360.0))
    # Mean Elongation of Moon from Sun.
    d = deg2rad(mod(297.85019547 + 16029616012.090 * w, 360.0))
    # Mean Longitude of Jupiter.
    elj = deg2rad(mod(34.35151874 + 109306899.89453 * w, 360.0))
    # Mean Longitude of Saturn.
    els = deg2rad(mod(50.07744430 + 44046398.47038 * w, 360.0))
    # TOPOCENTRIC TERMS:  Moyer 1981 and Murray 1983.
    wt = 0.00029e-10 * u * sin(tsol + elsun - els) +
        0.00100e-10 * u * sin(tsol - 2.0 * emsun) +
        0.00133e-10 * u * sin(tsol - d) +
        0.00133e-10 * u * sin(tsol + elsun - elj) -
        0.00229e-10 * u * sin(tsol + 2.0 * elsun + emsun) -
        0.02200e-10 * v * cos(elsun + emsun) +
        0.05312e-10 * u * sin(tsol - emsun) -
        0.13677e-10 * u * sin(tsol + 2.0 * elsun) -
        1.31840e-10 * v * cos(elsun) +
        3.17679e-10 * u * sin(tsol)

    # =====================
    # Fairhead et al. model
    # =====================

    # T**0
    w0 = 0.0
    for j in eachindex(fairhd0)
        @muladd w0 += fairhd0[j][1] * sin(fairhd0[j][2] * t + fairhd0[j][3])
    end
    # T**1
    w1 = 0.0
    for j in eachindex(fairhd1)
        @muladd w1 += fairhd1[j][1] * sin(fairhd1[j][2] * t + fairhd1[j][3])
    end
    # T**2
    w2 = 0.0
    for j in eachindex(fairhd2)
        @muladd w2 += fairhd2[j][1] * sin(fairhd2[j][2] * t + fairhd2[j][3])
    end
    # T**3
    w3 = 0.0
    for j in eachindex(fairhd3)
        @muladd w3 += fairhd3[j][1] * sin(fairhd3[j][2] * t + fairhd3[j][3])
    end
    # T**4
    w4 = 0.0
    for j in eachindex(fairhd4)
        @muladd w4 += fairhd4[j][1] * sin(fairhd4[j][2] * t + fairhd4[j][3])
    end

    # Multiply by powers of T and combine.
    wf = @evalpoly t w0 w1 w2 w3 w4

    # Adjustments to use JPL planetary masses instead of IAU.
    wj = 0.00065e-6 * sin(6069.776754 * t + 4.021194) +
        0.00033e-6 * sin(213.299095 * t + 5.543132) +
        (-0.00196e-6 * sin(6208.294251 * t + 5.696701)) +
        (-0.00173e-6 * sin(74.781599 * t + 2.435900)) +
        0.03638e-6 * t * t

    # TDB-TT in seconds.
    return -(wt + wf + wj)
end

"""
    getoffset(TT, TDB, second, fraction[, eop])

Return the offset between [`TT`](@ref) and [`TDB`](@ref) for the
current epoch (`second` after J2000 and `fraction`) for an observer on earth
in seconds.

### Arguments ###

- `second`, `fraction`: Current epoch
- `elong`: Longitude (east positive, radians)
- `u`: Distance from Earth's spin axis (km)
- `v`: Distance north of equatorial plane (km)

# Example

```jldoctest
julia> getoffset(TT, TDB, 0, 0.0, π, 6371.0, 0.0)
-9.92841981897215e-5
```

### References ###

- [ERFA](https://github.com/liberfa/erfa/blob/master/src/dtdb.c)
"""
@inline function getoffset(::TerrestrialTime, ::BarycentricDynamicalTime,
                           second, fraction, elong, u, v)
    tt1, tt2 = second, fraction
    tdb1, tdb2 = tt1, tt2
    offset = 0.0
    for _ in 1:3
        offset = -getoffset(TDB, TT, tdb1, tdb2, elong, u, v)
        tdb1, tdb2 = apply_offset(tt1, tt2, offset)
    end
    return offset
end


text:
struct AA
  name::String
  three_letter_code::String
  one_letter_code::String
  type::String
  polar::Bool
  hydrophobic::Bool
  charge::Int
end

const ALA = AA("Alanine"       , "ALA", "A",   "Aliphatic" , false, false,  0)
const ARG = AA("Arginine"      , "ARG", "R",   "Basic"     , true , false,  1)
const ASN = AA("Asparagine"    , "ASN", "N",   "Amide"     , true , false,  0)
const ASP = AA("Aspartic acid" , "ASP", "D",   "Acidic"    , true , false, -1)
const CYS = AA("Cysteine"      , "CYS", "C",   "Sulfuric"  , false, false,  0)
const GLN = AA("Glutamine"     , "GLN", "Q",   "Amide"     , true , false,  0)
const GLU = AA("Glutamic acid" , "GLU", "E",   "Acidic"    , true , false, -1)
const GLY = AA("Glycine"       , "GLY", "G",   "Aliphatic" , false, false,  0)
const HIS = AA("Histidine"     , "HIS", "H",   "Aromatic"  , true , false,  0)
const HSD = AA("Histidine"     , "HSD", "H",   "Aromatic"  , true , false,  1)
const HSE = AA("Histidine"     , "HSE", "H",   "Aromatic"  , true , false,  1)
const ILE = AA("Isoleucine"    , "ILE", "I",   "Aliphatic" , false, true,   0)
const LEU = AA("Leucine"       , "LEU", "L",   "Aliphatic" , false, true,   0)
const LYS = AA("Lysine"        , "LYS", "K",   "Basic"     , true , false,  1)
const MET = AA("Methionine"    , "MET", "M",   "Sulfuric"  , false, false,  0)
const PHE = AA("Phenylalanine" , "PHE", "F",   "Aromatic"  , false, true,   0)
const PRO = AA("Proline"       , "PRO", "P",   "Cyclic"    , false, false,  0)
const SER = AA("Serine"        , "SER", "S",   "Hydroxylic", true , false,  0)
const THR = AA("Threonine"     , "THR", "T",   "Hydroxylic", true , false,  0)
const TRP = AA("Tryptophan"    , "TRP", "W",   "Aromatic"  , false, true,   0)
const TYR = AA("Tyrosine"      , "TYR", "Y",   "Aromatic"  , true , false,  0)
const VAL = AA("Valine"        , "VAL", "V",   "Aliphatic" , false, true,   0)
const natural_aminoacids = [ ALA, ARG, ASN, ASP, CYS, GLN, GLU, GLY, HIS, HSD,
                             HSE, ILE, LEU, LYS, MET, PHE, PRO, SER, THR, TRP,
                             TYR, VAL ]


text:
function RLCore.Experiment(
    ::Val{:JuliaRL},
    ::Val{:DQN},
    ::Val{:SnakeGame},
    ::Nothing;
    save_dir = nothing,
    seed = 123,
)
    rng = Random.MersenneTwister(seed)

    SHAPE = (8, 8)
    inner_env = SnakeGameEnv(; action_style = FULL_ACTION_SET, shape = SHAPE, rng = rng)

    board_size = size(get_state(inner_env))
    N_FRAMES = 4

    env =
        inner_env |>
        StateOverriddenEnv(StackFrames(board_size..., N_FRAMES),) |>
        StateCachedEnv

    N_ACTIONS = length(get_actions(env))

    if isnothing(save_dir)
        t = Dates.format(now(), "yyyy_mm_dd_HH_MM_SS")
        save_dir = joinpath(pwd(), "checkpoints", "SnakeGame_$(t)")
    end

    lg = TBLogger(joinpath(save_dir, "tb_log"), min_level = Logging.Info)

    init = glorot_uniform(rng)

    update_freq = 4

    create_model() =
        Chain(
            x -> reshape(x, SHAPE..., :, size(x, ndims(x))),
            CrossCor(
                (3, 3),
                board_size[end] * N_FRAMES => 16,
                relu;
                stride = 1,
                pad = 1,
                init = init,
            ),
            CrossCor((3, 3), 16 => 32, relu; stride = 1, pad = 1, init = init),
            x -> reshape(x, :, size(x, ndims(x))),
            Dense(8 * 8 * 32, 256, relu; initW = init),
            Dense(256, N_ACTIONS; initW = init),
        ) |> cpu

    agent = Agent(
        policy = QBasedPolicy(
            learner = DQNLearner(
                approximator = NeuralNetworkApproximator(
                    model = create_model(),
                    optimizer = ADAM(0.001),
                ),
                target_approximator = NeuralNetworkApproximator(model = create_model()),
                update_freq = update_freq,
                γ = 0.99f0,
                update_horizon = 1,
                batch_size = 32,
                stack_size = nothing,
                min_replay_history = 20_000,
                loss_func = huber_loss,
                target_update_freq = 8_000,
                rng = rng,
            ),
            explorer = EpsilonGreedyExplorer(
                ϵ_init = 1.0,
                ϵ_stable = 0.01,
                decay_steps = 250_000,
                kind = :linear,
                rng = rng,
            ),
        ),
        trajectory = CircularCompactSALRTSALTrajectory(
            capacity = 500_000,
            state_type = Float32,
            state_size = (board_size..., N_FRAMES),
            legal_actions_mask_size = (N_ACTIONS,),
        ),
    )

    evaluation_result = []
    EVALUATION_FREQ = 100_000
    N_CHECKPOINTS = 3

    total_reward_per_episode = TotalRewardPerEpisode()
    time_per_step = TimePerStep()
    steps_per_episode = StepsPerEpisode()
    hook = ComposedHook(
        total_reward_per_episode,
        time_per_step,
        steps_per_episode,
        DoEveryNStep(update_freq) do t, agent, env
            with_logger(lg) do
                @info "training" loss = agent.policy.learner.loss log_step_increment =
                    update_freq
            end
        end,
        DoEveryNEpisode() do t, agent, env
            with_logger(lg) do
                @info "training" episode_length = steps_per_episode.steps[end] reward =
                    total_reward_per_episode.rewards[end] log_step_increment = 0
            end
        end,
    )

    N_TRAINING_STEPS = 1_000_000
    stop_condition = StopAfterStep(N_TRAINING_STEPS)
    description = """
    # Play Single Agent SnakeGame with DQN

    You can view the tensorboard logs with `tensorboard --logdir $(joinpath(save_dir, "tb_log"))`
    """
    Experiment(agent, env, stop_condition, hook, description)
end


text:
using Documenter, DemoCards
using ImageBase, ImageTransformations, ImageShow
using TestImages
using CoordinateTransformations, Interpolations, Rotations

# this is used to trigger artifact download and IO backend precompilation
testimage.(["cameraman", "lighthouse"])

examples, postprocess_cb, examples_assets = makedemos("examples")
format = Documenter.HTML(edit_link = "master",
                         prettyurls = get(ENV, "CI", nothing) == "true",
                         assets = [examples_assets,])

makedocs(modules  = [ImageTransformations, ImageBase],
         format   = format,
         sitename = "ImageTransformations",
         pages    = [
             "index.md",
             examples,
             "reference.md"])

postprocess_cb()

deploydocs(repo   = "github.com/JuliaImages/ImageTransformations.jl.git")


text:
using Base: @propagate_inbounds

const SplitView{T} = SubArray{T, 1, Vector{T}, Tuple{UnitRange{Int64}}, true}

struct SplitVector{T <: AbstractFloat} <: AbstractVector{T}

	#contiguous array of source data
	data::Vector{T}
	#array of data views of type Vector{T}, with indexing of
	#each view assumed to be over a contiguous range in the typedef
	views::Vector{SplitView{T}}
	split_by::AbstractConvexSet{T}   #sets that generated the views

	#constructor (composite set)
	function SplitVector{T}(x::Vector{T}, C::CompositeConvexSet{T}) where{T}

		#check for compatibility of vector to be split
		#and sum of all of the cone sizes
		@assert C.dim == length(x)

		views = Vector{SplitView{T}}(undef, num_subsets(C))

		# loop over the sets and create views
		sidx = 0
		for i = eachindex(C.sets)
			rng = (sidx + 1):(sidx + C.sets[i].dim)
			views[i] = view(x, rng)
			sidx += C.sets[i].dim
		end
		return new(x, views, C)
	end
	#constructor (non-composite)
	function SplitVector{T}(x::Vector{T}, C::AbstractConvexSet{T}) where{T}
		views = [view(x, 1:length(x))]
		return new(x, views, C)
	end
end

SplitVector(x::Vector{T}, C) where{T} = SplitVector{T}(x, C)

Base.size(A::SplitVector) = size(A.data)
Base.length(A::SplitVector) = length(A.data)
Base.IndexStyle(::Type{<:SplitVector}) = IndexLinear()
@propagate_inbounds Base.getindex(A::SplitVector, idx::Int) = getindex(A.data, idx)
@propagate_inbounds Base.setindex!(A::SplitVector, val, idx::Int) = setindex!(A.data, val, idx)

Base.iterate(A::SplitVector) = iterate(A.data)
Base.iterate(A::SplitVector,state) = iterate(A.data, state)
Base.firstindex(A::SplitVector) = 1
Base.lastindex(A::SplitVector)  = length(A.data)

Base.showarg(io::IO, A::SplitVector, toplevel) = print(io, typeof(A))