text:
import ModelingToolkit
import ModelingToolkit: build_function, substitute
import Catalyst
import Catalyst: ReactionSystem
using StaticArrays
import Base.show
using Setfield

struct ReactionSet
    rates::Vector{Float64}
    rstoich::Vector{Vector{Pair{Int64,Int64}}}
    nstoich::Vector{Vector{Pair{Int64,Int64}}}
    nspecies::Int
end

abstract type AbstractJumpRateAggregatorAlgorithm end

struct GillespieDirect <: AbstractJumpRateAggregatorAlgorithm end

struct DepGraphDirect <: AbstractJumpRateAggregatorAlgorithm end

abstract type AbstractJumpRateAggregator end

struct DirectAggregator{U} <: AbstractJumpRateAggregator
    "the sum of all reaction rates"
    sumrate::Float64

    "a vector of the current reaction rates"
    rates::Vector{Float64}

    "sum of group propensities"
    gsums::Vector{Float64}

    "maps reaction indices to group indices"
    ridtogroup::U

    "time span for aggregation"
    tspan::Tuple{Float64,Float64}

    "time of last recorded reaction"
    tprev::Float64

    "accumulated log-probability"
    weight::Float64
end

function build_aggregator(alg::GillespieDirect, reactions::ReactionSet, ridtogroup, tspan=(0.0, Inf64))
    ngroups = maximum(ridtogroup)
    nreactions = length(reactions.rates)
    DirectAggregator(0.0, zeros(nreactions), zeros(ngroups), ridtogroup, tspan, tspan[1], 0.0)
end

struct DepGraphAggregator{U,DepGraph} <: AbstractJumpRateAggregator
    "the sum of all reaction rates"
    sumrate::Float64

    "a vector of the current reaction rates"
    rates::Vector{Float64}

    "sum of group propensities"
    gsums::Vector{Float64}

    "maps reaction indices to group indices"
    ridtogroup::U

    "time span for aggregation"
    tspan::Tuple{Float64,Float64}

    "time of last recorded reaction"
    tprev::Float64

    "accumulated log-probability"
    weight::Float64

    "dependency graph"
    depgraph::DepGraph

    "previously executed reaction"
    prev_reaction::Int
end

function build_aggregator(alg::DepGraphDirect, reactions::ReactionSet, ridtogroup, tspan=(0.0, Inf64))
    ngroups = maximum(ridtogroup)
    nreactions = length(reactions.rates)
    depgraph = make_depgraph(reactions)
    DepGraphAggregator(0.0, zeros(nreactions), zeros(ngroups), ridtogroup, tspan, tspan[1], 0.0, depgraph, 0)
end

function initialize_aggregator(agg::AbstractJumpRateAggregator; tspan=(0.0, Inf64))
    agg = @set agg.tspan = tspan
    agg = @set agg.tprev = tspan[1]
    agg = @set agg.weight = 0.0
    agg = @set agg.sumrate = 0.0
    agg.rates .= 0.0
    agg.gsums .= 0.0
    agg
end

function initialize_aggregator(agg::DepGraphAggregator; tspan=(0.0, Inf64))
    agg = @set agg.tspan = tspan
    agg = @set agg.tprev = tspan[1]
    agg = @set agg.weight = 0.0
    agg = @set agg.sumrate = 0.0
    agg = @set agg.prev_reaction = 0
    agg.rates .= 0.0
    agg.gsums .= 0.0
    agg
end

function ReactionSet(js::ModelingToolkit.JumpSystem, p)
    parammap = map(Pair, ModelingToolkit.parameters(js), p)
    statetoid = Dict(ModelingToolkit.value(state) => i for (i, state) in enumerate(ModelingToolkit.states(js)))

    rates = Float64[]
    rstoich_vec = Vector{Pair{Int64,Int64}}[]
    nstoich_vec = Vector{Pair{Int64,Int64}}[]

    for eq in ModelingToolkit.equations(js)
        rate = ModelingToolkit.value(ModelingToolkit.substitute(eq.scaled_rates, parammap))
        rstoich = sort!([statetoid[ModelingToolkit.value(spec)] => stoich for (spec, stoich) in eq.reactant_stoch])
        nstoich = sort!([statetoid[ModelingToolkit.value(spec)] => stoich for (spec, stoich) in eq.net_stoch])

        push!(rates, rate)
        push!(rstoich_vec, rstoich)
        push!(nstoich_vec, nstoich)
    end

    ReactionSet(rates, rstoich_vec, nstoich_vec, length(ModelingToolkit.states(js)))
end

function species_to_dependent_reaction_map(reactions::ReactionSet)
    nspecies = reactions.nspecies
    # map from a species to reactions that depend on it
    spec_to_dep_rxs = [Vector{Int}() for n = 1:nspecies]
    for (rx, complex) in enumerate(reactions.rstoich)
        for (spec, stoch) in complex
            push!(spec_to_dep_rxs[spec], rx)
        end
    end

    foreach(s -> unique!(sort!(s)), spec_to_dep_rxs)
    spec_to_dep_rxs
end

function make_depgraph(reactions::ReactionSet)
    nreactions = length(reactions.rates)
    spec_to_dep_rxs = species_to_dependent_reaction_map(reactions)

    # create map from rx to reactions depending on it
    dep_graph = [Vector{Int}() for n = 1:nreactions]
    for rx in 1:nreactions
        # rx changes spec, hence rxs depending on spec depend on rx
        for (spec, stoch) in reactions.nstoich[rx]
            for dependent_rx in spec_to_dep_rxs[spec]
                push!(dep_graph[rx], dependent_rx)
            end
        end
    end

    foreach(deps -> unique!(sort!(deps)), dep_graph)
    dep_graph
end

struct TrajectoryDistribution{A}
    reactions::ReactionSet
    aggregator::A
end

function TrajectoryDistribution(reactions::ReactionSet, alg::AbstractJumpRateAggregatorAlgorithm, ridtogroup=1:length(reactions.rates))
    agg = build_aggregator(alg, reactions, ridtogroup)
    TrajectoryDistribution(reactions, agg)
end

distribution(rn::ReactionSystem, p; update_map=1:Catalyst.numreactions(rn), alg=GillespieDirect()) = TrajectoryDistribution(ReactionSet(convert(ModelingToolkit.JumpSystem, rn), p), alg, update_map)

@fastmath function Distributions.logpdf(dist::TrajectoryDistribution, trajectory)::Float64
    traj_iter = trajectory_iterator(trajectory)
    tprev = 0.0
    result = 0.0
    agg = initialize_aggregator(dist.aggregator)
    for (u, t, i) in traj_iter
        agg = update_rates(agg, (u, t, i), dist.reactions)

        dt = t - tprev
        result -= dt * agg.sumrate
        if i != 0
            @inbounds gid = agg.ridtogroup[i]
            gid != 0 && @inbounds result += log(agg.gsums[gid])
        end

        tprev = t
    end

    result
end

@fastmath function fold_logpdf(dist::TrajectoryDistribution, agg::AbstractJumpRateAggregator, (u, t, i))
    agg = update_rates(agg, (u, t, i), dist.reactions)
    dt = t - agg.tprev
    log_jump_prob = 0.0
    if i != 0
        @inbounds gid = agg.ridtogroup[i]
        gid != 0 && @inbounds log_jump_prob = log(agg.gsums[gid])
    end
    log_surv_prob = -dt * agg.sumrate
    agg = @set agg.weight = agg.weight + log_surv_prob + log_jump_prob
    agg = @set agg.tprev = t
end

function step_energy(dist::TrajectoryDistribution, agg::AbstractJumpRateAggregator, (u, t, i))
    if t <= agg.tspan[1]
        return agg
    end
    if t > agg.tspan[2]
        if agg.tprev >= agg.tspan[2]
            return agg
        end
        fold_logpdf(dist, agg, (u, agg.tspan[2], 0))
    else
        fold_logpdf(dist, agg, (u, t, i))
    end
end

function trajectory_energy(dist::TrajectoryDistribution, traj; tspan=(0.0, Inf64))
    agg = initialize_aggregator(dist.aggregator, tspan=tspan)
    traj_iter = trajectory_iterator(traj)
    agg = Base.foldl((acc, x) -> step_energy(dist, acc, x), traj_iter; init=agg)
    agg.weight
end

function cumulative_logpdf!(result::AbstractVector, dist::TrajectoryDistribution, traj, dtimes::AbstractVector)
    tspan = (first(dtimes), last(dtimes))
    agg = initialize_aggregator(dist.aggregator, tspan=tspan)
    result[1] = zero(eltype(result))
    traj_iter = trajectory_iterator(traj)
    result_agg, k = Base.foldl(traj_iter; init=(agg, 1)) do (agg, k), (u, t, i)
        if t <= agg.tspan[1]
            return agg, k
        end

        t = min(t, agg.tspan[2])
        agg = update_rates(agg, (u, t, i), dist.reactions)

        tprev = agg.tprev
        while k <= length(dtimes) && dtimes[k] < t
            result[k] -= (dtimes[k] - tprev) * agg.sumrate
            tprev = dtimes[k]
            k += 1
            result[k] = result[k-1]
        end
        result[k] -= (t - tprev) * agg.sumrate

        log_jump_prob = 0.0
        if i != 0
            @inbounds gid = agg.ridtogroup[i]
            gid != 0 && @inbounds log_jump_prob = log(agg.gsums[gid])
        end
        result[k] += log_jump_prob

        agg = @set agg.weight = agg.weight - (t - agg.tprev) * agg.sumrate + log_jump_prob
        agg = @set agg.tprev = t
        agg, k
    end

    result
end

cumulative_logpdf(dist::TrajectoryDistribution, trajectory, dtimes::AbstractVector) = cumulative_logpdf!(zeros(length(dtimes)), dist, trajectory, dtimes)

@inline @fastmath function evalrxrate(speciesvec::AbstractVector{T}, rxidx::Int64, rs::ReactionSet) where {T}
    val = Float64(1.0)
    @inbounds for specstoch in rs.rstoich[rxidx]
        @inbounds specpop = speciesvec[specstoch[1]]
        val *= Float64(specpop)
        @inbounds for k = 2:specstoch[2]
            specpop -= one(specpop)
            val *= Float64(specpop)
        end
    end

    @inbounds val * rs.rates[rxidx]
end

@inline function update_reaction_rate!(aggregator::AbstractJumpRateAggregator, u::AbstractVector, reactions, rx::Int)
    @inbounds gid = aggregator.ridtogroup[rx]
    if gid != 0
        rate = evalrxrate(u, rx, reactions)
        @inbounds aggregator.gsums[gid] += -aggregator.rates[rx] + rate
        @inbounds aggregator.rates[rx] = rate
    end
end

@fastmath function update_rates(aggregator::DirectAggregator, (u, t, i), reactions::ReactionSet)
    for rx in eachindex(reactions.rates)
        update_reaction_rate!(aggregator, u, reactions, rx)
    end
    aggregator = @set aggregator.sumrate = sum(aggregator.gsums)
end

@fastmath function update_rates(aggregator::DepGraphAggregator, (u, t, i), reactions::ReactionSet)
    if aggregator.prev_reaction == 0
        for rx in eachindex(reactions.rates)
            update_reaction_rate!(aggregator, u, reactions, rx)
        end
    else
        @inbounds for rx in aggregator.depgraph[aggregator.prev_reaction]
            update_reaction_rate!(aggregator, u, reactions, rx)
        end
    end
    aggregator = @set aggregator.sumrate = sum(aggregator.gsums)
    aggregator = @set aggregator.prev_reaction = i
end


text:
#tape interface
module TapeInterface

using ReverseDiffTape
import MathProgBase

const EPS=1e-5

export TapeSolver

###############################################################
#
#  TapeSolver
#
################################################################
type TapeSolver <: MathProgBase.AbstractMathProgSolver
    s::MathProgBase.AbstractMathProgSolver
end

type TapeData <: MathProgBase.AbstractMathProgModel
    evaluator   ##cannot be typed ? why
    m::MathProgBase.AbstractMathProgModel
end

function MathProgBase.NonlinearModel(solver::TapeSolver)
    # @show "TapeSolver model"
    return TapeData(nothing,MathProgBase.NonlinearModel(solver.s))
end

function MathProgBase.setwarmstart!(m::TapeData, x)
    return MathProgBase.setwarmstart!(m.m,x)
end

function MathProgBase.status(m::TapeData)
    e = m.evaluator
    @printf "hesslag_structure_time= %s \n" e.hesslag_structure_time
    @printf "eval_hesslag_time= %s \n" e.eval_hesslag_time
    return MathProgBase.status(m.m)
end

function MathProgBase.getobjval(m::TapeData)
    return MathProgBase.getobjval(m.m)
end

function MathProgBase.getsolution(m::TapeData)
    return MathProgBase.getsolution(m.m)
end

function MathProgBase.getreducedcosts(m::TapeData)
    return MathProgBase.getreducedcosts(m.m)
end

function MathProgBase.getconstrduals(m::TapeData)
    return MathProgBase.getconstrduals(m.m)
end

function MathProgBase.loadproblem!(m::TapeData, numVar, numConstr, l, u, lb, ub, sense, jd::MathProgBase.AbstractNLPEvaluator)
    # @show "TapeData - loadproblem"
    tape_evaluator = TapeNLPEvaluator(jd,numVar,numConstr)
    m.evaluator = tape_evaluator
    MathProgBase.loadproblem!(m.m,numVar,numConstr,l,u,lb,ub,sense,tape_evaluator)
end

function MathProgBase.optimize!(m::TapeData)
    MathProgBase.optimize!(m.m)
end

