text:
module CTMLEs

using Logging
if isdefined(Main, :CTMLE_LOG_LEVEL)
    @eval @Logging.configure(level=$(Main.CTMLE_LOG_LEVEL))
else
    @Logging.configure(level=WARNING)
end

export ctmle,

SearchStrategy,
ForwardStepwise,
PreOrdered,
OrderingStrategy,
LogisticOrdering,
PartialCorrOrdering,
HDPSOrdering,

StratifiedKfold,
StratifiedRandomSub,

flucinfo

using StatsBase, MLBase
using ..LReg, ..Common, ..Parameters, ..Qmodels

import ..Qmodels: fluctuate!, defluctuate!, Fluctuation

include("pcor.jl")
include("strategies.jl")

subset{P<:Parameter}(param::P, idx::AbstractVector{Int}) =
    P(map(d -> subset(d, idx), regimens(param))...)

subset(d::StaticRegimen, idx) = d
subset(d::DynamicRegimen, idx) = DynamicRegimen(d.a[idx])

"""
A `Qchunk` holds (a subset of) the data set and corresponding indexes along with a `Qmodel` which can compute predictions for
each observation, a `Parameter` (whos definition may depend on each observation) and the value of the emperical risk
for the current `Qmodel`.
"""
type Qchunk{T<:AbstractFloat}
    q::Qmodel{T}
    W::Matrix{T}
    A::Vector{T}
    Y::Vector{T}
    param::Parameter{T}
    idx::AbstractVector{Int}
    gbounds::Vector{T}
    penalize::Bool
    risk::T
end

function make_chunk_subset{T<:AbstractFloat}(logitQnA1::Vector{T}, logitQnA0::Vector{T},
                           W::Matrix{T}, A::Vector{T}, Y::Vector{T},
                           param::Parameter{T}, idx::AbstractVector{Int}, gbounds::Vector{T},
                           penalize::Bool)
    length(logitQnA1) == length(logitQnA0) == size(W, 1) == length(A) == length(Y) ||
        error(ArgumentError("Input sizes do not match"))
    q = Qmodel(logitQnA1[idx], logitQnA0[idx])
    param = subset(param, idx)
    W = W[idx, :]
    A = A[idx]
    Y = Y[idx]
    #Can't compute penalized risk if q isn't fluctuated yet,
    #so set to Inf
    q_risk = penalize? convert(T, Inf) : risk(q, A, Y, param, false)
    Qchunk(q, W, A, Y, param, idx, gbounds, penalize, q_risk)
end

StatsBase.nobs(chunk::Qchunk) = nobs(chunk.q)

"""
A `QCV` holds a training `Qchunk` and a validation `Qchunk`, a sequence of estimated gs corresponding to
each fluctuation of the `Qchunk`s, the `SearchStrategy` in use, and indexes of used and unused covariates.
"""
type QCV{T<:AbstractFloat}
    chunk_train::Qchunk{T}
    chunk_val::Qchunk{T}
    gseq::Vector{LR{T}}
    searchstrategy::SearchStrategy
    used_covars::IntSet
    unused_covars::IntSet
end

"""
Removes the last fluctuation from a `QCV`

** Arguments **

* `qcv` - A `QCV`

Returns last estimated g and last fluctuation from the training and validation chunks.
"""
function defluctuate!(qcv::QCV)
    (pop!(qcv.gseq),
     defluctuate!(qcv.chunk_train.q),
     defluctuate!(qcv.chunk_val.q))
end

"""
Adds a previously calculated fluctuation to a `QCV`

** Arguments **

* `qcv` - A `QCV`
* `g` - A `LR`, estimate of g
* `fluc_train` - fluctuation for the training set
* `fluc_val` - fluctuation for the validation set

Returns training and validation fluctuations.
"""
function fluctuate!(qcv::QCV, g::LR, fluc_train::Fluctuation, fluc_val::Fluctuation)
    push!(qcv.gseq, g)
    fluctuate!(qcv.chunk_train, fluc_train)
    fluctuate!(qcv.chunk_val, fluc_val)
    (fluc_train, fluc_val)
end

"""
Adds a fluctuation to a `QCV` given an estimated g

** Arguments **

* `qcv` - A `QCV`
* `g` - A `LR`, estimate of g

** Details **

A new fluctuation for the training chunk is calculated along with a corresponding flucutation for the validation
chunk, based on the estiamted `g`, then the `QCV` is updated.

Returns training and validation fluctuations.
"""
function fluctuate!(qcv, g::LR)
    train = qcv.train_chunk
    val = qcv.val_chunk
    gn1_train = bound!(predict(g, train.W), train.gbounds)
    gn1_val = bound!(predict(g, val.W), val.gbounds)
    fluc_train = computefluc(train.q, train.param, gn1_train, train.A, train.Y)
    fluc_val = valfluc(fluc_train, val.param, gn1_val, val.A)
    fluctuate!(qcv, g, fluc_train, fluc_val)
end

function makeQCV{T<:AbstractFloat}(logitQnA1::Vector{T}, logitQnA0::Vector{T},
                                   W::Matrix{T}, A::Vector{T}, Y::Vector{T},
                                   param::Parameter{T}, searchstrategy::SearchStrategy,
                                   idx_train::AbstractVector{Int}, gbounds::Vector{T},
                                   penalize::Bool)
    n, p = size(W)
    idx_train = sort(idx_train)
    idx_val = setdiff(1:n, idx_train)
    chunk_train = make_chunk_subset(logitQnA1, logitQnA0, W, A, Y, param, idx_train, gbounds, penalize)
    chunk_val= make_chunk_subset(logitQnA1, logitQnA0, W, A, Y, param, idx_val, gbounds, false)
    used_covars = IntSet()
    unused_covars = setdiff(IntSet(1:p), used_covars)
    gseq = LR{T}[]
    QCV(chunk_train, chunk_val, gseq, deepcopy(searchstrategy), used_covars, unused_covars)
end

immutable FluctuationInfo
    steps::Int
    covar_order::Vector{Vector{Int}}
    new_flucseq::Vector{Bool}
    added_each_fluc::Vector{Vector{Int}}
    function FluctuationInfo(steps::Int, covar_order::Vector{Vector{Int}},
                             new_flucseq::Vector{Bool},
                             added_each_fluc::Vector{Vector{Int}})
        steps == length(covar_order) == length(new_flucseq) || throw(ArgumentError("steps should equal length of covar_order and new_flucseq"))
        length(added_each_fluc) == sum(new_flucseq) || throw(ArgumentError("length(added_each_fluc) not equal to number of new flucs"))
        new(steps, covar_order, new_flucseq, added_each_fluc)
    end
end

function FluctuationInfo{T}(steps::Int, gseq::Vector{LR{T}}, new_flucseq::Vector{Bool})
    new_flucseq[1] || throw(ArgumentError("new_flucseq[1] should be true"))
    first_covars = Vector{Int}[convert(Vector{Int}, gseq[1].idx)]::Vector{Vector{Int}}
    added_covar_order = [setdiff(a.idx, b.idx)::Vector{Int} for (a, b) in zip(gseq[2:end], gseq[1:end-1])]::Vector{Vector{Int}}
    covar_order = [first_covars; added_covar_order]::Vector{Vector{Int}}
    start_new_fluc = find(new_flucseq)
    end_fluc = [find(new_flucseq)[2:end] - 1; length(new_flucseq)]
    added_each_fluc = [vcat(covar_order[st:en]...)::Vector{Int}
                       for (st, en) in zip(start_new_fluc, end_fluc)]
    FluctuationInfo(steps, covar_order, new_flucseq, added_each_fluc)
end

function Base.show(io::IO, flucinfo::FluctuationInfo)
    println(io, "Fluctuation info:")
    println(io, "Covariates added in $(flucinfo.steps) steps to $(length(flucinfo.added_each_fluc)) fluctuations.")
    for (i, new_covars) in enumerate(flucinfo.added_each_fluc)
        println(io, "Fluc $i covars added: $new_covars")
    end
end

type CTMLE{T<:AbstractFloat} <: AbstractScalarEstimate
    psi::T
    ic::Vector{T}
    n::Int
    estimand::AbstractString
    seasrchstrategy::SearchStrategy
    flucinfo::FluctuationInfo
    function CTMLE(psi, ic, n, estimand, searchstrategy, flucinfo)
        est = ScalarEstimate(psi, ic)
        new(est.psi, est.ic, n, estimand, searchstrategy, flucinfo)
    end
end

CTMLE{T<:AbstractFloat}(psi::T, ic::Vector{T}, n::Int, estimand::AbstractString,
                        searchstrategy::SearchStrategy, flucinfo::FluctuationInfo) =
    CTMLE{T}(psi, ic, n, estimand, searchstrategy, flucinfo)

"""
Performs colaborative targeted minimum loss-based estimation

** Arguments **

* `logitQnA1` - Vector of length n of initial estimates logit(\bar{Q}_n(1, W_i))
* `logitQnA0` - Vector of length n of initial estimates logit(\bar{Q}_n(0, W_i))
* `W` - Matrix of covariates to be potentially used to estimate g. n rows.
        The first column should be all ones for an intercept.
* `A` - Vector of length n of treatments, 0.0 or 1.0
* `Y` - Vector of length n of outcomes, bounded between 0 and 1.

** Keyword Arguments **

* `param` - Target parameter. See ?ATE or ?Mean. Default is ATE()
* `gbounds` - values for bounding g. Estimates of g will be bounded by `extrema(gbounds)`.
   To only bound away from 0 (or 1), set to e.g. `[0.01, 1.0]` (or `[0.0, 0.99]`).
   To turn off bounding all together, set to `[0.0, 1.0]`. Defaults to `[0.01, 0.99]`.
* `penalize_risk` - should a penalized risk be used when choosing covariates to estimate g?
   Defaults to false.
* `searchstrategy` - Strategy for adding covariates to estimates of g. See ?SearchStrategy.
* `cvplan` - An iterator of vectors of Int indexes of observations in each training set for
   cross validaiton. StratifiedKfold and StratifiedRandomSub from the MLBase package are useful
   here. Defaults to 10 fold CV stratifying by treatment for outcomes with more than 2 levels,
   or by treatment and outcome for binary outcomes.
* `patience` - For how many steps should CV continue after finding a local optimum? Defaults to typemax(Int)

"""
function ctmle{T<:AbstractFloat}(logitQnA1::Vector{T}, logitQnA0::Vector{T},
                                 W::Matrix{T}, A::Vector{T}, Y::Vector{T};
                                 param::Parameter{T}=ATE(),
                                 gbounds::Vector{T}=[0.01, 0.99],
                                 penalize_risk::Bool=false,
                                 searchstrategy::SearchStrategy = ForwardStepwise(),
                                 cvplan = StratifiedKfold(length(unique(Y))<3? zip(A, Y) : A, 10),
                                 patience::Int=typemax(Int)
                                )
    n,p = size(W)
    n == length(logitQnA1) ==
        length(logitQnA0) ==
        length(A) ==
        length(Y) ||
        throw(ArgumentError("logitQnA1, logitQnA0, A, and Y should have length equal to the number of rows of W."))

    all(W[:, 1] .== 1) || throw(ArgumentError("The first column of W should be all ones."))

    #create vector of QCV objects
    cvqs = [makeQCV(logitQnA1, logitQnA0, W, A, Y, param, searchstrategy,
                    idx_train, gbounds, penalize_risk)::QCV{T}
            for idx_train in cvplan]

    #using cross-validation, find best number of steps to take
    best_steps = find_steps(cvqs, patience=patience)

    #build a chunk with the full data set
    fullchunk = Qchunk(Qmodel(logitQnA1, logitQnA0), W, A, Y, param, 1:n, gbounds, penalize_risk, convert(T, Inf))
    gseq = LR{T}[]
    new_flucseq = Bool[]
    used_covars = IntSet()
    unused_covars = IntSet(1:p)
    #and add covariates steps # of times, saving g fits
    @info("Fitting final estimate of Q")
    for step in 1:best_steps
