text:
const OrderedEinCode{LT} = Union{NestedEinsum{DynamicEinCode{LT}}, SlicedEinsum{LT,NestedEinsum{DynamicEinCode{LT}}}}

# we implement the register interface because we want to use the operator system in Yao.
abstract type PEPS{T,NF,LT} <:AbstractRegister{NF} end

"""
    SimplePEPS{T,NF,LT<:Union{Int,Char},Ein<:OrderedEinCode} <: PEPS{T,NF,LT}

 a     b     c     d
 ┆     ┆     ┆     ┆ ← size = nflavor
 ●--α--●--β--●--γ--●
                ↑
            size ≤ Dmax

`α`, `β` and `γ` are virtual labels.
`a`, `b`, `c` and `d` are physical labels.
`aα`, `αbβ`, `βcγ` and `γd` are vertex (tensor) labels.

`NF` is the number of flavors, i.e. size of physical bonds. For spins, it is 2.

It has fields

* `physical_labels` is a vector of unique physical labels, should be a vector of integers.
* `virtual_labels` is a vector if unique virtual labels, should be a vector of integer 2-tuples.

* `vertex_labels` is a vector of vectors, i-th vector is the labels for i-th vertex tensor.
* `vertex_tensors` is a vector of tensors defined on vertices.
* `max_index` is the maximum index, used for creating new labels.

* `code_statetensor` is the optimized contraction code for obtaining state vector.
* `code_inner_product` is the optimized contraction code for obtaining the overlap.

* `Dmax` is the maximum virtual bond dimension.
* `ϵ` is useful in compression (e.g. with SVD), to determine the cutoff precision.
"""
struct SimplePEPS{T,NF,LT<:Union{Int,Char},Ein<:OrderedEinCode} <: PEPS{T,NF,LT}
    physical_labels::Vector{LT}
    virtual_labels::Vector{LT}

    vertex_labels::Vector{Vector{LT}}
    vertex_tensors::Vector{<:AbstractArray{T}}
    max_index::LT

    # optimized contraction codes
    code_statetensor::Ein
    code_inner_product::Ein

    Dmax::Int
    ϵ::Float64
end

function SimplePEPS{NF}(physical_labels::Vector{LT},
    virtual_labels::Vector{LT},

    vertex_labels::Vector{Vector{LT}},
    vertex_tensors::Vector{<:AbstractArray{T}},
    max_index::LT,

    code_statetensor::Ein,
    code_inner_product::Ein,

    Dmax::Int,
    ϵ::Float64
    ) where {LT,T,NF,Ein}
    SimplePEPS{T,NF,LT,Ein}(physical_labels, virtual_labels, vertex_labels, vertex_tensors, max_index,
        code_statetensor, code_inner_product, Dmax, ϵ)
end

function SimplePEPS{NF}(vertex_labels::AbstractVector{<:AbstractVector{LT}}, vertex_tensors::Vector{<:AbstractArray{T}},
        virtual_labels::AbstractVector{LT}, Dmax::Int, ϵ::Real, optimizer::CodeOptimizer, simplifier::CodeSimplifier) where {LT,T,NF}
    physical_labels = [vl[findall(∉(virtual_labels), vl)[]] for vl in vertex_labels]
    # check dimension of input tensors
    if !all(tv->all(tvi->(tvi[2] ∈ physical_labels ? tvi[1] == NF : tvi[1] <= Dmax), zip(size(tv[1]), tv[2])), zip(vertex_tensors,vertex_labels))
        error("Input tensor shape error, either: \n 1. physical degrees of freedom do not match, \n 2. or input tensor virtual degrees of freedom larger than `Dmax = $(Dmax)`.")
    end
    max_ind = max(maximum(physical_labels), maximum(virtual_labels))
    # optimal contraction orders
    optcode_statetensor, optcode_inner_product = _optimized_code(vertex_labels, physical_labels, virtual_labels,
        max_ind, NF, Dmax, optimizer, simplifier)
    SimplePEPS{NF}(physical_labels, virtual_labels, vertex_labels, vertex_tensors, max_ind,
        optcode_statetensor, optcode_inner_product, Dmax, ϵ)
end

function _optimized_code(alllabels, physical_labels::AbstractVector{LT}, virtual_labels, max_ind, nflavor, D, optimizer, simplifier) where LT
    code_statetensor = EinCode(alllabels, physical_labels)
    size_dict = Dict([[l=>nflavor for l in physical_labels]..., [l=>D for l in virtual_labels]...])
    optcode_statetensor = optimize_code(code_statetensor, size_dict, optimizer, simplifier)
    rep = [l=>max_ind+i for (i, l) in enumerate(virtual_labels)]
    merge!(size_dict, Dict([l.second=>D for l in rep]))
    code_inner_product = EinCode([alllabels..., [replace(l, rep...) for l in alllabels]...], LT[])
    optcode_inner_product = optimize_code(code_inner_product, size_dict, optimizer, simplifier)
    return optcode_statetensor, optcode_inner_product
end

# all labels for vertex tensors and bond tensors (if any)
alllabels(s::SimplePEPS) = s.vertex_labels
# all vertex tensors and bond tensors (if any)
alltensors(s::SimplePEPS) = s.vertex_tensors
Base.copy(peps::SimplePEPS{T,NF}) where {T,NF} = SimplePEPS{NF}(copy(peps.physical_labels), copy(peps.virtual_labels),
    copy(peps.vertex_labels), copy(peps.vertex_tensors), peps.max_index,
    peps.code_statetensor, peps.code_inner_product, peps.Dmax, peps.ϵ)

# ●----●----●----●   ← ⟨peps1|
# ┆    ┆    ┆    ┆
# ●----●----●----●   ← |peps2⟩
function inner_product(p1::PEPS, p2::PEPS)
    p1c = conj(p1)
    # we assume `p1` and `p2` have the same structure and virtual bond dimension.
    p1.code_inner_product(alltensors(p1c)..., alltensors(p2)...)
end

"""
    VidalPEPS{T, NF, LT<:Union{Int,Char}} <: PEPS{T,NF,LT}

Similar to `SimplePEPS` except it contains one extra field

 a     b     c     d
 ┆  α  ┆  β  ┆  γ  ┆ ← size = nflavor
 ●--◆--●--◆--●--◆--●
                ↑
            bond tensor

* `bond_tensors` is a vector of tensors on bonds.
"""
struct VidalPEPS{T, NF, LT<:Union{Int,Char}, Ein<:OrderedEinCode} <: PEPS{T,NF,LT}
    physical_labels::Vector{LT}
    virtual_labels::Vector{LT}

    vertex_labels::Vector{Vector{LT}}
    vertex_tensors::Vector{<:AbstractArray{T}}
    # bond labels are same as virtual_labels
    bond_tensors::Vector{<:AbstractVector{T}}
    max_index::LT

    # optimized contraction codes
    code_statetensor::Ein
    code_inner_product::Ein

    Dmax::Int
    ϵ::Float64
end
function VidalPEPS{NF}(vertex_labels::AbstractVector{<:AbstractVector{LT}}, vertex_tensors::Vector{<:AbstractArray{T}},
        virtual_labels::AbstractVector{LT}, bond_tensors::Vector{<:AbstractVector}, Dmax::Int, ϵ::Real,
        optimizer::CodeOptimizer, simplifier::CodeSimplifier) where {LT,T,NF}
    physical_labels = [vl[findall(∉(virtual_labels), vl)[]] for vl in vertex_labels]
    max_ind = max(maximum(physical_labels), maximum(virtual_labels))
    optcode_statetensor, optcode_inner_product = _optimized_code([vertex_labels..., [[l] for l in virtual_labels]...],
        physical_labels, virtual_labels, max_ind, NF, Dmax, optimizer, simplifier)
    VidalPEPS{NF}(physical_labels, virtual_labels, vertex_labels, vertex_tensors, bond_tensors, max_ind,
        optcode_statetensor, optcode_inner_product, Dmax, ϵ)
end

function VidalPEPS{NF}(physical_labels::Vector{LT},
    virtual_labels::Vector{LT},

    vertex_labels::Vector{Vector{LT}},
    vertex_tensors::Vector{<:AbstractArray{T}},
    bond_tensors::Vector{<:AbstractVector{T}},
    max_index::LT,

    code_statetensor::Ein,
    code_inner_product::Ein,

    Dmax::Int,
    ϵ::Float64
    ) where {LT,T,NF,Ein}
    VidalPEPS{T,NF,LT,Ein}(physical_labels, virtual_labels, vertex_labels, vertex_tensors, bond_tensors, max_index,
        code_statetensor, code_inner_product, Dmax, ϵ)
end

alllabels(peps::VidalPEPS) = [peps.vertex_labels..., [[l] for l in peps.virtual_labels]...]
alltensors(peps::VidalPEPS) = [peps.vertex_tensors..., peps.bond_tensors...]
# find bond tensor by virtual label
findbondtensor(peps::VidalPEPS, b) = peps.bond_tensors[findall(==(b), peps.virtual_labels)[]]
Base.copy(peps::VidalPEPS{T,NF}) where {T,NF} = VidalPEPS{NF}(copy(peps.physical_labels), copy(peps.virtual_labels),
    copy(peps.vertex_labels), copy(peps.virtual_tensors), copy(peps.bond_tensors), peps.max_index,
    peps.code_statetensor, peps.code_inner_product, peps.Dmax, peps.ϵ)
function Base.conj(peps::PEPS)
    replace_tensors(peps, conj.(alltensors(peps)))
end

# label APIs
getvlabel(peps::PEPS, i::Int) = peps.vertex_labels[i]  # vertex tensor labels
getphysicallabel(peps::PEPS, i::Int) = peps.physical_labels[i]  # physical label
newlabel(peps::PEPS, offset) = peps.max_index + offset  # create a new label
function virtualbonds(peps::PEPS)  # list all virtual bonds (a bond is a 2-tuple)
    bs = Tuple{Int,Int}[]
    for b in peps.virtual_labels
        i,j = findall(l->b∈l, peps.vertex_labels)
        push!(bs, (i, j))
    end
    return bs
end
nsite(peps::PEPS) = length(peps.physical_labels)
nflavor(::PEPS{T,NF}) where NF = NF
Dmax(peps::PEPS) = peps.Dmax
tensortype(peps::PEPS) = eltype(peps.vertex_tensors)
array_match_type(::Type{<:DiffTensor}, x::Array) = DiffTensor(x)
array_match_type(::Type{<:Array}, x::Array) = x

# all variables by flattening the tensors
variables(peps::PEPS) = vcat(vec.(alltensors(peps))...)
# load all variables to tensors
function load_variables!(peps::PEPS, variables)
    k = 0
    for t in alltensors(peps)
        t .= reshape(variables[k+1:k+length(t)], size(t))
        k += length(t)
    end
    return peps
end
function load_variables(peps::SimplePEPS, variables)  # for AD
    ats = alltensors(peps)
    ks = cumsum(length.(ats))
    tensors = map(1:length(ats)) do i
        t = ats[i]
        reshape(variables[(i>1 ? ks[i-1] : 0)+1:ks[i]], size(t))
    end
    return replace_tensors(peps, tensors)
end

function replace_tensors(peps::SimplePEPS{T,NF}, tensors) where {T,NF}
    SimplePEPS{NF}(peps.physical_labels, peps.virtual_labels,
        peps.vertex_labels, tensors, peps.max_index,
        peps.code_statetensor, peps.code_inner_product, peps.Dmax, peps.ϵ
    )
end
function replace_tensors(peps::VidalPEPS{T,NF}, tensors) where {T,NF}
    nv = nsite(peps)
    VidalPEPS{NF}(peps.physical_labels, peps.virtual_labels,
        peps.vertex_labels, tensors[1:nv], typeof(tensors[nv+1])[tensors[nv+1:end]...], peps.max_index,
        peps.code_statetensor, peps.code_inner_product, peps.Dmax, peps.ϵ
    )
end

function Base.show(io::IO, ::MIME"text/plain", peps::PEPS)
    println(io, typeof(peps), " (Dmax = $(peps.Dmax), ϵ = $(peps.ϵ))")
    println(io, " # of spins = $(nsite(peps))")
    print(io, " # of virtual degrees = $(length(peps.virtual_labels))")
end
Base.show(io::IO, peps::PEPS) = show(io, MIME"text/plain"(), peps)

# random and zero PEPSs
zero_vidalpeps(::Type{T}, g::SimpleGraph, D::Int; nflavor::Int=2, optimizer=GreedyMethod(), simplifier=MergeGreedy(), Dmax=D, ϵ=1e-12) where T = _peps_zero_state(Val(:Vidal), T, g, D, nflavor, Dmax, ϵ, optimizer, simplifier)
zero_simplepeps(::Type{T}, g::SimpleGraph, D::Int; nflavor::Int=2, optimizer=GreedyMethod(), simplifier=MergeGreedy(), Dmax=D, ϵ=1e-12) where T = _peps_zero_state(Val(:Simple), T, g, D, nflavor, Dmax, ϵ, optimizer, simplifier)
function _peps_zero_state(::Val{TYPE}, ::Type{T}, g::SimpleGraph, D::Int, nflavor::Int, Dmax::Int, ϵ::Real,
        optimizer::CodeOptimizer, simplifier::CodeSimplifier) where {TYPE, T}
    virtual_labels = collect(nv(g)+1:nv(g)+ne(g))
    vertex_labels = Vector{Int}[]
    vertex_tensors = Array{T}[]
    edge_map = Dict(zip(edges(g), virtual_labels))
    for i=1:nv(g)
        push!(vertex_labels, [i,[get(edge_map, SimpleEdge(i,nb), get(edge_map,SimpleEdge(nb,i),0)) for nb in neighbors(g, i)]...])
        t = zeros(T, nflavor, fill(D, degree(g, i))...)
        t[1] = 1
        push!(vertex_tensors, t)
    end
    if any(vl->any(iszero, vl), vertex_labels)
        error("incorrect input labels1")
    end
    if TYPE === :Vidal
        bond_tensors = [ones(T, D) for _=1:ne(g)]
        return VidalPEPS{nflavor}(vertex_labels, vertex_tensors, virtual_labels, bond_tensors, Dmax, ϵ, optimizer, simplifier)
    else
        return SimplePEPS{nflavor}(vertex_labels, vertex_tensors, virtual_labels, Dmax, ϵ, optimizer, simplifier)
    end
