text:
const _chapman_so_filename =
            string(
                    @__DIR__, "/chapman",
                    isapple ? "_macos" : islinux ? "_linux" : "_win",
                    ".so",
                )

_check_chapman() = _check_so(_chapman_so_filename)
_check_chapman()

# export atm8_chapman
atm8_chapman(xscale::Float64, chi::Float64)::Float64 =
    ccall((:atm8_chapman_, _chapman_so_filename),
            Float64,
            (Ref{Float64}, Ref{Float64}),
            xscale, chi,)

function test_chapman()
    ch1 = atm8_chapman(1000.0, 30.0)
    ch1
    # isapprox(ch1, 1.154317169142657)
end


text:
function coordinate_tag(n)
  str = replace("$n", ")" => "")
  str = replace(str, "(" => "")
  str = replace(str, " " => "")
  if length(n) > 1
    str = replace(str, "," => ".")
  else
    str = replace(str, "," => "")
  end
  return str
end

function link_tag(n1, n2)
  link_string = "$(coordinate_tag(n1))↔$(coordinate_tag(n2))"
  start_ind = nextind(link_string, 0, 1)
  stop_ind = min(ncodeunits(link_string), nextind(link_string, 0, 8))
  link_string = link_string[start_ind:stop_ind]
  return TagSet(link_string)
end

function ITensors.linkinds(lattice::HyperCubic; linkdims, addtags=ts"")
  dims = size(lattice)
  N = length(dims)
  linkinds_dict = Dict{Edge{N},Index{typeof(linkdims)}}()
  for n in sites(lattice), edge_n in incident_edges(lattice, n)
    l = Index(linkdims; tags=ITensors.addtags(link_tag(edge_n.edge...), addtags))
    get!(linkinds_dict, edge_n, l)
  end
  return linkinds_dict
end

function get_link_ind(linkinds_dict::Dict, edge::Edge, site::Tuple)
  l = linkinds_dict[edge]
  return is_in_edge(site, edge) ? dag(l) : l
end

# A network of link indices for a HyperCubic lattice, with
# no site indices.
function inds_network(dims::Int...; linkdims, kwargs...)
  site_inds = fill(Index{typeof(linkdims)}[], dims)
  return inds_network(site_inds; linkdims=linkdims, kwargs...)
end

function inds_network(site_inds::Array{<:Index,N}; kwargs...) where {N}
  return inds_network(map(x -> [x], site_inds); kwargs...)
end

# A network of link indices for a HyperCubic lattice, with
# site indices specified.
function inds_network(
  site_inds::Array{<:Vector{<:Index},N}; linkdims, addtags=ts""
) where {N}
  dims = size(site_inds)
  lattice = HyperCubic(dims)
  linkinds_dict = linkinds(lattice; linkdims=linkdims, addtags=addtags)
  inds = Array{Vector{Index{typeof(linkdims)}},N}(undef, dims)
  for n in sites(lattice)
    inds_n = [
      get_link_ind(linkinds_dict, edge_n, n) for edge_n in incident_edges(lattice, n)
    ]
    inds[n...] = append!(inds_n, site_inds[n...])
  end
  return inds
end

ext: jl


text:
# This file is a part of BAT.jl, licensed under the MIT License (MIT).

using BAT
using Test

using Random, Statistics, LinearAlgebra
using Distributions
using QuadGK, ForwardDiff
using StableRNGs

@testset "standard_normal" begin
    stblrng() = StableRNG(789990641)

    @test @inferred(BAT.LogUniform(0.1f0, 100)) isa BAT.LogUniform{Float32}
    @test_throws ArgumentError BAT.LogUniform(0, 100)
    @test_throws ArgumentError BAT.LogUniform(100, 10)

    d = BAT.LogUniform(0.1, 100)
    X = @inferred(rand(stblrng(), d, 10^5))

    @test params(d) == (0.1, 100)
    @test @inferred(minimum(d)) == 0.1
    @test @inferred(maximum(d)) == 100

    @test isapprox(@inferred(quantile(d, 0.3)), quantile(X, 0.3), rtol = 0.05)
    @test quadgk(x -> pdf(d, x), minimum(d), maximum(d))[1] ≈ 1
    @test cdf(d, minimum(d)) ≈ 0
    @test cdf(d, maximum(d)) ≈ 1

    @test isapprox(@inferred(mean(d)), mean(X), rtol = 0.05)
    @test isapprox(@inferred(median(d)), median(X), rtol = 0.05)
    @test @inferred(mode(d)) == minimum(d)
    @test isapprox(@inferred(var(d)), var(X), rtol = 0.05)
    @test isapprox(@inferred(std(d)), std(X), rtol = 0.05)

    for x in (minimum(d), maximum(d), (minimum(d) + maximum(d)) / 2.9, (minimum(d) + maximum(d)) / 5.2)
        @test @inferred(quantile(d, @inferred(cdf(d, x)))) ≈ x
        @test ForwardDiff.derivative(x -> cdf(d, x), x) ≈ @inferred(pdf(d, x))
        @test log(pdf(d, x)) ≈ @inferred(logpdf(d, x))
        @test log(cdf(d, x)) ≈ @inferred(logcdf(d, x))
    end
end


text:
using JuMP
using ParameterJuMP
using Polyhedra

mutable struct Conitope{T, VT <: AbstractVector{T}, D<:Polyhedra.FullDim}
    d::D
    points::Vector{VT}
    optimizer_constructor
    model::Union{Nothing, MOI.ModelLike}
    sum_con::Union{Nothing, MOI.ConstraintIndex{MOI.ScalarAffineFunction{T}, MOI.LessThan{T}}}
    z::Union{Nothing, Vector{MOI.ConstraintIndex{MOI.ScalarAffineFunction{T}, MOI.GreaterThan{T}}}}
    t_0::Union{Nothing, MOI.VariableIndex}
    function Conitope{T, VT, D}(
        d::Polyhedra.FullDim, points::Polyhedra.PointIt,
        optimizer_constructor=nothing) where {T, VT, D}
        new{T, VT, D}(Polyhedra.FullDim_convert(D, d),
                      Polyhedra.lazy_collect(points), optimizer_constructor, nothing, nothing, nothing)
    end
end
function Conitope(d::Polyhedra.FullDim, points::Polyhedra.PointIt, args...)
    return Conitope{Polyhedra.coefficient_type(eltype(points)), eltype(points),
                    typeof(d)}(d, points, args...)
end

Polyhedra.FullDim(v::Conitope) = v.d

function _sum(x::Vector{MOI.VariableIndex}, T::Type)
    MOI.ScalarAffineFunction([MOI.ScalarAffineTerm(one(T), xi) for xi in x],
                             zero(T))
end
function _convex_combination(t::Vector{MOI.VariableIndex}, points::Vector{Vector{T}}, d::Int) where T
    return [MOI.ScalarAffineFunction([
        MOI.ScalarAffineTerm(points[i][j], t[i]) for i in eachindex(t)
    ], zero(T)) for j in 1:d]
end
function _build_model(brp::Conitope{T}) where T
    @assert brp.model === nothing
    n = length(brp.points)
    brp.model = MOI.instantiate(brp.optimizer_constructor, with_bridge_type = T)
    t, ct = MOI.add_constrained_variables(brp.model, MOI.Nonnegatives(n))
    # sum t_i ≤ 1
    brp.sum_con = MOI.add_constraint(brp.model, _sum(t, T), MOI.LessThan(one(T)))
    return _convex_combination(t, brp.points, brp.d)
end
function _parametrized_model(brp::Conitope, v)
    hull = _build_model(brp)
    brp.z = [MOI.add_constraint(brp.model, hull[i], MOI.GreaterThan(v[i])) for i in 1:brp.d]
end
function _ratio_model(brp::Conitope{T}, v) where T
    hull = _build_model(brp)
    brp.t_0 = MOI.add_variable(brp.model)
    for i in 1:brp.d
        push!(hull[i].terms, MOI.ScalarAffineTerm(-v[i], brp.t_0))
    end
    brp.z = [MOI.add_constraint(brp.model, hull[i], MOI.GreaterThan(zero(T))) for i in 1:brp.d]
end
function _update_model(p::Conitope{T}, v) where T
    ts, ct = MOI.add_constrained_variables(p.model, MOI.Nonnegatives(1))
    t = ts[1]
    MOI.modify(p.model, p.sum_con, MOI.ScalarCoefficientChange(t, one(T)))
    for i in 1:p.d
        MOI.modify(p.model, p.z[i], MOI.ScalarCoefficientChange(t, v[i]))
    end
end

function _fix(brp::Conitope, v)
    for i in 1:brp.d
        if brp.t_0 === nothing
            MOI.set(brp.model, MOI.ConstraintSet(), brp.z[i], MOI.GreaterThan(v[i]))
        else
            MOI.modify(brp.model, brp.z[i], MOI.ScalarCoefficientChange(brp.t_0, -v[i]))
        end
    end
end

function _ei(i::Integer, d::Integer, T::Type)
    e = zeros(T, d)
    e[i] = one(T)
    return e
end
function Polyhedra.polyhedron(
    c::Conitope{T},
    lib::Polyhedra.Library=Polyhedra.default_library(Polyhedra.FullDim(c),
                                                     T)) where T
    n = fulldim(c)
    p = polyhedron(vrep([Ray(-_ei(i, n, T)) for i in 1:n]) + vrep(c.points),
                   lib)
    return p ∩ (Matrix(-LinearAlgebra.I, n, n) * p)
end


text:
function AbstractLambdaMethod()
    return AbstractLambdaMethod(())
end

function solve_ils(obj::AbstractLambdaMethod, arg0::jint, arg1::Vector{jdouble}, arg2::Vector{jint}, arg3::RealMatrix)
    return jcall(obj, "solveILS", Vector{IntegerLeastSquareSolution}, (jint, Vector{jdouble}, Vector{jint}, RealMatrix), arg0, arg1, arg2, arg3)
end


text:
function Qmat(N, m, lam)
    # Integration matrix
    Q = 1
    round(m) !== m ? (return Q = Qmat05(N, lam)) : nothing
    Q = Q*Qmat1(N, lam)
    return Q
end

function Qmat05(N, lam)
    # Half integral matrix
    if lam == 0.5
        v = (2/sqrt(pi)) ./ (2*collect(0:N-1).+1)
        Q = myspdiagm(spdiagm(0=>v, 1=>-v[2:end]), N)
    elseif lam == 1
        v = (sqrt(pi)/2) .* ones(N)
        Q = myspdiagm(spdiagm(-1=>v[2:end], 0=>v), N)
    end
    return Q
end

function Qmat1(N, lam)
    # Full integral matrix
    if lam == 0.5
        v = 1 ./ (2*collect(0:N-1).+1)
        Q = myspdiagm(spdiagm(-1 => v, 1 => -v[2:end]), N)
        Q[1, 1] = 1
    elseif lam == 1
        nn = collect(0:N)
        v = 1 ./ (2*nn.+1)
        v2 = 2 ./ (4*nn.^2 .-1)
        Q = myspdiagm(spdiagm(-1=>v[2:end], 0=>v2[2:end], 1=>-v[2:end-1]), N)
    end
    return Q
end

myspdiagm(M, N) = size(M, 1) !== N ? (return M[1:N, 1:N]) : (return M)


text:
export @dom_str, @js, @js_str

# adapted from Hiccup.jl
function cssparse(s)
    # e.g. s = "img#theid.class1.class2[src=image.jpg, alt=great pic]"
    # parse props first
    p = match(r"\[[^\]]+\]", s)
    props = Dict()
    if p != nothing
        m = strip(p.match, ['[',']'])
        props[:attributes] = Dict(map(x->Pair(split(x,r"\s*=\s*", limit=2)...), split(m, r",\s*")))
        # s is now just the "img#class1.c2.thirdclass"
        s = s[1:p.offset-1]
    end
    trimfirst(str) = str[2:end]
    id = match(r"#-?[_a-zA-Z][_a-zA-Z0-9-]*", s)
    id == nothing || (props[:id] = trimfirst(id.match))
    classes = collect(m.match for m in eachmatch(r"\.-?[_a-zA-Z][_a-zA-Z0-9-]*", s))
    isempty(classes) || (props[:className] = join(map(trimfirst, classes), " "))
    tagm = match(r"^[^\.#\[\]]+", s)
    tagm == nothing && error("Invalid tag syntax $s")
    tag = tagm.match
    return tag, props
end

function makedom(tag, props)
    d = if occursin(":", string(tag))
        ns, t = split(string(tag), ":")
        DOM(Symbol(ns), Symbol(t))
    else
        DOM(:html, tag)
    end
    function dom(args...; kwargs...)
        n = node(d, args...)(Dict(kwargs))
        isempty(props) ? n : n(props)
    end
end

"""
    dom"div.<class>#<id>[<prop>=<value>,...]"(x...; kw...)
"""
macro dom_str(sraw)
    str = Meta.parse(string('"', sraw, '"'))
    quote
        tagstr, props = WebIO.cssparse($(esc(str)))
        tag = Symbol(tagstr)
        WebIO.makedom(tag, props)
    end
end

struct JSString
    s::String
end

# Required to allow JSStrings to interpolate into normal strings correctly.
Base.print(io::IO, x::JSString) = print(io, x.s)

"""
    Interpolator(s)

An iterator that yields _interpolation chunks (_i.e._ a string if the chunk was a string literal and
an `Expr` if the chunk was an interpolation).

# Examples
"""
struct Interpolator{S<:AbstractString}
    s::S
end

Base.iterate(interp::Interpolator, args...) = iterate_interpolations(interp.s, args...)
Base.IteratorSize(::Interpolator) = Base.SizeUnknown()