#         @debug("used_covars: $used_covars")
#         @debug("unused_covars: $unused_covars")
        gfit, new_fluc = add_covars!(fullchunk, searchstrategy, used_covars, unused_covars)
        push!(gseq, gfit)
        push!(new_flucseq, new_fluc)
        @info(if length(gseq) == 1
                  "Step $step of $best_steps complete, covariate [1] added in new fluctuation."
              else
                  whichfluc = new_fluc? "in new" : "to current"
                  whichcovars = setdiff(gseq[end].idx, gseq[end-1].idx)
                  "Step $step of $best_steps complete, covariates $whichcovars added $whichfluc fluctuation."
              end)
    end

    #create flucinfo
    flucinfo = FluctuationInfo(best_steps, gseq, new_flucseq)

    q=fullchunk.q
    CTMLE(applyparam(param, q, A, Y)..., nobs(q), estimand(param), searchstrategy, flucinfo)

end

flucinfo(est::CTMLE) = est.flucinfo

"""
Determines the number of times covariates should be added using cross validation.

** Arguments **

* `qcvs` - A vector of one or more `QCV`s

** Keyword Arguments **

* `patience` - An integer. How many more steps should we try after finding a local minimum?
Defaults to a big number.

** Details **

Adds covariates to each `qcv` and computes the mean risk on validation chunks.
If minimum validation risk occured more than `patience` steps ago, we stop looking.
So for example if `patience` is `2`, we try for two more steps after finding a local minimum to see
if we can do better.
If `patience` is `0`, we stop as soon as risk starts increasing.
If `pateince` is negative, we stop after 1 step, which is probably not a good idea.

Returns best risk and number of steps.

"""
function find_steps{T<:AbstractFloat}(qcvs::Vector{QCV{T}}; patience::Int=typemax(Int))
    notdone = true
    steps = 0
    best_steps = 0
    best_risk = Inf

    @info("Choosing number of steps via cross validation")
    while notdone
        steps += 1
        map(add_covars!, qcvs)
        val_risk = mean(map(qcv -> qcv.chunk_val.risk, qcvs))
        if val_risk < best_risk
            best_risk = val_risk
            best_steps = steps
        end
        @info("Step $steps validation risk: $val_risk. Best risk so far: $best_risk at step $best_steps")
        if steps - best_steps > patience
            @info("No improvement in $patience steps, stopping at step $steps.")
            notdone = false
        end
        if all(map(qcv -> isempty(qcv.unused_covars), qcvs))
            @info("Crossvalidation completed after $steps steps.")
            notdone = false
        end
    end
    best_steps
end

"""
Adds next covariate(s) to a Qchunk based on a search strategy

** Arguments **

* `chunk` - a `Qchunk`
* `searchstrategy` - a `SearchStrategy`
* `used_covars` - indexes of covariates already included in g
* `unused_covars` - indexes of covariates not yet included in g

** Details **

If `used_covars` is emtpy (because this is the first time `add_covars!` has been called on this chunk,
a fluctuation based on an intercept only model is added to Q.
Otherwise, next covariates are added where the order is determined by `searchstrategy`

Returns most recent `gfit` object and a boolean which is `true` if the most recent covariates were added to a new fluctuation
or `false` if the previous fluctuation was updated with new covariates.

"""
function add_covars!(chunk::Qchunk, searchstrategy::SearchStrategy, used_covars, unused_covars)
    if isempty(used_covars)
#         @debug("Initial fluctuation")
        push!(used_covars, 1)
        delete!(unused_covars, 1)
        gfit = lreg(chunk.W, chunk.A, subset=1:1)
        fluc = computefluc(chunk.q, chunk.param, bound!(predict(gfit, chunk.W), chunk.gbounds), chunk.A, chunk.Y)
        fluctuate!(chunk.q, fluc)
        chunk.risk = risk(chunk.q,chunk.A, chunk.Y, chunk.param, chunk.penalize)
        return gfit, true
    else
#         @debug("Adding covariate(s)")
        abc = add_covars!(searchstrategy, chunk.q, chunk.param, chunk.W, chunk.A, chunk.Y,
                          used_covars, unused_covars, chunk.risk, chunk.gbounds, chunk.penalize)
        q_risk, gfit, new_fluc = abc
        chunk.risk=q_risk
        return gfit, new_fluc
    end
end

"""
Adds next covariate(s) to a QCV.

** Arguments **

* `qcv` - a `QCV`

** Details **
If `qcv` has no unused covariates, none are added.  New covariates are first added to the training chunk,
`qcv.train`. Then the validation chunk fluctuated using the newest `gfit`.
Returns updated `qcv`.

"""
function add_covars!(qcv::QCV)
    if isempty(qcv.unused_covars)
        @debug("No covariates to add")
        return qcv
    end
    train = qcv.chunk_train
    val = qcv.chunk_val
    #find next fluctuation for q_train
    # update q_val, which dpeends on g fit on q_train, and fluc

    #add next covar to training chunk and store gfit
    gfit, new_fluc = add_covars!(train, qcv.searchstrategy, qcv.used_covars, qcv.unused_covars)
    push!(qcv.gseq, gfit)

    #add newest fluctuation to val chunk and compute risk
    val = qcv.chunk_val
    gn1_val = bound!(predict(gfit, val.W), val.gbounds)
    if !new_fluc
        defluctuate!(val.q)
    end
    fluctuate!(val.q, valfluc(lastfluc(train.q), val.param, gn1_val, val.A))
    val.risk = risk(val.q, val.A, val.Y, val.param, false)

    qcv
end

function bound!{T<:AbstractFloat}(gn1::Vector{T}, gbounds::Vector{T})
    mn, mx = extrema(gbounds)
    mn < 0.5 || ArgumentError("minimum(gbounds) should be < 0.5")
    mx > 0.5 || ArgumentError("maximum(gbounds) should be > 0.5")

    mn > 0.0 || mx < 1.0 || return gn1

    for i in 1:length(gn1)
        @inbounds gn1[i] = max(mn, min(mx,gn1[i]))
    end
    gn1
end

function risk{T<:AbstractFloat}(q::Qmodel{T}, A::Vector{T}, Y::Vector{T}, param::Parameter{T}, penalize::Bool)
    total_loss = 0.0
    length(Y) == nobs(q) || error()
    lp = linpred(q, A)
    for i in 1:nobs(q)
      total_loss += LReg.loss(Y[i], lp[i])
    end
    q_risk = total_loss/nobs(q)

    if penalize
        ic = applyparam(param, q, A, Y)[2]
        pen = var(ic) / length(A)
    else
        pen = 0.0
    end
    q_risk + pen
end

end # module


text:
# ---
# title: 313. Super Ugly Number
# id: problem313
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Math, Heap
# link: <https://leetcode.com/problems/super-ugly-number/description/>
# hidden: true
# ---
#
# Write a program to find the `nth` super ugly number.
#
# Super ugly numbers are positive numbers whose all prime factors are in the
# given prime list `primes` of size `k`.
#
# **Example:**
#
#
#
#     Input: n = 12, primes = [2,7,13,19]
#     Output: 32
#     Explanation:[1,2,4,7,8,13,14,16,19,26,28,32] is the sequence of the first 12
#                  super ugly numbers given primes = [2,7,13,19] of size 4.
#
# **Note:**
#
#   * `1` is a super ugly number for any given `primes`.
#   * The given numbers in `primes` are in ascending order.
#   * 0 < `k` ≤ 100, 0 < `n` ≤ 106, 0 < `primes[i]` < 1000.
#   * The nth super ugly number is guaranteed to fit in a 32-bit signed integer.
#
#
## @lc code=start
using LeetCode

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

ext: jl


text:
using Expect
using Test

function relative(a::Float64, b::Float64)::Expected{Float64, ErrorException}
    if a == 0
        return ErrorException("cannot compute relative to 0")
    else
        return (b-a) / a
    end
end

x = relative(3.0, 2.0)
@test isassigned(x)
@test x[] == (2.0 - 3.0) / 3.0
@inferred x[]
y = relative(0.0, 2.0)
@test !isassigned(y)
@test_throws ErrorException y[]

function good_int_or_bad_int(a)::Expected{Int, Int}
    if a < 3
        return a
    else
        return unexpected(a)
    end
end

z = good_int_or_bad_int(1)
@test isassigned(z)
q = good_int_or_bad_int(4)
@test !isassigned(q)
@test_throws 4 q[]


text:
using Flux
using CuArrays

export NeuralNet

# This code is originally written by Avik Pal
struct ResidualBlock
    conv_layers
    norm_layers
    shortcut
end

Flux.@treelike ResidualBlock

# ResidualBlock Function allows us to define a Residual Block having any number of Convolution and Batch Normalization Layers
function ResidualBlock(filters, kernels::Array{Tuple{Int,Int}},
                       pads::Array{Tuple{Int,Int}},
                       strides::Array{Tuple{Int,Int}},
                       shortcut = identity)
    local conv_layers = []
    local norm_layers = []
    for i in 2:length(filters)
        @inbounds push!(conv_layers, Conv(kernels[i-1], filters[i-1]=>filters[i], pad = pads[i-1], stride = strides[i-1]) |> gpu)
        push!(norm_layers, BatchNorm(filters[i])|>gpu)
    end
    ResidualBlock(Tuple(conv_layers),Tuple(norm_layers),shortcut)
end

# Function converts the Array of scalar kernel, pad and stride values to tuples
function ResidualBlock(filters, kernels::Array{Int}, pads::Array{Int}, strides::Array{Int}, shortcut = identity)
    ResidualBlock(filters, [(i,i) for i in kernels], [(i,i) for i in pads], [(i,i) for i in strides], shortcut)
end

function (block::ResidualBlock)(input)
    local value = copy.(input)
    for i in 1:length(block.conv_layers)-1
        @inbounds value = relu.((block.norm_layers[i])((block.conv_layers[i])(value)))
    end
    @inbounds relu.(((block.norm_layers[end])((block.conv_layers[end])(value))) + block.shortcut(input))
end

struct NeuralNet
    base_net::Chain
    value::Chain
    policy::Chain
end

Flux.@treelike NeuralNet

function NeuralNet(N::Int = 9, tower_height::Int = 20, features::Int=256)

    res_block() = ResidualBlock([features,features,features], [3,3], [1,1], [1,1])
    # 9 residual blocks
    tower = [res_block() for i = 1:tower_height]

    # feature extraction
    base_net = Chain(Conv((3,3), 9=>features, pad=(1,1)),
                     BatchNorm(features, leakyrelu),
                     tower...) |> gpu

    # value head
    value = Chain(Conv((1,1), features=>1),
                  BatchNorm(1, relu), x->reshape(x, N^2, :),
                  Dense(N^2, features, leakyrelu),
                  Dense(features, 1, tanh)) |> gpu

    policy = Chain(Conv((1,1), features=>2),
                   BatchNorm(2, relu), x->reshape(x, 2N^2, size(x, 4)),
                   Dense(2N^2, 9), softmax) |> gpu

    NeuralNet(base_net, value, policy)