end

function rand_vidalpeps(::Type{T}, g::SimpleGraph, D::Int; nflavor::Int=2, optimizer=GreedyMethod(), simplifier=MergeGreedy(), Dmax::Int=D, ϵ=1e-12) where T
    randn!(zero_vidalpeps(T, g, D; nflavor, Dmax, ϵ, optimizer, simplifier))
end
function rand_simplepeps(::Type{T}, g::SimpleGraph, D::Int; nflavor::Int=2, optimizer=GreedyMethod(), simplifier=MergeGreedy(), Dmax::Int=D, ϵ=1e-12) where T
    randn!(zero_simplepeps(T, g, D; nflavor, Dmax, ϵ, optimizer, simplifier))
end
function Random.randn!(peps::PEPS)
    for t in alltensors(peps)
        randn!(t)
    end
    return peps
end

# multiply a peps with a constant
function LinearAlgebra.rmul!(peps::PEPS, c::Number)
    peps.vertex_tensors .*= c^(1/nsite(peps))
    return peps
end
Base.:*(c, peps::PEPS) = peps * c
function Base.:*(peps::PEPS, c::Number)   # to support AD? maybe not needed
    tensors = peps.vertex_tensors .* c^(1/nsite(peps))
    return replace_tensors(peps, tensors)
end
function Base.:*(peps::PEPS, c::AbstractArray{T,0}) where T   # to support AD
    tensors = [t * c .^ (1/nsite(peps)) for t in peps.vertex_tensors]
    return replace_tensors(peps, tensors)
end

# norm of a peps
#
# ●----●----●----●   ← ⟨peps|
# ┆    ┆    ┆    ┆
# ●----●----●----●   ← |peps⟩
LinearAlgebra.norm(peps::PEPS) = sqrt.(abs.(inner_product(peps, peps)))  # inner product returns a rank 0 array
function LinearAlgebra.normalize!(peps::PEPS)
    nm = sqrt.(abs.(inner_product(peps, peps)))
    return rmul!(peps, 1/nm[])
end

# contractor, the not cached version.
function direct_contract(code::EinCode, tensors)
    size_dict = OMEinsum.get_size_dict(OMEinsum.getixs(code), tensors)
    optcode = optimize_code(code, size_dict, GreedyMethod())
    optcode(tensors...)
end

# contract the peps and obtain the state vector
#
# ┆    ┆    ┆    ┆
# ●----●----●----●   ← |peps⟩
Base.vec(peps::PEPS) = vec(statetensor(peps))
function statetensor(peps::PEPS)
    peps.code_statetensor(alltensors(peps)...)
end

# apply a single site operator
#      ┆
#      ■  ← (operator)
# ┆    ┆    ┆    ┆
# ●----●----●----●   ← |peps⟩
function apply_onsite!(peps::PEPS{T,LT}, i, mat::AbstractMatrix) where {T,LT}
    @assert size(mat, 1) == size(mat, 2)
    ti = peps.vertex_tensors[i]
    old = getvlabel(peps, i)
    mlabel = [newlabel(peps, 1), getphysicallabel(peps, i)]
    peps.vertex_tensors[i] = EinCode([old, mlabel], replace(old, mlabel[2]=>mlabel[1]))(ti, mat)
    return peps
end

function apply_onsite(peps::PEPS{T,LT}, i, mat::AbstractMatrix) where {T,LT}  # to support AD
    @assert size(mat, 1) == size(mat, 2)
    old = getvlabel(peps, i)
    mlabel = [newlabel(peps, 1), getphysicallabel(peps, i)]
    tensors = map(1:length(peps.vertex_tensors)) do j
        tj = peps.vertex_tensors[j]
        j == i ? EinCode([old, mlabel], replace(old, mlabel[2]=>mlabel[1]))(tj, mat) : tj
    end
    return replace_tensors(peps, tensors)
end
# apply a two site operator
#      ┆    ┆
#      ■----■  ← (operator)
# ┆    ┆    ┆    ┆
# ●----●----●----●    ← |peps⟩
function apply_onbond!(peps::PEPS, i, j, mat::AbstractArray{T,4}) where T
    ti, tj = peps.vertex_tensors[i], peps.vertex_tensors[j]
    li, lj = getvlabel(peps, i), getvlabel(peps, j)
    shared_label = li ∩ lj; @assert length(shared_label) == 1
    only_left, only_right = setdiff(li, lj), setdiff(lj, li)
    lij = only_left ∪ only_right

    if peps isa VidalPEPS
        # multiple S tensors
        for b in li
            b == getphysicallabel(peps, i) && continue
            ti = _apply_vec(ti, li, sqrt.(findbondtensor(peps, b)), b)
        end
        for b in lj
            b == getphysicallabel(peps, j) && continue
            tj = _apply_vec(tj, lj, sqrt.(findbondtensor(peps, b)), b)
        end
    end

    # contract
    tij = EinCode([li, lj], lij)(ti, tj)
    lijkl = [newlabel(peps, 1), newlabel(peps, 2), getphysicallabel(peps, i), getphysicallabel(peps, j)]
    lkl = replace(lij, lijkl[3]=>lijkl[1], lijkl[4]=>lijkl[2])
    tkl = EinCode([lij, lijkl], lkl)(tij, mat)

    # SVD and truncate
    sl = [size(ti, findall(==(l), li)[]) for l in only_left]
    sr = [size(tj, findall(==(l), lj)[]) for l in only_right]
    U, S, Vt = svd(reshape(tkl, prod(sl), prod(sr)))
    Vt = conj.(Vt)
    D0 = findfirst(<(peps.ϵ), S)
    D = D0 === nothing ? min(peps.Dmax, size(U,2)) : min(D0-1, peps.Dmax)
    if D < size(U, 2)
        @info "truncation error is $(sum(S[D+1:end]))"
        S = S[1:D]
        U = U[:,1:D]
        Vt = Vt[:,1:D]
    end
    # reshape back
    ti_ = EinCode([[only_left..., shared_label[]]], li)(reshape(U, (sl..., size(U,2))))
    tj_ = EinCode([[only_right..., shared_label[]]], lj)(reshape(Vt, (sr..., size(Vt,2))))

    if peps isa VidalPEPS
        # devide S tensors
        for b in only_left
            b == getphysicallabel(peps, i) && continue
            ti_ = _apply_vec(ti_, li, safe_inv.(sqrt.(findbondtensor(peps, b))), b)  # assume S being positive
        end
        for b in only_right
            b == getphysicallabel(peps, j) && continue
            tj_ = _apply_vec(tj_, lj, safe_inv.(sqrt.(findbondtensor(peps, b))), b)
        end
        peps.bond_tensors[findall(==(shared_label[]), peps.virtual_labels)[]] = S
    else
        b = shared_label[]
        sqrtS = sqrt.(S)
        ti_ = _apply_vec(ti_, li, sqrtS, b)  # assume S being positive
        tj_ = _apply_vec(tj_, lj, sqrtS, b)
    end

    # update tensors
    peps.vertex_tensors[i] = ti_
    peps.vertex_tensors[j] = tj_
    return peps
end

function safe_inv(y::T) where T
    ϵ = 1e-10
    if abs(y) < ϵ
        return inv(real(y) >= 0 ? T(ϵ) : T(-ϵ))
    else
        return inv(y)
    end
end

function _apply_vec(t, l, v, b)
    return EinCode([l, [b]], l)(t, v)
end

#### Circuit interfaces ####
# NOTE: we should have a register type.
Yao.nqubits(peps::PEPS) = nsite(peps)
Yao.nactive(peps::PEPS) = nsite(peps)
Yao.statevec(peps::PEPS) = vec(peps)
for (APPLY, APPLY_ONSITE, MUL) in [(:unsafe_apply!, :apply_onsite!, :mul!),
        (:apply, :apply_onsite, :*)]
    @eval function YaoBlocks.$APPLY(peps::PEPS{T}, block::PutBlock{2,1}) where {T}
        $APPLY_ONSITE(peps, block.locs[1], array_match_type(tensortype(peps), Matrix{T}(block.content)))
    end
    @eval function YaoBlocks.$APPLY(peps::PEPS{T}, block::KronBlock{D,M,BT}) where {D,M,T,BT<:NTuple{M,AbstractBlock}}
        for (loc, g) in zip(block.locs, subblocks(block))
            peps = $APPLY_ONSITE(peps, loc[1], array_match_type(tensortype(peps), Matrix{T}(g)))
        end
        return peps
    end
    @eval function YaoBlocks.$APPLY(peps::PEPS{T}, block::ControlBlock{BT,1,1}) where {T,BT}
        # forward to PutBlock.
        YaoBlocks.$APPLY(peps, put(nqudits(block), (block.ctrl_locs[1], block.locs[1])=>control(2,1,2=>block.content)))
    end
    @eval function YaoBlocks.$APPLY(peps::PEPS{T}, block::Scale) where T
        $MUL(YaoBlocks.$APPLY(peps, content(block)), (1.0+0im)*Yao.factor(block))
    end
end
# no non-inplace version defined.
YaoBlocks.unsafe_apply!(peps::PEPS, block::PutBlock{2,2}) = invoke(YaoBlocks.unsafe_apply!, Tuple{PEPS,PutBlock{D,2} where D}, peps, block)
function YaoBlocks.unsafe_apply!(peps::PEPS{T}, block::PutBlock{D,2}) where {T,D}
    apply_onbond!(peps, block.locs..., reshape(Matrix{T}(block.content), 2, 2, 2, 2))
end

# compute the expectation value of a Hamiltonian
#
# ●----●----●----●   ← ⟨peps|
# ┆    ┆    ┆    ┆
# ■    ■    ■    ■   ← (product operator)
# ┆    ┆    ┆    ┆
# ●----●----●----●   ← |peps⟩
function Yao.expect(operator::Add, pa::PEPS, pb::PEPS)
    return sum(Yao.subblocks(operator)) do term
        expect(term, pa, pb)
    end
end
function Yao.expect(operator::AbstractBlock, pa::PEPS, pb::PEPS)
    opb = apply(pb, operator)
    inner_product(pa, opb)
end


text:
"
ReplacementOperator represents the method that replaces the previous generation
by a new one originated after crossover and mutation.

# Fields
- `needsComparison::Bool = true`: tells wether the user wants the system to compare
    between the individuals in order to select the ones that will replace the old
    population.
- `_eliteSize::UInt16`: number of fittest individuals to be maintained if there
    is elitism.
- `_method::Function`: method used for replacement.
- `_varArgs::Array{Any}`: arguments necessary for the replacement method.
"
struct ReplacementOperator
    _needsComparison::Bool
    _eliteSize::UInt16
    _method::Function
    _varArgs::Array{Any}
end # struct

"""
    needsComparison(replaceOp::ReplacementOperator)::Bool

Indicates wether the replacement method needs to compare between individuals
or not.
"""
needsComparison(replaceOp::ReplacementOperator)::Bool = replaceOp._needsComparison
# function

"""
    isElitist(replaceOp::ReplacementOperator)::Bool

Indicates wether the replacement method is elitist or not.
"""
isElitist(replaceOp::ReplacementOperator)::Bool = replaceOp._eliteSize > 0
# function

"""
    getEliteSize(replaceOp::ReplacementOperator)::UInt16

Obtains the size of the elite set to be maintained throughtout the generations.
"""
getEliteSize(replaceOp::ReplacementOperator)::UInt16 = replaceOp._eliteSize
# function

"""
    getDefaultReplacementOp()::Tuple{Function, Array}

Returns the default replacement operator, which is [`replaceAllPopulation`](@ref).
"""
function getDefaultReplacementOp()::Tuple{Function, Array}
    return replaceAllPopulation, []
end # function

"""
    selectFromSet(Ra::Array{Integer}, nLeft::Integer, population::Array{Individual},
                  compare::Function, rng::Random.AbstractRNG)

Selects `nLeft` individuals from the set `Ra` using a comparison method.
"""
function selectFromSet(Ra::Array, nLeft::Integer, population::Array{Individual},
                        compare::Function, rng::Random.AbstractRNG)

    subPop = population[Ra]

    for i=1:nLeft
        index = tournamentSelector(subPop[1:end-i+1], compare, rng, 4)
        subPop[index], subPop[end-i+1] = subPop[end-i+1], subPop[index]
        Ra[index], Ra[end-i+1] = Ra[end-i+1], Ra[index]
    end

    return Ra[end-nLeft+1:end]
end # function

"""
    fillNewPopulation(populationSelectedIndexes::Array{Integer},
                      offspringSelectedIndexes::Array{Integer},
                      population::Array{Individual},
                      offspring::Array{Individual},
                      compare::Function,
                      nLeft::Integer, rng::Random.AbstractRNG)

Fills a set of individuals with `nLeft` individuals using a comparison function.
"""
function fillNewPopulation(populationSelectedIndexes::Array,
                           offspringSelectedIndexes::Array,
                           population::Array{Individual},
                           offspring::Array{Individual},
                           compare::Function,
                           nLeft::Integer, rng::Random.AbstractRNG)

    lenPop = length(population)
    lenOff = length(offspring)

    left = Array{Integer}(undef, nLeft)

    belongs1(element) = findnext(x->x==element, offspringSelectedIndexes, 1) != nothing
    notSelectedOff = filter(!belongs1, collect(1:lenOff))

    newSelectedOff = []
    newSelectedPop = []
    if length(notSelectedOff) > nLeft
        newSelectedOff = selectFromSet(notSelectedOff, nLeft, offspring, compare, rng)
    else
        newSelectedOff = notSelectedOff
        if (nLeft - length(notSelectedOff)) != 0
            belongs2(element) = findnext(x->x==element, populationSelectedIndexes, 1) != nothing
            notSelectedPopulation = filter(!belongs2, collect(1:lenPop))
            newSelectedPop = selectFromSet(notSelectedPopulation, nLeft - length(notSelectedOff), population, compare, rng)
        end
    end

    return vcat(newSelectedPop, populationSelectedIndexes), vcat(newSelectedOff, offspringSelectedIndexes)
