text:
getshape(::FiniteElement{dim,shape}) where {dim,shape} = shape

function spatial_nodal_coordinate(mesh::AbstractPolytopalMesh{dim,T}, ci::Int,fe::FiniteElement{dim},n_point::Int) where {dim,T}
    x = getverticescoords(mesh,getcell(mesh,ci))
    n_base_funcs = getngeombasefunctions(fe)
    @assert length(x) == n_base_funcs
    vec = zero(Tensors.Vec{dim,T})
    @inbounds for i in 1:n_base_funcs
        vec += fe.M[i, n_point] * x[i]
    end
    return vec
end

"""
gradient_value(ip::FiniteElement{dim,shape<:Shape,order}, k::Int, ξ::Vec{dim,T}) where {dim,shape,order, T}
Compute value of Finite Element basis `j` derivative at point ξ
on the reference shape `shape`
"""
function gradient_value(ip::FiniteElement{dim}, k::Int, ξ::Tensors.Vec{dim,T}) where {dim, T}
    Tensors.gradient(ξ -> value(ip, k, ξ), ξ)
end

max_stars_repo_name: Paulms/jFEMToools

max_issues_repo_name: Paulms/jFEMToools

max_forks_repo_name: Paulms/jFEMToools


text:
include("src/PayGroups.jl")
using .PayGroups
igroupay()

ext: jl


text:
function triple(::Type{T}, x::BigFloat) where {T<:Base.IEEEFloat}
    prec = precision(BigFloat)
    setprecision(BigFloat, 768)
    hi = T(x)
    md = T(x - hi)
    lo = T(x - hi - md)
    setprecision(BigFloat, prec)
    return hi,md,lo
end

function triple(::Type{T}, x::String) where {T<:Base.IEEEFloat}
    prec = precision(BigFloat)
    setprecision(BigFloat, 768)
    z = parse(BigFloat, x)
    hi = T(z)
    md = T(z - hi)
    lo = T(z - hi - md)
    setprecision(BigFloat, prec)
    return hi,md,lo
end

triple64(x::BigFloat) = triple(Float64, x)
triple32(x::BigFloat) = triple(Float32, x)
triple16(x::BigFloat) = triple(Float16, x)

triple64(x::String) = triple(Float64, x)
triple32(x::String) = triple(Float32, x)
triple16(x::String) = triple(Float16, x)

function triple_inv(::Type{T}, x::BigFloat) where {T<:Base.IEEEFloat}
    prec = precision(BigFloat)
    setprecision(BigFloat, 768)
    x = inv(x)
    hi = T(x)
    md = T(x - hi)
    lo = T(x - hi - md)
    setprecision(BigFloat, prec)
    return hi,md,lo
end

triple64inv(x::BigFloat) = triple_inv(Float64, x)
triple32inv(x::BigFloat) = triple_inv(Float32, x)
triple16inv(x::BigFloat) = triple_inv(Float16, x)


text:
@testset "xavier_init" begin
    @test_nowarn a = xavier_init([100,10], Float32)
end

@testset "load_op" begin
    @test begin
        # ADCME.install_custom_op_dependency()
        ADCME.compile("SparseSolver")
        # somehow we cannot first call `load_op_and_grad` and then call `load_op`
        load_op("$(@__DIR__)/../deps/CustomOps/SparseSolver/build/libSparseSolver", "sparse_solver")
        load_op_and_grad("$(@__DIR__)/../deps/CustomOps/SparseSolver/build/libSparseSolver", "sparse_solver")
        true
    end
end

@testset "install" begin
    install("OTNetwork", force=true)
end

@testset "ae" begin
    n = ae_num([2,20,20,20,3])
    @test n==length(ae_init([2,20,20,20,3]))
end

@testset "register" begin
    forward = x->log(1+exp(x))
    backward = (dy, y, x)->dy*(1-1/(1+y))
    f = register(forward, backward)
    x_ = rand(10)
    x = constant(x_)
    @test run(sess, f(x)) ≈ @. log(1+exp(x_))
    @test_nowarn run(sess,gradients(sum(f(x)), x))

    function forward2(x, θ)
        f = (θ, y)->(y^3+1.0-θ -x, spdiag(3y^2))
        nr = newton_raphson(f, constant(ones(length(x))), θ, options=Dict("verbose"=>true))
        y = nr.x
        return y
    end

    function backward2(dy, y, x, θ)
        e = 1/3*(θ+x-1)^(-2/3)
        e, e
    end

    f2 = register(forward2, backward2)

    θ = constant(ones(5))
    x = constant(8ones(5))
    y = f2(x, θ)
    @test run(sess, y)≈ones(5)*2

    @test run(sess, gradients(sum(y), θ))≈0.08333333333333333*ones(5)
end


text:
#this file is to create omcis data for missing pattern (ii)
# for each omics feature, some random individuals have missing omics data
using Distributions, DelimitedFiles,Random,CSV,DataFrames
Random.seed!(123)

cd("/Users/tianjing/Box/Omicsdata/nonlinear/varmi_2/")

pct_all=[0.1,0.3,0.5,0.7,0.9]
repAll=20

omicsfile="data1/omics.csv"
omics_all= CSV.read(omicsfile,DataFrame,delim = ',',header=true,missingstrings=["NA"])
for i in 2:size(omics_all,2) #skip ID
    omics_all[!,i]= convert(Vector{Union{Missing,Float64}}, omics_all[!,i]);
end

@show size(omics_all)
@show omics_all[1:5,1:5];
@show omics_all[end-4:end,1:5];

Random.seed!(123)
for i in 1:length(pct_all)
    pct=pct_all[i]
    num_missing_per_gene=Int(floor(955*pct))  #1-955 are training individuals
    for rep in 1:repAll
        omics=copy(omics_all)
        for i in 2:1201
            omics[sample(1:955,num_missing_per_gene,replace=false),i].=missing;
        end
        @show [count(ismissing,col) for col in eachcol(omics)]
        CSV.write("data1/partial_missing_data/omics_partialmissing.pct$pct.rep$rep.csv",omics,missingstring="NA")
    end
end

max_stars_repo_name: zhaotianjing/NN-LMM

max_issues_repo_name: zhaotianjing/NN-LMM

max_forks_repo_name: zhaotianjing/NN-LMM


text:
"""
General state types
"""

abstract type AbstractState{T} end
abstract type AbstractAbsoluteState{T} <: AbstractState{T} end
abstract type AbstractUncertainState{T} <: AbstractState{T} end

const Covariance{T} = Array{T,2} where T

"""
    DiscreteState(x::Vector[, t::Int64])
    DiscreteState(x::AbstractFloat[, t::Int64])

Discrete-time state. Contains a state vector and an Integer-valued time step. If
not provided, the default value of `t` is 0.
"""
mutable struct DiscreteState{T<:AbstractFloat} <: AbstractAbsoluteState{T}
    x::Vector{T}
    t::Int64

    DiscreteState(x::Vector{T}, t::Int64) where {T} = new{T}(deepcopy(x), t)
    DiscreteState(x::Vector{T}) where {T} = new{T}(deepcopy(x), 0)
end
DiscreteState(x::AbstractFloat) = DiscreteState([x])
DiscreteState(x::AbstractFloat, t::Int64) = DiscreteState([x], t)

"""
    UncertainDiscreteState(x::Vector, P::Matrix[, t::Int64])
    UncertainDiscreteState(x::AbstractFloat, P::AbstractFloat[, t::Int64])

Uncertain discrete-time state. Contains a state vector, a covariance matrix
representing the uncertainty in that state, and an Integer-valued time step. If
not provided, the default value of `t` is 0.

Uncertain states are Gaussian distributioned random varaibles with mean \$x\$
and covariance \$P\$.
"""
mutable struct UncertainDiscreteState{T<:AbstractFloat} <:
    AbstractUncertainState{T}

    x::Vector{T}
    P::Covariance{T}
    t::Int64

    function UncertainDiscreteState(x::Vector{T}, P::Covariance{T}, t) where T
            if length(x) != size(P)[1] || length(x) != size(P)[2]
                throw(DimensionMismatch(
                    "Dimensions of covariance and state do not match."))
            end
            new{T}(deepcopy(x), deepcopy(P), t)
    end
    function UncertainDiscreteState(x::Vector{T}, P::Covariance{T}) where T
            if length(x) != size(P)[1] || length(x) != size(P)[2]
                throw(DimensionMismatch(
                    "Dimensions of covariance and state do not match."))
            end
            new{T}(deepcopy(x), deepcopy(P), 0)
    end
end
UncertainDiscreteState{T<:AbstractFloat}(x::T, P::T) =
    UncertainDiscreteState([x], reshape([P],1,1))
UncertainDiscreteState{T<:AbstractFloat}(x::T, P::T, t) =
    UncertainDiscreteState([x], reshape([P],1,1), t)

"""
Union of Discrete state types (absolute and uncertain).
"""
const UnionDiscrete{T} = Union{DiscreteState{T}, UncertainDiscreteState{T}} where T

"""
    ContinuousState(x::Vector[, t::AbstractFloat])
    ContinuousState(x::AbstractFloat[, t::AbstractFloat])

Continuous-time state. Contains a state vector and a floating-point time step.
If not provided, the default value of `t` is 0.0.
"""
mutable struct ContinuousState{T<:AbstractFloat} <: AbstractAbsoluteState{T}
    x::Vector{T}
    t::T

    ContinuousState(x::Vector{T}, t::T) where {T} = new{T}(deepcopy(x), t)
    ContinuousState(x::Vector{T}) where {T} = new{T}(deepcopy(x), 0.0)
end
ContinuousState(x::AbstractFloat) = ContinuousState([x])
ContinuousState{T<:AbstractFloat}(x::T, t::T) = ContinuousState([x], t)

"""
    UncertainContinuousState(x::Vector, P::Matrix[, t::AbstractFloat])
    UncertainContinuousState(x::AbstractFloat, P::AbstractFloat[, t::AbstractFloat])

Uncertain continuous-time state. Contains a state vector, a covariance matrix
representing the uncertainty in that state, and a floating-point time step. If
not provided, the default value for `t` is 0.0.

Uncertain states are Gaussian distributioned random varaibles with mean \$x\$
and covariance \$P\$.
"""
mutable struct UncertainContinuousState{T<:AbstractFloat} <:
    AbstractUncertainState{T}

    x::Vector{T}
    P::Covariance{T}
    t::T

    function UncertainContinuousState(x::Vector{T}, P::Covariance{T},
        t::T) where T
            if length(x) != size(P)[1] || length(x) != size(P)[2]
                throw(DimensionMismatch(
                    "Dimensions of covariance and state do not match."))
            end
            new{T}(deepcopy(x), deepcopy(P), t)
    end
    function UncertainContinuousState(x::Vector{T}, P::Covariance{T}) where T
            if length(x) != size(P)[1] || length(x) != size(P)[2]
                throw(DimensionMismatch(
                    "Dimensions of covariance and state do not match."))
            end
            new{T}(deepcopy(x), deepcopy(P), 0.0)
    end
end
UncertainContinuousState{T<:AbstractFloat}(x::T, P::T) =
    UncertainContinuousState([x], reshape([P],1,1))
UncertainContinuousState{T<:AbstractFloat}(x::T, P::T, t::T) =
    UncertainContinuousState([x], reshape([P],1,1), t)

"""
Union of Continuous state types (absolute and uncertain).
"""
const UnionContinuous{T} = Union{ContinuousState{T}, UncertainContinuousState{T}} where T

"""
    make_uncertain(state::AbstractUncertainState[, P::Matrix])

Create an uncertain copy of the input absolute state. If not provided, the
covariance matrix, `P`, is all zero.
"""
make_uncertain{T}(state::DiscreteState{T}) =
    UncertainDiscreteState(
        state.x, zeros(T, length(state.x), length(state.x)), state.t)
make_uncertain{T}(state::DiscreteState{T}, P::Covariance{T}) =
    UncertainDiscreteState(state.x, P, state.t)
make_uncertain{T}(state::ContinuousState{T}) =
    UncertainContinuousState(
        state.x, zeros(T, length(state.x), length(state.x)), state.t)
make_uncertain{T}(state::ContinuousState{T}, P::Covariance{T}) =
    UncertainContinuousState(state.x, P, state.t)

"""
    make_absolute(state::AbstractUncertainState)

Create an absolute copy of the input uncertain state.
"""
make_absolute(state::UncertainDiscreteState) = DiscreteState(state.x, state.t)
make_absolute(state::UncertainContinuousState) =
    ContinuousState(state.x, state.t)

"""
"""
absolute_type{T}(state::UncertainDiscreteState{T}) = DiscreteState{T}
absolute_type{T}(state::UncertainContinuousState{T}) = ContinuousState{T}

"""
"""
uncertain_type{T}(state::DiscreteState{T}) = UncertainDiscreteState{T}
uncertain_type{T}(state::ContinuousState{T}) = UncertainContinuousState{T}

