text:
@with_kw mutable struct DeepQLearningSolver{E<:ExplorationPolicy} <: Solver
    qnetwork::Any = nothing # intended to be a flux model
    learning_rate::Float32 = 1f-4
    max_steps::Int64 = 1000
    batch_size::Int64 = 32
    train_freq::Int64 = 4
    eval_freq::Int64 = 500
    target_update_freq::Int64 = 500
    num_ep_eval::Int64 = 100
    double_q::Bool = true
    dueling::Bool = true
    recurrence::Bool = false
    evaluation_policy::Any = basic_evaluation
    exploration_policy::E
    trace_length::Int64 = 40
    prioritized_replay::Bool = true
    prioritized_replay_alpha::Float32 = 0.6f0
    prioritized_replay_epsilon::Float32 = 1f-6
    prioritized_replay_beta::Float32 = 0.4f0
    buffer_size::Int64 = 1000
    max_episode_length::Int64 = 100
    train_start::Int64 = 200
    rng::AbstractRNG = MersenneTwister(0)
    logdir::Union{Nothing, String} = "log/"
    save_freq::Int64 = 3000
    log_freq::Int64 = 100
    verbose::Bool = true
end

function POMDPs.solve(solver::DeepQLearningSolver, problem::MDP)
    env = MDPCommonRLEnv{AbstractArray{Float32}}(problem) # ignores solver.rng because CommonRLEnv doesn't have rng support yet
    return solve(solver, env)
end

function POMDPs.solve(solver::DeepQLearningSolver, problem::POMDP)
    env = POMDPCommonRLEnv{AbstractArray{Float32}}(problem) # ignores solver.rng because CommonRLEnv doesn't have rng support yet
    return solve(solver, env)
end

function POMDPs.solve(solver::DeepQLearningSolver, env::AbstractEnv)
    action_map = collect(actions(env))
    action_indices = Dict(a=>i for (i, a) in enumerate(action_map))

    # check reccurence
    if isrecurrent(solver.qnetwork) && !solver.recurrence
        throw("DeepQLearningError: you passed in a recurrent model but recurrence is set to false")
    end
    replay = initialize_replay_buffer(solver, env, action_indices)
    if solver.dueling
        active_q = create_dueling_network(solver.qnetwork)
    else
        active_q = solver.qnetwork
    end
    policy = NNPolicy(env, active_q, action_map, length(obs_dimensions(env)))

    return dqn_train!(solver, env, policy, replay)
end

function dqn_train!(solver::DeepQLearningSolver, env::AbstractEnv, policy::AbstractNNPolicy, replay)
    if solver.logdir !== nothing
        logger = TBLogger(solver.logdir)
        solver.logdir = logger.logdir
    end
    active_q = getnetwork(policy) # shallow copy
    target_q = deepcopy(active_q)
    optimizer = ADAM(solver.learning_rate)
    # start training
    resetstate!(policy)
    reset!(env)
    obs = observe(env)
    done = false
    step = 0
    rtot = 0
    episode_rewards = Float64[0.0]
    episode_steps = Int64[]
    saved_mean_reward = -Inf
    scores_eval = -Inf
    model_saved = false
    eval_next = false
    save_next = false
    action_indices = Dict(a=>i for (i, a) in enumerate(actionmap(policy)))
    for t=1:solver.max_steps
        act = action(solver.exploration_policy, policy, t, obs)
        ai = action_indices[act]
        rew = act!(env, act)
        op = observe(env)
        done = terminated(env)
        exp = DQExperience(obs, ai, Float32(rew), op, done)
        if solver.recurrence
            add_exp!(replay, exp)
        elseif solver.prioritized_replay
            add_exp!(replay, exp, abs(exp.r))
        else
            add_exp!(replay, exp, 0f0)
        end
        obs = op
        step += 1
        episode_rewards[end] += rew
        if done || step >= solver.max_episode_length

            if eval_next # wait for episode to end before evaluating
                scores_eval, steps_eval, info_eval = evaluation(solver.evaluation_policy,
                policy, env,
                solver.num_ep_eval,
                solver.max_episode_length,
                solver.verbose)
                eval_next = false

                # only save after evaluation
                if save_next
                    model_saved, saved_mean_reward = save_model(solver, active_q, scores_eval, saved_mean_reward, model_saved)
                    save_next = false
                end

                if solver.logdir !== nothing
                    log_value(logger, "eval_reward", scores_eval, step = t)
                    log_value(logger, "eval_steps", steps_eval, step = t)
                    for (k, v) in info_eval
                        log_value(logger, k, v, step = t)
                    end
                end
            end

            reset!(env)
            obs = observe(env)
            resetstate!(policy)
            push!(episode_steps, step)
            push!(episode_rewards, 0.0)
            done = false
            step = 0
            rtot = 0
        end
        num_episodes = length(episode_rewards)
        avg100_reward = mean(episode_rewards[max(1, length(episode_rewards)-101):end])
        avg100_steps = mean(episode_steps[max(1, length(episode_steps)-101):end])
        if t%solver.train_freq == 0
            hs = hiddenstates(active_q)
            loss_val, grad_val = batch_train!(solver, env, policy, optimizer, target_q, replay)
            sethiddenstates!(active_q, hs)
        end

        if t%solver.target_update_freq == 0
            weights = Flux.params(active_q)
            Flux.loadparams!(target_q, weights)
        end

        if t % solver.eval_freq == 0
            eval_next = true
        end
        if t % solver.save_freq == 0
            save_next = true
        end

        if t % solver.log_freq == 0 && solver.logdir !== nothing
            nt = loginfo(solver.exploration_policy, t)
            for (k, v) in pairs(nt)
                log_value(logger, String(k), v, step=t)
            end
            if  solver.verbose
                @printf("%5d / %5d eps %0.3f |  avgR %1.3f | Loss %2.3e | Grad %2.3e | EvalR %1.3f \n",
                        t, solver.max_steps, nt[1], avg100_reward, loss_val, grad_val, scores_eval)
            end

            log_value(logger, "avg_reward", avg100_reward, step = t)
            log_value(logger, "loss", loss_val, step = t)
            log_value(logger, "grad_val", grad_val, step = t)
        end

    end # end training
    if model_saved
        if solver.verbose
            @printf("Restore model with eval reward %1.3f \n", saved_mean_reward)
            saved_model = BSON.load(joinpath(solver.logdir, "qnetwork.bson"))[:qnetwork]
            Flux.loadparams!(getnetwork(policy), saved_model)
        end
    end
    return policy
end

function initialize_replay_buffer(solver::DeepQLearningSolver, env::AbstractEnv, action_indices)
    # init and populate replay buffer
    if solver.recurrence
        replay = EpisodeReplayBuffer(env, solver.buffer_size, solver.batch_size, solver.trace_length)
    else
        replay = PrioritizedReplayBuffer(env, solver.buffer_size, solver.batch_size)
    end
    populate_replay_buffer!(replay, env, action_indices, max_pop=solver.train_start)
    return replay #XXX type unstable
end

function batch_train!(solver::DeepQLearningSolver,
                      env::AbstractEnv,
                      policy::AbstractNNPolicy,
                      optimizer,
                      target_q,
                      replay::PrioritizedReplayBuffer;
                      discount=default_discount(env)
                     )

    s_batch, a_batch, r_batch, sp_batch, done_batch, indices, importance_weights = sample(replay)

    active_q = getnetwork(policy)
    p = params(active_q)

    loss_val = nothing
    td_vals = nothing

    γ = convert(Float32, discount)
    if solver.double_q
        qp_values = active_q(sp_batch)
        target_q_values = target_q(sp_batch)
        best_a = [CartesianIndex(argmax(qp_values[:, i]), i) for i=1:solver.batch_size]
        q_sp_max = target_q_values[best_a]
    else
        q_sp_max = dropdims(maximum(target_q(sp_batch), dims=1), dims=1)
    end
    q_targets = r_batch .+ (1f0 .- done_batch) .* γ .* q_sp_max

    gs = Flux.gradient(p) do
        q_values = active_q(s_batch)
        q_sa = q_values[a_batch]
        td_vals = q_sa .- q_targets
        loss_val = sum(huber_loss, importance_weights.*td_vals)
        loss_val /= solver.batch_size
    end

    grad_norm = globalnorm(p, gs)
    Flux.Optimise.update!(optimizer, p, gs)

    if solver.prioritized_replay
        update_priorities!(replay, indices, td_vals)
    end

    return loss_val, grad_norm
end

# for RNNs
function batch_train!(solver::DeepQLearningSolver,
                      env::AbstractEnv,
                      policy::AbstractNNPolicy,
                      optimizer,
                      target_q,
                      replay::EpisodeReplayBuffer;
                      discount=default_discount(env)
                     )
    active_q = getnetwork(policy)
    s_batch, a_batch, r_batch, sp_batch, done_batch, trace_mask_batch = DeepQLearning.sample(replay)
    Flux.reset!(active_q)
    Flux.reset!(target_q)

    p = params(active_q)

    loss_val = nothing
    td_vals = nothing

    γ = convert(Float32, discount)
    q_targets = [zeros(Float32, solver.batch_size) for i=1:solver.trace_length]
    for i=1:solver.trace_length
        if solver.double_q
            qp_values = active_q(sp_batch[i])
            best_a = [CartesianIndex(argmax(qp_values[:, i]), i) for i=1:solver.batch_size]
            target_q_values = target_q(sp_batch[i])
            q_sp_max = target_q_values[best_a]
        else
            q_sp_max = dropdims(maximum(target_q(sp_batch[i]), dims=1), dims=1)
        end
        q_targets[i] .= r_batch[i] .+ (1f0 .- done_batch[i]) .* γ .* q_sp_max
    end

    Flux.reset!(active_q)

    gs = Flux.gradient(p) do
        loss_val = 0f0
        for i=1:solver.trace_length
            q_values = active_q(s_batch[i])
            q_sa = q_values[a_batch[i]]
            td_vals = q_sa .- q_targets[i]
            loss_val += sum(huber_loss, trace_mask_batch[i].*td_vals)/solver.batch_size
        end
        loss_val /= solver.trace_length
    end

    grad_norm = globalnorm(p, gs)
    Flux.Optimise.update!(optimizer, p, gs)
    return loss_val, grad_norm
end

function save_model(solver::DeepQLearningSolver, active_q, scores_eval::Float64, saved_mean_reward::Float64, model_saved::Bool)
    if scores_eval >= saved_mean_reward
        bson(joinpath(solver.logdir, "qnetwork.bson"), qnetwork=[w for w in params(active_q)])
        if solver.verbose
            @printf("Saving new model with eval reward %1.3f \n", scores_eval)
        end
        model_saved = true
        saved_mean_reward = scores_eval
    end
    return model_saved, saved_mean_reward
end

function restore_best_model(solver::DeepQLearningSolver, problem::MDP)
    env = convert(AbstractEnv, problem) # ignores solver.rng because CommonRLEnv doesn't have rng support yet
    restore_best_model(solver, env)
end

function restore_best_model(solver::DeepQLearningSolver, env::AbstractEnv)
    if solver.dueling
        active_q = create_dueling_network(solver.qnetwork)
    else
        active_q = solver.qnetwork
    end
    policy = NNPolicy(env, active_q, collect(actions(env)), length(obs_dimensions(env)))
    weights = BSON.load(solver.logdir*"qnetwork.bson")[:qnetwork]
    Flux.loadparams!(getnetwork(policy), weights)
    Flux.testmode!(getnetwork(policy))
    return policy
end

POMDPLinter.@POMDP_require POMDPs.solve(solver::DeepQLearningSolver, mdp::Union{MDP, POMDP}) begin
    P = typeof(mdp)
    S = POMDPs.statetype(P)
    A = POMDPs.actiontype(P)
    @req POMDPs.discount(::P)
    @req POMDPs.actions(::P)
    as = POMDPs.actions(mdp)
    @req length(::typeof(as))
    if isa(mdp, POMDP)
        O = obstype(mdp)
        @req POMDPs.convert_o(::Type{AbstractArray}, ::O, ::P)
    else
        @req POMDPs.convert_s(::Type{AbstractArray}, ::S, ::P)
    end
    @req POMDPs.reward(::P,::S,::A,::S)
end


text:
# Author: Philipp Witte, pwitte3@gatech.edu
# Date: January 2020

using InvertibleNetworks, LinearAlgebra, Test, Statistics

###############################################################################
# Test logdet implementation

# Input
nx = 4
ny = 4
nc = 3
batchsize = 1
X = rand(Float32, nx, ny, nc, batchsize)

# Actnorm and initialize
AN = ActNorm(nc; logdet=true)
AN.forward(X)

# Explicitely compute logdet of Jacobian through probing
# for small number of dimensions
J = zeros(Float32, Int(nx*ny*nc), Int(nx*ny*nc))
for i=1:nc
    count = 1
    for j=1:nx
        for k=1:ny
            E = zeros(Float32, nx, ny, nc, 1)
            E[k, j, i] = 1f0
            local Y = AN.forward(X)[1]
            J[:, (i-1)*nx*ny + count] = vec(AN.backward(E, Y)[1])
            count += 1
        end
    end