end # function

"""
    replacePopulation_(replaceOp::ReplacementOperator,
                       population::Array{Individual},
                       offspring::Array{Individual},
                       evaluator::Evaluator, rng::Random.AbstractRNG)

Performs the replacement of the population by the given replacement method using
the comparison function of the `Evaluator` to compare between individuals if necessary.
"""
function replacePopulation_(replaceOp::ReplacementOperator,
                            population::Array{Individual},
                            offspring::Array{Individual},
                            evaluator::Evaluator, rng::Random.AbstractRNG)

    if needsComparison(replaceOp)
        compare = getCompareFunction(evaluator)
        populationIndexes, offspringIndexes = replaceOp._method(population, offspring,
                                                                compare, replaceOp._varArgs...)
    else
        populationIndexes, offspringIndexes = replaceOp._method(population, offspring,
                                                                replaceOp._varArgs...)
    end

    nSelectedPop = length(populationIndexes)
    nSelectedOff = length(offspringIndexes)
    popSize = length(population)

    nLeft = popSize - (nSelectedPop + nSelectedOff)

    # Selected offspring is bigger or equal in size than popSize, which means we
    # only choose individuals among the user's selected offspring
    if nLeft < 0 && nSelectedOff >= popSize
        if !isElitist(replaceOp)
            offspringIndexes = selectFromSet(offspringIndexes, popSize, offspring,
                                             getCompareFunction(evaluator), rng)
            populationIndexes = []
        else
            offspringIndexes = selectFromSet(offspringIndexes, popSize-getEliteSize(replaceOp),
                                             offspring, getCompareFunction(evaluator), rng)
            populationIndexes = getBestIndividual(population,
                                                  getCompareFunction(evaluator),
                                                  getEliteSize(replaceOp))
        end

    # Selected offspring is less in size than popSize, but the sum of selected
    # individuals is bigger than popSize, which means we must choose the remaining
    # ones among the user's selected population plus all the user's selected offspring.
    elseif nLeft < 0
        if !isElitist(replaceOp)
            populationIndexes = selectFromSet(populationIndexes, popSize - nSelectedOff,
                                              population, getCompareFunction(evaluator), rng)
        else
            eliteSize = getEliteSize(replaceOp)

            if nSelectedOff + eliteSize > popSize
                function belongs1(element)
                    findnext(x->x==element, getWorstIndividual(offspring[offspringIndexes],
                                                               getCompareFunction(evaluator),
                                                               nSelectedOff + eliteSize - popSize), 1) != nothing
                end
                offspringIndexes = filter(!belongs1, offspringIndexes)

                populationIndexes = vcat(populationIndexes,
                                         getBestIndividual(population,
                                                           getCompareFunction(evaluator),
                                                           eliteSize))
            elseif nSelectedOff + eliteSize < popSize
                populationIndexes = getBestIndividual(population,
                                                      getCompareFunction(evaluator),
                                                      eliteSize)

                populationIndexes = vcat(populationIndexes,
                                         selectFromSet(populationIndexes,
                                                       popSize - (nSelectedOff + eliteSize),
                                                       population, getCompareFunction(evaluator),
                                                       rng))
            else
                populationIndexes = vcat(populationIndexes,
                                         getBestIndividual(population,
                                                           getCompareFunction(evaluator),
                                                           eliteSize))
            end
        end

    # The sum of selected individuals is less in size than popSize, we choose all
    # the user's selected individuals and we select new ones
    # through the remaining not selected offspring. If there is no more selected
    # offspring, we select individuals from the ones not selected of the population.
    elseif nLeft > 0
        if !isElitist(replaceOp)
            populationIndexes, offspringIndexes = fillNewPopulation(populationIndexes,
                                                                    offspringIndexes,
                                                                    population,
                                                                    offspring,
                                                                    getCompareFunction(evaluator), nLeft, rng)
        else
            eliteSize = getEliteSize(replaceOp)

            if nSelectedPop + nSelectedOff + eliteSize > popSize
                if !isempty(populationIndexes)
                    function belongs2(element)
                        findnext(x->x==element,
                                 getWorstIndividual(population[populationIndexes],
                                                    getCompareFunction(evaluator),
                                                    nSelectedPop+nSelectedOff+eliteSize-popSize),
                                 1) != nothing
                    end
                    populationIndexes = filter(!belongs2, populationIndexes)
                else
                    function belongs3(element)
                        findnext(x->x==element,
                                 getWorstIndividual(offspring[offspringIndexes],
                                                    getCompareFunction(evaluator),
                                                    nSelectedPop+nSelectedOff+eliteSize-popSize),
                                 1) != nothing
                    end
                    offspringIndexes = filter(!belongs3, offspringIndexes)
                end

                populationIndexes = vcat(populationIndexes,
                                         getBestIndividual(population,
                                                           getCompareFunction(evaluator),
                                                           eliteSize))
            elseif nSelectedPop + nSelectedOff + eliteSize < popSize
                populationIndexes = vcat(populationIndexes,
                                         getBestIndividual(population,
                                                           getCompareFunction(evaluator),
                                                           eliteSize))

                populationIndexes, offspringIndexes = fillNewPopulation(populationIndexes,
                                                                        offspringIndexes,
                                                                        population,
                                                                        offspring,
                                                                        getCompareFunction(evaluator),
                                                                        popSize-(nSelectedPop+nSelectedOff+eliteSize),
                                                                        rng)
            else
                populationIndexes = vcat(populationIndexes,
                                         getBestIndividual(population,
                                                           getCompareFunction(evaluator),
                                                           eliteSize))
            end
        end
    elseif nLeft == 0 && isElitist(replaceOp)
        eliteSize = getEliteSize(replaceOp)

        if !isempty(populationIndexes)
            function belongs4(element)
                findnext(x->x==element,
                         getWorstIndividual(population[populationIndexes],
                                            getCompareFunction(evaluator),
                                            eliteSize),
                         1) != nothing
            end
            populationIndexes = filter(!belongs4, populationIndexes)
        else
            function belongs5(element)
                findnext(x->x==element,
                         getWorstIndividual(offspring[offspringIndexes],
                                            getCompareFunction(evaluator),
                                            eliteSize),
                         1) != nothing
            end
            offspringIndexes = filter(!belongs5, offspringIndexes)
        end

        populationIndexes = vcat(populationIndexes,
                                 getBestIndividual(population,
                                                   getCompareFunction(evaluator),
                                                   eliteSize))
    end

    newPopulation = vcat(population[populationIndexes], offspring[offspringIndexes])

    return newPopulation
end # function

"""
    replaceAllPopulation(population::Array{Individual}, offspring::Array{Individual})

Performs a replacement of all the current population by the offspring.
"""
function replaceAllPopulation(population::Array{Individual}, offspring::Array{Individual})
    return [], collect(1:length(offspring))
end # function

"""
    replaceWorstIndividuals(population::Array{Individual}, offspring::Array{Individual}
                                 compareFunction::Function)

documentation
"""
function replaceWorstIndividuals(population::Array{Individual}, offspring::Array{Individual},
                                 compareFunction::Function)

    popSize = length(population)
    offSize = length(offspring)
    local populationIndexes
    local offspringIndexes

    if offSize >= popSize
        populationIndexes, offspringIndexes = [], collect(1:length(offspring))
    else

        function belongs(element)
            findnext(x->x==element,
                     getWorstIndividual(population,
                                        compareFunction,
                                        offSize),
                     1) != nothing
        end
        populationIndexes = filter(!belongs, collect(1:popSize))
        offspringIndexes = collect(1:offSize)
    end

    return populationIndexes, offspringIndexes
end # function

ext: jl

max_stars_repo_name: SergioGmezM/GeneticJulia

max_issues_repo_name: SergioGmezM/GeneticJulia

max_forks_repo_name: SergioGmezM/GeneticJulia


text:
module RocketDefaultIfEmptyOperatorTest

using Test
using Rocket

include("../test_helpers.jl")

@testset "operator: default_if_empty()" begin

    println("Testing: operator default_if_empty()")

    run_proxyshowcheck("DefaultIfEmpty", default_if_empty(0))
    run_proxyshowcheck("DefaultIfEmpty", default_if_empty(() -> 0))

    run_testset([
        (
            source = from(1:5) |> default_if_empty(0),
            values = @ts([ 1:5, c ])
        ),
        (
            source = completed(Int) |> default_if_empty(0),
            values = @ts([ 0, c ])
        ),
        (
            source      = completed(Int) |> default_if_empty("string"),
            values      = @ts([ "string", c ]),
            source_type = Union{Int, String}
        ),
        (
            source      = completed(Int) |> default_if_empty(0) |> default_if_empty("string"),
            values      = @ts([ 0, c ]),
            source_type = Union{Int, String}
        ),
        (
            source      = completed(Int) |> default_if_empty("string") |> default_if_empty(0) ,
            values      = @ts([ "string", c ]),
            source_type = Union{Int, String}
        ),
        (
            source = from(1:5) |> default_if_empty(() -> 0),
            values = @ts([ 1:5, c ]),
            source_type = Int,
        ),
        (
            source = completed(Int) |> default_if_empty(() -> 0),
            values = @ts([ 0, c ]),
            source_type = Int,
        ),
        (
            source      = completed(Int) |> default_if_empty(() -> 0) |> default_if_empty(() -> 1),
            values      = @ts([ 0, c ]),
            source_type = Int
        ),
        (
            source      = completed(Int) |> default_if_empty(() -> 1) |> default_if_empty(() -> 0) ,
            values      = @ts([ 1, c ]),
            source_type = Int
        ),
        (
            source      = completed(Int) |> default_if_empty(0) |> default_if_empty(() -> 1),
            values      = @ts([ 0, c ]),
            source_type = Int
        ),
        (
            source      = completed(Int) |> default_if_empty(1) |> default_if_empty(() -> 0) ,
            values      = @ts([ 1, c ]),
            source_type = Int
        ),
        (
            source      = completed(Int) |> default_if_empty(() -> 0) |> default_if_empty(1),
            values      = @ts([ 0, c ]),
            source_type = Int
        ),
        (
            source      = completed(Int) |> default_if_empty(() -> 1) |> default_if_empty(0) ,
            values      = @ts([ 1, c ]),
            source_type = Int
        ),
        (
            source      = completed(Int) |> default_if_empty(() -> error("42")) |> safe(),
            values      = @ts(e),
            source_type = Int
        ),
    ])

end

end


text:
@doc raw"""
    device_range(psi_container::PSIContainer,
                 range_data::Vector{DeviceRange},
                 cons_name::Symbol,
                 var_name::Symbol)

Constructs min/max range constraint from device variable.

# Constraints
If min and max within an epsilon width:

``` variable[name, t] == limits.max ```

Otherwise:

``` limits.min <= variable[name, t] <= limits.max ```

where limits in range_data.

# LaTeX

`` x = limits^{max}, \text{ for } |limits^{max} - limits^{min}| < \varepsilon ``

`` limits^{min} \leq x \leq limits^{max}, \text{ otherwise } ``

# Arguments
* psi_container::PSIContainer : the psi_container model built in PowerSimulations
* range_data::Vector{DeviceRange} : contains names and vector of min/max
* cons_name::Symbol : name of the constraint
* var_name::Symbol : the name of the continuous variable
"""
function device_range(
    psi_container::PSIContainer,
    range_data::Vector{DeviceRange},
    cons_name::Symbol,
    var_name::Symbol,
)
    time_steps = model_time_steps(psi_container)
    variable = get_variable(psi_container, var_name)
    ub_name = middle_rename(cons_name, PSI_NAME_DELIMITER, "ub")
    lb_name = middle_rename(cons_name, PSI_NAME_DELIMITER, "lb")
    names = (d.name for d in range_data)
    con_ub = add_cons_container!(psi_container, ub_name, names, time_steps)
    con_lb = add_cons_container!(psi_container, lb_name, names, time_steps)

    for data in range_data, t in time_steps
        expression_ub = JuMP.AffExpr(0.0, variable[data.name, t] => 1.0)
        for val in data.additional_terms_ub
            JuMP.add_to_expression!(
                expression_ub,
                get_variable(psi_container, val)[data.name, t],
            )
        end
        expression_lb = JuMP.AffExpr(0.0, variable[data.name, t] => 1.0)
        for val in data.additional_terms_lb
            JuMP.add_to_expression!(
                expression_lb,
                get_variable(psi_container, val)[data.name, t],
                -1.0,
            )
        end
        con_ub[data.name, t] =
            JuMP.@constraint(psi_container.JuMPmodel, expression_ub <= data.limits.max)
        con_lb[data.name, t] =
            JuMP.@constraint(psi_container.JuMPmodel, expression_lb >= data.limits.min)
    end

    return
end