"""
    iterate_interpolations(s[, i])

Advance the iterator to obtain the next string or expression to be interpolated. If no parts
remain (where a part is either a string literal or an interpolation), return nothing; otherwise,
return a tuple with either a string or an `Expr` and the new iteration state.

# Examples
See `test/syntax.jl`.
"""
function iterate_interpolations(s, i=firstindex(s))
    # Notation:
    # i0 - start index
    # l - last index (length)
    # c - the most recently iterated character
    # i - current index (while iterating character-by-character)
    #     note: this is actually the index of the character **after** c
    # j - the end of the current "chunk" of string literal (when set)
    # escape_next - true if we should escape the next character
    # escape_current - true if c was preceded by an escape

    i0 = i
    l = lastindex(s)

    # Return nothing if we're past the end.
    if i0 > l
        return nothing
    end

    # Get the initial character (we have to handle the first character specially).
    c, i = iterate(s, i0)

    # Special case: when the first character is a $, it's an interpolation, so we return the
    # `Expr` of the interpolated _thing_.
    # For example, if `s` is `$foo</script>` (and `i` is 1), we should return `:foo`.
    if c == '$'
        # Note: at this point, `i` is after the `$`, so we parse `foo</script>` and return `foo`
        return Meta.parse(s, i, greedy=false, raise=false)
    end

    escape_next = (c == '\\')
    escape_current = false

    # Iterate over characters, stopping when we hit a `$` or run iterate through end of the string.
    while (c != '$') && i <= l
        c, i = iterate(s, i)
        if escape_next
            escape_current = true
            break
        elseif c == '\\'
            escape_next = true
        end
    end

    # If we have `\$`, we **don't** want to interpolate.
    if c == '$' && escape_current
        # Suppose the string is `this.\$refs = $myvar`. We don't want to interpolate a variable
        # named refs, but we do want to interpolate `myvar`. We set up indices j and k as follows.
        # t h i s . \ $ r e f s ␣ = ␣ \ $ m y v a r
        # i0      j     i           k             l
        # We form the string as s[i0:j] * "$" * s[i:k].
        # We get k by continuing to iterate, starting at position i (which is past the $ and thus
        # treated like a literal string).
        j = prevind(s, i, 3)
        # If we're at the end of the string, `iterate_interpolations` returns `nothing`.
        next_s, k = i < l ? iterate_interpolations(s, i) : ("", i)
        return (s[i0:j] * "\$" * next_s), k
    end

    if escape_current
        # Julia's string literal escaping logic gets wonky when quotes are involved.
        # For example, `console.log("\\")` gets passed to us as `console.log("\")`
        # whereas `foo = '\\\\'` gets passed to us exactly as written. The issue seems to be that
        # whenever quotes are involved, it eagerly consumes the backslashes in a weird way; but the
        # upshot of this is that if a user ever writes `\"`, it gets passed to as as just `"`, so
        # if we ever see `\"`, we know that the user actually wrote `\\"`.
        if c == '"'
            # j is the index of the quote (_i.e._ s[j] = '"')
            j = prevind(s, i)
            # We get the next chunk of string...
            next_s, k = i < l ? iterate_interpolations(s, i) : ("", i)
            # And then concatenate it.
            return s[i0:j] * next_s, k
        end
        j = prevind(s, i, 3)
        return s[i0:j] * "$c", i
    end

    # We don't have a preceeding backslash, so we **do** want to interpolate the **next** piece.
    # We do this by returning the current chunk of string literal and setting the iteration pointer
    # to the `$` so that the next call to `iterate_interpolations` knows to return an `Expr`.
    if c == '$'
        # Suppose s is `alert($thing)`.
        # Indices:           j i
        # We get j by shifting back twice (once so that s[i] == '$', and again so that s[i] == '(').
        # We also shift the position of the iterator back one so that the next call will have the
        # string beginning with the $ and will thus start interpolating.
        j = prevind(s, i, 2)
        return s[i0:j], prevind(s, i)
    end

    # If we get here, we're at the end of the string and we didn't find any interpolations or
    # escapes; so just return the rest of the string.
    return s[i0:prevind(s, i)], i
end

"""
    tojs(x)

Returns a JSString object that constructs the same object as `x`
"""
tojs(x) = x

"""
    showjs(io, x)

Print Javascript code to `io` that constructs the equivalent of `x`.
"""
showjs(io, x::Any) = JSON.show_json(io, JSEvalSerialization(), x)
showjs(io, x::AbstractString) = write(io, JSON.json(x))

"""
    @js_str(s)

Create a `JSString` using a string literal and perform interpolations from Julia.

# Examples
```julia-repl
julia> mystr = "foo"; mydict = Dict("foo" => "bar", "spam" => "eggs");
julia> println(js"const myStr = \$mystr; const myObject = \$mydict;")
const myStr = "foo"; const myObject = {"spam":"eggs","foo":"bar"};
```
"""
macro js_str(s)
    writes = map(Interpolator(s)) do x
        if isa(x, AbstractString)
            # If x is a string, it was specified in the js"..." literal so let it
            # through as-is.
            :(write(io, $(x)))
        else
            # Otherwise, it's some kind of interpolation so we need to generate a
            # JavaScript representation of whatever it is/whatever it evaluates to.
            :(showjs(io, tojs($(esc(x)))))
       end
   end

   :(JSString(sprint(io->(begin; $(writes...) end))))
end

Base.string(s::JSString) = s.s
Base.:(==)(x::JSString, y::JSString) = x.s==y.s

JSON.lower(x::JSString) = JSON.lower(x.s)

const JSONContext = JSON.Writer.StructuralContext
const JSONSerialization = JSON.Serializations.CommonSerialization

struct JSEvalSerialization <: JSONSerialization end

const verbose_json = Ref(false)

# adapted (very slightly) from JSON.jl test/serializer.jl
function JSON.show_json(io::JSONContext, ::JSEvalSerialization, x::JSString)
    if verbose_json[]
        first = true
        for line in split(x.s, '\n')
            !first && JSON.indent(io)
            first = false
            Base.print(io, line)
        end
    else
        Base.print(io, x.s)
    end
end


text:
#Calculate the data transformations to the output
function y_transform(ys::Array{Float64},ml::ML_model)
    ml.yscale  = mean(ys); ys[:]=ys./ ml.yscale ; ys[:] = log.(ys);     # Scaling and log transform
 #   ml.yscale = std(ys,1); ml.ysmean  =mean(ys,1); ys[:]= (ys.-ml.ysmean)./ml.yscale; # Centering and scaling
    return ys
end

#Undo the data transformations to the output
function y_undo_transform(ys::Array{Float64}, ml::ML_model)
    ys[:] = exp.(ys); ys[:] = ys.*ml.yscale;    # Scaling and log transform
  #  ys[:]= ys.*ml.yscale.+ml.ysmean;     # Centering and scaling
    return ys
end

#Calculate the data transformations to the input x (feature vector)
function X_transform(X::Array{Float64}, ml::ML_model)
    ml.Xmean = mean(X,1);
  #  ml.Xstd = std(X,1);
  #  ind = (0.==ml.Xstd); ml.Xstd[ind] =1.0;  #replace 0 in std by 1 incase there is a constant feature
  #  X= (X.-ml.Xmean)./ml.Xstd; # Centering and scaling
    X= X.-ml.Xmean;
    sX = size(X);
    X = [X ones(sX[1],1)]; # Adding the bias term
    return X
end

#Apply the saved data transformations to a single data vector x (feature vector)
function x_apply_transform(x::Array{Float64}, ml::ML_model)
  #  x= (x.-ml.Xmean)./ml.Xstd; # Centering and scaling
    x =x.-ml.Xmean;
    x = [x 1];# Adding the bias term
   return x
end


text:
function set(str,var,val)
    !any(i->i==var, fieldnames(typeof(str))) && @warn "\n No field in struct $(typeof(str)) is named $(var)!\n Fields are $(fieldnames(typeof(str)))"
    dict = Dict{Symbol,Any}()
    for (i,name) in enumerate(fieldnames(typeof(str)))
        if typeof(getfield(str,name))<:Parameter && name == var
            dict[name] = Parameter(val,getfield(getfield(str,name),:description))
        elseif typeof(getfield(str,name))<:Parameter
            dict[name] = Parameter(getfield(getfield(str,name),:value),getfield(getfield(str,name),:description))
        elseif typeof(getfield(str,name))<:BoundedParameter && name == var
            dict[name] = BoundedParameter(val,getfield(getfield(str,name),:description),getfield(getfield(str,name),:lb),getfield(getfield(str,name),:ub))
        elseif typeof(getfield(str,name))<:BoundedParameter
            dict[name] = BoundedParameter(getfield(getfield(str,name),:value),getfield(getfield(str,name),:description),getfield(getfield(str,name),:lb),getfield(getfield(str,name),:ub))
        elseif name == var
            dict[name] = val
        else
            dict[name] = getfield(str,name)
        end
    end
    x = NamedTuple{Tuple(keys(dict))}(values(dict))
    str = typeof(str)(;x...)
end

function (←)(str,tpl::Tuple)
    if typeof(values(tpl)[1]) ≠ Symbol
        for value in values(tpl)
            str = set(str,value[1],value[2])
        end
    else
        str = set(str,tpl[1],tpl[2])
    end
    return str
end


text:
"""
julia dqn_with_target.jl

This example implements the DQN algorithm from
`Playing atari with deep reinforcement learning.`
Mnih, V., Kavukcuoglu, K., Silver, D., Graves, A., Antonoglou, I., Wierstra, D., & Riedmiller, M. (2013).
arXiv preprint arXiv:1312.5602. It also shows the usage of a target network for
the target q value predictions.
"""
module DQNWITHTARGET

using Gym, ArgParse, Knet

include("replay_buffer.jl")
include("mlp.jl")
include("piecewise_schedule.jl")

function loss(w, states, actions, targets; nh=1)
    qvals = predict_q(w, states; nh=1)
    nrows = size(qvals, 1)
    index = actions + nrows*(0:(length(actions)-1))
    qpred = reshape(qvals[index], size(targets)...)
    mse = sum(abs2, targets-qpred) / size(states, 2)
    return mse
end

lossgradient = gradloss(loss)

function train!(w, prms, states, actions, targets; nh=1)
    g, mse = lossgradient(w, states, actions, targets; nh=nh)
    update!(w, g, prms)
    return mse
end

function dqn_learn(w, opts, env, buffer, exploration, o)
    total = 0.0
    readytosave = 10000
    episode_rewards = Float32[]
    frames = Float32[]
    ob_t = reset!(env)

    target_w = deepcopy(w)

    for fnum=1:o["frames"]
        o["render"] && render(env)
        #process the raw ob
        ob_t_reshaped = reshape(ob_t, size(ob_t)..., 1)
        if !o["play"] && rand() < value(exploration, fnum)
            a = sample(env.action_space)
        else
            obses_t = encode_recent(buffer, ob_t_reshaped; stack=o["stack"])
            inp = convert(o["atype"], obses_t)
            qvals = predict_q(w, inp; nh=length(o["hiddens"]))
            a = findmax(vec(Array(qvals)))[2] - 1
        end

        ob_t, reward, done, _ = step!(env, a)
        total += reward

        if !o["play"]
            #process the raw ob
            ob_tp1_reshaped = reshape(ob_t, size(ob_t)..., 1)
            push!(buffer, ob_t_reshaped, a+1, reward, ob_tp1_reshaped, done)

            if can_sample(buffer, o["bs"])
                obses_t, actions, rewards, obses_tp1, dones = sample_batch(buffer, o["bs"]; stack=o["stack"])
                obses_tp1 = convert(o["atype"], obses_tp1)
                #predict next q values with the target network
                nextq = predict_q(target_w, obses_tp1; nh=length(o["hiddens"]))
                nextq = Array(nextq)
                maxs = maximum(nextq; dims=1)
                nextmax = sum(nextq .* (nextq.==maxs); dims=1)
                nextmax = reshape(nextmax, 1, length(nextmax))
                targets = reshape(rewards,1,length(rewards)) .+ (o["gamma"] .* nextmax .* dones)
                obses_t = convert(o["atype"], obses_t)
                targets = convert(o["atype"], targets)
                mse = train!(w, opts, obses_t, actions, targets; nh=length(o["hiddens"]))
            end

            if o["save"] != "" && fnum > readytosave
                save_model(w, o["save"])
                readytosave += 10000
            end

            if fnum % o["tupdate"] == 0
                target_w = deepcopy(w)
            end
        end

        if done
            ob_t = reset!(env)
            o["printinfo"] && println("Frame: $fnum , Total reward: $total, Exploration Rate: $(value(exploration, fnum))")
            push!(episode_rewards, total)
            push!(frames, fnum)
            total = 0.0
        end
    end
    return episode_rewards, frames
end

function main(args=ARGS)
    s = ArgParseSettings()
    s.description = "(c) Ozan Arkan Can, 2018. An implementation of the deep q network with a target network."
    @add_arg_table s begin
        ("--frames"; arg_type=Int; default=100; help="number of frames")
        ("--lr"; arg_type=Float64; default=0.001; help="learning rate")
        ("--gamma"; arg_type=Float64; default=0.99; help="discount factor")
        ("--hiddens"; arg_type=Int; nargs='+'; default=[32]; help="number of units in the hiddens for the mlp")
        ("--env_id"; default="CartPole-v0")
        ("--render"; action=:store_true)
        ("--memory"; arg_type=Int; default=1000; help="memory size")
        ("--bs"; arg_type=Int; default=32; help="batch size")
        ("--stack"; arg_type=Int; default=4; help="length of the frame history")
        ("--save"; default=""; help="model name")
        ("--load"; default=""; help="model name")
        ("--usegpu"; action=:store_true; help="use GPU or not")
        ("--play"; action=:store_true; help="only play")
        ("--printinfo"; action=:store_true; help="print the training messages")
        ("--tupdate"; arg_type=Int; default=500; help="update frequency for the target network")
    end
    isa(args, AbstractString) && (args=split(args))
    if in("--help", args) || in("-h", args)
        ArgParse.show_help(s; exit_when_done=false)
        return
    end

    o = parse_args(args, s)
    o["atype"] = !o["usegpu"] ? Array{Float32} : KnetArray{Float32}
    Knet.seed!(12345)
    env = GymEnv(o["env_id"])
    seed!(env, 12345)

    INPUT = env.observation_space.shape[1] * o["stack"]
    OUTPUT = env.action_space.n

    if o["load"] == ""
        w = init_weights(INPUT, o["hiddens"], OUTPUT, o["atype"])
    else
        w = load_model(o["load"], o["atype"])
    end

    opts = Dict()
    for k in keys(w)
        opts[k] = Rmsprop(lr=o["lr"])
    end

    buffer = ReplayBuffer(o["memory"])

    exploration = PiecewiseSchedule([(0, 1.0),
                                     (round(Int, o["frames"]/5), 0.1),
                                     (round(Int, o["frames"]/3.5), 0.1)])

    rewards, frames = dqn_learn(w, opts, env, buffer, exploration, o)