end
lgdet1 = log(abs(det(J)))
lgdet2 = AN.forward(X)[2]
@test isapprox((lgdet1 - lgdet2)/lgdet1, 0f0; atol=1f-6)

###############################################################################
# Initialization and invertibility

# Input
nx = 28
ny = 28
nc = 4
batchsize = 1
X = rand(Float32, nx, ny, nc, batchsize)

# Layer and initialization
AN = ActNorm(nc; logdet=false)
Y = AN.forward(X)

# Test initialization
@test isapprox(mean(Y), 0f0; atol=1f-6)
@test isapprox(var(Y), 1f0; atol=1f-3)

# Test invertibility
@test isapprox(norm(X - AN.inverse(AN.forward(X)))/norm(X), 0f0, atol=1f-6)
@test isapprox(norm(X - AN.forward(AN.inverse(X)))/norm(X), 0f0, atol=1f-6)

# Reversed layer (all combinations)
AN_rev = reverse(AN)

@test isapprox(norm(X - AN_rev.inverse(AN_rev.forward(X)))/norm(X), 0f0, atol=1f-6)
@test isapprox(norm(X - AN_rev.forward(AN_rev.inverse(X)))/norm(X), 0f0, atol=1f-6)

@test isapprox(norm(X - AN_rev.forward(AN.forward(X)))/norm(X), 0f0, atol=1f-6)
@test isapprox(norm(X - AN_rev.inverse(AN.inverse(X)))/norm(X), 0f0, atol=1f-6)

@test isapprox(norm(X - AN.forward(AN_rev.forward(X)))/norm(X), 0f0, atol=1f-6)
@test isapprox(norm(X - AN.inverse(AN_rev.inverse(X)))/norm(X), 0f0, atol=1f-6)

# Test with logdet enabled
AN = ActNorm(nc; logdet=true)
Y, lgdt = AN.forward(X)

# Test initialization
@test isapprox(mean(Y), 0f0; atol=1f-6)
@test isapprox(var(Y), 1f0; atol=1f-3)

# Test invertibility
@test isapprox(norm(X - AN.inverse(AN.forward(X)[1]))/norm(X), 0f0, atol=1f-6)
@test isapprox(norm(X - AN.forward(AN.inverse(X))[1])/norm(X), 0f0, atol=1f-6)

# Reversed layer (all combinations)
AN_rev = reverse(AN)

@test isapprox(norm(X - AN_rev.inverse(AN_rev.forward(X)[1]))/norm(X), 0f0, atol=1f-6)
@test isapprox(norm(X - AN_rev.forward(AN_rev.inverse(X))[1])/norm(X), 0f0, atol=1f-6)

@test isapprox(norm(X - AN_rev.forward(AN.forward(X)[1])[1])/norm(X), 0f0, atol=1f-6)
@test isapprox(norm(X - AN_rev.inverse(AN.inverse(X)))/norm(X), 0f0, atol=1f-6)

@test isapprox(norm(X - AN.forward(AN_rev.forward(X)[1])[1])/norm(X), 0f0, atol=1f-6)
@test isapprox(norm(X - AN.inverse(AN_rev.inverse(X)))/norm(X), 0f0, atol=1f-6)

###############################################################################
# Gradient Test

AN = ActNorm(nc; logdet=true)
X = randn(Float32, nx, ny, nc, batchsize)
X0 = randn(Float32, nx, ny, nc, batchsize)
dX = X - X0

# Forward pass
Y = AN.forward(X)[1]

function loss(AN, X, Y)

    # Forward pass
    Y_, lgdet = AN.forward(X)

    # Residual and function value
    ΔY = Y_ - Y
    f = .5f0/batchsize*norm(ΔY)^2
    AN.logdet == true && (f -= lgdet)

    # Back propagation
    ΔX, X_ = AN.backward(ΔY./batchsize, Y)

    # Check invertibility
    isapprox(norm(X - X_)/norm(X), 0f0, atol=1f-6)

    return f, ΔX, AN.s.grad, AN.b.grad
end

# Gradient test for X
maxiter = 6
print("\nGradient test actnorm\n")
f0, ΔX = loss(AN, X0, Y)[1:2]
h = .1f0
err1 = zeros(Float32, maxiter)
err2 = zeros(Float32, maxiter)
for j=1:maxiter
    f = loss(AN, X0 + h*dX, Y)[1]
    err1[j] = abs(f - f0)
    err2[j] = abs(f - f0 - h*dot(dX, ΔX))
    print(err1[j], "; ", err2[j], "\n")
    global h = h/2f0
end

@test isapprox(err1[end] / (err1[1]/2^(maxiter-1)), 1f0; atol=1f1)
@test isapprox(err2[end] / (err2[1]/4^(maxiter-1)), 1f0; atol=1f1)

# Gradient test for parameters
AN0 = ActNorm(nc; logdet=true); AN0.forward(randn(Float32, nx, ny, nc, batchsize))
AN_ini = deepcopy(AN0)
ds = AN.s.data - AN0.s.data
db = AN.b.data - AN0.b.data
maxiter = 6
print("\nGradient test actnorm\n")
f0, ΔX, Δs, Δb = loss(AN0, X, Y)
h = 1f0
err3 = zeros(Float32, maxiter)
err4 = zeros(Float32, maxiter)
for j=1:maxiter
    AN0.s.data = AN_ini.s.data + h*ds
    AN0.b.data = AN_ini.b.data + h*db
    f = loss(AN0, X, Y)[1]
    err3[j] = abs(f - f0)
    err4[j] = abs(f - f0 - h*dot(ds, Δs) - h*dot(db, Δb))
    print(err3[j], "; ", err4[j], "\n")
    global h = h/2f0
end

@test isapprox(err3[end] / (err3[1]/2^(maxiter-1)), 1f0; atol=1f1)
@test isapprox(err4[end] / (err4[1]/4^(maxiter-1)), 1f0; atol=1f1)

###############################################################################
# Gradient Test reversed layer

AN = reverse(ActNorm(nc; logdet=true))
X = randn(Float32, nx, ny, nc, batchsize)
X0 = randn(Float32, nx, ny, nc, batchsize)
dX = X - X0

# Forward pass
Y = AN.forward(X)[1]

function loss(AN, X, Y)

    # Forward pass
    Y_, lgdet = AN.forward(X)

    # Residual and function value
    ΔY = Y_ - Y
    f = .5f0/batchsize*norm(ΔY)^2
    AN.logdet == true && (f -= lgdet)

    # Back propagation
    ΔX, X_ = AN.backward(ΔY./batchsize, Y)

    # Check invertibility
    isapprox(norm(X - X_)/norm(X), 0f0, atol=1f-6)

    return f, ΔX, AN.s.grad, AN.b.grad
end

# Gradient test for X
maxiter = 6
print("\nGradient test actnorm reveres\n")
f0, ΔX = loss(AN, X0, Y)[1:2]
h = .1f0
err5 = zeros(Float32, maxiter)
err6 = zeros(Float32, maxiter)
for j=1:maxiter
    f = loss(AN, X0 + h*dX, Y)[1]
    err5[j] = abs(f - f0)
    err6[j] = abs(f - f0 - h*dot(dX, ΔX))
    print(err5[j], "; ", err6[j], "\n")
    global h = h/2f0
end

@test isapprox(err5[end] / (err5[1]/2^(maxiter-1)), 1f0; atol=1f1)
@test isapprox(err6[end] / (err6[1]/4^(maxiter-1)), 1f0; atol=1f1)

# Gradient test for parameters
AN0 = reverse(ActNorm(nc; logdet=true))
AN0.forward(randn(Float32, nx, ny, nc, batchsize))
AN_ini = deepcopy(AN0)

ds = randn(Float32, size(AN.s.data))
db = randn(Float32, size(AN.b.data))
maxiter = 6
print("\nGradient test actnorm reverse\n")
f0, ΔX, Δs, Δb = loss(AN0, X, Y)
h = 1f0
err7 = zeros(Float32, maxiter)
err8 = zeros(Float32, maxiter)
for j=1:maxiter
    AN0.s.data = AN_ini.s.data + h*ds
    AN0.b.data = AN_ini.b.data + h*db
    f = loss(AN0, X, Y)[1]
    err7[j] = abs(f - f0)
    err8[j] = abs(f - f0 - h*dot(ds, Δs) - h*dot(db, Δb))
    print(err7[j], "; ", err8[j], "\n")
    global h = h/2f0
end

@test isapprox(err7[end] / (err7[1]/2^(maxiter-1)), 1f0; atol=1f1)
@test isapprox(err8[end] / (err8[1]/4^(maxiter-1)), 1f0; atol=1f1)

ext: jl


text:
#
# Binary Operators
#

for fun in (:+, :-, :*, :/, :^, :min, :max)
    @eval begin
        # Konst eager algebra
        @inline function Base.$fun(k1::Konst{T1}, k2::Konst{T2}) where {T1<:Number, T2<:Number}
            return Konst(($fun)(k1.val, k2.val))
        end

        # Observables Binary operators
        @inline Base.$fun(o1::Observable, o2::Observable) = LiftObs2($fun, o1, o2)

        # Lift Number to Konst
        @inline Base.$fun(n::Number, o::Observable) = ($fun)(Konst(n), o)
        @inline Base.$fun(o::Observable, n::Number) = ($fun)(o, Konst(n))
    end
end

#
# Special functions
#

Base.round(k::Konst; digits=0) = Konst(round(k.val; digits=digits))
Base.trunc(k::Konst; digits=0) = Konst(trunc(k.val; digits=digits))
Base.iszero(k::Konst) = iszero(k.val)
Base.sqrt(o::Observable) = LiftObs(sqrt, o)
Base.exp(o::Observable) = LiftObs(exp, o)
Base.log(o::Observable) = LiftObs(log, o)
Base.exp(o::LiftObs{typeof(log)}) = o.o
Base.log(o::LiftObs{typeof(exp)}) = o.o

"""
    iszerokonst(observable) :: Bool

Returns `true` if `observable` is of type `Konst` and has a zero value.
Returns `false` otherwise.
"""
iszerokonst(k::Konst) = iszero(k)
iszerokonst(k::Observable) = false

#
# Associative ops
#

@inline function associative_op(f::F, k::Konst{N1}, o1::Konst{N2}, o2::Observable) where {F<:Function, N1<:Number, N2<:Number}
    return f(f(k, o1), o2)
end

@inline function associative_op(f::F, k::Konst{N1}, o1::Observable, o2::Konst{N2}) where {F<:Function, N1<:Number, N2<:Number}
    return associative_op(f, k, o2, o1)
end

@inline function associative_op(f::F, k::Konst, o1::Observable, o2::Observable) where {F<:Function}
    return LiftObs2(f, k, LiftObs2(f, o1, o2))
end

# to the left: k1 * (k2 * obs) == (k1 * k2) * obs
# to the left: k1 * (obs * k2) == (k1 * k2) * obs
# to the right: (k1 * obs) * k2 == (k1 * k2) * obs
# to the right: (obs * k1) * k2 == (k1 * k2) * obs
for fun in (:+, :*)
    @eval begin
        @inline function Base.$fun(k::Konst{N}, lo::LiftObs2{typeof($fun)}) where {N<:Number}
            associative_op($fun, k, lo.o1, lo.o2)
        end

        @inline function Base.$fun(lo::LiftObs2{typeof($fun)}, k::Konst{N}) where {N<:Number}
            associative_op($fun, k, lo.o1, lo.o2)
        end
    end
end

#
# algebra for minus (-)
#

@inline Base.:(-)(k::Konst{T}) where {T<:Number} = Konst(-k.val) # Konst eager algebra
@inline Base.:(-)(o::Observable) = LiftObs(-, o)

# -(a + b) == (-a) + (-b)
@inline function Base.:(-)(lo::LiftObs2{typeof(+)})
    return LiftObs2(+, -lo.o1, -lo.o2)
end

# -(a - b) == -a + b
@inline function Base.:(-)(lo::LiftObs2{typeof(-)})
    return LiftObs2(+, -lo.o1, lo.o2)
end

# k - (a +- b) == k + (-(a +- b))
@inline function Base.:(-)(k::Konst{N}, lo::LiftObs2{F}) where {N<:Number, F<:Union{typeof(+), typeof(-)}}
    return k + (-lo)
end

@inline function Base.:(-)(lo::LiftObs2{F}, k::Konst{N}) where {N<:Number, F<:Union{typeof(+), typeof(-)}}
    return lo + (-k)
end

# (k1 - obs) + k2 == (k1 + k2) - obs
# (obs - k1) + k2 == (k2 - k1) + obs
@inline function Base.:(+)(lo::LiftObs2{typeof(-)}, k::Konst{N}) where {N<:Number}
    associative_op(+, k, lo.o1, -lo.o2)
end

# k2 + (k1 - obs) == (k1 + k2) + (-obs)
@inline function Base.:(+)(k::Konst{N}, lo::LiftObs2{typeof(-)}) where {N<:Number}
    return lo + k