end

function Base.show(io::IO, nn::NeuralNet)
    compact = get(io, :compact, false)
    if !compact
        print(io, "NeuralNet()")
    end
end

function (nn::NeuralNet)(input::AbstractArray{Float32,4})
    features = nn.base_net(input)
    p = nn.policy(features)
    v = nn.value(features)
    return p, v
end


text:
export Trajectory,
    PrioritizedTrajectory,
    CircularArrayTrajectory,
    CircularVectorTrajectory,
    CircularArraySARTTrajectory,
    CircularArraySLARTTrajectory,
    CircularVectorSARTTrajectory,
    CircularVectorSARTSATrajectory,
    CircularArrayPSARTTrajectory,
    ElasticArrayTrajectory,
    ElasticSARTTrajectory,
    VectorTrajectory,
    VectorSATrajectory,
    VectorSARTTrajectory

using MacroTools: @forward
using ElasticArrays
using CircularArrayBuffers: CircularArrayBuffer, CircularVectorBuffer

#####
# Trajectory
#####

"""
    Trajectory(;[trace_name=trace_container]...)

A simple wrapper of `NamedTuple`.
Define our own type here to avoid type piracy with `NamedTuple`
"""
struct Trajectory{T} <: AbstractTrajectory
    traces::T
end

Trajectory(; kwargs...) = Trajectory(kwargs.data)

@forward Trajectory.traces Base.getindex, Base.keys

Base.merge(a::Trajectory, b::Trajectory) = Trajectory(merge(a.traces, b.traces))
Base.merge(a::Trajectory, b::NamedTuple) = Trajectory(merge(a.traces, b))
Base.merge(a::NamedTuple, b::Trajectory) = Trajectory(merge(a, b.traces))

#####

function CircularArrayTrajectory(; capacity, kwargs...)
    Trajectory(map(kwargs.data) do x
        CircularArrayBuffer{eltype(first(x))}(last(x)..., capacity)
    end)
end

function CircularVectorTrajectory(; capacity, kwargs...)
    Trajectory(map(kwargs.data) do x
        CircularVectorBuffer{x}(capacity)
    end)
end

#####

const CircularArraySARTTrajectory = Trajectory{
    <:NamedTuple{
        SART,
        <:Tuple{
            <:CircularArrayBuffer,
            <:CircularArrayBuffer,
            <:CircularArrayBuffer,
            <:CircularArrayBuffer,
        },
    },
}

CircularArraySARTTrajectory(;
    capacity::Int,
    state = Int => (),
    action = Int => (),
    reward = Float32 => (),
    terminal = Bool => (),
) = merge(
    CircularArrayTrajectory(; capacity = capacity + 1, state = state, action = action),
    CircularArrayTrajectory(; capacity = capacity, reward = reward, terminal = terminal),
)

const CircularArraySLARTTrajectory = Trajectory{
    <:NamedTuple{
        SLART,
        <:Tuple{
            <:CircularArrayBuffer,
            <:CircularArrayBuffer,
            <:CircularArrayBuffer,
            <:CircularArrayBuffer,
            <:CircularArrayBuffer,
        },
    },
}

CircularArraySLARTTrajectory(;
    capacity::Int,
    state = Int => (),
    legal_actions_mask,
    action = Int => (),
    reward = Float32 => (),
    terminal = Bool => (),
) = merge(
    CircularArrayTrajectory(;
        capacity = capacity + 1,
        state = state,
        legal_actions_mask = legal_actions_mask,
        action = action,
    ),
    CircularArrayTrajectory(; capacity = capacity, reward = reward, terminal = terminal),
)

#####

const CircularVectorSARTTrajectory = Trajectory{
    <:NamedTuple{
        SART,
        <:Tuple{
            <:CircularVectorBuffer,
            <:CircularVectorBuffer,
            <:CircularVectorBuffer,
            <:CircularVectorBuffer,
        },
    },
}

CircularVectorSARTTrajectory(;
    capacity::Int,
    state = Int,
    action = Int,
    reward = Float32,
    terminal = Bool,
) = merge(
    CircularVectorTrajectory(; capacity = capacity + 1, state = state, action = action),
    CircularVectorTrajectory(; capacity = capacity, reward = reward, terminal = terminal),
)

#####

const CircularVectorSARTSATrajectory = Trajectory{
    <:NamedTuple{
        SARTSA,
        <:Tuple{
            <:CircularVectorBuffer,
            <:CircularVectorBuffer,
            <:CircularVectorBuffer,
            <:CircularVectorBuffer,
            <:CircularVectorBuffer,
            <:CircularVectorBuffer,
        },
    },
}

CircularVectorSARTSATrajectory(;
    capacity::Int,
    state = Int,
    action = Int,
    reward = Float32,
    terminal = Bool,
    next_state = state,
    next_action = action,
) = CircularVectorTrajectory(;
    capacity = capacity,
    state = state,
    action = action,
    reward = reward,
    terminal = terminal,
    next_state = next_state,
    next_action = next_action,
)

#####

function ElasticArrayTrajectory(; kwargs...)
    Trajectory(map(kwargs.data) do x
        ElasticArray{eltype(first(x))}(undef, last(x)..., 0)
    end)
end

const ElasticSARTTrajectory = Trajectory{
    <:NamedTuple{SART,<:Tuple{<:ElasticArray,<:ElasticArray,<:ElasticArray,<:ElasticArray}},
}

function ElasticSARTTrajectory(;
    state = Int => (),
    action = Int => (),
    reward = Float32 => (),
    terminal = Bool => (),
)
    ElasticArrayTrajectory(;
        state = state,
        action = action,
        reward = reward,
        terminal = terminal,
    )
end

#####
# VectorTrajectory
#####

function VectorTrajectory(; kwargs...)
    Trajectory(map(kwargs.data) do x
        Vector{x}()
    end)
end

const VectorSARTTrajectory =
    Trajectory{<:NamedTuple{SART,<:Tuple{<:Vector,<:Vector,<:Vector,<:Vector}}}

function VectorSARTTrajectory(;
    state = Int,
    action = Int,
    reward = Float32,
    terminal = Bool,
)
    VectorTrajectory(; state = state, action = action, reward = reward, terminal = terminal)
end

const VectorSATrajectory =
    Trajectory{<:NamedTuple{(:state, :action),<:Tuple{<:Vector,<:Vector}}}

function VectorSATrajectory(; state = Int, action = Int)
    VectorTrajectory(; state = state, action = action)
end
#####

Base.@kwdef struct PrioritizedTrajectory{T,P} <: AbstractTrajectory
    traj::T
    priority::P
end

Base.keys(t::PrioritizedTrajectory) = (:priority, keys(t.traj)...)

Base.length(t::PrioritizedTrajectory) = length(t.priority)

Base.getindex(t::PrioritizedTrajectory, s::Symbol) =
    if s == :priority
        t.priority
    else
        getindex(t.traj, s)
    end

const CircularArrayPSARTTrajectory =
    PrioritizedTrajectory{<:SumTree,<:CircularArraySARTTrajectory}

CircularArrayPSARTTrajectory(; capacity, kwargs...) = PrioritizedTrajectory(
    CircularArraySARTTrajectory(; capacity = capacity, kwargs...),
    SumTree(capacity),
)

#####
# Common
#####

function Base.length(
    t::Union{
        CircularArraySARTTrajectory,
        CircularArraySLARTTrajectory,
        CircularVectorSARTSATrajectory,
        ElasticSARTTrajectory,
    },
)
    x = t[:terminal]
    size(x, ndims(x))
end

Base.length(t::VectorSARTTrajectory) = length(t[:terminal])
Base.length(t::VectorSATrajectory) = length(t[:action])


text:
using SequentialLOB
using Test

include("../src/main.jl")

@testset "All Tests" begin
    slob_model = SLOB()
    @testset "Reproducibility Non-Debug" begin
        slob_model_1 = SLOB(1, 200, 238.745, 500, 100.0, 2.0,
            1.0, 0.1, 20.0, SourceTerm(1.0, 0.5))
        slob_model_2 = SLOB(1, 200, 238.745, 500, 100.0, 2.0,
            1.0, 0.1, 20.0, SourceTerm(1.0, 0.5))
        @test all(slob_model_1(45) .== slob_model_2(45))

        @test all(slob_model_1(51) .== slob_model_2(51))

        @test all(slob_model_1(4565756745) .== slob_model_2(4565756745))
    end;

    @testset "Reproducibility Debug" begin
        slob_model_1 = SLOB(1, 200, 238.745, 500, 100.0, 2.0,
            1.0, 0.1, 20.0, SourceTerm(1.0, 0.5))
        slob_model_2 = SLOB(1, 200, 238.745, 500, 100.0, 2.0,
            1.0, 0.1, 20.0, SourceTerm(1.0, 0.5))
        @test all(slob_model_1(true, 45) .== slob_model_2(true, 45))

        @test all(slob_model_1(true, 51) .== slob_model_2(true, 51))

        @test all(slob_model_1(true, 4565756745) .== slob_model_2(true, 4565756745))
    end;

    @testset "Default Values" begin
        lob_densities, price_paths, mid_price_bars, P⁺s, P⁻s = slob_model(true)
        @test size(mid_price_bars) == (101,1)
    end;

    @testset "Command Line Parse Default Arguments" begin
        @test parse_commandline() == Dict("μ" => 0.5,
            "num_paths" => 1,
            "T" => 100,
            "p₀" => 100.0,
            "λ" => 1.0,
            "SEED" => 1,
            "nu" => 0.1,
            "σ" => 1.0,
            "D" => 4.0,
            "M" => 100,
            "L" => 50.0,
            "α" => 20.0)
    end

    @testset "Main Default Arguments" begin
        @test all(main("return") .== slob_model(1))
    end
end


text:
using vOptSpecific

n  =  3
C1 = [ 3  9  7 ; 16 10  6 ; 2  7 11]
C2 = [16 15  6 ;  5  7 13 ; 1  2 13]

myInstance = set2LAP( n , C1 , C2 )
algorithm = LAP_Przybylski2008( )
z1, z2, S = vSolve( myInstance , algorithm )

id = set2LAP(n, C1, C2)
z1,z2,solutions = vSolve(id)

for i = 1:length(z1)
    println("(" , z1[i], " | ", z2[i], ") : ", solutions[i,:])
end


text:
# (2017.07.21) copied from IJCNN2017 code base (file: test1_bd.jl)
#
# Test Bayesian Detector

include("./detectors/2bayes_detector_replica.jl")
# include("./lib/SignalGenerator.jl")
include("./lib/roc.jl")
# using SignalGenerator
using BayesChpOrig
using JSON
# using Metadata
using Roc

function run_sig_gen()
    ret = genSig(5, 100.0, 1.0, 2.0, 0.5, 1)
    sig = ret["sig"]
    chps= ret["changes"]
    writedlm(outdir*"testsig.txt", sig)
    writedlm(outdir*"testsigchps.txt", chps)
    print("run1")
end