@doc raw"""
    device_semicontinuousrange(psi_container::PSIContainer,
                                    range_data::Vector{DeviceRange},
                                    cons_name::Symbol,
                                    var_name::Symbol,
                                    binvar_name::Symbol)

Constructs min/max range constraint from device variable and on/off decision variable.

# Constraints
If device min = 0:

``` varcts[name, t] <= limits.max*varbin[name, t]) ```

``` varcts[name, t] >= 0.0 ```

Otherwise:

``` varcts[name, t] <= limits.max*varbin[name, t] ```

``` varcts[name, t] >= limits.min*varbin[name, t] ```

where limits in range_data.

# LaTeX

`` 0 \leq x^{cts} \leq limits^{max} x^{bin}, \text{ for } limits^{min} = 0 ``

`` limits^{min} x^{bin} \leq x^{cts} \leq limits^{max} x^{bin}, \text{ otherwise } ``

# Arguments
* psi_container::PSIContainer : the psi_container model built in PowerSimulations
* range_data::Vector{DeviceRange} : contains names and vector of min/max
* cons_name::Symbol : name of the constraint
* var_name::Symbol : the name of the continuous variable
* binvar_name::Symbol : the name of the binary variable
"""
function device_semicontinuousrange(
    psi_container::PSIContainer,
    range_data::Vector{DeviceRange},
    cons_name::Symbol,
    var_name::Symbol,
    binvar_name::Symbol,
)
    time_steps = model_time_steps(psi_container)
    varcts = get_variable(psi_container, var_name)
    varbin = get_variable(psi_container, binvar_name)
    ub_name = middle_rename(cons_name, PSI_NAME_DELIMITER, "ub")
    lb_name = middle_rename(cons_name, PSI_NAME_DELIMITER, "lb")
    names = (d.name for d in range_data)
    #MOI has a semicontinous set, but after some tests is not clear most MILP solvers support it.
    #In the future this can be updated
    con_ub = add_cons_container!(psi_container, ub_name, names, time_steps)
    con_lb = add_cons_container!(psi_container, lb_name, names, time_steps)

    for data in range_data, t in time_steps
        if JuMP.has_lower_bound(varcts[data.name, t])
            JuMP.set_lower_bound(varcts[data.name, t], 0.0)
        end
        expression_ub = JuMP.AffExpr(0.0, varcts[data.name, t] => 1.0)
        for val in data.additional_terms_ub
            JuMP.add_to_expression!(
                expression_ub,
                get_variable(psi_container, val)[data.name, t],
            )
        end
        expression_lb = JuMP.AffExpr(0.0, varcts[data.name, t] => 1.0)
        for val in data.additional_terms_lb
            JuMP.add_to_expression!(
                expression_lb,
                get_variable(psi_container, val)[data.name, t],
                -1.0,
            )
        end
        con_ub[data.name, t] = JuMP.@constraint(
            psi_container.JuMPmodel,
            expression_ub <= data.limits.max * varbin[data.name, t]
        )
        con_lb[data.name, t] = JuMP.@constraint(
            psi_container.JuMPmodel,
            expression_lb >= data.limits.min * varbin[data.name, t]
        )
    end

    return
end
#This function looks suspicious and repetitive. Needs verification
@doc raw"""
    reserve_device_semicontinuousrange(psi_container::PSIContainer,
                                    range_data::Vector{DeviceRange},
                                    cons_name::Symbol,
                                    var_name::Symbol,
                                    binvar_name::Symbol)

Constructs min/max range constraint from device variable and on/off decision variable.

# Constraints
If device min = 0:

``` varcts[name, t] <= limits.max * (1 - varbin[name, t]) ```

``` varcts[name, t] >= 0.0 ```

Otherwise:

``` varcts[name, t] <= limits.max * (1 - varbin[name, t]) ```

``` varcts[name, t] >= limits.min * (1 - varbin[name, t]) ```

where limits in range_data.

# LaTeX

`` 0 \leq x^{cts} \leq limits^{max} (1 - x^{bin} ), \text{ for } limits^{min} = 0 ``

`` limits^{min} (1 - x^{bin} ) \leq x^{cts} \leq limits^{max} (1 - x^{bin} ), \text{ otherwise } ``

# Arguments
* psi_container::PSIContainer : the psi_container model built in PowerSimulations
* range_data::Vector{DeviceRange} : contains names and vector of min/max
* cons_name::Symbol : name of the constraint
* var_name::Symbol : the name of the continuous variable
* binvar_name::Symbol : the name of the binary variable
"""
function reserve_device_semicontinuousrange(
    psi_container::PSIContainer,
    range_data::Vector{DeviceRange},
    cons_name::Symbol,
    var_name::Symbol,
    binvar_name::Symbol,
)

    time_steps = model_time_steps(psi_container)
    varcts = get_variable(psi_container, var_name)
    varbin = get_variable(psi_container, binvar_name)

    ub_name = middle_rename(cons_name, PSI_NAME_DELIMITER, "ub")
    lb_name = middle_rename(cons_name, PSI_NAME_DELIMITER, "lb")
    names = (d.name for d in range_data)
    #MOI has a semicontinous set, but after some tests is not clear most MILP solvers support it.
    #In the future this can be updated
    con_ub = add_cons_container!(psi_container, ub_name, names, time_steps)
    con_lb = add_cons_container!(psi_container, lb_name, names, time_steps)

    for data in range_data, t in time_steps
        if JuMP.has_lower_bound(varcts[data.name, t])
            JuMP.set_lower_bound(varcts[data.name, t], 0.0)
        end
        expression_ub = JuMP.AffExpr(0.0, varcts[data.name, t] => 1.0)
        for val in data.additional_terms_ub
            JuMP.add_to_expression!(
                expression_ub,
                get_variable(psi_container, val)[data.name, t],
            )
        end
        expression_lb = JuMP.AffExpr(0.0, varcts[data.name, t] => 1.0)
        for val in data.additional_terms_lb
            JuMP.add_to_expression!(
                expression_lb,
                get_variable(psi_container, val)[data.name, t],
                -1.0,
            )
        end
        con_ub[data.name, t] = JuMP.@constraint(
            psi_container.JuMPmodel,
            expression_ub <= data.limits.max * (1 - varbin[data.name, t])
        )
        con_lb[data.name, t] = JuMP.@constraint(
            psi_container.JuMPmodel,
            expression_lb >= data.limits.min * (1 - varbin[data.name, t])
        )
    end
    return
end


text:
# rmsd100
# root mean-squared difference between two arrays in %

export rmsd100

using Test: @test, @test_throws, @inferred
using LinearAlgebra: norm

"""
    rmsd = rmsd100(x, y ; mask)

Compute 100 * RMSD (root mean squared difference) between `x` and `y` within domain mask.

in
- `x` : array
- `y` : another array of same size

option:
- `mask::Array{Bool}` : domain over which to compute the RMSE; default `trues(size(x))`

out
- rmsd : rmsd of `x` vs `y` within `mask` in %
"""
function rmsd100(x::AbstractArray{<:Number}, y::AbstractArray{<:Number}
        ; mask::AbstractArray{Bool} = trues(size(x)))
    return 100 * norm(x[mask] - y[mask]) / sqrt(sum(mask))
end

"""
    rmsd100(:test)
self test
"""
function rmsd100(x::Symbol)
    x != :test && throw("Invalid symbol; use :test for testing functions.")
    a = [1 2; 4 5; 7 8]
    b = [3 2; 4 6; 8 7]
    c = [false true; true false; false false] #avg. 0
    d = [false false; false false; false true] #avg. 1
    e = [true true; true true; true false] # average 1,0,1,0,4 = sqrt(6/5)
    tests = cat(c, d, e, dims = 3)
    pred = [0, 100, 100*sqrt(6/5)]
    for i in 1:size(tests,1)
        @test rmsd100(a, b ; mask=tests[:,:,i]) == pred[i]
    end
    @test rmsd100([1, 0], [0, 1]) == 100 # no mask test

#=  not important to support or test.  user can pass mask=ig.mask if needed
    a = ones(300,400)
    b = ones(300,400)
    b[100:200,100:200] .= 0
    #so there's a blot of zeros from 100,100 to 200,200
    b[200:250,300:320] .= 3
    #and a blot of threes over here
    ig1 = image_geom(nx = 300,ny = 400,dx = 1,dy = 1)
    @test my_rmse(a,b,ig1) > 20
=#
    return true
end


text:
"""
    struct MATModel

Wrapper around the models loaded in dictionaries from the MATLAB representation.
"""
struct MATModel <: MetabolicModel
    mat::Dict{String,Any}
end

n_metabolites(m::MATModel)::Int = size(m.mat["S"], 1)
n_reactions(m::MATModel)::Int = size(m.mat["S"], 2)

"""
    reactions(m::MATModel)::Vector{String}

Extracts reaction names from `rxns` key in the MAT file.
"""
function reactions(m::MATModel)::Vector{String}
    if haskey(m.mat, "rxns")
        reshape(m.mat["rxns"], n_reactions(m))
    else
        "rxn" .* string.(1:n_reactions(m))
    end
end

"""
    _mat_has_squashed_coupling(mat)

Guesses whether C in the MAT file is stored in A=[S;C].
"""
_mat_has_squashed_coupling(mat) =
    haskey(mat, "A") && haskey(mat, "b") && length(mat["b"]) == size(mat["A"], 1)

"""
    metabolites(m::MATModel)::Vector{String}

Extracts metabolite names from `mets` key in the MAT file.
"""
function metabolites(m::MATModel)::Vector{String}
    nm = n_metabolites(m)
    if haskey(m.mat, "mets")
        reshape(m.mat["mets"], length(m.mat["mets"]))[begin:nm]
    else
        "met" .* string.(1:n_metabolites(m))
    end
end

"""
    stoichiometry(m::MATModel)

Extract the stoichiometry matrix, stored under key `S`.
"""
stoichiometry(m::MATModel) = sparse(m.mat["S"])

"""
    bounds(m::MATModel)

Extracts bounds from the MAT file, saved under `lb` and `ub`.
"""
bounds(m::MATModel) = (
    sparse(reshape(get(m.mat, "lb", fill(-Inf, n_reactions(m), 1)), n_reactions(m))),
    sparse(reshape(get(m.mat, "ub", fill(Inf, n_reactions(m), 1)), n_reactions(m))),
)

"""
    balance(m::MATModel)

Extracts balance from the MAT model, defaulting to zeroes if not present.
"""
function balance(m::MATModel)
    b = get(m.mat, "b", spzeros(n_metabolites(m), 1))
    if _mat_has_squashed_coupling(m.mat)
        b = b[1:n_metabolites(m), :]
    end
    sparse(reshape(b, n_metabolites(m)))
end

"""
    objective(m::MATModel)

Extracts the objective from the MAT model (defaults to zeroes).
"""
objective(m::MATModel) =
    sparse(reshape(get(m.mat, "c", zeros(n_reactions(m), 1)), n_reactions(m)))

"""
    coupling(m::MATModel)

Extract coupling matrix stored, in `C` key.
"""
coupling(m::MATModel) =
    _mat_has_squashed_coupling(m.mat) ? sparse(m.mat["A"][n_reactions(m)+1:end, :]) :
    sparse(get(m.mat, "C", zeros(0, n_reactions(m))))

"""
    coupling_bounds(m::MATModel)

Extracts the coupling constraints. Currently, there are several accepted ways to store these in MATLAB models; this takes the constraints from vectors `cl` and `cu`.
"""
function coupling_bounds(m::MATModel)
    nc = n_coupling_constraints(m)
    if _mat_has_squashed_coupling(m.mat)
        (
            sparse(fill(-Inf, nc)),
            sparse(reshape(m.mat["b"], length(m.mat["b"]))[n_reactions(m)+1:end]),
        )
    else
        (
            sparse(reshape(get(m.mat, "cl", fill(-Inf, nc, 1)), nc)),
            sparse(reshape(get(m.mat, "cu", fill(Inf, nc, 1)), nc)),
        )
    end
end

"""
    genes(m::MATModel)

Extracts the possible gene list from `genes` key.
"""
function genes(m::MATModel)
    x = get(m.mat, "genes", [])
    reshape(x, length(x))
end

"""
    reaction_gene_association(m::MATModel, rid::String)

Extracts the associations from `grRules` key, if present.
"""
function reaction_gene_association(m::MATModel, rid::String)
    if haskey(m.mat, "grRules")
        grr = m.mat["grRules"][findfirst(==(rid), reactions(m))]
        typeof(grr) == String ? _parse_grr(grr) : nothing
    else
        nothing
    end
end

"""
    metabolite_formula(m::MATModel, mid::String)

Extract metabolite formula from key `metFormula` or `metFormulas`.
"""
metabolite_formula(m::MATModel, mid::String) = _maybemap(
    x -> _parse_formula(x[findfirst(==(mid), metabolites(m))]),
    get(m.mat, "metFormula", get(m.mat, "metFormulas", nothing)),
)

"""
    metabolite_charge(m::MATModel, mid::String)

Extract metabolite charge from `metCharge` or `metCharges`.
"""
metabolite_charge(m::MATModel, mid::String) = _maybemap(
    x -> x[findfirst(==(mid), metabolites(m))],
    get(m.mat, "metCharge", get(m.mat, "metCharges", nothing)),
)

"""
    metabolite_compartment(m::MATModel, mid::String)

Extract metabolite compartment from `metCompartment` or `metCompartments`.
"""
metabolite_compartment(m::MATModel, mid::String) = _maybemap(
    x -> x[findfirst(==(mid), metabolites(m))],
    get(m.mat, "metCompartment", get(m.mat, "metCompartments", nothing)),
)

"""
    reaction_stoichiometry(model::MATModel, rxn_id::String)::Dict{String, Float64}

Return the reaction equation of reaction with id `rxn_id` in model. The reaction
equation maps metabolite ids to their stoichiometric coefficients.
"""
function reaction_stoichiometry(m::MATModel, rxn_id::String)::Dict{String,Float64}
    rxn_ind = first(indexin([rxn_id], m.mat["rxns"]))
    reaction_stoichiometry(m, rxn_ind)
end

"""
    reaction_stoichiometry(model::MATModel, rxn_index)::Dict{String, Float64}

Return the reaction equation of reaction with index `rxn_index` in model. The reaction
equation maps metabolite ids to their stoichiometric coefficients. Note, `rxn_index` can
be any suitable type that can index into an array.
"""
function reaction_stoichiometry(m::MATModel, rxn_ind)::Dict{String,Float64}
    met_inds = findall(m.mat["S"][:, rxn_ind] .!= 0.0)
    Dict(m.mat["mets"][met_ind] => m.mat["S"][met_ind, rxn_ind] for met_ind in met_inds)