"""
    sample(state::AbstractUncertainState)

Sample uncertain state. Returns absolute state of the corresponding time
format (i.e. UncertainDiscreteState => DiscreteState).
"""
function sample{T}(state::UncertainDiscreteState{T})
    if all(state.P .== 0)
        sampled_state = state.x
    else
        sampled_state = state.x+chol(Hermitian(state.P))'*randn(T, length(state.x))
    end
    return DiscreteState(sampled_state, state.t)
end
function sample{T}(state::UncertainContinuousState{T})
    if all(state.P .== 0)
        sampled_state = state.x
    else
        sampled_state = state.x+chol(Hermitian(state.P))'*randn(T, length(state.x))
    end
    return ContinuousState(sampled_state, state.t)
end

"""
    distance(state1::AbstractState, state2::AbstractState)

Computes the Euclidean distance between two states. Ignores uncertainty in
either of the states.

\$d(x,y) = \\left\\|{x-y}\\right\\|\$
"""
function distance{T}(state1::AbstractState{T}, state2::AbstractState{T})
    return sqrt(sum((state1.x .- state2.x).^2))
end

"""
    mahalanobis(x::AbstractAbsoluteState, y::AbstractUncertainState)
    mahalanobis(x::AbstractAbsoluteState, y::AbstractAbsoluteState, P::Matrix)

Compute the Mahalanobis distance between a point and a distribution or two
points in the same distribution. If `P` is the covariance of the distribution,

\$d(x,y) = \\sqrt{(x-y)^T P^{-1} (x-y)}\$
"""
function mahalanobis{T}(x::AbstractAbsoluteState{T},
                        y::AbstractUncertainState{T})
    diff_vec = x.x .- y.x
    return sqrt(diff_vec'*inv(y.P)*diff_vec)
end
function mahalanobis{T}(x::AbstractAbsoluteState{T},
                        y::AbstractAbsoluteState{T}, P::Array{T,2})
    diff_vec = x.x .- y.x
    return sqrt(diff_vec'*inv(P)*diff_vec)
end

"""
    block_diagonal(A1::Matrix, A2::Matrix)

Simple blkdiag alternative for two non-sparse matrices.
"""
function block_diagonal(A1, A2)
    m1, n1 = size(A1)
    m2, n2 = size(A2)
    A_out = zeros(m1+m2, n1+n2)
    A_out[1:m1, 1:n1] .= A1
    A_out[m1+1:end, n1+1:end] .= A2
    return A_out
end

"""
    expand_state(state::AbstractState, x[, P])

Increase the dimension of an AbstractState type by concatenating a new states.
For uncertain state types, a covariance is required. The covariances of the
existing and new states are joined in a block diagonal form. Works for both
scalar and array input.
"""
function expand_state(state::AbstractAbsoluteState, x)
    augmented_state = deepcopy(state)
    augmented_state.x = vcat(state.x, x)
    return augmented_state
end
function expand_state(state::AbstractUncertainState, x::Real, P::Real)
    return expand_state(state, [x], fill(P, (1,1)))
end
function expand_state(state::AbstractUncertainState, x::Vector, P::Matrix)
    n = length(x)
    if size(P) != (n,n)
        throw(DimensionMismatch("Dimensions of x and P are inconsistent."))
    end
    augmented_state = deepcopy(state)
    augmented_state.x = vcat(state.x, x)
    augmented_state.P = block_diagonal(state.P, P)
    return augmented_state
end

"""
    reduce_state(state::AbstractState, ndim::Integer)

Reduce the dimension of an AbstractState to a specified number of dimensions.
This method discards state data beyond the desired number of dimensions and,
for uncertain state types, discards the corresponding covariance data.
"""
function reduce_state(state::AbstractAbsoluteState, ndim::Integer)
    reduced_state = deepcopy(state)
    reduced_state.x = state.x[1:ndim]
    return reduced_state
end
function reduce_state(state::AbstractUncertainState, ndim::Integer)
    reduced_state = deepcopy(state)
    reduced_state.x = state.x[1:ndim]
    reduced_state.P = state.P[1:ndim, 1:ndim]
    return reduced_state
end

"""
State addition.
"""
function Base.:+(state::AbstractState, a)
    out = deepcopy(state)
    out.x .+= a
    return out
end
function Base.:+(a, state::AbstractState)
    out = deepcopy(state)
    out.x .+= a
    return out
end
function Base.broadcast!(::typeof(+), out::AbstractState, state::AbstractState,
                         a)
    out.x .= state.x .+ a
    return nothing
end
function Base.broadcast!(::typeof(+), out::AbstractState, a,
                         state::AbstractState)
    out.x .= state.x .+ a
    return nothing
end

"""
State subtraction.
"""
function Base.:-(state::AbstractState, a)
    out = deepcopy(state)
    out.x .-= a
    return out
end
function Base.:-(a, state::AbstractState)
    out = deepcopy(state)
    out.x .-= a
    return out
end
function Base.broadcast!(::typeof(-), out::AbstractState, state::AbstractState,
                         a)
    out.x .= state.x .- a
    return nothing
end
function Base.broadcast!(::typeof(-), out::AbstractState, a,
                         state::AbstractState)
    out.x .= state.x .- a
    return nothing
end

"""
State equality check.
"""
function Base.:(==){S<:AbstractAbsoluteState}(state1::S, state2::S)
    if state1.x == state2.x && state1.t == state2.t
        return true
    else
        return false
    end
end
function Base.:(==){S<:AbstractUncertainState}(state1::S, state2::S)
    if state1.x == state2.x && state1.P == state2.P && state1.t == state2.t
        return true
    else
        return false
    end
end
Base.:(==)(state1::AbstractState, state2::AbstractState) = false

"""
State multiplication.
"""
function Base.:*(state::AbstractState, a::Number)
    out = deepcopy(state)
    out .*= a
    return out
end
function Base.:*(a::Number, state::AbstractState)
    out = deepcopy(state)
    out .*= a
    return out
end
function Base.broadcast!(::typeof(*), out::AbstractAbsoluteState,
                         state::AbstractState, a::Number)
    out.x .= state.x .* a
    return nothing
end
function Base.broadcast!(::typeof(*), out::AbstractAbsoluteState, a::Number,
                         state::AbstractState)
    out.x .= state.x .* a
    return nothing
end
function Base.broadcast!(::typeof(*), out::AbstractUncertainState, a::Number,
                         state::AbstractUncertainState)
    out.x .= state.x .* a
    out.P .= state.P .* a^2
    return nothing
end
function Base.broadcast!(::typeof(*), out::AbstractUncertainState,
                         state::AbstractUncertainState, a::Number)
    out.x .= state.x .* a
    out.P .= state.P .* a^2
    return nothing
end

function Base.:*(A::Matrix, state::AbstractState)
    out = deepcopy(state)
    out .= A .* out
    return out
end
function Base.broadcast!(::typeof(*), out::AbstractAbsoluteState, A::Matrix,
                         state::AbstractState)
    out.x .= A * state.x
    return nothing
end
function Base.broadcast!(::typeof(*), out::AbstractUncertainState, A::Matrix,
                         state::AbstractUncertainState)
    out.x .= A * state.x
    out.P .= A * state.P * A'
    return nothing
end

"""
State assignment operations.
"""
function Base.broadcast!(::typeof(identity), out::DiscreteState,
                            state::DiscreteState)
    out.x .= state.x
    out.t = state.t
    return nothing
end
function Base.broadcast!(::typeof(identity), out::ContinuousState,
                            state::ContinuousState)
    out.x .= state.x
    out.t = state.t
    return nothing
end
function Base.broadcast!(::typeof(identity), out::UncertainDiscreteState,
                            state::UncertainDiscreteState)
    out.x .= state.x
    out.P .= state.P
    out.t = state.t
    return nothing
end
function Base.broadcast!(::typeof(identity), out::UncertainContinuousState,
                            state::UncertainContinuousState)
    out.x .= state.x
    out.P .= state.P
    out.t = state.t
    return nothing
end


text:
module TestDimensions

using Mimi
using Test

import Mimi:
    compdef, AbstractDimension, RangeDimension, Dimension, key_type, first_period, last_period,
    ComponentReference, ComponentPath, ComponentDef, time_labels

##
## Constants
##

dim_varargs = Dimension(:foo, :bar, :baz)   # varargs
dim_vec = Dimension([:foo, :bar, :baz]) # Vector
dim_range = Dimension(2010:2100)     # AbstractRange
rangedim = RangeDimension(2010:2100) # RangeDimension type
dim_vals = Dimension(4) # Same as 1:4

##
## Test a Bunch of Small Functionalities and Helpers
##

@test key_type(dim_varargs) == Symbol
@test key_type(dim_vec) == Symbol
@test key_type(dim_range) == Int
@test key_type(dim_vals) == Int

@test length(dim_varargs) == 3
@test length(dim_vec) == 3
@test length(dim_vals) == 4
@test length(dim_range) == 91
@test length(rangedim) == 91

# Test iteration
@test [x.first for x in dim_varargs] == collect(keys(dim_varargs))
@test [x.first for x in dim_range]   == collect(keys(dim_range))

@test lastindex(dim_varargs) == :baz
@test lastindex(dim_vec) == :baz
@test lastindex(dim_range) == 2100
@test lastindex(dim_vals) == 4

@test Base.keys(rangedim) == [2010:2100...]
@test Base.values(rangedim) == [1:91...]

for i = 1:length(dim_range)
    @test getindex(dim_range, 2010 + i - 1) == i
end
@test dim_range[:] == [1:91...]
@test dim_varargs[:bar] == 2
@test dim_varargs[:] == [1,2,3]

# @test rangedim[2011] == 2 # TODO: this errors..
@test get(dim_varargs, :bar, 999) == 2
@test get(dim_varargs, :new, 4) == 4 #adds a key/value pair
@test get(rangedim, 2010, 1) == 1
# @test get(rangedim, 2101, 92) == 92 # TODO:  this errors ...

#test iteratable
dim_vals2 = Dimension(2:2:8)
intkeys = [2,4,6,8]
# Work around new global scoping rules
# Without the `let`, `index` is unknown in the for loop
let index = 1
    for pair in dim_vals2
        @test dim_vals2[intkeys[index]] == index
        index += 1
    end
end

rangedim2 = RangeDimension(2:2:8)
# Work around new global scoping rules
# Without the `let`, `index` is unknown in the for loop
let index = 1
    for pair in rangedim2   # uses iterate()
        @test get(rangedim2, Base.keys(rangedim2)[index]) == index
        index += 1
    end
end

@test getindex(dim_varargs, :bar) == 2
@test getindex(dim_varargs, :) == [1,2,3]

# Test that d.time returns AbstracTimesteps that can be used as indexes
@defcomp bar begin
    v1 = Variable(index = [time])

    function run_timestep(p, v, d, t)
        if is_first(t)
            for i in d.time
                v.v1[i] = gettime(i)
            end
        end
    end
end

fixed_years = 2000:2010
variable_years = [2000, 2005, 2020, 2050, 2100]

m = Model()
set_dimension!(m, :time, fixed_years)
add_comp!(m, bar)
run(m)
@test m[:bar, :v1] == fixed_years

m = Model()
set_dimension!(m, :time, variable_years)
add_comp!(m, bar)
run(m)
@test m[:bar, :v1] == variable_years

# test variable_dimensions function
comp_path = ComponentPath((:bar,))
dims = [:time]
@test variable_dimensions(m, comp_path, :v1) == dims
@test variable_dimensions(m, :bar, :v1) == dims
@test variable_dimensions(m, (:bar,), :v1) == dims

##
## Test time dimension (and resetting it!)
##

@defcomp foo2 begin
    x = Parameter(index=[time])
    y = Variable(index=[4])
end

# build model, set dims, and add components
m = Model()
@test_throws ErrorException add_comp!(m, foo2) # cannot add a component before time dimension is set
set_dimension!(m, :time, 2000:2100)

@test_throws ErrorException add_comp!(m, foo2; first = 2000, last = 2105)   # 2105 cannot be found in the model's time dimension
@test_throws ErrorException add_comp!(m, foo2; first = 1950, last = 2100)   # 1950 cannot be found in the model's time dimension

foo2_ref1 = add_comp!(m, foo2)
foo2_ref2 = ComponentReference(m, :foo2)
@test foo2_ref1 === foo2_ref2
my_foo2 = compdef(foo2_ref1)

@test first_period(m.md) == first_period(m.md.namespace[:foo2]) == 2000
@test last_period(m.md) == last_period(m.md.namespace[:foo2]) == 2100

# Set Parameters
original_x_vals = collect(2000:2100)
@test_throws ErrorException set_param!(m, :foo2, :x, 1990:2200) # too long
@test_throws ErrorException set_param!(m, :foo2, :x, 2005:2095) # too short
set_param!(m, :foo2, :x, original_x_vals)

run(m)

# Reset Dimension
@test_throws ErrorException set_dimension!(m, :time, 2050:2100) # can't move time forward
@test_throws ErrorException set_dimension!(m, :time, 2105:2200) # can't move new first past old last
set_dimension!(m, :time, 1990:2050)

@test first_period(m.md) == 1990
@test last_period(m.md) == 2050
@test first_period(m.md.namespace[:foo2]) == 2000 # no change
@test last_period(m.md.namespace[:foo2]) == 2050 # trimmed with model

# check that parameters were padded properly
new_x_vals = m.md.external_params[:x].values.data
@test length(new_x_vals) == length(time_labels(m))
@test new_x_vals[11:end] == original_x_vals[1:51]
@test all(ismissing, new_x_vals[1:10])

run(m) # should still run because parameters were adjusted under the hood

# reset again with late end
set_dimension!(m, :time, 1990:2200)
new_x_vals = m.md.external_params[:x].values.data
@test length(new_x_vals) == length(time_labels(m))
@test all(ismissing, new_x_vals[1:10])
@test new_x_vals[11:61] == original_x_vals[1:51]
@test all(ismissing, new_x_vals[62:end])

run(m)
end #module


text:
function parse(instance::CompilerInstance)
    begin_diag(instance)
    try
        InitializeBuiltins(getPreprocessor(instance))
        ParseAST(getSema(instance))
    finally
        end_diag(instance)
    end
    return true
end

function parse(instance::CompilerInstance, parser::Parser, codegen::CodeGenerator)
    diag_csr = getDiagnosticClient(instance)

    preprocessor = getPreprocessor(instance)
    EnterMainSourceFile(preprocessor)
    InitializeBuiltins(preprocessor)

    begin_source_file(diag_csr, getLangOpts(instance), preprocessor)

    initialize(parser)

    sema = getSema(instance)
    ast_ctx = getASTContext(instance)

    if tryParseAndSkipInvalidOrParsedDecl(parser, codegen)
        processWeakTopLevelDecls(sema, codegen)
        HandleTranslationUnit(codegen, ast_ctx)
        end_source_file(diag_csr)
        return true
    else
        end_source_file(diag_csr)
        return false
    end
end

function parse_cxx_scope_spec(p::Parser, spec::CXXScopeSpec)
    sema = getSema(p)
    pp = getPreprocessor(p)
    tok = getCurToken(p)
    is_incremental(pp) && is_eof(tok) && ConsumeToken(p)
    while !is_eof(tok)
        if is_coloncolon(tok)  # `::foo::bar::baz`
            TryAnnotateCXXScopeToken(p) && error("failed to annotate token.")
        elseif is_identifier(tok)  # `foo::bar::baz`
            TryAnnotateCXXScopeToken(p) && error("failed to annotate token.")
            if !is_identifier(tok)
                continue
            else
                ConsumeToken(p)
            end
        elseif is_annot_cxxscope(tok) || is_annot_typename(tok)
            restore_nns_annotation(sema, tok, spec)
            break
        elseif is_annot_template_id(tok)
            TryAnnotateTypeOrScopeTokenAfterScopeSpec(p, spec)
            break
        else  # skip other cases which we may add support in the future
            ConsumeAnyToken(p)
        end
    end
    # parse to the end of file, ignore everything
    while !is_eof(tok)
        ConsumeAnyToken(p)
    end
    return nothing
end

function parse_cxx_scope_spec(p::Parser, fid::FileID, spec::CXXScopeSpec)
    pp = getPreprocessor(p)
    EnterSourceFile(pp, fid)
    try
        parse_cxx_scope_spec(p, spec)
    finally
        EndSourceFile(pp)
    end
    return nothing
end

function parse_cxx_scope_spec(ci::CompilerInstance, p::Parser, str::String, spec::CXXScopeSpec)
    src_mgr = getSourceManager(ci)
    begin_diag(ci)
    fid = FileID(src_mgr, get_buffer(str))
    try
        parse_cxx_scope_spec(p, fid, spec)
    finally
        dispose(fid)
        end_diag(ci)
    end
    return nothing
end


text:
###
### Conversion & Promotion
###
###
### Conversion methods for LongSequences.
###
### This file is a part of BioJulia.
### License is MIT: https://github.com/BioJulia/BioSequences.jl/blob/master/LICENSE.md

###
### Promotion
###
for alph in (DNAAlphabet, RNAAlphabet)
    @eval function Base.promote_rule(::Type{LongSequence{A}}, ::Type{LongSequence{B}}) where {A<:$alph,B<:$alph}
        return LongSequence{promote_rule(A, B)}
    end
end

###
### Conversion
###

# Create a 4 bit DNA/RNA sequence from a 2 bit DNA/RNA sequence, and vice-versa.
for (alpha, alphb) in [(DNAAlphabet{4}, DNAAlphabet{2}), # DNA to DNA
                       (DNAAlphabet{2}, DNAAlphabet{4}),
                       (RNAAlphabet{4}, RNAAlphabet{2}), # RNA to RNA
                       (RNAAlphabet{2}, RNAAlphabet{4}),
                       (DNAAlphabet{2}, RNAAlphabet{4}), # DNA to RNA
                       (DNAAlphabet{4}, RNAAlphabet{2}),
                       (DNAAlphabet{2}, RNAAlphabet{2}),
                       (DNAAlphabet{4}, RNAAlphabet{4}),
                       (RNAAlphabet{4}, DNAAlphabet{2}), # RNA to DNA
                       (RNAAlphabet{2}, DNAAlphabet{4}),
                       (RNAAlphabet{2}, DNAAlphabet{2}),
                       (RNAAlphabet{4}, DNAAlphabet{4})]

    @eval function Base.convert(::Type{LongSequence{$alpha}}, seq::LongSequence{$alphb})
        return LongSequence{$alpha}(seq)
    end
end

# Convert from a LongSequence to to a DNA or RNA vector
function Base.convert(::Type{LongSequence{A}}, seq::Vector) where A<:Alphabet
    return LongSequence{A}(seq)
end

Base.parse(::Type{LongSequence{A}}, seq::AbstractString) where A = LongSequence{A}(seq)


text:
module GenTransportMCMC

include("transport/transport.jl")
include("inference/inference.jl")

end

max_stars_repo_name: mvarble/gen-transport-mcmc

max_issues_repo_name: mvarble/gen-transport-mcmc

max_forks_repo_name: mvarble/gen-transport-mcmc


text:
# UnitCommitment.jl: Optimization Package for Security-Constrained Unit Commitment
# Copyright (C) 2020, UChicago Argonne, LLC. All rights reserved.
# Released under the modified BSD license. See COPYING.md for more details.

function _add_unit!(model::JuMP.Model, g::Unit, formulation::Formulation)
    if !all(g.must_run) && any(g.must_run)
        error("Partially must-run units are not currently supported")
    end
    if g.initial_power === nothing || g.initial_status === nothing
        error("Initial conditions for $(g.name) must be provided")
    end

    # Variables
    _add_production_vars!(model, g, formulation.prod_vars)
    _add_reserve_vars!(model, g)
    _add_startup_shutdown_vars!(model, g)
    _add_status_vars!(model, g, formulation.status_vars)

    # Constraints and objective function
    _add_min_uptime_downtime_eqs!(model, g)
    _add_net_injection_eqs!(model, g)
    _add_production_limit_eqs!(model, g, formulation.prod_vars)
    _add_production_piecewise_linear_eqs!(
        model,
        g,
        formulation.prod_vars,
        formulation.pwl_costs,
        formulation.status_vars,
    )
    _add_ramp_eqs!(
        model,
        g,
        formulation.prod_vars,
        formulation.ramping,
        formulation.status_vars,
    )
    _add_startup_cost_eqs!(model, g, formulation.startup_costs)
    _add_startup_shutdown_limit_eqs!(model, g)
    _add_status_eqs!(model, g, formulation.status_vars)
    return
end

_is_initially_on(g::Unit)::Float64 = (g.initial_status > 0 ? 1.0 : 0.0)

function _add_reserve_vars!(model::JuMP.Model, g::Unit)::Nothing
    reserve = _init(model, :reserve)
    reserve_shortfall = _init(model, :reserve_shortfall)
    for t in 1:model[:instance].time
        if g.provides_spinning_reserves[t]
            reserve[g.name, t] = @variable(model, lower_bound = 0)
        else
            reserve[g.name, t] = 0.0
        end
        reserve_shortfall[t] =
            (model[:instance].shortfall_penalty[t] >= 0) ?
            @variable(model, lower_bound = 0) : 0.0
    end
    return
end

function _add_reserve_eqs!(model::JuMP.Model, g::Unit)::Nothing
    reserve = model[:reserve]
    for t in 1:model[:instance].time
        add_to_expression!(expr_reserve[g.bus.name, t], reserve[g.name, t], 1.0)
    end
    return
end

function _add_startup_shutdown_vars!(model::JuMP.Model, g::Unit)::Nothing
    startup = _init(model, :startup)
    for t in 1:model[:instance].time
        for s in 1:length(g.startup_categories)
            startup[g.name, t, s] = @variable(model, binary = true)
        end
    end
    return
end

function _add_startup_shutdown_limit_eqs!(model::JuMP.Model, g::Unit)::Nothing
    eq_shutdown_limit = _init(model, :eq_shutdown_limit)
    eq_startup_limit = _init(model, :eq_startup_limit)
    is_on = model[:is_on]
    prod_above = model[:prod_above]
    reserve = model[:reserve]
    switch_off = model[:switch_off]
    switch_on = model[:switch_on]
    T = model[:instance].time
    for t in 1:T
        # Startup limit
        eq_startup_limit[g.name, t] = @constraint(
            model,
            prod_above[g.name, t] + reserve[g.name, t] <=
            (g.max_power[t] - g.min_power[t]) * is_on[g.name, t] -
            max(0, g.max_power[t] - g.startup_limit) * switch_on[g.name, t]
        )
        # Shutdown limit
        if g.initial_power > g.shutdown_limit
            eq_shutdown_limit[g.name, 0] =
                @constraint(model, switch_off[g.name, 1] <= 0)
        end
        if t < T
            eq_shutdown_limit[g.name, t] = @constraint(
                model,
                prod_above[g.name, t] <=
                (g.max_power[t] - g.min_power[t]) * is_on[g.name, t] -
                max(0, g.max_power[t] - g.shutdown_limit) *
                switch_off[g.name, t+1]
            )
        end
    end
    return
end

function _add_ramp_eqs!(
    model::JuMP.Model,
    g::Unit,
    formulation::RampingFormulation,
)::Nothing
    prod_above = model[:prod_above]
    reserve = model[:reserve]
    eq_ramp_up = _init(model, :eq_ramp_up)
    eq_ramp_down = _init(model, :eq_ramp_down)
    for t in 1:model[:instance].time
        # Ramp up limit
        if t == 1
            if _is_initially_on(g) == 1
                eq_ramp_up[g.name, t] = @constraint(
                    model,
                    prod_above[g.name, t] + reserve[g.name, t] <=
                    (g.initial_power - g.min_power[t]) + g.ramp_up_limit
                )
            end
        else
            eq_ramp_up[g.name, t] = @constraint(
                model,
                prod_above[g.name, t] + reserve[g.name, t] <=
                prod_above[g.name, t-1] + g.ramp_up_limit
            )
        end

        # Ramp down limit
        if t == 1
            if _is_initially_on(g) == 1
                eq_ramp_down[g.name, t] = @constraint(
                    model,
                    prod_above[g.name, t] >=
                    (g.initial_power - g.min_power[t]) - g.ramp_down_limit
                )
            end
        else
            eq_ramp_down[g.name, t] = @constraint(
                model,
                prod_above[g.name, t] >=
                prod_above[g.name, t-1] - g.ramp_down_limit
            )
        end
    end
end

function _add_min_uptime_downtime_eqs!(model::JuMP.Model, g::Unit)::Nothing
    is_on = model[:is_on]
    switch_off = model[:switch_off]
    switch_on = model[:switch_on]
    eq_min_uptime = _init(model, :eq_min_uptime)
    eq_min_downtime = _init(model, :eq_min_downtime)
    T = model[:instance].time
    for t in 1:T
        # Minimum up-time
        eq_min_uptime[g.name, t] = @constraint(
            model,
            sum(switch_on[g.name, i] for i in (t-g.min_uptime+1):t if i >= 1) <= is_on[g.name, t]
        )
        # Minimum down-time
        eq_min_downtime[g.name, t] = @constraint(
            model,
            sum(
                switch_off[g.name, i] for i in (t-g.min_downtime+1):t if i >= 1
            ) <= 1 - is_on[g.name, t]
        )
        # Minimum up/down-time for initial periods
        if t == 1
            if g.initial_status > 0
                eq_min_uptime[g.name, 0] = @constraint(
                    model,
                    sum(
                        switch_off[g.name, i] for
                        i in 1:(g.min_uptime-g.initial_status) if i <= T
                    ) == 0
                )
            else
                eq_min_downtime[g.name, 0] = @constraint(
                    model,
                    sum(
                        switch_on[g.name, i] for
                        i in 1:(g.min_downtime+g.initial_status) if i <= T
                    ) == 0
                )
            end
        end
    end
end

function _add_net_injection_eqs!(model::JuMP.Model, g::Unit)::Nothing
    expr_net_injection = model[:expr_net_injection]
    for t in 1:model[:instance].time
        # Add to net injection expression
        add_to_expression!(
            expr_net_injection[g.bus.name, t],
            model[:prod_above][g.name, t],
            1.0,
        )
        add_to_expression!(
            expr_net_injection[g.bus.name, t],
            model[:is_on][g.name, t],
            g.min_power[t],
        )
    end
end


text:
"""
Functions for interacting with the `CondaPkg.toml` file in the current project.

This is the main `Pkg`-like API, with functions like `add`, `rm`, `status`.
"""

function cur_deps_file()
    e = Base.load_path()[1]
    if !isdir(e)
        e = dirname(e)
        @assert isdir(e)
    end
    joinpath(e, "CondaPkg.toml")
end

function read_deps(; file=cur_deps_file())
    isfile(file) ? TOML.parsefile(file) : Dict{String,Any}()
end

function write_deps(toml; file=cur_deps_file())
    if isempty(toml)
        if isfile(file)
            Base.rm(file)
        end
    else
        open(file, "w") do io
            TOML.print(io, toml)
        end
    end
    return
end

function parse_deps(toml)
    # packages
    packages = PkgSpec[]
    if haskey(toml, "deps")
        deps = _convert(Dict{String,Any}, toml["deps"])
        for (name, dep) in deps
            version = ""
            channel = ""
            if dep isa AbstractString
                version = _convert(String, dep)
            elseif dep isa AbstractDict
                for (k, v) in _convert(Dict{String,Any}, dep)
                    if k == "version"
                        version = _convert(String, v)
                    elseif k == "channel"
                        channel = _convert(String, v)
                    else
                        error("deps keys must be 'version' or 'channel', got '$k'")
                    end
                end
            else
                error("deps must be String or Dict, got $(typeof(dep))")
            end
            pkg = PkgSpec(name, version=version, channel=channel)
            push!(packages, pkg)
        end
    end

    # channels
    channels = ChannelSpec[]
    if haskey(toml, "channels")
        chan_names = _convert(Vector{String}, toml["channels"])
        for name in chan_names
            push!(channels, ChannelSpec(name))
        end
    end

    # pip packages
    pip_packages = PipPkgSpec[]
    if haskey(toml, "pip")
        pip = _convert(Dict{String,Any}, toml["pip"])
        if haskey(pip, "deps")
            pip_deps = _convert(Dict{String,String}, pip["deps"])
            for (name, version) in pip_deps
                pkg = PipPkgSpec(name, version=version)
                push!(pip_packages, pkg)
            end
        end
    end

    # done
    return (packages=packages, channels=channels, pip_packages=pip_packages)
end

function current_packages()
    cmd = conda_cmd(`list -p $(envdir()) --json`)
    pkglist = JSON3.read(cmd)
    Dict(normalise_pkg(pkg.name) => pkg for pkg in pkglist)
end

function current_pip_packages()
    pkglist = withenv() do
        cmd = `$(which("pip")) list --format=json`
        JSON3.read(cmd)
    end
    Dict(normalise_pip_pkg(pkg.name) => pkg for pkg in pkglist)
end

"""
    status()

Show the status of the current environment.

This does not include dependencies from nested environments.
"""
function status(; io::IO=stderr)
    # collect information
    dfile = cur_deps_file()
    resolved = is_resolved()
    isnull = backend() == :Null
    pkgs, channels, pippkgs = parse_deps(read_deps(file=dfile))
    curpkgs = resolved && !isnull && !isempty(pkgs) ? current_packages() : nothing
    curpippkgs = resolved && !isnull && !isempty(pippkgs) ? current_pip_packages() : nothing
    blank = isempty(pkgs) && isempty(channels) && isempty(pippkgs)

    # print status
    printstyled(io, "CondaPkg Status", color=:light_green)
    printstyled(io, " ", dfile, bold=true)
    if blank
        print(io, " (empty)")
    end
    println(io)
    if !resolved
        printstyled(io, "Not Resolved", color=:yellow)
        println(io, " (resolve first for more information)")
    end
    if isnull
        printstyled(io, "Using the Null backend", color=:yellow)
        println(io, " (dependencies shown here are not being managed)")
    end
    if !isempty(pkgs)
        printstyled(io, "Packages", bold=true, color=:cyan)
        println(io)
        for pkg in pkgs
            print(io, "  ", pkg.name)
            if curpkgs !== nothing
                curpkg = get(curpkgs, pkg.name, nothing)
                if curpkg === nothing
                    printstyled(io, " uninstalled", color=:red)
                else
                    print(io, " v", curpkg.version)
                end
            end
            specparts = String[]
            pkg.version == "" || push!(specparts, pkg.version)
            pkg.channel == "" || push!(specparts, "channel=$(pkg.channel)")
            isempty(specparts) || printstyled(io, " (", join(specparts, ", "), ")", color=:light_black)
            println(io)
        end
    end
    if !isempty(channels)
        printstyled(io, "Channels", bold=true, color=:cyan)
        println(io)
        for chan in channels
            println(io, "  ", chan.name)
        end
    end
    if !isempty(pippkgs)
        printstyled(io, "Pip packages", bold=true, color=:cyan)
        println(io)
        for pkg in pippkgs
            print(io, "  ", pkg.name)
            if curpippkgs !== nothing
                curpkg = get(curpippkgs, pkg.name, nothing)
                if curpkg === nothing
                    printstyled(io, " uninstalled", color=:red)
                else
                    print(io, " v", curpkg.version)
                end
            end
            if !isempty(pkg.version)
                printstyled(io, " (", pkg.version, ")", color=:light_black)
            end
            println(io)
        end
    end
end

function add(pkgs::AbstractVector; resolve=true)
    toml = read_deps()
    for pkg in pkgs
        add!(toml, pkg)
    end
    write_deps(toml)
    STATE.resolved = false
    resolve && CondaPkg.resolve()
    return
end

add(pkg::Union{PkgSpec,PipPkgSpec,ChannelSpec}; resolve=true) = add([pkg], resolve=resolve)

function add!(toml, pkg::PkgSpec)
    deps = get!(Dict{String,Any}, toml, "deps")
    filter!(kv -> normalise_pkg(kv[1]) != pkg.name, deps)
    dep = Dict{String,Any}()
    if pkg.version != ""
        dep["version"] = pkg.version
    end
    if pkg.channel != ""
        dep["channel"] = pkg.channel
    end
    if issubset(keys(dep), ["version"])
        deps[pkg.name] = pkg.version
    else
        deps[pkg.name] = dep
    end
end

function add!(toml, channel::ChannelSpec)
    channels = get!(Vector{Any}, toml, "channels")
    push!(channels, channel.name)
end

function add!(toml, pkg::PipPkgSpec)
    pip = get!(Dict{String,Any}, toml, "pip")
    deps = get!(Dict{String,Any}, pip, "deps")
    filter!(kv -> normalise_pip_pkg(kv[1]) != pkg.name, deps)
    deps[pkg.name] = pkg.version
end

function rm(pkgs::AbstractVector; resolve=true)
    toml = read_deps()
    for pkg in pkgs
        rm!(toml, pkg)
    end
    write_deps(toml)
    STATE.resolved = false
    resolve && CondaPkg.resolve()
    return
end

rm(pkg::Union{PkgSpec,PipPkgSpec,ChannelSpec}; resolve=true) = rm([pkg], resolve=resolve)

function rm!(toml, pkg::PkgSpec)
    deps = get!(Dict{String,Any}, toml, "deps")
    n = length(deps)
    filter!(kv -> normalise_pkg(kv[1]) != pkg.name, deps)
    length(deps) < n || error("package not found: $(pkg.name)")
    isempty(deps) && delete!(toml, "deps")
end

function rm!(toml, channel::ChannelSpec)
    channels = get!(Vector{Any}, toml, "channels")
    n = length(channels)
    filter!(c -> normalise_channel(c) != channel.name, channels)
    length(channels) < n || error("channel not found: $(channel.name)")
    isempty(channels) && delete!(toml, "channels")
end

function rm!(toml, pkg::PipPkgSpec)
    pip = get!(Dict{String,Any}, toml, "pip")
    deps = get!(Dict{String,Any}, pip, "deps")
    n = length(deps)
    filter!(kv -> normalise_pip_pkg(kv[1]) != pkg.name, deps)
    length(deps) < n || error("pip package not found: $(pkg.name)")
    isempty(deps) && delete!(pip, "deps")
    isempty(pip) && delete!(toml, "pip")
end

"""
    add(pkg; version="", channel="", resolve=true)

Adds a dependency to the current environment.
"""
add(pkg::AbstractString; version="", channel="", resolve=true) = add(PkgSpec(pkg, version=version, channel=channel), resolve=resolve)

"""
    rm(pkg; resolve=true)

Removes a dependency from the current environment.
"""
rm(pkg::AbstractString; resolve=true) = rm(PkgSpec(pkg); resolve=resolve)

"""
    add_channel(channel; resolve=true)

Adds a channel to the current environment.
"""
add_channel(channel::AbstractString; resolve=true) = add(ChannelSpec(channel), resolve=resolve)

"""
    rm_channel(channel; resolve=true)

Removes a channel from the current environment.
"""
rm_channel(channel::AbstractString; resolve=true) = rm(ChannelSpec(channel), resolve=resolve)

"""
    add_pip(pkg; version="", resolve=true)

Adds a pip dependency to the current environment.

!!! warning

    Use conda dependencies instead if at all possible. Pip does not handle version
    conflicts gracefully, so it is possible to get incompatible versions.
"""
add_pip(pkg::AbstractString; version="", resolve=true) = add(PipPkgSpec(pkg, version=version), resolve=resolve)

"""
    rm_pip(pkg; resolve=true)

Removes a pip dependency from the current environment.
"""
rm_pip(pkg::AbstractString; resolve=true) = rm(PipPkgSpec(pkg), resolve=resolve)


text:
using MinimalRLCore

mutable struct Critterbot <: AbstractEnvironment
    num_steps::Int
    num_features::Int
    num_targets::Int
    sensors::Vector{Int}

    idx::Int
    data::Array{Float64}
end

function Critterbot(obs_sensors, target_sensors)
    all_sensors = vcat(obs_sensors, target_sensors)
    num_features = length(obs_sensors)
    return Critterbot(CritterbotUtils.numSteps(),
                      num_features, length(target_sensors), CritterbotUtils.getSensorIndices(all_sensors),
                      0, squish(CritterbotUtils.loadSensor(all_sensors)))
end

function Critterbot(obs_sensors, target_sensors, γs::AbstractArray)
    all_sensors = vcat(obs_sensors, target_sensors)
    num_features = length(obs_sensors)*length(γs)
    data = squish(vcat(CritterbotUtils.getReturns(obs_sensors, γs), CritterbotUtils.loadSensor(target_sensors)))

    return Critterbot(CritterbotUtils.numSteps(),
                      num_features, length(target_sensors), CritterbotUtils.getSensorIndices(all_sensors),
                      0, data)
end

Critterbot(obs_sensors, target_sensors, γ_str::String) =
    Critterbot(obs_sensors, target_sensors, eval(Meta.parse(γ_str)))

# Hack to use same features as targets; just duplicate the data in new cols
Critterbot(sensors::Vector{Int}) = Critterbot(sensors, sensors)
get_num_features(cb::Critterbot) = cb.num_features
get_num_targets(cb::Critterbot) = cb.num_targets

function MinimalRLCore.start!(cb::Critterbot)
    cb.idx = 1
    return MinimalRLCore.get_state(cb)
end

function MinimalRLCore.step!(cb::Critterbot)
    cb.idx += 1
    return MinimalRLCore.get_state(cb), MinimalRLCore.get_reward(cb)
end

# Data for each sensor in a row, so that we can access data for all sensors by col
MinimalRLCore.get_state(cb::Critterbot) = cb.data[1:cb.num_features, cb.idx]
MinimalRLCore.get_reward(cb::Critterbot) = cb.data[cb.num_features+1:end, cb.idx] #

mutable struct CritterbotTPC <: AbstractEnvironment
    num_steps::Int
    num_features::Int

    idx::Int
    obs_data::Array{Float64}
    rewards::Array{Float64}
    discounts::Array{Float64}
    the_all_seeing_eye::Array{Float64}
end

function CritterbotTPC(obs_sensors; γ=0.9875)
    # all_sensors = vcat(obs_sensors, target_sensors)

    volts = CritterbotUtils.loadSensor(["Motor$(i)0" for i in 0:2])
    cur = CritterbotUtils.loadSensor(["Motor$(i)2" for i in 0:2])

    rewards = sum(abs.(cur .* volts); dims=1)
    light3 = GVFN.CritterbotUtils.loadSensor("Light3")
    discounts = (light3 .< 1020) .* γ

    num_features = length(obs_sensors)
    feats = squish(CritterbotUtils.loadSensor(obs_sensors))

    the_all_seeing_eye = CritterbotUtils.loadSensor("powerToGoal-$(γ)")

    return CritterbotTPC(
        CritterbotUtils.numSteps(),
        num_features,
        0,
        feats,
        rewards,
        discounts,
        the_all_seeing_eye)
end

function CritterbotTPC(obs_sensors, γs::AbstractArray; γ=0.9875)
    # all_sensors = vcat(obs_sensors, target_sensors)
    num_features = length(obs_sensors)*length(γs)
    feats = squish(CritterbotUtils.getReturns(obs_sensors, γs))

    volts = CritterbotUtils.loadSensor(["Motor$(i)0" for i in 0:2])
    cur = CritterbotUtils.loadSensor(["Motor$(i)2" for i in 0:2])

    rewards = sum(abs.(cur .* volts); dims=1)
    light3 = GVFN.CritterbotUtils.loadSensor("Light3")
    discounts = (light3 .< 1020).*γ

    the_all_seeing_eye = CritterbotUtils.loadSensor("powerToGoal-$(γ)")

    return CritterbotTPC(
        CritterbotUtils.numSteps(),
        num_features,
        0,
        feats,
        rewards,
        discounts,
        the_all_seeing_eye)
end

CritterbotTPC(obs_sensors, γ_str::String) =
    CritterbotTPC(obs_sensors, eval(Meta.parse(γ_str)))

# Hack to use same features as targets; just duplicate the data in new cols
# CritterbotTPC(sensors::Vector{Int}) = Critterbot(sensors, sensors)
get_num_features(cb::CritterbotTPC) = cb.num_features
# get_num_targets(cb::CritterbotTPC) = 1

function MinimalRLCore.start!(cb::CritterbotTPC)
    cb.idx = 1
    return MinimalRLCore.get_state(cb)
end

function MinimalRLCore.step!(cb::CritterbotTPC)
    cb.idx += 1
    return MinimalRLCore.get_state(cb), MinimalRLCore.get_reward(cb)
end

# Data for each sensor in a row, so that we can access data for all sensors by col
MinimalRLCore.get_state(cb::CritterbotTPC) = cb.obs_data[:, cb.idx]
MinimalRLCore.get_reward(cb::CritterbotTPC) = [cb.rewards[cb.idx], cb.discounts[cb.idx]]
ground_truth(cb::CritterbotTPC) = [cb.the_all_seeing_eye[cb.idx]]


text:
const libjulia_path = joinpath(Sys.BINDIR, Base.LIBDIR, "libjulia.dylib")
@show isfile(libjulia_path)
ccall((:jl_init,libjulia_path),Cvoid,(Ptr{Cvoid},),C_NULL)

max_stars_repo_name: terasakisatoshi/juliaExer

max_issues_repo_name: terasakisatoshi/juliaExer

max_forks_repo_name: terasakisatoshi/juliaExer


text:
mountaincar = QuickMDP(
    function (s, a, rng)
        x, v = s
        vp = clamp(v + a*0.001 + cos(3*x)*-0.0025, -0.07, 0.07)
        xp = x + vp
        if xp > 0.5
            r = 100.0
        else
            r = -1.0
        end
        return (sp=(xp, vp), r=r)
    end,
    actions = [-1., 0., 1.],
    initialstate = Deterministic((-0.5, 0.0)),
    discount = 0.95,
    isterminal = s -> s[1] > 0.5,

    render = function (step)
        cx = step.s[1]
        cy = 0.45*sin(3*cx)+0.5
        car = (context(), circle(cx, cy+0.035, 0.035), fill("blue"))
        track = (context(), line([(x, 0.45*sin(3*x)+0.5) for x in -1.2:0.01:0.6]), stroke("black"))
        goal = (context(), star(0.5, 1.0, -0.035, 5), fill("gold"), stroke("black"))
        bg = (context(), rectangle(), fill("white"))
        ctx = context(0.7, 0.05, 0.6, 0.9, mirror=Mirror(0, 0, 0.5))
        return compose(context(), (ctx, car, track, goal), bg)
    end
)


text:
module FluxArchitectures

import Flux: reset!, trainable

using Reexport
@reexport using Flux

using BSON
using JuliennedArrays
using SliceMap
using LazyArtifacts

import Tables

include("./shared/StackedLSTM.jl")
include("./shared/Sequentialize.jl")
include("dataloader.jl")
include("DARNN.jl")
include("DSANet.jl")
include("LSTnet.jl")
include("TPALSTM.jl")

export DARNN, DSANet, LSTnet, TPALSTM, get_data, prepare_data

end # module

max_stars_repo_name: sdobber/FluxArchitectures

max_issues_repo_name: sdobber/FluxArchitectures

max_forks_repo_name: sdobber/FluxArchitectures


text:
using MathepiaOptimal
using Test

@testset "MathepiaOptimal.jl" begin
    # Write your tests here.
end


text:
using DiffEqBayes, CmdStan, DynamicHMC

using Distributions, BenchmarkTools
using OrdinaryDiffEq, RecursiveArrayTools, ParameterizedFunctions
using Plots

gr(fmt=:png)

f = @ode_def LotkaVolterraTest begin
    dx = a*x - b*x*y
    dy = -c*y + d*x*y
end a b c d

u0 = [1.0,1.0]
tspan = (0.0,10.0)
p = [1.5,1.0,3.0,1,0]

prob = ODEProblem(f,u0,tspan,p)
sol = solve(prob,Tsit5())

t = collect(range(1,stop=10,length=10))
sig = 0.49
data = convert(Array, VectorOfArray([(sol(t[i]) + sig*randn(2)) for i in 1:length(t)]))

scatter(t, data[1,:], lab="#prey (data)")
scatter!(t, data[2,:], lab="#predator (data)")
plot!(sol)

priors = [Truncated(Normal(1.5,0.5),0.5,2.5),Truncated(Normal(1.2,0.5),0,2),Truncated(Normal(3.0,0.5),1,4),Truncated(Normal(1.0,0.5),0,2)]

@btime bayesian_result_stan = stan_inference(prob,t,data,priors,num_samples=10_000,printsummary=false)

@btime bayesian_result_turing = turing_inference(prob,Tsit5(),t,data,priors,num_samples=10_000)

@btime bayesian_result_dynamichmc = dynamichmc_inference(prob,Tsit5(),t,data,priors,num_samples=10_000)

using DiffEqBenchmarks
DiffEqBenchmarks.bench_footer(WEAVE_ARGS[:folder],WEAVE_ARGS[:file])


text:
# vim:set ft=julia ts=4 sw=4 sts=4 autoindent:

# Header providing several conveniences for the testing scripts.
#
# Author:   Pontus Stenetorp    <pontus stenetorp se>
# Version:  2014-06-06

import Base.source_path

using Base.Test

srcpath = source_path()
srcdir  = dirname(srcpath)

repodir     = "$srcdir/../"
respath     = "$repodir/res/"
debugpath   = "$respath/debug.conllx"
talpath     = "$respath/talbanken/train.conllx"

# Inject the main source directory.
push!(LOAD_PATH, "$srcdir/../src")

max_stars_repo_name: ninjin/allen

max_issues_repo_name: ninjin/allen

max_forks_repo_name: ninjin/allen


text:
# This file is for Simulation object.

"""
    Simulation(model; simdir=tempdir(), simname=string(uuid4()), simprefix="Simulation-", logger=SimpleLogger())

Constructs a `Simulation` object for the simulation of `model`. The `Simulation` object is used to monitor the state of the simulation of the `model`. `simdir` is the path of the directory into which the simulation files(log, data files etc.) are recorded. `simname` is the name of the `Simulation` and `simprefix` is the prefix of the name of the `Simulation`. `logger` is used to log the simulation steps of the `model`. See also: [`Model`](@ref), [`Logging`](https://docs.julialang.org/en/v1/stdlib/Logging/)
"""
mutable struct Simulation{MD}
    model::MD
    path::String
    logger::Union{SimpleLogger, ConsoleLogger}
    state::Symbol
    retcode::Symbol
    name::String
    function Simulation(model; simdir=tempdir(), simname=string(uuid4()), simprefix="Simulation-",
        logger=SimpleLogger())
        name = simprefix * simname  # `get_instant()` may be used for time-based paths names.
        path = joinpath(simdir, name)
        ispath(path) || mkpath(path)
        check_writer_files(model, path, force=true)
        new{typeof(model)}(model, path, logger, :idle, :unknown, name)
    end
end

show(io::IO, sim::Simulation) = print(io, "Simulation(state:$(sim.state), retcode:$(sim.retcode), path:$(sim.path))")

##### Simulation checks

function check_writer_files(model, path; force=true)
    for node in filter(node-> isa(node.component, Writer), model.nodes)
        dirname(node.component.file.path) == path || mv(node.component, path, force=true)
    end
end

##### Simulation logging
"""
    setlogger(path, name; setglobal::Bool=true)

Returns a logger. `path` is the path and `name` is the name of the file of the logger. If `setglobal` is `true`, the returned logger is a global logger.

# Example
```julia
julia> logger = setlogger(tempdir(), "mylogger", setglobal=true)
Base.CoreLogging.SimpleLogger(IOStream(<file /tmp/mylogger>), Info, Dict{Any,Int64}())
```
"""
function setlogger(path::AbstractString, name::AbstractString; setglobal::Bool=true, loglevel::LogLevel=Logging.Info)
    io = open(joinpath(path, name), "w+")
    logger = SimpleLogger(io, loglevel)
    if setglobal
        global_logger(logger)
    end
    logger
end

"""
    closelogger(logger=global_logger())

Closes the `logger` the file of the `loggger`. See also: [`setlogger`](@ref)
"""
function closelogger(logger=global_logger())
    if isa(logger, AbstractLogger)
        close(logger.stream)
    end
end

##### Simulation reporting
"""
    report(simulation::Simulation)

Records the state of the `simulation` by writing all its fields into a data file. All the fields of the `simulation` is written into file. When the file is read back, the `simulation` object is constructed back. The data file is written under the path of the `simulation`.
"""
function report(simulation::Simulation)
    # Write simulation info.
    jldopen(joinpath(simulation.path, "report.jld2"), "w") do simreport
        simreport["name"] = simulation.name
        simreport["path"] = simulation.path
        simreport["state"] = simulation.state
        simreport["retcode"] = simulation.retcode

        # Save simulation model components.
        model = simulation.model
        components = [node.component for node in model.nodes]
        # foreach(unfasten!, filter(component->isa(component, AbstractSink), components))
        model_group = JLD2.Group(simreport, "model")
        model_group["id"] = string(simulation.model.id)
        model_group["name"] = string(simulation.model.name)
        model_group["clock"] = simulation.model.clock
        model_group["callbacks"] = simulation.model.callbacks
        model_blocks_group = JLD2.Group(simreport, "components")
        for component in filter(component->!isa(component, Writer), components)
            model_blocks_group[string(component.name)] = component
        end
    end
    # close(simreport)
end

##### SimulationError type
"""
    SimulationError(msg::String)

Thrown when an error occurs during a simulation.
"""
struct SimulationError <: Exception
    msg::String
end
Base.showerror(io::IO, err::SimulationError) = println(io, err.msg)


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

mutable struct ConnectionSharedKey <: SwaggerModel
    value::Any # spec type: Union{ Nothing, String } # spec name: value

    function ConnectionSharedKey(;value=nothing)
        o = new()
        validate_property(ConnectionSharedKey, Symbol("value"), value)
        setfield!(o, Symbol("value"), value)
        o
    end
end # type ConnectionSharedKey

const _property_map_ConnectionSharedKey = Dict{Symbol,Symbol}(Symbol("value")=>Symbol("value"))
const _property_types_ConnectionSharedKey = Dict{Symbol,String}(Symbol("value")=>"String")
Base.propertynames(::Type{ ConnectionSharedKey }) = collect(keys(_property_map_ConnectionSharedKey))
Swagger.property_type(::Type{ ConnectionSharedKey }, name::Symbol) = Union{Nothing,eval(Meta.parse(_property_types_ConnectionSharedKey[name]))}
Swagger.field_name(::Type{ ConnectionSharedKey }, property_name::Symbol) =  _property_map_ConnectionSharedKey[property_name]

function check_required(o::ConnectionSharedKey)
    (getproperty(o, Symbol("value")) === nothing) && (return false)
    true
end

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


text:
module KhepriFrame3DD
using KhepriBase

# functions that need specialization
include(khepribase_interface_file())

include("Frame3DD.jl")

function __init__()

#=  set_backend_family(default_wall_family(), autocad, acad_layer_family("Wall"))
  set_backend_family(default_slab_family(), autocad, acad_layer_family("Slab"))
  set_backend_family(default_roof_family(), autocad, acad_layer_family("Roof"))
  set_backend_family(default_beam_family(), autocad, acad_layer_family("Beam"))
  set_backend_family(default_column_family(), autocad, acad_layer_family("Column"))
  set_backend_family(default_door_family(), autocad, acad_layer_family("Door"))
  set_backend_family(default_panel_family(), autocad, acad_layer_family("Panel"))

  set_backend_family(default_table_family(), autocad, acad_layer_family("Table"))
  set_backend_family(default_chair_family(), autocad, acad_layer_family("Chair"))
  set_backend_family(default_table_chair_family(), autocad, acad_layer_family("TableChairs"))

  set_backend_family(default_curtain_wall_family(), autocad, acad_layer_family("CurtainWall"))
  set_backend_family(default_curtain_wall_family().panel, autocad, acad_layer_family("CurtainWall-Panel"))
  set_backend_family(default_curtain_wall_family().boundary_frame, autocad, acad_layer_family("CurtainWall-Boundary"))
  set_backend_family(default_curtain_wall_family().transom_frame, autocad, acad_layer_family("CurtainWall-Transom"))
  set_backend_family(default_curtain_wall_family().mullion_frame, autocad, acad_layer_family("CurtainWall-Mullion"))
  =#
  set_backend_family(default_truss_bar_family(),
     frame3dd,
     frame3dd_truss_bar_family(
       E=210e09,         # E (Young's modulus)
       G=81e09,          # G (Kirchoff's or Shear modulus)
       p=0.0,            # Roll angle
       d=7.701e3))       # RO (Density)

  add_current_backend(frame3dd)
end

end


text:
### matriz

A = [[1 0];[0 1]]
abs.(A)

for i = 1:100
    println(i)
end

f(x) = x^2

# ou

function quadrado(x)
    x = x^2
    return x
end

quadrado(2)

using Plots

x = randn(500)

plot(x)
scatter(x)

y = randn(200)
z = cumsum(y)

plot(z)

w = 0.5*x + randn(500)
d = zeros(500)
d[1] = 1

w = w + cumsum(randn(500))
vetor = ones(500)
vetor = cumsum(vetor)

y = sin(vetor) + randn(500)
plot(w)

plot!(x)
savefig("figura1.png")

using Optim

f(x) = 2*x[1]^2 - 5*x[1] + 40

solucao = optimize(f, ones(1), LBFGS())


text:
"""
    Statevector

Represents a pure state.
"""
struct Statevector <: AbstractArray{ComplexF64,1}
    state::Vector{ComplexF64}
    vec::Vector{ComplexF64}
    perm::Vector{Int}
    N::Int

    @doc """
        Statevector(v::Vector{ComplexF64})

    Create a `Statevector` object representing a pure quantum state.
    """
    function Statevector(v::Vector{ComplexF64})
        ispow2(length(v)) || error("Statevector must have length which is a power of 2")
        N = intlog2(length(v))
        new(v, zeros(ComplexF64, 2^N), collect(1:2^N), N)
    end
    @doc """
    Statevector(v::Vector{ComplexF64})

    Create a `Statevector` object representing a pure quantum state with 'N' qubits.
    """
    function Statevector(N::Int)
        new(zeros(ComplexF64, 2^N), zeros(ComplexF64, 2^N), collect(1:2^N), N)
    end
end

Base.setindex!(ψ::Statevector, v, i::Int64) = setindex!(ψ.state, v, i)
Base.getindex(ψ::Statevector, i::Int) = getindex(ψ.state, i)
Base.size(ψ::Statevector) = size(ψ.state)
Base.length(ψ::Statevector) = length(ψ.state)
Base.IndexStyle(ψ::Type{Statevector}) = IndexLinear()


text:
pages = [
    "Home" => "index.md",
    "Symbolic Modeling Tutorials" => Any["tutorials/ode_modeling.md",
                                         "tutorials/spring_mass.md",
                                         "tutorials/acausal_components.md",
                                         "tutorials/higher_order.md",
                                         "tutorials/tearing_parallelism.md",
                                         "tutorials/nonlinear.md",
                                         "tutorials/optimization.md",
                                         "tutorials/stochastic_diffeq.md",
                                         "tutorials/nonlinear_optimal_control.md",
                                         "tutorials/parameter_identifiability.md"],
    "ModelingToolkitize Tutorials" => Any["mtkitize_tutorials/modelingtoolkitize.md",
                                          "mtkitize_tutorials/modelingtoolkitize_index_reduction.md",
                                          "mtkitize_tutorials/sparse_jacobians.md"],
    "Basics" => Any["basics/AbstractSystem.md",
                    "basics/ContextualVariables.md",
                    "basics/Variable_metadata.md",
                    "basics/Composition.md",
                    "basics/Validation.md",
                    "basics/DependencyGraphs.md",
                    "basics/FAQ.md"],
    "System Types" => Any["systems/ODESystem.md",
                          "systems/SDESystem.md",
                          "systems/JumpSystem.md",
                          "systems/NonlinearSystem.md",
                          "systems/OptimizationSystem.md",
                          "systems/ControlSystem.md",
                          "systems/PDESystem.md"],
    "comparison.md",
    "internals.md",
]


text:
using FactorGraph
using SparseArrays
using Test

@testset "StaticGBP" begin
    ### Vanilla GBP with Damping and XLSX input
    gbp = continuousModel("data33_14.xlsx"; variance = 1e60)
    for iteration = 1:50
        messageFactorVariable(gbp)
        messageVariableFactor(gbp)
    end
    for iteration = 51:1000
        messageDampFactorVariable(gbp)
        messageVariableFactor(gbp)
    end
    marginal(gbp)
    exact = wls(gbp)
    displayData(gbp)
    @test maximum(gbp.inference.mean ./ exact.estimate) < 1.0000000009

    ### Kahan GBP with Damping and HDF5 input
    gbp = continuousModel("data33_14.h5"; variance = 1e60)
    for iteration = 1:50
        messageFactorVariableKahan(gbp)
        messageVariableFactorKahan(gbp)
    end
    for iteration = 51:1000
        messageDampFactorVariableKahan(gbp)
        messageVariableFactorKahan(gbp)
    end
    marginal(gbp)
    exact = wls(gbp)
    displayData(gbp)
    @test maximum(gbp.inference.mean ./ exact.estimate) < 1.0000000009

    ### Broadcast GBP with Damping and HDF5 input
    gbp = continuousModel("data33_14.h5"; variance = 1e10)
    for iteration = 1:50
        messageFactorVariableBroadcast(gbp)
        messageVariableFactorBroadcast(gbp)
    end
    for iteration = 51:1000
        messageDampFactorVariableBroadcast(gbp)
        messageVariableFactorBroadcast(gbp)
    end
    marginal(gbp)
    exact = wls(gbp)
    @test maximum(gbp.inference.mean ./ exact.estimate) < 1.00000000003

    ### Vanilla GBP with Damping with marginal computation in each iteration
    gbp = continuousModel("data33_14.h5"; variance = 1e60)
    for iteration = 1:50
        messageFactorVariable(gbp)
        messageVariableFactor(gbp)
        marginal(gbp)
    end
    for iteration = 51:1000
        messageDampFactorVariable(gbp)
        messageVariableFactor(gbp)
        marginal(gbp)
    end
    exact = wls(gbp)
    displayData(gbp, exact)
    @test maximum(gbp.inference.mean ./ exact.estimate) < 1.0000000009

    ### Vanilla GBP with Damping with error metrics
    gbp = continuousModel("data33_14.h5"; variance = 1e60)
    for iteration = 1:50
        messageFactorVariable(gbp)
        messageVariableFactor(gbp)
    end
    for iteration = 51:500
        messageDampFactorVariable(gbp)
        messageVariableFactor(gbp)
    end
    marginal(gbp)
    exact = wls(gbp)
    errors = errorMetric(gbp, exact)
    displayData(gbp, exact, errors)
    @test abs(errors.rmse[1] - exact.rmse[1]) < 1e-10
    @test abs(errors.mae[1] - exact.mae[1]) < 1e-10
    @test abs(errors.wrss[1] - exact.wrss[1]) < 1e-10
    @test errors.rmseGBPWLS[1] < 1e-10
    @test errors.maeGBPWLS[1] < 1e-10

    ### Vanilla GBP with Damping and passing arguments
    H = [3.0 2.0 -1.0; -2.0 2.0 1.0; 1.0 1.0 1.0]
    z = [6.0; 3.0; 4.0]
    v = [1.0; 1.0; 1.0]
    gbp = continuousModel(H, z, v)
    for iteration = 1:5
        messageFactorVariable(gbp)
        messageVariableFactor(gbp)
    end
    for iteration = 6:1000
        messageDampFactorVariable(gbp)
        messageVariableFactor(gbp)
    end
    marginal(gbp)
    @test round.(gbp.inference.mean, digits = 3) ≈ [1.0; 2.0; 1.0]

    ### Kahan GBP and passing arguments
    H = [1.0 0.0 0.0; 25 -25 0; -50 -40 90; 0 1 0]
    z = [0; 1.795; 1.966; -0.066]
    v = [1e-2; 1e-2; 1e-2; 1e-6]

    gbp = continuousModel(H, z, v, variance = 1e60)
    for iteration = 1:30
        messageFactorVariableKahan(gbp)
        messageVariableFactorKahan(gbp)
    end
    marginal(gbp)
    @test round.(gbp.inference.mean, digits = 5) ≈ [0.00579; -0.066; -0.00427]
end

@testset "DynamicGBP" begin
    H = [3.0 2.0 -1.0; -2.0 2.0 1.0; 1.0 1.0 1.0]
    z = [16.0; 13.0; 14.0]
    v = [10.0; 20.0; 30.0]

    ### Vanilla GBP with Damping and two dynamic updates
    gbp = continuousModel(H, z, v)
    for iteration = 1:9
        messageFactorVariable(gbp)
        messageVariableFactor(gbp)
    end
    dynamicFactor!(gbp; factor = 1, mean = 6, variance = 1)
    dynamicFactor!(gbp; factor = 3, mean = 4, variance = 1)
    for iteration = 10:99
        messageDampFactorVariable(gbp)
        messageVariableFactor(gbp)
    end
    dynamicFactor!(gbp; factor = 2, mean = 3, variance = 1)
    for iteration = 100:2000
        messageDampFactorVariable(gbp)
        messageVariableFactor(gbp)
    end
    marginal(gbp)
    @test round.(gbp.inference.mean, digits = 3) ≈ [1.0; 2.0; 1.0]

    ### Vanilla GBP with Damping and one dynamic updates
    gbp = continuousModel(H, z, v)
    for iteration = 1:9
        messageFactorVariable(gbp)
        messageVariableFactor(gbp)
        marginal(gbp)
    end
    dynamicFactor!(gbp; factor = 1, mean = 6, variance = 1)
    dynamicFactor!(gbp; factor = 3, mean = 4, variance = 1)
    dynamicFactor!(gbp; factor = 2, mean = 3, variance = 1)
    for iteration = 10:1000
        messageDampFactorVariable(gbp)
        messageVariableFactor(gbp)
        marginal(gbp)
    end
    @test round.(gbp.inference.mean, digits = 3) ≈ [1.0; 2.0; 1.0]
end

@testset "AgeingGBP" begin
    H = [1.0 0.0 0.0; 0.0 1.0 0.0; 0.0 0.0 1.0; 3.0 2.0 -1.0]
    z = [1.0; 2.0; 3.0; 11.0]
    v = [1.0; 1.0; 1.0; 1.0]

    ### Vanilla GBP with linear ageing
    gbp = continuousModel(H, z, v)
    for iteration = 1:100
        messageFactorVariable(gbp)
        messageVariableFactor(gbp)
    end
    for iteration = 1:900
        ageingVariance!(gbp; factor = 4, initial = 1, limit = 1e60, model = 1, a = 1e57, tau = iteration)
        messageFactorVariable(gbp)
        messageVariableFactor(gbp)
    end
    marginal(gbp)
    @test round.(gbp.inference.mean, digits = 3) ≈ [1.0; 2.0; 3.0]

    ### Vanilla GBP with logarithmic ageing
    gbp = continuousModel(H, z, v)
    for iteration = 1:100
        messageFactorVariable(gbp)
        messageVariableFactor(gbp)
    end
    for iteration = 1:1000
        ageingVariance!(gbp; factor = 4, initial = 1, limit = 1e60, model = 2, a = 1e57, b = 0.00002, tau = iteration)
        messageFactorVariable(gbp)
        messageVariableFactor(gbp)
    end
    marginal(gbp)
    @test round.(gbp.inference.mean, digits = 3) ≈ [1.0; 2.0; 3.0]

    ### Vanilla GBP with exponential ageing
    gbp = continuousModel(H, z, v)
    for iteration = 1:100
        messageFactorVariable(gbp)
        messageVariableFactor(gbp)
    end
    for iteration = 1:900
        ageingVariance!(gbp; factor = 4, initial = 1, limit = 1e60, model = 3, a = 0.08, b = 2, tau = iteration)
        messageFactorVariable(gbp)
        messageVariableFactor(gbp)
    end
    marginal(gbp)
    @test round.(gbp.inference.mean, digits = 3) ≈ [1.0; 2.0; 3.0]
end

@testset "GraphicalModel" begin
    ### Freeze factor node
    gbp = continuousModel("data33_14.h5")
    for i = 1:2
        messageFactorVariable(gbp)
        messageVariableFactor(gbp)
    end
    T = gbp.inference
    beforeFreeze = [T.fromFactor T.toVariable T.meanFactorVariable T.varianceFactorVariable]
    idx = T.fromFactor .== 10
    beforeFreeze = beforeFreeze[idx, :]
    freezeFactor!(gbp; factor = 10)
    for i = 1:20
        messageFactorVariable(gbp)
        messageVariableFactor(gbp)
    end
    T = gbp.inference
    afterFreeze = [T.fromFactor T.toVariable T.meanFactorVariable T.varianceFactorVariable]
    afterFreeze = afterFreeze[idx, :]
    @test sum(beforeFreeze[:, 3] - afterFreeze[:, 3]) == 0
    @test sum(beforeFreeze[:, 4] - afterFreeze[:, 4]) == 0

    ### Defreeze factor node
    defreezeFactor!(gbp; factor = 10)
    for i = 1:20
        messageFactorVariable(gbp)
        messageVariableFactor(gbp)
    end
    T = gbp.inference
    afterDefreeze = [T.fromFactor T.toVariable T.meanFactorVariable T.varianceFactorVariable]
    afterDefreeze = afterDefreeze[idx, :]
    @test sum(afterDefreeze[:, 3] - afterFreeze[:, 3]) != 0
    @test sum(afterDefreeze[:, 4] - afterFreeze[:, 4]) != 0

    ### Freeze variable node
    gbp = continuousModel("data33_14.h5")
    for i = 1:2
        messageFactorVariable(gbp)
        messageVariableFactor(gbp)
    end
    T = gbp.inference
    beforeFreeze = [T.fromVariable T.toFactor T.meanVariableFactor T.varianceVariableFactor]
    idx = T.fromVariable .== 9
    beforeFreeze = beforeFreeze[idx, :]
    freezeVariable!(gbp; variable = 9)
    for i = 1:20
        messageFactorVariable(gbp)
        messageVariableFactor(gbp)
    end
    T = gbp.inference
    afterFreeze = [T.fromVariable T.toFactor T.meanVariableFactor T.varianceVariableFactor]
    afterFreeze = afterFreeze[idx, :]
    @test sum(beforeFreeze[:, 3] - afterFreeze[:, 3]) == 0
    @test sum(beforeFreeze[:, 4] - afterFreeze[:, 4]) == 0

    ### Defreeze variable node
    defreezeVariable!(gbp; variable = 9)
    for i = 1:20
        messageFactorVariable(gbp)
        messageVariableFactor(gbp)
    end
    T = gbp.inference
    afterDefreeze = [T.fromFactor T.toVariable T.meanFactorVariable T.varianceFactorVariable]
    afterDefreeze = afterDefreeze[idx, :]
    @test sum(afterDefreeze[:, 3] - afterFreeze[:, 3]) != 0
    @test sum(afterDefreeze[:, 4] - afterFreeze[:, 4]) != 0

    ### Freeze and defreeze message variable node to factor node
    gbp = continuousModel("data33_14.h5")
    for iteration = 1:5
        messageFactorVariable(gbp)
        messageVariableFactor(gbp)
    end
    beforeFreeze6 = copy(gbp.inference.meanVariableFactor[6])
    beforeFreeze57 = copy(gbp.inference.meanVariableFactor[57])
    freezeVariableFactor!(gbp; variable = 9, factor = 29)
    freezeVariableFactor!(gbp; variable = 2, factor = 10)
    for iteration = 1:100
        messageFactorVariable(gbp)
        messageVariableFactor(gbp)
    end
    afterFreeze6 = copy(gbp.inference.meanVariableFactor[6])
    afterFreeze57 = copy(gbp.inference.meanVariableFactor[57])
    @test beforeFreeze6 == afterFreeze6
    @test beforeFreeze57 == afterFreeze57
    defreezeVariableFactor!(gbp; variable = 9, factor = 29)
    defreezeVariableFactor!(gbp; variable = 2, factor = 10)
    for iteration = 1:100
        messageFactorVariable(gbp)
        messageVariableFactor(gbp)
    end
    @test afterFreeze6 != gbp.inference.meanVariableFactor[6]
    @test afterFreeze57 != gbp.inference.meanVariableFactor[57]

    ### Freeze and defreeze message factor node to variable node
    gbp = continuousModel("data33_14.h5")
    for iteration = 1:5
        messageFactorVariable(gbp)
        messageVariableFactor(gbp)
    end
    beforeFreeze6 = copy(gbp.inference.meanFactorVariable[6])
    beforeFreeze57 = copy(gbp.inference.meanFactorVariable[57])
    freezeFactorVariable!(gbp; variable = 2, factor = 11)
    freezeFactorVariable!(gbp; variable = 14, factor = 25)
    for iteration = 1:100
        messageFactorVariable(gbp)
        messageVariableFactor(gbp)
    end
    afterFreeze6 = copy(gbp.inference.meanFactorVariable[6])
    afterFreeze57 = copy(gbp.inference.meanFactorVariable[57])
    @test beforeFreeze6 == afterFreeze6
    @test beforeFreeze57 == afterFreeze57
    defreezeFactorVariable!(gbp; variable = 2, factor = 11)
    defreezeFactorVariable!(gbp; variable = 14, factor = 25)
    for iteration = 1:100
        messageFactorVariable(gbp)
        messageVariableFactor(gbp)
    end
    @test afterFreeze6 != gbp.inference.meanFactorVariable[6]
    @test afterFreeze57 != gbp.inference.meanFactorVariable[57]

    ### Hide factor node
    H = [1.5 0.0 1.0;
         2.0 3.1 4.6;
         2.0 1.1 0.6;
         1.0 0.0 0.0;
         1.0 0.0 0.0]
    z = [0.5; 0.8; 4.1; 5.1; 10]
    v = [0.1; 1.0; 1.0; 0.0002;  0.0001]
    gbp = continuousModel(H,z,v)
    hideFactor!(gbp; factor = 2)
    for iteration = 1:200
        messageFactorVariable(gbp)
        messageVariableFactor(gbp)
    end
    marginal(gbp)
    exact = wls(gbp)
    @test round.(gbp.inference.mean, digits = 4) ≈ round.(exact.estimate, digits = 4)

    H = [1.5 0.0 1.0;
         2.0 3.1 4.6;
         2.0 1.1 0.6;
         1.0 0.0 0.0;
         1.0 0.0 0.0]
    z = [0.5; 0.8; 4.1; 5.1; 10]
    v = [0.1; 1.0; 1.0; 0.0002;  0.0001]
    gbp = continuousModel(H,z,v)
    hideFactor!(gbp; factor = 4)
    for iteration = 1:200
        messageFactorVariable(gbp)
        messageVariableFactor(gbp)
    end
    marginal(gbp)
    exact = wls(gbp)
    @test round.(gbp.inference.mean, digits = 4) ≈ round.(exact.estimate, digits = 4)

    ### Add factor node
    gbp = continuousModel("data33_14.h5")
    for iteration = 1:15
        messageFactorVariable(gbp)
        messageVariableFactor(gbp)
    end
    mean = [2.2; 3.1; 0.5]
    variance = [0.001; 0.001; 0.001]
    jacobian = zeros(3, 14)
    jacobian[1, 5] = 2.0
    jacobian[2, 1] = 0.45; jacobian[2, 5] = 0.23; jacobian[2, 14] = 0.1
    jacobian[3, 3] = 0.25; jacobian[3, 6] = 0.8
    addFactors!(gbp; mean = mean, variance = variance, jacobian = jacobian)
    for iteration = 1:200
        messageFactorVariable(gbp)
        messageVariableFactor(gbp)
    end
    marginal(gbp)
    exact = wls(gbp)
    @test round.(gbp.inference.mean, digits = 4) ≈ round.(exact.estimate, digits = 4)

    ### Hide and add factor node
    gbp = continuousModel("data33_14.h5")
    for iteration = 1:15
        messageFactorVariable(gbp)
        messageVariableFactor(gbp)
    end
    hideFactor!(gbp; factor = 2)
    mean = [2.2; 3.1; 0.5]
    variance = [0.001; 0.001; 0.001]
    jacobian = zeros(3, 14)
    jacobian[1, 5] = 2.0
    jacobian[2, 1] = 0.45; jacobian[2, 5] = 0.23; jacobian[2, 14] = 0.1
    jacobian[3, 3] = 0.25; jacobian[3, 6] = 0.8
    addFactors!(gbp; mean = mean, variance = variance, jacobian = jacobian)
    for iteration = 1:200
        messageFactorVariable(gbp)
        messageVariableFactor(gbp)
    end
    marginal(gbp)
    exact = wls(gbp)
    @test round.(gbp.inference.mean, digits = 4) ≈ round.(exact.estimate, digits = 4)
end

@testset "GraphicalModelTree" begin
    H = [2 3 0 0 0 0 0 0 0 0 0;
         0 4 6 8 2 3 0 0 0 0 0;
         0 0 0 5 0 0 0 0 0 0 0;
         0 0 0 0 2 0 0 0 0 0 0;
         0 0 0 0 0 8 0 2 0 1 0;
         0 0 0 0 0 4 0 0 0 0 0;
         0 0 0 0 0 3 7 0 2 0 0;
         0 0 0 0 0 0 0 5 0 0 0;
         0 0 0 0 0 0 0 3 0 0 0;
         2 0 0 0 0 0 0 0 0 0 0;
         2 0 0 0 0 0 0 0 0 0 0;
         1 0 0 0 0 0 0 0 0 0 2]
    z = [1; 2; 3; 4; 5; 6; 7; 8; 9; 10; 11; 12]
    v = [3; 4; 2; 5; 1; 8; 9; 1; 2; 3; 8; 4]

    ### isTree
    gbp = continuousTreeModel(H, z, v)
    tree = isTree(gbp)
    @test tree == true

    gbp = continuousModel("data33_14.h5")
    tree = isTree(gbp)
    @test tree == false

    # Forward-backward algorithm
    gbp = continuousModel(H, z, v)
    for iteration = 1:100
        messageFactorVariable(gbp)
        messageVariableFactor(gbp)
    end
    marginal(gbp)
    mean = copy(gbp.inference.mean)

    gbp = continuousTreeModel(H, z, v; root = 1)
    while gbp.graph.forward
        forwardVariableFactor(gbp)
        forwardFactorVariable(gbp)
    end
    while gbp.graph.backward
        backwardVariableFactor(gbp)
        backwardFactorVariable(gbp)
    end
    marginal(gbp)
    @test sum(round.(mean - gbp.inference.mean, digits = 4)) == 0

    # Forward-backward algorithm - new root
    gbp = continuousTreeModel(H, z, v; root = 8)
    while gbp.graph.forward
        forwardVariableFactor(gbp)
        forwardFactorVariable(gbp)
    end
    while gbp.graph.backward
        backwardVariableFactor(gbp)
        backwardFactorVariable(gbp)
    end
    marginal(gbp)
    @test sum(round.(mean - gbp.inference.mean, digits = 4)) == 0

    # Forward-backward algorithm - new root
    gbp = continuousTreeModel(H, z, v; root = 2)
    while gbp.graph.forward
        forwardVariableFactor(gbp)
        forwardFactorVariable(gbp)
    end
    while gbp.graph.backward
        backwardVariableFactor(gbp)
        backwardFactorVariable(gbp)
    end
    marginal(gbp)
    @test sum(round.(mean - gbp.inference.mean, digits = 4)) == 0
end

@testset "DicreteModelTree" begin
    ### isTree
    bp = discreteTreeModel("discrete6_4.xlsx")
    tree = isTree(bp)
    @test tree == true

    bp = discreteTreeModel("discrete6_4.h5")
    tree = isTree(bp)
    @test tree == true

    # Forward-backward algorithm
    bp = discreteTreeModel("discrete6_4.h5")
    while bp.graph.forward
        forwardVariableFactor(bp)
        forwardFactorVariable(bp)
    end
    while bp.graph.backward
        backwardVariableFactor(bp)
        backwardFactorVariable(bp)
    end
    marginal(bp)
end


text:
# ---
# title: 652. Find Duplicate Subtrees
# id: problem652
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Tree
# link: <https://leetcode.com/problems/find-duplicate-subtrees/description/>
# hidden: true
# ---
#
# Given the `root` of a binary tree, return all **duplicate subtrees**.
#
# For each kind of duplicate subtrees, you only need to return the root node of
# any **one** of them.
#
# Two trees are **duplicate** if they have the **same structure** with the
# **same node values**.
#
#
#
# **Example 1:**
#
# ![](https://assets.leetcode.com/uploads/2020/08/16/e1.jpg)
#
#
#
#     Input: root = [1,2,3,4,null,2,4,null,null,4]
#     Output: [[2,4],[4]]
#
#
# **Example 2:**
#
# ![](https://assets.leetcode.com/uploads/2020/08/16/e2.jpg)
#
#
#
#     Input: root = [2,1,1]
#     Output: [[1]]
#
#
# **Example 3:**
#
# ![](https://assets.leetcode.com/uploads/2020/08/16/e33.jpg)
#
#
#
#     Input: root = [2,2,2,3,null,3,null]
#     Output: [[2,3],[3]]
#
#
#
#
# **Constraints:**
#
#   * The number of the nodes in the tree will be in the range `[1, 10^4]`
#   * `-200 <= Node.val <= 200`
#
#
## @lc code=start
using LeetCode

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


text:
using BinaryBuilder

cuda_version = v"10.1.243"   # NOTE: could be less specific

sources_linux_x64 = [
    ArchiveSource("https://developer.nvidia.com/compute/machine-learning/cudnn/secure/8.0.4/10.1_20200923/cudnn-10.1-linux-x64-v8.0.4.30.tgz",
                  "eb4b888e61715168f57a0a0a21c281ada6856b728e5112618ed15f8637487715")
]
sources_linux_ppc64le = [
    ArchiveSource("https://developer.nvidia.com/compute/machine-learning/cudnn/secure/8.0.4/10.1_20200923/cudnn-10.1-linux-ppc64le-v8.0.4.30.tgz",
                  "690811bbf04adef635f4a6f480575fc2a558c4a2c98c85c7090a3a8c60dacea9")
]
sources_windows = [
    ArchiveSource("https://developer.nvidia.com/compute/machine-learning/cudnn/secure/8.0.4/10.1_20200923/cudnn-10.1-windows10-x64-v8.0.4.30.zip",
                  "f9a7e6632ac1e2058676eed141d05ce34d8b479df388cdc80334356dce63bac2")
]

include("../common.jl")


text:
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/FEMBase.jl/blob/master/LICENSE

using FEMBase
using FEMBase: IP
using Test

@testset "test integration point" begin
    a = sqrt(1.0/3.0)
    ip = IP(1, 1.0, (-a, a))
    strain = [1.0 2.0; 3.0 4.0]
    update!(ip, "strain", 0.0 => strain)
    @test isapprox(ip("strain", 0.0), strain)
    @test isapprox(ip[1], -a)
    for p in ip
        @test isapprox(abs.(p), a)
    end
    update!(ip, "strain", 1.0 => 2*strain)
    @test isapprox(ip("strain", 1.0), 2*strain)
    ip2 = IP(1, 1.0, [-a, a])
    @test isapprox(ip2[1], -a)
end


text:
# ---
# title: 108. Convert Sorted Array to Binary Search Tree
# id: problem108
# author: Tian Jun
# date: 2020-10-31
# difficulty: Easy
# categories: Tree, Depth-first Search
# link: <https://leetcode.com/problems/convert-sorted-array-to-binary-search-tree/description/>
# hidden: true
# ---
#
# Given an array where elements are sorted in ascending order, convert it to a
# height balanced BST.
#
# For this problem, a height-balanced binary tree is defined as a binary tree in
# which the depth of the two subtrees of _every_ node never differ by more than
# 1.
#
# **Example:**
#
#
#
#     Given the sorted array: [-10,-3,0,5,9],
#
#     One possible answer is: [0,-3,9,-10,null,5], which represents the following height balanced BST:
#
#           0
#          / \
#        -3   9
#        /   /
#      -10  5
#
#
#
## @lc code=start
using LeetCode

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


text:
#!/usr/bin/env julia

export pe3

# need Primes.jl
# https://github.com/JuliaMath/Primes.jl
using Primes

pe3(n::Integer=600851475143)=factor(Vector,n)|>maximum

if !haskey(ENV,"PROJECT_EULER_WITHOUT_RESULT")
    pe3()|>print
end

max_stars_repo_path: 3. Largest prime factor/3.jl

max_stars_repo_name: miRoox/ProjectEuler

max_issues_repo_path: 3. Largest prime factor/3.jl

max_issues_repo_name: miRoox/ProjectEuler

max_forks_repo_path: 3. Largest prime factor/3.jl

max_forks_repo_name: miRoox/ProjectEuler


text:
using Images, ImageView, ImageMorphology

geom_img = load("sample-images/geometrical-figures-and-noise.jpg");
geom_img_binary = Gray.(Gray.(geom_img) .> 0.5);

geom_img_binary_o = opening(geom_img_binary)
imshow(geom_img_binary_e)

geom_img_binary_c = closing(geom_img_binary)
imshow(geom_img_binary_c)


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

module VersionWeights

importall ....Base.Operators

export VersionWeight

struct HierarchicalValue{T}
    v::Vector{T}
    rest::T
end

HierarchicalValue{T}(v::Vector{T}) = HierarchicalValue{T}(v, zero(T))
HierarchicalValue(T::Type) = HierarchicalValue(T[])

Base.zero{T}(::Type{HierarchicalValue{T}}) = HierarchicalValue(T)

Base.typemin{T}(::Type{HierarchicalValue{T}}) = HierarchicalValue(T[], typemin(T))

for f in (:-, :+)
    @eval function Base.$f{T}(a::HierarchicalValue{T}, b::HierarchicalValue{T})
        av = a.v
        bv = b.v
        la = length(a.v)
        lb = length(b.v)
        l0 = min(la, lb)
        l1 = max(la, lb)
        ld = la - lb
        rv = Array{T}(l1)
        rf = ($f)(a.rest, b.rest)
        @inbounds for i = 1:l0
            rv[i] = ($f)(av[i], bv[i])
        end
        @inbounds for i = l0+1:l0+ld
            rv[i] = ($f)(av[i], b.rest)
        end
        @inbounds for i = l0+1:l0-ld
            rv[i] = ($f)(a.rest, bv[i])
        end
        return HierarchicalValue(rv, rf)
    end
end

Base.:-(a::HierarchicalValue) = HierarchicalValue(-a.v, -a.rest)

function Base.cmp{T}(a::HierarchicalValue{T}, b::HierarchicalValue{T})
    av = a.v
    bv = b.v
    la = length(a.v)
    lb = length(b.v)
    l0 = min(la, lb)
    l1 = max(la, lb)
    ld = la - lb
    @inbounds for i = 1:l0
        c = cmp(av[i], bv[i]); c != 0 && return c
    end
    @inbounds for i = l0+1:l0+ld
        c = cmp(av[i], b.rest); c != 0 && return c
    end
    @inbounds for i = l0+1:l0-ld
        c = cmp(a.rest, bv[i]); c != 0 && return c
    end
    return cmp(a.rest, b.rest)
end
Base.isless{T}(a::HierarchicalValue{T}, b::HierarchicalValue{T}) = cmp(a,b) < 0
Base.:(==){T}(a::HierarchicalValue{T}, b::HierarchicalValue{T}) = cmp(a,b) == 0

Base.abs{T}(a::HierarchicalValue{T}) = HierarchicalValue(T[abs(x) for x in a.v], abs(a.rest))

struct VWPreBuildItem
    nonempty::Int
    s::HierarchicalValue{Int}
    i::Int
end
VWPreBuildItem() = VWPreBuildItem(0, HierarchicalValue(Int), 0)
VWPreBuildItem(i::Int) = VWPreBuildItem(1, HierarchicalValue(Int), i)
VWPreBuildItem(s::String) = VWPreBuildItem(1, HierarchicalValue(Int[s...]), 0)

Base.zero(::Type{VWPreBuildItem}) = VWPreBuildItem()

Base.typemin(::Type{VWPreBuildItem}) = (x=typemin(Int); VWPreBuildItem(x, typemin(HierarchicalValue{Int}), x))

Base.:-(a::VWPreBuildItem, b::VWPreBuildItem) = VWPreBuildItem(a.nonempty-b.nonempty, a.s-b.s, a.i-b.i)
Base.:+(a::VWPreBuildItem, b::VWPreBuildItem) = VWPreBuildItem(a.nonempty+b.nonempty, a.s+b.s, a.i+b.i)

Base.:-(a::VWPreBuildItem) = VWPreBuildItem(-a.nonempty, -a.s, -a.i)

function Base.cmp(a::VWPreBuildItem, b::VWPreBuildItem)
    c = cmp(a.nonempty, b.nonempty); c != 0 && return c
    c = cmp(a.s, b.s); c != 0 && return c
    return cmp(a.i, b.i)
end
Base.isless(a::VWPreBuildItem, b::VWPreBuildItem) = cmp(a,b) < 0
Base.:(==)(a::VWPreBuildItem, b::VWPreBuildItem) = cmp(a,b) == 0

Base.abs(a::VWPreBuildItem) = VWPreBuildItem(abs(a.nonempty), abs(a.s), abs(a.i))

struct VWPreBuild
    nonempty::Int
    w::HierarchicalValue{VWPreBuildItem}
end

const _vwprebuild_zero = VWPreBuild(0, HierarchicalValue(VWPreBuildItem))

function VWPreBuild(ispre::Bool, desc::Tuple{Vararg{Union{Int,String}}})
    isempty(desc) && return _vwprebuild_zero
    desc == ("",) && return VWPreBuild(ispre ? -1 : 1, HierarchicalValue(VWPreBuildItem[]))
    nonempty = ispre ? -1 : 0
    w = Array{VWPreBuildItem}(length(desc))
    i = 1
    @inbounds for item in desc
        w[i] = VWPreBuildItem(item)
        i += 1
    end
    return VWPreBuild(nonempty, HierarchicalValue(w))
end
VWPreBuild() = _vwprebuild_zero

Base.zero(::Type{VWPreBuild}) = VWPreBuild()

const _vwprebuild_min = VWPreBuild(typemin(Int), typemin(HierarchicalValue{VWPreBuildItem}))
Base.typemin(::Type{VWPreBuild}) = _vwprebuild_min

function Base.:-(a::VWPreBuild, b::VWPreBuild)
    b === _vwprebuild_zero && return a
    a === _vwprebuild_zero && return -b
    VWPreBuild(a.nonempty-b.nonempty, a.w-b.w)
end
function Base.:+(a::VWPreBuild, b::VWPreBuild)
    b === _vwprebuild_zero && return a
    a === _vwprebuild_zero && return b
    VWPreBuild(a.nonempty+b.nonempty, a.w+b.w)
end

function Base.:-(a::VWPreBuild)
    a === _vwprebuild_zero && return a
    VWPreBuild(-a.nonempty, -a.w)
end

@inline function Base.cmp(a::VWPreBuild, b::VWPreBuild)
    a === _vwprebuild_zero && b === _vwprebuild_zero && return 0
    c = cmp(a.nonempty, b.nonempty); c != 0 && return c
    return cmp(a.w, b.w)
end
Base.isless(a::VWPreBuild, b::VWPreBuild) = cmp(a,b) < 0
Base.:(==)(a::VWPreBuild, b::VWPreBuild) = cmp(a,b) == 0

function Base.abs(a::VWPreBuild)
    a === _vwprebuild_zero && return a
    VWPreBuild(abs(a.nonempty), abs(a.w))
end

# The numeric type used to determine how the different
# versions of a package should be weighed
struct VersionWeight
    major::Int
    minor::Int
    patch::Int
    prerelease::VWPreBuild
    build::VWPreBuild
end
VersionWeight(major::Int, minor::Int, patch::Int, prerelease::VWPreBuild) = VersionWeight(major, minor, patch, prerelease, zero(VWPreBuild))
VersionWeight(major::Int, minor::Int, patch::Int) = VersionWeight(major, minor, patch, zero(VWPreBuild))
VersionWeight(major::Int, minor::Int) = VersionWeight(major, minor, 0)
VersionWeight(major::Int) = VersionWeight(major, 0)
VersionWeight() = VersionWeight(0)

VersionWeight(vn::VersionNumber) =
    VersionWeight(vn.major, vn.minor, vn.patch,
                  VWPreBuild(true, vn.prerelease), VWPreBuild(false, vn.build))

Base.zero(::Type{VersionWeight}) = VersionWeight()

Base.typemin(::Type{VersionWeight}) = (x=typemin(Int); y=typemin(VWPreBuild); VersionWeight(x, x, x, y, y))

Base.:-(a::VersionWeight, b::VersionWeight) =
    VersionWeight(a.major-b.major, a.minor-b.minor, a.patch-b.patch,
                  a.prerelease-b.prerelease, a.build-b.build)

Base.:+(a::VersionWeight, b::VersionWeight) =
    VersionWeight(a.major+b.major, a.minor+b.minor, a.patch+b.patch,
                  a.prerelease+b.prerelease, a.build+b.build)

Base.:-(a::VersionWeight) =
    VersionWeight(-a.major, -a.minor, -a.patch,
                  -a.prerelease, -a.build)

function Base.cmp(a::VersionWeight, b::VersionWeight)
    c = cmp(a.major, b.major); c != 0 && return c
    c = cmp(a.minor, b.minor); c != 0 && return c
    c = cmp(a.patch, b.patch); c != 0 && return c
    c = cmp(a.prerelease, b.prerelease); c != 0 && return c
    return cmp(a.build, b.build)
end
Base.isless(a::VersionWeight, b::VersionWeight) = cmp(a,b) < 0
Base.:(==)(a::VersionWeight, b::VersionWeight) = cmp(a,b) == 0

Base.abs(a::VersionWeight) =
    VersionWeight(abs(a.major), abs(a.minor), abs(a.patch),
                  abs(a.prerelease), abs(a.build))

end