end


text:
using Revise
using Koala
using KoalaRidge
using Base.Test
using DataFrames

## SYNTHETIC DATA TEST

# Define some linear, noise-free, synthetic data:
const bias = -42.0
const coefficients = Float64[1, 3, 7]
const A = randn(1000, 3)
const X = DataFrame(A)
const y = A*coefficients

# Train model on all data with no regularization and no
# standardization of target:
const ridge = RidgeRegressor(lambda=0.0)
ty = default_transformer_y(ridge)
ty.standardize = false
ridgeM = Machine(ridge, X, y, eachindex(y), transformer_y=ty)
fit!(ridgeM, eachindex(y))

# Training error:
E = err(ridgeM, eachindex(y))
@test E < 1e-10

# Get the true bias?
@test abs(ridgeM.predictor.bias) < 1e-10
@test norm(ridgeM.predictor.coefficients - coefficients) < 1e-10

# Train with standardization of target:
ridgeM1 = Machine(ridge, X, y, eachindex(y))
fit!(ridgeM1, eachindex(y))

# Get same predictions?
@test isapprox(predict(ridgeM, X), predict(ridgeM1, X))

## TEST OF OTHER METHODS ON REAL DATA

# Load some data and define train/test rows:
const X, y = load_ames();
y = log.(y) # log of the SalePrice
const train, test = split(eachindex(y), 0.7); # 70:30 split

# Instantiate a model:
ridge = RidgeRegressor(lambda=0.1)
ty = default_transformer_y(ridge)
ty.boxcox = false
showall(ty)
showall(ridge)

# Build a machine:
ridgeM = Machine(ridge, X, y, train, transformer_y=ty)

fit!(ridgeM, train)
showall(ridgeM)

# tune lambda using cross-validation
lambdas, rmserrors = @curve λ logspace(-3,1,100) begin
    ridge.lambda = λ
    mean(cv(ridgeM, train, n_folds=9, verbosity=0))
end

# set lambda to the optimal value and do final train:
ridge.lambda = lambdas[indmin(rmserrors)]
fit!(ridgeM, train)

score = err(ridgeM, test)
println("score = $score")
@test score > 0.12 && score < 0.14


text:
using UnicodeFun

@testset "Parity" begin
    @testset "Construction" begin
        @test p"even".p
        @test !p"odd".p
        @test convert(Parity, 1) == p"even"
        @test convert(Parity, -1) == p"odd"
        @test_throws ArgumentError convert(Parity, 5)
    end

    @testset "Boolean properties" begin
        @test iseven(p"even")
        @test !isodd(p"even")
        @test !iseven(p"odd")
        @test isodd(p"odd")
        @test p"odd" < p"even"
        @test p"even" ≥ p"odd"
    end

    @testset "Arithmetic" begin
        @test iseven(p"even"*p"even")
        @test isodd(p"even"*p"odd")
        @test isodd(p"odd"*p"even")
        @test iseven(p"odd"*p"odd")

        @test iseven(p"even"^0)
        @test iseven(p"even"^1)
        @test iseven(p"even"^2)
        @test iseven(p"even"^3)
        @test iseven(p"odd"^0)
        @test isodd(p"odd"^1)
        @test iseven(p"odd"^2)
        @test isodd(p"odd"^3)

        @test isodd(-p"even")
        @test iseven(-p"odd")
    end

    @testset "Pretty-printing" begin
        @test "$(p"even")" == "even"
        @test "$(p"odd")" == "odd"
        @test to_superscript(p"even") == ""
        @test to_superscript(p"odd") == "ᵒ"
    end
end


text:
module EFT
export two_sum, fast_two_sum, two_prod

using VectorizationBase: ifelse, AbstractSIMD

include("errorfree.jl")

end


text:
function init_V_gauss_G(strf, Ω, Ns, G2, A, alpha)
    Npoints = size(strf)[1]
    Nspecies = size(strf)[2]
    Vg = zeros(ComplexF64,Npoints)
    V = zeros(Float64,Npoints)
    for isp = 1:Nspecies
        pf = -A[isp]*(pi/alpha[isp])^1.5/Ω
        @printf("Gaussian parameters:\n")
        @printf("A[%2d]     = %18.10f\n", isp, A[isp])
        @printf("alpha[%2d] = %18.10f\n", isp, alpha[isp])
        for ig = 1:Npoints
            Vg[ig] = pf*exp(-0.25*G2[ig]/alpha[isp])*strf[ig,isp]
        end
        V = V + real( G_to_R(Ns,Vg) )*Npoints
    end
    return V
end


text:
using Plots, SparseArrays, JLD2, LinearAlgebra, Wavelets, MultiscaleGraphSignalTransforms

#img = load("test\\8.jpg")
#heatmap(img,ratio=1, yaxis =:flip, axis = false, color = :gray)

JLD2.@load "../spie_data.jld2"
matrix = vars["barbara"][1:4:512,1:4:512]

heatmap(matrix,ratio=1, yaxis =:flip, axis = false, color = :grays, clim = (0,1))

matrix_gt = matrix[:]

sigma = 0.01
m, n = size(matrix)

############### Get the neighbors for affinity matrix computing
r = 5; # specify radius of neighbors
l = 2*r + 1

temp_x, temp_y = fill(1,l^2), fill(1,l^2)
temp_xy = CartesianIndices((1:l,1:l))
for i = 1:l^2
    temp_x[i], temp_y[i] = temp_xy[i][1], temp_xy[i][2]
end
temp_ind = ((temp_x .- (r + 1)).^2 + (temp_y .- (r + 1)).^2).<= r^2
neighbor_x = temp_x[temp_ind] .- (r + 1)
neighbor_y = temp_y[temp_ind] .- (r + 1)
# for any index(x,y), (x + neighbor_x, y - neigbor_y) are the neighbors to calculate affinity

################ Create affinity matrix
W = fill(0., (m*n,m*n))
for i = 1:m*n
    cur = CartesianIndices((m,n))[i]
    for j = 1:length(neighbor_x)
        if 1 <= cur[1] + neighbor_x[j] <= m && 1 <= cur[2] + neighbor_y[j] <= n
            tempd = LinearIndices((m,n))[cur[1] + neighbor_x[j], cur[2] + neighbor_y[j]]
            W[i,tempd] = exp(-(matrix_gt[i] - matrix_gt[tempd])^2/sigma)
        end
    end
end
W = sparse(W)

############### Preprocess to get G and GP
G = GraphSig(W, f = reshape(matrix,(length(matrix_gt),1)))
GP = partition_tree_fiedler(G)
dmatrix = ghwt_analysis!(G, GP=GP)

G_tmp = Glevel(G,GP,4)
heatmap(reshape(G_tmp.f[:],(128,128)))
############# Construct or search the specific basis
############# Haar
BS_haar = bs_haar(GP)
dvec_haar = dmatrix2dvec(dmatrix, GP, BS_haar)

############# eGHWT
dvec_eghwt, BS_eghwt = ghwt_tf_bestbasis(dmatrix, GP)

############# Walsh
BS_walsh = bs_walsh(GP)
dvec_walsh = dmatrix2dvec(dmatrix, GP, BS_walsh)

############# GHWT_c2f
dvec_c2f, BS_c2f = ghwt_c2f_bestbasis(dmatrix, GP)

############# GHWT_f2c
dvec_f2c, BS_f2c = ghwt_f2c_bestbasis(dmatrix, GP)

###################################################################################
############## Visuallization of Top basis vector (multiple)#######################
###################################################################################
### visualization fo the first 12 vectors
### function to plot top 12 vectors
function top_vectors_plot(dvec::Array{Float64, 2}, BS::BasisSpec, GP::GraphPart; clims::Tuple{Float64,Float64} = (-0.02,0.02))
    sorted_ind = sortperm(abs.(dvec[:]), rev = true);
    plot(12,layout=(3,4),framestyle=:none, legend=false)
    for i=1:12
        dvecT = fill(0., size(dvec))
        #dvecT[sorted_ind[i]] = dvec_ghwt[sorted_ind[i]]
        dvecT[sorted_ind[i]] = 1
        f = ghwt_synthesis(dvecT, GP, BS)
        #print((maximum(f), minimum(f)))
        heatmap!(reshape(f, size(matrix)), subplot=i, ratio=1, yaxis=:flip, axis=false, color = :balance, clims = clims)
    end
    current()
end

### haar
top_vectors_plot(dvec_haar, BS_haar, GP)
### walsh
top_vectors_plot(dvec_walsh, BS_walsh, GP)
### c2f
top_vectors_plot(dvec_c2f, BS_c2f, GP)
### f2c
top_vectors_plot(dvec_f2c, BS_f2c, GP)
### eghwt
top_vectors_plot(dvec_eghwt, BS_eghwt, GP)

###################################################################################
############## Visuallization of Top basis vector (single)#########################
###################################################################################
### function to plot the first i-th vector
### Note that here only one vector is plotted instead of multiple vectors.
function top_vector_plot(dvec::Array{Float64, 2}, BS::BasisSpec, GP::GraphPart, G::GraphSig, i::Int64; clims::Tuple{Float64,Float64} = (-0.02,0.02))
    sorted_ind = sortperm(abs.(dvec[:]), rev = true);
    dvecT = fill(0., size(dvec))
    #dvecT[sorted_ind[i]] = dvec_ghwt[sorted_ind[i]]
    dvecT[sorted_ind[i]] = 1
    f, GS = ghwt_synthesis(dvecT, GP, BS, G)
    heatmap(reshape(f, size(matrix)), ratio=1, yaxis=:flip, axis=false, color = :grays, clims = clims)
end
########################## plot top i-th vector
i = 2
### haar
top_vector_plot(dvec_haar, BS_haar, GP, G, i)
### walsh
top_vector_plot(dvec_walsh, BS_walsh, GP, G, i)
### c2f
top_vector_plot(dvec_c2f, BS_c2f, GP, G, i)
### f2c
top_vector_plot(dvec_f2c, BS_f2c, GP, G, i)
### eghwt
top_vector_plot(dvec_eghwt, BS_eghwt, GP, G, i)

#######################################################################################
### generate figure used for spie paper
##########################################################################################
# Figure 8b
sorted_ind = sortperm(abs.(dvec_eghwt[:]), rev = true)
plot(4, layout=(2,2), framestyle=:none, legend=false)

for i = 1:4
    dvec_T = fill(0., size(dvec_eghwt))
    dvec_T[sorted_ind[i+1]] = 1
    f = ghwt_synthesis(dvec_T, GP, BS_eghwt)
    heatmap!(reshape(f, size(matrix)), subplot=i, ratio=1, yaxis=:flip, axis=false, color = :grays, clims = (-0.01,0.01))
end
current()

savefig("barbara_top4.pdf")

#########################
function approx_error(DVEC::Array{Array{Float64,1},1})
    plot(xaxis = "Fraction of Coefficients Retained", yaxis = "Relative Approximation Error")
    frac = 0:0.01:0.3
    T = ["Classical Haar transform", "eGHWT Haar basis", "eGHWT best basis"]
    L = [(:dashdot,:orange),(:dashdot,:blue),(:solid, :black)]
    for i = 1:3
        dvec = DVEC[i]
        N = length(dvec)
        dvec_norm = norm(dvec,2)
        dvec_sort = sort(dvec.^2, rev = true)

        er = fill(0., length(frac))

        for j = 1:length(frac)
            p = Int64(floor(frac[j]*N))
            er[j] = sqrt(dvec_norm^2 - sum(dvec_sort[1:p]))/dvec_norm
        end

        plot!(frac, er, yaxis=:log, xlims = (0.,0.3), label = T[i], line = L[i])
        print(er[26])
    end
end

#########################
function approx_error2(DVEC::Array{Array{Float64,1},1})
    plot(xaxis = "Fraction of Coefficients Retained", yaxis = "Relative Approximation Error")
    frac = 0.3
    T = ["Classical Haar transform", "eGHWT Haar basis", "eGHWT best basis"]
    L = [(:dashdot,:orange),(:dashdot,:blue),(:solid, :black)]
    for i = 1:3
        dvec = DVEC[i]
        N = length(dvec)
        dvec_norm = norm(dvec,2)
        dvec_sort = sort(dvec.^2) # the smallest first
        er = sqrt.(reverse(cumsum(dvec_sort)))/dvec_norm # this is the relative L^2 error of the whole thing, i.e., its length is N
        p = Int64(floor(frac*N)) + 1 # upper limit
        plot!(frac*(0:(p-1))/(p-1), er[1:p], yaxis=:log, xlims = (0.,frac), label = T[i], line = L[i])
    end
end

dvec_classichaar = dwt(matrix, wavelet(WT.haar))

### Figure 8(a)
#approx_error([dvec_classichaar[:], dvec_haar[:], dvec_eghwt[:]])
approx_error2([dvec_classichaar[:], dvec_haar[:], dvec_eghwt[:]])
current()


text:
using Test

TEST_PLOT_OUTPUT = !isempty(get(ENV, "CES_TEST_PLOT_OUTPUT", ""))