###############################################################
#
#  TapeEvaluator
#
################################################################
type TapeNLPEvaluator <: MathProgBase.AbstractNLPEvaluator
    jd::MathProgBase.AbstractNLPEvaluator
    numVar::Int
    numConstr::Int
    lag_start::Int
    pvals::Vector{Float64}
    lag_tt::Tape{Int,Float64}
    ttstarts::Vector{Int}
    ttends::Vector{Int}
    trends::Vector{Int}
    hs::HessStorage{Int,Float64}

    nl_idxes::Vector{Int}
    init::Int
    enable_timing_stats::Bool

    jac_I::Vector{Int}
    jac_J::Vector{Int}
    jac_nnz::Int
    jac_idxmap::Vector{Int}
    JJ::Vector{Float64}

    laghess_I::Vector{Int}
    laghess_J::Vector{Int}
    laghess_nnz::Int

    #working array
    stk::Vector{Float64}
    vals::Vector{Float64}
    imm::Vector{Float64}

    # timers
    tape_build::Float64
    eval_f_time::Float64
    eval_g_time::Float64
    eval_grad_f_time::Float64
    jac_structure_time::Float64
    eval_jac_g_time::Float64
    hesslag_structure_time::Float64
    eval_hesslag_time::Float64

    #jump nlp evaluator timer
    jd_init::Float64
end

function TapeNLPEvaluator(nlpe::MathProgBase.AbstractNLPEvaluator,numVar,numConstr;with_timing=false)
    e = TapeNLPEvaluator(nlpe,
        numVar,numConstr,0,
        Vector{Float64}(), #parameter values

        Tape{Int,Float64}(), #lag_tt
        Vector{Int}(),
        Vector{Int}(),
        Vector{Int}(),
        HessStorage{Int,Float64}(0,-1),

        Vector{Int}(), #nl_idxes
        -1, #init
        with_timing,

        Vector{Int}(), Vector{Int}(),-1, Vector{Int}(), Vector{Float64}(), #Jacobian
        Vector{Int}(), Vector{Int}(),-1, #Hessian

        Vector{Float64}(),Vector{Float64}(),Vector{Float64}(), #working array

        0.0, 0.0, 0.0, 0.0,    #timer
        0.0, 0.0, 0.0, 0.0
        ,0.0
        )
    finalizer(e, evaluator_cleanup)
    return e
end

function evaluator_cleanup(e::TapeNLPEvaluator)
    @timing e.enable_timing_stats begin
        f = open("tape_build.txt","w")
        writedlm(f,e.tape_build)
        close(f)

        f = open("eval_f_time.txt","w")
        writedlm(f, e.eval_f_time)
        close(f)

        f = open("eval_g_time.txt","w")
        writedlm(f, e.eval_g_time)
        close(f)

        f = open("eval_grad_f_time.txt","w")
        writedlm(f,e.eval_grad_f_time)
        close(f)

        f = open("jac_structure_time.txt","w")
        writedlm(f, e.jac_structure_time)
        close(f)

        f = open("eval_jac_g_time.txt","w")
        writedlm(f, e.eval_jac_g_time)
        close(f)

        f = open("hesslag_structure_time.txt","w")
        writedlm(f, e.hesslag_structure_time)
        close(f)

        f = open("eval_hesslag_time.txt","w")
        writedlm(f, e.eval_hesslag_time)
        close(f)
    end
end

function MathProgBase.initialize(d::TapeNLPEvaluator, requested_features::Vector{Symbol})
    @assert d.init == -1
    # @assert length(d.pvals) == length(d.ttstarts) == length(d.ttends) == length(d.trends) == 0  ## assume no values in pvals
    jd = d.jd
    gnvar = MathProgBase.numvar(jd.m)

    @timing d.enable_timing_stats tic()
    MathProgBase.initialize(jd,[:ExprGraph])
    @timing d.enable_timing_stats d.jd_init += toq()

    objexpr = MathProgBase.obj_expr(jd)

#indexing approach
    pholder = Vector{Int}()
    bigT = tapeBuilderNoHess(objexpr,d.pvals,gnvar)
    push!(d.ttstarts, 1)
    push!(d.ttends,length(bigT.tt))
    push!(d.trends,length(bigT.tr))
    appendMultParam(bigT,pholder)

    ntts = 1
    obj_nvnode = bigT.nvnode

    for i = 1:d.numConstr
        conexpr = MathProgBase.constr_expr(jd,i)
        if !MathProgBase.isconstrlinear(jd,i)
            push!(d.nl_idxes,i)
        end
        j = length(conexpr.args) == 3?1:3
        @timing d.enable_timing_stats tic()
        tt = tapeBuilderNoHess(conexpr.args[j],d.pvals, gnvar)
        @timing d.enable_timing_stats d.tape_build += toq()
        appendTapeMultParam(bigT, tt, d.pvals,d.ttstarts,d.ttends,d.trends,pholder)
        ntts +=1
    end
    @assert length(pholder) == 1+d.numConstr

    j = length(d.pvals)
    d.lag_start = j + 2
    for i=1:length(pholder)
        @assert bigT.tt[pholder[i]] == -2
        bigT.tt[pholder[i]] = j + i
    end
    append!(d.pvals,ones(Float64,ntts))
    @assert length(d.pvals) == (d.lag_start-1+d.numConstr)
    d.lag_tt = buildSumTape(bigT.tt, bigT.depth,ntts,gnvar,d.hs)
    @assert length(d.ttstarts) == length(d.ttends) == length(d.trends) == (1 + d.numConstr)
##### otherwise
    @show d.lag_tt.depth
    @assert (d.lag_start + d.numConstr - 1) == length(d.pvals)  "$(length(d.pvals)) $(d.numConstr) $(d.lag_start)"
    #end building single lag

    resize!(d.stk, d.lag_tt.stklen)
    resize!(d.vals, d.lag_tt.nnode)
    resize!(d.imm, d.lag_tt.immlen)

    # resize!(d.jac_I,(d.lag_tt.nvnode - obj_nvnode))
    # resize!(d.jac_J,(d.lag_tt.nvnode - obj_nvnode))
    resize!(d.JJ,(d.lag_tt.nvnode - obj_nvnode))
    resize!(d.jac_idxmap, (d.lag_tt.nvnode - obj_nvnode))
    fill!(d.jac_idxmap,0)

    @show "end initialize "
    d.init = 1  #set initialized
    nothing
end

MathProgBase.features_available(d::TapeNLPEvaluator) = [:Grad, :Jac, :Hess]

#evaluate the objective function given on iterate x
function MathProgBase.eval_f(d::TapeNLPEvaluator, x)
    @timing d.enable_timing_stats tic()
    tt = d.lag_tt.tt
    @inbounds ts = d.ttstarts[1]
    @inbounds te = d.ttends[1]
    v = feval(ts,te,tt,x,d.pvals,d.stk)
    @timing d.enable_timing_stats d.eval_f_time += toq()

    return v
end

#evaluate the objective gradient on given iterate x. Results is set to g
function MathProgBase.eval_grad_f(d::TapeNLPEvaluator, g, x)
    # @timing d.enable_timing_stats tic()
    tt = d.lag_tt.tt
    @inbounds ts = d.ttstarts[1]
    @inbounds te = d.ttends[1]
    grad_reverse_dense(ts,te,tt,x,d.pvals,d.vals,d.stk, d.imm, g)
    # @timing d.enable_timing_stats d.eval_grad_f_time += toq()

    return nothing
end

#constraint evaluation
function MathProgBase.eval_g(d::TapeNLPEvaluator, g, x)
    @timing d.enable_timing_stats tic()
    tt = d.lag_tt.tt
    @inbounds for i=1:d.numConstr
        # @inbounds tt = d.constr_tt[i].tt
        @inbounds ts = d.ttstarts[i+1]
        @inbounds te = d.ttends[i+1]
        @inbounds g[i]=feval(ts,te,tt,x,d.pvals,d.stk)
    end
    @timing d.enable_timing_stats d.eval_g_time += toq()

    return nothing
end

function MathProgBase.jac_structure(d::TapeNLPEvaluator)
    @assert d.jac_nnz == -1
    @assert length(d.jac_I) == length(d.jac_J) == 0

    nnz_repeat = length(d.JJ)
    II = Vector{Int}(nnz_repeat)
    JJ = Vector{Int}(nnz_repeat)

    # @timing d.enable_timing_stats tic()
    tt = d.lag_tt.tt
    start = 1
    for i=1:d.numConstr
        # @inbounds tt = d.constr_tt[i].tt
        @inbounds ts = d.ttstarts[i+1]
        @inbounds te = d.ttends[i+1]
        nz = grad_structure(ts, te, tt, start, JJ)
        @inbounds fill!(sub(II, start:(start+nz-1)),i)
        start += nz
    end
    @assert (start - 1) == nnz_repeat
    @assert length(II) == length(JJ) == length(d.JJ) == (start - 1) "$(tt_i.nzg) $(tt_i.nvnode) $(length(d.jac_J)) $(length(d.jac_I))"

    #postprocess non-repeated
    mergedindices = zeros(Int,nnz_repeat)
    @inbounds mergednnz = [0]
    mergedmap = zeros(Int,nnz_repeat)
    function combine(idx1,idx2)
        # @show idx1, idx2
        @inbounds if mergedmap[idx1] == 0 && mergedmap[idx2] != 0
            @assert false
            @inbounds mergednnz[1] += 1
            @inbounds mergedmap[idx1] = idx2
            @inbounds mergedindices[mergednnz[1]] = idx1
            return idx2
        else
            @inbounds @assert mergedmap[idx2] == 0
            @inbounds mergednnz[1] += 1
            @inbounds mergedmap[idx2] = idx1
            @inbounds mergedindices[mergednnz[1]] = idx2
            return idx1
        end
    end

    @assert length(d.jac_I) == length(d.jac_J) == 0
    Jmat = sparse(II, JJ, [i for i in 1:length(JJ)], d.numConstr, MathProgBase.numvar(d.jd.m), combine)

    nnz_non_repeat = length(Jmat.nzval)
    pre = length(d.jac_I)
    start = pre + 1
    append!(d.jac_I, zeros(Int,nnz_non_repeat))
    append!(d.jac_J, zeros(Int,nnz_non_repeat))

    @assert length(d.jac_idxmap) == nnz_repeat
    for col in 1:MathProgBase.numvar(d.jd.m)
        @inbounds for pos in Jmat.colptr[col]:(Jmat.colptr[col+1]-1)
            @inbounds row = Jmat.rowval[pos]
            @inbounds origidx = Jmat.nzval[pos] # this is the original index (on JJ) of this element
            @inbounds d.jac_idxmap[origidx] = pos
            @inbounds d.jac_I[start] = row
            @inbounds d.jac_J[start] = col
            start += 1
        end
    end

    @inbounds for k in 1:mergednnz[1]
        @inbounds origidx = mergedindices[k]
        @inbounds mergedwith = mergedmap[origidx]
        @inbounds @assert d.jac_idxmap[origidx] == 0
        @inbounds @assert d.jac_idxmap[mergedwith] != 0
        @inbounds d.jac_idxmap[origidx] = d.jac_idxmap[mergedwith]
    end

    for i in 1:nnz_repeat
        @inbounds @assert d.jac_idxmap[i] != 0
    end

    d.jac_nnz = nnz_non_repeat

    # @show "end jac_structure"
    return d.jac_I, d.jac_J
end

function MathProgBase.eval_jac_g(d::TapeNLPEvaluator, J, x)
    @assert d.jac_nnz != -1  #structure already computed
    @assert length(J) == d.jac_nnz

    tt = d.lag_tt.tt
    j = 1
    for i = 1:d.numConstr
        # @inbounds tt = d.constr_tt[i].tt
        @inbounds ts = d.ttstarts[i+1]
        @inbounds te = d.ttends[i+1]
        # @timing d.enable_timing_stats tic()
        nz = grad_reverse(ts,te,tt,x,d.pvals,d.vals,d.stk,d.imm,j,d.JJ)
        j += nz
        # @timing d.enable_timing_stats d.eval_jac_g_time += toq()
    end
    @assert (j-1) == length(d.JJ) "$j $(length(d.JJ)) $(d.jac_nnz)"

    nnz = length(d.JJ)
    fill!(J,0)
    for i = 1:nnz
        @inbounds J[d.jac_idxmap[i]] += d.JJ[i]
    end

    nothing
end

function MathProgBase.eval_hesslag_prod(
    d::TapeNLPEvaluator,
    h::Vector{Float64}, # output vector
    x::Vector{Float64}, # current solution
    v::Vector{Float64}, # rhs vector
    σ::Float64,         # multiplier for objective
    μ::Vector{Float64}) # multipliers for each constraint

    #@assert false
    MathProgBase.eval_hesslag_prod(d,h,x,v,σ,μ)
end

function MathProgBase.hesslag_structure(d::TapeNLPEvaluator)
    @assert d.laghess_nnz == -1
    @assert length(d.laghess_I) == length(d.laghess_J) == 0

    @timing true tic()

    #with single constraints block
    tt = d.lag_tt.tt
    tr = d.lag_tt.tr
    nz = hess_structure(1, length(tt), tt, length(tr), tr, d.hs, d.laghess_I, d.laghess_J)
    #end single constraint block

    d.laghess_nnz = length(d.laghess_I)
    @timing true d.hesslag_structure_time += toq()

    @assert d.laghess_nnz == length(d.laghess_J) == length(d.laghess_I) == nz
    @show "end hesslag_structure"
    return  d.laghess_I, d.laghess_J
end