end

# NOTE: There's no useful standard on how and where to store notes and
# annotations in MATLAB models. We therefore leave it very open for the users,
# who can easily support any annotation scheme using a custom wrapper.
# Even the (simple) assumptions about grRules, formulas and charges that we use
# here are very likely completely incompatible with >50% of the MATLAB models
# out there.

"""
    Base.convert(::Type{MATModel}, m::MetabolicModel)

Convert any metabolic model to `MATModel`.
"""
function Base.convert(::Type{MATModel}, m::MetabolicModel)
    if typeof(m) == MATModel
        return m
    end

    lb, ub = bounds(m)
    cl, cu = coupling_bounds(m)
    nr = n_reactions(m)
    nm = n_metabolites(m)
    return MATModel(
        Dict(
            "S" => stoichiometry(m),
            "rxns" => reactions(m),
            "mets" => metabolites(m),
            "lb" => Vector(lb),
            "ub" => Vector(ub),
            "b" => Vector(balance(m)),
            "c" => Vector(objective(m)),
            "C" => coupling(m),
            "cl" => Vector(cl),
            "cu" => Vector(cu),
            "genes" => genes(m),
            "grRules" =>
                _default.(
                    "",
                    _maybemap.(
                        x -> _unparse_grr(String, x),
                        reaction_gene_association.(Ref(m), reactions(m)),
                    ),
                ),
            "metFormulas" =>
                _default.(
                    "",
                    _maybemap.(
                        _unparse_formula,
                        metabolite_formula.(Ref(m), metabolites(m)),
                    ),
                ),
            "metCharges" => _default.(0, metabolite_charge.(Ref(m), metabolites(m))),
            "metCompartments" =>
                _default.("", metabolite_compartment.(Ref(m), metabolites(m))),
        ),
    )
end


text:
function parse_input()
    points = NTuple{4, Int}[]
    h = 0
    w = 0
    @inbounds for line in eachline("../data/input5.txt")
        p1, p2 = split(line, " -> ")
        x1, y1 = parse.(Int, split(p1, ","))
        x2, y2 = parse.(Int, split(p2, ","))
        w = max(x1, x2, w)
        h = max(y1, y2, h)
        push!(points, (x1, x2, y1, y2))
    end
    return points, w + 1, h + 1
end

function day5_part1()
    points, w, h = parse_input()
    occupied = zeros(Int, h+1, w+1)
    numcrowded = 0
    @inbounds for (x1, x2, y1, y2) in points
        if x1 == x2 || y1 == y2
            d = max(abs(x2 - x1), abs(y2 - y1))
            dx = sign(x2 - x1)
            dy = sign(y2 - y1)
            for step in 0:d
                xidx = x1 + dx * step + 1
                yidx = y1 + dy * step + 1
                occupied[yidx, xidx] += 1
                numcrowded += occupied[yidx, xidx] == 2
            end
        end
    end
    return numcrowded
end

function day5_part2()
    points, w, h = parse_input()
    occupied = zeros(Int, h+1, w+1)
    numcrowded = 0
    @inbounds for (x1, x2, y1, y2) in points
        d = max(abs(x2 - x1), abs(y2 - y1))
        dx = sign(x2 - x1)
        dy = sign(y2 - y1)
        for step in 0:d
            xidx = x1 + dx * step + 1
            yidx = y1 + dy * step + 1
            occupied[yidx, xidx] += 1
            numcrowded += occupied[yidx, xidx] == 2
        end
    end
    return numcrowded
end

max_stars_repo_name: lucaferranti/AoC

max_issues_repo_name: lucaferranti/AoC

max_forks_repo_name: lucaferranti/AoC


text:
using ContinuousDPs
using QuantEcon: PFI, VFI, solve
using BasisMatrices: Basis, ChebParams, SplineParams, LinParams
using Test

include("test_cdp.jl")


text:
# some model
mgrid = [range(0.0, stop=10.,step=0.05), range(0.0, stop=10.,step=0.05)];

model = Medium(mgrid);
vpb=[2100.,2200.];rhob=[2100., 2300.]
update!(model, [:vp,:rho], [vpb, rhob])
fill!(model)

nznx=prod(length.(model.mgrid))

for parameterization in [[:χKI, :χrhoI, :null],[:χKI, :null, :null],[:χvp, :χrho, :null],[:χvp, :null, :null]]
	δx1=randn(count(parameterization .≠ :null)*nznx)
	δxout1=zero(δx1)

	copyto!(δxout1, δx1, model, parameterization)

	δxout2=randn(size(δxout1));
	δx2=zero(δx1)
	GeoPhyInv.pert_chainrule!(δx2, δxout2, model, parameterization)

	@test LinearAlgebra.dot(δxout2,δxout1) ≈ LinearAlgebra.dot(δx1,δx2)
end


text:
@testset "461.hamming-distance.jl" begin
    @test hamming_distance(1, 4) == 2
    @test hamming_distance(235, 63626) == 8
end


text:
"""
    UniformBlockDiagonal{T}

Homogeneous block diagonal matrices.  `k` diagonal blocks each of size `m×m`
"""
struct UniformBlockDiagonal{T} <: AbstractMatrix{T}
    data::Array{T,3}
end

function Base.axes(A::UniformBlockDiagonal)
    m, n, l = size(A.data)
    return (Base.OneTo(m * l), Base.OneTo(n * l))
end

function Base.copyto!(dest::UniformBlockDiagonal{T}, src::UniformBlockDiagonal{T}) where {T}
    sdat = src.data
    ddat = dest.data
    size(ddat) == size(sdat) || throw(DimensionMismatch(""))
    copyto!(ddat, sdat)
    return dest
end

function Base.copyto!(dest::Matrix{T}, src::UniformBlockDiagonal{T}) where {T}
    size(dest) == size(src) || throw(DimensionMismatch(""))
    fill!(dest, zero(T))
    sdat = src.data
    m, n, l = size(sdat)
    @inbounds for k in axes(sdat, 3)
        ioffset = (k - 1) * m
        joffset = (k - 1) * n
        for j in axes(sdat, 2)
            jind = joffset + j
            for i in axes(sdat, 1)
                dest[ioffset + i, jind] = sdat[i, j, k]
            end
        end
    end
    return dest
end

function Base.getindex(A::UniformBlockDiagonal{T}, i::Int, j::Int) where {T}
    @boundscheck checkbounds(A, i, j)
    Ad = A.data
    m, n, l = size(Ad)
    iblk, ioffset = divrem(i - 1, m)
    jblk, joffset = divrem(j - 1, n)
    return iblk == jblk ? Ad[ioffset + 1, joffset + 1, iblk + 1] : zero(T)
end

function LinearAlgebra.Matrix(A::UniformBlockDiagonal{T}) where {T}
    return copyto!(Matrix{T}(undef, size(A)), A)
end

function Base.size(A::UniformBlockDiagonal)
    m, n, l = size(A.data)
    return (l * m, l * n)
end

"""
    BlockedSparse{Tv,S,P}

A `SparseMatrixCSC` whose nonzeros form blocks of rows or columns or both.

# Members
* `cscmat`: `SparseMatrixCSC{Tv, Int32}` representation for general calculations
* `nzasmat`: nonzeros of `cscmat` as a dense matrix
* `colblkptr`: pattern of blocks of columns

The only time these are created are as products of `ReMat`s.
"""
mutable struct BlockedSparse{T,S,P} <: AbstractMatrix{T}
    cscmat::SparseMatrixCSC{T,Int32}
    nzsasmat::Matrix{T}
    colblkptr::Vector{Int32}
end

function densify(A::BlockedSparse, threshold::Real=0.1)
    m, n = size(A)
    if nnz(A) / (m * n) ≤ threshold
        A
    else
        Array(A)
    end
end

Base.size(A::BlockedSparse) = size(A.cscmat)

Base.size(A::BlockedSparse, d) = size(A.cscmat, d)

Base.getindex(A::BlockedSparse, i::Integer, j::Integer) = getindex(A.cscmat, i, j)

LinearAlgebra.Matrix(A::BlockedSparse) = Matrix(A.cscmat)

SparseArrays.sparse(A::BlockedSparse) = A.cscmat

SparseArrays.nnz(A::BlockedSparse) = nnz(A.cscmat)

function Base.copyto!(L::BlockedSparse{T}, A::SparseMatrixCSC{T}) where {T}
    size(L) == size(A) && nnz(L) == nnz(A) ||
        throw(DimensionMismatch("size(L) ≠ size(A) or nnz(L) ≠ nnz(A"))
    copyto!(nonzeros(L.cscmat), nonzeros(A))
    return L
end

LinearAlgebra.rdiv!(A::BlockedSparse, B::Diagonal) = rdiv!(A.cscmat, B)

function LinearAlgebra.mul!(
    C::BlockedSparse{T,1,P},
    A::SparseMatrixCSC{T,Ti},
    adjB::Adjoint{T,BlockedSparse{T,P,1}},
    α,
    β,
) where {T,P,Ti}
    return mul!(C.cscmat, A, adjoint(adjB.parent.cscmat), α, β)
end


text:
using LoopVectorization, LinearAlgebra
BLAS.set_num_threads(1)

function jgemm_nkm!(C, A, B)
    C .= 0
    M, N = size(C); K = size(B,1)
    @inbounds for n ∈ 1:N, k ∈ 1:K
        @simd ivdep for m ∈ 1:M
            C[m,n] += A[m,k] * B[k,n]
        end
    end
end
function gemmavx!(C, A, B)
    @avx for i ∈ 1:size(A,1), j ∈ 1:size(B,2)
        Cᵢⱼ = 0.0
        for k ∈ 1:size(A,2)
            Cᵢⱼ += A[i,k] * B[k,j]
        end
        C[i,j] = Cᵢⱼ
    end
end
function jAtmulB!(C, A, B)
    @inbounds for n ∈ 1:size(C,2), m ∈ 1:size(C,1)
        Cₘₙ = zero(eltype(C))
        @simd ivdep for k ∈ 1:size(A,1)
            Cₘₙ += A[k,m] * B[k,n]
        end
        C[m,n] = Cₘₙ
    end
end
function jAtmulBavx!(C, A, B)
    @avx for n ∈ 1:size(C,2), m ∈ 1:size(C,1)
        Cₘₙ = zero(eltype(C))
        for k ∈ 1:size(A,1)
            Cₘₙ += A[k,m] * B[k,n]
        end
        C[m,n] = Cₘₙ
    end
end
function jdot(a, b)
    s = 0.0
    @inbounds @simd ivdep for i ∈ eachindex(a, b)
        s += a[i] * b[i]
    end
    s
end
function jdotavx(a, b)
    s = 0.0
    @avx for i ∈ eachindex(a, b)
        s += a[i] * b[i]
    end
    s
end
function jselfdot(a)
    s = 0.0
    @inbounds @simd ivdep for i ∈ eachindex(a)
        s += a[i] * a[i]
    end
    s
end
function jselfdotavx(a)
    s = 0.0
    @avx for i ∈ eachindex(a)
        s += a[i] * a[i]
    end
    s
end
function jdot3(x, A, y)
    M, N = size(A)
    s = zero(promote_type(eltype(x), eltype(A), eltype(y)))
    @inbounds for n ∈ 1:N
        @simd ivdep for m ∈ 1:M
            s += x[m] * A[m,n] * y[n]
        end
    end
    s
end
function jdot3avx(x, A, y)
    M, N = size(A)
    s = zero(promote_type(eltype(x), eltype(A), eltype(y)))
    @avx for m ∈ 1:M, n ∈ 1:N
        s += x[m] * A[m,n] * y[n]
    end
    s
end
function jvexp!(b, a)
    @inbounds for i ∈ eachindex(a)
        b[i] = exp(a[i])
    end
end
function jvexpavx!(b, a)
    @avx for i ∈ eachindex(a)
        b[i] = exp(a[i])
    end
end
function jsvexp(a)
    s = 0.0
    @inbounds for i ∈ eachindex(a)
        s += exp(a[i])
    end
    s
end
function jsvexpavx(a)
    s = 0.0
    @avx for i ∈ eachindex(a)
        s += exp(a[i])
    end
    s
end
function jgemv!(y, A, x)
    y .= 0.0
    for j ∈ eachindex(x)
        @simd ivdep for i ∈ eachindex(y)
            y[i] += A[i,j] * x[j]
        end
    end
end
function jgemvavx!(y, A, x)
    @avx for i ∈ eachindex(y)
        yᵢ = 0.0
        for j ∈ eachindex(x)
            yᵢ += A[i,j] * x[j]
        end
        y[i] = yᵢ
    end
end
function jvar!(s², A, x̄)
    @. s² = 0
    @inbounds for i ∈ 1:size(A,2)
        @simd for j ∈ eachindex(s²)
            δ = A[j,i] - x̄[j]
            s²[j] += δ*δ
        end
    end
end
function jvaravx!(s², A, x̄)
    @avx for j ∈ eachindex(s²)
        s²ⱼ = 0.0
        x̄ⱼ = x̄[j]
        for i ∈ 1:size(A,2)
            δ = A[j,i] - x̄ⱼ
            s²ⱼ += δ*δ
        end
        s²[j] = s²ⱼ
    end
end
japlucBc!(d, a, B, c) =      @. d = a + B * c';
japlucBcavx!(d, a, B, c) = @avx @. d = a + B * c';

function jOLSlp(y, X, β)
    lp = 0.0
    @inbounds @fastmath for i ∈ eachindex(y)
        δ = y[i]
        @simd for j ∈ eachindex(β)
            δ -= X[i,j] * β[j]
        end
        lp += δ * δ
    end
    lp
end
function jOLSlp_avx(y, X, β)
    lp = 0.0
    @avx for i ∈ eachindex(y)
        δ = y[i]
        for j ∈ eachindex(β)
            δ -= X[i,j] * β[j]
        end
        lp += δ * δ
    end
    lp
end


text:
abstract type BoundaryCondition end