if TEST_PLOT_OUTPUT
    using Plots
end

function include_test(_module)
    println("Starting tests for $_module")
    t = @elapsed include(joinpath(_module, "runtests.jl"))
    println("Completed tests for $_module, $(round(Int, t)) seconds elapsed")
    return nothing
end

@testset "EnsembleKalmanProcesses" begin
    all_tests = isempty(ARGS) || "all" in ARGS ? true : false

    function has_submodule(sm)
        any(ARGS) do a
            a == sm && return true
	    first(split(a, '/')) == sm && return true
	    return false
        end
    end

    for submodule in ["DataStorage",
                      "ParameterDistribution",
                      "Observations",
                      "EnsembleKalmanProcessModule"]

        if all_tests || has_submodule(submodule) || "EnsembleKalmanProcesses" in ARGS
            include_test(submodule)
        end
    end
end


text:
@testset "ee00" begin
    epsa =  0.4090789763356509900
    dpsi = -0.9630909107115582393e-5
    ee = ee00(2400000.5, 53736.0, epsa, dpsi)
    @test ee ≈ -0.8834193235367965479e-5 atol=1e-18
end

@testset "ee00a" begin
    ee = ee00a(2400000.5, 53736.0)
    @test ee ≈ -0.8834192459222588227e-5 atol=1e-18
end

@testset "ee00b" begin
    ee = ee00b(2400000.5, 53736.0)
    @test ee ≈ -0.8835700060003032831e-5 atol=1e-18
end

@testset "ee06a" begin
    ee = ee06a(2400000.5, 53736.0)
    @test ee ≈ -0.8834195072043790156e-5 atol=1e-15
end

@testset "eect00" begin
    ct = eect00(2400000.5, 53736.0)
    @test ct ≈ 0.2046085004885125264e-8 atol=1e-20
end

@testset "eo06a" begin
    eo = eo06a(2400000.5, 53736.0)
    @test eo ≈ -0.1332882371941833644e-2 atol=1e-15
end

@testset "eors" begin
    r = [0.9999989440476103608 -0.1332881761240011518e-2 -0.5790767434730085097e-3;
         0.1332858254308954453e-2 0.9999991109044505944 -0.4097782710401555759e-4;
         0.5791308472168153320e-3 0.4020595661593994396e-4 0.9999998314954572365]
    s = -0.1220040848472271978e-7
    eo = eors(r, s)
    @test eo ≈ -0.1332882715130744606e-2 atol=1e-15
    @test_throws ArgumentError eors(r[1:2,:], s)
end

@testset "eform" begin
    a, f = eform(ERFA.WGS84)
    @test a ≈ 6378137.0 atol=1e-10
    @test f ≈ 0.0033528106647474807 atol=1e-18
    a, f = eform(ERFA.GRS80)
    @test a ≈ 6378137.0 atol=1e-10
    @test f ≈ 0.0033528106811823189 atol=1e-18
    a, f = eform(ERFA.WGS72)
    @test a ≈ 6378135.0 atol=1e-10
    @test f ≈ 0.0033527794541675049 atol=1e-18
end

@testset "epb" begin
    b = epb(2415019.8135, 30103.18648)
    @test b ≈ 1982.418424159278580 atol=1e-12
end

@testset "epb2jd" begin
    dj0, dj1 = epb2jd(1957.3)
    @test dj0 ≈ 2400000.5 atol=1e-9
    @test dj1 ≈ 35948.1915101513 atol=1e-9
end

@testset "epj" begin
    j = epj(2451545, -7392.5)
    @test j ≈ 1979.760438056125941 atol=1e-12
end

@testset "epj2jd" begin
    dj0, dj1 = epj2jd(1996.8)
    @test dj0 ≈ 2400000.5 atol=1e-9
    @test dj1 ≈ 50375.7 atol=1e-9
end

@testset "epv00" begin
    pvh, pvb = epv00(2400000.5, 53411.52501161)
    @test pvh[1][1] ≈ -0.7757238809297706813 atol=1e-14
    @test pvh[1][2] ≈ 0.5598052241363340596 atol=1e-14
    @test pvh[1][3] ≈ 0.2426998466481686993 atol=1e-14
    @test pvh[2][1] ≈ -0.1091891824147313846e-1 atol=1e-15
    @test pvh[2][2] ≈ -0.1247187268440845008e-1 atol=1e-15
    @test pvh[2][3] ≈ -0.5407569418065039061e-2 atol=1e-15
    @test pvb[1][1] ≈ -0.7714104440491111971 atol=1e-14
    @test pvb[1][2] ≈ 0.5598412061824171323 atol=1e-14
    @test pvb[1][3] ≈ 0.2425996277722452400 atol=1e-14
    @test pvb[2][1] ≈ -0.1091874268116823295e-1 atol=1e-15
    @test pvb[2][2] ≈ -0.1246525461732861538e-1 atol=1e-15
    @test pvb[2][3] ≈ -0.5404773180966231279e-2 atol=1e-15
    @test_logs (:warn,) epv00(2.5245935e6, 0.0)
end

@testset "eqeq94" begin
    ee = eqeq94(2400000.5, 41234.0)
    @test ee ≈ 0.5357758254609256894e-4 atol=1e-17
end

@testset "era00" begin
    era = era00(2400000.5, 54388.0)
    @test era ≈ 0.4022837240028158102 atol=1e-12
end

@testset "eceq06" begin
    dr, dd = eceq06(2456165.5, 0.401182685, 5.1, -0.9)
    @test dr ≈ 5.533459733613627767 atol=1e-14
    @test dd ≈ -1.246542932554480576 atol=1e-14
end

@testset "eqec06" begin
    dl, db = eqec06(1234.5, 2440000.5, 1.234, 0.987)
    @test dl ≈ 1.342509918994654619 atol=1e-14
    @test db ≈ 0.5926215259704608132 atol=1e-14
end

@testset "ecm06" begin
    rm = ecm06(2456165.5, 0.401182685)
    @test rm[1,1] ≈ 0.9999952427708701137 atol=1e-14
    @test rm[1,2] ≈ -0.2829062057663042347e-2 atol=1e-14
    @test rm[1,3] ≈ -0.1229163741100017629e-2 atol=1e-14
    @test rm[2,1] ≈ 0.3084546876908653562e-2 atol=1e-14
    @test rm[2,2] ≈ 0.9174891871550392514 atol=1e-14
    @test rm[2,3] ≈ 0.3977487611849338124 atol=1e-14
    @test rm[3,1] ≈ 0.2488512951527405928e-5 atol=1e-14
    @test rm[3,2] ≈ -0.3977506604161195467 atol=1e-14
    @test rm[3,3] ≈ 0.9174935488232863071 atol=1e-14
end


text:
# Tests for Column and LazyColumn types, and associated methods

@testset "Column" begin
    @testset "Constructors" begin
        @test_throws DomainError Column(-1, [])

        a = Column(0, [1,1,1,1])
        b = Column(a)
        c = Column{Integer}(0, 4)
        d = Column(0, 4)
        @test a == b == c == d

        b.data[2] = 4
        @test a != b

        a = Column(3, [1.0, 4.0, 10.0, 20.0])
        b = Column(a)
        c = Column{Float64}(3, 4)
        d = Column(3, 4) # different data type
        @test a == b == c == d
    end

    @testset "Accessors" begin
        a = Column(5,4)
        @test colnumber(a) == 5
        @test value(a) == [1, 6, 21, 56]
    end

    @testset "Array functions" begin
        a = Column(7, 8)
        b = Column{Float64}(7, 8)
        @test size(a) == size(b) == (8,)

        @testset "Data access $i" for i ∈ 7:14
            @test a[i] == b[i]
        end

        @test_throws BoundsError a[6]
        @test_throws BoundsError a[15]

        @test firstindex(a) == firstindex(b) == 7
        @test lastindex(a) == lastindex(b) == 14
    end

    @testset "Checks" begin
        a = Column(0, 5)
        b = Column(4, 5)
        c = Column(4, [1, 5, 15, 35, 80])
        d = Column(4, [1.0, 5.0, 15.0, 35.0, 70.0])

        @test isfirst(a)
        @test !isfirst(b)

        @test PT.isvalid(a)
        @test !PT.isvalid(c)

        @test toarray(b) == toarray(d)
    end

    @testset "Movement" begin
        a = Column(0, 4)
        b = Column(8, 4)
        c = Column(9, 4)
        d = Column(10, 4)

        @test_throws OutOfBoundsError left!(a)
        @test_throws OutOfBoundsError left(a)
        @test_throws OutOfBoundsError prev!(a)
        @test_throws OutOfBoundsError prev(a)

        @test left(c) == prev(c) == b
        @test right(c) == next(c) == d

        e = Column(b)
        f = Column(d)

        right!(b)
        next!(e)
        left!(d)
        prev!(f)
        @test b == c == d == e == f
    end
end

@testset "LazyColumn" begin
    @testset "Constructors" begin
        @test_throws DomainError LazyColumn(-1, Dict{Int,Int}())

        a = LazyColumn(0, Dict(1 => 1))
        b = LazyColumn(a)
        c = LazyColumn{Integer}(0)
        d = LazyColumn(0)
        @test a == b == c == d

        b.data[2] = 4
        @test a != b

        a = LazyColumn(3, Dict(1 => 1.0))
        b = LazyColumn(a)
        c = LazyColumn{Float64}(3)
        d = LazyColumn(3) # different data type
        @test a == b == c == d
    end

    @testset "Accessors" begin
        a = LazyColumn(5)
        @test colnumber(a) == 5
    end

    @testset "Array functions" begin
        a = LazyColumn(7)
        b = LazyColumn{Float64}(7)

        @testset "Data access $i" for i ∈ 7:14
            @test a[i] == b[i]
        end

        @test_throws BoundsError a[-1]

        @test firstindex(a) == firstindex(b) == 7
    end

    @testset "Checks" begin
        a = LazyColumn(0)
        b = LazyColumn(4,
            Dict(1 => 1, 2 => 5, 3 => 15, 4 => 35, 5 => 70))
        c = LazyColumn(4,
            Dict(1 => 1, 2 => 5, 3 => 15, 4 => 35, 5 => 80))
        d = LazyColumn(4,
            Dict(1 => 1.0, 2 => 5.0, 3 => 15.0, 4 => 35.0, 5 => 70.0))

        @test isfirst(a)
        @test !isfirst(b)

        @test PT.isvalid(a)
        @test !PT.isvalid(c)

        @test toarray(b) == toarray(d)
    end

    @testset "Movement" begin
        a = LazyColumn(0)
        b = LazyColumn(8)
        c = LazyColumn(9)
        d = LazyColumn(10)

        @test_throws OutOfBoundsError left!(a)
        @test_throws OutOfBoundsError left(a)
        @test_throws OutOfBoundsError prev!(a)
        @test_throws OutOfBoundsError prev(a)

        @test left(c) == prev(c) == b
        @test right(c) == next(c) == d

        e = LazyColumn(b)
        f = LazyColumn(d)

        right!(b)
        next!(e)
        left!(d)
        prev!(f)
        @test b == c == d == e == f
    end
end


text:
module TextMining

using ASCIIPlots, LightXML

import Base: getindex, setindex!, isempty, keys, values,
             copy, length, haskey, display, show

abstract FeatureSpace

include("feature_vector.jl")
include("cluster.jl")
include("data_set.jl")
include("clustering.jl")
include("distribution.jl")
include("text_processing.jl")
include("naive_bayes.jl")
include("tree.jl")

export
  #types
  Cluster, DataSet, Distribution, FeatureSpace, FeatureVector,
  BinaryTree, BinaryTreeNode, EmptyTree,

  #functions

  #naive bayes
  get_probs_fv_in_clust, split_dataset, separate_by_class, naive_bayes,
  train_data, percentages,

  #classification
  knn,

  #cluster
  centroid, distance, dist_centroid, dist_matrix,

  #clustering
  random_init, max_min_init, kmeans, elbow_method, hclust,

  #data set

  #distribution
  get_num_texts_in_dist, get_num_texts_given_feature, probability, prob_fv,
  prob_cl, cond_prob_f_given_fv, cond_prob_f_given_fv, cond_prob_f_given_fv_in_clust,
  cond_prob_f_given_clust, cond_prob_fv_given_clust, prob_clust_in_dataset,
  prob_clust_given_feature, features, entropy, fv_entropy, clust_entropy,
  feature_entropy, clust_in_dataset_entropy, clust_given_feature_entropy, clust_info_gain,
  ld_info_gain, feature_info_gain, info_gain, perplexity, remove_smoothing!,
  delta_smoothing!, goodturing_smoothing!,

  #feature vector
  sanitize!, freq_list, find_common_type, add!, subtract!,
  multiply!, divide!, rationalize!, dist_cos, dist_zero,
  dist_zero_weighted, dist_taxicab, dist_euclidean, dist_infinite,

  #text processing
  clean, parse_xml, load_featurevector, load_cluster, load_dataset,
  get_metadata

end


text:
using LoopVectorization, Polyester