end

PROGRAM_FILE == "dqn_with_target.jl" && main(ARGS)

end


text:
using Revise
using mujoco_sim

target = ([0.5,0.,0.3], [1.,0.,0.,1.])
joint_names=["joint_1", "joint_2", "joint_3", "joint_4", "joint_5",
                    "joint_6"]
m = get_model("models/arm.xml")
d = get_data(m)
sim = MJSim(m,d)

q = inverse_kinematics!(sim, "gripsite", target[1], target[2],  joint_names)
1

# function ctrler!(s)
#         x = [s.d.qpos[1], s.d.qvel[1],s.d.qpos[2], s.d.qvel[2]]
#         if s.d.time < 0.25
#             u = [rand(Uniform(-0.5, 0.5))]
#         end
#         setaction!(s, u)
# end
# simulate( sim, controller = ctrler!; mode="passive")


text:
module KhepriFreeCAD
using KhepriBase
using Sockets

# functions that need specialization
include(khepribase_interface_file())
include("FreeCAD.jl")

function __init__()
  add_current_backend(freecad)
end
end


text:
#this file contains the functions in Mapping.jl

#the function naming convention is [filetype]_[year(s)], where all refers to all for the type
#examples:
# all_all, pf_all, pf_2003-2009_2011
# refers to all files, all private foundation files,
# and the 2003-2009 and 2011 private foundation files
# the first function is always all_all, which is the only one with no
# dataframe argument

#holds the 1:1 mapping

const MAPPING = Dict{Symbol, Dict{Symbol, Symbol}}(
  :all_all=>Dict(
    :ein=>:ein,
    :name => :name,
    :fisyr=>:fisyr),

  :pf_all=>Dict(
    :p1intrev => :cashinterest,
    :p1divid => :dividends,
    :p1ginvpf => :grprof,
    :p1othinc => :otherinc,
    :p1totrev => :totrev,
    :p1netinv => :returnamt,
    :p1tcont => :contributions,
    :p1contpd => :distributions,
    :p1texmex => :pfexpcharity,
    :p1totexp => :totexp,
    :p2gvtinv => :treasuries,
    :p2crpstk => :stock,
    :p2crpbnd => :bonds,
    :p2totast => :bookassets,
    :p2tliabl => :bookliabilities,
    ),

  :pf_1997_2015=>Dict(
    :p1invrv => :invgross,
    :p1nadjrv => :pfopinc,
    :p1grents => :rents,
    :p1tinvex => :invexpense,
    :p1tadjex => :pfexpincome,
    :p2cashbv => :cash,
    :p2mortg => :mortgages,
    :p2eyotin => :otherexcrealestate,
    :p3eytfnd => :pfnetassets,
    :p10tasst => :totnoncharassets,
    :taxper => :datefiscal,
    :afcd => :affiliationcode
    ),

  :pf_1998_2015=>Dict(
    :ntee1=>:category
    ),

  :pf_1992_2015=>Dict(
    :p10mir => :pfminreturn,
    :p2tasfmv => :marketassets
    ),

  :pf_2011_2015=>Dict(
    :naics=>:naics
    ),

  :pc_all=>Dict(
    :invinc => :totinv,
    :ntee1 => :category,
    :fundbal => :netassets,
    :totrev => :totrev,
    :secur => :secursales,
    :salesexp => :securexp,
    :exps => :totexp,
    :cont => :contributions,
    :grprof => :grprof,
    :ass_eoy => :bookassets,
    :liab_eoy => :bookliabilities,
    :taxper => :datefiscal,
    :afcd => :affiliationcode
    ),

  :pc_1991_2015=>Dict(
    :progrev => :progrev,
    :fundinc => :specrev,
    :othinc => :othinc
    ),

  :pc_1997_2015=>Dict(
    :salesecn => :securrealized,
    :saleothn => :othrealized,
    :netrent => :netrent
    ),

  :pc_2009_2015=>Dict(
    :naics=>:naics
    ),

  :co_all=>Dict(
    :ntee1 => :category
    ),
  :co_1997_2015=>Dict(
    :taxper => :datefiscal
    ),

  :co_2000_2015=>Dict(
    :invinc => :totinv,
    :fundbal => :netassets,
    :salesexp => :securexp,
    :salesecn => :securrealized,
    :saleothn => :othrealized,
    :totrev2 => :totrev,
    :exps => :totexp,
    :cont => :contributions,
    :progrev => :progrev,
    :grprof => :grprof,
    :fundinc => :specrev,
    :othinc => :othinc,
    :netrent => :netrent,
    :ass_eoy => :bookassets,
    :liab_eoy => :bookliabilities,
    :afcd => :affiliationcode
    ),

  :co_2011_2015=>Dict(
    :naics=>:naics
    ),
  :co_1989_1999=>Dict(
    :p1naseoy => :netassets,
    :grossrec => :totrev

    )

)

replacemissing(v::Union{Number, Missing} , default::Union{Number, Missing} = 0.) = ismissing(v) ? default : v

function all_all(data::NCCSFile)
  local outDF::DataFrame #the output dataframe
  local fields::Dict #contains the name mapping

  try
    fields = MAPPING[:all_all]
    outDF = data.df[collect(keys(fields))]
    rename!(outDF, fields)

    rename!(outDF, :fisyr=>:fisyrold)
    outDF[:fisyr] = ((s::String)->parse(Int, s)).(outDF[:fisyrold])
    deletecols!(outDF, :fisyrold)

    outDF[:fileyear] = data.meta.fileYear
  catch err
    @error "$err \nMetadata: $(data.meta)"
  end

  #println(outDF[1:3, :name])
  return outDF
end

#public charities
function pc_all(data::NCCSFile, outDF::DataFrame)
  fields::Dict = MAPPING[:pc_all]

  #1:1 mappings
  outDF = hcat(outDF, data.df[collect(keys(fields))])
  rename!(outDF, fields)

  outDF[ :charitytype] = fill(:pc, size(outDF, 1))
  #outDF[ :invrealized] =  outDF[ :totinv] .+ outDF[ :secursales] .-  outDF[ :securexp]
  outDF[ :netincome] =  outDF[ :totrev] .- outDF[:totexp]

  return outDF
end

function pc_1991_2015(data::NCCSFile, outDF::DataFrame)
  fields::Dict = MAPPING[:pc_1991_2015]

  #1:1 mappings
  outDF = hcat(outDF, data.df[collect(keys(fields))])
  rename!(outDF, fields)

  return outDF
end

function pc_1997_2015(data::NCCSFile, outDF::DataFrame)
  fields::Dict = MAPPING[:pc_1997_2015]

  #1:1 mappings
  outDF = hcat(outDF, data.df[collect(keys(fields))])
  rename!(outDF, fields)

  #For data quality reasons measure this two different ways
  outDF[:invrealizeddirect] =  outDF[:securrealized] .+  outDF[:othrealized]  .+
    outDF[:totinv] .+ (replacemissing).(outDF[:netrent])
  outDF[:invrealizedindirect] = outDF[:totrev]  .- outDF[:contributions] .- outDF[:grprof]
  outDF[:invrealizedindirect] .-=  (replacemissing).(outDF[:progrev]) .+
    (replacemissing).(outDF[:specrev]) .+ (replacemissing).(outDF[:othinc])

  return outDF
end

function pc_2009_2015(data::NCCSFile, outDF::DataFrame)
  fields::Dict = MAPPING[:pc_2009_2015]

  #1:1 mappings
  outDF = hcat(outDF, data.df[collect(keys(fields))])
  rename!(outDF, fields)

  return outDF
end

function co_all(data::NCCSFile, outDF::DataFrame)
  #get targeted fields
  fields::Dict = MAPPING[:co_all]

  #1:1 mappings
  outDF = hcat(outDF, data.df[collect(keys(fields))])
  rename!(outDF, fields) #make the appropriate names
  outDF[:charitytype] = fill(:co, size(outDF, 1))

  return outDF
end

function co_1997_2015(data::NCCSFile, outDF::DataFrame)
  #get targeted fields
  fields::Dict = MAPPING[:co_1997_2015]

  #1:1 mappings
  outDF = hcat(outDF, data.df[collect(keys(fields))])
  rename!(outDF, fields) #make the appropriate names

  return outDF
end

function co_2000_2015(data::NCCSFile, outDF::DataFrame)
  #get targeted fields
  fields::Dict = MAPPING[:co_2000_2015]

  #1:1 mappings
  outDF = hcat(outDF, data.df[collect(keys(fields))])
  rename!(outDF, fields) #make the appropriate names

  outDF[ :netincome] =  outDF[ :totrev] .- outDF[:totexp]
  #For data quality reasons measure this two different ways
  outDF[:invrealizeddirect] =  outDF[:securrealized] .+  outDF[:othrealized]  .+
    outDF[:totinv] .+ (replacemissing).(outDF[:netrent])
  outDF[:invrealizedindirect] = outDF[:totrev]  .- outDF[:contributions] .- outDF[:grprof]
  outDF[:invrealizedindirect] .-=  (replacemissing).(outDF[:progrev]) .+
    (replacemissing).(outDF[:specrev]) .+ (replacemissing).(outDF[:othinc])

  return outDF
end

function co_2011_2015(data::NCCSFile, outDF::DataFrame)
  #get targeted fields
  fields::Dict = MAPPING[:co_2011_2015]

  #1:1 mappings
  outDF = hcat(outDF, data.df[collect(keys(fields))])
  rename!(outDF, fields) #make the appropriate names

  return outDF
end

function co_1989_1999(data::NCCSFile, outDF::DataFrame)
  #get targeted fields
  fields::Dict = MAPPING[:co_1989_1999]

  #1:1 mappings
  outDF = hcat(outDF, data.df[collect(keys(fields))])
  rename!(outDF, fields) #make the appropriate names

  return outDF
end

#private foundations
function pf_all(data::NCCSFile, outDF::DataFrame)
  fields::Dict = MAPPING[:pf_all]

  #1:1 mappings
  outDF = hcat(outDF, data.df[collect(keys(fields))])
  rename!(outDF, fields)

  #transformations
  outDF[:returncheckamt] =  outDF[:totrev] .- outDF[:contributions] .- outDF[:grprof]

  outDF[:netassetsbook] = outDF[:bookassets] .- outDF[:bookliabilities]
  outDF[:netincome] =  outDF[ :totrev] .- outDF[:totexp]

  outDF[:charitytype] = fill(:pf, size(outDF, 1))

  return outDF
end
function pf_2011_2015(data::NCCSFile, outDF::DataFrame)
  #get targeted fields
  fields::Dict = MAPPING[:pf_2011_2015]

  #1:1 mappings
  outDF = hcat(outDF, data.df[collect(keys(fields))])
  rename!(outDF, fields) #make the appropriate names

  return outDF
end

#private foundation files 1997-2015
function pf_1998_2015(data::NCCSFile, outDF::DataFrame)
  #get targeted fields
  fields::Dict = MAPPING[:pf_1998_2015]

  #1:1 mappings
  outDF = hcat(outDF, data.df[collect(keys(fields))])
  rename!(outDF, fields) #make the appropriate names

  return outDF
end

#private foundation files 1997-2015
function pf_1997_2015(data::NCCSFile, outDF::DataFrame)
  #get targeted fields
  fields::Dict = MAPPING[:pf_1997_2015]

  #1:1 mappings
  outDF = hcat(outDF, data.df[collect(keys(fields))])
  rename!(outDF, fields) #make the appropriate names

  #other transformations
  outDF[:pfgrossfromnet] = outDF[:returnamt] .- outDF[:invexpense] # =outDF[:invgross]

  return outDF
end

#private foundation files 1992-2015
function pf_1992_2015(data::NCCSFile, outDF::DataFrame)
  #get targeted fields
  fields::Dict = MAPPING[:pf_1992_2015]

  #1:1 mappings
  outDF = hcat(outDF, data.df[collect(keys(fields))])
  rename!(outDF, fields) #make the appropriate names
  outDF[:netassets] = outDF[:marketassets] .- outDF[:bookliabilities]

  return outDF
end

max_stars_repo_name: clintonTE/CCA

max_issues_repo_name: clintonTE/CCA

max_forks_repo_name: clintonTE/CCA


text:
function theme_minimal()
    Theme(
        Axis = (
            backgroundcolor = :transparent,
            xgridvisible = false,
            ygridvisible = false,
            xminorgridvisible = false,
            yminorgridvisible = false,
            leftspinevisible = true,
            rightspinevisible = false,
            bottomspinevisible = true,
            topspinevisible = false,
            xminorticksvisible = false,
            yminorticksvisible = false,
            xticksvisible = false,
            yticksvisible = false,
            xlabelpadding = 3,
            ylabelpadding = 3
        ),
        Legend = (
            framevisible = false,
            padding = (0, 0, 0, 0),
        ),
        Axis3 = (
            xgridvisible = false,
            ygridvisible = false,
            zgridvisible = false,
            xspinesvisible = true,
            yspinesvisible = true,
            zspinesvisible = true,
            yzpanelvisible = false,
            xzpanelvisible = false,
            xypanelvisible = false,
            xticksvisible = false,
            yticksvisible = false,
            zticksvisible = false,
            xticklabelpad = 3,
            yticklabelpad = 3,
            zticklabelpad = 6,
            xspinecolor_2 = :transparent,
            xspinecolor_3 = :transparent,
            yspinecolor_2 = :transparent,
            yspinecolor_3 = :transparent,
            zspinecolor_2 = :transparent,
            zspinecolor_3 = :transparent,
        ),
        Colorbar = (
            ticksvisible = false,
            spinewidth = 0,
            ticklabelpad = 5,
        )
    )