function MathProgBase.eval_hesslag(
    d::TapeNLPEvaluator,
    H::Vector{Float64},         # Sparse hessian entry vector
    x::Vector{Float64},         # Current solution
    obj_factor::Float64,        # Lagrangian multiplier for objective
    lambda::Vector{Float64})    # Multipliers for each constraint

    @assert (length(d.pvals)-d.lag_start+1) == length(lambda)
    @assert length(H) == d.laghess_nnz

    @timing true tic()
    #with single constraint block
    @inbounds d.pvals[d.lag_start-1] = obj_factor

    j=1
    @simd for i=d.lag_start:length(d.pvals)
        @inbounds d.pvals[i] = lambda[j]
        j+=1
    end
    @assert (j-1) == d.numConstr
    tt = d.lag_tt.tt
    tr = d.lag_tt.tr
    hess_reverse(1, length(tt), tt, length(tr), tr ,x,d.pvals,d.stk, d.vals, d.imm, d.hs, 1,H)
    #end with single constraint block

    @timing true d.eval_hesslag_time += toq()
    return
end

end #end module
#########

ext: jl


text:
module FizzBuzzQuiz

"""文字列連結（多重再定義）"""
function Base.:*(n::Int, s2::AbstractString)::AbstractString
    str = String(UInt8[reinterpret(UInt8, [n]);codeunits(s2)])
    SubString(str, sizeof(Int)+1, lastindex(str))
end

function Base.:*(s1::SubString{String}, s2::AbstractString)::AbstractString
    str = String(UInt8[codeunits(s1.string);codeunits(s2)])
    SubString(str, s1.offset+1, lastindex(str))
end

"""文字列に埋め込んだ整数値を抽出（多重定義）"""
function choose_n(str::SubString{String})::Int
    reinterpret(Int, codeunits(str.string)[1:str.offset])[1]
end
choose_n(n::Int)::Int = n

"""fizz()"""
function fizz(val::Union{Int,SubString{String}})::Union{Int,AbstractString}
    choose_n(val) % 3 == 0 ? val * "Fizz" : val
end

"""buzz()"""
function buzz(val::Union{Int,SubString{String}})::Union{Int,AbstractString}
    choose_n(val) % 5 == 0 ? val * "Buzz" : val
end

"""pezz()"""
function pezz(val::Union{Int,SubString{String}})::Union{Int,AbstractString}
    choose_n(val) % 7 == 0 ? val * "Pezz" : val
end

"""hozz()"""
function hozz(val::Union{Int,SubString{String}})::Union{Int,AbstractString}
    3 ∈ digits(choose_n(val)) ? val * "Aho" : val
end

export fizz, buzz, pezz, hozz

end # module


text:
using Printf
using PrettyTables

const bad_compare_msg(obj,exp,got) = @sprintf("Expected %s to be %s, instead got %s",string(obj),string(exp),string(got))

@overload
function Concepts.check(object::Type{Val{:rank}},of::Matrix{T},is::Optional{Integer}=nothing) where T<:Number
    if isnothing(is)
        return LinearAlgebra.rank(of)
    end
    got_rank = LinearAlgebra.rank(of)
    if got_rank == is
        return true
    end
    @warn bad_compare_msg(:rank,is,got_rank)
    return false

end

@overload
function Concepts.check(object::Type{Val{:dimension}},of::Array{T,2},is::Optional{Tuple}=nothing) where T<:Any
    if isnothing(is)
        return size(of)
    end
    got_size = size(of)
    if got_size == is
        return true
    end
    @warn bad_compare_msg(:dimension,is,got_size)
    return false
end

@overload
function Concepts.check(object::Union{Type{Val{:l2difference}},Type{Val{:l2diff}}},
                        a::Union{T1,Array{T1}},b::Union{T2,Array{T2}},
                        against::Optional{S}=nothing) where {T1<:Number,T2<:Number,S<:Real}
    if isnothing(against)
        return LinearAlgebra.norm(a-b,2)
    end
    got_diff = LinearAlgebra.norm(a-b,2)
    if abs(got_diff - against) < 1e-5
        return true
    end
    @warn bad_compare_msg(:l2difference,against,got_diff)
    return false
end

@overload
function Base.zeros(like::Array{T}) where T<:Any
    return zeros(size(like))
end

# This is necessary due to compiler bug??
# @overload
# const Concepts.check(object::Symbol,arg1) = Concepts.check(Val{object},arg1)

# @overload
# const Concepts.check(object::Symbol,arg1,arg2) = Concepts.check(Val{object},arg1,arg2)

# @overload
# const Concepts.check(object::Symbol,arg1,arg2,arg3) =  Concepts.check(Val{object},arg1,arg2,arg3)


text:
# This file is auto-generated by AWSMetadata.jl
using AWS
using AWS.AWSServices: timestream_query
using AWS.Compat
using AWS.UUIDs