struct Dirichlet <: BoundaryCondition end
struct Periodic <: BoundaryCondition end

const dirichlet = Dirichlet()
const periodic = Periodic()

has(::Dirichlet) = true
has(::Periodic) = false


text:
function Base.isapprox(a::Epoch{S}, b::Epoch{S}; kwargs...) where S <: TimeScale
    return isapprox(a.fraction + a.second, b.fraction + b.second; kwargs...)
end

function Base.:(==)(a::Epoch, b::Epoch)
    a.second == b.second && a.fraction == b.fraction
end

Base.isless(ep1::Epoch, ep2::Epoch) = isless(value(ep1 - ep2), 0.0)

function Base.:+(ep::Epoch{S}, p::AstroPeriod) where {S}
    second, fraction, error = apply_offset(ep.second, ep.fraction, ep.error, p.second, p.fraction, p.error)
    return Epoch{S}(second, fraction, error)
end
Base.:-(ep::Epoch, p::AstroPeriod) = ep + (-p)

"""
    -(a::Epoch, b::Epoch)

Return the duration between epoch `a` and epoch `b`.

### Examples ###

```jldoctest; setup = :(using AstroTime)
julia> TAIEpoch(2018, 2, 6, 20, 45, 20.0) - TAIEpoch(2018, 2, 6, 20, 45, 0.0)
20.0 seconds
```
"""
function Base.:-(a::Epoch{S}, b::Epoch{S}) where S<:TimeScale
    second = a.second - b.second
    fraction = (a.error - b.error) + (a.fraction - b.fraction)
    return (fraction + second) * seconds
end

max_stars_repo_name: JuliaAstrodynamics/AstronomicalTime

max_issues_repo_name: JuliaAstrodynamics/AstronomicalTime

max_forks_repo_name: JuliaAstrodynamics/AstronomicalTime


text:
# include("C:/Users/WorkingDirectory/SoluSL.jl")

using JLD2
using Random
using LinearAlgebra
using Statistics
using CSV
using DataFrames
using FreqTables

function q1()
#################### Q1 ###################
######## a
Random.seed!(1234)
A = 15 .*rand(10,7) .- 5;
B = 15 .*randn(10,7) .- 2;
C = hcat(A[ 1:5 , 1:5 ] , B[ 1:5 , 6:7 ])

latD = copy(vec(A))

for i in 1:length(latD)
    if latD[i] > 0
        latD[i] = 0
    end
end

D = reshape(latD,10,7)

######## b
length(A)

######## c
length(unique(D))

######## d
E = reshape(B,length(B),1)
Ev = vec(B)

######## e
F = cat(A,B,dims=3)

######## f
F = permutedims(F, [3,1,2])
######## g

G = kron(B,C)
# Julia rise dimensions issues
######## h
@save "C:/Users/WorkingDirectory/matrixpractice.jld" A B C D E F G

######## i
@save "C:/Users/WorkingDirectory/firstmatrix.jld" A B C D

######## j
CSV.write("C:/Users/WorkingDirectory/Cmatrix.csv",DataFrame(C))

######## k
CSV.write("C:/Users/WorkingDirectory/Dmatrix.dat",DataFrame(D), delim="\t")

######## l

return(A,B,C,D)
end

A,B,C,D = q1()

function q2(A,B,C)

#################### Q2 ###################
######## a

AB = Array{Float64}(undef,10,7)

for i in 1:size(A,1), j = 1:size(B,2)
    AB[i,j] = A[i,j]*B[i,j]
end

AB2 = A .*B
######## b

MaxElementC = size(C,1)*size(C,2)
ll = 0
CprimeI = Array{Float64}(undef,MaxElementC,1)
for i in 1:size(C,1), j = 1:size(C,2)
    if  -5 <= C[i,j] <= 5
        ll = ll + 1
        CprimeI[ll] = C[i,j]
    end
end

Cprime = CprimeI[1:ll]
Cprime2 = C[abs.(C) .<= 5]
######## c

N, K, T = 15169, 6, 5
X = Array{Float64}(undef,N,K,T)
for t in 1:T
    X[:,1,t]  = ones(N)
    X[:,2,t]  = [Int64( rand() >= (0.75 *(6 - t)/5) ) for i in 1:N ]
    X[:,3,t]  = (5*(t-1)) .*randn(N) .+ (15 + t + 1)
    X[:,4,t]  = (1/exp(1)) .*randn(N) .+ pi*(6 - t)/3
    X[:,5,t]  = trunc.((2.19^3) .*rand(N) .+ 12)
    X[:,6,t]  = trunc.((20*0.5^2) .*rand(N) .+ (20*0.5))
end

######## d
Beta = Array{Float64}(undef,K,T)
[Beta[:,t] = [0.25*t + 0.75, log(t), -sqrt(t), exp(t) - exp(t + 1), t, t/3] for t in 1:T]

######## e
Y = Array{Float64}(undef,N,T)
[Y[:,t] =  X[:,:,t]*Beta[:,t] + 0.36*randn(N) for t in 1:T]

end

q2(A,B,C)

function q3()
#################### Q3 ###################
######## a
dff = CSV.read("C:/Users/WorkingDirectory/nlsw88.csv"; delim = ',') # in this data missing values are not ignored

df = CSV.read("C:/Users/WorkingDirectory/nlsw88.csv"; delim = ',', missingstring = "NA")
@save "C:/Users/WorkingDirectory/00-JuliaTips/nlsw88.jld" df

######## b
per_nm = mean(df.never_married)*100
per_cg = mean(df.collgrad)*100

######## c
per_race = prop(freqtable(df.race)).*100
######## d

sumstdf = describe(df, :mean, :median, :std, :min, :max, :nunique, :q75, :q25)
size(sumstdf)

m = size(sumstdf)[2]-2 # last col for interquartile computing
summarystats = Array{Any}(undef,size(sumstdf)[1],m)

#initiating matrix with the describe df values and replacing nothing values by missing which will allows computing

for i in 2:m
    summarystats[:,i-1] = replace(sumstdf[:,i], nothing => missing)
end

#computing interquartile
summarystats[:,m] = replace(sumstdf[:, m + 1], nothing => missing) - replace(sumstdf[:, m + 2], nothing => missing)

#from the describe(dff) command, there is 2 missing values
describe(dff)

######## e
freqtable(df.industry, df.occupation)[1:end .!=1 , 1:end .!=1 ]  # droping row and col of missing values

######## f
subdf = select(df, :wage, :occupation, :industry)
combine(groupby(subdf, All(:industry, :occupation)), :wage => mean)

end

function q4()
#################### Q4 ###################
######## a
@load "C:/Users/WorkingDirectory/firstmatrix.jld"
######## b
function matrixops(A,B)
    #compute and return a tulpe element-by-element product, matrix product and sum of elements of two matrix
    if(size(A)!=size(B))
        return("inputs must have the same size")
    end

    i = A .*B
    ii = A'B
    iii = sum(A + B)

    return(i,ii,iii)
end

######## c
# in the code above

######## d
matrixops(q1()[1],q1()[2])

######## e
# in the code above

######## f
matrixops(q1()[3],q1()[4])
# it rises the error message

######## g
convert(Array,df.ttl_exp)
convert(Array,df.wage)
matrixops(convert(Array,df.ttl_exp),convert(Array,df.wage))

######## h

end


text:
#All series by default
function smooth!(ec::EChart; smooth::Bool = true)

	[x.smooth = smooth for x in ec.series]

	return ec

end

#Single series
function smooth!(ec::EChart, series::Int; smooth::Bool = true)

	ec.series[series].smooth = smooth

	return ec

end

#Specified Series in Array
function smooth!(ec::EChart, series::Vector{Int}; smooth::Bool = true)

	[smooth!(ec, x, smooth = smooth) for x in series]

	return ec

end


text:
module SyringePumpJl
    using Measurements, Unitful,  Blink, Interact, LibSerialPort, Plots

    include("UI.jl")
    export InputField, UncertainInputField, value, Calibration,
            makeplot, launch_pump_interface

    include("RXTX.jl")
    export drive_pump, stop_pump, drive_physical
end # module

max_stars_repo_name: JuliaMakers/MakerProjects

max_issues_repo_name: JuliaMakers/MakerProjects

max_forks_repo_name: JuliaMakers/JuliaMakerProjects


text:
"""
Problem Overview:
A UAV chases a moving ground platform (the target) and intends to finally land on it.
Assume full knowledge of the UAV state with the only uncertainty lying on the target state.
A camera installed on the UAV receives intermittent observations of the target to estimate its state.
Assume a noisy measurement of the target location (2D) when it is inside the camera field-of-view,
and no measurement when it is outside.

"""

using POMDPs
using Random # for AbstractRNG
using POMDPModelTools # for Deterministic
using Parameters
using StaticArrays
using Distributions
using LinearAlgebra
using Base
using ParticleFilters
using POMCPOW
using MCTS
using Test
using POMDPSimulators
using POMDPPolicies
using Plots
# using PyPlot
using Reel
using ProgressMeter
# importall POMDPs

const Vec2 = SVector{2, Float64}
const Vec3 = SVector{3, Float64}

struct MDPState
    uavPose::Vec3
    uavHeading::Float64 # degree
    targetPose::Vec2
    targetPose_true::Vec2
end

struct MDPAction
    xy_speed::Float64 # m/s
    z_speed::Float64
    angle::Float64 # degree
end

struct POMDPAction
    look::Bool
    vel_steer::MDPAction
end

mutable struct targetObservationDistribution # not an actual distribution. Probably better to revise the name
    current_observation::Vec2
    if_observed::Bool
end
POMDPs.rand(rng::AbstractRNG, d::targetObservationDistribution) = d
# POMDPs.pdf(d::TMazeObservationDistribution, o::Int64) = o == d.current_observation ? (return 1.0) : (return 0.0)

@with_kw struct UAVchaseMDP <: MDP{MDPState, MDPAction}
    # mu::Float64          = 2.0
    target_velocity::Vec2 = SVector(0.2, 0.0)
    target_std::Float64     = 0.04# already very good: 0.01

    dt::Float64          = 0.1
    target_height::Float64 = 0.3 # height of the platform
    landing_radius::Float64  = 0.2 # old is 0.8
    ready_to_land_distance = 5.0 # ready for landing when target is this distance away from the origin
    UAV_xy_speed = 0.6
    UAV_z_speed = 0.2
    # initial conditions:
    init_UAVPose = SVector(0,0,0.5)
    init_UAVHeading = 0.0
    init_targetPose = SVector(3., -1.0)
    init_targetPose_true = SVector(2., 2.0)
    init_target_std = 5.
    # physical constraints:
    z_min::Float64 = 0.3
    z_max::Float64 = 5
    # rewards:
    r_undetected::Float64 = -800.
    r_outScene::Float64 = -100.0 # under the ground or above the ceiling
    r_action::Float64 = -10.0
    r_distance::Float64 = -300.

    discount::Float64    = 0.95
end

@with_kw struct UAVchasePOMDP <: POMDP{MDPState, MDPAction, targetObservationDistribution}
    mdp::UAVchaseMDP           = UAVchaseMDP()
    # meas_std::Float64          = 0.001 # camera measurement noise
    meas_std::Float64          = 0.1

    # angle to be within Field of View (FOV)
    FOV_angle1::Float64 = 40 # in degree
    FOV_angle2::Float64 = 40 # in degree
end

const UAVchaseProblem = Union{UAVchasePOMDP, UAVchaseMDP}
mdp(p::UAVchaseMDP) = p
mdp(p::UAVchasePOMDP) = p.mdp

POMDPs.discount(pp::UAVchasePOMDP) = mdp(pp).discount

function POMDPs.isterminal(pp::UAVchasePOMDP, s::MDPState)
    p = mdp(pp)
    target3Dpose = SVector(s.targetPose_true[1], s.targetPose_true[2], p.target_height)
    distance = norm(s.uavPose - target3Dpose)
    return (norm(target3Dpose) > p.ready_to_land_distance) && (distance < p.landing_radius)
end

mutable struct MDPStateDistribution
    current_uavPose::Vec3
    current_uavHeading::Float64
    current_targetPose_true::Vec2
    mean_targetPose::Vec2
    std_targetPose::Float64
end
POMDPs.sampletype(::Type{MDPStateDistribution}) = MDPState

function POMDPs.rand(rng::AbstractRNG, d::MDPStateDistribution)
    target_cov = Matrix(d.std_targetPose*Diagonal{Float64}(I, 2))
    rand_targetPose = d.mean_targetPose + rand(rng, MvNormal(target_cov))
    return MDPState(d.current_uavPose, d.current_uavHeading, rand_targetPose,  d.current_targetPose_true)
end

function POMDPs.initialstate_distribution(pp::UAVchaseProblem)
    p = mdp(pp)
    return MDPStateDistribution(p.init_UAVPose, p.init_UAVHeading, p.init_targetPose_true, p.init_targetPose, p.init_target_std)
end

function POMDPs.generate_s(pp::UAVchaseProblem, s::MDPState, a::MDPAction, rng::AbstractRNG)
    p = mdp(pp)
    # calculate target state
    target_dt_distance = p.dt*p.target_velocity
    ## sensor_noise = Base.rand(Normal(0, p.target_std), 2)
    target_dynamic_cov = Matrix(p.target_std*Diagonal{Float64}(I, 2))
    sensor_noise = rand(rng, MvNormal(target_dynamic_cov))
    ## sensor_noise = SVector(0,0)
    curr_targ_true = s.targetPose_true + target_dt_distance
    if if_in_FOV(pp, s)
        curr_targ = curr_targ_true + sensor_noise
    else
        curr_targ = s.targetPose + target_dt_distance + sensor_noise # next_target_pos(p, s.targetPose)
    end
    # calculate UAV state
    curr_angle = s.uavHeading + a.angle
    xy_dt_distance = p.dt*a.xy_speed*SVector(cosd(curr_angle), sind(curr_angle)) # careful
    z_dt_distance = p.dt*a.z_speed
    xyz_dt_distance = SVector(xy_dt_distance[1], xy_dt_distance[2], z_dt_distance)
    curr_pos = s.uavPose + xyz_dt_distance
    return MDPState(curr_pos, curr_angle, curr_targ, curr_targ_true)