"""
Simply run Bayesian detector on a single one dimensional signal.
Save detections and changes into files.
Plot using R.
"""
function run_detector()
    jcont = JSON.parsefile("./config.json")

    sig_dat = JSON.parsefile(jcont["input"]["art_sig1"])
    sig = sig_dat["sig"]
    #ret = genSig(10, 100.0, 10.0, 1.3, 0.75, 1)
    chps= sig_dat["changes"]

    matrix_max_vals, matrix = bayes_detector(sig, lambda = 120)
    detections = extract_changes(matrix_max_vals)

    max_delay = 30
    tpr = tpRate(detections, chps, max_delay)
    fpr = faRateScaled(detections, chps, max_delay, length(sig), tn_price_scale=0.1)
    println("... test_bd > TPrate = ", tpr, " / FArate = ", fpr)

    writedlm(jcont["output"]["bchp_test1"]["detections"], detections)
    writedlm(jcont["output"]["bchp_test1"]["signal"], sig)
    writedlm(jcont["output"]["bchp_test1"]["changes"], chps)
    writedlm(jcont["output"]["bchp_test1"]["matrix"], matrix, ",")
    writedlm(jcont["output"]["bchp_test1"]["matrix_max_vals"], matrix_max_vals, ",")
end
run_detector()

max_stars_repo_name: av-maslov/Recurrency

max_issues_repo_name: av-maslov/Recurrency

max_forks_repo_name: av-maslov/Recurrency


text:
module ArrayIteration

import Base: getindex, setindex!, iterate, length, eachindex, show, parent, isless
using Base: ReshapedArray, ReshapedIndex, IndexStyle, IndexLinear, IndexCartesian
using Base: Slice, OneTo
using Base.PermutedDimsArrays: PermutedDimsArray
using Base.Order

using SparseArrays

export Follower, inds, index, value, stored, each, sync

include("types.jl")
include("core.jl")
include("reshaped.jl")
include("sparse.jl")
include("sync_stored.jl")

end # module


text:
using GeometricAlgebra

v1, v2 = basis((1, 1))
v12 = v1*v2

x, y = basis((x=1, y=1))
xy = x*y

@testset "repr parseability" begin
	isparseable(a) = eval(Meta.parse(repr(a))) == a

	@test isparseable(4v1)
	@test isparseable(0v1)
	@test isparseable((1 + im)*v1)
	@test isparseable(v1 - 2v2)
	@test isparseable(4 + v12)

	@test isparseable(10x)
	@test isparseable(0*x)
	@test isparseable(0.0x)
	@test isparseable(1 + 2x + 3xy)
end

# just test against errors
for to_show ∈ (x, x + y, 1 + y, Any[100x, 2x*y, 1 + x, 0 + 0y], typeof(x), typeof(x + y))
	println()
	@test isnothing(show(stdout, MIME("text/plain"), to_show))
end

@testset "pretty-printed types" begin
	for sig ∈ [(1, 1, 1), (x=1, y=1, z=1)]
		pretty_sig = GeometricAlgebra.show_signature(sig)
		compare(a, b) = startswith(replace(a, r"\s+"=>""), replace(b, r"\s+"=>""))
		@test compare(sprint(show, Blade{sig,2,0b011,Float64}), "Blade{$pretty_sig, 2, 0b011, Float64}")
		@test compare(sprint(show, Blade{sig,2,:bits}), "Blade{$pretty_sig, 2, :bits")
		@test compare(sprint(show, Blade{sig,2,0b101,Float64} where sig), "Blade{sig,2, 0b101, Float64} where sig")
		@test compare(sprint(show, Multivector{sig,k,S} where {k,S}), "Multivector{$pretty_sig")
	end
	@test sprint(show, MixedMultivector) == "MixedMultivector"
end

@testset "blade printing" begin
	@test repr(Blade{(1,1,1)}(42, 0b11)) == "42v12"
	@test repr(Blade{(x=1,y=1,z=1)}(42, 0b11)) == "42xy"
	@test repr(Blade{OffsetSignature((-1,1,1,1),0:3)}(42, 0b11)) == "42v01"
	@test repr(Blade{OffsetSignature((t=-1,x=1,y=1,z=1),0:3)}(42, 0b11)) == "42tx"
end


text:
import Base: LineEdit, REPL

# original code from
# https://github.com/Keno/Cxx.jl/blob/master/src/CxxREPL/replpane.jl

function AddSwifterMode(key, repl)
    panel = LineEdit.Prompt("Swifter> ";
        prompt_prefix = "\e[0;3;5m",
        prompt_suffix = Base.text_colors[:white],
        on_enter = s->true)

    mirepl = isdefined(Base.active_repl,:mi) ? Base.active_repl.mi : Base.active_repl
    function linecall(line)
        try
            expr = parse(line)
            ex = isa(expr, Symbol) ? Expr(:block, expr) : expr
            Expr(:macrocall, Symbol("@query"), ex)
        catch e_
        end
    end
    panel.on_done = REPL.respond(linecall, repl, panel)

    main_mode = repl.interface.modes[1]
    push!(mirepl.interface.modes,panel)

    hp = main_mode.hist
    hp.mode_mapping[:Swifter] = panel
    panel.hist = hp

    const keymap = Dict{Any,Any}(
        key => function (s,args...)
            if isempty(s) || position(LineEdit.buffer(s)) == 0
                buf = copy(LineEdit.buffer(s))
                LineEdit.transition(s, panel) do
                    LineEdit.state(s, panel).input_buffer = buf
                end
            else
                LineEdit.edit_insert(s, key)
            end
        end
    )
    repl.interface = REPL.setup_interface(repl; extra_repl_keymap = keymap)

    search_prompt, skeymap = LineEdit.setup_search_keymap(hp)
    main_keymap = REPL.mode_keymap(main_mode)

    b = Dict{Any,Any}[skeymap, main_keymap, LineEdit.history_keymap, LineEdit.default_keymap, LineEdit.escape_defaults]
    panel.keymap_dict = LineEdit.keymap(b)

    main_mode.keymap_dict = LineEdit.keymap_merge(main_mode.keymap_dict, keymap)
    nothing
end

function RunSwifterREPL(; key=">")
    AddSwifterMode(key, Base.active_repl)
end


text:
using Libdl

const _DEPS_FILE = joinpath(dirname(@__FILE__), "deps.jl")
if isfile(_DEPS_FILE)
    rm(_DEPS_FILE)
end

const _CPX_VERS = [ # From oldest to most recent.
    "1210", "12100",
    "201", "2010", "20100",
]
const _BASE_ENV = "CPLEX_STUDIO_BINARIES"

function write_depsfile(path)
    open(_DEPS_FILE, "w") do f
        println(f, "const libcplex = \"$(escape_string(path))\"")
    end
end

function library_name(v)
    cpx_prefix = Sys.iswindows() ? "" : "lib"
    return "$(cpx_prefix)cplex$(v).$(Libdl.dlext)"
end

function default_installation_path(cplex_studio_path::AbstractString)
    if Sys.iswindows()
        return escape_string("C:\\Program Files\\IBM\\ILOG\\$cplex_studio_path\\cplex\\bin\\x64_win64\\")
    elseif Sys.isapple()
        return "/Applications/$cplex_studio_path/cplex/bin/x86-64_osx/"
    else
        return "/opt/$cplex_studio_path/cplex/bin/x86-64_linux/"
    end
end

function get_error_message_if_not_found()
    return """
    Unable to install CPLEX.jl.

    The versions of CPLEX supported by CPLEX.jl are:

    * 12.10
    * 20.1

    You must download and install one of these versions separately.

    You should set the `CPLEX_STUDIO_BINARIES` environment variable to point to
    the install location then try again. For example (updating the path to the
    correct location):

    ```
    ENV["CPLEX_STUDIO_BINARIES"] = "$(default_installation_path("CPLEX_Studio201"))"
    import Pkg
    Pkg.add("CPLEX")
    Pkg.build("CPLEX")
    ```

    See the CPLEX.jl README at https://github.com/jump-dev/CPLEX.jl for further
    instructions.
    """
end

function check_cplex_in_libnames(libnames)
    for l in libnames
        d = Libdl.dlopen_e(l)
        if d == C_NULL
            continue
        end
        return l
    end
    return nothing
end

function check_cplex_in_environment_variables()
    # Find CPLEX in the CPLEX environment variables.
    libnames = String[]

    for v in reverse(_CPX_VERS)
        name = library_name(v)

        # Library name is not always using the same suffix as in _CPX_VERS.
        # E.g., on Windows, mix between 201 and 2010 for 20.1:
        # C:\Program Files\IBM\ILOG\CPLEX_Studio201\opl\bin\x64_win64\cplex2010.dll
        for env in [_BASE_ENV, _BASE_ENV * v, [_BASE_ENV * v2 for v2 in reverse(_CPX_VERS) if v2 != v]...]
            if !haskey(ENV, env)
                continue
            end

            for d in split(ENV[env], ';')
                if isdir(d) && isfile(joinpath(d, name))
                    push!(libnames, joinpath(d, name))
                end
            end
        end
    end

    return libnames
end

function check_cplex_in_default_paths()
    # Find CPLEX in the default installation locations, based on the platform.
    libnames = String[]
    for v in reverse(_CPX_VERS)
        name = library_name(v)
        for product in ["CPLEX_Studio$v", "CPLEX_Enterprise_Server$v/CPLEX_Studio"]
            path = default_installation_path(product)
            if isdir(path)
                guessed_file = joinpath(path, name)
                if isfile(guessed_file)
                    push!(libnames, guessed_file)
                end
            end
        end
    end
    return libnames
end

function try_local_installation()
    # Iterate through a series of places where CPLEX could be found: either
    # from an environment variable, in the path (directly the callable library
    # or the CPLEX executable), or in a default install location, in that
    # order. Indeed, some software packages propose a version of CPLEX in the
    # PATH that is not useable from Julia.
    for libnames in [check_cplex_in_environment_variables(), check_cplex_in_default_paths()]
        found_cplex_lib = check_cplex_in_libnames(libnames)
        if found_cplex_lib !== nothing
            write_depsfile(Libdl.dlpath(found_cplex_lib))
            @info("Using CPLEX found in location `$(found_cplex_lib)`")
            return
        end
    end

    error(get_error_message_if_not_found())
end

function try_ci_installation()
    CPLEX_VERSION = ENV["CPLEX_VERSION"]
    url = ENV["SECRET_CPLEX_URL_" * CPLEX_VERSION]
    local_filename = joinpath(@__DIR__, "libcplex.so")
    download(url, local_filename)
    write_depsfile(local_filename)
end

if get(ENV, "JULIA_REGISTRYCI_AUTOMERGE", "false") == "true"
    # We need to be able to install and load this package without error for
    # Julia's registry AutoMerge to work. Just write a fake libcplex path.
    write_depsfile("julia_registryci_automerge")
elseif get(ENV, "SECRET_CPLEX_URL_12100", "") != ""
    try_ci_installation()
else
    try_local_installation()
end


text:
using ReconTPCF, Plots
using BenchmarkTools

#dims, C2, S2, philen = get_C2_S2("test/uniform_circles3.png")
#guess, S2n, C2n, S2_BN1, C2_BN1, SN1 = histrecon(dims, C2, S2, philen);
#heatmap(guess, color=:grays, aspect_ratio=1)
#save("test38.png", colorview(Gray, guess))
#imshow(guess)
function test()
    t = vec(CartesianIndices(zeros((100,100))))
    return ReconTPCF.blas_stat_st2(t, 1.0, 100.0)
end

@btime test()