end


text:
### ALSMMALA state subtypes

## MuvALSMMALAState holds the internal state ("local variables") of the ALSMMALA sampler for multivariate parameters

type MuvALSMMALAState <: MuvPSMMALAState
  pstate::ParameterState{Continuous, Multivariate}
  tune::MCTunerState
  sqrttunestep::Real
  ratio::Real
  μ::RealVector
  newinvtensor::RealMatrix
  oldinvtensor::RealMatrix
  cholinvtensor::RealLowerTriangular
  newfirstterm::RealVector
  oldfirstterm::RealVector
  presentupdatetensor::Bool
  pastupdatetensor::Bool
  count::Integer
  updatetensorcount::Integer

  function MuvALSMMALAState(
    pstate::ParameterState{Continuous, Multivariate},
    tune::MCTunerState,
    sqrttunestep::Real,
    ratio::Real,
    μ::RealVector,
    newinvtensor::RealMatrix,
    oldinvtensor::RealMatrix,
    cholinvtensor::RealLowerTriangular,
    newfirstterm::RealVector,
    oldfirstterm::RealVector,
    presentupdatetensor::Bool,
    pastupdatetensor::Bool,
    count::Integer,
    updatetensorcount::Integer
  )
    if !isnan(ratio)
      @assert 0 < ratio < 1 "Acceptance ratio should be between 0 and 1"
      @assert sqrttunestep > 0 "Square root of tuned drift step is not positive"
    end
    new(
      pstate,
      tune,
      sqrttunestep,
      ratio,
      μ,
      newinvtensor,
      oldinvtensor,
      cholinvtensor,
      newfirstterm,
      oldfirstterm,
      presentupdatetensor,
      pastupdatetensor,
      count,
      updatetensorcount
    )
  end
end

MuvALSMMALAState(pstate::ParameterState{Continuous, Multivariate}, tune::MCTunerState=PSMMALAMCTune()) =
  MuvALSMMALAState(
  pstate,
  tune,
  NaN,
  NaN,
  Array(eltype(pstate), pstate.size),
  Array(eltype(pstate), pstate.size, pstate.size),
  Array(eltype(pstate), pstate.size, pstate.size),
  RealLowerTriangular(Array(eltype(pstate), pstate.size, pstate.size)),
  Array(eltype(pstate), pstate.size),
  Array(eltype(pstate), pstate.size),
  false,
  false,
  0,
  0
)

### Metropolis-adjusted Langevin Algorithm (ALSMMALA)

immutable ALSMMALA <: PSMMALA
  driftstep::Real
  identitymala::Bool
  update!::Function
  transform::Union{Function, Void}
  initupdatetensor::Tuple{Bool,Bool} # The tuple ordinates refer to (sstate.presentupdatetensor, sstate.pastupdatetensor)

  function ALSMMALA(
    driftstep::Real,
    identitymala::Bool,
    update!::Function,
    transform::Union{Function, Void},
    initupdatetensor::Tuple{Bool,Bool}
    )
    @assert driftstep > 0 "Drift step is not positive"
    new(driftstep, identitymala, update!, transform, initupdatetensor)
  end
end

ALSMMALA(
  driftstep::Real=1.;
  identitymala::Bool=false,
  update::Function=rand_update!,
  transform::Union{Function, Void}=nothing,
  initupdatetensor::Tuple{Bool,Bool}=(false, false)
) =
  ALSMMALA(driftstep, identitymala, update, transform, initupdatetensor)

### Initialize ALSMMALA sampler

## Initialize parameter state

function initialize!(
  pstate::ParameterState{Continuous, Multivariate},
  parameter::Parameter{Continuous, Multivariate},
  sampler::ALSMMALA,
  outopts::Dict
)
  parameter.uptotensorlogtarget!(pstate)
  if sampler.transform != nothing
    pstate.tensorlogtarget = sampler.transform(pstate.tensorlogtarget)
  end
  pstate.diagnosticvalues[1] = false
  @assert isfinite(pstate.logtarget) "Log-target not finite: initial value out of support"
  @assert all(isfinite(pstate.gradlogtarget)) "Gradient of log-target not finite: initial values out of support"
  @assert all(isfinite(pstate.tensorlogtarget)) "Tensor of log-target not finite: initial values out of support"

  if !isempty(outopts[:diagnostics])
    pstate.diagnostickeys = copy(outopts[:diagnostics])
    pstate.diagnosticvalues = Array(Any, length(pstate.diagnostickeys))
  end
end

## Initialize ALSMMALA state

tuner_state(sampler::ALSMMALA, tuner::PSMMALAMCTuner) =
  PSMMALAMCTune(
    BasicMCTune(NaN, 0, 0, tuner.smmalatuner.period),
    BasicMCTune(NaN, 0, 0, tuner.malatuner.period),
    BasicMCTune(sampler.driftstep, 0, 0, tuner.totaltuner.period)
  )

function sampler_state(
  parameter::Parameter{Continuous, Multivariate},
  sampler::ALSMMALA,
  tuner::MCTuner,
  pstate::ParameterState{Continuous, Multivariate},
  vstate::VariableStateVector
)
  sstate = MuvALSMMALAState(generate_empty(pstate, parameter.diffmethods, parameter.diffopts), tuner_state(sampler, tuner))
  sstate.sqrttunestep = sqrt(sampler.driftstep)
  sstate.oldinvtensor = inv(pstate.tensorlogtarget)
  sstate.cholinvtensor = ctranspose(chol(Hermitian(sstate.oldinvtensor)))
  sstate.oldfirstterm = sstate.oldinvtensor*pstate.gradlogtarget
  sstate.presentupdatetensor, sstate.pastupdatetensor = sampler.initupdatetensor
  sstate
end

### Reset ALSMMALA sampler

## Reset parameter state

function reset!(
  pstate::ParameterState{Continuous, Multivariate},
  x::RealVector,
  parameter::Parameter{Continuous, Multivariate},
  sampler::ALSMMALA
)
  pstate.value = copy(x)
  parameter.uptotensorlogtarget!(pstate)
  pstate.diagnosticvalues[1] = false
end

## Reset sampler state

function reset!(
  sstate::PSMMALAState,
  pstate::ParameterState{Continuous, Multivariate},
  parameter::Parameter{Continuous, Multivariate},
  sampler::ALSMMALA,
  tuner::MCTuner
)
  reset!(sstate.tune, sampler, tuner)
  sstate.sqrttunestep = sqrt(sampler.driftstep)
  sstate.oldinvtensor = inv(pstate.tensorlogtarget)
  sstate.cholinvtensor = chol(sstate.oldinvtensor, Val{:L})
  sstate.oldfirstterm = sstate.oldinvtensor*pstate.gradlogtarget
  sstate.presentupdatetensor, sstate.pastupdatetensor = sampler.initupdatetensor
end

Base.show(io::IO, sampler::ALSMMALA) = print(io, "ALSMMALA sampler: drift step = $(sampler.driftstep)")

Base.show(io::IO, ::MIME"text/plain", sampler::ALSMMALA) = show(io, sampler)


text:
using .MathLibSignatures
using ..Utilities.FastEigen

using LinearAlgebra

function MathLibSignatures.project(to::SemidefiniteCone, mat::Array{Float64, 2};
                                   eigs_implementation = KrylovMethods())
  if isnothing(rank(to))
    # do a full projection
    # @warn("Doing full eigen projection, could be costly!")
    eigDecomposition    = eigen(mat);
    posEigenValuesIndex = findall(x -> real(x) > 0,eigDecomposition.values);
    posEigenValues      = eigDecomposition.values[posEigenValuesIndex];
    posEigenVectors     = eigDecomposition.vectors[:,posEigenValuesIndex];
    projectedMatrix     = posEigenVectors * diagm(0 => posEigenValues) *posEigenVectors';
    return projectedMatrix;
  end
  # we are computing the full projection
  Λ, X = eigs(eigs_implementation, mat, nev = to.rank)
  # return X * diagm(0 => Λ) * X'
  id = findall(x -> real(x) > 0, Λ)
  return X[:, id] * diagm(0 => Λ[id]) * (X[:,id])'
end

function MathLibSignatures.project(to::ClosedInterval, x::AutoboxedArray{Float64})
  return max.(to.ll, min.(x, to.rr))
end

function MathLibSignatures.project!(to::ClosedInterval, x::Array{Float64})
  @. x = max.(to.ll, min.(x, to.rr))
end


text:
# TODO: needs more work later
"""
    transmittance(θ_in::FT, n::FT) where {FT <:AbstractFloat}

Transmittance of S and P polarizations for a plane dielectric surface for light incident, given
- `θ_i` Incoming light angle
- `n` Relative index of refraction
"""
function transmittance(θ_i::FT, n::FT) where {FT <:AbstractFloat}
    _cos_θ_i = cosd(θ_i);
    _cos_θ_t = sqrt(1 - sind(θ_i)^2 / n^2);

    # reflection of s and p polarizations
    _r_s = ((_cos_θ_i - n * _cos_θ_t) / (_cos_θ_i + n * _cos_θ_t)) ^ 2;
    _r_p = ((n * _cos_θ_i - _cos_θ_t) / (n * _cos_θ_i + _cos_θ_t)) ^ 2;

    return 1 - _r_s, 1 - _r_p
end

# TODO: needs more work later
###############################################################################
#
# Compute transmission of isotropic radiation
#
###############################################################################
"""
    calctav(α::FT, nr::FT) where {FT<:AbstractFloat}

Computes transmission of isotropic radiation across an interface between two
    dielectrics (Stern F., 1964; Allen W.A., 1973)). From calctav.m in
    PROSPECT-D
- `α` angle of incidence
- `nr` Index of refraction
"""
function calctav(α::FT, nr::FT) where {FT<:AbstractFloat}
    a   = ((nr+1) ^ 2) / 2;
    a3  = a  ^ 3;
    n2  = nr ^ 2;
    n4  = nr ^ 4;
    n6  = nr ^ 6;
    np  = n2 + 1;
    np2 = np ^ 2;
    np3 = np ^ 3;
    nm2 = (n2 - 1) ^2;
    k   = ((n2-1) ^ 2) / -4;
    k2  = k  ^ 2;
    sa2 = sind(α) ^ 2;

    _b1 = (α==90 ? 0 : sqrt( (sa2 - np/2)^2 + k ));
    _b2 = sa2 - np/2;
    b   = _b1 - _b2;
    b3  = b ^ 3;
    ts  = ( k2 / (6*b3) + k/b - b/2 ) - ( k2 / (6*a3) + k/a - a/2 );
    tp1 = -2 * n2 * (b-a) / (np2);
    tp2 = -2 * n2 * np * log(b/a) / (nm2);
    tp3 = n2 * (1/b - 1/a) / 2;
    tp4 = 16 * n4 * (n4+1) * log((2*np*b - nm2) / (2*np*a - nm2)) / (np3*nm2);
    tp5 = 16 * n6 * (1 / (2*np*b-nm2) - 1 / (2*np*a-nm2)) / (np3);
    tp  = tp1 + tp2 + tp3 + tp4 + tp5;
    tav = (ts + tp) / (2 * sa2);

    return tav
end


text:
module ServoWithRamp

using ModiaLang
using Unitful

# ModiaLang models
include("$(ModiaLang.path)/models/Blocks.jl")
include("$(ModiaLang.path)/models/Electric.jl")
include("$(ModiaLang.path)/models/Rotational.jl")

Controller = Model(
    # Interface
    refLoadAngle     = input,
    actualMotorAngle = input,
    actualMotorSpeed = input,
    refTorque        = output,

    # Components
    gainOuter     = Gain,
    angleFeedback = Feedback,
    gainInner     = Gain,
    speedFeedback = Feedback,
    PI = PI,

    connect = :[
        (refLoadAngle, gainOuter.u)
        (gainOuter.y, angleFeedback.u1)
        (actualMotorAngle, angleFeedback.u2)
        (angleFeedback.y, gainInner.u)
        (gainInner.y, speedFeedback.u1)
        (actualMotorSpeed, speedFeedback.u2)
        (speedFeedback.y, PI.u)
        (PI.y, refTorque)
    ]
)

Drive = Model(
    # Interface
    refTorque   = input,
    actualSpeed = output,
    actualAngle = output,
    flange      = Flange,
    # Components
    torque       = UnitlessTorque,
    speedSensor  = UnitlessSpeedSensor,
    angleSensor  = UnitlessAngleSensor,
    motorInertia = Inertia,
    idealGear    = IdealGear,

    connect = :[
        (refTorque, torque.tau)
        (torque.flange, speedSensor.flange, angleSensor.flange, motorInertia.flange_a)
        (motorInertia.flange_b, idealGear.flange_a)
        (idealGear.flange_b, flange)
        (speedSensor.w, actualSpeed)
        (angleSensor.phi, actualAngle)
    ]
)

Servo = Model(
    #gearRatio = 1.0,
    refLoadAngle = input,
    flange = Flange,
    ppi   = Controller, # | Map(gearRatio=:gearRatio),
    drive = Drive, # | Map(gearRatio=:gearRatio),

    connect = :[
        # inputs of ppi
        (refLoadAngle, ppi.refLoadAngle)
        (ppi.actualMotorSpeed, drive.actualSpeed)
        (ppi.actualMotorAngle, drive.actualAngle)
        # output of ppi --> input of motor
        (ppi.refTorque, drive.refTorque)
        # Flange of drive
        (drive.flange, flange)
    ]
)