"""
    cancel_query(query_id)
    cancel_query(query_id, params::Dict{String,<:Any})

 Cancels a query that has been issued. Cancellation is guaranteed only if the query has not
completed execution before the cancellation request was issued. Because cancellation is an
idempotent operation, subsequent cancellation requests will return a CancellationMessage,
indicating that the query has already been canceled.

# Arguments
- `query_id`:  The id of the query that needs to be cancelled. QueryID is returned as part
  of QueryResult.

"""
function cancel_query(QueryId; aws_config::AbstractAWSConfig=global_aws_config())
    return timestream_query(
        "CancelQuery",
        Dict{String,Any}("QueryId" => QueryId);
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function cancel_query(
    QueryId, params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()
)
    return timestream_query(
        "CancelQuery",
        Dict{String,Any}(mergewith(_merge, Dict{String,Any}("QueryId" => QueryId), params));
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

"""
    describe_endpoints()
    describe_endpoints(params::Dict{String,<:Any})

DescribeEndpoints returns a list of available endpoints to make Timestream API calls
against. This API is available through both Write and Query. Because Timestream’s SDKs
are designed to transparently work with the service’s architecture, including the
management and mapping of the service endpoints, it is not recommended that you use this
API unless:   Your application uses a programming language that does not yet have SDK
support   You require better control over the client-side implementation   For detailed
information on how to use DescribeEndpoints, see The Endpoint Discovery Pattern and REST
APIs.

"""
function describe_endpoints(; aws_config::AbstractAWSConfig=global_aws_config())
    return timestream_query(
        "DescribeEndpoints"; aws_config=aws_config, feature_set=SERVICE_FEATURE_SET
    )
end
function describe_endpoints(
    params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()
)
    return timestream_query(
        "DescribeEndpoints", params; aws_config=aws_config, feature_set=SERVICE_FEATURE_SET
    )
end

"""
    query(query_string)
    query(query_string, params::Dict{String,<:Any})

 Query is a synchronous operation that enables you to execute a query. Query will timeout
after 60 seconds. You must update the default timeout in the SDK to support a timeout of 60
seconds. The result set will be truncated to 1MB. Service quotas apply. For more
information, see Quotas in the Timestream Developer Guide.

# Arguments
- `query_string`:  The query to be executed by Timestream.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"ClientToken"`:  Unique, case-sensitive string of up to 64 ASCII characters that you
  specify when you make a Query request. Providing a ClientToken makes the call to Query
  idempotent, meaning that multiple identical calls have the same effect as one single call.
  Your query request will fail in the following cases:    If you submit a request with the
  same client token outside the 5-minute idepotency window.     If you submit a request with
  the same client token but a change in other parameters within the 5-minute idempotency
  window.     After 4 hours, any request with the same client token is treated as a new
  request.
- `"MaxRows"`:  The total number of rows to return in the output. If the total number of
  rows available is more than the value specified, a NextToken is provided in the command's
  output. To resume pagination, provide the NextToken value in the starting-token argument of
  a subsequent command.
- `"NextToken"`:  A pagination token passed to get a set of results.
"""
function query(QueryString; aws_config::AbstractAWSConfig=global_aws_config())
    return timestream_query(
        "Query",
        Dict{String,Any}("QueryString" => QueryString, "ClientToken" => string(uuid4()));
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function query(
    QueryString,
    params::AbstractDict{String};
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return timestream_query(
        "Query",
        Dict{String,Any}(
            mergewith(
                _merge,
                Dict{String,Any}(
                    "QueryString" => QueryString, "ClientToken" => string(uuid4())
                ),
                params,
            ),
        );
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end


text:
#=
    The File contains types and interfaces definitions
=#

# All interfaces defined here must provide one-based indexing,
# yet internal indexing is implementation-defined

#------------------------------------------------------------------------------

include("policiks.jl")

#------------------------------------------------------------------------------

import Base: ==, !=, +, -, *
import Nemo: gfp_elem, gfp_fmpz_elem
import AbstractAlgebra: Field

#------------------------------------------------------------------------------

# AbstractSparseObject<T>
#
# `T` is the type of the base_ring field
#
# Base interface for anything Sparse implemented in this project
abstract type AbstractSparseObject{T} end

#------------------------------------------------------------------------------

# AbstractSparseVector<T>
#
# `T` is the type if the base_ring field
#
# Base interface for Sparse vectors family
# An object which implements `AbstractSparseVector` can be treated
# as an element of Rⁿ vectorspace. Thus, some vectorspace operations and
# left matrix multiplication are supported for such objects
abstract type AbstractSparseVector{T} <: AbstractSparseObject{T} end

#------------------------------------------------------------------------------

# AbstractSparseMatrix<T>
#
# `T` is the type if the base_ring field
#
# Base interface for Sparse matrices family
# An object which implements `AbstractSparseMatrix` can be treated
# as an element from a Matrix algebra
# An object which implements `AbstractSparseMatrix` can be treated
# as a linear operator in Rⁿ
#
#
# Implementing some functions for other sparse formats
# rather than DOK is tricky,
# probably, this interface it too strong to be applied to every class
abstract type AbstractSparseMatrix{T} <: AbstractSparseObject{T} end

#------------------------------------------------------------------------------

function base_ring(::AbstractSparseMatrix) end

function first_nonzero(::T) where {T<:AbstractSparseMatrix} end
function getindex(::AbstractSparseMatrix, i::Int) end

function scale(::AbstractSparseMatrix, ::C) where {C} end
function scale!(::AbstractSparseMatrix, ::C) where {C} end

# To Be Changed sometime
# !!!
# A tender for a better name!
function Base.reduce(::T, ::T, ::C) where {T<:AbstractSparseMatrix{F}} where {F, C} end
function reduce!(::T, ::T, ::C) where {T<:AbstractSparseMatrix{F}} where {F, C} end

function inner(::AbstractSparseMatrix, ::AbstractSparseMatrix) end

# hmmm
function Base.prod(::T, ::T) where {T<:AbstractSparseMatrix{F}} where {F} end

function Base.zero(::AbstractSparseMatrix) end
function Base.isempty(::AbstractSparseMatrix) end
function Base.iszero(::AbstractSparseMatrix) end

function dim(::AbstractSparseMatrix) end
function Base.length(::AbstractSparseMatrix) end
function Base.size(::AbstractSparseMatrix) end
function Base.size(::AbstractSparseMatrix, i::Int) end

function density(::AbstractSparseMatrix) end

function ==(::T, ::T) where {T<:AbstractSparseMatrix{F}} where {F} end
function !=(::T, ::T) where {T<:AbstractSparseMatrix{F}} where {F} end

function Base.valtype(::AbstractSparseMatrix) end
function Base.eltype(::AbstractSparseMatrix) end

function Base.show(::AbstractSparseMatrix) end
function print_sparsik(::AbstractSparseMatrix) end

function modular_reduction(::AbstractSparseMatrix, ::T) where {T} end
function rational_reconstruction(::AbstractSparseMatrix) end

function Base.iterate(::AbstractSparseMatrix, state) end

function unit_vector(::Int, ::Int, ::T) where {T} end

function random_sparsik(::Tuple{Int}, ::T) where {T} end

#------------------------------------------------------------------------------

function Base.convert(::Type{Int}, x::gfp_elem)
    return Int(x.data)
end

function Base.convert(::Type{BigInt}, x::gfp_fmpz_elem)
    return BigInt(x.data)
end

#------------------------------------------------------------------------------


text:
# *********************************************************************************
# REopt, Copyright (c) 2019-2020, Alliance for Sustainable Energy, LLC.
# All rights reserved.
#
# Redistribution and use in source and binary forms, with or without modification,
# are permitted provided that the following conditions are met:
#
# Redistributions of source code must retain the above copyright notice, this list
# of conditions and the following disclaimer.
#
# Redistributions in binary form must reproduce the above copyright notice, this
# list of conditions and the following disclaimer in the documentation and/or other
# materials provided with the distribution.
#
# Neither the name of the copyright holder nor the names of its contributors may be
# used to endorse or promote products derived from this software without specific
# prior written permission.
#
# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
# ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
# WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
# IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT,
# INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,
# BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
# DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
# LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE
# OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED
# OF THE POSSIBILITY OF SUCH DAMAGE.
# *********************************************************************************
function add_flexible_hvac_constraints(m, p::REoptInputs; _n="")

    binFlexHVAC = @variable(m, binary = true)
    (N, J) = size(p.s.flexible_hvac.input_matrix)
    dvTemperature = @variable(m, [1:N, p.time_steps])
    @variable(m, lower_comfort_slack[p.time_steps] >= 0)
    @variable(m, upper_comfort_slack[p.time_steps] >= 0)

    # initialize space temperatures
    @constraint(m, dvTemperature[:, 1] .== p.s.flexible_hvac.initial_temperatures)

    # TODO time scaling for dt?
    input_vec = zeros(N)
    input_vec[p.s.flexible_hvac.control_node] = 1

    if !isempty(p.techs.heating) && !isempty(p.techs.cooling)
        # space temperature evolution based on state-space model
        @constraint(m, [n in 1:N, ts in 2:length(p.time_steps)],
            binFlexHVAC => { dvTemperature[n, ts] == dvTemperature[n, ts-1] +
                sum(p.s.flexible_hvac.system_matrix[n, i] * dvTemperature[i, ts-1] for i=1:N) +
                sum(p.s.flexible_hvac.input_matrix[n, j] * p.s.flexible_hvac.exogenous_inputs[j, ts-1] for j=1:J) +
                input_vec[n] * p.s.flexible_hvac.input_matrix[n, p.s.flexible_hvac.control_node] * (
                    sum(m[Symbol("dvThermalProduction"*_n)][t, ts-1] for t in p.techs.heating) -
                    sum(m[Symbol("dvThermalProduction"*_n)][t, ts-1] for t in p.techs.cooling)
                )}
        )
        @constraint(m, [ts in p.time_steps],
            p.s.flexible_hvac.temperature_lower_bound_degC - lower_comfort_slack[ts] <=
            dvTemperature[p.s.flexible_hvac.control_node, ts]
        )
        @constraint(m, [ts in p.time_steps],
            dvTemperature[p.s.flexible_hvac.control_node, ts] <=
            p.s.flexible_hvac.temperature_upper_bound_degC + upper_comfort_slack[ts]
        )

    elseif !isempty(p.techs.heating)

        @constraint(m, [n in 1:N, ts in 2:length(p.time_steps)],
            binFlexHVAC => { dvTemperature[n, ts] == dvTemperature[n, ts-1] +
            sum(p.s.flexible_hvac.system_matrix[n, i] * dvTemperature[i, ts-1] for i=1:N) +
            sum(p.s.flexible_hvac.input_matrix[n, j] * p.s.flexible_hvac.exogenous_inputs[j, ts-1] for j=1:J) +
            input_vec[n] * p.s.flexible_hvac.input_matrix[n, p.s.flexible_hvac.control_node] * (
                sum(m[Symbol("dvThermalProduction"*_n)][t, ts-1] for t in p.techs.heating)
            )}
        )
        @constraint(m, [ts in p.time_steps],
            p.s.flexible_hvac.temperature_lower_bound_degC - lower_comfort_slack[ts] <=
            dvTemperature[p.s.flexible_hvac.control_node, ts]
        )
        # when only heating the upper temperature limit is the highest temperature seen naturally
        @constraint(m, [ts in p.time_steps],
            dvTemperature[p.s.flexible_hvac.control_node, ts] <=
            maximum(p.s.flexible_hvac.bau_hvac.temperatures[p.s.flexible_hvac.control_node, :]) +
            upper_comfort_slack[ts]
        )

    elseif !isempty(p.techs.cooling)

        @constraint(m, [n in 1:N, ts in 2:length(p.time_steps)],
            binFlexHVAC => { dvTemperature[n, ts] == dvTemperature[n, ts-1] +
            sum(p.s.flexible_hvac.system_matrix[n, i] * dvTemperature[i, ts-1] for i=1:N) +
            sum(p.s.flexible_hvac.input_matrix[n, j] * p.s.flexible_hvac.exogenous_inputs[j, ts-1] for j=1:J) -
            input_vec[n] * p.s.flexible_hvac.input_matrix[n, p.s.flexible_hvac.control_node] * (
                sum(m[Symbol("dvThermalProduction"*_n)][t, ts-1] for t in p.techs.cooling)
            )}
        )
        # when only cooling the lower temperature limit is the lowest temperature seen naturally
        @constraint(m, [ts in p.time_steps],
            minimum(p.s.flexible_hvac.bau_hvac.temperatures[p.s.flexible_hvac.control_node, :]) -
            lower_comfort_slack[ts] <=
            dvTemperature[p.s.flexible_hvac.control_node, ts]
        )
        @constraint(m, [ts in p.time_steps],
            dvTemperature[p.s.flexible_hvac.control_node, ts] <=
            p.s.flexible_hvac.temperature_upper_bound_degC + upper_comfort_slack[ts]
        )
    end

    dvComfortLimitViolationCost = @expression(m,
        1e9 * sum(lower_comfort_slack[ts] + upper_comfort_slack[ts] for ts in p.time_steps)
    )
    # TODO convert dvThermalProduction units? to ? shouldn't the conversion be in input_matrix coef? COP in Xiang's test is 4-5, fan_power_ratio = 0, hp prod factor generally between 1 and 2
    ## TODO check eigen values / stability of system matrix?

    @variable(m, dvFlexHVACcost >= 0)
    @constraint(m, binFlexHVAC => { dvFlexHVACcost >= p.s.flexible_hvac.installed_cost})
    m[:TotalTechCapCosts] += dvFlexHVACcost

    # If not buying FlexibleHVAC then the BAU (deadband) thermal loads must be met
    # TODO account for different tech efficiencies in following?

    if !isempty(p.techs.heating)
        @constraint(m, [ts in p.time_steps],
            !binFlexHVAC => { sum(m[Symbol("dvThermalProduction"*_n)][t, ts] for t in p.techs.heating) == p.s.flexible_hvac.bau_hvac.existing_boiler_kw_thermal[ts]
            }
        )
    end
    if !isempty(p.techs.cooling)
        @constraint(m, [ts in p.time_steps],
            !binFlexHVAC => { sum(m[Symbol("dvThermalProduction"*_n)][t, ts] for t in p.techs.cooling) == p.s.flexible_hvac.bau_hvac.existing_chiller_kw_thermal[ts]
            }
        )
    end

    m[Symbol("binFlexHVAC"*_n)] = binFlexHVAC
    m[Symbol("dvTemperature"*_n)] = dvTemperature
    m[Symbol("dvComfortLimitViolationCost"*_n)] = dvComfortLimitViolationCost
    nothing
end

"""
    function add_flexible_hvac_constraints(m, p::REoptInputs{BAUScenario}; _n="")

For the BAU scenario we enforce the deadband control pattern using the values from BAU_HVAC.
"""
function add_flexible_hvac_constraints(m, p::REoptInputs{BAUScenario}; _n="")

    # If not buying FlexibleHVAC then the BAU (deadband) thermal loads must be met
    # TODO account for different tech efficiencies in following?

    if !isempty(p.techs.heating)
        @constraint(m, [ts in p.time_steps],
            sum(m[Symbol("dvThermalProduction"*_n)][t, ts] for t in p.techs.heating) ==
            p.s.flexible_hvac.existing_boiler_kw_thermal[ts]
        )
    end
    if !isempty(p.techs.cooling)
        @constraint(m, [ts in p.time_steps],
            sum(m[Symbol("dvThermalProduction"*_n)][t, ts] for t in p.techs.cooling) ==
            p.s.flexible_hvac.existing_chiller_kw_thermal[ts]
        )
    end

    m[Symbol("binFlexHVAC"*_n)] = 0
    m[Symbol("dvTemperature"*_n)] = p.s.flexible_hvac.temperatures
    m[Symbol("dvComfortLimitViolationCost"*_n)] = 0.0
    nothing
end

max_stars_repo_name: NREL/REopt

max_issues_repo_name: NREL/REopt

max_forks_repo_name: NREL/REopt


text:
chars_only = vcat(collect('A':'Z'),collect('a':'z'),'"')
digits_only = ['0', '1','2','3','4','5','6','7','8','9']
numbers =  vcat(digits_only,['+','-','.'])

chars_numbers = vcat(chars_only,numbers)
separators = ['[', ']','{','}','(',')',';',',','-','*','/',' ','<','>','=',' ']
valid_chars = vcat(chars_numbers,separators,':')

function checkisint(lexem::String)
	try
		parse(Int,lexem)
		if (contains(==,lexem,'.'))
			return false
		end
			return true
	catch
		return false
	end
end;

function checkisfloat(lexem::String)
	try
		parse(Float64,lexem)
		return true
	catch
		return false
	end
end;

checkisnumber(lexem::String) = checkisfloat(lexem)

function matchlexem(item::String)

	if (item == "!=")
		return Int(OPRLR_REL)
	end
	if (item == ">")
		return Int(OPRLR_REL)
	end

	#Checks if the lexem has invalid characters
	invalidpart =
	filter(item) do x
			!contains(==,vcat(valid_chars),x)
	end

	if (length(invalidpart) > 0)
		return Int(LEX_ERR)
	end

	if (contains(==,digits_only,item[1]) && !checkisnumber(item)) #if the lexem starts with a number, but the lexem itself is not a number
		return Int(LEX_ERR)
	end

	if checkisnumber(item)
		if (checkisint(item))
			return Int(CT_INT)
		else
			return Int(CT_FLOAT)
		end
	end

	if item[1] == '"' && item[length(item)] == '"' #if it starts and ends with "
		return Int(CT_STRING)
	end

	#TODO ISSO AQUI TÁ PERMITINDO CHAR DE TAMANHO >1
	if item[1] == '\'' && item[length(item)] == '\'' #if it starts and ends with \)
		return Int(CT_CHAR)
	end

	#From this point we are sure the token is not a number neither is a string

	return	@match item begin
				"for"=>Int(FOR)
				"while"=>Int(WHILE)
				"return"=>Int(RETURN)
				"~" => Int(OPR_UN_NEG)
				"-" => Int(OPR_PM)
				"+" => Int(OPR_PM)
				"++" => Int(OPR_CONCAT)
				"*" => Int(OPR_DM)
				"/" => Int(OPR_DM)
				"[" => Int(O_SQRBRCKT)
				"]" => Int(C_SQRBRCKT)
				"(" => Int(O_BRCKT)
				")" => Int(C_BRCKT)
				"," => Int(COMMA)
				";" => Int(SMCL)
				"{" => Int(O_C_BRCKT)
				"}" => Int(C_C_BRCKT)
				"==" => Int(OPRLR_REL)
				"=" => Int(OPR_ATR)
				">" => Int(OPRLR_REL)
				"<" => Int(OPRLR_REL)
				">=" => Int(OPRLR_REL)
				"<=" => Int(OPRLR_REL)
				"!=" => Int(OPRLR_REL)
				"if" => Int(IF)
				"int" => Int(IDT_INT)
				"float" => Int(IDT_FLOAT)
				"void" => Int(IDT_VOID)
				"char" => Int(IDT_CHAR)
				"string" => Int(IDT_STRING)
				"print" => Int(ID)
				"main" => Int(ID)
				"AND" => Int(OPRLR_AND)
				"OR" => Int(OPRLR_OR)
				"NOT" => Int(OPRL_NOT)
				"true" => Int(CT_TRUE)
				"false" => Int(CT_FALSE)
				"#" => Int(CMNT_LN)
				"::" => Int(VEC_IN)
				"EOF" => Int(EOF)
				 _ => false
			end
end

max_stars_repo_name: NaelsonDouglas/namelesscompiler

max_issues_repo_name: NaelsonDouglas/namelesscompiler

max_forks_repo_name: NaelsonDouglas/namelesscompiler


text:
const libldap = joinpath(@__DIR__, "..", "libldap", "lib", "libldap.so")

function check_deps()
    global libldap
    if !isfile(libldap)
        error("\$(libldap) does not exist.")
    end

    if Libdl.dlopen_e(libldap) in (C_NULL, nothing)
        error("\$(libldap) cannot be opened.")
    end
end

#int ldap_is_ldap_url( const char *url )
function ldap_is_ldap_url(url::AbstractString)
    ccall((:ldap_is_ldap_url, libldap), Cint, (Cstring,), url)
end

#=
int ldap_initialize(ldp, uri)
       LDAP **ldp;
       char *uri;
=#
function ldap_initialize(ldp_handle_ref::Ref{Ptr{Cvoid}}, uri::AbstractString)
    ccall((:ldap_initialize, libldap), Cint, (Ref{Ptr{Cvoid}}, Cstring), ldp_handle_ref, uri)
end

#char *ldap_err2string( int err );
function ldap_err2string(err::Integer)
    ccall((:ldap_err2string, libldap), Cstring, (Cint,), err)
end

# int ldap_url_parse( const char *url, LDAPURLDesc **ludpp )
function ldap_url_parse(url::AbstractString, ldap_url_desc_handle_ref::Ref{Ptr{LDAPURLDesc}})
    ccall((:ldap_url_parse, libldap), Cint, (Cstring, Ref{Ptr{LDAPURLDesc}}), url, ldap_url_desc_handle_ref)
end

# void ldap_free_urldesc( LDAPURLDesc *ludp );
function ldap_free_urldesc(ldap_url_desc_handle::Ptr{LDAPURLDesc})
    ccall((:ldap_free_urldesc, libldap), Cvoid, (Ptr{LDAPURLDesc},), ldap_url_desc_handle)
end

# int ldap_get_option(LDAP *ld, int option, void *outvalue);
function ldap_get_option(ldp_handle::Ptr{Cvoid}, option::LDAPOption, out_value_ref::Ref{Cint})
    opt_int = Cint(option)
    ccall((:ldap_get_option, libldap), Cint, (Ptr{Cvoid}, Cint, Ref{Cint}), ldp_handle, opt_int, out_value_ref)
end

# int ldap_set_option(LDAP *ld, int option, const void *invalue);
function ldap_set_option(ldp_handle::Ptr{Cvoid}, option::LDAPOption, in_value::T) where {T}
    opt_int = Cint(option)
    in_value_ref = Ref(Cint(in_value))
    ccall((:ldap_set_option, libldap), Cint, (Ptr{Cvoid}, Cint, Ref{Cint}), ldp_handle, opt_int, in_value_ref)
end

# int ldap_simple_bind_s(LDAP *ld, const char *who, const char *passwd);
function ldap_simple_bind_s(ldp_handle::Ptr{Cvoid}, who::AbstractString, password::AbstractString)
    ccall((:ldap_simple_bind_s, libldap), Cint, (Ptr{Cvoid}, Cstring, Cstring), ldp_handle, who, password)
end

# int ldap_unbind_s(LDAP *ld);
function ldap_unbind_s(ldp_handle::Ptr{Cvoid})
    ccall((:ldap_unbind_s, libldap), Cint, (Ptr{Cvoid},), ldp_handle)
end


text:
abstract type AbstractAbundanceTable{T <: Real,
                                     F <: AbstractFeature,
                                     S <: AbstractSample} <: EcoBase.AbstractAssemblage{T, F, S}
end

"""
    CommunityProfile{T, F, S} <: AbstractAbundanceTable{T, F, S}

An `AbstractAssemblage` from [EcoBase.jl](https://github.com/EcoJulia/EcoBase.jl)
that uses an `AxisArray` of a `SparseMatrixCSC` under the hood.

`CommunityProfile`s are tables with `AbstractFeature`-indexed rows and
`AbstractSample`-indexed columns.
Note - we can use the `name` of samples and features to index.
"""
mutable struct CommunityProfile{T, F, S} <: AbstractAbundanceTable{T, F, S}
    aa::NamedAxisArray

    function CommunityProfile(aa::NamedAxisArray)
        @assert dimnames(aa) == (:features, :samples)
        T = eltype(parent(aa))
        F = eltype(keys(axes(aa, 1)))
        S = eltype(keys(axes(aa, 2)))
        return new{T, F, S}(aa)
    end
end

function CommunityProfile(tab::SparseMatrixCSC{<:Real},
                          features::AbstractVector{<:AbstractFeature},
                          samples::AbstractVector{<:AbstractSample})
    return CommunityProfile(NamedAxisArray(tab, features=features, samples=samples))
end

function CommunityProfile{T, F, S}(tab::SparseMatrixCSC{<:T},
                                   features::AbstractVector{F},
                                   samples::AbstractVector{S}) where {T, F, S}
    return CommunityProfile(tab, features, samples)
end

function CommunityProfile(tab::AbstractMatrix,
                          features::AbstractVector{<:AbstractFeature},
                          samples::AbstractVector{<:AbstractSample})
    return CommunityProfile(sparse(tab), features, samples)
end
## -- Convienience functions -- ##

function ==(p1::CommunityProfile, p2::CommunityProfile)
    return abundances(p1) == abundances(p2) &&
           samples(p1)    == samples(p2) &&
           features(p1)   == features(p2)
end

"""
    features(at::AbstractAbundanceTable)

Returns features in `at`. To get featurenames instead, use [`featurenames`](@ref).
"""
features(at::AbstractAbundanceTable) = axes(at.aa, 1) |> keys

"""
    samples(at::AbstractAbundanceTable)

Returns samples in `at`. To get samplenames instead, use [`samplenames`](@ref).
"""
samples(at::AbstractAbundanceTable) = axes(at.aa, 2) |> keys

"""
    samples(at::AbstractAbundanceTable, name::AbstractString)

Returns sample in `at` with name `name`.
"""
function samples(at::AbstractAbundanceTable, name::AbstractString)
    idx = findall(==(name), samplenames(at))
    length(idx) == 0 && throw(IndexError("No samples called $name"))
    length(idx) > 1 && throw(IndexError("More than one sample matches name $name"))
    return samples(at)[axes(at.aa, 2)][first(idx)]
end

profiletype(at::AbstractAbundanceTable) = eltype(features(at))
ranks(at::AbstractAbundanceTable) = taxrank.(features(at))

Base.size(at::AbstractAbundanceTable, dims...) = size(at.aa, dims...)

Base.copy(at::AbstractAbundanceTable) = typeof(at)(copy(abundances(at)), copy(features(at)), deepcopy(samples(at)))

# -- Indexing -- #

function _index_profile(at, idx, inds)
    # single value - return that value
    ndims(idx) == 0 && return idx
    # another table - return a new CommunityProfile with that table
    ndims(idx) == 2 && return CommunityProfile(idx)
    # a row or a column, figure out which, and make it 2D
    if ndims(idx) == 1
        dn = dimnames(idx)[1]
        # if it's a row...
        if dn == :samples
            return at[[inds[1]], inds[2]]
        # if it's a column
        elseif dn == :features
            return at[inds[1], [inds[2]]]
        end
    end
end

function _toinds(arr, inds::AbstractVector{Regex})
    return findall(a-> any(ind-> contains(a, ind), inds), arr)
end

function _toinds(arr, inds::AbstractVector{<: Union{AbstractSample, AbstractFeature, AbstractString}})
    return findall(a-> any(==(a), inds), arr)
end

# fall back ↑
_toinds(arr, ind::Union{AbstractSample, AbstractFeature, AbstractString, Regex}) = _toinds(arr, [ind])

# if inds are integers, just return them
_toinds(_, ind::Int) = ind
_toinds(_, inds::AbstractVector{Int}) = inds

function Base.getindex(at::CommunityProfile, inds...)
    idx = at.aa[inds...]

    _index_profile(at, idx, inds)
end

function Base.getindex(at::CommunityProfile, rowind::Union{T, AbstractVector{<:T}} where T<:Union{AbstractString,Regex}, colind)
    rows = _toinds(featurenames(at), rowind)
    idx = at.aa[rows, colind]

    _index_profile(at, idx, (rows, colind))
end

function Base.getindex(at::CommunityProfile, rowind, colind::Union{T, AbstractVector{<:T}} where T<:Union{AbstractString,Regex})
    cols = _toinds(samplenames(at), colind)
    idx = at.aa[rowind, cols]

    _index_profile(at, idx, (rowind, cols))
end

function Base.getindex(at::CommunityProfile, rowind::Union{T, AbstractVector{<:T}} where T<:Union{AbstractString,Regex},
                                             colind::Union{S, AbstractVector{<:S}} where S<:Union{AbstractString,Regex})
    rows = _toinds(featurenames(at), rowind)
    at[rows, colind]
end

## -- EcoBase Translations -- ##
# see src/ecobase.jl for Microbiome function names
# thing => feature
# place => sample
# occurrences => abundances (or totals)

EcoBase.thingnames(at::AbstractAbundanceTable) = name.(features(at))
EcoBase.placenames(at::AbstractAbundanceTable) = name.(samples(at))
EcoBase.occurrences(at::AbstractAbundanceTable) = parent(parent(at.aa)) # first parent is the unnamed AxisArray
EcoBase.nthings(at::AbstractAbundanceTable) = size(at, 1)
EcoBase.nplaces(at::AbstractAbundanceTable) = size(at, 2)
## todo
# EcoBase.thingoccurrences(at::AbstractAbundanceTable, things) = nothing
# EcoBase.placeoccurrences(at::AbstractAbundanceTable, places) = nothing

# for custom printing
EcoBase.thingkind(asm::AbstractAbundanceTable) = "feature"
EcoBase.placekind(asm::AbstractAbundanceTable) = "sample"
## not needed for now
# EcoBase.thingkindplural(asm::AbstractAbundanceTable) = "$(thingkind(asm))s"
# EcoBase.placekindplural(asm::AbstractAbundanceTable) = "$(placekind(asm))s"

"""
    featuretotals(at::AbstractAbundanceTable)

Returns sum of each row (feature) in `at`.
Note, return value is a nfeatures x 1 `Matrix`, not a `Vector`.
If you need 1D `Vector`, use `vec(featuretotals(at))`.
"""
featuretotals(at::AbstractAbundanceTable) = sum(abundances(at), dims=2)

"""
    sampletotals(at::AbstractAbundanceTable)

Returns sum of each row (feature) in `at`.
Note, return value is a 1 x nsamples `Matrix`, not a `Vector`.
If you need 1D `Vector`, use `vec(sampletotals(at))`.
"""
sampletotals(at::AbstractAbundanceTable) = sum(abundances(at), dims=1)

## -- Tables Interface -- ##

Tables.istable(::AbstractAbundanceTable) = true
Tables.columnaccess(::AbstractAbundanceTable) = true
Tables.rowaccess(::AbstractAbundanceTable) = true

Tables.getcolumn(at::AbstractAbundanceTable, i::Int) = i == 1 ? featurenames(at) : abundances(at[:, i-1])
Tables.getcolumn(at::AbstractAbundanceTable, i::AbstractString) = i == "features" ? features(at) : abundances(at[:, i])
Tables.getcolumn(at::AbstractAbundanceTable, i::Symbol) = Tables.getcolumn(at, string(i))

Tables.columnnames(at::AbstractAbundanceTable) = [:features, Symbol.(samplenames(at))...]

function Tables.schema(at::AbstractAbundanceTable)
    elt = eltype(abundances(at))
    coltypes = [eltype(features(at)), (elt for _ in 1:nsamples(at))...]
    return Tables.Schema(Tables.columnnames(at), coltypes)
end

Tables.columns(at::AbstractAbundanceTable) = (; (col => Tables.getcolumn(at, col) for col in Tables.columnnames(at))...)

function _makerow(row::AbstractAbundanceTable)
    size(row, 1) == 1 || error("Can't make row from table of size $(size(row))")
    NamedTuple{(:features, Symbol.(samplenames(row))...)}((first(features(row)), abundances(row)...))
end

Tables.rows(at::AbstractAbundanceTable) = (_makerow(at[i, :]) for i in 1:nfeatures(at))

## -- Methods for absolute and relative abundances -- ##

# """
#     filterabund(abun::AbstractComMatrix, n::Int=minimum(10, nfeatures(abun)))

# Filter an abundance table to the top `n` features accross all samples

# This function also adds a row for "other", which sums the abundances of the
# remaining features.
# """
# function filterabund(abun::AbstractAbundanceTable, n::Int=minimum(10, nfeatures(abun)))
#     # TODO: add prevalence filter

#     totals = featuretotals(abun)

#     srt = sortperm(totals, rev=true)

#     newabun = getfeature(abun, srt[1:n])

#     remainder = [sum(occurrences(abun)[srt[n+1:end], i]) for i in 1:size(abun, 2)]'
#     newabun = vcat(newabun, remainder)
#     newrows = cat(featurenames(abun)[srt[1:n]], ["other"], dims=1)

#     return abundancetable(newabun, samplenames(abun), newrows)
# end

"""
    relativeabundance!(a::AbstractAbundanceTable; kind::Symbol=:fraction)

Normalize each sample in AbstractAbundanceTable to the sum of the sample.

By default, columns sum to 1.0.
Use `kind=:percent` for columns to sum to 100.
"""
function relativeabundance!(at::AbstractAbundanceTable; kind::Symbol=:fraction)
    in(kind, [:percent, :fraction]) || throw(ArgumentError("Invalid kind: $kind"))
    eltype(abundances(at)) <: AbstractFloat || throw(ArgumentError("relativeabundance! requires profile to have AbstractFloat eltype. Try relativeabundance instead"))
    abund = abundances(at)
    abund ./= sampletotals(at)
    kind == :percent && (abund .*= 100)
    dropzeros!(abund) # shouldn't need dropzeros - https://github.com/JuliaLang/julia/issues/39018
    return at
end

"""
    relativeabundance(at::AbstractAbundanceTable, kind::Symbol=:fraction)

Like [`relativeabundance!`](@ref), but does not mutate original.
"""
function relativeabundance(at::AbstractAbundanceTable, kind::Symbol=:fraction)
    comm = typeof(at)(float.(abundances(at)), deepcopy(features(at)), deepcopy(samples(at)))
    relativeabundance!(comm)
end

"""
    present(t::Union{Real, Missing}, minabundance::Real=0.0)
    present(at::AbstractAbundanceTable, minabundance::Real=0.0)

Check if a given (non-zero) value is greater than or equal to a minimum value.
If the minimum abundance is 0, just checks if value is non-zero.

If used on an `AbstractAbundanceTable`, returns a sparse boolean matrix of the same size.
"""
function present(t::Real, minabundance::Real=0.0)
    (minabundance >= 0 && t >= 0) || throw(DomainError("Only defined for positive values"))
    t == 0 ? false : t >= minabundance
end

present(::Missing, m) = missing

function present(at::AbstractAbundanceTable, minabundance::Real=0.0)
    mat = spzeros(Bool, Ssize(at)...)
    for i in eachindex(mat)
        mat[i] = present(at[i], minabundance)
    end
    return mat
end

"""
    prevalence(a::AbstractArray{<:Real}, minabundance::Real=0.0)
    prevalence(at::AbstractAbundanceTable, minabundance::Real=0.0)

Return the fraction of values that are greater than or equal to a minimum.
If the minimum abundance is 0, returns the fraction of non-zero values.

If used on an `AbstractAbundanceTable`,
returns a prevalence value for each `feature` accross the `sample`s.
"""
prevalence(a::AbstractArray{<:Real}, minabundance::Real=0.0) = mean(x-> present(x, minabundance), a)

# makes it work for any iterable
prevalence(a, minabundance::Real=0.0) = mean(x-> present(x, minabundance), (y for y in a))

function prevalence(at::AbstractAbundanceTable, minabundance::Real=0.0)
    mean(x-> present(x, minabundance), abundances(at), dims=2)
end

"""
    prevalence_filter(comm::AbstractAbundanceTable; minabundance=0.0; minprevalence=0.05, renorm=false)

Return a filtered `CommunityProfile` where features with prevalence lower than `minprevalence` are removed.
By default, a feature is considered "present" if > 0, but this can be changed by setting `minabundance`.

Optionally, set `renorm = true` to calculate relative abundances after low prevalence features are removed.
"""
function prevalence_filter(comm::AbstractAbundanceTable; minabundance=0.0, minprevalence=0.05, renorm=false)
    comm = comm[vec(prevalence(comm, minabundance) .>= minprevalence), :]
    return renorm ? relativeabundance(comm) : comm
end

"""
    rankfilter(comm::AbstractAbundanceTable, cl::Union{Symbol, Int}; keepempty=false)

Return a copy of `comm`, where only rows that have `taxrank(feature) == cl` are kept.
Use `keepempty = true` to also keep features that don't have a `rank` (eg "UNIDENTIFIED").
"""
function rankfilter(comm::AbstractAbundanceTable, cl::Symbol; keepempty=false)
    in(cl, keys(_ranks)) ||  error("Invalid rank $cl, must be one of $(keys(_ranks))")
    if keepempty
        return filter(f-> !hasrank(f) || taxrank(f) == cl, comm)
    else
        return filter(f-> hasrank(f) && taxrank(f) == cl, comm)
    end
end

function rankfilter(comm::AbstractAbundanceTable, rank::Int; keepempty=false)
    0 <= rank <= 9 ||  error("Invalid rank $rank, must be one of $_ranks")
    return rankfilter(comm, keys(_ranks)[rank+1]; keepempty)
end

## Metadata

"""
    metadata(commp::CommunityProfile)

Returns iterator of `NamedTuple` per sample, where keys are `:sample`
and each metadata key found in `commp`.
Samples without given metadata are filled with `missing`.

Returned values can be passed to any Tables.rowtable - compliant type,
eg `DataFrame`.
"""
function metadata(commp::CommunityProfile)
    ss = samples(commp)
    cols = unique(reduce(vcat, collect.(keys.(metadata.(samples(commp))))))
    return Tables.rowtable(merge((; sample=name(s)),
                     NamedTuple(c => get(s, c, missing) for c in cols)
                    ) for s in ss)
end

"""
    set!(commp::CommunityProfile, sample::AbstractString, prop::Symbol, val)
    set!(commp::CommunityProfile, sample::AbstractString, md::Union{AbstractDict, NamedTuple})

Update or insert a value `val` to the metadata of `sample` in the CommunityProfile `commp` using a Symbol `prop`.
If you want an error to be thrown if the value already exists, use [`insert!`](@ref).

Can also pass a Dictionary or NamedTuple containing key=> value pairs,
all of which will be `set!`.
"""
function set!(commp::CommunityProfile, sample::AbstractString, prop::Symbol, val)
    sample = samples(commp, sample)
    prop in _restricted_fields(sample) && error("Cannot set! $prop for $(typeof(sample)).")
    set!(sample.metadata, prop, val)
    return sample
end

function set!(commp::CommunityProfile, sample::AbstractString, md::Dictionary)
    for (key, value) in pairs(md)
        set!(commp, sample, key, value)
    end
    return nothing
end

function set!(commp::CommunityProfile, sample::AbstractString, md::Union{AbstractDict, NamedTuple})
    md = Dictionary(md)
    set!(commp, sample, md)
    return nothing
end

"""
    set!(cp::CommunityProfile, md; namecol=:sample)

Add metadata (in the form of a `Tables.jl` table) a `CommunityProfile`.
One column (`namecol`) should contain sample names that exist in `commp`,
and other columns should contain metadata that will be added to the metadata of each sample.
"""
function set!(commp::CommunityProfile, md; namecol=:sample)
    Tables.istable(md) || throw(ArgumentError("Metadata must be a Tables.table"))
    sns = Set(samplenames(commp))
    md = filter(row-> row[namecol] in sns, md)
    for row in Tables.rows(md)
        sample = samples(commp, row[namecol])
        ks = filter(!=(namecol), keys(first(md)))
        for k in ks
            set!(sample, k, row[k])
        end
    end
    return nothing
end

"""
    insert!(commp::CommunityProfile, sample::AbstractString, prop::Symbol, val)

Insert a value `val` to the metadata of `sample` in a CommunityProfile `commp` using a Symbol `prop`,
and it will throw an error if `prop` exists.
If you don't want an error to be thrown if the value exists, use [`set!`](@ref).
"""
function insert!(commp::CommunityProfile, sample::AbstractString, prop::Symbol, val)
    sample = samples(commp, sample)
    prop in _restricted_fields(sample) && error("Cannot insert! $prop for $(typeof(sample)).")
    insert!(sample.metadata, prop, val)
    return sample
end

function insert!(commp::CommunityProfile, sample::AbstractString, md::Dictionary)
    isempty(Set(keys(commp, sample)) ∩ Set(keys(md))) || throw(IndexError("Duplicate keys found. Use `set!` to overwrite"))
    for (key, value) in pairs(md)
        insert!(commp, sample, key, value)
    end
    return nothing
end

function insert!(commp::CommunityProfile, sample::AbstractString, md::Union{<:AbstractDict, NamedTuple})
    md = Dictionary(md)
    insert!(commp, sample, md)
    return nothing
end

"""
    insert!(cp::CommunityProfile, md; namecol=:sample)

Add metadata (in the form of a `Tables.jl` table) a `CommunityProfile`.
One column (`namecol`) should contain sample names that exist in `commp`,
and other columns should contain metadata that will be added to the metadata of each sample.

Before starting, this will check that every value in every row is `insert!`able,
and will throw an error if not.
This requires iterating over the metadata table twice, which may be slow.
If performance matters, you can use `set!` instead,
though this will overwrite existing data.
"""
function insert!(commp::CommunityProfile, md; namecol=:sample, careful=true)
    Tables.istable(md) || throw(ArgumentError("Metadata must be a Tables.table"))
    sns = Set(samplenames(commp))
    md = filter(row-> row[namecol] in sns, md)
    for row in Tables.rows(md)
        sample = samples(commp, row[namecol])
        ks = filter(!=(namecol), keys(first(md)))
        for k in ks
            haskey(sample, k) && throw(IndexError("Duplicate metadata detected. Use `set!` to force overwrite."))
        end
    end
    for row in Tables.rows(md)
        sample = samples(commp, row[namecol])
        ks = filter(!=(namecol), keys(first(md)))
        for k in ks
            set!(sample, k, row[k])
        end
    end
    return nothing
end

"""
unset!(commp::CommunityProfile, sample::AbstractString, prop::Symbol)

Delete a metadata entry in `sample` from CommunityProfile `commp` using the Symbol `prop`.
If you want an error to be thrown if the value does not exist, use [`delete!`](@ref).
"""
function unset!(commp::CommunityProfile, sample::AbstractString, prop::Symbol)
    sample = samples(commp, sample)
    prop in _restricted_fields(sample) && error("Cannot unset! $prop for $(typeof(sample)).")
    unset!(sample.metadata, prop)
    return sample
end

"""
    delete!(commp::CommunityProfile, sample::AbstractString, prop::Symbol)

Delete a metadata entry in `sample` from CommunityProfile `commp` using the Symbol `prop` if it exists, or throw an error otherwise.
If you don't want an error to be thrown if the value does not exist, use [`unset!`](@ref).
"""
function delete!(commp::CommunityProfile, sample::AbstractString, prop::Symbol)
    sample = samples(commp, sample)
    prop in _restricted_fields(sample) && error("Cannot delete! $prop for $(typeof(sample)).")
    delete!(sample.metadata, prop)
    return sample
end

"""
    keys(commp::CommunityProfile, sample::AbstractString)

Return an iterator over all keys of the metadata attached to `sample` in a CommunityProfile `commp`.
`collect(keys(commp, sample))` returns an array of keys.
"""
Base.keys(commp::CommunityProfile, sample::AbstractString) = keys(metadata(samples(commp, sample)))

"""
    haskey(commp::CommunityProfile, sample::AbstractString, key::Symbol)

Determine whether the metadata of `sample` in a CommunityProfile `commp` has a mapping for a given `key`.
Use `!haskey` to determine whether a `sample` in a CommunityProfile doesn't have a mapping for a given `key`
"""
Base.haskey(commp::CommunityProfile, sample::AbstractString, key::Symbol) = in(key, keys(samples(commp, sample)))

"""
    get(commp::CommunityProfile, sample::AbstractString, key::Symbol, default)

Return the value of the metadata in a `sample` stored for the given `key`, or the given `default` value if no mapping for the key is present.
"""
Base.get(commp::CommunityProfile, sample::AbstractString, key::Symbol, default) = get(metadata(samples(commp, sample)), key, default)

"""
    filter(f, comm::CommunityProfile)

Apply `f` to the features of `comm`,
and return a copy where `f(feature)` is `true`.
"""
function Base.filter(f::Function, commp::CommunityProfile)
    ridx = findall(f, features(commp))
    isempty(ridx) && error("Can't return empty profile")
    return copy(commp[ridx, :])
end


text:
include("utils.jl")

@testset "Rotations" begin
    M = Manifolds.Rotations(2)
    @test repr(M) == "Rotations(2)"
    @test representation_size(M) == (2, 2)
    @test injectivity_radius(M) == π * sqrt(2.0)
    @test injectivity_radius(M, [1.0 0.0; 0.0 1.0]) == π * sqrt(2.0)
    @test injectivity_radius(M, ExponentialRetraction()) == π * sqrt(2.0)
    @test injectivity_radius(M, [1.0 0.0; 0.0 1.0], ExponentialRetraction()) ==
          π * sqrt(2.0)
    @test injectivity_radius(M, PolarRetraction()) ≈ π / sqrt(2)
    @test injectivity_radius(M, [1.0 0.0; 0.0 1.0], PolarRetraction()) ≈ π / sqrt(2)
    types = [Matrix{Float64}]
    TEST_FLOAT32 && push!(types, Matrix{Float32})
    TEST_STATIC_SIZED && push!(types, MMatrix{2,2,Float64,4})
    retraction_methods = [Manifolds.PolarRetraction(), Manifolds.QRRetraction()]

    inverse_retraction_methods =
        [Manifolds.PolarInverseRetraction(), Manifolds.QRInverseRetraction()]

    basis_types = (DefaultOrthonormalBasis(), ProjectedOrthonormalBasis(:svd))

    @testset "vee/hat" begin
        M = Manifolds.Rotations(2)
        v = [1.23]
        x = Matrix{Float64}(I, 2, 2)
        V = Manifolds.hat(M, x, v)
        @test isa(V, AbstractMatrix)
        @test norm(M, x, V) / sqrt(2) ≈ norm(v)
        @test Manifolds.vee(M, x, V) == v

        V = project(M, x, randn(2, 2))
        v = Manifolds.vee(M, x, V)
        @test isa(v, AbstractVector)
        @test Manifolds.hat(M, x, v) == V
    end

    for T in types
        angles = (0.0, π / 2, 2π / 3, π / 4)
        pts = [convert(T, [cos(ϕ) -sin(ϕ); sin(ϕ) cos(ϕ)]) for ϕ in angles]
        test_manifold(
            M,
            pts;
            test_reverse_diff=false,
            test_injectivity_radius=false,
            test_project_tangent=true,
            test_musical_isomorphisms=true,
            retraction_methods=retraction_methods,
            inverse_retraction_methods=inverse_retraction_methods,
            point_distributions=[Manifolds.normal_rotation_distribution(M, pts[1], 1.0)],
            tvector_distributions=[Manifolds.normal_tvector_distribution(M, pts[1], 1.0)],
            basis_types_to_from=basis_types,
        )

        @testset "log edge cases" begin
            v = Manifolds.hat(M, pts[1], [Float64(π)])
            x = exp(M, pts[1], v)
            @test isapprox(x, exp(M, pts[1], log(M, pts[1], x)))
        end

        v = log(M, pts[1], pts[2])
        @test norm(M, pts[1], v) ≈ (angles[2] - angles[1]) * sqrt(2)

        v14_polar = inverse_retract(M, pts[1], pts[4], Manifolds.PolarInverseRetraction())
        p4_polar = retract(M, pts[1], v14_polar, Manifolds.PolarRetraction())
        @test isapprox(M, pts[4], p4_polar)

        v14_qr = inverse_retract(M, pts[1], pts[4], Manifolds.QRInverseRetraction())
        p4_qr = retract(M, pts[1], v14_qr, Manifolds.QRRetraction())
        @test isapprox(M, pts[4], p4_qr)
    end

    @testset "Distribution tests" begin
        usd_mmatrix =
            Manifolds.normal_rotation_distribution(M, (@MMatrix [1.0 0.0; 0.0 1.0]), 1.0)
        @test isa(rand(usd_mmatrix), MMatrix)

        usd1_mmatrix =
            Manifolds.normal_rotation_distribution(Rotations(1), (@MMatrix [1.0]), 1.0)
        @test isa(rand(usd1_mmatrix), MMatrix)
        @test rand(usd1_mmatrix) == @MMatrix [1.0]

        gtsd_mvector =
            Manifolds.normal_tvector_distribution(M, (@MMatrix [1.0 0.0; 0.0 1.0]), 1.0)
        @test isa(rand(gtsd_mvector), MMatrix)
    end

    Random.seed!(42)
    for n in (3, 4, 5)
        @testset "Rotations: SO($n)" begin
            SOn = Manifolds.Rotations(n)
            ptd = Manifolds.normal_rotation_distribution(SOn, Matrix(1.0I, n, n), 1.0)
            tvd = Manifolds.normal_tvector_distribution(SOn, Matrix(1.0I, n, n), 1.0)
            pts = [rand(ptd) for _ in 1:3]
            test_manifold(
                SOn,
                pts;
                test_forward_diff=n == 3,
                test_reverse_diff=false,
                test_injectivity_radius=false,
                test_musical_isomorphisms=true,
                test_mutating_rand=true,
                retraction_methods=retraction_methods,
                inverse_retraction_methods=inverse_retraction_methods,
                point_distributions=[ptd],
                tvector_distributions=[tvd],
                basis_types_to_from=basis_types,
                exp_log_atol_multiplier=20,
                retraction_atol_multiplier=12,
            )

            @testset "vee/hat" begin
                x = Matrix(1.0I, n, n)
                v = randn(manifold_dimension(SOn))
                V = Manifolds.hat(SOn, x, v)
                @test isa(V, AbstractMatrix)
                @test norm(SOn, x, V) / sqrt(2) ≈ norm(v)
                @test Manifolds.vee(SOn, x, V) == v

                V = project(SOn, x, randn(n, n))
                v = Manifolds.vee(SOn, x, V)
                @test isa(v, AbstractVector)
                @test Manifolds.hat(SOn, x, v) ≈ V
            end

            if n == 4
                @testset "exp/log edge cases" begin
                    vs = [
                        [0, 0, π, 0, 0, π],  # θ = (π, π)
                        [0, 0, π, 0, 0, 0],  # θ = (π, 0)
                        [0, 0, π / 2, 0, 0, π],  # θ = (π, π/2)
                        [0, 0, π, 0, 0, 0] ./ 2,  # θ = (π/2, 0)
                        [0, 0, π, 0, 0, π] ./ 2,  # θ = (π/2, π/2)
                        [0, 0, 0, 0, 0, 0],  # θ = (0, 0)
                        [0, 0, 1, 0, 0, 1] .* 1e-100, # α = β ≈ 0
                        [0, 0, 1, 0, 0, 1] .* 1e-6, # α = β ⩰ 0
                        [0, 0, 10, 0, 0, 1] .* 1e-6, # α ⪆ β ⩰ 0
                        [0, 0, π / 4, 0, 0, π / 4 - 1e-6], # α ⪆ β > 0
                    ]
                    for v in vs
                        @testset "rotation vector $v" begin
                            V = Manifolds.hat(SOn, Matrix(1.0I, n, n), v)
                            x = exp(V)
                            @test x ≈ exp(SOn, one(x), V)
                            @test ForwardDiff.derivative(t -> exp(SOn, one(x), t * V), 0) ≈
                                  V
                            x2 = exp(log(SOn, one(x), x))
                            @test isapprox(x, x2; atol=1e-6)
                        end
                    end
                end
            end

            v = Matrix(
                Manifolds.hat(SOn, pts[1], π * normalize(randn(manifold_dimension(SOn)))),
            )
            x = exp(SOn, pts[1], v)
            v2 = log(SOn, pts[1], x)
            @test x ≈ exp(SOn, pts[1], v2)
        end
    end
    @testset "Test Manifold Point and Tangent Vector checks" begin
        M = Manifolds.Rotations(2)
        for x in [1, [2.0 0.0; 0.0 1.0], [1.0 0.5; 0.0 1.0]]
            @test_throws DomainError is_manifold_point(M, x, true)
            @test !is_manifold_point(M, x)
        end
        x = one(zeros(2, 2))
        @test is_manifold_point(M, x)
        @test is_manifold_point(M, x, true)
        for v in [1, [0.0 1.0; 0.0 0.0]]
            @test_throws DomainError is_tangent_vector(M, x, v, true)
            @test !is_tangent_vector(M, x, v)
        end
        v = [0.0 1.0; -1.0 0.0]
        @test is_tangent_vector(M, x, v)
        @test is_tangent_vector(M, x, v, true)
    end
    @testset "Project point" begin
        M = Manifolds.Rotations(2)
        x = Matrix{Float64}(I, 2, 2)
        x1 = project(M, x)
        @test is_manifold_point(M, x1, true)

        M = Manifolds.Rotations(3)
        x = collect(reshape(1.0:9.0, (3, 3)))
        x2 = project(M, x)
        @test is_manifold_point(M, x2, true)

        rng = MersenneTwister(44)
        x3 = project(M, randn(rng, 3, 3))
        @test is_manifold_point(M, x3, true)
    end
    @testset "Edge cases of Rotations" begin
        @test_throws OutOfInjectivityRadiusError inverse_retract(
            Rotations(2),
            [1.0 0.0; 0.0 1.0],
            [1.0 0.0; 0.0 -1.0],
            PolarInverseRetraction(),
        )
    end
end


text:
module QuantumStateBase
    using LinearAlgebra
    using ClassicalOrthogonalPolynomials
    using OhMyArtifacts

    const my_artifacts = Ref{String}()

    function __init__()
        my_artifacts[] = @my_artifacts_toml!()
        return
    end

    # polynomial
    include("polynomial.jl")

    # state
    include("operator.jl")
    include("state.jl")

    # wigner
    include("wigner_util.jl")
    include("wigner.jl")
end


text:
module ITensorParallel

using Distributed
using MPI
using ITensors
using ITensors.NDTensors

import ITensors: product, position!, noiseterm, lproj, rproj

include("partition.jl")
include("threaded_projmposum.jl")
include("distributed_projmposum.jl")
include("mpi_projmposum.jl")

export ThreadedProjMPOSum, partition, DistributedSum, MPISum

end


text:
## Handler for reading outputs of crowdsource processing on DECaPS
module decam

export inject_rename #
export read_decam #
export read_crowdsource #
export load_psfmodel_cs #
export save_fxn #
export get_catnames #

export proc_ccd
export proc_all

export decam_dir

using CloudCovErr
using PyCall
using FITSIO
import ImageFiltering
import Distributions
import StatsBase
using Random
using LinearAlgebra
import Conda

"""
    __int__()

Builds the required python crowdsource dependency and exports the required load
function for obtaining the position dependent psf to the python namespace.
"""
function __init__()
    if !haskey(Conda._installed_packages_dict(),"crowdsourcephoto")
        Conda.add("crowdsourcephoto",channel="conda-forge")
    end
    decam_dir = dirname(@__FILE__)*"/decam_dir"
    py"""
    import sys
    import crowdsource.psf as psfmod
    from crowdsource.decam_proc import read_data
    from astropy.io import fits
    from functools import partial
    import os

    # default decam_dir at Harvard
    if 'DECAM_DIR' not in os.environ:
        os.environ['DECAM_DIR'] = $decam_dir

    def load_psfmodel(outfn, ccd, filter, pixsz=9):
        f = fits.open(outfn)
        psfmodel = psfmod.linear_static_wing_from_record(f[ccd+"_PSF"].data[0],filter=filter)
        f.close()
        psfmodel.fitfun = partial(psfmod.fit_linear_static_wing, filter=filter, pixsz=pixsz)
        return psfmodel
    """
end

## File specific read/write functions

"""
    inject_rename(fname) -> ifname

Convenience renaming file paths to read data from injection tests
which are stored separately from data obtained on DECam for data
provenance purposes.

# Arguments:
- `fname`: file name for exposure data from DECam

# Output:
- `ifname`: corresponding file name for injection tests into that exposure
"""
function inject_rename(fname)
    splitname = split(fname,"/")
    splitname[end-1] = "decapsi"
    return chop(join(splitname,"/"),tail=7)*"I.fits.fz"
end

"""
    read_decam(base,date,filt,vers,ccd;corrects7=true) -> ref_im, d_im

Read in raw image files associated with exposures obtain on the DarkEnergyCamera.
Returns the image and a quality flag mask image. See [NOAO
handbook](http://ast.noao.edu/sites/default/files/NOAO_DHB_v2.2.pdf) for more details
on what is contained in each file and how they are obtained.

# Arguments:
- `base`: parent directory and file name prefix for exposure files
- `date`: date_time of the exposure
- `filt`: optical filter used to take the exposure
- `vers`: NOAO community processing version number
- `ccd`: which ccd we are pulling the image for

# Keywords:
- `corrects7`: use `crowdsource` load to read ccd "S7" to correct for floating amplifier on half of the chip (default true)

# Output:
- `ref_im`: image of photoelectron counts from observation on DECam
- `d_im`: quality flag mask image from NOAO community pipeline

# Example
```julia
ref_im, d_im = read_decam("/n/fink2/decaps/c4d_","170420_040428","g","v1","N14")
```
"""
function read_decam(base,date,filt,vers,ccd;corrects7=true)
    ifn = base*date*"_ooi_"*filt*"_"*vers*".fits.fz"
    wfn = base*date*"_oow_"*filt*"_"*vers*".fits.fz"
    dfn = base*date*"_ood_"*filt*"_"*vers*".fits.fz"
    if last(ccd,1) == "I"
        ifn = inject_rename(ifn)
        wfn = inject_rename(wfn)
        dfn = inject_rename(dfn)
    end
    if corrects7 .& ((ccd == "S7") .| (ccd == "S7I"))
        # a little wasteful to throw away load times in python for w_im and d_im
        # but we need the crowdsource formatting for s7 correction and it is easiest
        # to keep the CloudCovErr formatting consistent
        py_ref_im, py_w_im, py_d_im, nebprob = py"read_data"(ifn,wfn,dfn,ccd,maskdiffuse=false)
        py_w_im = nothing
        py_d_im = nothing
        nebprob = nothing
        ref_im = py_ref_im'
    else
        f = FITS(ifn)
        ref_im = read(f[ccd])
        close(f)
    end
    f = FITS(dfn)
    d_im = read(f[ccd])
    close(f)
    return ref_im, d_im
end

"""
    read_crowdsource(basecat,date,filt,vers,ccd) -> x_stars, y_stars, flux_stars, decapsid, gain, mod_im, sky_im, wcol, w

Read in outputs of crowdsource, a photometric pipeline. To pair with an arbitrary
photometric pipeline, an analogous read in function should be created. The relevant
outputs are the model image (including the sources) so that we can produce the
residual image, the sky/background model (no sources), and the coordinates of the stars.
The survey id number is also readout of the pipeline solution file to help
cross-validate matching of the CloudCovErr outputs and the original sources. The empirical
gain is read out of the header (for other photometric pipelines which don't perform this estiamte,
the gain from DECam is likely sufficient). All columns from the photometric catalogue are
also read in at this point to be rexported with the CloudCovErr outputs.

# Arguments:
- `basecat`: parent directory of the cat directory holding all of the single-epoch crowdsource catalogue files
- `date`: date_time of the exposure
- `filt`: optical filter used to take the exposure
- `vers`: NOAO community processing version number
- `ccd`: which ccd we are pulling the image for

# Keywords:
- `corrects7`: use `crowdsource` load to read ccd "S7" to correct for floating amplifier on half of the chip (default true)

# Output:
- `x_stars`: list of source x-coordinates (accounting for indexing order and start point)
- `y_stars`: list of source y-coordinates (accounting for indexing order and start point)
- `flux_stars`: list of stellar fluxes in ADU
- `decapsid`: list of survey id number for each detection
- `gain`: gain of detector to convert from photon count noise to detector noise
- `mod_im`: model image (including the sources) from photometric pipeline
- `sky_im`: sky image (background estimate) from photometric pipeline
- `wcol`: list of all column names in photometric catalogue
- `w`: list of all column values in photometric catalogue
"""
function read_crowdsource(basecat,date,filt,vers,ccd)
    f = FITS(basecat*"cat/c4d_"*date*"_ooi_"*filt*"_"*vers*".cat.fits")
    x_stars = read(f[ccd*"_CAT"],"x")
    y_stars = read(f[ccd*"_CAT"],"y")
    flux_stars = read(f[ccd*"_CAT"],"flux")
    decapsid = read(f[ccd*"_CAT"],"decapsid")
    gain = read_key(f[ccd*"_HDR"],"GAINCRWD")[1]
    wcol = String[]
    w = []
    for col in FITSIO.colnames(f[ccd*"_CAT"])
        push!(wcol,col)
        push!(w,read(f[ccd*"_CAT"],col))
    end
    close(f)

    if gain < 0
        println("WARNING: Gain is negative! Do not trust the input from this exposure!")
        gain = abs(gain)
    end

    f = FITS(basecat*"mod/c4d_"*date*"_ooi_"*filt*"_"*vers*".mod.fits")
    mod_im = read(f[ccd*"_MOD"])
    sky_im = read(f[ccd*"_SKY"])
    close(f)
    #switch x, y order and compensate for 0 v 1 indexing between Julia and python
    return y_stars.+1, x_stars.+1, flux_stars, decapsid, gain, mod_im, sky_im, wcol, w
end

"""
    load_psfmodel_cs(base,date,filt,vers,ccd) -> psfmodel

Julia wrapper function for the PyCall that reads the position dependent psfmodel
produced by crowdsource from the catalogue file for a given exposure and ccd. The
returned psfmodel takes an x- and y-position for the source location and the size
of the desired psfstamp (the stamps are square and required to be odd).

# Arguments:
- `base`: parent directory and file name prefix for catalogue files
- `date`: date_time of the exposure
- `filt`: optical filter used to take the exposure
- `vers`: NOAO community processing version number
- `ccd`: which ccd we are pulling the image for

# Output:
- `psfmodel`: function that returns PSF stamp from parametric PSF model that is a function of position
"""
function load_psfmodel_cs(base,date,filt,vers,ccd)
    psfmodel_py = py"load_psfmodel"(base*"cat/c4d_"*date*"_ooi_"*filt*"_"*vers*".cat.fits",ccd,filt)
    function psfmodel_jl(x,y,sz)
        # accounts for x, y ordering and 0 v 1 indexing between python and Julia
        return psfmodel_py(y.-1,x.-1,sz)'
    end
    return psfmodel_jl
end

"""
    save_fxn(wcol,w,basecat,date,filt,vers,ccd)

Saves CloudCovErr.jl outputs and initial photometric catalogue outputs to a new
single-epoch catalogue. Massages types of columns to reduce data storage size.
Converts the native CloudCovErr.jl output of the bias offset value into a `cflux`
corrected flux column for the ease of catalogue users.

# Arguments:
- `wcol`: list of all column names in photometric catalogue
- `w`: list of all column values in photometric catalogue
- `basecat`: parent directory of the cat directory holding all of the single-epoch crowdsource catalogue files
- `date`: date_time of the exposure
- `filt`: optical filter used to take the exposure
- `vers`: NOAO community processing version number
- `ccd`: which ccd we are pulling the image for
"""
function save_fxn(wcol,w,basecat,date,filt,vers,ccd)
    for i=1:length(wcol)
        if (wcol[i] == "passno") | (wcol[i] == "dnt")
            w[i] .= convert.(Int8,w[i])
        elseif (wcol[i] == "kcond0") | (wcol[i] == "kcond") | (wcol[i] == "kpred")
            w[i] .= convert.(Int32,w[i])
        elseif (wcol[i] == "rchi2") | (wcol[i] == "spread_model") | (wcol[i] == "dspread_model")
            w[i] .= convert.(Float32,w[i])
        elseif wcol[i] == "fdb_tot"
            for j=1:length(wcol)
                if wcol[j] == "flux"
                    w[i].+=w[j]
                    wcol[i] = "cflux"
                end
            end
        end
    end
    f = FITS(basecat*"cer/c4d_"*date*"_ooi_"*filt*"_"*vers*".cat.cer.fits","r+")
    write(f,wcol,w,name=ccd*"_CAT")
    close(f)
end

"""
    get_catnames(f) -> extnames

Reads list of extension names from an open FITS file to determine which
CCDs have completed photometric catalogues and are eligible for CloudCovErr.jl.

# Arguments:
- `f`: an open FITS file handle containing `crowdsource` catalogues

# Outputs:
- `extnames`: list of CCDs that have photometric catalogues in f

"""
function get_catnames(f)
    nhdu = length(f)
    extnames = String[]
    for i = 1:nhdu-1
        FITSIO.fits_movabs_hdu(f.fitsfile, i+1)
        extname = something(FITSIO.fits_try_read_extname(f.fitsfile), "")
        if last(extname,3) == "CAT"
            push!(extnames,chop(extname,tail=4))
        end
    end
    return extnames
end

## Main run functions for processing exposures
# combines all of the functions in the repo for
# an acutal implementation

"""
    proc_ccd(base,date,filt,vers,basecat,ccd;thr=20,outthr=20000,Np=33,corrects7=true,widx=129,widy=widx,tilex=1,tiley=tilex,ftype::Int=32)

Primary run function for a given CCD image of a larger exposure.

# Arguments:
- `base`: parent directory and file name prefix for exposure files
- `date`: date_time of the exposure
- `filt`: optical filter used to take the exposure
- `vers`: NOAO community processing version number
- `basecat`: parent directory of the cat directory holding all of the single-epoch crowdsource catalogue files
- `ccd`: which ccd we are pulling the image for

# Keywords:
- `thr`: threshold used for flux-dependent masking (default 20)
- `outthr`: threshold for residual-based masking (default 20000)
- `Np`: size of local covariance matrix in pixels (default 33)
- `corrects7`: use `crowdsource` load to read ccd "S7" to correct for floating amplifier on half of the chip (default true)
- `widx`: width of boxcar window in x which determines size of region used for samples for the local covariance estimate (default 129)
- `widy`: width of boxcar window in y which determines size of region used for samples for the local covariance estimate (default 129)
- `tilex`: total number of tile divisions along x (default 1)
- `tiley`: total number of tile divisions along y (default tilex)
- `ftype::Int`: determine the Float precision, 32 is Float32, otherwise Float64
"""
function proc_ccd(base,date,filt,vers,basecat,ccd;thr=20,outthr=20000,Np=33,corrects7=true,widx=129,widy=widx,tilex=1,tiley=tilex,ftype::Int=32)
    println("Started $ccd")
    flush(stdout)

    if ftype == 32
        T = Float32
    else
        T = Float64
    end

    # loads from disk
    ref_im, d_im = CloudCovErr.read_decam(base,date,filt,vers,ccd,corrects7=corrects7)
    bmaskd = (d_im .!= 0)
    d_im = nothing
    (sx0, sy0) = size(ref_im)
    x_stars, y_stars, flux_stars, decapsid, gain, mod_im, sky_im, wcol, w = CloudCovErr.read_crowdsource(basecat,date,filt,vers,ccd)
    (Nstars,) = size(x_stars)

    if Nstars > 0
        psfmodel = CloudCovErr.load_psfmodel_cs(basecat,date,filt,vers,ccd)
        psfstatic511 = psfmodel(sx0÷2,sy0÷2,511)
        psfstatic33 = psfmodel(sx0÷2,sy0÷2,Np)

        # mask bad camera pixels/cosmic rays, then mask out star centers
        CloudCovErr.gen_mask_staticPSF2!(bmaskd, psfstatic511, psfstatic33, x_stars, y_stars, flux_stars; thr=thr)

        testim = mod_im .- ref_im
        bimage = zeros(T,sx0,sy0)
        bimageI = zeros(Int64,sx0,sy0)
        testim2 = zeros(T,sx0,sy0)
        bmaskim2 = zeros(Bool,sx0,sy0)
        goodpix = zeros(Bool,sx0,sy0)

        bmaskd .|= (abs.(testim) .> outthr)
        prelim_infill!(testim,bmaskd,bimage,bimageI,testim2,bmaskim2,goodpix,ccd;widx=19,widy=19,ftype=ftype)
        testim .= mod_im .- ref_im #fixes current overwrite for 0 infilling
        ref_im = nothing
        bimage = nothing
        bimageI = nothing
        bmaskim2 = nothing
        goodpix = nothing

        ## calculate the star farthest outside the edge of the image in x and y
        cx = round.(Int,x_stars)
        cy = round.(Int,y_stars)
        px0 = outest_bounds(cx,sx0)
        py0 = outest_bounds(cy,sy0)

        ## these have to be allocating to get the noise model right
        Δx = (widx-1)÷2
        Δy = (widy-1)÷2
        padx = Np+Δx+px0
        pady = Np+Δy+py0
        in_image = ImageFiltering.padarray(testim2,ImageFiltering.Pad(:reflect,(padx+2,pady+2)));
        in_image_raw = ImageFiltering.padarray(testim,ImageFiltering.Pad(:reflect,(padx+2,pady+2)));
        in_sky_im = ImageFiltering.padarray(sky_im,ImageFiltering.Pad(:reflect,(padx+2,pady+2)));
        in_stars_im = ImageFiltering.padarray(abs.(mod_im.-sky_im)./gain,ImageFiltering.Pad(:reflect,(padx+2,pady+2)));
        in_bmaskd = ImageFiltering.padarray(bmaskd,ImageFiltering.Fill(true,(padx+2,pady+2)));

        # exposure datetime based seed
        rndseed = parse(Int,date[1:6])*10^6 + parse(Int,date[8:end])
        CloudCovErr.add_sky_noise!(in_image,in_bmaskd,in_sky_im,gain;seed=rndseed)

        ## iterate over all star positions and compute errorbars/debiasing corrections
        star_stats = zeros(T,10,Nstars)
        star_k = zeros(Int32,10,Nstars)

        # preallocate the cov and μ per star variables
        cov = zeros(T,Np*Np,Np*Np)
        μ = zeros(T,Np*Np)

        # compute a radial mask for reduced num cond pixels
        circmask = kstar_circle_mask(Np,rlim=256)

        # some important global sizes for the loop
        cntStar0 = 0
        stepx = (sx0+2) ÷ tilex
        stepy = (sy0+2) ÷ tiley

        # precallocate the image subblocks
        #GC.gc(false)
        in_subimage = zeros(T,stepx+2*padx,stepy+2*pady)
        ism = zeros(T,stepx+2*padx,stepy+2*pady)
        bimage = zeros(T,stepx+2*padx-2*Δx,stepy+2*pady-2*Δy)
        bism = zeros(T,stepx+2*padx-2*Δx,stepy+2*pady-2*Δy,2*Np-1, Np);
        for jx=1:tilex, jy=1:tiley
            xrng, yrng, star_ind = im_subrng(jx,jy,cx,cy,sx0+2,sy0+2,px0,py0,stepx,stepy,padx,pady,tilex,tiley)
            cntStar = length(star_ind)
            if cntStar > 0
                in_subimage .= in_image[xrng,yrng]
                cov_avg!(bimage, ism, bism, in_subimage, widx=widx, widy=widy)
                offx = padx-Δx-(jx-1)*stepx
                offy = pady-Δy-(jy-1)*stepy
                for i in star_ind
                    build_cov!(cov,μ,cx[i]+offx,cy[i]+offy,bimage,bism,Np,widx,widy)
                    data_in, stars_in, kmasked2d = stamp_cutter(cx[i],cy[i],in_image_raw,in_stars_im,in_bmaskd;Np=Np)
                    psft, kstar, kpsf2d, kcond0, kcond, kpred, dnt = gen_pix_mask(kmasked2d,psfmodel,circmask,x_stars[i],y_stars[i],flux_stars[i];Np=Np,thr=thr)
                    try
                        star_stats[:,i] .= [condCovEst_wdiag(cov,μ,kstar,kpsf2d,data_in,stars_in,psft,diag_on=true)[1]..., kcond0, kcond, kpred, dnt]
                    catch
                        star_stats[:,i] .= [NaN, NaN, NaN, NaN, NaN, NaN, kcond0, kcond, kpred, dnt]
                    end
                end
            end
            cntStar0 += cntStar
            println("Finished $cntStar stars in tile ($jx, $jy) of $ccd")
            flush(stdout)
        end
        # prepare for export by appending to cat vectors
        # if doing the ops float64, might want to do a final convert to float32 before saving
        for (ind,col) in enumerate(["dcflux","dcflux_diag","fdb_tot","fdb_res","fdb_pred","cchi2","kcond0","kcond","kpred","dnt"])
            push!(wcol,col)
            push!(w,star_stats[ind,:])
        end
        CloudCovErr.save_fxn(wcol,w,basecat,date,filt,vers,ccd)
        println("Saved $ccd processing $cntStar0 of $Nstars stars")
        flush(stdout)
        pdefer = count(isnan.(star_stats[1,:]))
        if (pdefer > 0)
            println("There were posDef errors: $pdefer")
            flush(stdout)
        end

        ## manually call out memory as needed to be collected
        x_stars=nothing
        y_stars=nothing
        flux_stars = nothing
        cx = nothing
        cy = nothing
        wcol = nothing
        w = nothing
        star_stats = nothing
    else
        # prepare for export by appending to cat vectors
        # if doing the ops float64, might want to do a final convert to float32 before saving
        for (ind,col) in enumerate(["dcflux","dcflux_diag","fdb_tot","fdb_res","fdb_pred","cchi2","kcond0","kcond","kpred","dnt"])
            push!(wcol,col)
            push!(w,Float32[])
        end
        CloudCovErr.save_fxn(wcol,w,basecat,date,filt,vers,ccd)
        println("Saved $ccd processing 0 of $Nstars stars")
        flush(stdout)
    end
    return
end

"""
    proc_all(base,date,filt,vers,basecat;ccdlist=String[],resume=false,corrects7=true,thr=20,outthr=20000,Np=33,widx=129,widy=widx,tilex=1,tiley=tilex,ftype::Int=32)

Exposure level run function the manages which ccds to run and calls proc_ccd serially.

# Arguments:
- `base`: parent directory and file name prefix for exposure files
- `date`: date_time of the exposure
- `filt`: optical filter used to take the exposure
- `vers`: NOAO community processing version number
- `basecat`: parent directory of the cat directory holding all of the single-epoch crowdsource catalogue files

# Keywords:
- `ccdlist`: run only ccds in this list
- `resume`: if the exposure is partially complete, resume running from where it left off (default false)
- `corrects7`: use `crowdsource` load to read ccd "S7" to correct for floating amplifier on half of the chip (default true)
- `thr`: threshold used for flux-dependent masking (default 20)
- `outthr`: threshold for residual-based masking (default 20000)
- `Np`: size of local covariance matrix in pixels (default 33)
- `widx`: width of boxcar window in x which determines size of region used for samples for the local covariance estimate (default 129)
- `widy`: width of boxcar window in y which determines size of region used for samples for the local covariance estimate (default 129)
- `tilex`: total number of tile divisions along x (default 1)
- `tiley`: total number of tile divisions along y (default tilex)
- `ftype::Int`: determine the Float precision, 32 is Float32, otherwise Float64
"""
function proc_all(base,date,filt,vers,basecat;ccdlist=String[],resume=false,corrects7=true,thr=20,outthr=20000,Np=33,widx=129,widy=widx,tilex=1,tiley=tilex,ftype::Int=32)
    infn = basecat*"cat/c4d_"*date*"_ooi_"*filt*"_"*vers*".cat.fits"
    outfn = basecat*"cer/c4d_"*date*"_ooi_"*filt*"_"*vers*".cat.cer.fits"
    println("Starting to process "*infn)

    f = FITS(infn)
    extnames=get_catnames(f)

    # write output file or read it in and check which ccds are completed
    if ((!resume) | (!isfile(outfn)))
        prihdr = read_header(f[1])
        close(f)
        f = FITS(outfn,"w")
        write(f,[0], header=prihdr)
        close(f)
        extnamesdone = String[]
    else
        close(f)
        f = FITS(outfn,"r+")
        extnamesdone=get_catnames(f)
        close(f)
    end

    ## prepare list of ccds to run
    # remove the ones already done
    if length(extnamesdone) !=0
        alreadydone = intersect(extnames,extnamesdone)
        extnames = setdiff(extnames,extnamesdone)
        println("Skipping ccds already completed: ", alreadydone)
    end
    # limit to the ccds in the list not yet run
    if length(ccdlist) != 0
        extnames = intersect(extnames,ccdlist)
        println("Only running unfinished ccds in ccdlist: ", extnames)
    end

    # main loop over ccds
    for ccd in extnames
        proc_ccd(base,date,filt,vers,basecat,ccd;thr=thr,outthr=outthr,Np=Np,corrects7=corrects7,widx=widx,widy=widx,tilex=tilex,tiley=tilex,ftype=ftype)
        GC.gc()
    end
end

end