text:
# Zigui writes the code to remove missing value and duplicate rows. And try to deal with the white space. But the code is too redundant and needs improvement.
function remove_missing_duplicate(file; delim = ',', header = 0,miss = missing)
    """
    Remove missing value and duplicate
    """
    dat = CSV.read(file, delim = delim, header = header)
    dat = unique!(dat) #remove duplicate
    nCol = ncol(dat)
    dat = dropmissing(dat)
    #Drop the missing value defined in julia

    if ismissing(miss)
        #Check the if the user defined missing value is the missing
        return dat
        #The user defined missing value is the missing, already removed

    else
    #The user defined missing value is not the missing, we need to remove the
    #rows containing user defined missing value
        for i in 1:nCol
            dat = dat[dat[:,i].!=miss,:] #Remove rows containing user defined missing value
        end
    end
    return dat

end

#Testing
dat = remove_missing_duplicate("pedigree_clean.txt",miss = missing)

#function can deal with the white space,it works but too redundant. Need imporvement.
function remove_missing_duplicate_2(file; delim = ',', header = 0,miss = missing)
    """
    Remove missing value and duplicate and can deal with the case of white space.
    But too redundant, needs improvement
    """
    dat = CSV.read(file, delim = delim, header = header) # no header
    dat = unique!(dat) #remove duplicate
    nCol = ncol(dat)
    nRow = nrow(dat)
    dat = dropmissing(dat) #Drop the julia missing value (missing)

    if ismissing(miss)
        #Check the if the user defined missing value is the missing
        return dat
        #The user defined missing value is the missing, already removed

    else
    #The user defined missing value is not the missing, we need to remove the
    #rows containing user defined missing value
        for i in 1:nCol
            row_record = []
            for j in 1:nrow(dat)
                dat[j,i] = strip(dat[j,i])
                if isempty(dat[j,i])
                    push!(row_record,j)
                end
            end
            dat = deleterows!(dat, row_record)
            dat = dat[dat[:,i].!=miss,:] #Remove rows containing user defined missing value
        end
    end
    return dat

end

#Testing
dat = remove_missing_duplicate_2("pedigree_clean.txt",miss = "0")

#Jiayi deal with the question that same individuals having different parents and remove missing value. And add three warnings which can be seen in the function documentation.
function ErrorDetection(file; delim = ',', header = 0, missingType = missing)
    """
    Input: pedigree.txt file
    Output: DataFrame
    Remove missing values given input missing Type, and give corresponding warnings
    Remove leading and trailing whitespace between values
    Remove rows that same individuals having different parents, and give corresponding warnings
    """
    dat = CSV.read(file, delim = delim, header = header)
    nCol = ncol(dat)

    # remove missing value
    if ismissing(missingType) # if missing values are just left empty in the pedigree
        MissingIndex = findall(ismissing, convert(Array,dat)) # an array of CartesianIndex
        indexForMissing = [value[1] for (i, value) in enumerate(MissingIndex)]
        dropmissing!(dat)
    else
        indexForMissing =  Array{Int,1}(undef, 0)
        for cols in 1:nCol
            Mindex = findall(x -> x == missingType, dat[cols]) # find all the rows that contain missing values
            append!(indexForMissing, Mindex)
        end
        unique!(indexForMissing) # get the unique combination of rows to delete
        deleterows!(dat, indexForMissing)
    end

   if !isempty(indexForMissing)
        sort!(indexForMissing)
        println("Warning: Missing values are found in row $indexForMissing and corresponding rows were deleted!")
   end

    # remove leading and trailing whitespace between values

    nCol = ncol(dat)
    for cols in 1:nCol
        if any(typeof.(dat[cols]) .== String)
        dat[cols] = strip.(dat[cols])
        end
    end

    # remove duplicate rows
    unique!(dat)

    index = collect(1:1:nrow(dat))

    # remove rows that same individuals having different parents
    SameInd = Array{Int,1}(undef, 0)
    Duplicate_index = dat[nonunique(dat, 1), 1]
    unique!(Duplicate_index)
    for ID in Duplicate_index
        IDindex = findall(x -> x == ID, dat[1])
        append!(SameInd, IDindex)
    end

    Duplicate_index=String.(Duplicate_index)
    println("Warning: Same individuals with different parents are found for individuals $Duplicate_index and corresponding rows were deleted!")
    deleterows!(dat, SameInd)
    return dat
end

#Tianjing re-write the code to deal with the problem that individual appears in both sire and dam.
function remove_dup_son(dat) #HAO: GOOD, SEE LINE 37
    """
    Goal:   remove duplicate son (this step is after removing all duplicate)
    Input:  a4 a1 a2
            a4 a9 a10
            a4 a6 a2
    Output: a4 a1 a2
    """
    dup_ele = nonunique(dat, 1)
    dup_index=findall(x -> x==true, dup_ele)
    deleterows!(dat, dup_index)

    return dat
end

#Testing
dat=remove_dup_son(dat)

function same_sire_dam(dat, root_signal = "0")
    """
    Goal:   remove element in both dam and sire
    Input:  a8  a4 a6
            a9  a4 a6
            a13 a6 a2
    Output: empty
    """
    dup = intersect(dat[2], dat[3])
    filter!(e -> e ≠ root_signal, dup)  #remove the dam and sire of root, eg. "0"
    dup_both = vcat(findall(x->x in dup, dat[2]),findall(x->x in dup, dat[3])) #find index of dup

    unique!(dup_both)
    sort!(dup_both)

    deleterows!(dat, dup_both)

    return dat
end

#Testing
dat = same_sire_dam(dat)

function quality_control(file; delim = ',', header = 0,miss = missing)
    #Jiayi
    dat = ErrorDetection(file; delim = delim, header = header, missingType = miss)

    #zigui
    #dat = remove_missing_duplicate_2(file; delim = delim, header = header,miss = miss) #Remove missing and duplicate

    #Tianjing
    dat = same_sire_dam(dat) #individual appear in sire and dam
    dat = remove_dup_son(dat) #individual has different father and mother
    return dat
end

#Testing
quality_control("pedigree_clean.txt",miss = missing)


text:
const PyJuliaSetValue_Type__ref = Ref(PyNULL)
PyJuliaSetValue_Type() = begin
    ptr = PyJuliaSetValue_Type__ref[]
    if isnull(ptr)
        c = []
        base = PyJuliaAnyValue_Type()
        isnull(base) && return PyNULL
        t = fill(
            PyType_Create(
                c,
                name = "juliaaa.SetValue",
                base = base,
                methods = [
                    (name = "add", flags = Py_METH_O, meth = pyjlset_add),
                    (name = "clear", flags = Py_METH_NOARGS, meth = pyjlset_clear),
                    (name = "copy", flags = Py_METH_NOARGS, meth = pyjlset_copy),
                    (name = "difference", flags = Py_METH_O, meth = pyjlset_difference),
                    (
                        name = "difference_update",
                        flags = Py_METH_O,
                        meth = pyjlset_difference_update,
                    ),
                    (name = "discard", flags = Py_METH_O, meth = pyjlset_discard),
                    (name = "intersection", flags = Py_METH_O, meth = pyjlset_intersection),
                    (
                        name = "intersection_update",
                        flags = Py_METH_O,
                        meth = pyjlset_intersection_update,
                    ),
                    (name = "isdisjoint", flags = Py_METH_O, meth = pyjlset_isdisjoint),
                    # (name="issubset", flags=Py_METH_O, meth=pyjlset_issubset),
                    # (name="issuperset", flags=Py_METH_O, meth=pyjlset_issuperset),
                    (name = "pop", flags = Py_METH_NOARGS, meth = pyjlset_pop),
                    (name = "remove", flags = Py_METH_O, meth = pyjlset_remove),
                    (
                        name = "symmetric_difference",
                        flags = Py_METH_O,
                        meth = pyjlset_symmetric_difference,
                    ),
                    (
                        name = "symmetric_difference_update",
                        flags = Py_METH_O,
                        meth = pyjlset_symmetric_difference_update,
                    ),
                    (name = "union", flags = Py_METH_O, meth = pyjlset_union),
                    (name = "update", flags = Py_METH_O, meth = pyjlset_update),
                ],
                # as_number = (
                #     or = pyjlset_or,
                #     and = pyjlset_and,
                #     xor = pyjlset_xor,
                #     subtract = pyjlset_sub,
                #     inplace_or = pyjlset_ior,
                #     inplace_and = pyjlset_iand,
                #     inplace_xor = pyjlset_ixor,
                #     inplace_subtract = pyjlset_isub,
                # ),
            ),
        )
        ptr = PyPtr(pointer(t))
        err = PyType_Ready(ptr)
        ism1(err) && return PyNULL
        abc = PyMutableSetABC_Type()
        isnull(abc) && return PyNULL
        ism1(PyABC_Register(ptr, abc)) && return PyNULL
        PYJLGCCACHE[ptr] = push!(c, t)
        PyJuliaSetValue_Type__ref[] = ptr
    end
    ptr
end

PyJuliaSetValue_New(x::AbstractSet) = PyJuliaValue_New(PyJuliaSetValue_Type(), x)
PyJuliaValue_From(x::AbstractSet) = PyJuliaSetValue_New(x)

pyjlset_add(xo::PyPtr, vo::PyPtr) = begin
    x = PyJuliaValue_GetValue(xo)::AbstractSet
    ism1(PyObject_Convert(vo, eltype(x))) && return PyNULL
    v = takeresult(eltype(x))
    try
        push!(x, v)
        PyNone_New()
    catch err
        PyErr_SetJuliaError(err)
        PyNULL
    end
end

pyjlset_discard(xo::PyPtr, vo::PyPtr) = begin
    x = PyJuliaValue_GetValue(xo)::AbstractSet
    r = PyObject_TryConvert(vo, eltype(x))
    r == -1 && return PyNULL
    r == 0 && return PyNone_New()
    v = takeresult(eltype(x))
    try
        delete!(x, v)
        PyNone_New()
    catch err
        PyErr_SetJuliaError(err)
        PyNULL
    end
end

pyjlset_clear(xo::PyPtr, _::PyPtr) = begin
    x = PyJuliaValue_GetValue(xo)::AbstractSet
    try
        empty!(x)
        PyNone_New()
    catch err
        PyErr_SetJuliaError(err)
        PyNULL
    end
end

pyjlset_copy(xo::PyPtr, _::PyPtr) = begin
    x = PyJuliaValue_GetValue(xo)::AbstractSet
    try
        PyObject_From(copy(x))
    catch err
        PyErr_SetJuliaError(err)
        PyNULL
    end
end

pyjlset_pop(xo::PyPtr, _::PyPtr) = begin
    x = PyJuliaValue_GetValue(xo)::AbstractSet
    try
        if isempty(x)
            PyErr_SetString(PyExc_KeyError(), "pop from an empty set")
            PyNULL
        else
            PyObject_From(pop!(x))
        end
    catch err
        PyErr_SetJuliaError(err)
        PyNULL
    end
end

pyjlset_remove(xo::PyPtr, vo::PyPtr) = begin
    x = PyJuliaValue_GetValue(xo)::AbstractSet
    r = PyObject_TryConvert(vo, eltype(x))
    if r == -1
        return PyNULL
    elseif r == 0
        PyErr_SetObject(PyExc_KeyError(), vo)
        return PyNULL
    end
    v = takeresult(eltype(x))
    try
        if v in x
            delete!(x, v)
            PyNone_New()
        else
            PyErr_SetObject(PyExc_KeyError(), vo)
            PyNULL
        end
    catch err
        PyErr_SetJuliaError(err)
        PyNULL
    end
