text:
@everywhere using LinearAlgebra
#TRPO Specifics:
#-------------------------------------------------------------------------------
#(1) Define Fisher-Matrix as Fisher-Vector-Product / Or calculate directly:
#-------------------------------------------------------------------------------

function Fisher_Vector_Product(all_traj::Traj_Data, θ, params::Params, v::Array{Float64,1})
    #Determine data to average from:
    n_traj = all_traj.n_traj
    shuffle_arr = shuffle(1:n_traj)
    l_traj = params.l_traj
    cg_damping = params.cg_damping
    data_size = Int(floor(n_traj*0.2)) + 1

    #Constants:
    FVP = net_to_vec(θ)*0 #Preallocation
    data_len = 0

    #Go over 25% of all collected data:
    for i_true = 1:data_size
        #shuffle the trajectory taken:
        i_traj = shuffle_arr[i_true]
        l_traj = all_traj.traj[i_traj].len
        data_len += l_traj

        for i_t = 1:l_traj
            #Take current state and build derivative
            s_t = all_traj.traj[i_traj].s_t[i_t,:] #Current state
            a_t = all_traj.traj[i_traj].a_t[i_t,:] #Current state
            Jt  = ∇logπ(θ,s_t,a_t,params)
            Jt_flat = net_to_vec(Jt)
            #Fisher Vector Product
            FVP += Jt_flat*dot(Jt_flat, v)
        end
    end

    FVP /= (data_len) #Average over collected data
    FVP += cg_damping * v      #Add the damping factor
    return FVP
end

#For small Networks: get all the Training data and collect the Fisher approximation!
function update_Fisher(all_traj::Traj_Data, θ, params::Params, Fisher_old)
    n_traj = all_traj.n_traj
    data_size = n_traj#Int(floor(n_traj/5)) + 1
    fac = 0.8 #Building the weighted averge of the FIM for increased performance!

    len = Int(sqrt(length(Fisher_old)))
    FIM = zeros(len,len)
    data_len = 0
    #Go over 10% of all collected data:
    for i_traj = 1:data_size
        l_traj = all_traj.traj[i_traj].len
        data_len += l_traj
        for i_t = 1:l_traj
            #Take current state and build derivative
            s_t = all_traj.traj[i_traj].s_t[i_t,:] #Current state
            a_t = all_traj.traj[i_traj].a_t[i_t,:]
            Jt  = ∇logπ(θ,s_t,a_t,params)
            Jt_flat = net_to_vec(Jt)

            #FIM += Jt_flat*D_KL_II*Jt_flat'
            Base.LinAlg.BLAS.syr!('U',1.0,Jt_flat,FIM)
        end
    end

    #Rebuild the whole Matrix from the Symmetric one:
    FIM = (Symmetric(FIM) / data_len) #Whole Matrix

    #Keep Running Average of FIM
    FIM = fac*FIM + (1-fac)*(Fisher_old)

    return FIM
end

#Build the logarithmic gradient
@everywhere function logπ(θ::Any, s_t::Array{Float64,1},a_t::Array{Float64,1},params::Params)
    #Only for segway:
    s_mod = s_t
    #=-----------
    s_mod = zeros(8)
    s_mod[1:2]   = s_t[1:2]
    s_mod[3]     = sin(s_t[3])
    s_mod[4]     = cos(s_t[3])
    s_mod[5:end] = s_t[4:end]
    #-----------=#

    out_t   = mlp_nlayer_policy_norm(θ,s_mod,params) #.+ 0.5#Evaluate the Network
    prob_pi = sum(logpdf(a_t,out_t,params))
    return prob_pi
end

@everywhere ∇logπ = grad(logπ,1)

#-------------------------------------------------------------------------------
#(2) Define Search direction F⁻¹g -> as result of Conjugate Gradient algorithm
#-------------------------------------------------------------------------------

#Conjugate Gradient Algorithm:
function conjugate_gradient(f_Ax::Function, b::Array{Float64,1}, all_traj::Traj_Data, θ, params::Params, cg_iters=10, residual_tol=1e-10)
    # Initialise
    d = copy(b)
    r = copy(b)
    x = copy(b)*0
    rTr = r'r
    #Do the Gradient Steps:
    for iter_conj = 1: cg_iters
        z = f_Ax(all_traj, θ, params, d) #Build Fisher_Vector_Product
        α = rTr / d'z
        x += α * d
        r -= α * z
        new_rTr = r'r
        β = new_rTr / rTr
        d = r + β * d
        rTr = new_rTr
        if rTr < residual_tol
            break
        end
    end
    return x #Is Inv(A)*g => NaturalGradient direction!
end