end

function if_in_FOV(p::UAVchasePOMDP, s::MDPState)
    # pose = s.uavPose
    range1 = s.uavPose[3]*tand(p.FOV_angle1)
    range2 = s.uavPose[3]*tand(p.FOV_angle2)
    # within_x = (s.targetPose[1]>(s.uavPose[1]-range1))||(s.targetPose[1]<(s.uavPose[1]+range1))
    if !((s.targetPose_true[1]>(s.uavPose[1]-range1))&&(s.targetPose_true[1]<(s.uavPose[1]+range1)))
        return false # if target is NOT within FOV in x direction
    end
    if (s.targetPose_true[2]>(s.uavPose[2]-range2))&&(s.targetPose_true[2]<(s.uavPose[2]+range2))
        return true # if target is within FOV in y direction
    else
        return false
    end
end

function if_in_FOV_fake(p::UAVchasePOMDP, s::MDPState)
    # pose = s.uavPose
    range1 = s.uavPose[3]*tand(p.FOV_angle1)
    range2 = s.uavPose[3]*tand(p.FOV_angle2)
    # within_x = (s.targetPose[1]>(s.uavPose[1]-range1))||(s.targetPose[1]<(s.uavPose[1]+range1))
    if !((s.targetPose[1]>(s.uavPose[1]-range1))&&(s.targetPose[1]<(s.uavPose[1]+range1)))
        return false # if target is NOT within FOV in x direction
    end
    if (s.targetPose[2]>(s.uavPose[2]-range2))&&(s.targetPose[2]<(s.uavPose[2]+range2))
        return true # if target is within FOV in y direction
    else
        return false
    end
end
# I only count the case of POMDP
function POMDPs.reward(pp::UAVchasePOMDP, a::MDPAction, sp::MDPState)
    p = mdp(pp)
    # action_reward = p.r_action * ((a.xy_speed != 0)||(a.z_speed != 0))
    if ((a.xy_speed != 0)||(a.z_speed != 0))
        action_reward = p.r_action
    else
        action_reward = 0.0
    end
    # detection_reward = p.r_detect * if_in_FOV(pp, sp)
    if if_in_FOV(pp, sp)
        detection_reward = 0.0
    else
        detection_reward = p.r_undetected
    end
    # outScene_reward = p.r_outScene * ((sp.uavPose[3] < p.z_min)||(sp.uavPose[3] > p.z_max))
    if ((sp.uavPose[3] < p.z_min)||(sp.uavPose[3] > p.z_max))
        outScene_reward = p.r_outScene
    else
        outScene_reward = 0.0
    end
    # distance reward:
    target3Dpose = SVector(sp.targetPose[1], sp.targetPose[2], p.target_height)
    distance_reward = p.r_distance*(norm(sp.uavPose - target3Dpose)+0.00001)
    return action_reward + detection_reward + outScene_reward + distance_reward
end

function POMDPs.generate_sr(pp::UAVchasePOMDP, s::MDPState, a::MDPAction, rng::AbstractRNG)
    sp = generate_s(pp, s, a, rng)
    return sp, reward(pp,a, sp)
end

function POMDPs.generate_o(pp::UAVchasePOMDP, s::MDPState, a::MDPAction, sp::MDPState, rng::AbstractRNG)
    p = mdp(pp)
    if if_in_FOV(pp, sp) # if target is IN camera field of view (FOV) <=> having a measurement
        meas_cov = Matrix(pp.meas_std*Diagonal{Float64}(I, 2))
        measurement_noise = rand(rng, MvNormal(meas_cov))
        targetPose_true_pred = s.targetPose_true + p.dt*p.target_velocity
        meas_location = targetPose_true_pred + measurement_noise # measurement is the target location + noise
        return targetObservationDistribution(meas_location, true)
    else # not in FOV <=> no measurement available
        targetPose_pred = s.targetPose + p.dt*p.target_velocity # a naive predictor when no measurement is available when taget is NOT within FOV
        return targetObservationDistribution(targetPose_pred, false)
    end
end

function POMDPs.observation(pp::UAVchasePOMDP, s::MDPState, a::MDPAction, sp::MDPState, rng::AbstractRNG)
    return generate_o(pp, s, a, sp, rng) # not sure if this function is needed for anything
end

function POMDPModelTools.obs_weight(pp::UAVchasePOMDP, s::MDPState, a::MDPAction, sp::MDPState, o::targetObservationDistribution)
     p = mdp(pp)
     if if_in_FOV(pp, sp)
         if o.if_observed
             y = o.current_observation # true measurement
             y_hat = s.targetPose + p.dt*p.target_velocity # not sure if prediction should be sp
             sigma_square = pp.meas_std
             return 1.0/(2*pi*sigma_square)*exp(-0.5*( (y[1]-y_hat[1])^2 + (y[2]-y_hat[2])^2 )/sigma_square) # assume normal distribution
         else

             return 0.0
         end
     else # if not in FOV
         if o.if_observed
             return 0.0
         else
             return 1.0
         end
     end
 end

function POMDPs.actions(pp::UAVchasePOMDP)
    p = mdp(pp)
    xy = p.UAV_xy_speed
    z = p.UAV_z_speed
    action_list= [MDPAction(xy,0.,0.), MDPAction(xy,0.,15.), MDPAction(xy,0.,30.), MDPAction(xy,0.,45.),
                  MDPAction(0.,0.,0.), MDPAction(xy,0.,-15.), MDPAction(xy,0.,-30.), MDPAction(xy,0.,-45.),
                  MDPAction(0.,z*2/3,0.), MDPAction(0.,-z, 0.), MDPAction(xy,0.,90.), MDPAction(xy,0.,-90.),
                  MDPAction(2*xy,0.,0.), MDPAction(2*xy,0.,15.), MDPAction(2*xy,0.,30.), MDPAction(2*xy,0.,45.),
                  MDPAction(0.,-3*z,0.), MDPAction(2*xy,0.,-15.), MDPAction(2*xy,0.,-30.), MDPAction(2*xy,0.,-45.),
                  MDPAction(0.,-2*z, 0.), MDPAction(2*xy,0.,90.), MDPAction(2*xy,0.,-90.),
                  MDPAction(-xy/2,0.,0.),MDPAction(0.,10*z,0.), MDPAction(0.,20*z,0.), MDPAction(0.,50*z,0.),
                  MDPAction(xy/2,0.,0.), MDPAction(xy/2,0.,-15.), MDPAction(xy/2,0.,30.), MDPAction(xy/2,0.,-45.),
                  MDPAction(xy/3,0.,0.), MDPAction(xy/3,0.,15.), MDPAction(xy/3,0.,-30.), MDPAction(xy/3,0.,45.),
                  MDPAction(xy/6,0.,0.), MDPAction(xy/6,0.,-30.), MDPAction(xy/6,0.,30.),MDPAction(-xy/6,0.,0.)]
    return action_list
end


text:
using Dates

mutable struct NPLCitation
    id::String
    num::Int64
    text::String
    external_ids::Vector{String}
    NPLCitation() = new("", 0, "", String[])
end

mutable struct Title
    text::String
    lang::String
    Title() = new("", "")
end

mutable struct Abstract
    text::String
    lang::String
    Abstract() = new("", "")
end

lang(t::Title) = t.lang
text(t::Title) = t.text

lang(t::Abstract) = t.lang
text(t::Abstract) = t.text

mutable struct Classification
    symbol::String
    title::String
    Classification() = new("", "")
end

code(c::Classification) = c.symbol
subgroup(c::Classification) = code(c)
maingroup(c::Classification) = match(r"[^/]*", code(c)).match
subclass(c::Classification) = first(code(c), 4)
class(c::Classification) = first(code(c), 3)
section(c::Classification) = first(code(c), 1)

==(c1::Classification, c2::Classification) = code(c1) == code(c2)
hash(class::Classification, h::UInt) = hash(code(class), h)

mutable struct ApplicationID
    source::String
    id::String
    jurisdiction::String
    doc_number::String
    kind::String
    date::Date
    ApplicationID() = new("", "", "", "", "", Date(9999))
    ApplicationID(s, i, j, d, k ,dt) = new(s, i, j, d, k, dt)
end

source(a::ApplicationID) = a.source
id(a::ApplicationID) = a.id
jurisdiction(a::ApplicationID) = a.jurisdiction
docnr(a::ApplicationID) = a.doc_number
kind(a::ApplicationID) = a.kind
date(a::ApplicationID) = a.date

==(a1::ApplicationID, a2::ApplicationID) = id(a1) == id(a2)
hash(app::ApplicationID, h::UInt) = hash(id(app), h)

mutable struct Application
    id::ApplicationID
    publication_type::String
    inventors::Vector{String}
    applicants::Vector{String}
    title::Vector{Title}
    abstract::Vector{Abstract}
    cpc::Vector{Classification}
    siblings_simple::Vector{ApplicationID}
    cites_patents::Vector{ApplicationID}
    cites_npl::Vector{NPLCitation}
    cited_by::Vector{ApplicationID}
    family_size::Int
    cites_patents_count::Int
    cites_npl_count::Int
    cited_by_count::Int
    Application() = new(ApplicationID(), "", [""], [""], Title[], Abstract[],
                        Classification[], ApplicationID[], ApplicationID[],
                        NPLCitation[], ApplicationID[], 0, 0, 0, 0)

    Application(id, typ, inv, app, tit, abs, cpc, fam, cit, npl,
                citby, fsize, citcount, citcountnpl, citbycount) = begin

        new(id, typ, inv, app, tit, abs, cpc, fam, cit, npl,
            citby, fsize, citcount, citcountnpl, citbycount
        )
    end
end

id(a::Application) = id(a.id)
jurisdiction(a::Application) = jurisdiction(a.id)
docnr(a::Application) = docnr(a.id)
kind(a::Application) = kind(a.id)
date(a::Application) = date(a.id)
status(a::Application) = a.publication_type

title(a::Application) = a.title
abstract(a::Application) = a.abstract

inventors(a::Application) = a.inventors
applicants(a::Application) = a.applicants

siblings(a::Application) = a.siblings_simple
familysize(a::Application) = a.family_size

classification(a::Application) = a.cpc

cites(a::Application) = a.cites_patents
cites_npl(a::Application) = a.cites_npl
cites_count(a::Application) = a.cites_patents_count
cites_count_npl(a::Application) = a.cites_npl_count

citedby(a::Application) = a.cited_by
citedby_count(a::Application) = a.cited_by_count

==(a1::Application, a2::Application) = id(a1) == id(a2)
hash(app::Application, h::UInt) = hash(id(app), h)

function Base.show(io::IO, a::ApplicationID)
    print(io, "$(id(a)) | $(date(a)) | $(jurisdiction(a) * docnr(a) * kind(a))" )
end

function Base.show(io::IO, a::Application)
    print(io, "$(id(a)) | $(date(a)) | $(jurisdiction(a) * docnr(a) * kind(a))" )
end

Base.show(io::IO, t::Title) = print(io, "$(text(t))" )
Base.show(io::IO, a::Abstract) = print(io, "$(text(a))" )


text:
# the dependecy scheduler
module ArbiterSchedulers

import ..ArbiterGraphs: ArbiterGraph, NodeSet, ImmutableNodeSet
import ..ArbiterGraphs: roots, nodes, prune!
import ..ArbiterGraphs: Strategy, Orphan, Promote, Remove
import ..ArbiterTasks: ArbiterTask

export  Scheduler,
        completed,
        failed,
        running,
        runnable,
        isfinished,
        add_task!,
        start_task!,
        end_task!,
        remove_unrunnable!,
        fail_remaining!

type Scheduler
    _graph::ArbiterGraph
    _tasks::Dict{Symbol, ArbiterTask}
    _running::NodeSet
    _completed::NodeSet
    _failed::NodeSet
end

"""
A dependency scheduler
"""
function Scheduler(tasks=NodeSet(), completed=NodeSet(), failed=NodeSet())
    scheduler = Scheduler(
        ArbiterGraph(),
        Dict{Symbol, ArbiterTask}(),
        NodeSet(),
        completed,
        failed,
    )

    for task in tasks
        add_task!(scheduler, task)
    end

    return scheduler
end

"""
A copy of the set of successfully completed tasks.
"""
completed(s::Scheduler) = ImmutableNodeSet(s._completed)

"""
A copy of the set of failed tasks.
"""
failed(s::Scheduler) = ImmutableNodeSet(s._failed)

"""
A copy of the set of running tasks.
"""
running(s::Scheduler) = ImmutableNodeSet(s._running)

"""
Get the set of tasks that are currently runnable.
"""
runnable(s::Scheduler) = setdiff(roots(s._graph), s._running)

"""
Have all runnable tasks completed?
"""
isfinished(s::Scheduler) = isempty(roots(s._graph)) && isempty(s._running)

"""
Add a task to the scheduler.

task: The task to add.
"""
function add_task!(scheduler::Scheduler, task::ArbiterTask)
    scheduler._tasks[task.name] = task

    incomplete_dependencies = NodeSet()
    failed = false

    for dependency in task.dependencies
        if dependency in scheduler._failed
            cascade_failure!(scheduler, task.name)
            failed = true
            break
        end

        if !(dependency in scheduler._completed)
            push!(incomplete_dependencies, dependency)
        end
    end

    if !failed
        try
            push!(scheduler._graph, task.name, incomplete_dependencies)
        catch
            cascade_failure!(scheduler, task.name)
        end
    end

    return nothing
end