mutable struct Layer
    input_size::Int64
    layer_size::Int64

    weights::Array{Float32}
    voltage::Array{Float32}
    output::Array{Float32}

    fire_rec::Array{Int64}

    Pᵢ::Array{Float32}
    M::Array{Float32}

    guide::Bool
    record::Bool
    sample_output::Array{Float32}

    activate
    reset

    function Layer(;input_size::Int64, layer_size::Int64)
        weights = rand(0.0f0:1.0f-5:4.0f-1, layer_size,input_size)
        voltage = zeros(Float32, layer_size)
        output = zeros(Float32, layer_size)

        fire_rec = zeros(Int64, layer_size)

        Pᵢ = zeros(Float32, layer_size,input_size)
        M = zeros(Float32, layer_size)

        new(input_size, layer_size, weights, voltage, output, fire_rec, Pᵢ, M, false, false, Float32[], activate, reset)
    end

    function activate(self::Layer, input::Array{Float32}, η, A₊, A₋, γᵥ, γₚ, γₘ, w₊, w₋, θ, t, r, first_spike, table)
        @batch for l in 1:self.layer_size
            if t-self.fire_rec[l]>=r
                s = 0
                @avx for i in 1:self.input_size
                    s += self.weights[l,i] * input[i]
                    self.Pᵢ[l,i] *= γₚ
                    self.Pᵢ[l,i] += input[i] * A₊
                    self.weights[l,i] += η * self.M[l] * (w₋-self.weights[l,i])
                end
                self.voltage[l] *= γᵥ
                self.voltage[l] += s
            end
        end

        if self.guide || self.record
            if !first_spike
                if maximum(self.voltage)>θ
                    if self.record
                        first_spike = true
                        winner = argmax(self.voltage)
                        @avxt for i in 1:self.layer_size
                            table[winner, i] += self.sample_output[i]
                        end
                    end
                    if self.guide
                        @avxt for l in 1:self.layer_size
                            self.voltage[l] *= self.sample_output[l]
                        end
                    end
                end
            else
                if self.guide
                    @avxt for l in 1:self.layer_size
                        self.voltage[l] *= self.sample_output[l]
                    end
                end
            end
            if self.guide
                if self.voltage[argmax(self.sample_output)]<θ && t!=0 && t%r==0
                    self.voltage[argmax(self.sample_output)] += θ
                end
            end
        end

        @batch for l in 1:self.layer_size
            self.M[l] *= γₘ
            if self.voltage[l]>=θ
                self.voltage[l] = 0
                self.output[l] = 1
                self.fire_rec[l] = t
                @avx for i in 1:self.input_size
                    self.weights[l,i] += η * self.Pᵢ[l,i] * (w₊-self.weights[l,i])
                    self.Pᵢ[l,i] = 0
                end
                self.M[l] -= A₋
            elseif self.output[l]!=0
                self.output[l] = 0
            end
        end

        return first_spike
    end

    function reset(self::Layer, r)
        @avxt self.voltage .= 0.0f0
        @avxt self.Pᵢ .= 0.0f0
        @avxt self.M .= 0.0f0
        @avxt self.fire_rec .= -r
    end
end


text:
### Day 02 - Problem 2
## Bathroom Security
## Author: Thanasis Georgiou

# Load puzzle input
f = open("bathroom-code.txt")
input = readlines(f)
close(f)

# Pad buttons
pad = [
    ' ' ' ' '1' ' ' ' ';
    ' ' '2' '3' '4' ' ';
    '5' '6' '7' '8' '9';
    ' ' 'A' 'B' 'C' ' ';
    ' ' ' ' 'D' ' ' ' '
]

# Finger location
type Point
    x::Int
    y::Int
end

finger = Point(2,2)

# Bathroom PIN
PIN = ""

# Parse instructions
for line in input
    for char in line
        # Create test finger
        testFinger = Point(finger.x, finger.y)

        # Move test finger
        if char == 'U'
            testFinger.x -= 1
        elseif char == 'D'
            testFinger.x += 1
        elseif char == 'L'
            testFinger.y -= 1
        elseif char == 'R'
            testFinger.y += 1
        end

        # Clamp within array bounds
        testFinger.x = clamp(testFinger.x, 1, 5)
        testFinger.y = clamp(testFinger.y, 1, 5)

        # Check if the finger landed on a button
        if pad[testFinger.x, testFinger.y] != ' '
            # Move the actual finger
            finger.x = testFinger.x
            finger.y = testFinger.y
        end
    end

    # Add key under finger to passcode
    PIN = "$(PIN)$(string(pad[finger.x, finger.y]))"
end

# Output result
println("PIN: $PIN")


text:
function NNPredict(θ,LayerStruct)
 k=size(LayerStruct,1)-1           #Number of additional layers ( It does not include the first layer)

    ΘNparams=ΘNparams=Array{Int64,1}(undef,k)

    for i=1:k
       ΘNparams[i]=(LayerStruct[i]+1)*LayerStruct[i+1]
    end
    p= function (X)
        m=size(X,1)
                    Θ=Array{Array{Float64,2},1}(undef, k)   #Initialize the shape of the Θ Array of Arrays
                    a=Array{Array{Float64,2},1}(undef, k+1)
                    z=Array{Array{Float64,2},1}(undef, k)
                    a[1]=[ones(m) X]
                    r=zeros(k)
                    c=0
                    for i=1:k
                        Θ[i]=reshape(θ[c+1:c+ΘNparams[i]],(LayerStruct[i+1],LayerStruct[i]+1))
                        c=ΘNparams[i]
                        z[i]=a[i]*Θ[i]'
                        a[i+1]=Sigmoid(z[i])
                        a[i+1]=[ones(m) a[i+1]]
                    end

        s=zeros(m)
        for i=1:m
            s[i]=argmax(a[end][i,2:end])
        end
        return s
    end
        return p
end


text:
struct Geometric <: Distribution{Int} end

"""
    geometric(p::Real)

Sample an `Int` from the Geometric distribution with parameter `p`.
"""
const geometric = Geometric()

function Gen.logpdf(::Geometric, x::Int, p::Real)
    Distributions.logpdf(Distributions.Geometric(p), x)
end

function Gen.logpdf_grad(::Geometric, x::Int, p::Real)
    if x >= 0
        p_grad = 1/p - (1/(1-p) * x)
    else
        p_grad = 0.0
    end
    (nothing, p_grad)
end

function Gen.random(::Geometric, p::Real)
    rand(Distributions.Geometric(p))
end

is_discrete(::Geometric) = true

(::Geometric)(p) = random(Geometric(), p)

Gen.has_output_grad(::Geometric) = false
Gen.has_argument_grads(::Geometric) = (true,)

export geometric


text:
using Oceananigans.BoundaryConditions: ContinuousBoundaryFunction

function test_boundary_condition(arch, FT, topo, side, field_name, boundary_condition)
    grid = RegularCartesianGrid(FT, size=(1, 1, 1), extent=(1, π, 42), topology=topo)

    boundary_condition_kwarg = Dict(side => boundary_condition)
    field_boundary_conditions = TracerBoundaryConditions(grid; boundary_condition_kwarg...)
    bcs = NamedTuple{(field_name,)}((field_boundary_conditions,))

    model = IncompressibleModel(grid=grid, architecture=arch, float_type=FT,
                                boundary_conditions=bcs)

    success = try
        time_step!(model, 1e-16, euler=true)
        true
    catch err
        @warn "test_boundary_condition errored with " * sprint(showerror, err)
        false
    end

    return success
end

function test_flux_budget(arch, FT, fldname)
    N, κ, Lz = 16, 1, 0.7
    grid = RegularCartesianGrid(FT, size=(N, N, N), extent=(1, 1, Lz))

    bottom_flux = FT(0.3)
    flux_bc = BoundaryCondition(Flux, bottom_flux)

    if fldname == :u
        field_bcs = UVelocityBoundaryConditions(grid, bottom=flux_bc)
    elseif fldname == :v
        field_bcs = VVelocityBoundaryConditions(grid, bottom=flux_bc)
    else
        field_bcs = TracerBoundaryConditions(grid, bottom=flux_bc)
    end

    model_bcs = NamedTuple{(fldname,)}((field_bcs,))

    closure = IsotropicDiffusivity(FT, ν=κ, κ=κ)
    model = IncompressibleModel(grid=grid, closure=closure, architecture=arch, tracers=(:T, :S),
                                float_type=FT, buoyancy=nothing, boundary_conditions=model_bcs)

    field = get_model_field(fldname, model)
    @. field.data = 0

    τκ = Lz^2 / κ   # Diffusion time-scale
    Δt = 1e-6 * τκ  # Time step much less than diffusion time-scale
    Nt = 10         # Number of time steps

    for n in 1:Nt
        time_step!(model, Δt, euler= n==1)
    end

    # budget: Lz*∂<ϕ>/∂t = -Δflux = -top_flux/Lz (left) + bottom_flux/Lz (right)
    # therefore <ϕ> = bottom_flux * t / Lz
    return mean(interior(field)) ≈ bottom_flux * model.clock.time / Lz
end

function fluxes_with_diffusivity_boundary_conditions_are_correct(arch, FT)
    Lz = 1
    κ₀ = FT(exp(-3))
    bz = FT(π)
    flux = - κ₀ * bz

    grid = RegularCartesianGrid(FT, size=(16, 16, 16), extent=(1, 1, Lz))

    buoyancy_bcs = TracerBoundaryConditions(grid, bottom=BoundaryCondition(Gradient, bz))
    κₑ_bcs = DiffusivityBoundaryConditions(grid, bottom=BoundaryCondition(Value, κ₀))
    model_bcs = (b=buoyancy_bcs, κₑ=(b=κₑ_bcs,))

    model = IncompressibleModel(
        grid=grid, architecture=arch, float_type=FT, tracers=:b, buoyancy=BuoyancyTracer(),
        closure=AnisotropicMinimumDissipation(), boundary_conditions=model_bcs
    )

    b₀(x, y, z) = z * bz
    set!(model, b=b₀)

    b = model.tracers.b
    mean_b₀ = mean(interior(b))

    τκ = Lz^2 / κ₀  # Diffusion time-scale
    Δt = 1e-6 * τκ  # Time step much less than diffusion time-scale
    Nt = 10         # Number of time steps

    for n in 1:Nt
        time_step!(model, Δt, euler= n==1)
    end

    # budget: Lz*∂<ϕ>/∂t = -Δflux = -top_flux/Lz (left) + bottom_flux/Lz (right)
    # therefore <ϕ> = bottom_flux * t / Lz
    #
    # Use an atol of 1e-6 so test passes with Float32 as there's a big cancellation
    # error due to buoyancy order of magnitude.
    #
    # Float32:
    # mean_b₀ = -1.5707965f0
    # mean(interior(b)) = -1.5708286f0
    # mean(interior(b)) - mean_b₀ = -3.20673f-5
    # (flux * model.clock.time) / Lz = -3.141593f-5
    #
    # Float64
    # mean_b₀ = -1.5707963267949192
    # mean(interior(b)) = -1.57082774272148
    # mean(interior(b)) - mean_b₀ = -3.141592656086267e-5
    # (flux * model.clock.time) / Lz = -3.141592653589793e-5

    return isapprox(mean(interior(b)) - mean_b₀, flux * model.clock.time / Lz, atol=1e-6)
end

test_boundary_conditions(C, FT, ArrayType) = (integer_bc(C, FT, ArrayType),
                                              float_bc(C, FT, ArrayType),
                                              irrational_bc(C, FT, ArrayType),
                                              array_bc(C, FT, ArrayType),
                                              simple_function_bc(C, FT, ArrayType),
                                              parameterized_function_bc(C, FT, ArrayType),
                                              field_dependent_function_bc(C, FT, ArrayType),
                                              parameterized_field_dependent_function_bc(C, FT, ArrayType),
                                              discrete_function_bc(C, FT, ArrayType),
                                              parameterized_discrete_function_bc(C, FT, ArrayType)
                                             )