servoParameters = Map(gearRato = 156.0,
                      ppi = Map(
                          gainOuter = Map(k=156.0),
                          gainInner = Map(k=50.0),
                          PI        = Map(k=0.1, T=1.0)
                      ),
                      drive = Map(
                          motorInertia = Map(J= 0.0000135 + 0.000000409),
                          idealGear = Map(ratio=156.0)
                      )
                  )

TestServo = Model(
    ramp  = Ramp | Map(duration=1.18u"s", height=2.95),
    servo = Servo | servoParameters,
    load  = Inertia | Map(J=0.32u"kg*m^2"),
    equations =:[load.flange_b.tau = 0u"N*m"],
    connect = :[
        (ramp.y, servo.refLoadAngle)
        (servo.flange, load.flange_a)
    ]
)

testServo = @instantiateModel(TestServo, unitless=true, logCode=false, log=false)

stopTime = 4.0
tolerance = 1e-6
requiredFinalStates = [-0.004445925511355035, 2.9499994183801284, 5.816286501116783e-7]
simulate!(testServo, stopTime=stopTime, tolerance=tolerance, log=true,logStates=true, requiredFinalStates=requiredFinalStates)

@usingModiaPlot
plotVariables = [("ramp.y", "load.phi"); "load.w"; "servo.ppi.PI.x"; "servo.ppi.refTorque"]
plot(testServo, plotVariables, figure=1)

end


text:
module StateSpaceModels

# package code goes here
include("types.jl")
include("filters.jl")

end # module


text:
import DiffEqJump
using Transducers
using Transducers: @next, complete, R_, inner, wrapping, xform, start, next, wrap, unwrap, Unseen

struct SSAIter{F,uType,tType,P,S,CB,SA,OPT,TS}
    integrator::DiffEqJump.SSAIntegrator{F,uType,tType,P,S,CB,SA,OPT,TS}
end

Base.IteratorSize(::Type{SSAIter{F,uType,tType,P,S,CB,SA,OPT,TS}}) where {F,uType,tType,P,S,CB,SA,OPT,TS} = Base.SizeUnknown()
Base.eltype(::Type{SSAIter{F,uType,tType,P,S,CB,SA,OPT,TS}}) where {F,uType,tType,P,S,CB,SA,OPT,TS} = Tuple{uType,tType,Int}

function Base.iterate(iter::SSAIter)
    integrator = iter.integrator
    end_time = integrator.sol.prob.tspan[2]
    if integrator.t >= end_time
        return nothing
    end
    aggregator = integrator.cb.condition
    next_jump_time = integrator.tstop > integrator.t ? integrator.tstop : typemax(integrator.tstop)
    if !isempty(integrator.tstops) &&
            integrator.tstops_idx <= length(integrator.tstops) &&
            integrator.tstops[integrator.tstops_idx] < next_jump_time
        next_jump_time = integrator.tstops[integrator.tstops_idx]
    end
    t = min(next_jump_time, end_time)
    i = t == end_time ? 0 : aggregator.next_jump
    return (integrator.u, t, i), ()
end

function Base.iterate(iter::SSAIter, state::Tuple{})
    integrator = iter.integrator
    if DiffEqJump.should_continue_solve(integrator)
        end_time = integrator.sol.prob.tspan[2]
        step!(integrator)
        aggregator = integrator.cb.condition
        i = aggregator.next_jump

        next_jump_time = integrator.tstop > integrator.t ? integrator.tstop : typemax(integrator.tstop)
        if !isempty(integrator.tstops) &&
                integrator.tstops_idx <= length(integrator.tstops) &&
                integrator.tstops[integrator.tstops_idx] < next_jump_time
            next_jump_time = integrator.tstops[integrator.tstops_idx]
            i = 0
        end
        t = min(next_jump_time, end_time)
        if t == end_time
            i = 0
        end

        return (integrator.u, t, i), ()
    end
    nothing
end

function sub_trajectory(traj, indices)
    traj |> Map((u,t,i)::Tuple -> ((@view u[indices]), t, i)) |> Thin() |> collect_trajectory
end
struct Chain{V,T <: AbstractVector,U <: AbstractVector} <: AbstractVector{V}
    head::T
    tail::U

    function Chain(head::AbstractVector{X}, tail::AbstractVector{Y}) where {X, Y}
        new{promote_type(X, Y), typeof(head), typeof(tail)}(head, tail)
    end
end

Base.IndexStyle(::Type{<:Chain}) = IndexLinear()
Base.size(ch::Chain) = size(ch.head) .+ size(ch.tail)
Base.getindex(ch::Chain{V}, i::Int) where {V} = convert(V, i > length(ch.head) ? ch.tail[i - length(ch.head)] : ch.head[i])
Base.setindex!(ch::Chain, v, i::Int) = i > length(ch.head) ? ch.tail[i - length(ch.head)] = v : ch.head[i] = v

# merge more than 2 trajectories using recursion
merge_trajectories(traj) = traj
function merge_trajectories(traj1, traj2)
    traj1 |> MergeWith(traj2)
end
function merge_trajectories(traj1, traj2, other_trajs...)
    merge12 = merge_trajectories(traj1, traj2)
    merge_trajectories(merge12, other_trajs...)
end

struct Thin <: Transducers.AbstractFilter end

Transducers.start(rf::R_{Thin}, result) = wrap(rf, Unseen(), start(inner(rf), result))

function Transducers.complete(rf::R_{Thin}, result)
    result = wrapping(rf, result) do prev, iresult
        if prev isa Unseen
            return prev, iresult
        end
        return prev, next(inner(rf), iresult, prev)
    end
    complete(inner(rf), unwrap(rf, result)[2])
end

@inline Transducers.next(rf::R_{Thin}, result, (u, t, i)) =
    wrapping(rf, result) do prev, iresult
        if prev isa Unseen
            return (copy(u), t, i), iresult
        end
        if prev[1] != u
            return (copy(u), t, i), next(inner(rf), iresult, prev)
        else
            return (prev[1], t, i), iresult
        end
    end

struct MergeWith{uType,tType} <: Transducer
    traj::Trajectory{uType,tType}
    start_index::Int
end

MergeWith(traj::Trajectory) = MergeWith(traj, 1)

function Transducers.start(rf::R_{MergeWith}, result)
    private_state = xform(rf).start_index
    return wrap(rf, private_state, start(inner(rf), result))
end

function Transducers.next(rf::R_{MergeWith}, result, (u, t, i))
    f = let u=u, t=t, i=i, rf=rf
        function (index::Int, iresult)
            merge_traj = xform(rf).traj
            ri = i
            @inbounds while index <= length(merge_traj.t) && merge_traj.t[index] < t
                merge_i = merge_traj.i[index]
                iresult = next(inner(rf), iresult, (Chain(u, merge_traj.u[index]), merge_traj.t[index], merge_i))
                index += 1
            end

            if index > length(merge_traj.t)
                return index, reduced(iresult)
            end

            iresult = next(inner(rf), iresult, (Chain(u, (@inbounds merge_traj.u[index])), t, ri))

            @inbounds if index <= length(merge_traj.t) && merge_traj.t[index] == t
                index += 1
            end

            index, iresult
        end
    end

    wrapping(f, rf, result)
end

function Transducers.complete(rf::R_{MergeWith}, result)
    _private_state, inner_result = unwrap(rf, result)
    return complete(inner(rf), inner_result)
end