end

pyjlset_ibinop_named(xo, yo, op, skip) = begin
    x = PyJuliaValue_GetValue(xo)::AbstractSet
    y = PyIterable_Collect(yo, eltype(x), skip)
    isempty(y) && PyErr_IsSet() && return PyNULL
    try
        op(x, y)
        PyNone_New()
    catch err
        PyErr_SetJuliaError(err)
        PyNULL
    end
end
pyjlset_update(xo, yo) = pyjlset_ibinop_named(xo, yo, union!, false)
pyjlset_difference_update(xo, yo) = pyjlset_ibinop_named(xo, yo, setdiff!, true)
pyjlset_symmetric_difference_update(xo, yo) = pyjlset_ibinop_named(xo, yo, symdiff!, false)
pyjlset_intersection_update(xo, yo) = pyjlset_ibinop_named(xo, yo, intersect!, true)

# pyjlset_ibinop_operator(xo, yo, op, skip) = begin
#     r = PySetABC_Check(yo)
#     r == -1 && return PyNULL
#     r ==  0 && return PyNotImplemented_New()
#     x = PyJuliaValue_GetValue(xo)::AbstractSet
#     y = PyIterable_Collect(yo, eltype(x), skip)
#     isempty(y) && PyErr_IsSet() && return PyNULL
#     try
#         PyObject_From(op(x, y))
#     catch err
#         PyErr_SetJuliaError(err)
#         PyNULL
#     end
# end
# pyjlset_ior(xo, yo) = pyjlset_ibinop_operator(xo, yo, union!, false)
# pyjlset_isub(xo, yo) = pyjlset_ibinop_operator(xo, yo, setdiff!, true)
# pyjlset_ixor(xo, yo) = pyjlset_ibinop_operator(xo, yo, symdiff!, false)
# pyjlset_iand(xo, yo) = pyjlset_ibinop_operator(xo, yo, intersect!, true)

pyjlset_binop_generic(xo::PyPtr, yo::PyPtr, op) = begin
    xo2 = PySet_New(xo)
    isnull(xo2) && return PyNULL
    yo2 = PySet_New(yo)
    isnull(yo2) && (Py_DecRef(xo2); return PyNULL)
    ro = op(xo2, yo2)
    Py_DecRef(xo2)
    Py_DecRef(yo2)
    ro
end
pyjlset_or_generic(xo, yo) = pyjlset_binop_generic(xo, yo, PyNumber_InPlaceOr)
pyjlset_xor_generic(xo, yo) = pyjlset_binop_generic(xo, yo, PyNumber_InPlaceXor)
pyjlset_and_generic(xo, yo) = pyjlset_binop_generic(xo, yo, PyNumber_InPlaceAnd)
pyjlset_sub_generic(xo, yo) = pyjlset_binop_generic(xo, yo, PyNumber_InPlaceSubtract)

pyjlset_binop_special(xo::PyPtr, yo::PyPtr, op) = begin
    x = PyJuliaValue_GetValue(xo)::AbstractSet
    y = PyJuliaValue_GetValue(yo)::AbstractSet
    try
        PyObject_From(op(x, y))
    catch err
        PyErr_SetJuliaError(err)
        PyNULL
    end
end
pyjlset_or_special(xo, yo) = pyjlset_binop_special(xo, yo, union)
pyjlset_xor_special(xo, yo) = pyjlset_binop_special(xo, yo, symdiff)
pyjlset_and_special(xo, yo) = pyjlset_binop_special(xo, yo, intersect)
pyjlset_sub_special(xo, yo) = pyjlset_binop_special(xo, yo, setdiff)

# pyjlset_binop_operator(xo::PyPtr, yo::PyPtr, gop, sop) = begin
#     t = PyJuliaSetValue_Type()
#     isnull(t) && return PyNULL
#     if (r=PyObject_IsInstance(xo, t); ism1(r) && return PyNULL; r!=0)
#         if (r=PyObject_IsInstance(yo, t); ism1(r) && return PyNULL; r!=0)
#             sop(xo, yo)
#         elseif (r=PySetABC_Check(yo); ism1(r) && return PyNULL; r!=0)
#             gop(xo, yo)
#         else
#             PyNotImplemented_New()
#         end
#     elseif (r=PyObject_IsInstance(yo, t); ism1(r) && return PyNULL; r!=0)
#         if (r=PySetABC_Check(xo); ism1(r) && return PyNULL; r!=0)
#             gop(xo, yo)
#         else
#             PyNotImplemented_New()
#         end
#     else
#         PyNotImplemented_New()
#     end
# end
# pyjlset_or(xo, yo) = pyjlset_binop_operator(xo, yo, pyjlset_or_generic, pyjlset_or_special)
# pyjlset_xor(xo, yo) = pyjlset_binop_operator(xo, yo, pyjlset_xor_generic, pyjlset_xor_special)
# pyjlset_and(xo, yo) = pyjlset_binop_operator(xo, yo, pyjlset_and_generic, pyjlset_and_special)
# pyjlset_sub(xo, yo) = pyjlset_binop_operator(xo, yo, pyjlset_sub_generic, pyjlset_sub_special)

pyjlset_binop_named(xo, yo, gop, sop) = begin
    t = PyJuliaSetValue_Type()
    isnull(t) && return PyNULL
    if (r = PyObject_IsInstance(yo, t); ism1(r) && return PyNULL; r != 0)
        sop(xo, yo)
    elseif (r = PySetABC_Check(yo); ism1(r) && return PyNULL; r != 0)
        gop(xo, yo)
    else
        PyNotImplemented_New()
    end
end
pyjlset_union(xo, yo) = pyjlset_binop_named(xo, yo, pyjlset_or_generic, pyjlset_or_special)
pyjlset_symmetric_difference(xo, yo) =
    pyjlset_binop_named(xo, yo, pyjlset_xor_generic, pyjlset_xor_special)
pyjlset_intersection(xo, yo) =
    pyjlset_binop_named(xo, yo, pyjlset_and_generic, pyjlset_and_special)
pyjlset_difference(xo, yo) =
    pyjlset_binop_named(xo, yo, pyjlset_sub_generic, pyjlset_sub_special)

pyjlset_isdisjoint(x::PyPtr, y::PyPtr) =
    _pyjlset_isdisjoint(PyJuliaValue_GetValue(x)::AbstractSet, y)
_pyjlset_isdisjoint(x::AbstractSet, yso::PyPtr) = begin
    isempty(x) && return PyObject_From(true)
    r = PyIterable_Map(yso) do yo
        r = PyObject_TryConvert(yo, eltype(x))
        r == -1 && return -1
        r == 0 && return 1
        y == takeresult(eltype(x))
        y in x ? 0 : 1
    end
    r == -1 ? PyNULL : r == 0 ? PyObject_From(false) : PyObject_From(true)
end


text:
using Documenter

push!(LOAD_PATH,"../src/")

using BerghAlgorithms

DocMeta.setdocmeta!(BerghAlgorithms, :DocTestSetup, :(using BerghAlgorithms); recursive = true)

makedocs(
    modules=[BerghAlgorithms],
    sitename="Bergh Documentation",
    pages = [
        "Home" => "index.md",
        "Functions" => "functions.md",
        "Contributors" => "contributors.md"
    ]
)

deploydocs(;
    repo="github.com/nbenabla/destackification",
)

max_stars_repo_name: nbenabla/destackification

max_issues_repo_name: nbenabla/destackification

max_forks_repo_name: nbenabla/destackification


text:
# # Event counting with Negative Binomial

# We load all the necessary packages
using AbstractGPs
using ApproximateGPs
using AugmentedGPLikelihoods
using Distributions
using LinearAlgebra
# Plotting libraries
using Plots

# We create some random data (sorted for plotting reasons)
N = 100
x = range(-10, 10; length=N)
kernel = with_lengthscale(SqExponentialKernel(), 2.0)
gp = GP(kernel)
lik = NegativeBinomialLikelihood(NBParamFailure(15))
lf = LatentGP(gp, lik, 1e-6)
f, y = rand(lf(x));
# We plot the sampled data
plt = scatter(x, y; label="Data")
plot!(plt, x, f; color=:red, label="Latent GP")
# ## CAVI Updates
# We write our CAVI algorithmm
function u_posterior(fz, m, S)
    return posterior(SparseVariationalApproximation(Centered(), fz, MvNormal(m, S)))
end

function cavi!(fz::AbstractGPs.FiniteGP, x, y, m, S, qΩ; niter=10)
    K = ApproximateGPs._chol_cov(fz)
    for _ in 1:niter
        post_u = u_posterior(fz, m, S)
        post_fs = marginals(post_u(x))
        aux_posterior!(qΩ, lik, y, post_fs)
        S .= inv(Symmetric(inv(K) + Diagonal(only(expected_auglik_precision(lik, qΩ, y)))))
        m .= S * (only(expected_auglik_potential(lik, qΩ, y)) + K \ mean(fz))
    end
    return m, S
end;
# Now we just initialize the variational parameters
m = zeros(N)
S = Matrix{Float64}(I(N))
qΩ = init_aux_posterior(lik, N)
fz = gp(x, 1e-8);
# And visualize the current posterior
x_te = -10:0.01:10
plot!(
    plt, x_te, u_posterior(fz, m, S); color=:blue, alpha=0.3, label="Initial VI Posterior"
)
# We run CAVI for 3-4 iterations
cavi!(fz, x, y, m, S, qΩ; niter=4);
# And visualize the obtained variational posterior
plot!(
    plt,
    x_te,
    u_posterior(fz, m, S);
    color=:darkgreen,
    alpha=0.3,
    label="Final VI Posterior",
)
# ## ELBO
# How can one compute the Augmented ELBO?
# Again AugmentedGPLikelihoods provides helper functions
# to not have to compute everything yourself
function aug_elbo(lik, u_post, x, y)
    qf = marginals(u_post(x))
    qΩ = aux_posterior(lik, y, qf)
    return expected_logtilt(lik, qΩ, y, qf) - aux_kldivergence(lik, qΩ, y) -
           kldivergence(u_post.approx.q, u_post.approx.fz)     # approx.fz is the prior and approx.q is the posterior
end

aug_elbo(lik, u_posterior(fz, m, S), x, y)
# ## Gibbs Sampling
# We create our Gibbs sampling algorithm (we could do something fancier with
# AbstractMCMC)
function gibbs_sample(fz, f, Ω; nsamples=200)
    K = ApproximateGPs._chol_cov(fz)
    Σ = zeros(length(f), length(f))
    μ = zeros(length(f))
    return map(1:nsamples) do _
        aux_sample!(Ω, lik, y, f)
        Σ .= inv(Symmetric(inv(K) + Diagonal(only(auglik_precision(lik, Ω, y)))))
        μ .= Σ * (only(auglik_potential(lik, Ω, y)) + K \ mean(fz))
        rand!(MvNormal(μ, Σ), f)
        return copy(f)
    end
end;
# We initialize our random variables
f = randn(N)
Ω = init_aux_variables(lik, N);
# Run the sampling for default number of iterations (200)
fs = gibbs_sample(fz, f, Ω);
# And visualize the samples overlapped to the variational posterior
# that we found earlier.
for f in fs
    plot!(plt, x, f; color=:black, alpha=0.07, label="")
end
plt