@testset "Boundary condition integration tests" begin
    @info "Testing boundary condition integration into IncompressibleModel..."

    @testset "Boundary condition regularization" begin
        @info "  Testing boundary condition regularization in IncompressibleModel constructor..."

        FT = Float64
        arch = first(archs)

        # We use Periodic, Bounded, Bounded here because triply Bounded domains don't work on the GPU
        # yet.
        grid = RegularCartesianGrid(FT, size=(1, 1, 1), extent=(1, π, 42), topology=(Periodic, Bounded, Bounded))

        u_boundary_conditions = UVelocityBoundaryConditions(grid;
                                                            bottom = simple_function_bc(Value),
                                                            top    = simple_function_bc(Value),
                                                            north  = simple_function_bc(Value),
                                                            south  = simple_function_bc(Value))

        v_boundary_conditions = VVelocityBoundaryConditions(grid;
                                                            bottom = simple_function_bc(Value),
                                                            top    = simple_function_bc(Value),
                                                            north  = simple_function_bc(NormalFlow),
                                                            south  = simple_function_bc(NormalFlow))

        w_boundary_conditions = VVelocityBoundaryConditions(grid;
                                                            bottom = simple_function_bc(NormalFlow),
                                                            top    = simple_function_bc(NormalFlow),
                                                            north  = simple_function_bc(Value),
                                                            south  = simple_function_bc(Value))

        T_boundary_conditions = TracerBoundaryConditions(grid;
                                                         bottom = simple_function_bc(Value),
                                                         top    = simple_function_bc(Value),
                                                         north  = simple_function_bc(Value),
                                                         south  = simple_function_bc(Value))

        boundary_conditions = (u=u_boundary_conditions,
                               v=v_boundary_conditions,
                               w=w_boundary_conditions,
                               T=T_boundary_conditions)

        model = IncompressibleModel(architecture = arch,
                                    grid = grid,
                                    float_type = FT,
                                    boundary_conditions = boundary_conditions)

        @test location(model.velocities.u.boundary_conditions.bottom.condition) == (Face, Cell, Nothing)
        @test location(model.velocities.u.boundary_conditions.top.condition)    == (Face, Cell, Nothing)
        @test location(model.velocities.u.boundary_conditions.north.condition)  == (Face, Nothing, Cell)
        @test location(model.velocities.u.boundary_conditions.south.condition)  == (Face, Nothing, Cell)

        @test location(model.velocities.v.boundary_conditions.bottom.condition) == (Cell, Face, Nothing)
        @test location(model.velocities.v.boundary_conditions.top.condition)    == (Cell, Face, Nothing)
        @test location(model.velocities.v.boundary_conditions.north.condition)  == (Cell, Nothing, Cell)
        @test location(model.velocities.v.boundary_conditions.south.condition)  == (Cell, Nothing, Cell)

        @test location(model.velocities.w.boundary_conditions.bottom.condition) == (Cell, Cell, Nothing)
        @test location(model.velocities.w.boundary_conditions.top.condition)    == (Cell, Cell, Nothing)
        @test location(model.velocities.w.boundary_conditions.north.condition)  == (Cell, Nothing, Face)
        @test location(model.velocities.w.boundary_conditions.south.condition)  == (Cell, Nothing, Face)

        @test location(model.tracers.T.boundary_conditions.bottom.condition) == (Cell, Cell, Nothing)
        @test location(model.tracers.T.boundary_conditions.top.condition)    == (Cell, Cell, Nothing)
        @test location(model.tracers.T.boundary_conditions.north.condition)  == (Cell, Nothing, Cell)
        @test location(model.tracers.T.boundary_conditions.south.condition)  == (Cell, Nothing, Cell)
    end

    @testset "Boudnary condition time-stepping works" begin
        for arch in archs, FT in (Float64,) #float_types
            @info "  Testing that time-stepping with boundary conditions works [$(typeof(arch)), $FT]..."

            topo = arch isa CPU ? (Bounded, Bounded, Bounded) : (Periodic, Bounded, Bounded)

            for C in (Gradient, Flux, Value), boundary_condition in test_boundary_conditions(C, FT, array_type(arch))
                arch isa CPU && @test test_boundary_condition(arch, FT, topo, :east, :T, boundary_condition)

                @test test_boundary_condition(arch, FT, topo, :south, :T, boundary_condition)
                @test test_boundary_condition(arch, FT, topo, :top, :T, boundary_condition)
            end

            for boundary_condition in test_boundary_conditions(NormalFlow, FT, array_type(arch))
                 arch isa CPU && @test test_boundary_condition(arch, FT, topo, :east, :u, boundary_condition)

                @test test_boundary_condition(arch, FT, topo, :south, :v, boundary_condition)
                @test test_boundary_condition(arch, FT, topo, :top, :w, boundary_condition)
            end
        end
    end

    @testset "Budgets with Flux boundary conditions" begin
        for arch in archs, FT in float_types
            @info "  Testing budgets with Flux boundary conditions on u, v, T [$(typeof(arch)), $FT]..."
            for field_name in (:u, :v, :T)
                @test test_flux_budget(arch, FT, field_name)
            end
        end
    end

    @testset "Custom diffusivity boundary conditions" begin
        for arch in archs, FT in (Float64,) #float_types
            @info "  Testing flux budgets with diffusivity boundary conditions [$(typeof(arch)), $FT]..."
            @test fluxes_with_diffusivity_boundary_conditions_are_correct(arch, FT)
        end
    end
end


text:
"""
    KNN(; kwargs...)

Imputation using k-Nearest Neighbor algorithm.

# Keyword Arguments
* `k::Int`: number of nearest neighbors
* `dist::MinkowskiMetric`: distance metric suppports by `NearestNeighbors.jl` (Euclidean, Chebyshev, Minkowski and Cityblock)
* `threshold::AbstractFloat`: threshold for missing neighbors

# Reference
* Troyanskaya, Olga, et al. "Missing value estimation methods for DNA microarrays." Bioinformatics 17.6 (2001): 520-525.
"""
struct KNN{M} <: Imputor where M <: NearestNeighbors.MinkowskiMetric
    k::Int
    threshold::AbstractFloat
    dist::M
end

# TODO : Support Categorical Distance (NearestNeighbors.jl support needed)
function KNN(; k=1, threshold=0.5, dist=Euclidean())
    k < 1 && throw(ArgumentError("The number of nearset neighbors should be greater than 0"))

    !(0 < threshold < 1) && throw(ArgumentError("Missing neighbors threshold should be within 0 to 1"))

    # to exclude missing value itself
    KNN(k + 1, threshold, dist)
end