"""
Start a task.

Returns the task that was started (or None if no task has been
    started).

name: (optional, None) The task to start. If a name is given,
    Scheduler will attempt to start the task (and raise an
    exception if the task doesn't exist or isn't runnable). If
    no name is given, a task will be chosen arbitrarily
"""
function start_task!(scheduler::Scheduler, name::Symbol)
    if !(name in scheduler._graph._roots) || name in scheduler._running
        error(name)
    end

    push!(scheduler._running, name)

    return Nullable(scheduler._tasks[name])
end

start_task!(scheduler::Scheduler, name) = start_task!(scheduler, symbol(name))

function start_task!(scheduler::Scheduler)
    found = false
    local name

    for possibility in scheduler._graph._roots
        if !(possibility in scheduler._running)
            name = possibility
            found = true
            break
        end
    end

    if found
        push!(scheduler._running, name)
        task = Nullable(scheduler._tasks[name])
    else
        # all tasks blocked/running/completed/failed
        task = Nullable{ArbiterTask}()
    end

    return task
end

"""
End a running task. Raises an exception if the task isn't
running.

name: The name of the task to complete.
success: (optional, True) Whether the task was successful.
"""
function end_task!(scheduler::Scheduler, name::Symbol, success=true)
    pop!(scheduler._running, name)

    if success
        push!(scheduler._completed, name)
        pop!(scheduler._graph, name; strategy=Orphan)
    else
        cascade_failure!(scheduler, name)
    end

    return nothing
end

end_task!(scheduler::Scheduler, name, success=true) = end_task(scheduler, symbol(name), success)

"""
Remove any tasks that are dependent on non-existent tasks.
"""
function remove_unrunnable!(scheduler::Scheduler)
    union!(scheduler._failed, prune!(scheduler._graph))
end

"""
Mark all unfinished tasks (including currently running ones) as
failed.
"""
function fail_remaining!(scheduler::Scheduler)
    union!(scheduler._failed, nodes(scheduler._graph))
    scheduler._graph = ArbiterGraph()
    scheduler._running = NodeSet()
end

"""
Mark a task (and anything that depends on it) as failed.

name: The name of the offending task
"""
function cascade_failure!(scheduler::Scheduler, name)
    if name in scheduler._graph
        union!(scheduler._failed, pop!(scheduler._graph, name; strategy=Remove))
    else
        push!(scheduler._failed, name)
    end

    return nothing
end

"""
Remove all unrunnable tasks and call a function. When
the function is exited, all non-complete tasks will be
failed.
"""
function Scheduler(f::Function, args...)
    scheduler = Scheduler(args...)

    remove_unrunnable!(scheduler)
    try
        f(scheduler)
    finally
        fail_remaining!(scheduler)
    end

    return nothing
end

end


text:
include("meta.jl")

mutable struct Diffable #<: Real

  value::Float64
  precursors::Vector{Tuple{Float64,Diffable}} #gradient weight + the Diffable node that is dependent during a forward pass of the computation graph (antecedent from a backwards perspective, subsequent from a forward perspective)
  grad

  Diffable(a) = new(a, Vector{Tuple{Float64,Diffable}}(), 0.0)
end

function clear!(x::Diffable)

  x.precursors = similar(x.precursors, 0)
end

"""sadly 8x faster than a cooler one-liner"""
function backwards!(y::Diffable)

    if isempty(y.precursors) # terminal node
      y.grad = 1
    else
      y.grad = 0 # remove this to accumulate
      for (weight,precursor) in y.precursors
        y.grad += weight * backwards!(precursor)
      end
    end

  return y.grad
end

"""
generate conversion functions for a binary operator
e.g. `@create_binaryfns +`
yields
`+(a::Diffable, b) = +(a, Diffable(b))`
`+(a, b::Diffable) = +(Diffable(a), b)`
"""
macro create_conversions(operator)

  quote
    $(esc(operator))(a::Diffable, b::T) where T <: Real = $(esc(operator))(a, Diffable(b))
    $(esc(operator))(a::T, b::Diffable) where T <: Real = $(esc(operator))(Diffable(a), b)
  end
end

"""
generate base functions for a binary operator in terms of arguments `a` and `b`

calculate the result and push the result (weight) and partial derivative into the computation graph

example output:
```
function *(a::Diffable, b::Diffable)

  result = Diffable(a.value * b.value)
  push!(a.precursors, (b.value, result))
  push!(b.precursors, (a.value, result))

  return result
end
```
"""
macro create_binaryfn(operator, ∂a, ∂b, a = :a, b = :b)

  ∂a = replace_expr!(∂a, a => :($a.value), b => :($b.value))
  ∂b = replace_expr!(∂b, a => :($a.value), b => :($b.value))

  quote
    function $(esc(operator))(a::Diffable, b::Diffable)

      result = Diffable($(esc(operator))($a.value, $b.value))
      push!($(:(a.precursors)), ($∂a, result))
      push!($(:(b.precursors)), ($∂b, result))

      return result
    end

    @create_conversions($(esc(operator)))
  end
end

"""unary version"""
macro create_unaryfn(fn, ∂a, a = :a)

  ∂a = replace_expr!(∂a, a => :($a.value))

  quote
    function $(esc(fn))(a::Diffable)

      result = Diffable($(esc(fn))($a.value))
      push!($(:(a.precursors)), ($∂a, result))

      return result
    end
  end

end

# actually define the adjoints:

# binary functions
@create_binaryfn(+, 1.0, 1.0)
@create_binaryfn(-, 1.0, -1.0)
@create_binaryfn(*, b, a)
@create_binaryfn(^, b * a^(b-1), a^b * log(a)) #x^a, a^x

/(a::Diffable, b::Diffable) = a * (b^-1)
@create_conversions /

# unary functions
-(a::Diffable) = -1 * a

@create_unaryfn(inv, -a * a^-2) #used internally for calls to a^-1
@create_unaryfn(exp, result.value) #short-circuit recalculating exp(a)
@create_unaryfn(log, 1/a)
@create_unaryfn(sin, cos(a))
@create_unaryfn(cos, -sin(a))
@create_unaryfn(abs, sign(a)) #equivalent to |x|/x, but handles 0 to avoid random errors at the expense of correctness.

relu(x) = x > 0 ? x : 0
@create_unaryfn(relu, a > 0)


text:
using MeasureTheory
using ConcreteStructs

import MeasureTheory: logdensity

@concrete terse struct Chain <: AbstractMeasure
    κ
    μ
end

evolve(mc::MarkovChain, μ) =  μ ⋅ mc.κ
evolve(mc::MarkovChain, ::Nothing) =  mc.μ

function dyniterate(mc::MarkovChain, u::Sample)
    xnew =  rand(u.rng, evolve(mc, Dirac(u.x)))
    xnew, Sample(u.rng, xnew)
end

function logdensity(mc::Chain, x)
    μ = mc.μ
    ℓ = 0.0
    for xj in x
        ℓ += logdensity(μ, xj)
        μ = mc.κ(xj)
    end
    return ℓ
end

mc = Chain(Normal()) do x Normal(μ=x) end

logdensity(mc, randn(100))

using Soss

hmm = @model begin
    ε ~ Exponential() #  transition
    σ ~ Exponential() # Observation noise
    x ~ Chain(Normal()) do xj
        Normal(xj, ε)
    end

    y ~ For(x) do xj
        Normal(xj, σ)
    end
end

using Soss

mbind = @model μ,next begin
    x ~ μ
    y ~ next(x)
    return y
end;

mbind2 = @model μ,next begin
    x ~ μ
    y ~ next(last(x))
    return (x...,y)
end;

⋅(μ,next) = mbind(μ,next)

⊛(μ,next) = mbind2(μ,next);

next(x) = Normal(x+1,1);

d =  Normal() ⋅ next ⋅ next ⋅ next

μ0 = @model d begin
    x ~ d
    return (x,)
end

d2 = (Normal()^1) ⊛ next ⊛ next ⊛ next

rand(d)
rand(d2)

simulate(d)

t = xform(d);
t(randn(4))

# julia> d =  Cauchy() ⋅ (x -> Normal(μ=x)) ⋅ (x -> Normal(μ=x)) ⋅ (x -> Normal(μ=x))
# ConditionalModel given
#     arguments    (:μ, :κ)
#     observations ()
# @model (μ, κ) begin
#         x ~ μ
#         y ~ κ(x)
#         return y
#     end

# julia> rand(d)
# -3.0414465047589037

# julia> t = xform(d)
# TransformVariables.TransformTuple{NamedTuple{(:x, :y), Tuple{TransformVariables.TransformTuple{NamedTuple{(:x, :y), Tuple{TransformVariables.TransformTuple{NamedTuple{(:x, :y), Tuple{TransformVariables.Identity, TransformVariables.Identity}}}, TransformVariables.Identity}}}, TransformVariables.Identity}}}((x = TransformVariables.TransformTuple{NamedTuple{(:x, :y), Tuple{TransformVariables.TransformTuple{NamedTuple{(:x, :y), Tuple{TransformVariables.Identity, TransformVariables.Identity}}}, TransformVariables.Identity}}}((x = TransformVariables.TransformTuple{NamedTuple{(:x, :y), Tuple{TransformVariables.Identity, TransformVariables.Identity}}}((x = asℝ, y = asℝ), 2), y = asℝ), 3), y = asℝ), 4)

# julia> t(randn(4))
# (x = (x = (x = -0.24259286698966315, y = 0.278190893626807), y = -1.361907586870645), y = 0.05914265096096323)

# julia> simulate(d)
# (value = 5.928939554009484, trace = (x = (value = 5.307006358072237, trace = (x = (value = 3.2023770380851797, trace = (x = 3.677550124255551, y = 3.2023770380851797)), y = 5.307006358072237)), y = 5.928939554009484))


text:
# Test Knitro with MINLPTests
#
using MINLPTests, JuMP, Test

using KNITRO

const SOLVER = JuMP.with_optimizer(KNITRO.Optimizer, outlev=0, opttol=1e-8)

const NLP_SOLVERS = [SOLVER]
const MINLP_SOLVERS = [SOLVER]
const POLY_SOLVERS = []
const MIPOLY_SOLVERS = []

@testset "JuMP Model Tests" begin
    @testset "$(solver.constructor): nlp" for solver in NLP_SOLVERS
        MINLPTests.test_nlp(solver, exclude = [
            "005_010",
            "005_011",  # Uses the function `\`
        ])
        # For 005_010, Knitro founds a different solution, close
        # to those of MINLPTests.
        MINLPTests.nlp_005_010(solver, 1e-5, 1e-5, 1e-5)
        MINLPTests.test_nlp_cvx(solver)
    end
    @testset "$(solver.constructor): nlp_mi" for solver in MINLP_SOLVERS
        MINLPTests.test_nlp_mi(solver, exclude = [
            "005_010",
            "005_011",  # Uses the function `\`
            "006_010"   # Bug in Juniper - handling of user-defined functions.
        ])
        MINLPTests.test_nlp_mi_cvx(solver)
    end
end


text:
module GraphGenerators
export path_graph, cycle_graph, complete_graph, star_graph, wheel_graph

using ...CSetDataStructures, ..BasicGraphs
using ...CSetDataStructures: hom

""" Path graph on ``n`` vertices.
"""
function path_graph(::Type{T}, n::Int; V=(;), E=(;)) where T <: ACSet
  g = T()
  add_vertices!(g, n; V...)
  add_edges!(g, 1:(n-1), 2:n; E...)
  g
end

""" Cycle graph on ``n`` vertices.

When ``n = 1``, this is a loop graph.
"""
function cycle_graph(::Type{T}, n::Int; V=(;), E=(;)) where T <: ACSet
  g = T()
  add_vertices!(g, n; V...)
  add_edges!(g, 1:n, circshift(1:n, -1); E...)
  g
end

""" Complete graph on ``n`` vertices.
"""
function complete_graph(::Type{T}, n::Int; V=(;)) where T <: ACSet
  g = T()
  add_vertices!(g, n; V...)
  for i in vertices(g), j in vertices(g)
    if i != j && (is_directed(T) || i < j)
      add_edge!(g, i, j)
    end
  end
  g
end

""" Star graph on ``n`` vertices.

In the directed case, the edges point outward.
"""
function star_graph(::Type{T}, n::Int) where T <: ACSet
  g = T()
  add_vertices!(g, n)
  add_edges!(g, fill(n,n-1), 1:(n-1))
  g
end

""" Wheel graph on ``n`` vertices.

In the directed case, the outer cycle is directed and the spokes point outward.
"""
function wheel_graph(::Type{T}, n::Int) where T <: ACSet
  g = cycle_graph(T, n-1)
  add_vertex!(g)
  add_edges!(g, fill(n,n-1), 1:(n-1))
  g
end

# Should this be exported from `BasicGraphs`?
@generated function is_directed(::Type{T}) where {S, Ts, T<:StructACSet{S,Ts}}
  !(:inv ∈ hom(S))
end

end


text:
# This is a rewrite of SimJulia.jl example of
# Ross S. (2012): Simulation, 5th edition, Section 7.7, p. 124-126

using EventSimulation
using Distributions
using Random

const RUNS = 5
const N = 10
const S = 3
const LAMBDA = 100
const MU = 1
const SEED = 150

const F = Exponential(LAMBDA)
const G = Exponential(MU)

mutable struct RepairState <: AbstractState
    spare::Int
    repair::Int
end

function breakdown(s)
    if s.state.spare == 0
        println("At time $(s.now): No more spares!")
        terminate!(s)
    else
        s.state.spare -= 1
        register!(s, breakdown, rand(F))
        s.state.repair += 1
        s.state.repair == 1 && register!(s, repaired, rand(G))
    end
end

function repaired(s)
    s.state.spare += 1
    s.state.repair -= 1
    s.state.repair > 0 && register!(s, repaired, rand(G))
end

function exec()
    state = RepairState(S, 0)
    s = Scheduler(state)
    for i in 1:N
        register!(s, breakdown, rand(F))
    end
    go!(s)
    s.now
end

Random.seed!(SEED)
results = Float64[]
for i in 1:RUNS push!(results, exec()) end
println("Average crash time: ", sum(results)/RUNS)