#Conjugate Gradient Algorithm:
function conjugate_gradient_with_Fisher(b::Array{Float64,1}, FIM, cg_iters=10, residual_tol=1e-10)
    # Initialise
    d = copy(b)
    r = copy(b)
    x = copy(b)*0
    rTr = r'r
    cg_damping = params.cg_damping

    #Do the Gradient Steps:
    for iter_conj = 1: cg_iters
        z = (FIM*d) + (cg_damping * d)   #Build Fisher_Vector_Product
        α = rTr / (d'z)
        x += α * d
        r -= α * z
        new_rTr = r'r
        β = new_rTr / rTr
        d = r + β * d
        rTr = new_rTr
        if rTr < residual_tol
            break
        end
    end
    return x #Is Inv(A)*g => NaturalGradient direction!
end

#-------------------------------------------------------------------------------
#(3) Perform line search along the new space to determine improvement
#-------------------------------------------------------------------------------
function linesearch(loss_all::Function,full_step, θa, all_traj::Traj_Data, params::Params)
    accept_ratio = .1
    max_backtracks = 10
    Loss_old = loss_all(θa,all_traj,params)
    stepfrac = 1.0

    for steps = 1: max_backtracks
        θ_new = θa - stepfrac * full_step
        Loss_new = loss_all(θ_new,all_traj,params)

        actual_improve = Loss_new - Loss_old

        #Break if u have improved:
        if actual_improve <= 0.0
            return θ_new
        end

        #Decrease the Stepsize
        stepfrac *= 0.6
    end

    return θa
end

#Calculate the loss on all trajectories / or at least a part of it..
function loss_all(θ,all_traj::Traj_Data,params::Params)
    n_traj  = all_traj.n_traj
    σ       = params.σ
    δ       = params.δ

    KL_mean  = 0
    L_theta  = 0
    data_len = 0

    #Size to be evaluating from:
    data_size = Int(floor(n_traj/2)) + 1 #here only on 10% of data...
    for i_traj = 1: data_size
        # As we regard costs => a negative Advantage is desired
        # Therefore we have a smaller loss is a negative advantage occurs
        l_traj = all_traj.traj[i_traj].len
        data_len += l_traj
        A_t = all_traj.traj[i_traj].A_t

        for i_t = 1: l_traj
            #Collect from struct
            s_t = all_traj.traj[i_traj].s_t[i_t,:]
            a_t = all_traj.traj[i_traj].a_t[i_t,:]
            #The two outputs:
            mu_old = all_traj.traj[i_traj].out_t[i_t,:]

            #Only for segway:
            s_mod = s_t
            #=-----------
            s_mod = zeros(8)
            s_mod[1:2]   = s_t[1:2]
            s_mod[3]     = sin(s_t[3])
            s_mod[4]     = cos(s_t[3])
            s_mod[5:end] = s_t[4:end]
            #-----------=#

            mu_new = mlp_nlayer_policy_norm(θ,s_mod,params)  #Evaluate the Network
            #The two probabilities
            p_old = prod(all_traj.traj[i_traj].p_a_t[i_t,:])
            p_new = prod(get_π_prob(mu_new,σ,a_t))
            #Update the losses
            KL_mean += p_new * log(p_new / p_old)#kl(mu_old,mu_new,σ)
            L_theta += (p_new / p_old) * A_t[i_t]
        end
    end
    #Average
    KL_mean /= (data_len)
    L_theta /= (data_len)

    #Check if the KL divergence is below δ
    if KL_mean < δ
        return L_theta
    end

    return Inf

end

#KL divergence between two gaussian distributions
function kl(μ1::Array{Float64,1},μ2::Array{Float64,1},σ::Array{Float64,1})
    return 0.5 * (μ1 - μ2)' * ((μ1 - μ2) .* (1 ./ (σ.*σ)))
end

#-------------------------------------------------------------------------------
# Performjm the complete TRPO update:
#-------------------------------------------------------------------------------

#Update without explicity calculating FIM
function do_TRPO_update(ĝ_a, all_traj::Traj_Data, θa, params::Params)
     # Flatten the gradient
     ĝ_v = net_to_vec(ĝ_a)

     #1. Determine the step direction with conjugate Gradient and FVP:
     stepdir = conjugate_gradient(Fisher_Vector_Product_parallel, ĝ_v, all_traj, θa, params)

     #2. Calculate the max step length β
     sTAs = stepdir' * Fisher_Vector_Product_parallel(all_traj, θa, params, stepdir)
     β = sqrt( (2*params.δ) / sTAs)
     full_step = vec_to_net(stepdir * β, θa)

     #3. Perform Line-Search along Step-direction until Constraint is fullfilled:
     θ_new = linesearch(loss_all, full_step, θa, all_traj, params)
     #θ_new = θa - full_step

     return θ_new
end

function do_TNPG_update(ĝ_a, all_traj::Traj_Data, θa, params::Params)
     # Flatten the gradient
     ĝ_v = net_to_vec(ĝ_a)

     #1. Determine the step direction with conjugate Gradient and FVP:
     stepdir = conjugate_gradient(Fisher_Vector_Product_parallel, ĝ_v, all_traj, θa, params)

     #2. Calculate the max step length β
     sTAs = stepdir' * Fisher_Vector_Product_parallel(all_traj, θa, params, stepdir)
     β = sqrt( (2*params.δ) / sTAs)
     full_step = vec_to_net(stepdir * β, θa)

     #3. Perform Line-Search along Step-direction until Constraint is fullfilled:
     θ_new = θa - full_step

     return θ_new
end

#Update with explicit calculation of FIM
function do_TRPO_update_small(ĝ, all_traj::Traj_Data, θ, params::Params, Fisher)
    ĝ_v = net_to_vec(ĝ)

    FIM = update_Fisher(all_traj, θ, params, Fisher)
    #@time FIM = update_Fisher_parallel(all_traj, θ, params, Fisher)

    stepdir = conjugate_gradient_with_Fisher(ĝ_v, FIM)

    sTAs = stepdir' * FIM * stepdir
    β = 0.9*sqrt( (2*params.δ) / sTAs)

    full_step = vec_to_net(stepdir * β, θ)

    θ_new = full_step#linesearch(loss_all, full_step, θ, all_traj, params)

    return θ_new

end

"""
Function for the parallel-update of the FIM
"""
function Fisher_Vector_Product_parallel(all_traj::Traj_Data, θ::Any, params::Params, v::Array{Float64,1})
    data_percent = 0.33 #Take only partial data for increased performance
    n_traj = all_traj.n_traj
    shuffle_arr = shuffle(1:n_traj) #shuffle trajectories taken

    l_traj = params.l_traj
    cg_damping = params.cg_damping
    data_size = Int(floor(n_traj*data_percent)) + 1

    FVP = net_to_vec(θ)*0 #Preallocation
    data_len = 0

    # Create Channels storing results and starting valus
    traj_channel = RemoteChannel(()->Channel{Trajectory}(length(workers())+1)) # will contain s_init
    FVP_channel  = RemoteChannel(()->Channel{Array{Float64,1}}(length(workers())+1)) # will contain traj and gradients

    # Create problem on every worker
    begin
        for p in workers()
            remote_do(distributed_fisher_VP, p,  θ, params, traj_channel, FVP_channel, v)
        end
    end

    # Start optimizations
    for i_true = 1:data_size
        i_traj = shuffle_arr[i_true]
        data_len += all_traj.traj[i_traj].len
        loc_traj = all_traj.traj[i_traj]
        @async put!(traj_channel,loc_traj)
    end

    # Collect the resulting data
    for i_true = 1:data_size
        #Take from Channel:
        fvp_loc = take!(FVP_channel)
        FVP += fvp_loc
    end

    # Tell workers that its done
    NaN_mat = zeros(length(params.l_traj))
    fill!(NaN_mat,NaN)
    ending_traj = init_traj(params,params.l_traj)
    ending_traj.r_t = NaN_mat

    for j = 1:length(workers())
            @async put!(traj_channel, ending_traj)
    end

    #Close Channels
    sleep(0.1)
    close(traj_channel)
    close(FVP_channel)

    #Get the final Fisher result
    FVP /= data_len
    FVP  = FVP + cg_damping * v
    #Return the Gradients:
    return FVP
end

# Function running on the worker calculating the FIsher Matrix
@everywhere function distributed_fisher_VP(θ::Any, params::Params ,traj_channel::RemoteChannel, FVP_channel::RemoteChannel, v::Array{Float64,1})
    loc_traj = init_traj(params,params.l_traj)
    # Start loop where work is done
    while true
        # break loop if NaN state received
        try
            loc_traj = copy_traj(take!(traj_channel),params)
            if (sum(isnan, loc_traj.r_t) > 0)
                break
            end
        catch
            println("Some problem occured, worker $(myid()) exiting trajectory optimization function.")
            break
        end

        #execute the calculation
        loc_FVP = local_Fisher_VP(θ,loc_traj,params,v)

        put!(FVP_channel, loc_FVP)
    end

    return nothing
end

#calculate the local Fisher-Matrix
@everywhere function local_Fisher_VP(θ::Any,loc_traj::Trajectory,params::Params, v::Array{Float64,1})
    FVP = copy(v)*0
    l_traj = loc_traj.len

    for i_t = 1:l_traj
        #Take current state and build derivative
        s_t = loc_traj.s_t[i_t,:] #Current state
        a_t = loc_traj.a_t[i_t,:]
        Jt  = ∇logπ(θ,s_t,a_t,params)
        Jt_flat = net_to_vec(Jt)
        #Fisher Vector Product
        FVP += Jt_flat*dot(Jt_flat, v)
    end

    return FVP
end

#calculate the local Fisher-Matrix- Version 2
@everywhere function local_Fisher_VP2(θ::Any,traj::Trajectory,params::Params, v::Array{Float64,1})

    Jt      = ∇logπ_traj(θ, traj, params)
    Jt_flat = net_to_vec(Jt)
    FVP     = Jt_flat*dot(Jt_flat, v)

    return FVP
end

#Build the logarithmic gradient
@everywhere function logπ_traj(θ::Any, traj::Trajectory, params::Params)
    l_traj  = traj.len
    prob_pi = 0

    for i_t = 1:l_traj
        s_t = traj.s_t[i_t,:] #Current state
        a_t = traj.a_t[i_t,:] #Current action

        out_t    = mlp_nlayer_policy_norm(θ,s_t,params) #Evaluate the Network
        prob_pi += sum(logpdf(a_t,out_t,params))
    end

    return prob_pi / l_traj
end

@everywhere ∇logπ_traj = grad(logπ_traj,1)

function kl_all(θ,all_traj::Traj_Data,params::Params)
    n_traj  = all_traj.n_traj
    σ       = params.σ
    data_len = 0

    KL_mean  = 0

    #Size to be evaluating from:
    data_size = Int(floor(n_traj/2)) + 1 #here only on 10% of data...
    for i_traj = 1: data_size
        # As we regard costs => a negative Advantage is desired
        # Therefore we have a smaller loss is a negative advantage occurs
        l_traj = all_traj.traj[i_traj].len
        data_len += l_traj
        A_t = all_traj.traj[i_traj].A_t

        for i_t = 1: l_traj
            #Collect from struct
            s_t = all_traj.traj[i_traj].s_t[i_t,:]
            a_t = all_traj.traj[i_traj].a_t[i_t,:]
            #The two outputs:
            mu_old = all_traj.traj[i_traj].out_t[i_t,:]

            #Only for segway:
            s_mod = s_t
            mu_new = mlp_nlayer_policy_norm(θ,s_mod,params)  #Evaluate the Network
            #The two probabilities
            p_old = prod(all_traj.traj[i_traj].p_a_t[i_t,:])
            p_new = prod(get_π_prob(mu_new,σ,a_t))
            #Update the losses
            #KL_mean += p_new * log(p_new / p_old)#kl(mu_old,mu_new,σ)
            KL_mean += kl(mu_old,mu_new,σ)
        end
    end
    #Average
    KL_mean /= (data_len)

    return KL_mean
end

max_stars_repo_name: NikonPic/REINFORCEMENT-JULIA

max_issues_repo_name: NikonPic/REINFORCEMENT-JULIA

max_forks_repo_name: NikonPic/REINFORCEMENT-JULIA


text:
using QuasiGeostrophy, Plots, FFTW, BenchmarkTools
using LinearAlgebra
import Plots: plot
import QuasiGeostrophy: compute

##
Ωxy = S¹(0,2π) × S¹(0,2π)
Nx = 2^8; Ny = 2^8;
fourier_grid = create_grid((Nx, Ny), Ωxy)
x, y = fourier_grid.grid
fourier_transform = Transform(fourier_grid)

names = ("ϕ1", "ϕ2", "ϕ3", "ϕ4")
create_fields(names = names,
              grid = fourier_grid,
              transform = fourier_transform)
f1 = @. sin(x) + 0im * y
f2 = @. sin(y) + 0im * x
f3 = @. sin(x) * sin(y) # product
f4 = @. sin(x) + sin(y) # sum
f1 = fourier_transform.forward * f1
f2 = fourier_transform.forward * f2
f3 = fourier_transform.forward * f3
f4 = fourier_transform.forward * f4
ϕ1.data .= f1
ϕ2.data .= f2
ϕ3.data .= f3
ϕ4.data .= f4
## Check Algebra
norm((ϕ1 * ϕ2 - ϕ3).data)/norm((ϕ3).data)

f_ϕ1 = Field(ϕ1, BasicMetaData("ϕ1"))
f_ϕ2 = Field(ϕ2, BasicMetaData("ϕ2"))
tt =  2 * f_ϕ1 + f_ϕ2 * f_ϕ1 + 2 + tanh(f_ϕ1) + 2*2
compute(a::FourierField) = a
compute(tt)
evaluate(tt)

## Check Calculus
kx, ky = fourier_grid.wavenumbers
∂x = FourierOperator(im .* kx, FourierOperatorMetaData("∂x"))
∂y = FourierOperator(im .* ky, FourierOperatorMetaData("∂y"))

∂x(ϕ)

∂x(ϕ) + 1 + (∂x^2 + ∂y^2)(ϕ)
##
dmd = DerivativeMetaData(FourierOperator(im .* kx, FourierOperatorMetaData("∂x")), "x")
∂x = Operator(nothing, dmd)
∂x(f_ϕ1)
dmd = DerivativeMetaData(FourierOperator(im .* ky, FourierOperatorMetaData("∂y")), "y")
∂y = Operator(nothing, dmd)
∂ʸϕ1 = ∂y(f_ϕ2) * ∂x(f_ϕ1)
evaluate(∂ʸϕ1)
tmp = compute(∂ʸϕ1)
p1 = plot(tmp)

a1 = compute(∂x)^2
a2 = compute(∂x^2)
norm(a1.op - a2.op)

Δ1 = compute(∂x)^2 + compute(∂y)^2
Δ2 = compute(∂x^2+∂y^2)
norm(Δ1.op - Δ2.op)

ext: jl

max_stars_repo_name: upiterbarg/QuasiGeostrophy

max_issues_repo_name: upiterbarg/QuasiGeostrophy

max_forks_repo_name: upiterbarg/QuasiGeostrophy


text:
""" FS::CholeskySpecial{T,MT}

Special Cholesky type to allow for "removal" of rows and columns from a Cholesky factorization.

F0 = cholesky(A)
F = CholeskySpecial(F0)

Will destroy F0!

"Removing" row and column generates the factorization corresponding to
`A[:,i] = 0; A[i,:] = 0; A[i,i] = 1`

`deleterowcol!(F::CholeskySpecial, i)`
    Delete row and column i

`addrowcol!(F::CholeskySpecial, i)`
    Reverse the deletion row and column i

Supports `ldiv!`

Requires roughly n^2 + 3n extra memory
Update of row k requires O(n^2) computations
Solving is equivalent to a normal cholesky

"""
struct CholeskySpecial{T,MT} <: AbstractCholeskySpecial{T,MT}
    F::Cholesky{T,MT}
    idx::Set{Int}               # Indices of deleted rows and column
    tmp::Vector{T}              # tmp storage in deleterowcol
    tmp2::Vector{T}             # tmp storage in deleterowcol
    M::MT
end

function CholeskySpecial(F::Cholesky{T,MT}, M = F.U'F.U) where {T,MT}
    if F.uplo != 'U'
        error("Only implemented for U type")
    end
    CholeskySpecial{T,MT}(F, Set{Int}(),                    # F and idx
        fill(zero(T), size(F,1)), fill(zero(T), size(F,1)), # tmp, tmp2
        M)                                                  # M
end

"""
ldiv!(F::CholeskySpecial, b)

Solve Mx=b where F is Cholesky factorization if M with some rows and columns set to identity
"""
LinearAlgebra.ldiv!(F::CholeskySpecial{T,MT}, b::AbstractVector{T}) where {T,MT} =
    ldiv!(F.F, b)

Base.size(F::CholeskySpecial) = size(F.F)
Base.size(F::CholeskySpecial, i::Int) = size(F.F)[i]

function Base.getproperty(F::CholeskySpecial{T,MT}, s::Symbol) where {T,MT}
    if s == :U
        return F.F.U
    else
        getfield(F,s)
    end
end

""" `deleterowcol!(F::CholeskySpecial{T}, i)`
    Given F=U*U', update factorization corresponding to setting
    F[i, :] = 0, F[:, i] = 0
    F[i,i] = 1.
"""
function deleterowcol!(F::CholeskySpecial{T,MT}, i) where {T,MT}
    if i in F.idx
        error("Can not delete row and column that is already deleted")
    end
    n = size(F.F,1)
    # Create vector for rank update
    F.tmp[1:i] .= zero(T)
    F.tmp[(i+1):n] .= F.F.U[i, (i+1):n]
    # Record that  this row/column is "removed"
    push!(F.idx, i)
    # Do rank update
    lowrankupdate!(F.F, F.tmp)
    # Zero out
    for i in F.idx
        F.F.U[1:(i-1),i] .= zero(T)
        F.F.U[i,(i+1):end] .= zero(T)
        F.F.U[i,i] = one(T)
    end
    return
end

"""  Reverse deleterowcol!
"""
function addrowcol!(F::CholeskySpecial{T,MT}, i) where {T,MT}
    if !(i in F.idx)
        error("Can not add row and column that has not been deleted")
    end
    n = size(F.F,1)
    b = F.tmp
    b .= F.M[:,i]
    b[i:end] .= zero(T)
    for j in F.idx
        #This column should be zero in M
        b[j] = zero(T)
    end
    ldiv!(F.F.U', b) # b is now [S12; junk]
    S12 = view(b, 1:(i-1))
    # Scalar result
    S22 = sqrt(F.M[i,i] - S12'S12)

    # Get row from M (under assumption that some rows/cols are identity)
    M23 = view(F.tmp2, (i+1):n)
    M23 .= F.M[i,(i+1):end]
    for j in F.idx
        if j > i
            #This column should be zero in M
            F.tmp2[j] = zero(T)
            # Equivalent to
            # M23[j-i] = zero(T)
        end
    end
    # TODO see if view on factor matrix directly is faster
    # Now set the new row to
    # S23 = (M23 - F.U[1:(i-1),(i+1):end]'S12)/S22 # Transposed
    S23 = view(F.F.U, i, (i+1):n) # This is where we write
    S13 = view(F.F.U, 1:(i-1), (i+1):n)
    mul!(S23, S13', S12)
    S23 .= (M23 .- S23)./S22

    # Set the diagonal
    F.F.U[i,i] = S22
    # And the column
    F.F.U[1:(i-1),i] .= S12

    # No longer need S12, so reuse tmp for rank update
    F.tmp[1:i] .= zero(T)
    F.tmp[(i+1):n] .= S23

    #Update the triangular part S33
    lowrankdowndate!(F.F, F.tmp)

    # Rememeber that i is added back
    pop!(F.idx, i)
    return
end

max_stars_repo_name: mfalt/QPDAS

max_issues_repo_name: mfalt/QPDAS

max_forks_repo_name: mfalt/QPDAS


text:
const observables = [
    "Phosphorylated_Akt",
    "Phosphorylated_ERK",
    "Phosphorylated_c-Myc",
]

function observables_index(observable_name::String)::Union{Int,Nothing}
    if !(observable_name in observables)
        error("$observable_name is not defined in observables.")
    end
    return findfirst(isequal(observable_name), observables)
end


text:
array = zeros(Int64, 3, 3, 3)

array[:,:,1] = [1, 2, 3]
println(array)


text:
module CoreToggle

using ..Ahorn, Maple

const placements = Ahorn.PlacementDict(
    "Core Mode Toggle (Fire)" => Ahorn.EntityPlacement(
        Maple.CoreFlag,
        "point",
        Dict{String, Any}(
            "onlyFire" => true
        )
    ),
    "Core Mode Toggle (Ice)" => Ahorn.EntityPlacement(
        Maple.CoreFlag,
        "point",
        Dict{String, Any}(
            "onlyIce" => true
        )
    ),
    "Core Mode Toggle (Both)" => Ahorn.EntityPlacement(
        Maple.CoreFlag
    ),
)

function switchSprite(entity::Maple.CoreFlag)
    onlyIce = get(entity.data, "onlyIce", false)
    onlyFire = get(entity.data, "onlyFire", false)

    if onlyIce
        return "objects/coreFlipSwitch/switch13.png"

    elseif onlyFire
        return "objects/coreFlipSwitch/switch15.png"

    else
        return "objects/coreFlipSwitch/switch01.png"
    end
end

function Ahorn.selection(entity::Maple.CoreFlag)
    x, y = Ahorn.position(entity)
    sprite = switchSprite(entity)

    return Ahorn.getSpriteRectangle(sprite, x, y)
end

function Ahorn.render(ctx::Ahorn.Cairo.CairoContext, entity::Maple.CoreFlag, room::Maple.Room)
    sprite = switchSprite(entity)

    Ahorn.drawSprite(ctx, sprite, 0, 0)
end

end


text:
@testset "combinedims" begin
    #m = [1 2 3 4; 5 6 7 8; 9 10 11 12]
    m = [1 5 9; 2 6 10; 3 7 11; 4 8 12]
    vv = [[1,2,3,4],[5,6,7,8],[9,10,11,12]]
    @test combinedims(vv) == m

    @test combinedimsview(vv) == m
end


text:
# ---
# title: 1518. Water Bottles
# id: problem1518
# author: Tian Jun
# date: 2020-10-31
# difficulty: Easy
# categories: Greedy
# link: <https://leetcode.com/problems/water-bottles/description/>
# hidden: true
# ---
#
# Given `numBottles` full water bottles, you can exchange `numExchange` empty
# water bottles for one full water bottle.
#
# The operation of drinking a full water bottle turns it into an empty bottle.
#
# Return the **maximum** number of water bottles you can drink.
#
#
#
# **Example 1:**
#
# **![](https://assets.leetcode.com/uploads/2020/07/01/sample_1_1875.png)**
#
#
#
#     Input: numBottles = 9, numExchange = 3
#     Output: 13
#     Explanation: You can exchange 3 empty bottles to get 1 full water bottle.
#     Number of water bottles you can drink: 9 + 3 + 1 = 13.
#
#
# **Example 2:**
#
# ![](https://assets.leetcode.com/uploads/2020/07/01/sample_2_1875.png)
#
#
#
#     Input: numBottles = 15, numExchange = 4
#     Output: 19
#     Explanation: You can exchange 4 empty bottles to get 1 full water bottle.
#     Number of water bottles you can drink: 15 + 3 + 1 = 19.
#
#
# **Example 3:**
#
#
#
#     Input: numBottles = 5, numExchange = 5
#     Output: 6
#
#
# **Example 4:**
#
#
#
#     Input: numBottles = 2, numExchange = 3
#     Output: 2
#
#
#
#
# **Constraints:**
#
#   * `1 <= numBottles <= 100`
#   * `2 <= numExchange <= 100`
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
function load_job(job_dir::AbstractString)
    orig_dir = job_dir
    s = Server("localhost")
    if !isabspath(job_dir)
        job_dir = joinpath(s, job_dir)
    end
    scriptpath = joinpath(job_dir, "job.tt")
    if ispath(scriptpath)
        version = Jobs.version(job_dir)
    else
        error("No valid job found in $job_dir.")
    end
    metadata = ispath(joinpath(job_dir, ".metadata.jld2")) ? JLD2.load(joinpath(job_dir, ".metadata.jld2"))["metadata"] : Dict{Symbol, Any}()

    name, tcalcs, header, environment = FileIO.read_job_script(scriptpath)
    calculations, structure = FileIO.parse_calculations(tcalcs)
    job = Job(name         = name,
              dir          = orig_dir,
              version      = version,
              metadata     = metadata,
              structure    = structure,
              calculations = calculations,
              header       = header,
              environment  = environment)
    return job
end

job_versions(args...) = Jobs.job_versions(args...)
registered_jobs(args...) = Jobs.registered_jobs(args...)
running_jobs(args...) = Jobs.running_jobs(args...)

function workflow_logger(job::Job)
    return TeeLogger(MinLevelLogger(FileLogger(joinpath(job, ".workflow", "info.log");
                                               append = true), Logging.Info),
                     MinLevelLogger(FileLogger(joinpath(job, ".workflow", "info.log");
                                               append = true), Logging.Warn),
                     MinLevelLogger(FileLogger(joinpath(job, ".workflow", "error.log");
                                               append = true), Logging.Error))
end

queued_dir(job::Job, args...) = joinpath(job, ".workflow/queued", args...)
finished_dir(job::Job, args...) = joinpath(job, ".workflow/finished", args...)

function queue_steps(job::Job, funcs)
    qd = queued_dir(job)
    if !ispath(qd)
        mkpath(qd)
    end
    prev_steps = readdir(qd)
    last_step = !isempty(prev_steps) ? parse(Int, splitext(prev_steps[end])[1]) : 0
    for (i, f) in enumerate(funcs)
        write(joinpath(qd, "$(i + last_step).jl"), @code_string f(job, Dict{Symbol,Any}()))
    end
end

function write_workflow_files(job::Job)
    all = string.(values(Base.loaded_modules))
    valid = String[]
    ks = keys(Pkg.project().dependencies)
    for p in all
        if p ∈ ks
            push!(valid, p)
        end
    end
    JLD2.jldsave(joinpath(job, ".workflow/environment.jld2"); modules = valid,
                        project = Base.current_project())
    return JLD2.jldsave(joinpath(job, ".workflow/ctx.jld2"); ctx = Dict{Symbol,Any}())
end

function clear_queue!(job::Job)
    qd = queued_dir(job)
    if !ispath(qd)
        return
    end
    for f in readdir(qd)
        rm(joinpath(qd, f))
    end
end

queued(job::Job) = readdir(queued_dir(job))
finished(job::Job) = readdir(finished_dir(job))

# function submit_workflow(job::Job, funcs, d::Daemon = init_daemon())
#     queue_steps(job, funcs)
#     write_workflow_files(job)
#     while !is_started(d)
#         sleep(1)
#     end
#     return runexpr("""
#                DFControl.spawn_worker(DAEMON, Job("$(job.dir)"))
#                DFControl.save(DAEMON)
#                """; port = d.port)
# end

mods_test() = Base.loaded_modules

"""
    save(job::Job)

Saves the job's calculations and `job.tt` submission script in `job.dir`.
Some sanity checks will be performed on the validity of flags, execs, pseudopotentials, etc.
The job will also be registered for easy retrieval at a later stage.

If a previous job is present in the job directory (indicated by a valid job script),
it will be copied to the `.versions` sub directory as the previous version of `job`,
and the version of `job` will be incremented.
"""
function save(jobdir::String, files; kwargs...)

    #Since at this stage we know the job will belong to the current localhost we change the server
    # Here we find the main directory, needed for if a job's local dir is a .versions one
    dir = Jobs.main_job_dir(jobdir)
    if ispath(joinpath(dir, "job.tt"))
        tj = load_job(dir)
        cp(tj, joinpath(tj, Jobs.VERSION_DIR_NAME, "$(tj.version)"); force = true)
    end
    if !occursin(jobdir, dir)
        # We know for sure it was a previously saved job
        # Now that we have safely stored it we can clean out the directory to then fill
        # it with the files from the job.version
        clean_dir!(dir)
        for f in readdir(jobdir)
            cp(f, dir; force = true)
        end
    end

    # Needs to be done so the inputs `dir` also changes.
    mkpath(dir)

    version = Jobs.last_job_version(dir) + 1
    for (name, f) in files
        write(joinpath(dir, name), f)
    end

    Jobs.maybe_register_job(jobdir)
    return version
end

"""
    submit(job::Job; kwargs...)

First saves the job, then tries to submit the job script through `sbatch job.tt` or `bash job.tt` if the former command fails.
`kwargs...` get passed to `save(job; kwargs...)`.
"""
function submit(job_dir::String)
    open(PENDING_JOBS_FILE, "a", lock=true) do f
        return write(f, job_dir * "\n")
    end
end

# The actual submitting
function submit(s::Server, job_dir::String)
    if s.scheduler == Servers.Bash
        return bash_submit(job_dir)
    elseif s.scheduler == Servers.Slurm
        return slurm_submit(job_dir)
    end
end

"""
    last_running_calculation(path::String)

Returns the last `Calculation` for which an output file was created.
"""
function last_running_calculation(path::String)
    scrpath = joinpath(path, "job.tt")
    job = load_job(path)
    t = mtime(Jobs.scriptpath(job))
    times = map(x -> (o = joinpath(job, x.outfile); ispath(o) ? mtime(o) : 0.0), job.calculations)
    return isempty(times) ? nothing : findmax(times)[2]
end

"""
    state(job_dir::String)

Returns the job state of the job in `job_dir`.
"""
state(job_dir::String) =
    haskey(JOB_QUEUE[], job_dir) ? JOB_QUEUE[][job_dir][2] : Jobs.Unknown

function dirsize(path::String)
    totsize = 0.0
    for (root, dirs, files) in walkdir(path)
        for file in files
            totsize += filesize(root, file)
        end
    end
    return totsize
end

"""
    filesize(job::Job)

Total filesize on disk for a job and all its versions.
"""
Base.filesize(job::Job) = dirsize(job.dir)

"""
    cleanup(job::Job)

Cleanup `job.dir` interactively.
"""
function cleanup(job::Job)
    labels = String[]
    paths = String[]
    for v in versions(job)
        vpath = version_dir(job, v)
        s = round(dirsize(vpath) / 1e6; digits = 3)
        push!(labels, "Version $v:  $s Mb")
        push!(paths, vpath)
        opath = joinpath(vpath, Jobs.TEMP_CALC_DIR)
        if ispath(opath)
            s_out = round(dirsize(opath) / 1e6; digits = 3)
            push!(labels, "Version $v/outputs:  $s_out Mb")
            push!(paths, opath)
        end
    end
    menu = MultiSelectMenu(labels)
    choices = request("Select job files to delete:", menu)
    for i in choices
        if ispath(paths[i]) # Could be that outputs was already deleted
            @info "Deleting $(paths[i])"
            rm(paths[i]; recursive = true)
        end
    end
end

function save_metadata(job)
    return jldsave(joinpath(job, ".metadata.jld2"); metadata = job.metadata,
                   version = job.version)
end

timestamp(job::Job) = timestamp(job.dir)
has_timestamp(job) = haskey(job.metadata, :timestamp)

function clean_dir!(dir::AbstractString)
    for f in readdir(dir)
        if f == Jobs.TEMP_CALC_DIR ||
           f == Jobs.VERSION_DIR_NAME ||
           splitext(f)[end] == ".jl"
            continue
        end
        rm(joinpath(dir, f); recursive = true)
    end
end

exists_job(d::AbstractString) = ispath(d) && ispath(joinpath(d, "job.tt"))

"Finds the output files for each of the calculations of a job, and groups all found data into a dictionary."
function outputdata(jobdir::String, calculations::Vector{String})
    job = load_job(jobdir)
    calculations = isempty(calculations) ? map(x->x.name, job.calculations) : calculations
    respath = joinpath(job, "results.jld2")
    if ispath(respath)
        datadict = JLD2.load(respath, "outputdata")
    else
        datadict = Dict{String,Dict{Symbol,Any}}()
    end
    stime = isempty(datadict) ? 0.0 : mtime(respath)
    new_data = false
    for c in calculations
        calculation = job[c]
        p = joinpath(job, calculation.outfile)
        if mtime(p) > stime
            tout = outputdata(calculation, p)
            if !isempty(tout)
                datadict[calculation.name] = tout
                new_data = true
            end
        end
    end
    if new_data
        JLD2.jldsave(respath; outputdata=datadict)
        return respath
    elseif ispath(respath)
        return respath
    else
        return nothing
    end
end

rm_version!(jobdir::String, version::Int) = Jobs.rm_version!(load_job(job), version)

add_environment(env::Environment, name::AbstractString) = Jobs.save(env, name)
function get_environment(name::AbstractString)
    out = Jobs.load_environment(name)
    if out === nothing
        error("No Environment found with name $name")
    end
    return out
end

rm_environment!(args...) = Jobs.rm_environment!(args...)

queue(args...) = queue!(Dict{String, Tuple{Int, Jobs.JobState}}(), args...)

function queue!(q, s::Server, init=false)
    if s.scheduler == Servers.Bash
        return bash_queue!(q, init)
    elseif s.scheduler == Servers.Slurm
        return slurm_queue!(q, init)
    end
end


text:
#Run analysis with non-neutral data with 1 clone, input parameters:
# mu = 20.0
# cm = 200
# cellularity = 0.8
println()
println("####################")
println("Checking inference on non-neutral simulated data with known input parameters returns ground truth...")
Random.seed!(1)
@time out = fitABCmodels("data/oneclone.txt",
  "oneclone",
  read_depth = 150,
  resultsdirectory = "output",
  nparticles = 100,
  maxiterations =  10^4,
  Nmax = 10^3,
  maxclones = 2,
  verbose = true,
  firstpass = false,
  save = true,
  progress = true);

# check if we get the correct model
println("\tTesting posterior model probability returns one clone as most probable model...")
@test argmax(out.ModelProb[:Probability]) == 2

println("\tTesting 2 clone model did not die out...")
@test out.ModelProb[:Probability][3] > 0.0

# extract parameters
mu = out.Posterior[2].Parameters[:mu]
clonalmutations = out.Posterior[2].Parameters[:clonalmutations]
cellularity = out.Posterior[2].Parameters[:cellularity]

println("\tChecking true parameters are within the 95% credible interval range...")
# check if parameters are within 80% credible interval
@test quantile(mu, 0.025) < 20.0 < quantile(mu, 0.975)
@test quantile(clonalmutations, 0.025) < 200.0 < quantile(clonalmutations, 0.975)
@test quantile(cellularity, 0.025) < 0.8 < quantile(cellularity, 0.975)

println("\tChecking plotting functions work and are saved...")
saveallplots(out)
@test isfile("output/oneclone/finalpopulation/plots/oneclone-histogram-0clone.pdf")
@test isfile("output/oneclone/finalpopulation/posterior/oneclone-histogram-clone0.csv")
rm("output", recursive = true)

println("All tests passed for single clone.")


text:
module LMDB_CUR
    using LMDB
    using Base.Test

    const dbname = "testdb"
    key = 10
    val = "key value is "

    # Create dir
    mkdir(dbname)

    # Procedural style
    env = create()
    try
        open(env, dbname)
        txn = start(env)
        dbi = open(txn)
        commit(txn)

        txn = start(env)
        cur = open(txn, dbi)
        try
            @test 0 == put!(cur, key+1, val*string(key+1))
            @test 0 == put!(cur, key, val*string(key))
        finally
            close(cur)
            commit(txn)
        end
        @test !isopen(cur)
        @test !isopen(txn)
    finally
        close(env)
    end
    @test !isopen(env)

    # Remove db dir
    rm(dbname, recursive=true)
end


text:
@test ExactDiag.permutation_parity(Int[]) == 0
@test ExactDiag.permutation_parity([1]) == 0
@test ExactDiag.permutation_parity([1,2]) == 0
@test ExactDiag.permutation_parity([2,1]) == 1
@test ExactDiag.permutation_parity([1,2,3]) == 0
@test ExactDiag.permutation_parity([3,1,2]) == 0
@test ExactDiag.permutation_parity([1,3,2]) == 1

@test ExactDiag.site_spinflip(0) == 0
@test ExactDiag.site_spinflip(1) == 2
@test ExactDiag.site_spinflip(2) == 1
@test ExactDiag.site_spinflip(3) == 3

function test_1d_hubbard_hamiltonian(lattice)
    apply_hamiltonian = hubbard_hamiltonian(t=1, U=3, ϵ_total_spin=pi/1000, ϵ_total_pseudospin=MathConstants.e/1000)
    hs = HubbardHilbertSpace(lattice)
    seed_state!(hs, N_up=div(length(lattice), 2), N_dn=div(length(lattice), 2))
    mat = operator_matrix(hs, apply_hamiltonian)
    @test ishermitian(mat)
    zzz = HilbertSpaceTranslationCache(hs, 1)
    for j in 1:length(hs.indexer)
        i, η = translateη(zzz, j)
        # FIXME: test value of η
        @test hs.indexer[i] == my_1d_translate(hs.indexer[j])
        debug && println("$(hs.indexer[j])\t$(hs.indexer[i])\t$η")
    end
    # FIXME: test GS energy
end
test_1d_hubbard_hamiltonian(ChainLattice([8]))

# With abelian symmetries (this one always assumes half filling)
function test_1d_hubbard_symmetries(lattice, symmetries::Vector{Tuple{F,Int}} where {F<:Function})
    L = length(lattice)
    hs = HubbardHilbertSpace(lattice)
    seed_state!(hs, N_up=div(L, 2), N_dn=div(L, 2))
    apply_hamiltonian = hubbard_hamiltonian(t=1, U=2)
    rst = RepresentativeStateTable(hs, apply_hamiltonian, symmetries)

    full_ham = operator_matrix(hs, apply_hamiltonian)

    processed_length = 0
    symmbounds = (repeat([2], inner=[length(symmetries)])...,)
    for k_idx in eachmomentumindex(hs.lattice)
        for symm_idx in 1:2*length(symmetries)
            symm = [Tuple(CartesianIndices(symmbounds)[symm_idx])...] .- 1
            diagsect = DiagonalizationSector(rst, 1, k_idx, symm)
            length(diagsect) != 0 || continue
            processed_length += length(diagsect)
            reduced_ham = Matrix(construct_reduced_hamiltonian(diagsect))
            @test norm(reduced_ham - reduced_ham') < 1e-5
            fact = eigen(Hermitian((reduced_ham + reduced_ham') / 2))
            evals, evecs = fact.values, fact.vectors
            for i in 1:length(evals)
                eval = evals[i]
                evec = evecs[:,i]
                ψ = get_full_psi(diagsect, evec)
                @test eigenstate_badness(full_ham, eval, ψ) ≈ 0 atol=1e-8
            end
        end
    end
    @test processed_length == length(hs.indexer)
end

let
    for L in 2:2:6
        # Take advantage of spinflip symmetry
        test_1d_hubbard_symmetries(ChainLattice([L]), [spinflip_symmetry])
        test_1d_hubbard_symmetries(ChainLattice([L], diagm(0 => [L]), [1//2]), [spinflip_symmetry])
        test_1d_hubbard_symmetries(ChainLattice([L], diagm(0 => [L]), [1//3]), [spinflip_symmetry])
        test_1d_hubbard_symmetries(ChainLattice([L], diagm(0 => [0])), [spinflip_symmetry])

        # Particle-hole symmetry
        test_1d_hubbard_symmetries(ChainLattice([L]), [particlehole_symmetry])
        test_1d_hubbard_symmetries(ChainLattice([L], diagm(0 => [L]), [1//2]), [particlehole_symmetry])
        test_1d_hubbard_symmetries(ChainLattice([L], diagm(0 => [0])), [particlehole_symmetry])

        # Spinflip + particle-hole
        test_1d_hubbard_symmetries(ChainLattice([L]), [spinflip_symmetry, particlehole_symmetry])
        test_1d_hubbard_symmetries(ChainLattice([L], diagm(0 => [L]), [1//2]), [spinflip_symmetry, particlehole_symmetry])
        test_1d_hubbard_symmetries(ChainLattice([L], diagm(0 => [0])), [spinflip_symmetry, particlehole_symmetry])
    end
end

let
    lattice = ChainLattice([6], diagm(0 => [6]), [1//7])
    ltrc = LatticeTranslationCache(lattice, 1)
    hs = HubbardHilbertSpace(lattice)
    push!(hs.indexer, [0, 0, 3, 1, 0, 2])
    t1, η1 = translateη(hs, ltrc, 1)
    @test t1 == findfirst(isequal([2, 0, 0, 3, 1, 0]), hs.indexer) == 2
    @test η1 == 5//14
    t2, η2 = translateη(hs, ltrc, 2)
    @test t2 == findfirst(isequal([0, 2, 0, 0, 3, 1]), hs.indexer) == 3
    @test η2 == 0//1
    t3, η3 = translateη(hs, ltrc, 3)
    @test t3 == findfirst(isequal([1, 0, 2, 0, 0, 3]), hs.indexer) == 4
    @test η3 == 5//14
    t4, η4 = translateη(hs, ltrc, 4)
    @test t4 == findfirst(isequal([3, 1, 0, 2, 0, 0]), hs.indexer) == 5
    @test η4 == -2//7
    t5, η5 = translateη(hs, ltrc, 5)
    @test t5 == findfirst(isequal([0, 3, 1, 0, 2, 0]), hs.indexer) == 6
    @test η5 == 0//1
    t6, η6 = translateη(hs, ltrc, 6)
    @test t6 == 1
    @test η6 == 0//1
end

# explicit Slater determinant for ground state
let
    L = 6
    lattice = ChainLattice([L])
    apply_hamiltonian = hubbard_hamiltonian(t=1)
    hs = HubbardHilbertSpace(lattice)
    seed_state!(hs, N_up=div(L, 2), N_dn=div(L, 2))

    rst = RepresentativeStateTable(hs, apply_hamiltonian, [spinflip_symmetry])

    # Find the ground state wavefunction via ExactDiag
    local gs_eval, gs_evec
    gs_eval = maxintfloat(Float64) # fixme
    for k_idx in eachmomentumindex(hs.lattice)
        for spinflip_idx in 0:1
            diagsect = DiagonalizationSector(rst, 1, k_idx, [spinflip_idx])
            length(diagsect) != 0 || continue
            reduced_ham = Matrix(construct_reduced_hamiltonian(diagsect))
            @test norm(reduced_ham - reduced_ham') < 1e-5
            fact = eigen(Hermitian((reduced_ham + reduced_ham') / 2))
            evals, evecs = fact.values, fact.vectors
            eval = evals[1]
            evec = evecs[:,1]
            ψ = get_full_psi(diagsect, evec)
            if eval < gs_eval
                gs_eval = eval
                gs_evec = ψ
            end
        end
    end
    debug && @show gs_eval

    # Construct the Slater determinant wavefunction
    slater = Array{ComplexF64}(undef, length(hs.indexer))
    for (i, state) in enumerate(hs.indexer)
        # First figure out the positions of the particles
        pn = 3
        r_up = findall(state) do x
            x & 1 != 0
        end
        @assert length(r_up) == pn
        r_dn = findall(state) do x
            x & 2 != 0
        end
        @assert length(r_dn) == pn

        # Now construct and calculate the two determinants.
        mat1 = Array{ComplexF64}(undef, pn, pn)
        mat2 = Array{ComplexF64}(undef, pn, pn)
        for j in 1:3
            mat1[j, 1] = exp(im * pi * r_up[j] / 3.)
            mat1[j, 2] = 1
            mat1[j, 3] = exp(-im * pi * r_up[j] / 3.)
            mat2[j, 1] = exp(im * pi * r_dn[j] / 3.)
            mat2[j, 2] = 1
            mat2[j, 3] = exp(-im * pi * r_dn[j] / 3.)
        end

        # Normalize properly
        mat1 /= √L
        mat2 /= √L

        slater[i] = det(mat1) * det(mat2)
    end

    overlap = dot(gs_evec, slater)
    debug && @show overlap
    @test abs(overlap) ≈ 1
end

# Slater determinants for all eigenstates
function degenerate_ranges(v::AbstractVector{T}, tol::T) where {T<:Real}
    @assert issorted(v)
    diffs = v[2:end] - v[1:end-1]
    indices = findall(diffs) do d
        d > tol
    end
    return UnitRange{Int}[x:y for (x, y) in zip([1; indices .+ 1], [indices; length(v)])]
end

@test degenerate_ranges([1,2,3,4.0], 0.5) == UnitRange{Int}[1:1, 2:2, 3:3, 4:4]
@test degenerate_ranges([1,2,3,4.0], 1.2) == UnitRange{Int}[1:4]
@test degenerate_ranges([1,2,3,5.0], 1.2) == UnitRange{Int}[1:3, 4:4]
@test degenerate_ranges([0,2,3,4.0], 1.2) == UnitRange{Int}[1:1, 2:4]
@test degenerate_ranges([0,2,3,5.0], 1.2) == UnitRange{Int}[1:1, 2:3, 4:4]

function test_slater_determinants(lattice::AbstractLattice, N_up::Int, N_dn::Int, ϵ)
    apply_hamiltonian = hubbard_hamiltonian(t=1)
    hs = HubbardHilbertSpace(lattice)
    seed_state!(hs, N_up=N_up, N_dn=N_dn)

    rst = RepresentativeStateTable(hs, apply_hamiltonian)

    # Determine the energy and momentum of each possible band filling.  Store
    # everything in the `band_fillings` dict, which is keyed by total momentum.
    @assert isbravais(lattice)
    band_fillings = Dict{Vector{Rational{Int}}, Vector{Tuple{Float64, Vector{Int}, Vector{Int}}}}()
    for filled_k_up_indices in combinations(eachmomentumindex(lattice), N_up)
        for filled_k_dn_indices in combinations(eachmomentumindex(lattice), N_dn)
            filled_orbitals = map(k_idx -> momentum(lattice, k_idx),
                                  [filled_k_up_indices; filled_k_dn_indices])
            energy = mapreduce(ϵ, +, filled_orbitals)
            total_momentum = rem.(reduce(+, filled_orbitals), 1)
            v = get!(band_fillings, total_momentum) do
                Tuple{Float64, Vector{Int}, Vector{Int}}[]
            end
            push!(v, (energy, filled_k_up_indices, filled_k_dn_indices))
        end
    end
    @assert length(band_fillings) == nmomenta(lattice)

    # Sort by energy at each momentum.  This will help us soon determine
    # nearly-adjacent energy levels.
    for v in values(band_fillings)
        sort!(v, by=(x -> x[1]))
    end
    @assert mapreduce(length, +, values(band_fillings)) == length(hs.indexer)

    total_momenta = [rem.(momentum(lattice, k_idx, N_up + N_dn), 1) for k_idx in eachmomentumindex(lattice)]

    root_V = sqrt(length(lattice))

    # Construct each Slater determinant wavefunction
    slater_wfs = [Array{ComplexF64}(undef, length(hs.indexer), length(band_fillings[k])) for k in total_momenta]
    for (i, state) in enumerate(hs.indexer)
        # First figure out the positions of the particles
        r_up = findall(x -> x & 1 != 0, state)
        r_dn = findall(x -> x & 2 != 0, state)
        @assert length(r_up) == N_up
        @assert length(r_dn) == N_dn

        # Calculate each Slater determinant at these positions
        for (total_k, slater_wfs_k) in zip(total_momenta, slater_wfs)
            for (j, (energy, filled_k_up_indices, filled_k_dn_indices)) in enumerate(band_fillings[total_k])
                # Construct and calculate the two determinants with proper normalization.
                up_mat = [cis(kdotr(lattice, k, r)) / root_V for r in r_up, k in filled_k_up_indices]
                dn_mat = ComplexF64[cis(kdotr(lattice, k, r)) / root_V for r in r_dn, k in filled_k_dn_indices]
                @test eltype(up_mat) == eltype(dn_mat) == ComplexF64
                slater_wfs_k[i, j] = det(up_mat) * det(dn_mat)
            end
        end
    end

    # Diagonalize each momentum sector, and check the energies and overlaps
    for k_idx in eachmomentumindex(hs.lattice)
        slater_energies = [energy for (energy, ku, kd) in band_fillings[total_momenta[k_idx]]]
        diagsect = DiagonalizationSector(rst, 1, k_idx)
        length(diagsect) != 0 || continue
        reduced_ham = Matrix(construct_reduced_hamiltonian(diagsect))
        @test norm(reduced_ham - reduced_ham') < 1e-5
        fact = eigen(Hermitian((reduced_ham + reduced_ham') / 2))
        evals, evecs = fact.values, fact.vectors

        # Check the energies
        @test evals ≈ slater_energies atol=1e-10

        # Check the overlaps
        i_check = 0
        for degenerate_range in degenerate_ranges(slater_energies, 1e-8)
            for i in degenerate_range
                i_check += 1
                @assert i == i_check

                ψ = get_full_psi(diagsect, @view evecs[:,i])
                # Project ψ into the subspace of Slater determinant eigenstates
                # at this energy.  It should remain unchanged by this
                # projection.
                ϕ = zero(ψ)
                for j in degenerate_range
                    slater_wf = @view slater_wfs[k_idx][:,j]
                    ϕ .+= slater_wf .* dot(slater_wf, ψ)
                end
                @test ϕ ≈ ψ atol=1e-10
            end
        end
        @test i_check == length(evals)
    end
end

let
    hypercubic_ϵ(k) = -2 * sum(cos(2π * kα) for kα in k)

    test_slater_determinants(ChainLattice([6]), 2, 0, hypercubic_ϵ)
    test_slater_determinants(ChainLattice([6]), 3, 3, hypercubic_ϵ)
    test_slater_determinants(ChainLattice([6], diagm(0 => [6]), [1//2]), 1, 0, hypercubic_ϵ)
    test_slater_determinants(ChainLattice([6], diagm(0 => [6]), [1//5]), 1, 0, hypercubic_ϵ)
    test_slater_determinants(ChainLattice([6], diagm(0 => [6]), [1//2]), 3, 3, hypercubic_ϵ)
    test_slater_determinants(ChainLattice([6], diagm(0 => [6]), [1//5]), 3, 3, hypercubic_ϵ)
    test_slater_determinants(ChainLattice([5]), 2, 3, hypercubic_ϵ)
    test_slater_determinants(SquareLattice([2, 3]), 3, 3, hypercubic_ϵ)

    #test_slater_determinants(TriangularLattice([2, 3]), 3, 3)
    #test_slater_determinants(TriangularLattice([2, 3], diagm(0 => [2,3]), [1//2, 1//3]), 3, 3)
end

function test_hubbard_abelian_spinflip(lattice, N_updn; t=1, U=2, kwargs...) # does not assume half filling
    hs = HubbardHilbertSpace(lattice)
    seed_state!(hs, N_up=N_updn, N_dn=N_updn)
    apply_hamiltonian = hubbard_hamiltonian(t=t, U=U; kwargs...)
    rst = RepresentativeStateTable(hs, apply_hamiltonian, [spinflip_symmetry])

    full_ham = operator_matrix(hs, apply_hamiltonian)

    processed_length = 0
    for k_idx in eachmomentumindex(hs.lattice)
        for spinflip_idx in 0:1
            diagsect = DiagonalizationSector(rst, 1, k_idx, [spinflip_idx])
            length(diagsect) != 0 || continue
            processed_length += length(diagsect)
            reduced_ham = Matrix(construct_reduced_hamiltonian(diagsect))
            @test norm(reduced_ham - reduced_ham') < 1e-5
            fact = eigen(Hermitian((reduced_ham + reduced_ham') / 2))
            evals, evecs = fact.values, fact.vectors
            for i in 1:length(evals)
                eval = evals[i]
                evec = evecs[:,i]
                ψ = get_full_psi(diagsect, evec)

                @test eigenstate_badness(full_ham, eval, ψ) ≈ 0 atol=1e-8

                if i == 1
                    let L = length(lattice)
                        for L_A in 0:div(L, 2)
                            ent_cut1 = entanglement_entropy(Tracer(hs, 1:L_A), ψ)
                            ent_cut2 = entanglement_entropy(Tracer(hs, 1:L-L_A), ψ)
                            if ndimensions(hs.lattice) == 1
                                @test ent_cut1 ≈ ent_cut2 atol=1e-8
                            end
                        end
                    end
                end
            end
        end
    end
    @test processed_length == length(hs.indexer)
end

test_hubbard_abelian_spinflip(ChainLattice([6]), 2)
test_hubbard_abelian_spinflip(ChainLattice([6]), 3)
test_hubbard_abelian_spinflip(ChainLattice([2], diagm(0 => [2]), [1//3]), 1)
test_hubbard_abelian_spinflip(ChainLattice([4], diagm(0 => [4]), [1//3]), 2)
test_hubbard_abelian_spinflip(ChainLattice([6], diagm(0 => [6]), [1//2]), 3)
test_hubbard_abelian_spinflip(ChainLattice([6], diagm(0 => [6]), [1//5]), 3)
test_hubbard_abelian_spinflip(ChainLattice([6], diagm(0 => [6]), [1//5]), 3, t2=0.2)
test_hubbard_abelian_spinflip(ChainLattice([4], diagm(0 => [4]), [1//5]), 2, t2=0.2, J=0.3)

test_hubbard_abelian_spinflip(SquareLattice([2,3]), 3)
test_hubbard_abelian_spinflip(TriangularLattice([2,3]), 3)
# NOTE: The cubic lattice does not (yet) implement next-nearest
# neighbors, so this should fail.
@test_throws MethodError test_hubbard_abelian_spinflip(CubicLattice([2,3,2]), 3, t2=0.2)


text:
_inner_type(l::SimpleSDMResponse{T}) where {T <: Any} = T
_inner_type(l::SimpleSDMPredictor{T}) where {T <: Any} = T

"""
    mask!(l1::T1, l2::T2) where {T1 <: SimpleSDMLayer, T2 <: SimpleSDMLayer}

Changes the second layer so that the positions for which the first layer is zero
(of the appropriate type) or `nothing` are set to `nothing`. This is mostly
useful in cases where you have a `Bool` layer.
"""
function mask!(l1::T1, l2::T2) where {T1 <: SimpleSDMLayer, T2 <: SimpleSDMLayer}
    _itype = _inner_type(l1)
    dropfunc = (x) -> isnothing(x) || (x == zero(_itype))
    todrop = findall(dropfunc, l1.grid)
    l2.grid[todrop] .= nothing
end

"""
    mask(l1::T1, l2::T2) where {T1 <: SimpleSDMLayer, T2 <: SimpleSDMLayer}

Returns a copy of the second layer in which the positions for which the first
layer is zero (of the appropriate type) or `nothing` are set to `nothing`. This
is mostly useful in cases where you have a `Bool` layer.
"""
function mask(l1::T1, l2::T2) where {T1 <: SimpleSDMLayer, T2 <: SimpleSDMLayer}
    l3 = copy(l2)
    mask!(l1, l3)
    return l3
end


text:
# Stoyan & Stoyan, 1994, 15.11

"""
    Ktrans(xy, r, window)

Return K-function estimate with translation edge correction.
"""
function Ktrans(pp, r)
  n = npoints(pp)
  cumsum(disthisttrans(pp.data, step(r), length(r), window(pp))) / (n*(n-1)/area(window(pp)))
end

"""
   Kmark(xy, r, window, f)

Return mark-weighted K-function with translation edge correction.

This is not properly normalised.
"""
function Kmark(pp, r, f)
  n = npoints(pp)
  cumsum(disthisttransf(pp.data, f, step(r), length(r), window(pp))) / (n*(n-1)/area(window(pp)))
end

function K12(x, y, r)
  @assert window(x) == window(y)
  n1 = npoints(x)
  n2 = npoints(y)
  cumsum(disthisttrans12(x.data, y.data, step(r), length(r), window(x))) / (n1*n2/area(window(x)))
end
function L12(x, y, r)
  sqrt.(K12(x, y, r) ./ pi)
end
Lest(pp, r) = sqrt.(Kest(pp, r) ./ pi)


text:
using Coverage
a = analyze_malloc(".")
map(println, a[end-9:end]);
getbytes(x) = x.bytes
println(sum(getbytes.(a)))


text:
# -----------------------------------------------
# Types
# -----------------------------------------------

# using Optim

import Base.show

# Import types from Optim.jl
const AbsOptim = Optim.AbstractOptimizer
const AbsRes   = Optim.OptimizationResults

# #########   REGRESSION TYPES   #########

"""
## Regression Types

Every regression has a specific type inherited by the abstract type `AbstractRegressionType`.
"""
abstract type AbstractRegressionType end

# ---- LINEAR REGRESSION ----

"""
## Linear Regression

The abstract `LinearRegressionType` is the supertype for all *linear* regressions.
"""
abstract type LinearRegressionType<:AbstractRegressionType end

"""
### OLS (Ordinary Least Squares) Regression

Type that defines the simple OLS regression.
"""
struct OLS<:LinearRegressionType
    intercept::Bool
    robust::Any

    function OLS(; intercept::Bool=true, robust=false)
        if robust != false && !(typeof(robust) <: AbstractVarianceMatrix)
            error("Robust option is unrecognized")
        end
        new(intercept, robust)
    end
end

"""
### GLS (Generalized Least Squares) Regression

Type that defines the GLS estimator. Since this is a theoretical estimator, the matrix Ω must be provided.
"""
struct GLS<:LinearRegressionType
    omega
    intercept::Bool
    robust::Any

    function GLS{T<:Number}(omega::Array{T}; intercept::Bool=true, robust=false)
        if robust != false && !(typeof(robust) <: AbstractVarianceMatrix)
            error("Robust option is unrecognized")
        end
        new(omega, intercept, robust)
    end
end

"""
### FGLS (Feasible Generalized Least Squares) Regression

Type that defines the FGLS estimator: Ω is estimated by means of the residuals from an OLS regression.
"""
struct FGLS<:LinearRegressionType
    intercept::Bool
    robust::Any

    function FGLS(; intercept::Bool=true, robust=false)
        if robust != false && !(typeof(robust) <: AbstractVarianceMatrix)
            error("Robust option is unrecognized")
        end
        new(intercept, robust)
    end
end

"""
### Iterated FGLS Regression

Type that defines the *Iterated* FGLS estimator: in this case, the process of estimation for Ω through residuals is iterated until convergence.
"""
struct IteratedFGLS<:LinearRegressionType
    iterations::Int
    intercept::Bool
    robust::Any

    function IteratedFGLS(iterations::Int=1; intercept::Bool=true, robust=false)
        if robust != false && !(typeof(robust) <: AbstractVarianceMatrix)
            error("Robust option is unrecognized")
        end
        if iterations <= 0
            error("Number of iterations must be => 1")
        end
        new(iterations, intercept, robust)
    end
end

# ---- GENERAL LINEAR MODEL REGRESSION ----

abstract type GLMRegressionType<:AbstractRegressionType end

struct Logit<:GLMRegressionType
    intercept::Bool
    robust::Any

    function Logit(; intercept::Bool=true, robust=false)
        if robust != false && !(typeof(robust) <: AbstractVarianceMatrix)
            error("Robust option is unrecognized")
        end
        new(intercept, robust)
    end
end

# ---- NONLINEAR REGRESSION ----

"""
## Non-Linear Regression

The abstract `NonLinearRegressionType` is the supertype for all *non-linear* regressions.
"""
abstract type NonLinearRegressionType<:AbstractRegressionType end

"""
### 'Optimize' Regression

Type that defines a regression where the relationship between response and regressors is a function. In this case, the package *Optim.jl* is used to carry out the minimization of the sum of squared residuals.
"""
struct Optimize<:NonLinearRegressionType
    method::AbsOptim
    modelfun
    beta0
    intercept::Bool
    robust::Any

    function Optimize{T<:Number}(modelfun, beta0::Array{T}, intercept::Bool=true, robust=false; method::AbsOptim=NewtonTrustRegion())
        if robust != false && !(typeof(robust) <: AbstractVarianceMatrix)
            error("Robust option is unrecognized")
        end
        new(method, modelfun, beta0, intercept, robust)
    end
end

# #########   FITTING/RESULT TYPES   #########

"""
## Fitting Types

Every regression produces a collection of results/estimates.
Every output is thus packaged in a precise type, depending on the originating regression type, and inherits from the abstract type `AbstractFittingResult`.
"""
abstract type AbstractFittingResult end

"""
### Generic Fitting Result

Type that implements a neutral fitting result.
Should not be used in actual implementations.
"""
struct genericfitresult<:AbstractFittingResult
    callertype
    N
    beta
    residuals
    mse
    variance
    tstat
    tsig
    confint

    function genericfitresult{T<:Number}(callertype::Any, beta::Array{T}, residuals::Array{T}, variance::Array{T})
        N = length(residuals)
        mse = sum(residuals.^2)
        tstat = squeeze(beta ./ sqrt.(diag(variance)), 2)
        tsig = 2 * CDF(-abs.(tstat), (N - length(beta)), Student)
        confint = [beta - 2*sqrt.(diag(variance)), beta + 2*sqrt.(diag(variance))]
        new(callertype, N, beta, residuals, mse, variance, tstat, tsig, confint)
    end
end

"""
### Linear Fitting Result

Type that wraps results for outputs of a `LinearRegressionType` regression.
"""
struct linearfitresult<:AbstractFittingResult
    callertype
    N
    beta
    residuals
    mse
    variance
    tstat
    tsig
    confint
    robust::Any
    iterations::Int

    function linearfitresult{T<:Number}(callertype::Any, X::Array{T}, Y::Array{T}, beta::Array{T}, residuals::Array{T}, variance::Array{T}, robust::Any, iterations::Int=0)
        N = length(residuals)
        mse = sum(residuals.^2)
        tstat = squeeze(beta ./ sqrt.(diag(variance)), 2)
        tsig = 2 * CDF(-abs.(tstat), (N - length(beta)), Student)
        confint = [beta - 2*sqrt.(diag(variance)), beta + 2*sqrt.(diag(variance))]
        if robust != false && !(typeof(robust) <: AbstractVarianceMatrix)
            error("Robust option is unrecognized")
        end
        new(callertype, N, beta, residuals, mse, variance, tstat, tsig, confint, robust, iterations)
    end
end

"""
### GLM Fitting Result

Type that wraps results for outputs of a `GLMRegressionType` regression.
"""

struct glmfitresult<:AbstractFittingResult
    callertype
    N
    beta
    residuals
    fit
    variance
    tstat
    tsig
    confint
    robust::Any
    converged::Bool
    iterations::Int

    function glmfitresult{T<:Number}(callertype::Logit, X::Array{T}, Y::Array{T}, beta::Array{T}, residuals::Array{T}, variance::Array{T}, robust::Any, converged::Bool, iterations::Int=0)
        N = length(residuals)
        fit = -2*logLike(X, Y, beta, Logistic)
        tstat = beta ./ sqrt.(diag(variance))
        if length(size(tstat)) > 1
            tstat = squeeze(tstat, 2)
        end
        tsig = 2 * CDF(-abs.(tstat), (N - length(beta)), Student)
        confint = [beta - 2*sqrt.(diag(variance)), beta + 2*sqrt.(diag(variance))]
        if robust != false && !(typeof(robust) <: AbstractVarianceMatrix)
            error("Robust option is unrecognized")
        end
        new(callertype, N, beta, residuals, fit, variance, tstat, tsig, confint, robust, converged, iterations)
    end
end

"""
### Non-Linear Fitting Result

Type that wraps results for outputs of a `NonLinearRegressionType` regression.
"""
struct nlinfitresult<:AbstractFittingResult
    callertype
    N
    beta
    residuals
    mse
    variance
    tstat
    tsig
    confint
    robust::Any
    converged::Bool
    iterations::Int

    function nlinfitresult{T<:Number}(callertype::Any, X::Array{T}, Y::Array{T}, beta::Array{T}, residuals::Array{T}, variance::Array{T}, robust::Any, converged::Bool, iterations::Int=0)
        """
        beta = Optim.minimizer(optimresult)
        residuals = (Y - modelfun(X, beta))
        mse = Optim.minimum(optimresult)
        converged = Optim.converged(optimresult)
        """
        N = length(residuals)
        mse = sum(residuals.^2)
        tstat = beta ./ sqrt.(diag(variance))
        tsig = 2 * CDF(-abs.(tstat), (N - length(beta)), Student)
        confint = [beta - 2*sqrt.(diag(variance)), beta + 2*sqrt.(diag(variance))]
        if robust != false && !(typeof(robust) <: AbstractVarianceMatrix)
            error("Robust option is unrecognized")
        end
        new(callertype, N, beta, residuals, mse, variance, tstat, tsig, confint, robust, converged, iterations)
    end
end

# #########   ESTIMATION TYPES   #########

"""
## Estimation Types

Additional types that help with the estimation process of the regressions.

### Variance Matrices

Variance-covariance matrices are key elements in regression analysis, and require care.
Different types of *standard* and *robust* variance matrices are provided, and can be used standalone, if required.
All variance matrices types inherit from the abstract type `AbstractVarianceMatrix`.
"""
abstract type AbstractVarianceMatrix end

"""
### Basic-Naive Variance Matrix

Type that implements the naive variance-covariance matrix.
"""
struct BasicVariance<:AbstractVarianceMatrix

    function BasicVariance()
        new()
    end

    function BasicVariance{T<:Number}(X::Array{T}, residuals::Array{T}, k::Number)
        s2 = sum(residuals.^2) / (length(residuals) - k)
        return s2 .* inv(X' * X)
    end
end

"""
### HCE (Eicker-White) Heteroskedastic Consistent Variance Matrix

Type that implements the HCE variance-covariance matrix à la Eicker-White.
"""
struct HCEVariance<:AbstractVarianceMatrix

    function HCEVariance()
        new()
    end

    function HCEVariance{T<:Number}(X::Array{T}, residuals::Array{T})
        return inv(X' * X) * (X' * diagm(squeeze(residuals, 2).^2) * X) * inv(X' * X)
    end
end


text:
module Minibatch_GD
    using LoopVectorization, Polyester

    function fit(;model::Any, input_data::Array{Float32}, output_data::Array{Float32}, loss_function::Any, monitor::Any, α::Float64=0.01, epochs::Int64=20, batch::Real=32, mini_batch::Int64=5)
        model.initialize(model, mini_batch)
        input_shape = size(input_data)[1:end-1]
        output_shape = size(output_data)[1:end-1]
        examples = size(input_data)[end]
        batch_size = ceil(Int64, examples/batch)*mini_batch

        current_input_data = zeros(Float32, input_shape..., mini_batch)
        current_output_data = zeros(Float32, output_shape..., mini_batch)

        for e in 1:epochs
            print("Epoch ", e, "\n[")
            @time begin
                for t in 1:mini_batch:batch_size-mini_batch+1
                    @batch for i in 1:mini_batch
                        index = rand(1:examples)
                        selectdim(current_input_data, length(input_shape)+1, i) .= selectdim(input_data, length(input_shape)+1, index)
                        selectdim(current_output_data, length(output_shape)+1, i) .= selectdim(output_data, length(output_shape)+1, index)
                    end

                    if ((t+mini_batch-1)÷mini_batch)%ceil(batch_size/(50*mini_batch))==0
                        print("=")
                    end
                    model.update(model, current_input_data, current_output_data, loss_function, monitor, "Minibatch_GD", α)
                end
                print("] with loss ", model.loss, "\nTime usage: ")
                model.loss = 0.0
            end
        end
    end
end


text:
"""
This function creates a PowerSystems ThermalStandard unit.
"""
function create_PSY_generator(gen::ThermalGenEMIS{<: BuildPhase}, sys::PSY.System)

    tech = get_tech(gen)
    base_power = get_maxcap(gen)

    buses = PSY.get_components(PSY.Bus, sys)
    project_bus = PSY.Bus[]

    for b in buses
        if PSY.get_number(b) == get_bus(tech)
            push!(project_bus, b)
        end
    end

    PSY_gen =  PSY.ThermalStandard(
        get_name(gen), # name
        true,   # available
        true,   # status
        project_bus[1], # bus
        get_maxcap(gen) / base_power,    # active power
        1.0,                # reactive power
        get_maxcap(gen) / base_power,    # rating
        (min = get_mincap(gen) / base_power, max = get_maxcap(gen) / base_power), # active power limits
        nothing,      # reactive power limits
        (up = get_ramp_limits(tech)[:up] / (base_power * 60), down = get_ramp_limits(tech)[:down] / (base_power * 60)), # ramp limits
        get_operation_cost(tech), # operation cost
        base_power, # base power
        get_time_limits(tech), # up and down time limits
        PSY.PrimeMovers.Symbol(get_type(tech)), # primemover
        PSY.ThermalFuels.Symbol(get_fuel(tech)), # fuel type
    )
    return PSY_gen
end

"""
This function creates a PowerSystems RenewableDispatch unit.
"""
function create_PSY_generator(gen::RenewableGenEMIS{<: BuildPhase}, sys::PSY.System)
    tech = get_tech(gen)
    base_power = get_maxcap(gen)

    if get_type(tech) == "WT"
        primemover = PSY.PrimeMovers.WT
    elseif get_type(tech) == "PVe"
        primemover = PSY.PrimeMovers.PVe
    end

    buses = PSY.get_components(PSY.Bus, sys)
    project_bus = PSY.Bus[]

    for b in buses
        if PSY.get_number(b) == get_bus(tech)
            push!(project_bus, b)
        end
    end

    PSY_gen =  PSY.RenewableDispatch(
        get_name(gen),  # name
        true,           # available
        project_bus[1], # bus
        get_maxcap(gen) / base_power, # active power
        1.0,             # reactive power
        get_maxcap(gen) / base_power,    # rating
        primemover,     # primemover
        nothing,        # reactivepower limits
        1.0,            # power factor
        get_operation_cost(tech),
        base_power, # base power
    )
    return PSY_gen
end

"""
This function creates a PowerSystems GenericBattery unit.
"""
function create_PSY_generator(gen::BatteryEMIS{<: BuildPhase}, sys::PSY.System)
    tech = get_tech(gen)
    base_power = get_maxcap(gen)

    buses = PSY.get_components(PSY.Bus, sys)
    project_bus = PSY.Bus[]

    for b in buses
        if PSY.get_number(b) == get_bus(tech)
            push!(project_bus, b)
        end
    end

    PSY_gen =  PSY.GenericBattery(
        get_name(gen),  # name
        true,           # available
        project_bus[1], # bus
        PSY.PrimeMovers.Symbol(get_type(tech)), # primemover
        get_soc(tech) / (get_maxcap(gen) * base_power), #initial state of charge
        (min = get_storage_capacity(tech)[:min] / (get_maxcap(gen) * base_power), max = get_storage_capacity(tech)[:max] / (get_maxcap(gen) * base_power)), # state of charge limits
        get_maxcap(gen) / base_power, # rating
        get_maxcap(gen) / base_power, # active power
        (min = get_input_active_power_limits(tech)[:min] / base_power, max = get_input_active_power_limits(tech)[:max] / base_power), # input active power limits
        (min = get_output_active_power_limits(tech)[:min] / base_power, max = get_output_active_power_limits(tech)[:max] / base_power), # output active power limits
        get_efficiency(tech), # efficiency
        1.0,             # reactive power
        nothing,      # reactive power limits
        base_power, # base power
    )
    return PSY_gen
end


text:
#
# Define observations, or "truth data".
#
# For observational time-series, we specify
# the time-range that the observations are defined, which may be different from the
# data contained in a file.
#
# If the observations are one-dimensional, then forward maps are constructed that use
# a single 3D Oceananigans model to evaluate an ensemble of column models.
#
# If the observations are two-dimensional, then forward maps are constructed that use
# a single 3D Oceananigans model to evaluate an ensemble of "slice" 2D models.
#

free_convection_path = "/path/to/data"
free_convection_observation = OneDimensionalTimeSeries(free_convection_path, normalization=Variance(), field_names=(:u, :v, :b, :e), time_range=range(2hours, stop=2days, step=4hours))

# normalize(::Variance, field) =

# Some other possibilities

# 1. Construct a batch of observations:
# wind_stress_observation = OneDimensionalTimeSeries(wind_stress_path, time_range=range(2hours, stop=2days, step=4hours))

file_paths = [free_convection_path, wind_stress_path]
observations = [OneDimensionalTimeSeries(path, field_names=(:u, :v, :b, :e)) for path in file_paths]

# 2. Construct two-dimensional observations at a single snapshot:
# observation = TwoDimensionalSnapshot(data_path, iteration=LastIteration())

# EKI solves
#
# y = G(θ) + η
#
# where G(θ) is the forward map and y is the observations.

@free_parameters StabilityFnParameters CᴷRiʷ CᴷRiᶜ Cᴷu⁻ Cᴷuʳ Cᴷc⁻ Cᴷcʳ Cᴷe⁻ Cᴷeʳ
@free_parameters StabilityFnParameters CᴷRiʷ CᴷRiᶜ Cᴷuʳ Cᴷc⁻ Cᴷcʳ Cᴷe⁻ Cᴷeʳ

stability_fn_parameters_priors = StabilityFnParameters(
    CᴷRiʷ = Normal(0.25, 0.05) |> logify,
    CᴷRiᶜ = Normal(0.25, 0.5) |> logify,
    Cᴷu⁻ = Uniform(0.0, 10.0),
    Cᴷuʳ = Uniform(0.0, 10.0),
    Cᴷc⁻ = Uniform(0.0, 10.0),
    Cᴷcʳ = Uniform(0.0, 10.0),
    Cᴷe⁻ = Uniform(0.0, 10.0),
    Cᴷeʳ = Uniform(0.0, 10.0)
)

free_parameters = FreeParameters(parameters = StabilityFnParameters,
                                 prior = stability_fn_parameters_priors)

free_parameters = FreeParameters(parameters = BoundaryLayerMesoscaleCombinedParameters)

# Tricky question: where do we build the model?
#
# 1. Build the model in ForwardMap, dispatching on the type of observations
# 2. Just ask users to "pull up their pants" and build the model in their scripts.

# Define the forward map.
#
# With forward_map = ForwardMap(), G(θ) outputs a concatenated vector of all fields
# at all times.
#
# With forward_map = TimeAveragedForwardMap(If forwardloss_function = TimeAveragedLossFunction() (for example), then G(θ) outputs a _scalar_,
# representing some time-averaged measurement of discrepency between model and data.
#
# Questions:
#
# * How do we handle NaNs

struct InverseProblem
    observations
    simulation
    transformation
    parameters
end

simple_calibration_problem = InverseProblem(simulation, observation, transformation=nothing)
time_averaged_calibration_problem = InverseProblem(simulation, observation, transformation=AccumulatedTimeAverage())

calibrate!(calibration_problem)


text:
using QCMaterial
using PyCall
using LinearAlgebra
using HDF5
import Random
import PyCall: pyimport

#make taus list file.h5
#length=8
function taus_list(file_name)
    x = Float64[]
    for i in 1:10
        tau = 0.005*(i-1) + exp((0.1*(i-1))) -1
        push!(x,tau)
    end
    fid = h5open(file_name,"w")
    fid["/test/data"] = x
    close(fid)
end

#read taus list file.h5
function read_taus_list(file_name)
    fid = h5open(file_name,"r")
    x = fid["/test/data"][:]
    close(fid)
    return x
end

#read taus list file.txt
function read_and_parse_float(file_name)
    x = Float64[]
    open(file_name, "r") do fp
        num_elm = parse(Int64,readline(fp))
        for i in 1:num_elm
            push!(x, parse(Float64, readline(fp)))
        end
    end
    return x
end

function generate_impurity_ham_with_1imp_3bath_dmft(U::Float64, μ::Float64,  nsite::Integer)
    spin_index_functions = [up_index, down_index]
    so_idx(iorb, ispin) = spin_index_functions[ispin](iorb)

    ham_op1 = FermionOperator()
    ham_op1 += FermionOperator("$(up_index(1))^ $(down_index(1))^ $(up_index(1)) $(down_index(1))", -U)

    #chemical potential
    for ispin in [1, 2]
        ham_op1 += FermionOperator("$(so_idx(1, ispin))^ $(so_idx(1, ispin))", -μ)
    end

    #bath energy level
    ε = [0.0, 1.11919, 0.00000, -1.11919]
    for ispin in [1, 2]
        for i in 2:nsite
            ham_op1 += FermionOperator("$(so_idx(i, ispin))^ $(so_idx(i, ispin))", ε[i])
        end
    end

    #hybridization
    V = [0.0, -1.26264, 0.07702, -1.26264]

    for ispin in [1, 2]
        for i in 2:nsite
            ham_op1 += FermionOperator("$(so_idx(1, ispin))^ $(so_idx(i, ispin))", V[i])
            ham_op1 += FermionOperator("$(so_idx(i, ispin))^ $(so_idx(1, ispin))", V[i])
        end
    end
    ham_op1
end

function hev(n_qubit, depth)
    circuit = QulacsParametricQuantumCircuit(n_qubit)
    for d in 1:depth
        for i in 1:n_qubit
            add_parametric_RY_gate!(circuit, i, 0.0)
            add_parametric_RZ_gate!(circuit, i, 0.0)
        end
        for i in 1:n_qubit÷2
            add_CNOT_gate!(circuit, 2*i-1, 2*i)
        end
        for i in 1:n_qubit÷2-1
            add_CNOT_gate!(circuit, 2*i, 2*i+1)
        end
    end
    for i in 1:n_qubit
        add_parametric_RY_gate!(circuit, i, 0.0)
        add_parametric_RZ_gate!(circuit, i, 0.0)
    end
    circuit = QulacsVariationalQuantumCircuit(circuit)
    circuit
end

#beta = 1000 (T=0.001)
nsite = 4
n_qubit　= 2*nsite
U = 4.0
μ = U/2
d_theta = 1e-5
verbose = QCMaterial.MPI_rank == 0
Random.seed!(100)

taus = read_and_parse_float(ARGS[3])
println("taus=",taus)

#Hamiltonian
ham_op1 = generate_impurity_ham_with_1imp_3bath_dmft(U, μ, nsite)
n_electron_gs = 4
@assert mod(n_electron_gs, 2) == 0
sparse_mat = get_number_preserving_sparse_operator(ham_op1, n_qubit, n_electron_gs);
enes_ed = eigvals(sparse_mat.toarray());

#N+1 system
n_electron_ex = 5
sparse_mat_ex = get_number_preserving_sparse_operator(ham_op1, n_qubit, n_electron_ex);
enes_ed_ex = eigvals(sparse_mat_ex.toarray());
#debug
println("Ground energy_ED=",minimum(enes_ed))
println("Ground energy_N+1_ED=",minimum(enes_ed_ex))

#1 system
n_electron_1 = 1
sparse_mat_enum_1 = get_number_preserving_sparse_operator(ham_op1, n_qubit, n_electron_1);
enes_ed_enum_1 = eigvals(sparse_mat_enum_1.toarray());
#debug
#println("Ground energy_ED=",minimum(enes_ed_enum_1))
println("Ground energy_1_ED=",minimum(enes_ed_enum_1))

#2 system
n_electron_2 = 2
sparse_mat_enum_2 = get_number_preserving_sparse_operator(ham_op1, n_qubit, n_electron_2);
enes_ed_enum_2 = eigvals(sparse_mat_enum_2.toarray());
#debug
#println("Ground energy_ED=",minimum(enes_ed_enum_2))
println("Ground energy_2_ED=",minimum(enes_ed_enum_2))

#3 system
n_electron_3 = 3
sparse_mat_enum_3 = get_number_preserving_sparse_operator(ham_op1, n_qubit, n_electron_3);
enes_ed_enum_3 = eigvals(sparse_mat_enum_3.toarray());
#debug
#println("Ground energy_ED=",minimum(enes_ed_enum_2))
println("Ground energy_3_ED=",minimum(enes_ed_enum_3))

#ansatz
state0 = create_hf_state(n_qubit, n_electron_gs)
vc = uccgsd(n_qubit, nocc = 2, orbital_rot=true, uccgsd = true, p_uccgsd = false)
#depth = n_qubit*2
#vc = hev(n_qubit, depth)
theta_init = rand(num_theta(vc))

#Perform VQE
cost_history, thetas_opt =
QCMaterial.solve_gs(jordan_wigner(ham_op1), vc, state0, theta_init=theta_init, verbose=true,
    comm=QCMaterial.MPI_COMM_WORLD
)

#debug
println("Ground energy_VQE=",cost_history[end])

#c^{dag},c
up1 = up_index(1)
down1 = down_index(1)
up2 = up_index(2)
down2 = down_index(2)

if ARGS[1] == "plus_true"
    right_op = FermionOperator("$(up1)^ ", 1.0)
    right_op = jordan_wigner(right_op)
    left_op = FermionOperator("$(up1) ", 1.0)
    left_op = jordan_wigner(left_op)

    n_electron_ex = 5

    sign = 1
end

if ARGS[2] == "minus_true"
    right_op = FermionOperator("$(up1) ", 1.0)
    right_op = jordan_wigner(right_op)
    left_op = FermionOperator("$(up1)^ ", 1.0)
    left_op = jordan_wigner(left_op)

    n_electron_ex = 3

    sign = -1
end

vc_ex = uccgsd(n_qubit, nocc = 2, orbital_rot=true, uccgsd = true, p_uccgsd = false)
#vc_ex = hev(n_qubit, depth)
#println(vc_ex)
state_gs = create_hf_state(n_qubit, n_electron_gs)
update_circuit_param!(vc, thetas_opt)
update_quantum_state!(vc, state_gs)

state0_ex = create_hf_state(n_qubit, n_electron_ex)

#taus = collect(range(0.0, 0.02, length = 2))

#make taus list file
#taus_list("dimer_plus.h5")
#generate taus list file
#taus = read_taus_list("dimer_plus.h5")
#println("taus=",taus)

Gfunc_ij_list, norm =  compute_gtau(
    jordan_wigner(ham_op1),
    left_op,
    right_op,
    vc_ex,
    state_gs,
    state0_ex,
    taus,
    d_theta,
    verbose = verbose,
    algorithm = "vqs",
    recursive= true
)

Gfunc_ij_list *= sign

println("Gfunc_ij_list_plus=", Gfunc_ij_list)
println("norm=", norm)

function write_to_txt(file_name, x, y)
    open(file_name, "w") do fp
        for i = 1:length(x)
            println(fp, x[i], " ", real(y[i]))
        end
    end
end

function write_to_txt_(file_name, x, y)
    open(file_name, "w") do fp
        for i = 1:length(x)
            println(fp, 1000-x[i], " ", -real(y[i]))
        end
    end
end
write_to_txt("gf_4site_plus_relative_noise_100000.txt", taus, Gfunc_ij_list)
write_to_txt("norm_4site_plus_relative_noise_100000", taus, norm)
#write_to_txt_("gf_4site_minus_recursive_direct_doubles_.txt", taus, Gfunc_ij_list)
println("done!")

max_stars_repo_name: sakurairihito/QCMaterialNew

max_issues_repo_name: sakurairihito/QCMaterialNew

max_forks_repo_name: sakurairihito/QCMaterialNew


text:
##############################################################################
##
## References
##
##############################################################################

# - https://en.bitcoin.it/wiki/WIF

##############################################################################
##
## Wallet interchange format function definitions
##
##############################################################################

# Convert private key to WIF.
# TODO: turn keys into objects to hold metadata
# TODO: Update to use byte arrays
# - network_id: which network to use; 0x80 for mainnet, 0xef for testnet
# - compression: 01 if private key corresponds to compressed public key
function private2wif(private_key; network_id = "80", compression = "")
  private_key = string(network_id, private_key, compression)

  hashed = Crypto.digest("SHA256", private_key, is_hex=true)
  hashed = Crypto.digest("SHA256", hashed)

  checksum = Crypto.oct2hex(hashed[1:4])

  private_key = string(private_key, checksum)

  return encode58(parseint(BigInt, private_key, 16))
end

function wif2private(wif; is_compressed=false)
  private_key = hex(decode58(wif))

  # Drop the first two characters (network identifier) and last
  # 8 (checksum)
  private_key = private_key[3:end-8]

  if is_compressed
    private_key = private_key[1:end-2]
  end

  return private_key
end

function wif_check_sum(wif)
  result = hex(decode58(wif))
  return get_checksum(result[1:end-8], is_hex=true) == result[end-8+1:end]
end

function pub2base58(pub_key::String; network_id = "00")
  pub_key_length = div(length(pub_key), 2)

  # If public key is elliptic curve coordinate, hash with SHA-256
  if pub_key_length == 65
    pub_key = Crypto.digest("SHA256", pub_key, is_hex = true)
    pub_key = Crypto.oct2hex(pub_key)
    pub_key_length = div(length(pub_key), 2)
  end

  # If public key has been SHA-256 hashed, hash with RIPEMD-160
  if pub_key_length == 32
    pub_key = Crypto.digest("RIPEMD160", pub_key, is_hex = true)
    pub_key = Crypto.oct2hex(pub_key)
    pub_key_length = div(length(pub_key), 2)
  end

  # If public key has been RIPEMD-160 hashed, add network id
  if pub_key_length == 20
    pub_key = string(network_id, pub_key)
    pub_key_length = div(length(pub_key), 2)
  end

  # If public key has network id added, add checksum
  if pub_key_length == 21
    checksum = get_checksum(pub_key, is_hex = true)
    pub_key = string(pub_key, checksum)
    pub_key_length = div(length(pub_key), 2)
  end

  # If public key has checksum added
  if pub_key_length == 25
    pub_key = encode58(Crypto.hex2oct(pub_key))
  end

  return pub_key
end


text:
using Test
using MPI

@test !MPI.Initialized()
MPI.Init()
@test MPI.Initialized()

@test 0 <= MPI.Comm_rank(MPI.COMM_WORLD) < MPI.Comm_size(MPI.COMM_WORLD)

@test !MPI.Finalized()
MPI.Finalize()
@test MPI.Finalized()


text:
__precompile__()

module NODAL
    using Distributed, Statistics, Random, Printf

    # New Methods for Base Functions
    import Base: convert,
                 show,
                 getindex,
                 setindex!,
                 put!,
                 take!,
                 fetch,
                 isready

    # Types
    include("core/parameters.jl")
    include("core/configuration.jl")
    include("core/run.jl")

    # Methods
    include("core/perturb!.jl")
    include("core/perturb_elements!.jl")
    include("core/update!.jl")
    include("core/neighbor!.jl")
    include("core/unit_value.jl")

    # Measurement
    include("core/measurement/measure.jl")
    include("core/measurement/results.jl")

    # Search Blocks
    include("core/search/blocks/first_improvement.jl")
    include("core/search/blocks/probabilistic_improvement.jl")
    include("core/search/blocks/greedy_construction.jl")
    include("core/search/blocks/random_walk.jl")

    # Techniques
    include("core/search/techniques/technique.jl")
    include("core/search/techniques/iterative_first_improvement.jl")
    include("core/search/techniques/iterative_probabilistic_improvement.jl")
    include("core/search/techniques/randomized_first_improvement.jl")
    include("core/search/techniques/iterative_greedy_construction.jl")
    include("core/search/techniques/simulated_annealing.jl")
    include("core/search/techniques/iterated_local_search.jl")

    # Tools
    include("core/search/tools/initialize_search_tasks!.jl")
    include("core/search/tools/end_search_tasks!.jl")
    include("core/search/tools/get_new_best.jl")
    include("core/search/tools/stopping_criteria.jl")
    include("core/search/tools/reporting_criteria.jl")
    include("core/search/tools/temperatures.jl")

    # Optimize
    include("core/search/optimize.jl")

    # Utilities
    include("util/random.jl")
    include("util/show.jl")
    include("util/resultchannel.jl")

    # Types
    export Parameter,
           NumberParameter,
           IntegerParameter,
           FloatParameter,
           PermutationParameter,
           EnumParameter,
           StringParameter,
           BoolParameter,
           Configuration,
           Result,
           ResultChannel,
           Run

    # Methods
    export perturb!,
           perturb_elements!,
           neighbor!,
           optimize,
           test_names,
           test_A,
           update!,
           unit_value,
           unit_value!

    # Measurement Tools
    export measure_mean!,
           sequential_measure_mean!

    # Search Building Blocks
    export first_improvement,
           probabilistic_improvement,
           greedy_construction,
           random_walk

    # Search Techniques
    export technique,
           simulated_annealing,
           iterative_gredy_construction,
           iterative_first_improvement,
           randomized_first_improvement,
           iterative_probabilistic_improvement,
           iterated_local_search

    # Search Tools
    export initialize_search_tasks!,
           end_search_tasks!,
           get_new_best,
           elapsed_time_criterion,
           iterations_criterion,
           iterations_reporting_criterion,
           elapsed_time_reporting_criterion,
           log_temperature
end


text:
module BaytesMCMC

################################################################################
#Import modules
using BaytesCore:
    BaytesCore,
    AbstractAlgorithm,
    AbstractTune,
    AbstractConfiguration,
    AbstractDiagnostics,
    AbstractKernel,
    AbstractConstructor,
    Updater,
    Iterator,
    UpdateBool,
    UpdateTrue,
    UpdateFalse,
    logsumexp,
    logaddexp,
    logmeanexp,
    update,
    AcceptStatistic,
    SampleDefault,
    to_Tuple

import BaytesCore:
    BaytesCore,
    update,
    update!,
    infer,
    results,
    init,
    init!,
    propose,
    propose!,
    propagate,
    propagate!,
    result!,
    get_result,
    generate_showvalues

using ModelWrappers:
    ModelWrappers,
    ModelWrapper,
    Tagged,
    Objective,
    DiffObjective,
    AbstractDifferentiableTune,
    ℓObjectiveResult,
    ℓDensityResult,
    ℓGradientResult,
    checkfinite,
    AutomaticDiffTune,
    AnalyticalDiffTune,
    sample,
    sample!

import ModelWrappers: ModelWrappers, predict, generate, checkfinite

using Random: Random, AbstractRNG, GLOBAL_RNG, randexp
using LinearAlgebra:
    LinearAlgebra, Symmetric, Diagonal, diag, cholesky, inv, UniformScaling, dot, det, I
using Statistics: Statistics, mean, median, std, sqrt, quantile, var, cor, cov, middle
using Distributions: Distributions, logpdf, MvNormal

using DocStringExtensions:
    DocStringExtensions, TYPEDEF, TYPEDFIELDS, FIELDS, SIGNATURES, FUNCTIONNAME
using ArgCheck: ArgCheck, @argcheck
using UnPack: UnPack, @unpack, @pack!

################################################################################
#Abstract types to be dispatched in Examples section
abstract type MCMCKernel <: AbstractKernel end
abstract type MCMCKernelDiagnostics <: AbstractDiagnostics end
function get_acceptrate() end
function move() end

include("Core/Core.jl")
include("Kernels/Kernels.jl")

################################################################################
export
    # BaytesCore
    UpdateBool,
    UpdateTrue,
    UpdateFalse,
    propose,
    propose!,
    propagate,
    propagate!,
    results,
    update!,
    SampleDefault,

    # ModelWrappers
    checkfinite,
    init,
    init!,
    predict,
    generate,

    #MCMC
    MCMCKernel,
    MCMCKernelDiagnostics,
    MCMCDiagnostics,
    get_acceptrate,
    move

end


text:
"""
    splitup(data::Vector{UInt8}, offsets, T::Type; skipbytes=0)

Given the `offsets` and `data` return by `array(...; raw = true)`, reconstructed the actual
array (with custome struct, can be jagged as well).
"""
function splitup(data::Vector{UInt8}, offsets, T::Type; skipbytes=0, jagged=true)
    packedsize = packedsizeof(T)

    if jagged
        out = Vector{Vector{T}}()
    else
        out = Vector{T}()
    end
    sizehint!(out, length(offsets))

    io = IOBuffer(data)
    for l in diff(offsets)
        skip(io, skipbytes)
        if jagged
            n = (l - skipbytes) / packedsize
            push!(out, [readtype(io, T) for _ in 1:n])
        else
            # n != 1 && warning("The packed size of the entry does not match the data size.")
            push!(out, readtype(io, T))
        end
    end
    out
end

# Custom struct interpretation
abstract type CustomROOTStruct end

# TLorentzVector
const LVF64 = LorentzVector{Float64}
Base.show(io::IO, lv::LorentzVector) = print(io, "LV(x=$(lv.x), y=$(lv.y), z=$(lv.z), t=$(lv.t))")
function Base.reinterpret(::Type{LVF64}, v::AbstractVector{UInt8}) where T
    # first 32 bytes are TObject header we don't care
    # x,y,z,t in ROOT
    v4 = ntoh.(reinterpret(Float64, @view v[1+32:end]))
    # t,x,y,z in LorentzVectors.jl
    LVF64(v4[4], v4[1], v4[2], v4[3])
end

"""
    interped_data(rawdata, rawoffsets, ::Type{Vector{LorentzVector{Float64}}}, ::Type{Offsetjagg})

The `interped_data` method specialized for `LorentzVector`. This method will get called by
[`basketarray`](@ref) instead of the default method for `TLorentzVector` branch.
"""
function interped_data(rawdata, rawoffsets, ::Type{Vector{LVF64}}, ::Type{Offsetjagg})
    _size = 64 # needs to account for 32 bytes header
    dp = 0 # book keeping for copy_to!
    lr = length(rawoffsets)
    offset = Vector{Int32}(undef, lr)
    offset[1] = 0
    @views @inbounds for i in 1:lr-1
        start = rawoffsets[i]+10+1
        stop = rawoffsets[i+1]
        l = stop-start+1
        if l > 0
            unsafe_copyto!(rawdata, dp+1, rawdata, start, l)
            dp += l
            offset[i+1] = offset[i] + l
        else
            offset[i+1] = offset[i]
        end
    end
    resize!(rawdata, dp)
    real_data = interped_data(rawdata, offset, LVF64, Nojagg)
    offset .÷= _size
    offset .+= 1
    VectorOfVectors(real_data, offset)
end
function interped_data(rawdata, rawoffsets, ::Type{LVF64}, ::Type{J}) where {T, J <: JaggType}
    # even with rawoffsets, we know each TLV is destinied to be 64 bytes
    [
     reinterpret(LVF64, x) for x in Base.Iterators.partition(rawdata, 64)
    ]
end
# TLorentzVector ends

# KM3NeT
struct _KM3NETDAQHit <: CustomROOTStruct
    dom_id::Int32
    channel_id::UInt8
    tdc::Int32
    tot::UInt8
end
function readtype(io::IO, T::Type{_KM3NETDAQHit})
    T(readtype(io, Int32), read(io, UInt8), read(io, Int32), read(io, UInt8))
end
function interped_data(rawdata, rawoffsets, ::Type{Vector{_KM3NETDAQHit}}, ::Type{J}) where {T, J <: UnROOT.JaggType}
    UnROOT.splitup(rawdata, rawoffsets, _KM3NETDAQHit, skipbytes=10)
end

struct _KM3NETDAQTriggeredHit
    dom_id::Int32
    channel_id::UInt8
    tdc::Int32
    tot::UInt8
    trigger_mask::UInt64
end
packedsizeof(::Type{_KM3NETDAQTriggeredHit}) = 24  # incl. cnt and vers
function readtype(io::IO, T::Type{_KM3NETDAQTriggeredHit})
    dom_id = readtype(io, Int32)
    channel_id = read(io, UInt8)
    tdc = read(io, Int32)
    tot = read(io, UInt8)
    skip(io, 6)
    trigger_mask = readtype(io, UInt64)
    T(dom_id, channel_id, tdc, tot, trigger_mask)
end

function UnROOT.interped_data(rawdata, rawoffsets, ::Type{Vector{_KM3NETDAQTriggeredHit}}, ::Type{J}) where {T, J <: UnROOT.JaggType}
    UnROOT.splitup(rawdata, rawoffsets, _KM3NETDAQTriggeredHit, skipbytes=10)
end

struct _KM3NETDAQEventHeader
    detector_id::Int32
    run::Int32
    frame_index::Int32
    UTC_seconds::UInt32
    UTC_16nanosecondcycles::UInt32
    trigger_counter::UInt64
    trigger_mask::UInt64
    overlays::UInt32
end
packedsizeof(::Type{_KM3NETDAQEventHeader}) = 76

function readtype(io::IO, T::Type{_KM3NETDAQEventHeader})
    skip(io, 18)
    detector_id = readtype(io, Int32)
    run = readtype(io, Int32)
    frame_index = readtype(io, Int32)
    skip(io, 6)
    UTC_seconds = readtype(io, UInt32)
    UTC_16nanosecondcycles = readtype(io, UInt32)
    skip(io, 6)
    trigger_counter = readtype(io, UInt64)
    skip(io, 6)
    trigger_mask = readtype(io, UInt64)
    overlays = readtype(io, UInt32)
    T(detector_id, run, frame_index, UTC_seconds, UTC_16nanosecondcycles, trigger_counter, trigger_mask, overlays)
end

function UnROOT.interped_data(rawdata, rawoffsets, ::Type{_KM3NETDAQEventHeader}, ::Type{J}) where {T, J <: UnROOT.JaggType}
    UnROOT.splitup(rawdata, rawoffsets, _KM3NETDAQEventHeader, jagged=false)
end


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

## hashing a single value ##

hash(x::Any) = hash(x, zero(UInt))
hash(w::WeakRef, h::UInt) = hash(w.value, h)

## hashing general objects ##

hash(x::ANY, h::UInt) = hash_uint(3h - object_id(x))

## core data hashing functions ##

function hash_64_64(n::UInt64)
    local a::UInt64 = n
    a = ~a + a << 21
    a =  a ⊻ a >> 24
    a =  a + a << 3 + a << 8
    a =  a ⊻ a >> 14
    a =  a + a << 2 + a << 4
    a =  a ⊻ a >> 28
    a =  a + a << 31
    return a
end

function hash_64_32(n::UInt64)
    local a::UInt64 = n
    a = ~a + a << 18
    a =  a ⊻ a >> 31
    a =  a * 21
    a =  a ⊻ a >> 11
    a =  a + a << 6
    a =  a ⊻ a >> 22
    return a % UInt32
end

function hash_32_32(n::UInt32)
    local a::UInt32 = n
    a = a + 0x7ed55d16 + a << 12
    a = a ⊻ 0xc761c23c ⊻ a >> 19
    a = a + 0x165667b1 + a << 5
    a = a + 0xd3a2646c ⊻ a << 9
    a = a + 0xfd7046c5 + a << 3
    a = a ⊻ 0xb55a4f09 ⊻ a >> 16
    return a
end

if UInt === UInt64
    hash_uint64(x::UInt64) = hash_64_64(x)
    hash_uint(x::UInt)     = hash_64_64(x)
else
    hash_uint64(x::UInt64) = hash_64_32(x)
    hash_uint(x::UInt)     = hash_32_32(x)
end

## symbol & expression hashing ##

if UInt === UInt64
    hash(x::Expr, h::UInt) = hash(x.args, hash(x.head, h + 0x83c7900696d26dc6))
else
    hash(x::Expr, h::UInt) = hash(x.args, hash(x.head, h + 0x96d26dc6))
end

hash(x::QuoteNode, h::UInt) = hash(x.value, hash(QuoteNode, h))

# hashing ranges by component at worst leads to collisions for very similar ranges
const hashr_seed = UInt === UInt64 ? 0x80707b6821b70087 : 0x21b70087
function hash(r::Range, h::UInt)
    h += hashr_seed
    h = hash(first(r), h)
    h = hash(step(r), h)
    h = hash(last(r), h)
end

max_stars_repo_name: caodi/julia

max_issues_repo_name: caodi/julia

max_forks_repo_name: caodi/julia