function impute!(data::AbstractMatrix{Union{T, Missing}}, imp::KNN; dims=nothing) where T<:Real
    d = dim(data, dims)

    # KDTree expects data of the form dims x n
    X = d == 1 ? data : transpose(data)

    # Get mask array first
    missing_mask = ismissing.(X)

    # Fill missing value as mean value
    impute!(X, Substitute(); dims=1)

    # Disallow `missings` for NearestNeighbors
    X = disallowmissing(X)

    # Our search points are just observations containing `missing`s
    points = X[:, vec(any(missing_mask; dims=1))]

    # Contruct our KDTree over the entire dataset
    kdtree = KDTree(X, imp.dist)

    # Query for neighbors to our missing observations
    # NOTES:
    # 1. It's generally faster to query for all points at once
    # 2. We wanted the results sorted so that the first idx is our data points
    #   location in the original dataset.
    for (idxs, dists) in zip(NearestNeighbors.knn(kdtree, points, imp.k, true)...)
        # Our closest neighbor should always be our input data point (distance of zero)
        @assert iszero(first(dists))

        # Location of point to impute
        j = first(idxs)

        # Update each missing value in this point
        for i in axes(points, 1)
            # Skip non-missing elements
            missing_mask[i, j] || continue

            # Grab our neighbor mask to excluding neighbor values that were also missing.
            neighbor_mask = missing_mask[i, idxs]

            # Skip if there are too many missing neighbor values
            (count(neighbor_mask) / imp.k) > imp.threshold && continue

            # Weight valid neighbors based on inverse distance
            neighbor_dists = dists[.!neighbor_mask]
            wv = weights(1.0 ./ neighbor_dists)

            # Only fill with the weighted mean of neighbors if the sum of the weights are
            # non-zero and finite.
            if isfinite(sum(wv)) && !iszero(sum(wv))
                neighbor_vals = X[i, idxs[.!neighbor_mask]]
                X[i, j] = mean(neighbor_vals, wv)
            end
        end
    end

    # for type stability
    return allowmissing(d == 1 ? X : X')
end

impute!(data::AbstractMatrix{Missing}, imp::KNN; kwargs...) = data

function impute(data::AbstractMatrix{Union{T, Missing}}, imp::KNN; kwargs...) where T<:Real
    return impute!(trycopy(data), imp; kwargs...)
end

impute(data::AbstractMatrix{Missing}, imp::KNN; kwargs...) = trycopy(data)


text:
@inline function ExaTronKernels.daxpy(n::Int,da::Float64,
                       dx::CuDeviceArray{Float64,1},incx::Int,
                       dy::CuDeviceArray{Float64,1},incy::Int)
    tx = threadIdx().x

    if tx <= n
        @inbounds dy[tx] = dy[tx] + da*dx[tx]
    end
    CUDA.sync_threads()

    return
end


text:
using Test

include("./max-total.jl")

@testset "4 rows" begin
  input = """
  3
  7 4
  2 4 6
  8 5 9 3
  """

  @test calc_max_path(input; show_path=true) == (23, [(1, 1), (2, 1), (3, 2), (4, 3)])
end

@testset "15 rows" begin
  input = """
  75
  95 64
  17 47 82
  18 35 87 10
  20 04 82 47 65
  19 01 23 75 03 34
  88 02 77 73 07 63 67
  99 65 04 28 06 16 70 92
  41 41 26 56 83 40 80 70 33
  41 48 72 33 47 32 37 16 94 29
  53 71 44 65 25 43 91 52 97 51 14
  70 11 33 28 77 73 17 78 39 68 17 57
  91 71 52 38 17 14 91 43 58 50 27 29 48
  63 66 04 68 89 53 67 30 73 16 69 87 40 31
  04 62 98 27 23 09 70 98 73 93 38 53 60 04 23
  """
  @test calc_max_path(input; show_path=true) == (1074,
                                                 [(1, 1), (2, 2), (3, 3), (4, 3), (5, 3), (6, 4), (7, 4), (8, 4),
                                                  (9, 5), (10, 6), (11, 7), (12, 8), (13, 9), (14, 9), (15, 10)])
  # 1074 ≡ Σ [75, 64, 82, 87, 82, 75, 73, 28, 83, 32, 91, 78, 58, 73, 93]
end

@testset "100 rows" begin
  @test calc_max_path("./p067_triangle.txt")[1] == 7273
end

max_stars_repo_name: pascal-p/Euler

max_issues_repo_name: pascal-p/Euler

max_forks_repo_name: pascal-p/Euler


text:
"""
Stores results data for one DecisionModel
"""
mutable struct DecisionModelStore <: AbstractModelStore
    duals::Dict{ConstraintKey, OrderedDict{Dates.DateTime, DataFrames.DataFrame}}
    parameters::Dict{ParameterKey, OrderedDict{Dates.DateTime, DataFrames.DataFrame}}
    variables::Dict{VariableKey, OrderedDict{Dates.DateTime, DataFrames.DataFrame}}
    aux_variables::Dict{AuxVarKey, OrderedDict{Dates.DateTime, DataFrames.DataFrame}}
    expressions::Dict{ExpressionKey, OrderedDict{Dates.DateTime, DataFrames.DataFrame}}
    optimizer_stats::OrderedDict{Dates.DateTime, OptimizerStats}
end

function DecisionModelStore()
    return DecisionModelStore(
        Dict{
            ConstraintKey,
            Dict{ConstraintKey, OrderedDict{Dates.DateTime, DataFrames.DataFrame}},
        }(),
        Dict{
            ParameterKey,
            Dict{ParameterKey, OrderedDict{Dates.DateTime, DataFrames.DataFrame}},
        }(),
        Dict{
            VariableKey,
            Dict{VariableKey, OrderedDict{Dates.DateTime, DataFrames.DataFrame}},
        }(),
        Dict{
            AuxVarKey,
            Dict{ConstraintKey, OrderedDict{Dates.DateTime, DataFrames.DataFrame}},
        }(),
        Dict{
            AuxVarKey,
            Dict{ExpressionKey, OrderedDict{Dates.DateTime, DataFrames.DataFrame}},
        }(),
        OrderedDict{Dates.DateTime, OptimizerStats}(),
    )
end

function initialize_storage!(
    store::DecisionModelStore,
    container::OptimizationContainer,
    params::ModelStoreParams,
)
    num_of_executions = get_num_executions(params)
    time_steps_count = get_time_steps(container)[end]
    initial_time = get_initial_time(container)
    model_interval = get_interval(params)
    for type in STORE_CONTAINERS
        field_containers = getfield(container, type)
        results_container = getfield(store, type)
        for (key, field_container) in field_containers
            !should_write_resulting_value(key) && continue
            @debug "Adding $(encode_key_as_string(key)) to DecisionModelStore" _group =
                LOG_GROUP_MODEL_STORE
            results_container[key] = OrderedDict{Dates.DateTime, DataFrames.DataFrame}()
            column_names = get_column_names(key, field_container)
            for timestamp in
                range(initial_time, step=model_interval, length=num_of_executions)
                results_container[key][timestamp] = DataFrames.DataFrame(
                    OrderedDict(c => fill(NaN, time_steps_count) for c in column_names),
                )
            end
        end
    end
end

function write_result!(
    store::DecisionModelStore,
    name::Symbol,
    key::OptimizationContainerKey,
    index::DecisionModelIndexType,
    update_timestamp::Dates.DateTime,
    array::AbstractArray,
)
    df = axis_array_to_dataframe(array, key)
    write_result!(store, name, key, index, update_timestamp, df)
    return
end

function write_result!(
    store::DecisionModelStore,
    ::Symbol,
    key::OptimizationContainerKey,
    index::DecisionModelIndexType,
    update_timestamp::Dates.DateTime,
    df::Union{DataFrames.DataFrame, DataFrames.DataFrameRow},
)
    container = getfield(store, get_store_container_type(key))
    container[key][index] = df
    return
end

function read_results(
    store::DecisionModelStore,
    key::OptimizationContainerKey;
    index::Union{DecisionModelIndexType, Nothing}=nothing,
)
    container = getfield(store, get_store_container_type(key))
    data = container[key]
    if isnothing(index)
        @assert length(data) == 1
        index = first(keys(data))
    end

    # Return a copy because callers may mutate it.
    return copy(data[index], copycols=true)
end

function write_optimizer_stats!(
    store::DecisionModelStore,
    stats::OptimizerStats,
    index::DecisionModelIndexType,
)
    @assert !(index in keys(store.optimizer_stats))
    store.optimizer_stats[index] = stats
    return
end

function read_optimizer_stats(store::DecisionModelStore)
    stats = [to_namedtuple(x) for x in values(store.optimizer_stats)]
    df = DataFrames.DataFrame(stats)
    DataFrames.insertcols!(df, 1, :DateTime => keys(store.optimizer_stats))
    return df
end


text:
"""
    BayesianLinearRegressor{Tmw, TΛw}

A Bayesian Linear Regressor is a distribution over linear functions given by
```julia
w ~ Normal(mw, Λw)
f(x) = dot(x, w)
```
where `mw` and `Λw` are the mean and precision of `w`, respectively.
"""
struct BayesianLinearRegressor{Tmw<:AbstractVector,TΛw<:AbstractMatrix} <: AbstractGP
    mw::Tmw
    Λw::TΛw
end

const FiniteBLR = FiniteGP{<:BayesianLinearRegressor}

# All code below implements the primary + secondary AbstractGPs.jl APIs.

x_as_colvecs(fx::FiniteBLR) = x_as_colvecs(fx.x)

x_as_colvecs(x::ColVecs) = x

x_as_colvecs(x::RowVecs) = ColVecs(x.X')

function x_as_colvecs(x::T) where {T<:AbstractVector}
    return error(
        "$T is not a subtype of AbstractVector that is known. Please provide either a",
        "ColVecs or RowVecs.",
    )
end

AbstractGPs.mean(fx::FiniteBLR) = x_as_colvecs(fx).X' * fx.f.mw

function AbstractGPs.cov(fx::FiniteBLR)
    α = _cholesky(fx.f.Λw).U' \ x_as_colvecs(fx).X
    return Symmetric(α' * α + fx.Σy)
end

function AbstractGPs.var(fx::FiniteBLR)
    α = _cholesky(fx.f.Λw).U' \ x_as_colvecs(fx).X
    return vec(sum(abs2, α; dims=1)) .+ diag(fx.Σy)
end

AbstractGPs.mean_and_cov(fx::FiniteBLR) = (mean(fx), cov(fx))

AbstractGPs.mean_and_var(fx::FiniteBLR) = (mean(fx), var(fx))

function AbstractGPs.rand(rng::AbstractRNG, fx::FiniteBLR, samples::Int)
    X = x_as_colvecs(fx).X
    w = fx.f.mw .+ _cholesky(fx.f.Λw).U \ randn(rng, size(X, 1), samples)
    return X' * w .+ _cholesky(fx.Σy).U' * randn(rng, size(X, 2), samples)
end

function AbstractGPs.logpdf(fx::FiniteBLR, y::AbstractVector{<:Real})
    _, Bt, δy, logpdf_δy, Λεy = __compute_inference_quantities(fx, y)
    return -(logdet(Λεy) - sum(abs2, Λεy.U' \ (Bt'δy))) / 2 + logpdf_δy
end

function AbstractGPs.posterior(fx::FiniteBLR, y::AbstractVector{<:Real})
    Uw, Bt, δy, _, Λεy = __compute_inference_quantities(fx, y)

    # Compute posterior over decorrelated weights.
    mεy = Λεy \ (Bt'δy)

    # Compute posterior over weights.
    T = Λεy.U * Uw
    return BayesianLinearRegressor(fx.f.mw + Uw \ mεy, Symmetric(T'T))
end

# Computation utilised in both `logpdf` and `posterior`.
function __compute_inference_quantities(fx::FiniteBLR, y::AbstractVector{<:Real})
    X = x_as_colvecs(fx).X
    length(y) == size(X, 2) || throw(error("length(y) != size(fx.x.X, 2)"))
    blr = fx.f
    N = length(y)

    Uw = _cholesky(blr.Λw).U
    Σy = _cholesky(fx.Σy)

    Bt = Σy.U' \ (Uw' \ X)'
    δy = Σy.U' \ (y - mean(fx))

    logpdf_δy = -(N * log(2π) + logdet(Σy) + sum(abs2, δy)) / 2

    Λεy = _cholesky(Symmetric(Bt'Bt + I))

    return Uw, Bt, δy, logpdf_δy, Λεy
end

# Random function sample generation
# Following the Random API: https://docs.julialang.org/en/v1/stdlib/Random/#Hooking-into-the-Random-API
struct BLRFunctionSample{Tw<:AbstractVector}
    w::Tw
end

(s::BLRFunctionSample)(X::AbstractMatrix{<:Real}) = X's.w
(s::BLRFunctionSample)(X::ColVecs) = X.X's.w
(s::BLRFunctionSample)(X::RowVecs) = X.X * s.w

function Random.Sampler(
    ::Type{<:AbstractRNG}, blr::BayesianLinearRegressor, ::Random.Repetition
)
    return blr
end

function Random.rand(rng::AbstractRNG, blr::BayesianLinearRegressor)
    w = blr.mw .+ _cholesky(blr.Λw).U \ randn(rng, size(blr.mw))
    return BLRFunctionSample(w)
end

function Random.rand(rng::AbstractRNG, blr::BayesianLinearRegressor, dims::Dims)
    ws = blr.mw .+ _cholesky(blr.Λw).U \ randn(rng, (only(size(blr.mw)), prod(dims)))
    bs = [BLRFunctionSample(w) for w in eachcol(ws)]
    return reshape(bs, dims)
end

function Random.rand!(
    rng::AbstractRNG, A::AbstractArray{<:BLRFunctionSample}, blr::BayesianLinearRegressor
)
    ws = blr.mw .+ _cholesky(blr.Λw).U \ randn(rng, (only(size(blr.mw)), prod(size(A))))
    for i in LinearIndices(A)
        @inbounds A[i] = BLRFunctionSample(ws[:, i])
    end
    return A
end


text:
using BinaryBuilder

# Collection of pre-build pandoc binaries
name = "pandoc_crossref"

# Actually, this should be v"0.3.6.4".
version = v"0.3.6"
crossref_ver = "0.3.6.4"
pandoc_ver = "2.9.2.1"

url_prefix = "https://github.com/lierdakil/pandoc-crossref/releases/download/v$(crossref_ver)/pandoc-crossref"
sources = [
    ArchiveSource("$(url_prefix)-Linux-$(pandoc_ver).tar.xz", "57e8f71d46c401daf0a3247b7405b49503b836e93893730e6e5907f6cb2c0885"; unpack_target = "x86_64-linux-gnu"),
    ArchiveSource("$(url_prefix)-macOS-$(pandoc_ver).tar.xz", "b63ac830d7164279eb80041cd5f793e6b611ce1a701a86308ce5b2a3a99d33b2"; unpack_target = "x86_64-apple-darwin14"),
    FileSource("$(url_prefix)-Windows-$(pandoc_ver).7z", "3e943fcdd3f91951fe18f9900e136ccb9ea810b2582f4bd929760357eaf83ef4"; filename = "x86_64-w64-mingw32"),
    FileSource("https://raw.githubusercontent.com/lierdakil/pandoc-crossref/65858c01a76f75990e7e30bcdb571cd84a69d47c/LICENSE", "39db8f9acf036595a2566ea3fe560bc7bd65d8749f088e0f4a4ef2f8a6cb4b34"),
]

# Bash recipe for building across all platforms
script = raw"""
cd ${WORKSPACE}/srcdir/
mkdir -p "${bindir}"
subdir="${target}"
if [[ "${target}" == *-mingw* ]]; then
    7z x "${target}"
    subdir="windows-build"
fi
install -m 755 "${subdir}/pandoc-crossref${exeext}" "${bindir}/pandoc-crossref${exeext}"
install_license LICENSE
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = [
    Platform("x86_64", "linux"),
    Platform("x86_64", "macos"),
    Platform("x86_64", "windows"),
]

# The products that we will ensure are always built
products = [
    ExecutableProduct("pandoc-crossref", :pandoc_crossref),
]

# Dependencies that must be installed before this package can be built
dependencies = [
    HostBuildDependency("p7zip_jll"),
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies)

max_stars_repo_name: samtkaplan/Yggdrasil

max_issues_repo_name: samtkaplan/Yggdrasil

max_forks_repo_name: samtkaplan/Yggdrasil


text:
Base.@kwdef struct Spec{T} <: AbstractElement
    analyses::Tuple=()
    pkeys::NamedTuple=NamedTuple()
    style::Style=Style()
    options::NamedTuple=NamedTuple()
end

Spec(ctx::AbstractContext) = Spec{Any}(style=Style(ctx))
Spec(style::Style) = Spec{Any}(style=style)
Spec(s::Spec) = s

spec(args...; kwargs...) = Spec{Any}(options=namedtuple(args...; kwargs...))
spec(T::Union{Type, Symbol}, args...; kwargs...) = Spec{T}(options=namedtuple(args...; kwargs...))

plottype(::Spec{T}) where {T} = T

function Base.merge(t1::Spec{T1}, t2::Spec{T2}) where {T1, T2}
    T = T2 === Any ? T1 : T2
    analyses = (t1.analyses..., t2.analyses...)
    pkeys = merge(t1.pkeys, t2.pkeys)
    style = merge(t1.style, t2.style)
    options = merge(t1.options, t2.options)
    return Spec{T}(analyses, pkeys, style, options)
end

Base.:*(a1::AbstractElement, a2::AbstractElement) = merge(Spec(a1), Spec(a2))

layers(g::AbstractElement) = AlgebraicList([Spec(g)])
layers(g::AlgebraicList) = g

Base.:*(s1::AbstractElement, s2::AlgebraicList) = layers(s1) * s2
Base.:*(s1::AlgebraicList, s2::AbstractElement) = s1 * layers(s2)

const ElementOrList = Union{AbstractElement, AlgebraicList}

Base.:+(s1::ElementOrList, s2::ElementOrList) = layers(s1) + layers(s2)

# pipeline

# Expand pairs and run the analyses
function expand(sp::Spec{T}) where {T}
    analyses, pkeys, style, options = sp.analyses, sp.pkeys, sp.style, sp.options
    @assert isempty(pkeys)
    v = [Spec{T}(style=val, pkeys=key, options=options) for (key, val) in pairs(style)]
    list = AlgebraicList(v)
    return foldl((ls, an) -> apply(an, ls), analyses, init=list)
end

global_options(f, d::AlgebraicList) = NamedTuple()

# default fallback to apply a callable to a vector of key => value pairs
# if customized, it must return a vector of key => value pairs
function apply(f, d::AlgebraicList)
    global_kwargs = global_options(f, d)
    v = map(parent(d)) do layer
        analyses, pkeys, style, options = layer.analyses, layer.pkeys, layer.style, layer.options
        T = plottype(layer)
        args, kwargs = split(style.value)
        res = f(args...; global_kwargs..., kwargs...) * Spec{T}(analyses=analyses, options=options, pkeys=pkeys)
        return parent(layers(res))
    end
    return AlgebraicList(reduce(vcat, v))
end

# Expand styles, apply analyses, compute scales, and return vector of traces
function run_pipeline(s::ElementOrList)
    nested = [parent(expand(layer)) for layer in layers(s)]
    computescales(reduce(vcat, nested))
end


text:
if occursin(r"RvSpectMLEcoSystem$", pwd())   cd("EXPRES-StellarSignals")   end
using Pkg
Pkg.activate(".")

verbose = true
 make_plots = true
 if verbose && !isdefined(Main,:RvSpectML)  println("# Loading RvSpecML")    end
 using RvSpectMLBase
 using RvSpectML
 using EchelleInstruments, EchelleInstruments.EXPRES
if verbose   println("# Loading other packages")    end
 using CSV, DataFrames, Query, StatsBase, Statistics, Dates

 target_subdir = "26965"   # USER: Replace with directory of your choice
  fits_target_str = "26965"
  output_dir = "examples/output/"
  paths_to_search_for_param = [pwd(),"examples",joinpath(pkgdir(RvSpectMLBase),"..","RvSpectML","examples"), "/gpfs/group/ebf11/default/ebf11/expres"]
  # NOTE: make_manifest does not update its paths_to_search when default_paths_to_search is defined here, so if you change the line above, you must also include "paths_to_search=default_paths_to_search" in the make_manifest() function call below
  pipeline_plan = PipelinePlan()
  dont_make_plot!(pipeline_plan, :movie)

 reset_all_needs!(pipeline_plan)
 if need_to(pipeline_plan,:read_spectra)
    if verbose println("# Finding what data files are avaliable.")  end
    eval(read_data_paths(paths_to_search=paths_to_search_for_param))
    @assert isdefined(Main,:expres_data_path)
    df_files = make_manifest(expres_data_path, target_subdir, EXPRES )

    if verbose println("# Reading in customized parameters from param.jl.")  end
    eval(code_to_include_param_jl(paths_to_search=paths_to_search_for_param))

    if verbose println("# Reading in ", size(df_files_use,1), " FITS files.")  end
    @time all_spectra = map(row->EXPRES.read_data(row,store_min_data=true, store_tellurics=true, normalization=:continuum),eachrow(df_files_use))
    GC.gc()
    dont_need_to!(pipeline_plan,:read_spectra)
 end

#max_orders = 12:83
lfc_orders = 43:72
  #order_list_timeseries = extract_orders(all_spectra,pipeline_plan, orders_to_use=max_orders, recalc=true )
  order_list_timeseries = extract_orders(all_spectra,pipeline_plan, orders_to_use=lfc_orders, recalc=true )

line_list_filename = "G9.espresso.mas"   # or K2?
 line_width_50_default = 10.0e3
 lsf_width = 2.0e3
 linelist_for_ccf_fn_w_path = joinpath(pkgdir(EchelleCCFs),"data","masks",line_list_filename)
 line_list_espresso = prepare_line_list(linelist_for_ccf_fn_w_path, all_spectra, pipeline_plan, v_center_to_avoid_tellurics=ccf_mid_velocity,
      Δv_to_avoid_tellurics = 2*RvSpectMLBase.max_bc+3*lsf_width, recalc=true, verbose=true)

(ccfs, v_grid) = ccf_total(order_list_timeseries, line_list_espresso, pipeline_plan,  mask_scale_factor=4.0, ccf_mid_velocity=ccf_mid_velocity, recalc=true)

using Plots
ccfs_norm = EchelleCCFs.calc_normalized_ccfs( v_grid, ccfs, v_mid=ccf_mid_velocity, dv_min=RvSpectMLBase.max_bc,dv_max= 2*RvSpectMLBase.max_bc)
ccf_template = EchelleCCFs.calc_ccf_template(ccfs_norm, assume_normalized=true)
plot(v_grid,ccfs, labels=:none)
plot(v_grid,ccfs_norm , labels=:none)
plot(v_grid,(ccfs_norm.-ccf_template)[:,range(1,size(ccfs,2),step=10)], labels=:none)
heatmap(v_grid,1:size(ccfs_norm,2),ccfs_norm'.-ccf_template')
xlabel!("v (m/s)")
ylabel!("Observation ID")
title!("26965: CCF - <CCF>, ESPRESSO mask + SNR weights")
#savefig("ccf_heatmaps_26965.png")
line_width_50 = RvSpectMLBase.calc_line_width(v_grid,view(ccfs,:,1),frac_depth=0.5)
line_width_05 = RvSpectMLBase.calc_line_width(v_grid,view(ccfs,:,1),frac_depth=0.05)

line_list_espresso_old = copy(line_list_espresso)
 line_list_espresso = prepare_line_list(linelist_for_ccf_fn_w_path, all_spectra, pipeline_plan, v_center_to_avoid_tellurics=ccf_mid_velocity,
      Δv_to_avoid_tellurics = 1*RvSpectMLBase.max_bc+5*line_width_50+2*default_ccf_mask_v_width(EXPRES1D()), recalc=true, verbose=true)
 CSV.write(joinpath("pennstate","26965","26965_line_list_espresso.csv"), line_list_espresso)

msf = 4.0; fwtf = 2.0
 println("# mask_scale_factor = ", msf, ", frac_width_to_fit =  ", fwtf, ". ")
 ((ccfs_espresso, ccf_vars_espresso), v_grid) = ccf_total(order_list_timeseries, line_list_espresso, pipeline_plan,  mask_scale_factor=msf, range_no_mask_change=line_width_05, ccf_mid_velocity=ccf_mid_velocity, v_step=100,
   #v_max=RvSpectMLBase.max_bc,
   v_max= 0*RvSpectMLBase.max_bc+5*line_width_50,
   calc_ccf_var=true, recalc=true)
 fwtf = 1.5; println("# mask_scale_factor = ", msf, ", frac_width_to_fit =  ", fwtf, ". ")
 alg_fit_rv =  EchelleCCFs.MeasureRvFromCCFGaussian(frac_of_width_to_fit=fwtf)
 rvs_ccf_espresso = RvSpectML.calc_rvs_from_ccf_total(ccfs_espresso, ccf_vars_espresso, pipeline_plan, v_grid=v_grid, times = order_list_timeseries.times, recalc=true, bin_nightly=true, alg_fit_rv=alg_fit_rv)
 σ_rvs_espresso = deepcopy(read_cache(pipeline_plan,:σ_rvs_ccf_total))

ccfs_norm = EchelleCCFs.calc_normalized_ccfs( v_grid, ccfs_espresso, v_mid=ccf_mid_velocity, dv_min=2.5*line_width_50,dv_max= 2*RvSpectMLBase.max_bc)
ccf_template = EchelleCCFs.calc_ccf_template(ccfs_norm, assume_normalized=true)
plot(v_grid,ccfs_espresso, labels=:none)
plot(v_grid,ccfs_norm , labels=:none)
plot(v_grid,(ccfs_norm.-ccf_template)[:,range(1,size(ccfs,2),step=10)], labels=:none)
xlims!(ccf_mid_velocity-3*line_width_50,ccf_mid_velocity+3*line_width_50)
cols_to_fit = findlast(x->x<ccf_mid_velocity-3*line_width_50,v_grid):findfirst(x->x>ccf_mid_velocity+3*line_width_50,v_grid)
heatmap(view(v_grid,cols_to_fit),1:size(ccfs_norm,2),view(ccfs_norm,cols_to_fit,:)'.-view(ccf_template,cols_to_fit,:)')
#heatmap(v_grid,1:size(ccfs_norm,2),ccfs_norm'.-ccf_template')
xlabel!("v (m/s)")
ylabel!("Observation ID")
title!("26965: CCF - <CCF>, ESPRESSO mask + SNR weights")
#xlims!(ccf_mid_velocity-RvSpectMLBase.max_bc,ccf_mid_velocity+RvSpectMLBase.max_bc)
#xlims!(ccf_mid_velocity-3*line_width_50,ccf_mid_velocity+3*line_width_50)
savefig("ccf_heatmaps_26965_espresso_msf=4.0.png")

#EchelleCCFs.write_each_ccf_fits(map(s->s.metadata,all_spectra), v_grid, ccfs_espresso, ccf_vars_espresso,   fits_hdr=EchelleCCFs.make_ccf_fits_header(all_spectra[1].metadata),
#            total_vel=rvs_ccf_espresso, sigma_total_vel=σ_rvs_espresso )

max_orders = 12:83
  order_list_timeseries = extract_orders(all_spectra,pipeline_plan, orders_to_use=max_orders, recalc=true )

linelist_for_ccf_fn_w_path = joinpath(pkgdir(EchelleCCFs),"data","masks",line_list_filename)
 line_list_espresso = prepare_line_list(linelist_for_ccf_fn_w_path, all_spectra, pipeline_plan, v_center_to_avoid_tellurics=ccf_mid_velocity,
      Δv_to_avoid_tellurics = 1*RvSpectMLBase.max_bc+3*line_width_50+2*default_ccf_mask_v_width(EXPRES1D()),
      orders_to_use=max_orders, recalc=true, verbose=true)

(order_ccfs, order_ccf_vars, v_grid_order_ccfs) = ccf_orders(order_list_timeseries, line_list_espresso, pipeline_plan, mask_scale_factor=4.0, ccf_mid_velocity=ccf_mid_velocity,
 v_step=100, v_max= 0*RvSpectMLBase.max_bc+5*line_width_50, orders_to_use=max_orders, calc_ccf_var=true,recalc=true)

obs_id = 1
heatmap(v_grid, max_orders, (order_ccfs[:,:,obs_id]./mean(order_ccfs[:,:,obs_id],dims=1))')
xlabel!("v (m/s)")
ylabel!("Order index")
title!("26965: CCF_" *string(obs_id) * ", ESPRESSO mask + SNR weights")
#xlims!(ccf_mid_velocity-RvSpectMLBase.max_bc,ccf_mid_velocity+RvSpectMLBase.max_bc)
xlims!(ccf_mid_velocity-3*line_width_50,ccf_mid_velocity+3*line_width_50)

alg_fit_rv2 =  EchelleCCFs.MeasureRvFromCCFGaussian(frac_of_width_to_fit=fwtf)
order_rvs_g = zeros(length(all_spectra), length(max_orders))
order_rv_std_g = zeros(length(all_spectra), length(max_orders))
order_rvs_t = zeros(length(all_spectra), length(max_orders))
order_rv_std_t = zeros(length(all_spectra), length(max_orders))
cols_to_fit = findlast(x->x<ccf_mid_velocity-3*line_width_50,v_grid):findfirst(x->x>ccf_mid_velocity+3*line_width_50,v_grid)
for (i,ord) in enumerate(max_orders)
   println("# Order: ", i)
   if sum(view(order_ccfs,:,i,1)) > 0
      rvs_order_ccf = RvSpectML.calc_rvs_from_ccf_total(view(order_ccfs,cols_to_fit,i,:), view(order_ccf_vars,cols_to_fit,i,:), pipeline_plan, v_grid=view(v_grid,cols_to_fit), times = order_list_timeseries.times, recalc=true, bin_nightly=true, alg_fit_rv=alg_fit_rv2)
      alg_fit_rv3 =  EchelleCCFs.MeasureRvFromCCFTemplate(v_grid=view(v_grid,cols_to_fit), frac_of_width_to_fit=fwtf, template = vec(mean(view(order_ccfs,cols_to_fit,i,:),dims=2)) )
      order_rvs_g[:,i] .= read_cache(pipeline_plan, :rvs_ccf_total )
      order_rv_std_g[:,i] .= read_cache(pipeline_plan, :σ_rvs_ccf_total)
      rvs_order_ccf = RvSpectML.calc_rvs_from_ccf_total(view(order_ccfs,cols_to_fit,i,:), view(order_ccf_vars,cols_to_fit,i,:), pipeline_plan, v_grid=view(v_grid,cols_to_fit), times = order_list_timeseries.times, recalc=true, bin_nightly=true, alg_fit_rv=alg_fit_rv3)
      order_rvs_t[:,i] .= read_cache(pipeline_plan, :rvs_ccf_total )
      order_rv_std_t[:,i] .= read_cache(pipeline_plan, :σ_rvs_ccf_total)
   end
end

EchelleCCFs.write_each_ccf_fits(map(s->s.metadata,all_spectra), 100.0 .* v_grid, ccfs_espresso, ccf_vars_espresso,   fits_hdr=EchelleCCFs.make_ccf_fits_header(all_spectra[1].metadata),
                  line_list_filename=line_list_filename, orders=161 .-max_orders,ccf_orders=order_ccfs,ccf_orders_var=order_ccf_vars,
                  total_vel=100.0 .* rvs_ccf_espresso, sigma_total_vel= 100.0 .* σ_rvs_espresso,
                  order_vels = 100.0 .* order_rvs_g, sigma_order_vels = 100.0 .* order_rv_std_g )

|
|


text:
using JetPack, Jets, Test

@testset "JopTranslation, dot product test" begin
    n1,n2 = 100,120
    A = JopTranslation(rand(n1,n2), rand(n1,n2))
    lhs, rhs = dot_product_test(A, rand(domain(A)), rand(range(A)))
    @test lhs ≈ rhs
end

@testset "JopTranslation, Linearity" begin
    n1,n2 = 100,120
    A = JopTranslation(randn(n1,n2), randn(n1,n2))
    lhs,rhs = linearity_test(A)
    @test lhs ≈ rhs
end

@testset "JopTranslation, Correctness" begin
    n1,n2 = 100,120

    m = [1:n1;] * ones(n2)'
    A = JopTranslation(ones(n1,n2), zeros(n1,n2))
    d = A*m

    @test m[2:end-2,2:end-1] ≈ d[1:end-3,2:end-1]
end

@testset "JopTranslation, Correctness half width in t" begin
    n1,n2 = 100,120

    m = [1:n1;] * ones(n2)'
    A = JopTranslation(0.5*ones(n1,n2), 0.0*ones(n1,n2))
    d = A*m

    @test m[2:end-1,2:end-1] .+ 0.5 ≈ d[2:end-1,2:end-1]
end

@testset "JopTranslation, Correctness half width in x" begin
    n1,n2 = 100,120

    m = ones(n1)*[1:n2;]'
    A = JopTranslation(0.0*ones(n1,n2), 0.5*ones(n1,n2))
    d = A*m

    @test m[2:end-1,2:end-1] .+ 0.5 ≈ d[2:end-1,2:end-1]
end