function collect_trajectory(xf::Transducers.Transducer, itr)
    rf = Transducers.ProductRF(Map(copy)'(Transducers.push!!), Transducers.push!!, Transducers.push!!)
    (u, t, i) = foldxl(rf, xf, itr; init=(Transducers.Empty(Vector), Transducers.Empty(Vector), Transducers.Empty(Vector)))
    Trajectory(u, t, i)
end

collect_trajectory(ed::Transducers.Eduction) = collect_trajectory(Transducers.extract_transducer(ed)...)


text:
"""
A [`User`](@ref)'s voice connection status.
More details [here](https://discordapp.com/developers/docs/resources/voice#voice-state-object).
"""
struct VoiceState
    guild_id::Union{Snowflake, Missing}
    channel_id::Union{Snowflake, Nothing}
    user_id::Snowflake
    member::Union{Member, Missing}
    session_id::String
    deaf::Bool
    mute::Bool
    self_deaf::Bool
    self_mute::Bool
    suppress::Bool
end
@boilerplate VoiceState :dict :docs :lower :merge


text:
# # Basic Usage

using Distributions
using HMMBase
using PyPlot
using Seaborn

rc("axes", xmargin = 0) # hide
set_style("whitegrid")  # hide

# ### Model Specification

# B can contains any probability distribution from the `Distributions` package

a = [0.6, 0.4]
A = [0.9 0.1; 0.1 0.9]
B = [MvNormal([0.0, 5.0], ones(2)*1), MvNormal([0.0, 5.0], ones(2)*3)]
hmm = HMM(a, A, B)
size(hmm) # (number of states, observations dimension)

# ### Sampling

z, y = rand(hmm, 500, seq = true)

# Let's plot the observations and the hidden state sequence:

_, axes = subplots(nrows = 2, figsize = (9,3))
axes[1].plot(y)
axes[2].plot(z, linestyle = :steps)
gcf() # hide

# We can also drop the time dimension and plot the data in the plane:

_, axes = subplots(ncols = 2, figsize = (9,3))
axes[1].scatter(y[:,1], y[:,2], s = 3.0)
axes[2].scatter(y[:,1], y[:,2], s = 3.0, c = z, cmap = "tab10")
axes[1].set_title("Observations")
axes[2].set_title("Observations and hidden states")
gcf() # hide

# ### Inference

α, logtot = forward(hmm, y)
β, logtot = backward(hmm, y)

γ = posteriors(hmm, y) # or
γ = posteriors(α, β)

size(α), size(β), size(γ)
#-

figure(figsize = (9,2)) # hide
plot([α[:,1] β[:,1] γ[:,1]])
legend(["Forward", "Backward", "Posteriors"], loc = "upper right")
gcf() # hide
#-

_, axes = subplots(ncols = 3, figsize = (9,3))
for (ax, probs, title) in zip(axes, [α, β, γ], ["Forward", "Backward", "Posteriors"])
    ax.scatter(y[:,1], y[:,2], s = 3.0, c = probs[:,1], cmap = "Reds")
    ax.set_title(title)
end
gcf() # hide
#-

z_map = [z.I[2] for z in argmax(γ, dims = 2)][:]
z_viterbi = viterbi(hmm, y)

figure(figsize = (9,2)) # hide
plot([z z_map z_viterbi])
legend(["True sequence", "MAP", "Viterbi"], loc = "upper right")
gcf() # hide
#-

_, axes = subplots(ncols = 2, figsize = (9,3))
for (ax, seq, title) in zip(axes, [z_map, z_viterbi], ["MAP", "Viterbi"])
    ax.scatter(y[:,1], y[:,2], s = 3.0, c = seq, cmap = "Reds_r")
    ax.set_title(title)
end
gcf() # hide

# ### Parameters Estimation

hmm = HMM(randtransmat(2), [MvNormal(rand(2), ones(2)), MvNormal(rand(2), ones(2))])
hmm, hist = fit_mle(hmm, y, display = :iter, init = :kmeans)
hmm
#-

figure(figsize = (4,3)) # hide
plot(hist.logtots)
gcf() # hide

max_stars_repo_name: Shu-Lea-Lai/module-MultiHMMBase


text:
using Test

@testset "MOI_wrapper" begin
    include("MOI_wrapper.jl")
end


text:
module adapters_odbc_postgresql_repo_test

using Test # @test
using Octo.Adapters.ODBC # Repo Schema

Repo.debug_sql()

Repo.connect(
    adapter  = Octo.Adapters.ODBC,
    database = Octo.DBMS.SQL,
    dsn      = "PgSQL-test",
    username = "postgres",
    password = "",
)

struct Employee
end
Schema.model(Employee, table_name="Employee", primary_key="ID")

Repo.execute([DROP TABLE IF EXISTS Employee])
Repo.execute(Raw("""
    CREATE TABLE Employee (
        ID SERIAL,
        Name VARCHAR(255),
        Salary FLOAT(8),
        PRIMARY KEY (ID)
    )"""))

changes = (Name="John", Salary=10000.50)
Repo.insert!(Employee, changes)
changes = (Name="Cloris", Salary=20000.50)
Repo.insert!(Employee, changes)
df = Repo.query(Employee)
@test size(df) == (2,)
@test df == [(id=1, name="John", salary=10000.50),
             (id=2, name="Cloris", salary=20000.50)]

Repo.disconnect()

end # module adapters_odbc_postgresql_repo_test


text:
function primitiven{T<:Integer}(m::T)
    1 < m || return T[]
    m != 2 || return T[1]
    !isprime(m) || return T[2:2:m-1]
    rp = trues(m-1)
    if isodd(m)
        rp[1:2:m-1] = false
    end
    for p in keys(factor(m))
        rp[p:p:m-1] = false
    end
    T[1:m-1][rp]
end

function pythagoreantripcount{T<:Integer}(plim::T)
    primcnt = 0
    fullcnt = 0
    11 < plim || return (primcnt, fullcnt)
    for m in 2:plim
        p = 2m^2
        p+2m <= plim || break
        for n in primitiven(m)
            q = p + 2m*n
            q <= plim || break
            primcnt += 1
            fullcnt += div(plim, q)
        end
    end
    return (primcnt, fullcnt)
end

println("Counting Pythagorian Triplets within perimeter limits:")
println("    Limit          All   Primitive")
for om in 1:10
    (pcnt, fcnt) = pythagoreantripcount(10^om)
    println(@sprintf "    10^%02d  %11d   %9d" om fcnt pcnt)
end

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
module Test_ode_save_func_callback

#=
StackOverFlowError with two usages of save_func #548

https://github.com/JuliaDiffEq/DifferentialEquations.jl/issues/548

Issue is fixed by setting:

   ArrayInterface.ismutable(::Type{Model}) = true

=#

using DifferentialEquations, StaticArrays

mutable struct Model{R<:Real}
    m::R  # mass
    L::R  # length
    d::R  # damping
    g::R  # gravity
    ϕ::R  # angle
    ω::R  # der(angle)
    α::R  # der(der(angle))
    r::SVector{2,R}   # absolute position of end point

    function Model{R}(m::R, L::R, d::R, g::R) where {R<:Real}
        new(m, L, d, g, 0, 0, 0, SVector{2,R}(0,0))
    end
end

# Due to StackOverFlowError https://github.com/JuliaDiffEq/DifferentialEquations.jl/issues/548
using ArrayInterface
ArrayInterface.ismutable(::Type{Model}) = true

Model(;m::R=0, L::R=0, d::R=0, g::R =0) where {R<:Real} = Model{R}(m,L,d,g)

"simplePendulum: Compute m::Model at actual time instant"
function simplePendulum!(m::Model, t, x)::Nothing
    m.ϕ = x[1]
    m.ω = x[2]
    m.α = (-m.m*m.g*m.L*sin(m.ϕ) - m.d*m.ω) / (m.m*m.L^2)
    m.r = @SVector [m.L*sin(m.ϕ), -m.L*cos(m.ϕ)]
    return nothing
end

"SimplePendulum1: Called by integrator"
function simplePendulum1!(der_x, x, m::Model, t)
    simplePendulum!(m, t, x)
    der_x[1] = m.ω
    der_x[2] = m.α
end

"SimplePendulum2: Called at every communication point"
function simplePendulum2(x, t, integrator)
    m::Model = integrator.p
    simplePendulum!(m, t, x)
    return deepcopy(m)
end

t_start = 0.0
t_inc   = 0.01
t_end   = 7.0
tspan   = (t_start, t_end)
tspan2  = t_start:t_inc:t_end

model1 = Model(m=1.0, L=1.0, d=0.2, g=9.81)
x1₀ = [π/2, 0.0]

saved_values1=SavedValues(typeof(x1₀[1]), Model)
cb1 = SavingCallback(simplePendulum2, saved_values1, saveat=tspan2)

prob1 = ODEProblem(simplePendulum1!, x1₀, tspan, model1)
sol1 = solve(prob1, Tsit5(), reltol = 1e-6, saveat=t_inc, callback=cb1)

# Triggers stack-overflow error if ArrayInterface.ismutable(::Type{Model}) = true is not set
sol2 = solve(prob1, Tsit5(), reltol = 1e-6, saveat=t_inc, callback=cb1)

end


text:
# Julia wrapper for header: deps/RustDylib/include/myheader.h
# Automatically generated using Clang.jl wrap_c, version 0.0.0

function rustdylib_abs_f32(i::Cfloat)
    ccall((:rustdylib_abs_f32, librustdylib), Cfloat, (Cfloat,), i)
end

function rustdylib_abs_f64(i::Cdouble)
    ccall((:rustdylib_abs_f64, librustdylib), Cdouble, (Cdouble,), i)
end

function rustdylib_abs_i32(i::Int32)
    ccall((:rustdylib_abs_i32, librustdylib), Int32, (Int32,), i)
end

function rustdylib_abs_i64(i::Int64)
    ccall((:rustdylib_abs_i64, librustdylib), Int64, (Int64,), i)
end

function rustdylib_free_rust_owned_string(s::Cstring)
    ccall((:rustdylib_free_rust_owned_string, librustdylib), Cvoid, (Cstring,), s)
end

function rustdylib_generate_rust_owned_string()
    ccall((:rustdylib_generate_rust_owned_string, librustdylib), Cstring, ())
end

function rustdylib_inspect_string(cstring)
    ccall((:rustdylib_inspect_string, librustdylib), Cvoid, (Cstring,), cstring)
end

function rustdylib_is_true_bool(b::Bool)
    ccall((:rustdylib_is_true_bool, librustdylib), Bool, (Bool,), b)
end

function rustdylib_printhello()
    ccall((:rustdylib_printhello, librustdylib), Cvoid, ())
end


text:
# Test constraint deletion
@testset "JuMP.delete (Constraints)" begin
    # initialize model and references
    m = InfiniteModel()
    @infinite_parameter(m, 0 <= par <= 1)
    @infinite_variable(m, inf(par))
    @point_variable(m, inf(0.5), pt)
    @hold_variable(m, x)
    data = DiscreteMeasureData(par, [1], [1])
    meas = measure(inf + par - x, data)
    rv = ReducedInfiniteVariableRef(m, 42)
    @constraint(m, cref, par - inf + pt + 2x - rv + meas <= 1)
    # test normal deletion
    @test isa(delete(m, cref), Nothing)
    @test !is_valid(m, cref)
    @test !haskey(m.param_to_constrs, JuMP.index(par))
    @test !haskey(m.var_to_constrs, JuMP.index(inf))
    @test !haskey(m.var_to_constrs, JuMP.index(pt))
    @test !haskey(m.var_to_constrs, JuMP.index(x))
    @test !haskey(m.reduced_to_constrs, JuMP.index(rv))
    @test !haskey(m.meas_to_constrs, JuMP.index(meas))
    @test !haskey(m.constrs, JuMP.index(cref))
    @test !haskey(m.constr_to_name, JuMP.index(cref))
    @test !haskey(m.constr_in_var_info, JuMP.index(cref))
    # test assertion error
    @test_throws AssertionError delete(m, cref)
end

# Test deleting variable information
@testset "Variable Information" begin
    # initialize model and references
    m = InfiniteModel()
    @infinite_parameter(m, 0 <= par <= 1)
    @infinite_variable(m, 0 <= inf1(par) <= 1, Bin)
    @infinite_variable(m, inf2(par) == 1, Int)
    @point_variable(m, inf1(0.5), pt1)
    @point_variable(m, inf2(0.5), pt2)
    @hold_variable(m, 0 <= gb1 <= 1, Bin)
    @hold_variable(m, gb2 == 1, Int)
    # test delete_lower_bound
    @testset "JuMP.delete_lower_bound" begin
        # test with infinite variable
        @test has_lower_bound(inf1)
        @test isa(delete_lower_bound(inf1), Nothing)
        @test !has_lower_bound(inf1)
        @test_throws ErrorException delete_lower_bound(inf2)
        # test with point variable
        @test has_lower_bound(pt1)
        @test isa(delete_lower_bound(pt1), Nothing)
        @test !has_lower_bound(pt1)
        @test_throws ErrorException delete_lower_bound(pt2)
        # test with hold variable
        @test has_lower_bound(gb1)
        @test isa(delete_lower_bound(gb1), Nothing)
        @test !has_lower_bound(gb1)
        @test_throws ErrorException delete_lower_bound(gb2)
    end
    # test delete_upper_bound
    @testset "JuMP.delete_upper_bound" begin
        # test with infinite variable
        @test has_upper_bound(inf1)
        @test isa(delete_upper_bound(inf1), Nothing)
        @test !has_upper_bound(inf1)
        @test_throws ErrorException delete_upper_bound(inf2)
        # test with point variable
        @test has_upper_bound(pt1)
        @test isa(delete_upper_bound(pt1), Nothing)
        @test !has_upper_bound(pt1)
        @test_throws ErrorException delete_upper_bound(pt2)
        # test with hold variable
        @test has_upper_bound(gb1)
        @test isa(delete_upper_bound(gb1), Nothing)
        @test !has_upper_bound(gb1)
        @test_throws ErrorException delete_upper_bound(gb2)
    end
    # test unfix
    @testset "JuMP.unfix" begin
        # test with infinite variable
        @test is_fixed(inf2)
        @test isa(unfix(inf2), Nothing)
        @test !is_fixed(inf2)
        @test_throws ErrorException unfix(inf1)
        # test with point variable
        @test is_fixed(pt2)
        @test isa(unfix(pt2), Nothing)
        @test !is_fixed(pt2)
        @test_throws ErrorException unfix(pt1)
        # test with hold variable
        @test is_fixed(gb2)
        @test isa(unfix(gb2), Nothing)
        @test !is_fixed(gb2)
        @test_throws ErrorException unfix(gb1)
    end
    # test unset_binary
    @testset "JuMP.unset_binary" begin
        # test with infinite variable
        @test is_binary(inf1)
        @test isa(unset_binary(inf1), Nothing)
        @test !is_binary(inf1)
        @test_throws ErrorException unset_binary(inf2)
        # test with point variable
        @test is_binary(pt1)
        @test isa(unset_binary(pt1), Nothing)
        @test !is_binary(pt1)
        @test_throws ErrorException unset_binary(pt2)
        # test with hold variable
        @test is_binary(gb1)
        @test isa(unset_binary(gb1), Nothing)
        @test !is_binary(gb1)
        @test_throws ErrorException unset_binary(gb2)
    end
    # test unset_integer
    @testset "JuMP.unset_integer" begin
        # test with infinite variable
        @test is_integer(inf2)
        @test isa(unset_integer(inf2), Nothing)
        @test !is_integer(inf2)
        @test_throws ErrorException unset_integer(inf1)
        # test with point variable
        @test is_integer(pt2)
        @test isa(unset_integer(pt2), Nothing)
        @test !is_integer(pt2)
        @test_throws ErrorException unset_integer(pt1)
        # test with hold variable
        @test is_integer(gb2)
        @test isa(unset_integer(gb2), Nothing)
        @test !is_integer(gb2)
        @test_throws ErrorException unset_integer(gb1)
    end
end

# Test parameter deletion
@testset "Parameters" begin
    # initialize model and references
    m = InfiniteModel()
    @infinite_parameter(m, 0 <= par <= 1)
    @infinite_parameter(m, 0 <= pars[1:2] <= 1, container = SparseAxisArray)
    @infinite_parameter(m, 0 <= par2 <= 1)
    @infinite_parameter(m, 0 <= par3 <= 1)
    @infinite_variable(m, inf(par))
    @infinite_variable(m, inf2(par, par2))
    @infinite_variable(m, inf3(par, pars))
    @infinite_variable(m, inf4(par, par2, pars))
    @point_variable(m, inf(0.5), pt)
    pt2 = @point_variable(m, inf2(0.5, 0.5))
    @point_variable(m, inf3(0, [0, 0]), pt3)
    @hold_variable(m, x)
    m.reduced_info[-1] = ReducedInfiniteInfo(inf4, Dict(2 => 0.5))
    m.infinite_to_reduced[JuMP.index(inf4)] = [-1]
    rv = ReducedInfiniteVariableRef(m, -1)
    data = DiscreteMeasureData(par, [1], [1])
    data2 = DiscreteMeasureData(pars, [1], [[1, 1]])
    meas = measure(inf + par - x + rv, data)
    meas2 = measure(inf3 - x + pt2, data2)
    meas3 = measure(inf2 + par2, data)
    meas4 = add_measure(m, Measure(par3, data))
    meas5 = measure(pars[2] + par, data)
    @constraint(m, con, inf2 + inf4 - par2 <= 0)
    m.constrs[-1] = ScalarConstraint(par3, MOI.GreaterThan(0.))
    m.param_to_constrs[JuMP.index(par3)] = [-1]
    con2 = InfiniteConstraintRef(m, -1, JuMP.shape(m.constrs[-1]))
    set_name(con2, "")

    # test _check_param_in_data (scalar)
    @testset "_check_param_in_data (Scalar)" begin
        @test_throws ErrorException InfiniteOpt._check_param_in_data(par, data)
        @test isa(InfiniteOpt._check_param_in_data(par2, data), Nothing)
    end
    # test _check_param_in_data (array)
    @testset "_check_param_in_data (Array)" begin
        @test_throws ErrorException InfiniteOpt._check_param_in_data(pars[1],
                                                                     data2)
        @test isa(InfiniteOpt._check_param_in_data(par, data2), Nothing)
    end
    # test _check_param_in_data (fallback)
    @testset "_check_param_in_data (Fallback)" begin
        @test_throws ErrorException InfiniteOpt._check_param_in_data(par, BadData())
    end
    # test _contains_pref (scalar)
    @testset "_contains_pref (Scalar)" begin
        @test InfiniteOpt._contains_pref(par, par)
        @test !InfiniteOpt._contains_pref(par, par2)
    end
    # test _contains_pref (array)
    @testset "_contains_pref (Array)" begin
        @test InfiniteOpt._contains_pref(pars, pars[1])
        @test !InfiniteOpt._contains_pref(pars, par)
    end
    # test _remove_parameter
    @testset "_remove_parameter" begin
        # test removing single element
        @test InfiniteOpt._remove_parameter((par,), par) == ((), (1,))
        @test InfiniteOpt._remove_parameter((par2, par), par) == ((par2,), (2,))
        # test removing element from array
        new_pars = JuMP.Containers.SparseAxisArray(Dict((2,) => pars[2]))
        @test InfiniteOpt._remove_parameter((pars,),
                                            pars[1]) == ((new_pars,), (1, (1,)))
        @test InfiniteOpt._remove_parameter((par, pars),
                                        pars[1]) == ((par, new_pars), (2, (1,)))
        # test removing from array that becomes empty
        @test InfiniteOpt._remove_parameter((new_pars,), pars[2]) == ((), (1,))
        @test InfiniteOpt._remove_parameter((par, new_pars),
                                            pars[2]) == ((par,), (2,))
    end
    # test _update_infinite_variable
    @testset "_update_infinite_variable" begin
        # test updating when used by measure
        @test isa(InfiniteOpt._update_infinite_variable(inf3, (pars,)), Nothing)
        @test parameter_refs(inf3) == (pars,)
        @test name(inf3) == "inf3(pars)"
        @test name(meas2) == "measure(inf3(pars) - x + inf2(0.5, 0.5))"
        # Undo the changes
        @test isa(InfiniteOpt._update_infinite_variable(inf3, (par, pars)),
                  Nothing)
        @test parameter_refs(inf3) == (par, pars)
        @test name(inf3) == "inf3(par, pars)"
        @test name(meas2) == "measure(inf3(par, pars) - x + inf2(0.5, 0.5))"
    end
    # test _remove_parameter_values
    @testset "_remove_parameter_values" begin
        # test remove single element
        @test InfiniteOpt._remove_parameter_values((0, 1), (1,)) == (1,)
        @test InfiniteOpt._remove_parameter_values((0, 1), (2,)) == (0,)
        @test InfiniteOpt._remove_parameter_values(([0, 0],), (1,)) == ()
        # test remove array element
        arr = convert(JuMP.Containers.SparseAxisArray, [0, 1, 0])
        arr2 = JuMP.Containers.SparseAxisArray(Dict((1,) => 0, (2,) => 1,
                                                    (3,) => 0))
        filter!(x -> x.first != (2,), arr2.data)
        @test InfiniteOpt._remove_parameter_values((arr, 3), (1, (2,))) == (arr2, 3)
        arr3 = JuMP.Containers.SparseAxisArray(filter(x -> x.first != (1,),
                                                      arr2.data))
        @test InfiniteOpt._remove_parameter_values((arr2,), (1, (1,))) == (arr3,)
    end
    # test _update_point_variable
    @testset "_update_point_variable" begin
        # test point variable with single parameter
        @test isa(InfiniteOpt._update_point_variable(pt, ()), Nothing)
        @test parameter_values(pt) == ()
        @test name(pt) == "pt"
        # Undo changes
        @test isa(InfiniteOpt._update_point_variable(pt, (0.5,)), Nothing)
        @test parameter_values(pt) == (0.5,)
        @test name(pt) == "pt"
        # test point variable with multiple parameters
        @test isa(InfiniteOpt._update_point_variable(pt2, (0.5,)), Nothing)
        @test parameter_values(pt2) == (0.5,)
        @test name(pt2) == "inf2(0.5)"
        @test name(meas2) == "measure(inf3(par, pars) - x + inf2(0.5))"
        # Undo changes
        @test isa(InfiniteOpt._update_point_variable(pt2, (0.5, 0.5)), Nothing)
        @test parameter_values(pt2) == (0.5, 0.5)
        @test name(pt2) == "inf2(0.5, 0.5)"
        @test name(meas2) == "measure(inf3(par, pars) - x + inf2(0.5, 0.5))"
    end
    # test _update_reduced_variable
    @testset "_update_reduced_variable" begin
        # test removing single parameter that is not reduced
        @test isa(InfiniteOpt._update_infinite_variable(inf4, (par2, pars)),
                                                        Nothing)
        @test isa(InfiniteOpt._update_reduced_variable(rv, (1, )), Nothing)
        @test m.reduced_info[JuMP.index(rv)].infinite_variable_ref == inf4
        @test m.reduced_info[JuMP.index(rv)].eval_supports == Dict(1 => 0.5)
        @test name(rv) == "inf4(0.5, pars)"
        @test name(meas) == "measure(inf(par) + par - x + inf4(0.5, pars))"
        # Undo changes
        m.reduced_info[-1] = ReducedInfiniteInfo(inf4, Dict(2 => 0.5))
        @test isa(InfiniteOpt._update_infinite_variable(inf4, (par, par2, pars)),
                                                        Nothing)
        # test removing single parameter that is reduced
        @test isa(InfiniteOpt._update_infinite_variable(inf4, (par, pars)),
                                                        Nothing)
        @test isa(InfiniteOpt._update_reduced_variable(rv, (2, )), Nothing)
        @test m.reduced_info[JuMP.index(rv)].infinite_variable_ref == inf4
        @test m.reduced_info[JuMP.index(rv)].eval_supports == Dict{Int,
                               Union{Number, JuMP.Containers.SparseAxisArray}}()
        @test name(rv) == "inf4(par, pars)"
        @test name(meas) == "measure(inf(par) + par - x + inf4(par, pars))"
        # Undo changes
        @test isa(InfiniteOpt._update_infinite_variable(inf4, (par, par2, pars)),
                                                        Nothing)
        # test removing a different single parameter
        m.reduced_info[-1] = ReducedInfiniteInfo(inf4, Dict(1 => 0.5, 2 => 0.5))
        @test isa(InfiniteOpt._update_infinite_variable(inf4, (par, pars)),
                                                        Nothing)
        @test isa(InfiniteOpt._update_reduced_variable(rv, (2, )), Nothing)
        @test m.reduced_info[JuMP.index(rv)].infinite_variable_ref == inf4
        @test m.reduced_info[JuMP.index(rv)].eval_supports == Dict(1 => 0.5)
        @test name(rv) == "inf4(0.5, pars)"
        @test name(meas) == "measure(inf(par) + par - x + inf4(0.5, pars))"
        # Undo changes
        m.reduced_info[-1] = ReducedInfiniteInfo(inf4, Dict(2 => 0.5))
        @test isa(InfiniteOpt._update_infinite_variable(inf4, (par, par2, pars)),
                                                        Nothing)
        # prepare for removing array element
        supp = convert(JuMP.Containers.SparseAxisArray, [0, 0])
        m.reduced_info[-1] = ReducedInfiniteInfo(inf4, Dict(2 => 0.5, 3 => supp))
        new_supp = JuMP.Containers.SparseAxisArray(Dict((1,) => 0, (2,) => 0))
        filter!(x -> x.first != (1,), new_supp.data)
        # test removing array element
        @test isa(InfiniteOpt._update_reduced_variable(rv, (3, (1,))), Nothing)
        @test m.reduced_info[JuMP.index(rv)].infinite_variable_ref == inf4
        @test m.reduced_info[JuMP.index(rv)].eval_supports == Dict(2 => 0.5, 3 => new_supp)
        @test name(rv) == "inf4(par, 0.5, " * string(new_supp) * ")"
        @test name(meas) == "measure(inf(par) + par - x + " * name(rv) * ")"
        # Undo changes
        m.reduced_info[-1] = ReducedInfiniteInfo(inf4, Dict(2 => 0.5))
        @test isa(InfiniteOpt._update_infinite_variable(inf4, (par, par2, pars)),
                                                        Nothing)
    end
    # test JuMP.delete for parameters
    @testset "JuMP.delete (Parameters)" begin
        # test normal usage
        idx = JuMP.index(par2)
        @test isa(delete(m, par2), Nothing)
        @test name(meas3) == "measure(inf2(par))"
        @test !haskey(m.param_to_meas, idx)
        @test parameter_refs(inf2) == (par,)
        @test parameter_refs(inf4) == (par, pars)
        @test parameter_values(pt2) == (0.5,)
        @test parameter_refs(rv) == (par, pars)
        @test name(inf2) == "inf2(par)"
        @test name(inf4) == "inf4(par, pars)"
        @test name(pt2) == "inf2(0.5)"
        @test name(rv) == "inf4(par, pars)"
        @test !haskey(m.param_to_vars, idx)
        @test string(m.constrs[JuMP.index(con)].func) == "inf2(par) + " *
                                                         "inf4(par, pars)"
        @test !haskey(m.param_to_constrs, idx)
        @test !haskey(m.params, idx)
        @test !haskey(m.param_to_name, idx)
        @test !haskey(m.param_to_group_id, idx)
        # test special measure case with single parameter (possible through
        # multiple deletions) and with single parameter in constraint
        idx = JuMP.index(par3)
        @test isa(delete(m, par3), Nothing)
        @test name(meas4) == "measure(0)"
        @test !haskey(m.param_to_meas, idx)
        @test string(m.constrs[JuMP.index(con2)].func) == "0"
        @test !haskey(m.param_to_constrs, idx)
        @test !haskey(m.params, idx)
        @test !haskey(m.param_to_name, idx)
        @test !haskey(m.param_to_group_id, idx)
        # test invalid parameter
        @test_throws AssertionError delete(m, par2)
        @test_throws AssertionError delete(m, par3)
        # test measure data check
        @test_throws ErrorException delete(m, par)
        @test_throws ErrorException delete(m, pars[1])
        # test array element deletion
        delete!(m.measures, JuMP.index(meas2))
        m.param_to_meas[JuMP.index(pars[2])] = [JuMP.index(meas5)]
        delete!(m.var_to_meas, JuMP.index(inf3))
        idx = JuMP.index(pars[2])
        @test isa(delete(m, pars[2]), Nothing)
        @test name(meas5) == "measure(par)"
        @test !haskey(m.param_to_meas, idx)
        new_pars = convert(JuMP.Containers.SparseAxisArray, [pars[1]])
        @test parameter_refs(inf3) == (par, new_pars)
        @test parameter_refs(inf4) == (par, new_pars)
        new_vals = convert(JuMP.Containers.SparseAxisArray, [0])
        @test parameter_values(pt3) == (0, new_vals)
        @test parameter_refs(rv) == (par, new_pars)
        @test name(inf3) == "inf3(par, pars)"
        @test name(inf4) == "inf4(par, pars)"
        @test name(pt3) == "pt3"
        @test name(rv) == "inf4(par, pars)"
        @test !haskey(m.param_to_vars, idx)
        @test string(m.constrs[JuMP.index(con)].func) == "inf2(par) + " *
                                                         "inf4(par, pars)"
        @test !haskey(m.param_to_constrs, idx)
        @test !haskey(m.params, idx)
        @test !haskey(m.param_to_name, idx)
        @test !haskey(m.param_to_group_id, idx)
    end
 end

 # Test reduced variable deletion
 @testset "JuMP.delete (Reduced Variables)" begin
     # intialize the model
     m = InfiniteModel()
     @infinite_parameter(m, 0 <= par <= 1)
     @infinite_parameter(m, 0 <= par2 <= 1)
     @infinite_variable(m, inf(par, par2))
     @point_variable(m, inf(0.5, 0.5), pt)
     @hold_variable(m, x)
     m.reduced_info[-1] = ReducedInfiniteInfo(inf, Dict(2 => 0.5))
     m.infinite_to_reduced[JuMP.index(inf)] = [-1]
     rv = ReducedInfiniteVariableRef(m, -1)
     m.reduced_info[-2] = ReducedInfiniteInfo(inf, Dict(2 => 0.5))
     m.infinite_to_reduced[JuMP.index(inf)] = [-2]
     rv2 = ReducedInfiniteVariableRef(m, -2)
     data = DiscreteMeasureData(par, [1], [1])
     meas = measure(inf + par - x + rv, data)
     meas2 = measure(rv2, data)
     @constraint(m, con, x + rv <= 0)
     m.constrs[-1] = ScalarConstraint(rv2, MOI.GreaterThan(0.))
     m.reduced_to_constrs[JuMP.index(rv2)] = [-1]
     con2 = InfiniteConstraintRef(m, -1, JuMP.shape(m.constrs[-1]))
     set_name(con2, "")
     # test normal deletion
     @test isa(delete(m, rv), Nothing)
     @test name(meas) == "measure(inf(par, par2) + par - x)"
     @test !haskey(m.reduced_to_meas, JuMP.index(rv))
     @test string(m.constrs[JuMP.index(con)].func) == "x"
     @test !haskey(m.reduced_to_constrs, JuMP.index(rv))
     @test m.infinite_to_reduced[JuMP.index(inf)] == [JuMP.index(rv2)]
     @test !haskey(m.reduced_info, JuMP.index(rv))
     # test deletion of special cases
     @test isa(delete(m, rv2), Nothing)
     @test name(meas2) == "measure(0)"
     @test !haskey(m.reduced_to_meas, JuMP.index(rv2))
     @test string(m.constrs[JuMP.index(con2)].func) == "0"
     @test !haskey(m.reduced_to_constrs, JuMP.index(rv2))
     @test !haskey(m.infinite_to_reduced, JuMP.index(inf))
     @test !haskey(m.reduced_info, JuMP.index(rv2))
     # test error
     @test_throws AssertionError delete(m, rv)
     @test_throws AssertionError delete(m, rv2)
 end

 # Test variable deletion
 @testset "JuMP.delete (Hold Variables)" begin
     # intialize the model
     m = InfiniteModel()
     @infinite_parameter(m, 0 <= par <= 1)
     @hold_variable(m, 0 <= x <= 1, Bin)
     @hold_variable(m, y == 1, Int)
     data = DiscreteMeasureData(par, [1], [1])
     meas1 = measure(x + y + par, data)
     meas2 = add_measure(m, Measure(y, data))
     @constraint(m, con1, x + y + par <= 0)
     con2 = add_constraint(m, ScalarConstraint(y, MOI.LessThan(0.)))
     @objective(m, Min, x + y)
     # test deletion of x
     @test isa(delete(m, x), Nothing)
     @test !haskey(m.var_to_lower_bound, JuMP.index(x))
     @test !haskey(m.var_to_upper_bound, JuMP.index(x))
     @test !haskey(m.var_to_zero_one, JuMP.index(x))
     @test name(meas1) == "measure(y + par)"
     @test !haskey(m.var_to_meas, JuMP.index(x))
     @test string(m.constrs[JuMP.index(con1)].func) == "y + par"
     @test !haskey(m.var_to_constrs, JuMP.index(x))
     @test string(objective_function(m)) == "y"
     @test !haskey(m.var_in_objective, JuMP.index(x))
     @test !haskey(m.vars, JuMP.index(x))
     @test !haskey(m.var_to_name, JuMP.index(x))
     # test deletion of y
     set_objective_function(m, y)
     @test isa(delete(m, y), Nothing)
     @test !haskey(m.var_to_fix, JuMP.index(y))
     @test !haskey(m.var_to_integrality, JuMP.index(y))
     @test name(meas1) == "measure(par)"
     @test name(meas2) == "measure(0)"
     @test !haskey(m.var_to_meas, JuMP.index(y))
     @test string(m.constrs[JuMP.index(con1)].func) == "par"
     @test string(m.constrs[JuMP.index(con2)].func) == "0"
     @test !haskey(m.var_to_constrs, JuMP.index(y))
     @test string(objective_function(m)) == "0"
     @test !haskey(m.var_in_objective, JuMP.index(y))
     @test !haskey(m.vars, JuMP.index(y))
     @test !haskey(m.var_to_name, JuMP.index(y))
     # test errors
     @test_throws AssertionError delete(m, x)
     @test_throws AssertionError delete(m, y)
 end

 # Test variable deletion
 @testset "JuMP.delete (Point Variables)" begin
     # intialize the model
     m = InfiniteModel()
     @infinite_parameter(m, 0 <= par <= 1)
     @infinite_variable(m, inf(par))
     @point_variable(m, inf(0), 0 <= x <= 1, Bin)
     @point_variable(m, inf(1), y == 1, Int)
     data = DiscreteMeasureData(par, [1], [1])
     meas1 = measure(x + y + par, data)
     meas2 = add_measure(m, Measure(y, data))
     @constraint(m, con1, x + y + par <= 0)
     con2 = add_constraint(m, ScalarConstraint(y, MOI.LessThan(0.)))
     @objective(m, Min, x + y)
     # test deletion of x
     @test isa(delete(m, x), Nothing)
     @test !haskey(m.var_to_lower_bound, JuMP.index(x))
     @test !haskey(m.var_to_upper_bound, JuMP.index(x))
     @test !haskey(m.var_to_zero_one, JuMP.index(x))
     @test name(meas1) == "measure(y + par)"
     @test !haskey(m.var_to_meas, JuMP.index(x))
     @test string(m.constrs[JuMP.index(con1)].func) == "y + par"
     @test !haskey(m.var_to_constrs, JuMP.index(x))
     @test string(objective_function(m)) == "y"
     @test m.infinite_to_points[JuMP.index(inf)] == [JuMP.index(y)]
     @test !haskey(m.var_in_objective, JuMP.index(x))
     @test !haskey(m.vars, JuMP.index(x))
     @test !haskey(m.var_to_name, JuMP.index(x))
     # test deletion of y
     set_objective_function(m, y)
     @test isa(delete(m, y), Nothing)
     @test !haskey(m.var_to_fix, JuMP.index(y))
     @test !haskey(m.var_to_integrality, JuMP.index(y))
     @test name(meas1) == "measure(par)"
     @test name(meas2) == "measure(0)"
     @test !haskey(m.var_to_meas, JuMP.index(y))
     @test string(m.constrs[JuMP.index(con1)].func) == "par"
     @test string(m.constrs[JuMP.index(con2)].func) == "0"
     @test !haskey(m.var_to_constrs, JuMP.index(y))
     @test string(objective_function(m)) == "0"
     @test !haskey(m.infinite_to_points, JuMP.index(inf))
     @test !haskey(m.var_in_objective, JuMP.index(y))
     @test !haskey(m.vars, JuMP.index(y))
     @test !haskey(m.var_to_name, JuMP.index(y))
     # test errors
     @test_throws AssertionError delete(m, x)
     @test_throws AssertionError delete(m, y)
 end

 # Test variable deletion
 @testset "JuMP.delete (Infinite Variables)" begin
     # intialize the model
     m = InfiniteModel()
     @infinite_parameter(m, 0 <= par <= 1)
     @infinite_variable(m, 0 <= x(par) <= 1, Bin)
     @infinite_variable(m, y(par) == 1, Int)
     @point_variable(m, x(0), x0)
     m.reduced_info[-1] = ReducedInfiniteInfo(x, Dict(1 => 0.5))
     m.infinite_to_reduced[JuMP.index(x)] = [-1]
     rv = ReducedInfiniteVariableRef(m, -1)
     data = DiscreteMeasureData(par, [1], [1])
     meas1 = measure(x + y + par, data)
     meas2 = add_measure(m, Measure(y, data))
     @constraint(m, con1, x + y + par <= 0)
     con2 = add_constraint(m, ScalarConstraint(y, MOI.LessThan(0.)))
     # test deletion of x
     @test isa(delete(m, x), Nothing)
     @test !haskey(m.var_to_lower_bound, JuMP.index(x))
     @test !haskey(m.var_to_upper_bound, JuMP.index(x))
     @test !haskey(m.var_to_zero_one, JuMP.index(x))
     @test name(meas1) == "measure(y(par) + par)"
     @test !haskey(m.var_to_meas, JuMP.index(x))
     @test string(m.constrs[JuMP.index(con1)].func) == "y(par) + par"
     @test !haskey(m.var_to_constrs, JuMP.index(x))
     @test m.param_to_vars[JuMP.index(par)] == [JuMP.index(y)]
     @test !is_valid(m, x0)
     @test !haskey(m.infinite_to_points, JuMP.index(x))
     @test !is_valid(m, rv)
     @test !haskey(m.infinite_to_reduced, JuMP.index(x))
     @test !haskey(m.var_in_objective, JuMP.index(x))
     @test !haskey(m.vars, JuMP.index(x))
     @test !haskey(m.var_to_name, JuMP.index(x))
     # test deletion of y
     @test isa(delete(m, y), Nothing)
     @test !haskey(m.var_to_fix, JuMP.index(y))
     @test !haskey(m.var_to_integrality, JuMP.index(y))
     @test name(meas1) == "measure(par)"
     @test name(meas2) == "measure(0)"
     @test !haskey(m.var_to_meas, JuMP.index(y))
     @test string(m.constrs[JuMP.index(con1)].func) == "par"
     @test string(m.constrs[JuMP.index(con2)].func) == "0"
     @test !haskey(m.var_to_constrs, JuMP.index(y))
     @test !haskey(m.param_to_vars, JuMP.index(par))
     @test !haskey(m.var_in_objective, JuMP.index(y))
     @test !haskey(m.vars, JuMP.index(y))
     @test !haskey(m.var_to_name, JuMP.index(y))
     # test errors
     @test_throws AssertionError delete(m, x)
     @test_throws AssertionError delete(m, y)
 end

 # Test variable deletion
 @testset "JuMP.delete (Measures)" begin
     # intialize the model
     m = InfiniteModel()
     @infinite_parameter(m, 0 <= par <= 1)
     @infinite_parameter(m, 0 <= par2 <= 1)
     @infinite_parameter(m, 0 <= pars[1:2] <= 1)
     @infinite_variable(m, w(pars))
     @infinite_variable(m, x(par))
     @infinite_variable(m, y(par2))
     @point_variable(m, x(0), x0)
     m.reduced_info[-1] = ReducedInfiniteInfo(x, Dict(1 => 0.5))
     m.infinite_to_reduced[JuMP.index(x)] = [-1]
     rv = ReducedInfiniteVariableRef(m, -1)
     data = DiscreteMeasureData(par, [1], [1])
     data2 = DiscreteMeasureData(pars, [1], [[1, 1]])
     meas = measure(x, data)
     meas1 = measure(x + x0 + rv + par + meas + y + par2, data)
     meas2 = measure(x + x0 + rv + par + y + par2, data)
     meas3 = measure(meas1 + x0, data)
     meas4 = measure(meas2, data)
     meas5 = measure(w, data2)
     @constraint(m, con1, x0 + meas1 <= 0)
     con2 = add_constraint(m, ScalarConstraint(meas2, MOI.LessThan(0.)))
     @objective(m, Min, meas1 + x0)
     # test deletion of meas1
     @test isa(delete(m, meas1), Nothing)
     @test name(meas3) == "measure(x0)"
     @test !haskey(m.meas_to_meas, JuMP.index(meas1))
     @test string(m.constrs[JuMP.index(con1)].func) == "x0"
     @test !haskey(m.meas_to_constrs, JuMP.index(meas1))
     @test string(objective_function(m)) == "x0"
     @test m.var_to_meas[JuMP.index(x)] == [JuMP.index(meas), JuMP.index(meas2)]
     @test m.var_to_meas[JuMP.index(y)] == [JuMP.index(meas2)]
     @test m.param_to_meas[JuMP.index(par)] == [1, 3, 4, 5]
     @test m.param_to_meas[JuMP.index(par2)] == [JuMP.index(meas2)]
     @test !haskey(m.meas_to_meas, JuMP.index(meas))
     @test m.reduced_to_meas[JuMP.index(rv)] == [JuMP.index(meas2)]
     @test !haskey(m.meas_in_objective, JuMP.index(meas1))
     @test !haskey(m.measures, JuMP.index(meas1))
     @test !haskey(m.meas_to_name, JuMP.index(meas1))
     # test deletion of meas2
     set_objective_function(m, meas2)
     @test isa(delete(m, meas2), Nothing)
     @test name(meas4) == "measure(0)"
     @test !haskey(m.meas_to_meas, JuMP.index(meas1))
     @test string(m.constrs[JuMP.index(con2)].func) == "0"
     @test !haskey(m.meas_to_constrs, JuMP.index(meas2))
     @test string(objective_function(m)) == "0"
     @test m.var_to_meas[JuMP.index(x)] == [JuMP.index(meas)]
     @test !haskey(m.var_to_meas, JuMP.index(y))
     @test m.param_to_meas[JuMP.index(par)] == [1, 4, 5]
     @test !haskey(m.param_to_meas, JuMP.index(par2))
     @test !haskey(m.reduced_to_meas, JuMP.index(rv))
     @test !haskey(m.meas_in_objective, JuMP.index(meas2))
     @test !haskey(m.measures, JuMP.index(meas2))
     @test !haskey(m.meas_to_name, JuMP.index(meas2))
     # test deletion of meas5
     @test isa(delete(m, meas5), Nothing)
     @test !haskey(m.param_to_meas, JuMP.index(pars[1]))
     @test !haskey(m.param_to_meas, JuMP.index(pars[2]))
     # test errors
     @test_throws AssertionError delete(m, meas1)
     @test_throws AssertionError delete(m, meas2)
 end


text:
module neural_network

using NPZ

σ(x::AbstractFloat) = 1 / (1 + exp(-x))
σ_prime(y::AbstractFloat) = 0.25 + y * (1 - y)

struct NeuralNetwork
    neuron_counts::Vector{Int}
    weights::Vector{Matrix{Float64}}
    biases::Vector{Vector{Float64}}
    layers::Vector{Vector{Float64}}
end

@doc """
Construct a neural network with the provided neuron counts, weights and biases.
"""
NeuralNetwork(neuron_counts::Vector{Int}, weights::Vector{Matrix{Float64}}, biases::Vector{Vector{Float64}}) = NeuralNetwork(
    neuron_counts,
    weights,
    biases,
    [Vector{Float64}(undef, neuron_count) for neuron_count ∈ neuron_counts]
)

@doc """
Construct a neural network with random weights and biases.
"""
NeuralNetwork(neuron_counts::Vector{Int}) = NeuralNetwork(
    neuron_counts,
    [
        Array{Float64}(undef, 0, 0),
        [
            (rand(Float64, (next_layer_neuron_count, prev_layer_neuron_count)) .* 10) .- 5
            for (next_layer_neuron_count, prev_layer_neuron_count) ∈ zip(neuron_counts[2:end], neuron_counts)
        ]...
    ],
    [
        Array{Float64}(undef, 0),
        [rand(neuron_count) .* 20 .- 10 for neuron_count ∈ neuron_counts[2:end]]...
    ]
)

@doc """
Construct a neural network from params stored in the provided directory path.
"""
NeuralNetwork(neuron_counts::Vector{Int}, directory::String) = NeuralNetwork(
    neuron_counts,
    [npzread("$(directory)/weight$(i).npy") for i in 1:length(neuron_counts)],
    [npzread("$(directory)/bias$(i).npy") for i in 1:length(neuron_counts)]
)

@doc """
Run the neural network on the provided input image.
"""
function feed_forward!(network::NeuralNetwork, image::AbstractVector{Float64})
    network.layers[1] .= image
    @inbounds for i ∈ 2:length(network.neuron_counts)
        network.layers[i] .= σ.((network.weights[i] * network.layers[i - 1]) .+ network.biases[i])
    end
end

@doc """
Given a network run on some input, the target output for that input and the learning rate, backpropagate and rectify weights and biases.
"""
function back_propagate!(network::NeuralNetwork, target::Vector{Float64}, η::Float64)
    delc_dela::Vector{Float64} = 2 .* (network.layers[end] .- target)
    @inbounds for L ∈ length(network.neuron_counts):-1:2
        delc_delz::Vector{Float64} = delc_dela .* σ_prime.(network.layers[L])

        delc_delw::Matrix{Float64} = [j * k for j ∈ delc_delz, k ∈ network.layers[L-1]]
        delc_delb::Vector{Float64} = delc_delz
        delc_dela = vec(sum((.*).(network.weights[L], delc_delz), dims=1))

        # update weights and biases at this layer
        network.biases[L] .+= -η * delc_delb
        network.weights[L] .+= -η * delc_delw
    end
end

@doc """
Save the weights and biases of the network on the disk.
"""
function save_params(network::NeuralNetwork)
    mkpath("src/params")
    @info "saving weights"
    for (i, weight) in enumerate(network.weights)
        npzwrite("src/params/weight$(i).npy", weight)
    end
    @info "saved weights"

    @info "saving biases"
    for (i, bias) in enumerate(network.biases)
        npzwrite("src/params/bias$(i).npy", bias)
    end
    @info "saved biases"
end

@doc """
Train the neural network with a single image and its label.
"""
function train!(network::NeuralNetwork, image::AbstractVector{Float64}, label::Int, η::Float64 = 0.05)
    feed_forward!(network, image)

    prediction = findmax(network.layers[end])
    @debug "" predict = (prediction[2] - 1) actual=label confidence_actual=network.layers[end][label + 1] confidence_predict=prediction[1]

    target = [i == label ? 1. : 0. for i ∈ 0:9]
    back_propagate!(network, target, η)
end

export NeuralNetwork, feed_forward!, back_propagate!, train!, save_params
end

max_stars_repo_name: ShakyaMajumdar/neural-network

max_issues_repo_name: ShakyaMajumdar/neural-network

max_forks_repo_name: ShakyaMajumdar/neural-network