text:
# To run only selected tests:
# $julia -e 'include("test/runtests.jl")' TEST1...

# To run during development (from REPL) include this file INSIDE A MODULE

using Pseudospectra, Test

tests = ["basic","arpack1","medcplx","medrect","smrect","spdirect",
          "projection_and_transient","radius_abscissa","negproj","eigsnoconv",
          "numrange", "linmap", "big", "power_transient", "swmethod", "zoom"
          ]

if length(ARGS) > 0
    tests = ARGS
end

include("testplotter.jl")

for f in tests
    println("Testing $f functionality")
    fnam = f*".jl"
    include(fnam)
end


text:
# an example script on how to run regression tests
# best do it manually for now, until we learn how to use
# this well and then setup an automated workflow

using ACE, PkgBenchmark

results = benchmarkpkg(ACE)
writeresults(@__DIR__() * "/results.json", results)

results_main = benchmarkpkg(ACE, "main")
writeresults(@__DIR__() * "/results_main.json", results)

results = readresults(@__DIR__() * "/results.json")
results_base = readresults(@__DIR__() * "/results_main.json")
judgement = judge(results, results_main)
export_markdown(@__DIR__() * "/judge.md", judgement)


text:
"""
    solve_hadc(model, experiment) -> NamedTuple

Compute ADCs using the homogenized ADC model (HADC).
"""
function solve_hadc(model, experiment)

    # Measure function evalutation time
    starttime = Base.time()

    # Extract input parameters
    @unpack mesh, D, T₂, ρ = model
    @unpack directions, sequences = experiment.gradient
    @unpack odesolver, reltol, abstol = experiment.hadc

    # Deduce sizes
    ncompartment, nboundary = size(mesh.facets)
    npoint_cmpts = size.(mesh.points, 2)
    inds_cmpts = cumsum([0; npoint_cmpts])
    ndirection = size(directions, 2)
    nsequence = length(sequences)

    # Number of points in each compartment
    npoint_cmpts = size.(mesh.points, 2)

    # Assemble finite element matrices
    M = []
    S = []
    G = []
    volumes = zeros(ncompartment)
    for icmpt = 1:ncompartment
        points = mesh.points[icmpt]
        elements = mesh.elements[icmpt]
        fevolumes, = get_mesh_volumes(points, elements)
        volumes[icmpt] = sum(fevolumes)

        # Assemble finite element matrices
        push!(M, assemble_mass_matrix(elements', fevolumes))
        push!(S, assemble_stiffness_matrix(elements', points', D[icmpt]))

        # Compute surface integrals
        push!(G, zeros(npoint_cmpts[icmpt], 3))
        for iboundary = 1:nboundary
            facets = mesh.facets[icmpt, iboundary]
            if !isempty(facets)
                # Get facet normals
                _, _, normals = get_mesh_surfacenormals(points, elements, facets)

                # Surface normal weighted flux matrix (in each canonical direction)
                for dim = 1:3
                    Q = assemble_flux_matrix(facets', points', normals[dim, :])
                    G[icmpt][:, dim] += sum(Q, dims = 2)
                end
            end
        end
    end

    # Initial conditions
    ω₀ = zeros.(npoint_cmpts)

    # Time dependent ODE function
    function Mdω!(dω, ω, p, t)
        # @show t
        @unpack mS, f, surfint = p
        mul!(dω, mS, ω)
        dω .+= integral(f, t) .* surfint
    end

    # Allocate output arrays
    adc = zeros(ncompartment, nsequence, ndirection)
    adc_allcmpts = zeros(nsequence, ndirection)
    itertimes = zeros(ncompartment, nsequence, ndirection)

    # Iterate over compartments and gradient sequences and directions
    for idir = 1:ndirection, iseq = 1:nsequence, icmpt = 1:ncompartment

        # Measure iteration time
        itertime = time()

        # Extract parameters for iteration
        g = directions[:, idir]
        f = sequences[iseq]
        TE = echotime(f)

        # Free diffusivity in direction g
        D₀ = g' * D[icmpt] * g

        # Compute surface integrals in gradient direction
        surfint = G[icmpt] * (D[icmpt] * g)

        # Display state of iterations
        println("Solving HADC model of size $(sum(npoint_cmpts)):")
        println("  Direction   $idir of $ndirection: g = $g")
        println("  Sequence    $iseq of $nsequence: f = $f")
        println("  Compartment $icmpt of $ncompartment")

        # Create ODE function and Jacobian from matrices
        jac = (J, _, p, t) -> (J .= p.mS)
        p = (; mS = -S[icmpt], f, surfint)
        odefunction =
            ODEFunction(Mdω!, jac = jac, jac_prototype = p.mS, mass_matrix = M[icmpt])
        odeproblem = ODEProblem(odefunction, ω₀[icmpt], (0, TE), p, progress = false)

        # Solve ODE, keep all time steps (for integral)
        sol = solve(odeproblem, odesolver, reltol = reltol, abstol = abstol)

        # Integral over compartment boundary
        a, err = quadgk(t -> integral(f, t) * (surfint' * sol(t)), 0, TE)

        # HADC (free diffusivity minus correction)
        adc[icmpt, iseq, idir] = D₀ - a / volumes[icmpt] / bvalue_no_q(f)

        # Computational time
        itertimes[icmpt, iseq, idir] = time() - itertime
    end

    # Compute total HADC (weighted sum over compartments)
    weights = ρ .* volumes
    weights = weights / sum(weights)
    adc_allcmpts[:] = sum(weights .* adc, dims = 1)

    # Create output structure
    totaltime = time() - starttime

    (; adc, adc_allcmpts, totaltime, itertimes)
end


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct IoK8sApiBatchV1JobStatus <: SwaggerModel
    active::Any # spec type: Union{ Nothing, Int32 } # spec name: active
    completionTime::Any # spec type: Union{ Nothing, IoK8sApimachineryPkgApisMetaV1Time } # spec name: completionTime
    conditions::Any # spec type: Union{ Nothing, Vector{IoK8sApiBatchV1JobCondition} } # spec name: conditions
    failed::Any # spec type: Union{ Nothing, Int32 } # spec name: failed
    startTime::Any # spec type: Union{ Nothing, IoK8sApimachineryPkgApisMetaV1Time } # spec name: startTime
    succeeded::Any # spec type: Union{ Nothing, Int32 } # spec name: succeeded

    function IoK8sApiBatchV1JobStatus(;active=nothing, completionTime=nothing, conditions=nothing, failed=nothing, startTime=nothing, succeeded=nothing)
        o = new()
        validate_property(IoK8sApiBatchV1JobStatus, Symbol("active"), active)
        setfield!(o, Symbol("active"), active)
        validate_property(IoK8sApiBatchV1JobStatus, Symbol("completionTime"), completionTime)
        setfield!(o, Symbol("completionTime"), completionTime)
        validate_property(IoK8sApiBatchV1JobStatus, Symbol("conditions"), conditions)
        setfield!(o, Symbol("conditions"), conditions)
        validate_property(IoK8sApiBatchV1JobStatus, Symbol("failed"), failed)
        setfield!(o, Symbol("failed"), failed)
        validate_property(IoK8sApiBatchV1JobStatus, Symbol("startTime"), startTime)
        setfield!(o, Symbol("startTime"), startTime)
        validate_property(IoK8sApiBatchV1JobStatus, Symbol("succeeded"), succeeded)
        setfield!(o, Symbol("succeeded"), succeeded)
        o
    end
end # type IoK8sApiBatchV1JobStatus

const _property_map_IoK8sApiBatchV1JobStatus = Dict{Symbol,Symbol}(Symbol("active")=>Symbol("active"), Symbol("completionTime")=>Symbol("completionTime"), Symbol("conditions")=>Symbol("conditions"), Symbol("failed")=>Symbol("failed"), Symbol("startTime")=>Symbol("startTime"), Symbol("succeeded")=>Symbol("succeeded"))
const _property_types_IoK8sApiBatchV1JobStatus = Dict{Symbol,String}(Symbol("active")=>"Int32", Symbol("completionTime")=>"IoK8sApimachineryPkgApisMetaV1Time", Symbol("conditions")=>"Vector{IoK8sApiBatchV1JobCondition}", Symbol("failed")=>"Int32", Symbol("startTime")=>"IoK8sApimachineryPkgApisMetaV1Time", Symbol("succeeded")=>"Int32")
Base.propertynames(::Type{ IoK8sApiBatchV1JobStatus }) = collect(keys(_property_map_IoK8sApiBatchV1JobStatus))
Swagger.property_type(::Type{ IoK8sApiBatchV1JobStatus }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_IoK8sApiBatchV1JobStatus[name]))}
Swagger.field_name(::Type{ IoK8sApiBatchV1JobStatus }, property_name::Symbol) =  _property_map_IoK8sApiBatchV1JobStatus[property_name]

function check_required(o::IoK8sApiBatchV1JobStatus)
    true
end

function validate_property(::Type{ IoK8sApiBatchV1JobStatus }, name::Symbol, val)
end


text:
function part1(input)
    fields, your_ticket, nearby_tickets = parse_input(input)
    valid_values = find_valid_values(fields)
    return sum(value
               for value in reduce(vcat, nearby_tickets)
               if value ∉ valid_values)
end

function part2(input)
    fields, your_ticket, nearby_tickets = parse_input(input)
    valid_values = find_valid_values(fields)
    filter!(numbers -> all(number in valid_values for number in numbers),
            nearby_tickets)
    possible_assignments = Dict(key => find_field_number(nearby_tickets, value)
                                for (key, value) in fields)
    assignments = Dict{String, Int}()
    while !isempty(possible_assignments)
        assigned = nothing
        for (key, value) in possible_assignments
            if length(value) == 1
                assigned = only(value)
                assignments[key] = assigned
                delete!(possible_assignments, key)
                break
            end
        end
        @assert !isnothing(assigned)
        for value in values(possible_assignments)
            filter!(!=(assigned), value)
        end
    end

    return prod(Int[your_ticket[value]
                    for (key, value) in assignments
                    if startswith(key, "departure")])
end

function find_field_number(nearby_tickets, ranges)
    possible_fields = Set(1:length(first(nearby_tickets)))
    for ticket in nearby_tickets
        bad_fields = Set{Int}()
        for field in possible_fields
            if !any(ticket[field] in range for range in ranges)
                push!(bad_fields, field)
            end
            possible_fields = setdiff(possible_fields, bad_fields)
        end
    end
    return possible_fields
end

function find_valid_values(fields)
    valid_values = Set{Int}()
    for range in reduce(vcat, values(fields))
        for i in range
            i in valid_values || push!(valid_values, i)
        end
    end
    return valid_values
end

function parse_input(input)
    fields = Dict{String, Vector{UnitRange{Int}}}()
    field_data, your_ticket, nearby_tickets = split(read(input, String), "\n\n")
    for line in split(field_data, "\n")
        name, a, b, c, d = match(r"([\w ]+): (\d+)-(\d+) or (\d+)-(\d+)",
                                 line).captures
        fields[name] = [parse(Int, a):parse(Int, b),
                        parse(Int, c):parse(Int, d)]
    end

    your = get_numbers(last(split(your_ticket, "\n")))

    nearby = get_numbers.(split(nearby_tickets, "\n", keepempty = false)[2:end])

    return fields, your, nearby
end

get_numbers(line) = parse.(Int, split(line, ","))


text:
using Drawing
using ColorTypes

ignore = false  # when re-generating

r1, r2, r3, n, w = 1, 1, 0.5, 200, 0.5
d1, d2, rotn = 120, -140, 40
sat, val, alpha = 0.7, 0.5, 0.1
side = 600

width = 2*pi*r1 / n * w
srand(0)

with(
     PNG("colour-ring.png", side, side),
     Axes(negative=true),
     Translate(0.9, 0.9), Scale(3),
     Paper("black"),
     Pen(width; cap="round")) do
     for x in shuffle!([i / n for i in 1:n])
         hue = x * 360
         theta1 = deg2rad(rotn + x*360)
         theta2 = theta1 + deg2rad(d1)
         theta3 = theta1 + deg2rad(d2)
         draw(Ink(HSVA(hue, sat, val + (1-val)*rand(), alpha))) do
             move(r1 * cos(-theta1), r1  * sin(-theta1))
             line(r2 * cos(-theta2), r2  * sin(-theta2))
             move(r1 * cos(-theta1), r1  * sin(-theta1))
             line(r3 * cos(-theta3), r3  * sin(-theta3))
         end
     end
end
ignore || compare("colour-ring.png")


text:
# Model
include_model("particle")

# Horizon
T = 201

# Time step
tf = 1.0
h = tf / (T-1)

# Bounds
_uu = Inf * ones(model.m)
_uu[model.idx_u] .= Inf
_ul = zeros(model.m)
_ul[model.idx_u] .= -Inf
ul, uu = control_bounds(model, T, _ul, _uu)

# Circle
θ = range(0, stop = 2 * π, length = T)
x = cos.(θ)
y = sin.(θ)
plot(x, y, ratio = :equal)

# Initial and final states
q_ref = [[x[t]; y[t]; 1.0] for t = 1:T]
x_ref = configuration_to_state(q_ref)

xl, xu = state_bounds(model, T,
    x1 = [q_ref[1]; q_ref[1]],
    xT = [q_ref[T]; q_ref[T]])

# Objective
obj_penalty = PenaltyObjective(1.0e5, model.m)
obj_track = quadratic_tracking_objective(
    [Diagonal(1000.0 * ones(model.n)) for t = 1:T],
    [Diagonal([1.0e-1 * ones(model.nu); zeros(model.m - model.nu)]) for t = 1:T-1],
    [x_ref[t] for t = 1:T],
    [zeros(model.m) for t = 1:T])
obj = MultiObjective([obj_penalty, obj_track])

# Constraints
include_constraints("contact")
con_contact = contact_constraints(model, T)

# Problem
prob = trajectory_optimization_problem(model,
               obj,
               T,
               h = h,
               xl = xl,
               xu = xu,
               ul = ul,
               uu = uu,
               con = con_contact)

# Trajectory initialization
x0 = deepcopy(x_ref) #linear_interpolation(x1, x1, T) # linear interpolation on state
u0 = [1.0e-5 * rand(model.m) for t = 1:T-1] # random controls

# Pack trajectories into vector
z0 = pack(x0, u0, prob)

#NOTE: may need to run examples multiple times to get good trajectories
# Solve nominal problem
@time z̄ , info = solve(prob, copy(z0), tol = 1.0e-4, c_tol = 1.0e-4, mapl = 5)

check_slack(z̄, prob)
x̄, ū = unpack(z̄, prob)

Q = [Diagonal(ones(model.n)) for t = 1:T]
R = [Diagonal(1.0 * ones(model.m)) for t = 1:T-1]
x_proj, u_proj = lqr_projection(model, x̄, ū, h, Q, R)
q_proj = state_to_configuration(x_proj)

using Plots
plot(hcat(ū...)[1:3, :]',
    linetype = :steppost,
    label = "",
    color = :red,
    width = 2.0)

plot!(hcat(u_proj...)[1:3, :]',
    linetype = :steppost,
    label = "", color = :black)

plot(hcat(state_to_configuration(x̄)...)',
    color = :red,
    width = 2.0,
	label = "")

plot!(hcat(state_to_configuration(x_proj)...)',
    color = :black,
	label = "")

q̄ = state_to_configuration(x̄)

include(joinpath(pwd(), "models/visualize.jl"))
vis = Visualizer()
render(vis)
visualize!(vis, model, state_to_configuration(x̄), Δt = h)

ee_positions = hcat(q_ref...)
points = collect(eachcol(ee_positions))
material = LineBasicMaterial(color=colorant"orange", linewidth=6.0)
setobject!(vis["ee path"], Object(PointCloud(points), material, "Line"))

# Simulator
function dynamics(model::Particle, v1, q1, q2, q3, u, λ, b, w, h)
	(M_func(model, q1) * v1
	- M_func(model, q2) * (SVector{3}(q3) - SVector{3}(q2)) / h
	+ h * (transpose(B_func(model, q3)) * SVector{3}(u)
	+ transpose(N_func(model, q3)) * SVector{1}(λ)
	+ transpose(P_func(model, q3)) * SVector{4}(b)
	- G_func(model, q2)))
end

function maximum_dissipation(model::Particle, q2, q3, ψ, η, h)
	ψ_stack = ψ[1] * ones(model.nb)
	return P_func(model, q3) * (q3 - q2) / h + ψ_stack - η
end

function no_slip(model::Particle, q2, q3, λ, h)
	return (λ' * _P_func(model, q3) * (q3 - q2) / h)[1]
end

function friction_cone(model::Particle, λ, b)
	return @SVector [model.μ * λ[1] - sum(b)]
end

mutable struct MOISimulator <: MOI.AbstractNLPEvaluator
    num_var::Int                 # number of decision variables
    num_con::Int                 # number of constraints
    primal_bounds
    constraint_bounds

	model
	slack_penalty
	v1
	q1
	q2
	u
	w
	h
end

function simulator_problem(model, v1, q1, q2, u, w, h; slack_penalty = 1.0e5)
	num_var = model.nq + model.nc + model.nb + model.nc + model.nb + model.ns
	num_con = model.nq + model.nc + model.nc + model.nb +  3

	zl = zeros(num_var)
	zl[1:model.nq] .= -Inf
	zu = Inf * ones(num_var)

	cl = zeros(num_con)
	cu = zeros(num_con)
	cu[model.nq .+ (1:model.nc)] .= Inf
	cu[model.nq + model.nc .+ (1:model.nc)] .= Inf
	cu[model.nq + model.nc + model.nc + model.nb + 1] = Inf
	cu[model.nq + model.nc + model.nc + model.nb + 2] = Inf
	cu[model.nq + model.nc + model.nc + model.nb + 3] = Inf

	MOISimulator(num_var,
		num_con,
		(zl, zu),
		(cl, cu),
		model,
		slack_penalty,
		v1,
		q1,
		q2,
		u,
		w,
		h)
end
prob_sim = simulator_problem(model, ones(nq), ones(nq), ones(nq), zeros(nu), ones(nq), h)
z0 = rand(prob_sim.num_var)

function MOI.eval_objective(prob::MOISimulator, x)
    return prob.slack_penalty * x[prob.num_var]
end

MOI.eval_objective(prob_sim, z0)

function MOI.eval_objective_gradient(prob::MOISimulator, grad_f, x)
    grad_f .= 0.0
	grad_f[prob.num_var] = prob.slack_penalty
	return nothing
end

∇obj = zeros(prob_sim.num_var)
MOI.eval_objective_gradient(prob_sim, ∇obj, z0)

function MOI.eval_constraint(prob::MOISimulator, g, x)
	model = prob.model
	q3 = view(x, 1:model.nq)
	λ = view(x, model.nq .+ (1:model.nc))
	b = view(x, model.nq + model.nc .+ (1:model.nb))
	ψ = view(x, model.nq + model.nc + model.nb .+ (1:model.nc))
	η = view(x, model.nq + model.nc + model.nb + model.nc .+ (1:model.nb))
	s = view(x, model.nq + model.nc + model.nb + model.nc + model.nb .+ (1:model.ns))
	# λ = 0
	# b = 0
    g[1:model.nq] = dynamics(model,
		prob.v1, prob.q1, prob.q2, q3, prob.u, λ, b, prob.w, prob.h)
	g[model.nq .+ (1:model.nc)] = ϕ_func(model, q3)
	g[model.nq + model.nc .+ (1:model.nc)] = friction_cone(model, λ, b)
	g[model.nq + model.nc + model.nc .+ (1:model.nb)] = maximum_dissipation(model,
		prob.q2, q3, ψ, η, prob.h)
	g[model.nq + model.nc + model.nc + model.nb + 1] = s[1] - (λ' * ϕ_func(model, q3))[1]
	g[model.nq + model.nc + model.nc + model.nb + 2] = s[1] - (ψ' * friction_cone(model, λ, b))[1]
	g[model.nq + model.nc + model.nc + model.nb + 3] = s[1] - (η' * b)[1]

    return nothing
end
prob_sim.h
c0 = zeros(prob_sim.num_con)
MOI.eval_constraint(prob_sim, c0, z0)

function MOI.eval_constraint_jacobian(prob::MOISimulator, jac, x)
    con!(g, z) = MOI.eval_constraint(prob, g, z)
	ForwardDiff.jacobian!(reshape(jac, prob.num_var, prob.num_con), con!, zeros(prob.num_con), x)
    return nothing
end

∇c0 = vec(zeros(prob_sim.num_con, prob_sim.num_var))
MOI.eval_constraint_jacobian(prob_sim, ∇c0, z0)

function sparsity_jacobian(prob::MOISimulator)
    row = []
    col = []

    r = (1:prob.num_con)
    c = (1:prob.num_var)

    row_col!(row, col, r, c)

    return collect(zip(row, col))
end

sparsity_jacobian(prob_sim)

@time solve(prob_sim, z0)

function step_contact(model, v1, q1, q2, u, w, h)
    prob = simulator_problem(model, v1, q1, q2, u, w, h)
    z0 = [copy(q2); 1.0e-5 * ones(model.nc + model.nb + model.nc + model.nb + model.ns)]
    @time z , info = solve(prob, copy(z0), tol = 1.0e-6, c_tol = 1.0e-6, mapl = 0)

    @assert z[end] < 1.0e-5

    return z[1:model.nq]
end

q_sim = [q_proj[1], q_proj[2]]
v_sim = [(q_proj[2] - q_proj[1]) / h]

for t = 1:T-1
	# d x rate
	_q_sim = [q_sim[end-1], q_sim[end]]
	_v_sim = [v_sim[end]]

	d = 10
	for i = 1:d
		_h = h / convert(Float64, d)
		_q = step_contact(model,
			_v_sim[end], _q_sim[end-1], _q_sim[end], u_proj[t][model.idx_u], zeros(model.d), _h)
		_v = (_q - _q_sim[end]) / _h

		push!(_q_sim, _q)
		push!(_v_sim, _v)
	end

	push!(q_sim, _q_sim[end])
	push!(v_sim, _v_sim[end])

	# same rate
	# push!(q_sim, step_contact(model,
	# 	v_sim[end], q_sim[end-1], q_sim[end], ū[t][model.idx_u], zeros(model.d), h))
	# push!(v_sim, (q_sim[end] - q_sim[end-1]) / h)
end

plot(hcat(q_proj...)',
    color = :red,
    width = 2.0,
	label = "")
plot!(hcat(q_sim...)',
    color = :black,
    width = 1.0,
	label = "")

plot(hcat(v_sim...)',
    color = :black,
    width = 1.0,
	label = "")