text:
include("../src/GeneralizedStateSpace.jl")

using DifferentialEquations
using Plots: plot, plot!, theme
using Interact: @manipulate, settheme!

theme(:juno)
settheme!(:nativehtml)

# ż = -2z + u
# ẍ = -x + z
# y = 3x + z
#
# ż = -2z + u
# ẋ = v
# v̇ = -x + z
# y = z + 3x
#
# dx[1] = -2*x[1]               + u
# dx[2] =                  x[3]
# dx[3] =    x[1] - x[2]

function f1!(dx, x, u, t)
    dx[1] = -2*x[1] + u
    return nothing
end
function f2!(dx, x, u, t)
    dx[1] = x[2]
    dx[2] = u[1] - x[1]
    return nothing
end
function f!(dx, x, u, t)
    dx[1] = -2*x[1] + u
    dx[2] = x[3]
    dx[3] = x[1] - x[2]
    return nothing
end

h1(x, t) = x[1]
h2(x, u, t) = u[1] + 3*x[1]
h(x, t) = x[1] + 3*x[2]

sys_base = NLStateSpace(StateFunction(f!, 1), OutputFunction(h, 1))
sys1 = NLStateSpace(StateFunction(f1!, 1), OutputFunction(h1, 1, 1))
sys2 = NLStateSpace(StateFunction(f2!, 1), OutputFunction(h2, 1, 1))
sys_comp = sys2 * sys1

# Pulse width modulation
pwm(t, period, duty) = mod(t, period) ≤ duty*period ? 1 : 0
pwm(period, duty) = t -> pwm(t, period, duty)

# Trigger times for pwm to be passed into the solver tstops keyword argument
vzip(vect1, vect2) = vcat([vcat(v1, v2) for (v1, v2) in zip(vect1, vect2)]...)
function pwm_triggers(stoptime, period, duty)
    rise = period:period:stoptime
    fall = rise .+ duty*period
    return vzip(rise, fall)
end

x₀ = [0., 0., 0.]
tspan = (0., 20.)

# Pulse width modulation
# period = 0.1
# duty = 1/3
# u(t) = pwm(t, period, duty)
u = 1

@manipulate for duty = 0.01:0.01:1-0.01, period=0.01:0.01:3
    u_const = duty
    tstops = pwm_triggers(tspan[2], period, duty)

    affect!(integrator) = (integrator.p = 1 - integrator.p)
    time_choice(integrator) = integrator.t + period*(integrator.p==1 ? duty : (1-duty))
    cb = IterativeCallback(time_choice, affect!)

    prob_base = ODEProblem(ODEFunction{true}(sys_base.f), x₀, tspan, u_const)
    prob_comp = ODEProblem(ODEFunction{true}(sys_base.f), x₀, tspan, u)

    base = solve(prob_base, Tsit5())
    comp = solve(prob_comp, Tsit5(), callback=cb)

    time = comp.t
    detailed_time = LinRange(tspan..., 2000)

    plot(time, repeat([u_const], length(time)), ylim=[0, 3.5], label="duty cycle")
    plot!(detailed_time, pwm.(detailed_time, period, duty), label="input pwm")
    plot!(time, sys_comp.h.(comp.(time), u, time), label="composed")
    plot!(time, sys_base.h.(base.(time), u, time), label="baseline")
end

max_stars_repo_name: jonniedie/ComposableDiffEqs

max_issues_repo_name: jonniedie/ComposableDiffEqs

max_forks_repo_name: jonniedie/ComposableDiffEqs


text:
# ------------------------------------------------------------------
# Licensed under the MIT License. See LICENSE in the project root.
# ------------------------------------------------------------------

"""
    AbstractData

An abstract type to aid with custom spatial data types. Users can
subtype their domains from this type, and implement the methods in
`geotraits/data.jl`.
"""
abstract type AbstractData end

"""
    sdata₁ == sdata₂

Tells whether or not the spatial `sdata₁` and `sdata₂` are equal.
"""
==(sdata₁::AbstractData, sdata₂::AbstractData) =
  domain(sdata₁) == domain(sdata₂) &&
  values(sdata₁) == values(sdata₂)

# -----------
# TABLES API
# -----------

Tables.istable(::Type{<:AbstractData}) = true
Tables.materializer(sdata::AbstractData) = Tables.materializer(values(sdata))
Tables.columnaccess(sdata::AbstractData) = Tables.columnaccess(values(sdata))
Tables.rowaccess(sdata::AbstractData) = Tables.rowaccess(values(sdata))
Tables.schema(sdata::AbstractData) = Tables.schema(values(sdata))
Tables.columns(sdata::AbstractData) = Tables.columns(values(sdata))
Tables.columnnames(sdata::AbstractData) = Tables.columnnames(values(sdata))
Tables.getcolumn(sdata::AbstractData, c::Symbol) = Tables.getcolumn(values(sdata), c)
Tables.rows(sdata::AbstractData) = Tables.rows(values(sdata))

# -------------
# VARIABLE API
# -------------

function variables(sdata::AbstractData)
  s = Tables.schema(sdata)
  ns, ts = s.names, s.types
  @. Variable(ns, nonmissing(ts))
end

Base.getindex(sdata::AbstractData, var::Symbol) =
  Tables.getcolumn(sdata, var)

# ---------
# VIEW API
# ---------

Base.view(sdata::AbstractData, inds::AbstractVector{Int}) =
  DataView(sdata, inds, collect(name.(variables(sdata))))
Base.view(sdata::AbstractData, vars::AbstractVector{Symbol}) =
  DataView(sdata, 1:nelms(sdata), vars)
Base.view(sdata::AbstractData, inds, vars) =
  DataView(sdata, inds, vars)

#------------------
# IMPLEMENTATIONS
#------------------
"""
    SpatialData(domain, data)

Tabular `data` georeferenced in a spatial `domain`.
"""
struct SpatialData{𝒟,𝒯} <: AbstractData
  domain::𝒟
  table::𝒯

  function SpatialData{𝒟,𝒯}(domain, table) where {𝒟,𝒯}
    ne = nelms(domain)
    nr = length(Tables.rows(table))
    @assert ne == nr "number of table rows ≠ number of mesh elements"
    new(domain, table)
  end
end

SpatialData(domain::𝒟, table::𝒯) where {𝒟,𝒯} =
  SpatialData{𝒟,𝒯}(domain, table)

domain(sdata::SpatialData) = sdata.domain
values(sdata::SpatialData) = sdata.table

# ------------
# IO methods
# ------------
function Base.show(io::IO, sdata::SpatialData)
  N = ncoords(sdata)
  T = coordtype(sdata)
  n = nelms(sdata)
  print(io, "$n SpatialData{$T,$N}")
end

function Base.show(io::IO, ::MIME"text/plain", sdata::SpatialData)
  𝒟 = domain(sdata)
  𝒯 = values(sdata)
  s = Tables.schema(𝒯)
  vars = zip(s.names, s.types)
  println(io, 𝒟)
  println(io, "  variables")
  varlines = ["    └─$var ($V)" for (var,V) in vars]
  print(io, join(sort(varlines), "\n"))
end

ext: jl


text:
# https://www.philipzucker.com/egraph-2/
# https://github.com/philzook58/EGraphs.jl/blob/main/src/matcher.jl
# https://www.hpl.hp.com/techreports/2003/HPL-2003-148.pdf
# TODO support destructuring

# ematching seems to be faster without spawning tasks

# we keep a pair of EClass, Any in substitutions because
# when evaluating dynamic rules we also want to know
# what was the value of a matched literal
const Sub = Base.ImmutableDict{Any, Tuple{EClass, Any}}
const SubBuf = Vector{Sub}

# https://www.hpl.hp.com/techreports/2003/HPL-2003-148.pdf
# page 48
"""
From [https://www.hpl.hp.com/techreports/2003/HPL-2003-148.pdf](https://www.hpl.hp.com/techreports/2003/HPL-2003-148.pdf)
The iterator `ematchlist` matches a list of terms `t` to a list of E-nodes by first finding
all substitutions that match the first term to the first E-node, and then extending
each such substitution in all possible ways that match the remaining terms to
the remaining E-nodes. The base case of this recursion is the empty list, which
requires no extension to the substitution; the other case relies on Match to find the
substitutions that match the first term to the first E-node.
"""
function ematchlist(e::EGraph, t::AbstractVector{Pattern}, v::AbstractVector{Int64}, sub::Sub, buf::SubBuf)::SubBuf
    if length(t) == 0 || length(v) == 0
        if !isempty(sub)
            push!(buf, sub)
        end
    else
        for sub1 in ematchstep(e, t[1], v[1], sub, SubBuf(), nothing)
            ematchlist(e, (@view t[2:end]), (@view v[2:end]), sub1, buf)
        end
    end
    return buf
end

# Tries to match on a pattern variable
function ematchstep(g::EGraph, t::PatVar, v::Int64, sub::Sub, buf::SubBuf, lit=nothing)::SubBuf
    if haskey(sub, t)
        if find(g, first(sub[t])) == find(g, v)
            push!(buf, sub)
        end
    else
        push!(buf, Base.ImmutableDict(sub, t => (geteclass(g, find(g, v)), lit)))
    end
    return buf
end

# Tries to match on literals
function ematchstep(g::EGraph, t::PatLiteral, v::Int64, sub::Sub, buf::SubBuf, lit=nothing)::SubBuf
    ec = geteclass(g, v)
    for n in ec
        if arity(n) == 0 && t.val == n.head
            if haskey(sub, t)
                if find(g, first(sub[t])) == ec.id
                    push!(buf, sub)
                end
            else
                push!(buf, Base.ImmutableDict(sub, t => (ec, n.head)))
            end
        end
    end
    return buf
end

# tries to match on type assertions
function ematchstep(g::EGraph, t::PatTypeAssertion, v::Int64, sub::Sub, buf::SubBuf, lit=nothing)::SubBuf
    ec = geteclass(g, v)
    for n in ec
        if arity(n) == 0
            !(typeof(n.head) <: t.type) && continue
            ematchstep(g, t.var, v, sub, buf, n.head)
            continue
        end
    end
    return buf
end

# PATEQUIV mechanism
function ematchstep(g::EGraph, t::PatEquiv, v::Int64, sub::Sub, buf::SubBuf, lit=nothing)::SubBuf
    buf1 = SubBuf()
    buf2 = SubBuf()

    for sub1 ∈ ematchstep(g, t.left, v, sub, buf1)
        ematchstep(g, t.right, v, sub1, buf2)
    end

    if !isempty(buf1) && !isempty(buf2)
        for sub ∈ vcat(buf1, buf2)
            push!(buf, sub)
        end
    end
    return buf
end

# tries to match on composite expressions
function ematchstep(g::EGraph, t::PatTerm, v::Int64, sub::Sub, buf::SubBuf, lit=nothing)::SubBuf
    ec = geteclass(g, v)
    for n in ec
        (arity(n) > 0) && n.head == t.head && arity(t) == arity(n) || continue
         ematchlist(g, t.args, n.args, sub, buf)
    end
    return buf
end

function ematchstep(g::EGraph, t::PatAllTerm, v::Int64, sub::Sub, buf::SubBuf, lit=nothing)::SubBuf
    ec = geteclass(g, v)
    for n in ec
        (arity(n) > 0) && arity(t) == arity(n) || continue
        if haskey(sub, t.head)
            @show sub[t.head]
            if find(g, first(sub[t.head])) == find(g, v)
                ematchlist(g, t.args, n.args, sub, buf)
            end
        else
            nsub = Base.ImmutableDict(sub, t.head => (geteclass(g, find(g, v)), lit))
            ematchlist(g, t.args, n.args, nsub, buf)
        end
    end
    return buf
end

function ematch(g::EGraph, pat::Pattern, id::Int64, sub=Sub(), buf=SubBuf())
    ematchstep(g, pat, id, sub, buf)
    # @show pat
    # println.(buf)
    return buf
end


text:
#= -------------------------------------------------------------------------
# @file problem1.jl
#
# @date 12/22/15 21:56:50
# @author Martin Noblia
# @email martin.noblia@openmailbox.org
#
# @brief
#
# @detail
#
  Licence:
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or (at
# your option) any later version.
#
# This program is distributed in the hope that it will be useful, but
# WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
# General Public License for more details.
#
# You should have received a copy of the GNU General Public License

---------------------------------------------------------------------------=#
# Enun
#=
If we list all the natural numbers below 10 that are multiples of 3 or 5,
we get 3, 5, 6 and 9. The sum of these multiples is 23.
Find the sum of all the multiples of 3 or 5 below 1000.
=#

"""
Function for collect the integers multiples of some integer
Input:
N(Int) Maximum integer
n(Int) integer multiple
Output:
collection of integers multiples of n
"""
mult_of(N::Int,n::Int) = find(collect(1:N) .% n .== 0)
s = sum(unique(vcat(mult_of(999, 3), mult_of(999, 5))));
print("The sum of all the multiples of 3 or 5 bellow 1000 is: $s")

max_stars_repo_name: elsuizo/Euler

max_issues_repo_name: elsuizo/Euler

max_forks_repo_name: elsuizo/Euler


text:
export savelpconstraints, MatrixConstraintSet

type MatrixConstraintSet
    components::Vector{Symbol}
    names::Vector{Symbol}
    lowers::Vector{Float64}
    uppers::Vector{Float64}
    f::Vector{Float64}
    A::SparseMatrixCSC{Float64, Int64}
    b::Vector{Float64}
end

function lpconstraints(model::Model, components::Vector{Symbol}, names::Vector{Symbol}, objective::Function, constraints::Vector{Function}; verbose=false)
    if verbose
        println("Making zero point.")
    end
    initial = make0(model, names)

    if model.numberType == Number
        myobjective = unaryobjective(model, components, names, objective)
        f = ForwardDiff.gradient(myobjective, initial)

        myconstraints = Function[]
        for constraint in constraints
            myconstraints = [myconstraints; unaryobjective(model, components, names, constraint)]
        end

        # Copied from makematrix.jl
        # Could be made more efficient to get b values from gradient
        A = Float64[]
        b = Float64[]
        for myconstraint in myconstraints
            A = [A; ForwardDiff.gradient(myconstraint, initial)]
            b = [b; myconstraint(initial)] # Not the final values for 'b' yet!
        end

        A = reshape(A, (length(initial), div(length(A), length(initial))))'
    else
        if verbose
            println("Adding to gradient collection.")
        end
        rg = RegisterGradient(model, components, names, 1.0)
        addfunction!(rg, objective)
        for constraint in constraints
            addfunction!(rg, constraint)
        end

        vb, fA = getgradients(rg, initial, verbose=verbose, usesparse=false)
        f = vec(fA[1, :])
        A = fA[2:end, :]
        b = vb[2:end]
    end

    b = A * initial - b

    f, b, A
end

function savelpconstraints(model::Model, components::Vector{Symbol}, names::Vector{Symbol}, lowers::Vector{Float64}, uppers::Vector{Float64}, objective::Function, constraints::Vector{Function}=Function[]; verbose=false)
    f, b, A = lpconstraints(model, components, names, objective, constraints, verbose=verbose)

    sA = sparse(A)
    MatrixConstraintSet(components, names, lowers, uppers, f, sA, b)
end

function combineconstraints(f::Vector{Float64}, b::Vector{Float64}, A::Matrix{Float64}, model::Model, components::Vector{Symbol}, names::Vector{Symbol}, constraints::Vector{MatrixConstraintSet})
    for constraint in constraints
        f = [f; constraint.f] # these need to be variables never seen before!
        b = [b; constraint.b]
        A, components, names = addnamedcols(A, model, components, names, constraint.components, constraint.names)
        A = addnamedrows(A, model, components, names, constraint.components, constraint.names, full(constraint.A))
    end

    f, b, A, components, names
end

function addnamedcols(A::Matrix{Float64}, model::Model, components::Vector{Symbol}, names::Vector{Symbol}, pluscomponents::Vector{Symbol}, plusnames::Vector{Symbol})
    found = Bool[]
    for ii in 1:length(pluscomponents)
        found = [found; any((components == pluscomponents[ii]) & (names == plusnames[ii]))]
    end

    # Add 0's for any unfound
    for (ii, len, isscalar) in @task nameindexes(model, plusnames[!found])
        A = [A zeros(size(A)[1], len)]
    end

    A, [components; pluscomponents[!found]], [names; plusnames[!found]]
end

function addnamedrows(A::Matrix{Float64}, model::Model, components::Vector{Symbol}, names::Vector{Symbol}, pluscomponents::Vector{Symbol}, plusnames::Vector{Symbol}, plusA::Matrix{Float64})
    rows = zeros(size(plusA)[1], size(A)[2])

    # Get the indexes into plusA
    plusAindexes = Int64[]
    startindex = 1
    for (ii, len, isscalar) in @task nameindexes(model, plusnames)
        plusAindexes = [plusAindexes; startindex]
        startindex += len
    end
    plusAindexes = [plusAindexes; startindex]

    # Filter these into the right spots in rows
    startindex = 1
    for (ii, len, isscalar) in @task nameindexes(model, names)
        if names[ii] in plusnames
            which = find(plusnames .== names[ii])[1]
            println(startindex:startindex+len-1)
            println(plusAindexes[which])
            println(plusAindexes[which+1])
            rows[:, startindex:startindex+len-1] = plusA[:, plusAindexes[which]:plusAindexes[which+1]-1]
        end
    end

    [A; rows]
end

function combinelimits(exlowers::Vector{Float64}, exuppers::Vector{Float64}, model::Model, components::Vector{Symbol}, names::Vector{Symbol}, constraints::Vector{MatrixConstraintSet})
    for constraint in constraints
        for (ii, len, isscalar) in @task nameindexes(model, constraint.names)
            append!(exlowers, [constraint.lowers[ii] for jj in 1:len])
            append!(exuppers, [constraint.uppers[ii] for jj in 1:len])
        end
    end

    exlowers, exuppers
end


text:
function normalized_substrate_flux!(F_c::Vector{T}, S::Vector{T}, K::Matrix{T}) where {T<:Real}
    (I,J) = size(K)
    @assert I == size(K,1)
    @assert J == size(K,2)
    @inbounds for l in 1:I, j in 1:J
            val = S[l]/K[l,j]
            F_c[j] = F_c[j] + val[1]
    end
    return F_c
end # normalized_substrate_flux!

function normalized_substrate_flux!(F_c::Vector{T}, S::Vector{T}, K::Matrix{T}, α_kin::Matrix{T}) where {T<:Real}
    (I,J) = size(K)
    @assert I == size(K,1)
    @assert J == size(K,2)
    @inbounds for j in 1:J, l in 1:I
            val = α_kin[l,j]*S[l]/K[l,j]
            F_c[j] = F_c[j] + val[1]
    end
    return F_c
end # normalized_substrate_flux!

function conjugate_substrate_flux!(F_r::Vector{T}, E::Vector{T}, K::Matrix{T}) where {T<:Real}
    (I,J) = size(K)
    @assert J == size(E,1)
    @assert I == size(F_r,1)
    @inbounds for i in 1:I, l in 1:J
        val = E[l]/K[i,l]
        F_r[i] = F_r[i] + val[1]
    end
    return F_r
end

function conjugate_substrate_flux!(F_r::Vector{T}, E::Vector{T}, K::Matrix{T}, N_SB::Matrix{T}) where {T<:Real}
    (I,J) = size(K)
    @assert J == size(E,1)
    @assert I == size(F_r,1)
    @inbounds for i in 1:I, l in 1:J
        val = N_SB[i,l]*E[l]/K[i,l]
        F_r[i] = F_r[i] + val[1]
    end
    return F_r
end

function ECA_kinetics!(ECA::Matrix{T}, S::Vector{T}, E::Vector{T}, K::Matrix{T}, k2p::Matrix{T}) where {T<:Real}
    (I,J) = size(K)
    @assert I == size(S,1)
    @assert J == size(E,1)
    @assert (I,J) == size(k2p)
    F_c = zeros(J)
    F_r = zeros(I)
    F_r = conjugate_substrate_flux!(F_r, E, K)
    F_c = normalized_substrate_flux!(F_c, S, K)
    @inbounds for i in 1:I, j in 1:J
            val = k2p[i,j]*E[j]*(S[i]/K[i,j])/(1 + F_r[i] + F_c[j])
            ECA[i,j] = val[1]
    end

    any(x->x==true, isnan.(ECA)) ? throw(DomainError("NaN in DEBJulia.ECA_kinetics!")) : return ECA
end

function ECA_kinetics!(ECA::Matrix{T}, S::Vector{T}, E::Vector{T}, K::Matrix{T}, k2p::Matrix{T}, N_SB::Matrix{T}) where {T<:Real}
    (I,J) = size(ECA)
    @assert I == size(S,1)
    @assert J == size(E,1)
    @assert (I,J) == size(k2p)
    F_c = zeros(J)
    F_r = zeros(I)
    F_r = conjugate_substrate_flux!(F_r, E, K, N_SB)
    F_c = normalized_substrate_flux!(F_c, S, K)
    @inbounds for i in 1:I, j in 1:J
            val = k2p[i,j]*N_SB[i,j]*E[j]*(S[i]/K[i,j])/(1 + F_r[i] + F_c[j])
            ECA[i,j] = val[1]
    end

    #any(x->x==true, isnan.(ECA)) ? throw(DomainError("NaN in DEBJulia.ECA_kinetics!")) : return ECA
    any(x->x==true, isnan.(ECA)) ? throw(DomainError("NaN in DEBJulia.ECA_kinetics!"))  : return ECA
end

function ECA_kinetics!(ECA::Matrix{T}, S::Vector{T}, E::Vector{T}, K::Matrix{T}, k2p::Matrix{T}, α_kin::Matrix{T}, t::String) where {T<:Real}
    (I,J) = size(ECA)
    @assert I == size(S,1)
    @assert J == size(E,1)
    @assert (I,J) == size(k2p)
    F_c = zeros(J)
    F_r = zeros(I)
    F_r = conjugate_substrate_flux!(F_r, E, K)
    F_c = normalized_substrate_flux!(F_c, S, K, α_kin)
    @inbounds for i in 1:I, j in 1:J
            val = k2p[i,j]*E[j]*(α_kin[i,j]*S[i]/K[i,j])/(1 + F_r[i] + F_c[j])
            ECA[i,j] = val[1]
    end

    any(x->x==true, isnan.(ECA)) ? throw(DomainError("NaN in DEBJulia.ECA_kinetics!")) : return ECA
end


text:
module CAPI

const libmine = joinpath(dirname(dirname(@__FILE__)), joinpath("c",
                        @static if is_linux() "libmine.so"
                                elseif is_osx() "limine.dylab"
                                else error("this package only support osx and linux")
                                end))
function __init__()
    println(libmine)
end

macro fncall(fn, return_type)
    f = eval(fn)
    quote
        ccall(($(Meta.quot(f)), $libmine), $return_type,
              ())
    end
end

macro fncall(fn, return_type, argtypes, args...)
    f = eval(fn)
    args = map(esc, args)
    quote
        ccall(($(Meta.quot(f)), $libmine), $return_type,
              $argtypes, $(args...))
    end
end

include("type.jl")
include("function.jl")

end


text:
"""
    ReluVal(max_iter::Int64, tree_search::Symbol)

ReluVal combines symbolic reachability analysis with iterative interval refinement to
minimize over-approximation of the reachable set.

# Problem requirement
1. Network: any depth, ReLU activation
2. Input: Hyperrectangle
3. Output: LazySet

# Return
`CounterExampleResult`

# Method
Symbolic reachability analysis and iterative interval refinement (search).
- `max_iter` default `10`.
- `tree_search` default `:DFS` - depth first search.

# Property
Sound but not complete.

# Reference
[S. Wang, K. Pei, J. Whitehouse, J. Yang, and S. Jana,
"Formal Security Analysis of Neural Networks Using Symbolic Intervals,"
*CoRR*, vol. abs/1804.10829, 2018. arXiv: 1804.10829.](https://arxiv.org/abs/1804.10829)

[https://github.com/tcwangshiqi-columbia/ReluVal](https://github.com/tcwangshiqi-columbia/ReluVal)
"""
@with_kw struct ReluVal <: Solver
    max_iter::Int64     = 10
    tree_search::Symbol = :DFS # only :DFS/:BFS allowed? If so, we should assert this.
end

function solve(solver::ReluVal, problem::Problem)
    isbounded(problem.input) || throw(UnboundedInputError("ReluVal can only handle bounded input sets."))

    reach_list = []
    interval = problem.input
    for i in 1:solver.max_iter
        if i > 1
            interval = select!(reach_list, solver.tree_search)
        end
        reach = forward_network(solver, problem.network, init_symbolic_mask(interval))
        result = check_inclusion(reach.sym, problem.output, problem.network)

        if result.status === :violated
            return result
        elseif result.status === :unknown
            intervals = bisect_interval_by_max_smear(problem.network, reach)
            append!(reach_list, intervals)
        end
        isempty(reach_list) && return CounterExampleResult(:holds)
    end
    return CounterExampleResult(:unknown)
end

function bisect_interval_by_max_smear(nnet::Network, reach::SymbolicIntervalMask)
    LG, UG = get_gradient_bounds(nnet, reach.LΛ, reach.UΛ)
    feature, monotone = get_max_smear_index(nnet, reach.sym.domain, LG, UG) #monotonicity not used in this implementation.
    return collect(split_interval(reach.sym.domain, feature))
end

function select!(reach_list, tree_search)
    if tree_search == :BFS
        reach = popfirst!(reach_list)
    elseif tree_search == :DFS
        reach = pop!(reach_list)
    else
        throw(ArgumentError(":$tree_search is not a valid tree search strategy"))
    end
    return reach
end

function check_inclusion(reach::SymbolicInterval{<:Hyperrectangle}, output, nnet::Network)
    reachable = Hyperrectangle(low = low(reach), high = high(reach))

    issubset(reachable, output) && return CounterExampleResult(:holds)

    # Sample the middle point
    middle_point = center(domain(reach))
    y = compute_output(nnet, middle_point)
    y ∈ output || return CounterExampleResult(:violated, middle_point)

    return CounterExampleResult(:unknown)
end

# Symbolic forward_linear
function forward_linear(solver::ReluVal, L::Layer, input::SymbolicIntervalMask)
    output_Low, output_Up = interval_map(L.weights, input.sym.Low, input.sym.Up)
    output_Up[:, end] += L.bias
    output_Low[:, end] += L.bias
    sym = SymbolicInterval(output_Low, output_Up, domain(input))
    return SymbolicIntervalGradient(sym, input.LΛ, input.UΛ)
end

# Symbolic forward_act
function forward_act(::ReluVal, L::Layer{ReLU},  input::SymbolicIntervalMask)
    output_Low, output_Up = copy(input.sym.Low), copy(input.sym.Up)
    n_node = n_nodes(L)
    LΛᵢ, UΛᵢ = falses(n_node), trues(n_node)

    for j in 1:n_node
        # If the upper bound of the upper bound is negative, set
        # the generators and centers of both bounds to 0, and
        # the gradient mask to 0
        if upper_bound(upper(input), j) <= 0
            LΛᵢ[j], UΛᵢ[j] = 0, 0
            output_Low[j, :] .= 0
            output_Up[j, :] .= 0

        # If the lower bound of the lower bound is positive,
        # the gradient mask should be 1
        elseif lower_bound(lower(input), j) >= 0
            LΛᵢ[j], UΛᵢ[j] = 1, 1

        # if the bounds overlap 0, concretize by setting
        # the generators to 0, and setting the new upper bound
        # center to be the current upper-upper bound.
        else
            LΛᵢ[j], UΛᵢ[j] = 0, 1
            output_Low[j, :] .= 0
            if lower_bound(upper(input), j) < 0
                output_Up[j, :] .= 0
                output_Up[j, end] = upper_bound(upper(input), j)
            end
        end
    end

    sym = SymbolicInterval(output_Low, output_Up, domain(input))
    LΛ = push!(input.LΛ, LΛᵢ)
    UΛ = push!(input.UΛ, UΛᵢ)
    return SymbolicIntervalGradient(sym, LΛ, UΛ)
end

# Symbolic forward_act
function forward_act(::ReluVal, L::Layer{Id}, input::SymbolicIntervalMask)
    n_node = size(input.sym.Up, 1)
    LΛ = push!(input.LΛ, trues(n_node))
    UΛ = push!(input.UΛ, trues(n_node))
    return SymbolicIntervalGradient(input.sym, LΛ, UΛ)
end

function get_max_smear_index(nnet::Network, input::Hyperrectangle, LG::Matrix, UG::Matrix)

    smear(lg, ug, r) = sum(max.(abs.(lg), abs.(ug))) * r

    ind = argmax(smear.(eachcol(LG), eachcol(UG), input.radius))
    monotone = all(>(0), LG[:, ind] .* UG[:, ind]) # NOTE should it be >= 0 instead?

    return ind, monotone
end


text:
using FaultTree
using DD
using Test

@testset "FaultTree1" begin
    x = ftevent(:x) & ftevent(:y)
    println(x)
    @test x.op == Symbol(:AND)
end

@testset "FaultTree1" begin
    x = ftevent(:x)
    y = ftevent(:y)
    z = ftevent(:z)
    expr = ~(x * y + z)
    println(expr)
end

@testset "todot1" begin
    x = ftevent(:x)
    y = ftevent(:y)
    z = ftevent(:z)
    expr = ~((x & y) | z | x)
    expr = ftkofn(2, expr, expr, expr)
    println(todot(expr))
end

@testset "bdd1" begin
    x = ftevent(:x)
    y = ftevent(:y)
    z = ftevent(:z)
    expr = ~((x & y) | z | x)
    forest = BDDForest{Int,Int,Int}(FullyReduced())
    ft = bdd!(forest, expr)
    println(todot(forest, ft))
end

@testset "bdd2" begin
    x = ftevent(:x)
    y = ftevent(:y)
    z = ftevent(:z)
    expr = ftkofn(2, x, y, z)
    forest = BDDForest{Int,Int,Int}(FullyReduced())
    ft = bdd!(forest, expr)
    println(todot(forest, ft))
end

@testset "prob1" begin
    x = ftevent(:x)
    y = ftevent(:y)
    z = ftevent(:z)
    expr = ftand(x, y, z)
    forest = BDDForest{Int,Int,Int}(FullyReduced())
    ft = bdd!(forest, expr)
    env = Dict([:x => 0.1, :y => 0.5, :z => 0.8])
    cache = Dict{AbstractDDNode{Int,Int},Float64}()
    result = fteval!(ft, env, cache)
    println(result)
    @test isapprox(result, 0.1*0.5*0.8)
end

@testset "prob2" begin
    x = ftevent(:x)
    y = ftevent(:y)
    z = ftevent(:z)
    expr = ftor(x, y, z)
    ft, = bdd(expr)
    env = Dict([:x => 0.1, :y => 0.5, :z => 0.8])
    result = fteval(ft, env)
    println(result)
    @test isapprox(result, 1-(1-0.1)*(1-0.5)*(1-0.8))
end

@testset "prob3" begin
    x = ftevent(:x)
    y = ftevent(:y)
    z = ftevent(:z)
    expr = ftkofn(2, x, y, z)
    ft, = bdd(expr)
    env = Dict([:x => 0.1, :y => 0.5, :z => 0.8])
    result = fteval(ft, env)
    println(result)
    @test isapprox(result, 0.1*0.5*(1-0.8) + (1-0.1)*0.5*0.8 + 0.1*(1-0.5)*0.8 + 0.1*0.5*0.8)
end

@testset "mcs1" begin
    x = ftevent(:x)
    y = ftevent(:y)
    z = ftevent(:z)
    expr = ftand(x, y, z)
    forest = BDDForest{Int,Int,Int}(FullyReduced())
    ft = bdd!(forest, expr)
    cache = Dict{AbstractDDNode{Int,Int},Vector{Vector{Symbol}}}()
    result = ftmcs!(ft, cache)
    println(result)
    @test result == [[:x, :y, :z]]
end

@testset "mcs2" begin
    x = ftevent(:x)
    y = ftevent(:y)
    z = ftevent(:z)
    expr = ftor(x, y, z)
    forest = BDDForest{Int,Int,Int}(FullyReduced())
    ft = bdd!(forest, expr)
    cache = Dict{AbstractDDNode{Int,Int},Vector{Vector{Symbol}}}()
    result = ftmcs!(ft, cache)
    println(result)
    @test result == [[:x], [:y], [:z]]
end

@testset "mcs3" begin
    x = ftevent(:x)
    y = ftevent(:y)
    z = ftevent(:z)
    expr = ftkofn(2, x, y, z)
    forest = BDDForest{Int,Int,Int}(FullyReduced())
    ft = bdd!(forest, expr)
    cache = Dict{AbstractDDNode{Int,Int},Vector{Vector{Symbol}}}()
    result = ftmcs!(ft, cache)
    println(result)
    @test result == [[:x, :y], [:x, :z], [:y, :z]]
end

@testset "mcs4" begin
    x = [ftevent(Symbol("x$i")) for i = 1:10]
    expr = (x[9] | (x[2] & x[5] & x[6]) | (x[1] & x[10]) | (x[2] & x[6]))
    forest = BDDForest{Int,Int,Int}(FullyReduced())
    ft = bdd!(forest, expr)
    # println(todot(forest, ft))
    cache = Dict{AbstractDDNode{Int,Int},Vector{Vector{Symbol}}}()
    result = ftmcs!(ft, cache)
    println(result)
    @test result == [[:x1, :x10], [:x2, :x6], [:x9]]
end

@testset "probx" begin
    x = ftevent(:x)
    y = ftevent(:y)
    z = ftevent(:z)
    expr = ftand(x, y, z)
    forest = BDDForest{Int,Int,Int}(FullyReduced())
    ft = bdd!(forest, expr)
    env = Dict([:x => 0.0, :y => 0.0, :z => 0.0])
    cache = Dict{AbstractDDNode{Int,Int},Float64}()
    result = fteval!(ft, env, cache)
    println(result)
    @test isapprox(result, 0.0)
end

@testset "dprob1" begin
    x = ftevent(:x)
    y = ftevent(:y)
    z = ftevent(:z)
    expr = ftand(x, y, z)
    forest = BDDForest{Int,Int,Int}(FullyReduced())
    ft = bdd!(forest, expr)
    env = Dict([:x => 0.1, :y => 0.5, :z => 0.8])
    denv = Dict([:x => 1.0, :y => 0.0, :z => 0.0])
    cache = Dict{AbstractDDNode{Int,Int},Float64}()
    dcache = Dict{AbstractDDNode{Int,Int},Float64}()
    result = fteval!(ft, env, denv, cache, dcache)
    println(result)
    @test isapprox(result, 0.5*0.8)
end

@testset "dprob2" begin
    x = ftevent(:x)
    y = ftevent(:y)
    z = ftevent(:z)
    expr = ftand(x, y, z)
    forest = BDDForest{Int,Int,Int}(FullyReduced())
    ft = bdd!(forest, expr)
    env = Dict([:x => 0.1, :y => 0.5, :z => 0.8])
    denv = Dict([:x => 0.0, :y => 1.0, :z => 0.0])
    cache = Dict{AbstractDDNode{Int,Int},Float64}()
    dcache = Dict{AbstractDDNode{Int,Int},Float64}()
    result = fteval!(ft, env, denv, cache, dcache)
    println(result)
    @test isapprox(result, 0.1*0.8)
end

@testset "fteval21b" begin
    x = ftevent(:x)
    y = ftevent(:y)
    z = ftevent(:z)
    expr = ftand(x, y, z)
    forest = BDDForest{Int,Int,Int}(FullyReduced())
    ft = bdd!(forest, expr)
    lam1 = 0.001
    lam2 = 0.01
    env = Dict(:x => exp(-lam1), :y => exp(-lam1), :z => exp(-lam2))
    cache = Dict{AbstractDDNode{Int,Int},Float64}()
    result = fteval!(ft, env, cache)
    @test isapprox(result, exp(-2*lam1-lam2))
end

@testset "fteval2b" begin
    x = ftevent(:x)
    y = ftevent(:y)
    z = ftevent(:z)
    expr = ftand(x, y, z)
    forest = BDDForest{Int,Int,Int}(FullyReduced())
    ft = bdd!(forest, expr)
    lam1 = 0.001
    lam2 = 0.01
    env0 = Dict(
        :x => Float64[1-exp(-lam1)  0            0            0;
                      exp(-lam1)    1-exp(-lam1) 0            0;
                      0             0            1-exp(-lam1) 0;
                      0             0            exp(-lam1)   1-exp(-lam1)],
        :y => Float64[1-exp(-lam1)  0            0            0;
                      exp(-lam1)    1-exp(-lam1) 0            0;
                      0             0            1-exp(-lam1) 0;
                      0             0            exp(-lam1)   1-exp(-lam1)],
        :z => Float64[1-exp(-lam2)  0            0            0;
                      0             1-exp(-lam2) 0            0;
                      exp(-lam2)    0            1-exp(-lam2) 0;
                      0             exp(-lam2)   0            1-exp(-lam2)],
                      )
    env1 = Dict(
        :x => Float64[exp(-lam1)  0            0            0;
                      -exp(-lam1)    exp(-lam1) 0            0;
                      0             0            exp(-lam1) 0;
                      0             0            -exp(-lam1)   exp(-lam1)],
        :y => Float64[exp(-lam1)  0            0            0;
                      -exp(-lam1)    exp(-lam1) 0            0;
                      0             0            exp(-lam1) 0;
                      0             0            -exp(-lam1)   exp(-lam1)],
        :z => Float64[exp(-lam2)  0            0            0;
                      0             exp(-lam2) 0            0;
                      -exp(-lam2)    0            exp(-lam2) 0;
                      0             -exp(-lam2)   0            exp(-lam2)],
                      )
    cache = Dict{AbstractDDNode{Int,Int},Vector{Float64}}()
    cache[ddval!(forest, 0)] = Float64[0,0,0,0]
    cache[ddval!(forest, 1)] = Float64[1,0,0,0]
    result = ftevalgen!(ft, env0, env1, cache)
    println(result)
    @test isapprox(result, [exp(-2*lam1-lam2), -2*exp(-2*lam1-lam2), -exp(-2*lam1-lam2), 2*exp(-2*lam1-lam2)])
end

@testset "mcs5" begin
    x = [ftevent(Symbol("x", i)) for i = 1:10]
    expr = (x[9] | (x[2] & x[5] & x[6]) | (x[1] & x[10]) | (x[2] & x[6]))
    result = ftmcs(expr)
    println(result)
    @test result == [[:x1, :x10], [:x2, :x6], [:x9]]
end

@testset "prob4" begin
    x = ftevent(:x)
    y = ftevent(:y)
    z = ftevent(:z)
    expr = ftkofn(2, x, y, z)
    env = Dict([:x => 0.1, :y => 0.5, :z => 0.8])
    result = fteval(expr, env)
    println(result)
    @test isapprox(result, 0.1*0.5*(1-0.8) + (1-0.1)*0.5*0.8 + 0.1*(1-0.5)*0.8 + 0.1*0.5*0.8)
end

@testset "mcs6" begin
    x = ftevent(:x)
    y = ftevent(:y)
    z = ftevent(:z)
    expr = x * y + z
    expr = ftkofn(2, expr, expr, expr)
    result = ftmcs(expr)
    println(result)
    @test result == [[:x, :y], [:z]]
end

@testset "ftevent1" begin
    x = [ftevent("x", i) for i = 1:3]
    y = [ftevent("y", i) for i = 1:3]
    z = [ftevent("z", i) for i = 1:3]
    ft = @. x * y + z
    expr = ftkofn(2, ft[1], ft[2], ft[3])
    top, forest = bdd(expr)
    # println(todot(forest, top))
    result = ftmcs(expr)
    println(result)
    @test result == [[:x1, :x2, :y1, :y2], [:x1, :x3, :y1, :y3], [:x1, :y1, :z2], [:x1, :y1, :z3], [:x2, :x3, :y2, :y3], [:x2, :y2, :z1], [:x2, :y2, :z3], [:x3, :y3, :z1], [:x3, :y3, :z2], [:z1, :z2], [:z1, :z3], [:z2, :z3]]
end

@testset "ftevent2" begin
    x = [ftevent("x", i) for i = 1:3]
    y = [ftevent("y", i) for i = 1:2]
    z = ftevent("z")
    s1 = x[1] + z
    s2 = x[2] + z
    s3 = x[3] + z
    ft = ftkofn(2, s1, s2, s3) + y[1] * y[2]
    top, forest = bdd(ft)
    # println(todot(forest, top))
    result = ftmcs(ft)
    println(result)
    @test result == [[:x1, :x2], [:x1, :x3], [:x2, :x3], [:y1, :y2], [:z]]
end

@testset "ftevent3" begin
    x = [ftevent("motor failure ", i) for i = ["arm1", "arm2", "arm3"]] # 3 motors
    y = [ftevent("sensor failure ", i) for i = ["left", "right"]] # 2 sensors
    z = ftevent("battery failure") # battery failure
    s1 = x[1] + z
    s2 = x[2] + z
    s3 = x[3] + z
    ft = ftkofn(2, s1, s2, s3) + y[1] * y[2]
    top, forest = bdd(ft)
    # println(todot(forest, top))
    result = ftmcs(ft)
    println(result)
    @test result == [[Symbol("battery failure")], [Symbol("motor failure arm1"), Symbol("motor failure arm2")], [Symbol("motor failure arm1"), Symbol("motor failure arm3")], [Symbol("motor failure arm2"), Symbol("motor failure arm3")], [Symbol("sensor failure left"), Symbol("sensor failure right")]]
end


text:
module CodeTools

using Lazy, LNR

include("utils.jl")
include("eval.jl")
include("module.jl")
include("summaries.jl")

include("completions.jl")
include("doc.jl")

end # module


text:
using Test
export checkgrad

#=
macro checkgrad(f, params, atol=1e-3, cuda=true)
    quote
        checkgrad(() -> $(esc(f)), $(map(esc,params)...), atol=$atol, cuda=$cuda)
    end
end
=#

function checkgrad(f, inputs::Var...; atol=1e-3, cuda=true)
    params = tuple(filter(isparam, [inputs...])...)
    y = f()
    zerograd!.(params)
    gradient!(y)
    gxs1 = gradient.(params)
    gxs2 = map(params) do p
        x = p.data
        gx2 = similar(x)
        for k = 1:length(x)
            xk = x[k]
            x[k] = xk + 1e-3
            y1 = copy(f().data) # In case y == x
            x[k] = xk - 1e-3
            y2 = copy(f().data)
            x[k] = xk
            gx2[k] = sum(y1-y2) / 2e-3
        end
        gx2
    end

    for (gx1,gx2) in zip(gxs1,gxs2)
        @test maximum(abs,gx1-gx2) <= atol
    end

    if cuda && CUDA.AVAILABLE
        foreach(x -> todevice!(x,0), inputs)
        d_y = f()
        zerograd!.(params)
        gradient!(d_y)
        d_gxs = gradient.(params)

        @test maximum(abs,y.data-Array(d_y.data)) <= atol
        for (gx,d_gx) in zip(gxs1,d_gxs)
            @test maximum(abs,gx-Array(d_gx)) <= atol
        end
        foreach(x -> todevice!(x,-1), inputs)
    end
end


text:
@testset "MPC: Policy for Quadruped" begin
    # ## Simulation
    s = ContactImplicitMPC.get_simulation("quadruped", "flat_2D_lc", "flat");
    model = s.model
    env = s.env

    # ## Reference Trajectory
    ref_traj = deepcopy(get_trajectory(s.model, s.env,
        joinpath(module_dir(), "src/dynamics/quadruped/gaits/gait2.jld2"),
        load_type = :split_traj_alt))
    update_friction_coefficient!(ref_traj, model, env);

    H = ref_traj.H
    h = ref_traj.h

    # ## MPC setup
    N_sample = 5
    H_mpc = 10
    h_sim = h / N_sample
    H_sim = 1000
    κ_mpc = 2.0e-4

    obj = TrackingObjective(model, env, H_mpc,
        q = [Diagonal(1e-2 * [1.0; 0.02; 0.25; 0.25 * ones(model.nq-3)]) for t = 1:H_mpc],
        u = [Diagonal(3e-2 * ones(model.nu)) for t = 1:H_mpc],
        γ = [Diagonal(1.0e-100 * ones(model.nc)) for t = 1:H_mpc],
        b = [Diagonal(1.0e-100 * ones(model.nc * friction_dim(env))) for t = 1:H_mpc]);

    p = ci_mpc_policy(ref_traj, s, obj,
        H_mpc = H_mpc,
        N_sample = N_sample,
        κ_mpc = κ_mpc,
        mode = :configuration,
        n_opts = NewtonOptions(
            solver = :lu_solver,
            r_tol = 3e-4,
            max_iter = 5,
            ),
        mpc_opts = CIMPCOptions(),
        ip_opts = InteriorPointOptions(
                        undercut = 5.0,
                        γ_reg = 0.1,
                        κ_tol = κ_mpc,
                        r_tol = 1.0e-8,
                        diff_sol = true,
                        solver = :empty_solver,
                        ),
        )

    # ## Initial conditions
    q1_sim, v1_sim = initial_conditions(ref_traj);

    # ## Simulator
    sim = simulator(s, H_sim, h=h_sim, policy=p);

    # ## Simulate
    status = simulate!(sim, q1_sim, v1_sim);

    @test status
    qerr, uerr, γerr, berr = ContactImplicitMPC.tracking_error(ref_traj, sim.traj, N_sample, idx_shift=[1])
    @test qerr < 0.0201 * 1.5 # 0.0201
    @test uerr < 0.0437 * 1.5 # 0.0437
    @test γerr < 0.374 * 1.5 # 0.374
    @test berr < 0.0789 * 1.5 # 0.0789
    qerr > 0.0201 * 1.2 && @warn "mild regression on q tracking: current tracking error = $qerr, nominal tracking error = 0.0201"
    uerr > 0.0437 * 1.2 && @warn "mild regression on u tracking: current tracking error = $uerr, nominal tracking error = 0.0437"
    γerr > 0.3740 * 1.2 && @warn "mild regression on γ tracking: current tracking error = $γerr, nominal tracking error = 0.374"
    berr > 0.0789 * 1.2 && @warn "mild regression on b tracking: current tracking error = $berr, nominal tracking error = 0.0789"
end

# @testset "MPC quadruped: long trajectory with structured newton solver" begin
	# s = get_simulation("quadruped", "flat_2D_lc", "flat")
	# model = s.model
	# env = s.env

	# ref_traj = deepcopy(ContactImplicitMPC.get_trajectory(s.model, s.env,
	#     joinpath(module_dir(), "src/dynamics/quadruped/gaits/gait2.jld2"),
	#     load_type = :split_traj_alt))

	# # time
	# H = ref_traj.H
	# h = ref_traj.h
	# N_sample = 5
	# H_mpc = 10
	# h_sim = h / N_sample
	# H_sim = 1000

	# # barrier parameter
	# κ_mpc = 1.0e-4

	# obj_mpc = ContactImplicitMPC.quadratic_objective(model, H_mpc,
	#     q = [Diagonal(1e-2 * [1.0; 0.02; 0.25; 0.25 * ones(model.nq-3)]) for t = 1:H_mpc+2],
	#     v = [Diagonal(0.0 * ones(model.nq)) for t = 1:H_mpc],
	#     u = [Diagonal(3e-2 * ones(model.nu)) for t = 1:H_mpc-1])
	# obj_mpc = TrackingObjective(model, env, H_mpc,
	#     q = [Diagonal(1e-2 * [1.0; 0.02; 0.25; 0.25 * ones(model.nq-3)]) for t = 1:H_mpc],
	#     u = [Diagonal(3e-2 * ones(model.nu)) for t = 1:H_mpc],
	#     γ = [Diagonal(1.0e-100 * ones(model.nc)) for t = 1:H_mpc],
	#     b = [Diagonal(1.0e-100 * ones(model.nc * friction_dim(env))) for t = 1:H_mpc])

	# p = ci_mpc_policy(ref_traj, s, obj_mpc,
	#     H_mpc = H_mpc,
	#     N_sample = N_sample,
	#     κ_mpc = κ_mpc,
	# 	mode = :configuration,
	# 	newton_mode = :structure,
	#     n_opts = NewtonOptions(
	# 		solver = :lu_solver,
	# 		r_tol = 3e-4,
	# 		max_iter = 5,
	# 		),
	#     mpc_opts = CIMPCOptions(),
	# 	ip_opts = InteriorPointOptions(
	# 		max_iter = 100,
	# 		verbose = false,
	# 		r_tol = 1.0e-4,
	# 		κ_tol = 1.0e-4,
	# 		diff_sol = true,
	# 		solver = :empty_solver,
	# 		),
	#     )
    # # initial configurations
    # q1_sim = ref_traj.q[2]
    # v1_sim = (ref_traj.q[2] - ref_traj.q[1]) ./ h

    # sim = ContactImplicitMPC.simulator(s, H_sim, h=h_sim, policy=p)

    # # simulator
    # @test status = ContactImplicitMPC.simulate!(sim, q1_sim, v1_sim);
    # ref_traj = deepcopy(ref_traj)
	# qerr, uerr, γerr, berr = tracking_error(ref_traj, sim.traj, N_sample, idx_shift = [1])
	# @test qerr < 0.0202 * 1.5 # 0.0202
	# @test uerr < 0.0437 * 1.5 # 0.0437
	# @test γerr < 0.3780 * 1.5 # 0.3780
	# @test berr < 0.0799 * 1.5 # 0.0799
	# qerr > 0.0201 * 1.2 && @warn "mild regression on q tracking: current tracking error = $qerr, nominal tracking error = 0.0202"
	# uerr > 0.0437 * 1.2 && @warn "mild regression on u tracking: current tracking error = $uerr, nominal tracking error = 0.0437"
	# γerr > 0.3790 * 1.2 && @warn "mild regression on γ tracking: current tracking error = $γerr, nominal tracking error = 0.3780"
	# berr > 0.0798 * 1.2 && @warn "mild regression on b tracking: current tracking error = $berr, nominal tracking error = 0.0799"

# end


text:
#=
Copyright (c) 2018-2022 Chris Coey, Lea Kapelevich, and contributors

This Julia package Hypatia.jl is released under the MIT license; see LICENSE
file in the root directory or at https://github.com/chriscoey/Hypatia.jl

step distance search helpers for s,z
=#

mutable struct StepSearcher{T <: Real}
    min_prox::T
    prox_bound::T
    use_sum_prox::Bool
    alpha_sched::Vector{T}

    skzk::Vector{T}
    nup1::T
    cone_times::Vector{Float64}
    cone_order::Vector{Int}
    prev_sched::Int

    function StepSearcher{T}(
        model::Models.Model{T};
        min_prox::T = T(0.01),
        prox_bound::T = T(0.99),
        use_sum_prox::Bool = false,
        alpha_sched::Vector{T} = default_alpha_sched(T),
        ) where {T <: Real}
        cones = model.cones
        searcher = new{T}()
        searcher.min_prox = min_prox
        searcher.prox_bound = prox_bound
        searcher.use_sum_prox = use_sum_prox
        searcher.alpha_sched = alpha_sched
        searcher.skzk = zeros(T, length(cones))
        searcher.nup1 = T(model.nu + 1)
        searcher.cone_times = zeros(length(cones))
        searcher.cone_order = collect(1:length(cones))
        searcher.prev_sched = 0
        return searcher
    end
end

default_alpha_sched(T::Type{<:Real}) = T[
    0.9999, 0.999, 0.99, 0.97, 0.95, 0.9, 0.85, 0.8, 0.7, 0.6, 0.5,
    0.3, 0.1, 0.05, 0.01, 0.005, 0.001, 0.0005]

# backwards search on alphas in alpha schedule
function search_alpha(
    point::Point{T},
    model::Models.Model{T},
    stepper::Stepper{T};
    sched::Int = start_sched(stepper, stepper.searcher)
    ) where {T <: Real}
    searcher = stepper.searcher

    while sched <= length(searcher.alpha_sched)
        alpha = searcher.alpha_sched[sched]
        # update ztsk only in stepper.temp
        update_stepper_points(alpha, point, stepper, true)

        # NOTE updates cone points and grad
        if check_cone_points(model, stepper)
            searcher.prev_sched = sched
            return alpha
        end
        sched += 1
    end

    searcher.prev_sched = sched
    return zero(T)
end

# fallback starts at first alpha in schedule
start_sched(stepper::Stepper, searcher::StepSearcher) = 1

function check_cone_points(
    model::Models.Model{T},
    stepper::Stepper{T};
    ) where {T <: Real}
    searcher = stepper.searcher
    cand = stepper.temp
    skzk = searcher.skzk
    cones = model.cones
    prox_bound = searcher.prox_bound
    min_prox = searcher.min_prox
    use_sum_prox = searcher.use_sum_prox

    taukap_cand = cand.tau[] * cand.kap[]
    (min(cand.tau[], cand.kap[], taukap_cand) < eps(T)) && return false

    for k in eachindex(cones)
        skzk[k] = dot(cand.primal_views[k], cand.dual_views[k])
        (skzk[k] < eps(T)) && return false
    end
    mu_cand = (sum(skzk) + taukap_cand) / searcher.nup1

    taukap_rel = taukap_cand / mu_cand
    taukap_prox = abs(taukap_rel - 1)
    if (mu_cand < eps(T)) || (taukap_rel < min_prox) || (taukap_prox > prox_bound)
        return false
    end
    prox_bound_sqr = abs2(prox_bound)
    for k in eachindex(cones)
        nu_k = Cones.get_nu(cones[k])
        skzkmu = skzk[k] / mu_cand
        if (skzkmu < min_prox * nu_k) || (abs2(skzkmu - nu_k) > prox_bound_sqr * nu_k)
            return false
        end
    end

    # order the cones by how long it takes to check proximity condition and
    # iterate in that order, to improve efficiency
    cone_order = searcher.cone_order
    sortperm!(cone_order, searcher.cone_times, initialized = true) # stochastic

    sum_prox = taukap_prox
    rtmu = sqrt(mu_cand)
    irtmu = inv(rtmu)
    for k in cone_order
        cone_k = cones[k]
        start_time = time()
        Cones.load_point(cone_k, cand.primal_views[k], irtmu)
        Cones.load_dual_point(cone_k, cand.dual_views[k])
        Cones.reset_data(cone_k)

        in_prox_k = false
        if Cones.is_feas(cone_k) && Cones.is_dual_feas(cone_k) &&
            Cones.check_numerics(cone_k)
            prox_k = Cones.get_proximity(cone_k, rtmu, use_sum_prox)
            if !isnan(prox_k) && (prox_k < prox_bound)
                sum_prox += prox_k
                if !use_sum_prox || (sum_prox < prox_bound)
                    in_prox_k = true
                end
            end
        end
        searcher.cone_times[k] = time() - start_time
        in_prox_k || return false
    end

    return true
end


text:
#=
    "low level" way i.e. without Platform
    Scenario constructor, loading measurements from file
    `sim` for single scenario, scenario array, pairs array
    `fit` for scenario array, pairs array
    `mc` for single scenario, scenario array, pairs array
    statistics
=#
using HetaSimulator, Plots

################################## Model Upload ###########################################
# heta_update_dev()

platform = load_platform("$HetaSimulatorDir/cases/story_1", rm_out = false);
model = platform.models[:nameless]
# parameters(model) # def_parameters
# events_active(model) # def_events_active
# events_save(model) # def_events_save
# observables(model) # def_observations

################################## Single Simulation ######################################

Scenario(model; tspan = (0., 200.)) |> sim |> plot
scn0 = Scenario(model; tspan = (0., 200.))
sim(scn0, parameters_upd = [:k1=>0.01]) |> plot
Scenario(model; saveat = 0:10:100) |> sim |> plot
Scenario(model; saveat = 0:10:100, tspan = (0., 50.)) |> sim |> plot
Scenario(model; saveat = 0:10:100, tspan = (0., 500.)) |> sim |> plot
Scenario( # throw error
    model;
    tspan = (0., 500.),
    events_active=[:sw1=>false, :ss1 => false],
    events_save=[:sw1=>(true,true), :ss1=>(true,true)]
    ) |> sim |> plot
Scenario(
    model;
    tspan = (0., 500.),
    events_active=[:sw1=>false],
    events_save=[:sw1=>(true,true)]
    ) |> sim |> plot
Scenario(
    model;
    tspan = (0., 500.),
    events_active=[:sw1=>true],
    events_save=[:sw1=>(false,false)]
    ) |> sim |> plot
sim(Scenario(model, tspan = (0., 10.), parameters=[:k1=>1e-3]), parameter_upd = [:k1=>1e-3])

### single scenario sim()
scn1 = Scenario(model; tspan = (0., 200.), saveat = [0.0, 150., 250.]);
sim(scn1) |> plot
sim(scn1; parameters_upd=[:k1=>0.01]) |> plot

scn2 = Scenario(
    model;
    tspan = (0., 200.),
    events_active=[:sw1=>false],
    parameters = [:k2 => 0.001, :k3 => 0.02]
    );
sim(scn2) |> plot

scn3 = Scenario(
    model;
    tspan = (0., 250.),
    events_active=[:sw1=>false],
    parameters = [:k2 => 0.1]
    );
sim(scn3) |> plot

### sim sequentially
sim.([scn1, scn2, scn3]) .|> plot
### sim together
sim([scn1, scn2, scn3]) |> plot
sim([:x => scn1, :y=>scn2, :z=>scn3]) |> plot
sim([:x => scn1, :y=>scn2, :z=>scn3]; parameters_upd=[:k1=>0.01]) |> plot

### load measurements from CSV
#measurements_csv = read_measurements("$HetaSimulatorDir/cases/story_1/measurements.csv")
measurements_csv = read_measurements("$HetaSimulatorDir/cases/story_1/measurements_no_scope.csv")
measurements_xlsx = read_measurements("$HetaSimulatorDir/cases/story_1/measurements.xlsx")
scn4 = Scenario(model; parameters = [:k2=>0.001, :k3=>0.04], saveat = [0.0, 50., 150., 250.]);
add_measurements!(scn4, measurements_csv; subset = [:scenario => :dataone])

### fit many scenarios
res1 = fit([:x=>scn2, :y=>scn3, :z=>scn4], [:k1=>0.1,:k2=>0.2,:k3=>0.3])
res2 = fit([scn2, scn3, scn4], [:k1=>0.1,:k2=>0.2,:k3=>0.3])
sim(scn3, parameters_upd = optim(res2))

# sim all scenarios
sol = sim([:c1=>scn1, :c2=>scn2, :c3=>scn3, :c4=>scn4]);
plot(sol)

# plot selected observables
plot(sol; vars=[:a,:c])
# plot selected observables without data
plot(sol; vars=[:a,:c], show_measurements=false)

################################## Monte-Carlo Simulations  #####################

mc_scn1 = Scenario(
    model;
    tspan = (0., 200.),
    parameters = [:k1=>0.01],
    saveat = [50., 80., 150.]
    );

mc_scn2 = Scenario(
    model;
    tspan = (0., 200.),
    parameters = [:k1=>0.02],
    saveat = [50., 100., 200.]
    );

mc_scn3 = Scenario(
    model;
    tspan = (0., 200.),
    parameters = [:k1=>0.03],
    events_active=[:sw1 => false]
    );

# single MC Simulation
mcsim1 = mc(mc_scn1, [:k1=>Uniform(1e-3,1e-2), :k2=>Normal(1e-3,1e-4), :k3=>Normal(1e-4,1e-5)], 1000)
plot(mcsim1)

# multi MC Simulation
mcsim2 = mc(
    [:mc1=>mc_scn1,:mc2=>mc_scn2,:mc3=>mc_scn3],
    [:k1=>0.01, :k2=>Normal(1e-3,1e-4), :k3=>Uniform(1e-4,1e-2)],
    1000
  )
plot(mcsim2)

mcsim3 = mc(
    [mc_scn1, mc_scn2, mc_scn3],
    [:k1=>0.01, :k2=>LogNormal(1e-3,1e-4), :k3=>Uniform(1e-4,1e-2)],
    1000
  )
plot(mcsim3)

################################## Monte-Carlo Statistics  #####################

# mean
timestep_mean(mcsim1,2)
timepoint_mean(mcsim1,80)

# median
timestep_median(mcsim1,2)
timepoint_median(mcsim1,80)

# meanvar
timestep_meanvar(mcsim1,2)
timepoint_meanvar(mcsim1,80)

# meancov
timestep_meancov(mcsim1,2,3)
timepoint_meancov(mcsim1,80.,150.)

# meancor
timestep_meancor(mcsim1,2,3)
timepoint_meancor(mcsim1,80.,150.)

# !!!quantile
timestep_quantile(mcsim1,0.95,2)
timepoint_quantile(mcsim1,0.95,80.)

# full time steps statistics
timeseries_steps_mean(mcsim1) # Computes the mean at each time step
timeseries_steps_median(mcsim1) # Computes the median at each time step
timeseries_steps_quantile(mcsim1,0.95) # Computes the quantile q at each time step
timeseries_steps_meanvar(mcsim1) # Computes the mean and variance at each time step
timeseries_steps_meancov(mcsim1) # Computes the covariance matrix and means at each time step
timeseries_steps_meancor(mcsim1) # Computes the correlation matrix and means at each time step

# Ensemble Summary
ens = EnsembleSummary(mcsim1;quantiles=[0.05,0.95])
plot(ens)

################################## Monte-Carlo Parallel #####################
#=
using Distributed
addprocs(2)
@everywhere using HetaSimulator

mcscn1 = Scenario(model; tspan = (0., 200.), parameters = [:k1=>0.01], saveat = [50., 80., 150.]);
mcsim0 = mc(mcscn1, [:k2=>Normal(1e-3,1e-4), :k3=>Normal(1e-4,1e-5)], 20)
mcsim1 = mc(mcscn1, [:k2=>Normal(1e-3,1e-4), :k3=>Normal(1e-4,1e-5)], 150, parallel_type=EnsembleDistributed())
=#


text:
# This file defines the elements required to create an app with PackageCompiler
# - module
# - julia_main function called at the execution of the app

module Triangulation

include("lib.jl")

# disable logging for levels Info and below
disable_logging(Logging.Info)

function julia_main()::Cint
    try
        main(ARGS)
        return 0
    catch
        return 1
    end
end

end # module

# CLI call from a terminal: call the main function with CLI args
if abspath(PROGRAM_FILE) == @__FILE__
    Triangulation.julia_main()
end
;


text:
println("\n\n\nStarting runtests.jl ...")

using ExtremelyRandomizedTrees
using FactCheck

doPlot = false

if doPlot
	using PyPlot
end

doDummyClassification = true
doDummyRegression = true
doClassification = true
doRegression1to1 = true
doRegression1to2 = true
doRegression2to2 = true
testIndivProbPrediction = true

if doDummyClassification
    data = [-2 -1 1 2]
    labels = [1 1 2 2]
    a = ExtraTrees(float(data), labels)

    testdata = linspace(-5,5,100)'

    result, votes = predict(a, testdata, returnvotes = true)
    if doPlot
        figure()
        plot(testdata', [votes' result'])
    end
end

if doDummyRegression
    trainingData = collect(-10:10)'
    trainingLabels = trainingData+100
    testData = trainingData

    a = ExtraTrees(trainingData, trainingLabels; regression = true)
    result = predict(a,testData)

    if doPlot
        figure()
        plot(testData', result', "b-")
    end
end

if doClassification
    trainingData = hcat(rand(2,100), rand(2,100)+1, rand(2,100)+2)
    trainingLabels = hcat(ones(1,100), 2*ones(1,100), 3*ones(1,100))

    a = ExtraTrees(trainingData, trainingLabels)

    X = [x for x in -5:0.1:10, y in -5:0.1:10]
    Y = [y for x in -5:0.1:10, y in -5:0.1:10]

    result = predict(a, [vec(X) vec(Y)]');

    if doPlot
        figure()
        imshow(reshape(result, size(X)))
    end
end

####################################
# univariate regression
# 1 feature, 1 target
####################################

if doRegression1to1
    trainingData = linspace(0,2*pi,100)'
    testData = linspace(0,2*pi,1000)'
    trainingLabels = 10 + sin(trainingData) + (0.3*rand(size(trainingData))-0.15)

    a = ExtraTrees(trainingData, trainingLabels; regression = true)
    result = predict(a,testData)

	if doPlot
        figure()
		plot(trainingData, trainingLabels,"b.", hold = true)
		plot(testData', result',"g-", hold = true)
		title("univariate regression - 1 feature, 1 target")
	end
end

if testIndivProbPrediction
    trainingData = hcat(rand(2,100), rand(2,100)+1, rand(2,100)+2)
    trainingLabels = hcat(ones(1,100), 2*ones(1,100), 3*ones(1,100))

    a = ExtraTrees(trainingData, trainingLabels)

    X = [x for x in -5:0.1:10, y in -5:0.1:10]
    Y = [y for x in -5:0.1:10, y in -5:0.1:10]

    result, votes = predict(a, [vec(X) vec(Y)]', returnvotes = true)
    indivvotes = Array(Any,3)
    for i = 1:3
        indivvotes[i] = predict(a, [vec(X) vec(Y)]'; votesfor=i)
        #assert(all(votes[i,:]/sum(votes[i,:]) == indivvotes[i]))
    end
	if doPlot
        figure()
		for i = 1:3
			subplot(1,3,i); imshow(reshape(indivvotes[i], size(X)))
		end
	end
end

####################################
# multivariate regression
# 1 feature, 2 target values
####################################

if doRegression1to2
    trainingData = linspace(0,2*pi,100)'
    testData = linspace(0,2*pi,1000)'
    trainingLabels = vcat(3 + sin(trainingData) + (0.3*rand(size(trainingData))-0.15),
        10 + trainingData.^2 + (0.3*rand(size(trainingData))-0.15))

    a = ExtraTrees(trainingData, trainingLabels; regression = true)
    result = predict(a,testData)

	if doPlot
        figure()
		plot(trainingData, trainingLabels[1,:],"b.", hold = true)
		plot(trainingData, trainingLabels[2,:],"b.", hold = true)
		plot(testData', result[1,:]',"g-", hold = true)
		plot(testData', result[2,:]',"g-", hold = true)
		title("multivariate regression - 1 feature, 2 target values")
	end
end

####################################
# multivariate regression
# 2 features, 2 target values
####################################

if doRegression2to2
    trainingData = linspace(0,2*pi,100)'
    testData = linspace(0,2*pi,1000)'
    trainingTargets = vcat(10 + sin(trainingData) + (0.3*rand(size(trainingData))-0.15),
    10 + cos(trainingData) + (0.3*rand(size(trainingData))-0.15));
    testTargets = vcat(10 + sin(testData),
    10 + cos(testData));

    a = ExtraTrees(trainingData, trainingTargets; regression = true)
    result = predict(a,testData)

    assert(maximum(abs(result-testTargets))<0.3)
end

println("  ... runtests.jl done!")


text:
# To be implemented by user.

"""
    ASTInterface.Simulation

Parent type for user simulation.
"""
abstract type Simulation <: AbstractSimulation end

"""
    ASTInterface.reset!

Resets simulation.
"""
function reset!(sim::Simulation)::Nothing end

"""
    ASTInterface.environment

Returns Environment object constructed in simulation.
Environment is an alias of Dict{Symbol, Sampleable}.
"""
function environment(sim::Simulation)::Environment end

"""
    ASTInterface.observe

Returns quasi-normalized observation of simulation.
"""
function observe(sim::Simulation)::Vector{<:Real} end

"""
    ASTInterface.step!

Steps simulation given an EnvironmentValue object.
EnvironmentValue is an alias of Dict{Symbol, Any}.
"""
function step!(sim::Simulation, x::EnvironmentValue)::Nothing end

"""
    ASTInterface.isterminal

Checks whether simulation has finished due to time limit or terminal state, independent of event status.
"""
function isterminal(sim::Simulation)::Bool end

"""
    ASTInterface.isevent

Checks whether simulation is in an event state.
"""
function isevent(sim::Simulation)::Bool end

"""
    ASTInterface.distance

Returns custom metric of distance to event. For best results, metric should depend only on current state.
"""
function distance(sim::Simulation)::Real end

"""
    ASTInterface.flatten

Flattens environment variable to quasi-normalized array.
"""
function flatten(distribution::Any, value::Any)::Vector{<:Real} end

"""
    ASTInterface.unflatten

Reconstructs environment variable from quasi-normalized array.
"""
function unflatten(distribution::Any, array::Vector{<:Real})::Any end


text:
macro test_da_approx_eq(da1, da2)
    quote
        v1 = $(esc(da1))
        v2 = $(esc(da2))
        na = ismissing.(v1)
        @test na == ismissing.(v2)
        defined = (!).(na)
        if any(defined)
            @test isapprox(v1[defined], v2[defined], nans = true)
        end
    end
end

@testset "Reducedim" begin
    # Test for fast unit stride BitArray functions
    function test_any()
        bits = bitrand(64*5)
        for i = 1:length(bits), j = i:length(bits)
            v = false
            for k = i:j
                v |= bits[k]
            end
            Base.Test.@test DataArrays._any(bits, i, j) == v
        end
    end
    test_any()

    function test_count()
        bits = bitrand(64*5)
        for i = 1:length(bits), j = i:length(bits)
            v = 0
            for k = i:j
                v += bits[k]
            end
            Base.Test.@test DataArrays._count(bits, i, j) == v
        end
    end
    test_count()

    # mapslices from Base, hacked to work for these cases
    function safe_mapslices(f::Function, A::AbstractArray{T}, region, skipmissing) where T
        dims = intersect(region, 1:ndims(A))
        if isempty(dims)
            if skipmissing
                naval = f(Missings.T(T)[], 1)
                A = copy(A)
                A[ismissing.(A)] = isempty(naval) ? missing : naval[1]
            end
            return A
        end

        if isempty(dims)
            return map(f,A)
        end

        dimsA = collect(size(A))
        ndimsA = ndims(A)
        alldims = collect(1:ndimsA)

        otherdims = setdiff(alldims, dims)

        idx = Vector{Any}(ndimsA)
        fill!(idx, 1)
        Asliceshape = tuple(dimsA[dims]...)
        itershape   = tuple(dimsA[otherdims]...)
        for d in dims
            idx[d] = 1:size(A,d)
        end

        r1 = f(copy(reshape(A[idx...], Asliceshape)); skipmissing=skipmissing)

        # determine result size and allocate
        Rsize = copy(dimsA)
        # TODO: maybe support removing dimensions
        if !isa(r1, AbstractArray) || ndims(r1) == 0
            r2 = similar(A, T, 1)
            r2[1] = r1
            r1 = r2
        end
        Rsize[dims] = [size(r1)...; ones(Int,max(0,length(dims)-ndims(r1)))]
        R = similar(r1, tuple(Rsize...))

        ridx = Vector{Any}(ndims(R))
        fill!(ridx, 1)
        for d in dims
            ridx[d] = 1:size(R,d)
        end

        R[ridx...] = r1

        first = true

        for idxs = CartesianIndices(itershape)
            if first
                first = false
            else
                ia = [idxs.I...]
                idx[otherdims] = ia
                ridx[otherdims] = ia
                try
                    R[ridx...] = f(copy(reshape(A[idx...], Asliceshape)); skipmissing=skipmissing)
                catch e
                    if (isa(e, ErrorException) && e.msg == "Reducing over an empty array is not allowed.") || (isa(e, ArgumentError) && e.msg == "reducing over an empty collection is not allowed")

                        R[ridx...] = missing
                    else
                        println(typeof(e))
                        rethrow(e)
                    end
                end
            end
        end

        return R
    end

    myvarzm(x; skipmissing::Bool=false) = var(x; mean=0, skipmissing=skipmissing)
    myvar1m(x; skipmissing::Bool=false) = var(x; mean=1, skipmissing=skipmissing)

    for Areduc in (DataArray(rand(3, 4, 5, 6)),
                   DataArray(rand(3, 4, 5, 6), rand(3, 4, 5, 6) .< 0.2))
        for skipmissing = (false, true)
            for region in Any[
                1, 2, 3, 4, 5, (1, 2), (1, 3), (1, 4), (2, 3), (2, 4), (3, 4),
                (1, 2, 3), (1, 3, 4), (2, 3, 4), (1, 2, 3, 4)]
                # println("region = $region, skipmissing = $skipmissing")

                outputs = Any[DataArray(fill(NaN, length.(Base.reduced_indices(indices(Areduc), region))))]
                hasmissings = any(ismissing, Areduc)
                if hasmissings && !skipmissing
                    # Should throw an error reducing to non-DataArray
                    @test_throws MissingException sum!(outputs[1].data, Areduc; skipmissing=skipmissing)
                else
                    # Should be able to reduce to non-DataArray
                    push!(outputs, outputs[1].data)
                end

                for r in outputs
                    @test_da_approx_eq sum!(r, Areduc; skipmissing=skipmissing) safe_mapslices(sum, Areduc, region, skipmissing)
                    @test_da_approx_eq prod!(r, Areduc; skipmissing=skipmissing) safe_mapslices(prod, Areduc, region, skipmissing)
                    if !hasmissings
                        @test_da_approx_eq maximum!(r, Areduc; skipmissing=skipmissing) safe_mapslices(maximum, Areduc, region, skipmissing)
                        @test_da_approx_eq minimum!(r, Areduc; skipmissing=skipmissing) safe_mapslices(minimum, Areduc, region, skipmissing)
                    end
                    @test_da_approx_eq mean!(r, Areduc; skipmissing=skipmissing) safe_mapslices(mean, Areduc, region, skipmissing)
                end

                @test_da_approx_eq sum(Areduc, region; skipmissing=skipmissing) safe_mapslices(sum, Areduc, region, skipmissing)
                @test_da_approx_eq prod(Areduc, region; skipmissing=skipmissing) safe_mapslices(prod, Areduc, region, skipmissing)
                @test_da_approx_eq maximum(Areduc, region; skipmissing=skipmissing) safe_mapslices(maximum, Areduc, region, skipmissing)
                @test_da_approx_eq minimum(Areduc, region; skipmissing=skipmissing) safe_mapslices(minimum, Areduc, region, skipmissing)
                @test_da_approx_eq mean(Areduc, region; skipmissing=skipmissing) safe_mapslices(mean, Areduc, region, skipmissing)

                if region != 5
                    @test_da_approx_eq var(Areduc, region; skipmissing=skipmissing) safe_mapslices(var, Areduc, region, skipmissing)
                    @test_da_approx_eq var(Areduc, region; mean=0, skipmissing=skipmissing) safe_mapslices(myvarzm, Areduc, region, skipmissing)
                    for r in outputs
                        @test_da_approx_eq var(Areduc, region; mean=fill!(r, 1), skipmissing=skipmissing) safe_mapslices(myvar1m, Areduc, region, skipmissing)
                    end
                end
            end
        end
    end

    # Test missing-skipping behavior for maximum
    a = @data([missing missing; 3 4])
    @test isequal(maximum(a, 1; skipmissing=true), [3 4])
    @test isequal(maximum!(zeros(1, 2), a; skipmissing=true), [3 4])

    # Maximum should give an missing in the output if all values along dimension are missing
    @test isequal(maximum(a, 2; skipmissing=true), @data([missing 4])')
    # Maximum should refuse to reduce to a non-DataArray
    @test_throws MissingException maximum!(zeros(2, 1), a; skipmissing=true)
end


text:
module JuliaGrid

export runpf
export runopf
export runmg
export runse

using SparseArrays, LinearAlgebra, SuiteSparse
using PrettyTables, Printf
using HDF5, CSV, XLSX
using Random
using JuMP, Ipopt, GLPK
using LightGraphs, SimpleWeightedGraphs

### Includes
include("system/input.jl")
include("system/routine.jl")
include("system/results.jl")
include("system/headers.jl")

include("flow/flowdc.jl")
include("flow/flowac.jl")
include("flow/optimaldc.jl")
include("flow/measurements.jl")
include("flow/flowfunctions.jl")

include("estimation/estimatedc.jl")
include("estimation/estimatepmu.jl")
include("estimation/estimateac.jl")
include("estimation/estimatefunctions.jl")

### Power Flow

function runpf(
    args...;
    max::Int64 = 100,
    stop::Float64 = 1.0e-8,
    reactive::Int64 = 0,
    solve::String = "",
    save::String = "",
)
    path = loadpath(args)
    system, num, info = loadsystem(path)
    settings = pfsettings(args, max, stop, reactive, solve, save, system, num)

    if settings.algorithm == "dc"
        results = rundcpf(system, num, settings, info)
    else
        results = runacpf(system, num, settings, info)
    end

    return results, system, info
end

### Optimal Power Flow
function runopf(
    args...;
    save::String = "",
)

    path = loadpath(args)
    system, num, info = loadsystem(path)
    settings = opfsettings(args, save, system, num)

    if settings.algorithm == "dc"
        results = rundcopf(system, num, settings, info)
    else

    end

    return results, system, info
end

### Measurement Generator
function runmg(
    args...;
    runflow::Int64 = 1,
    max::Int64 = 100,
    stop::Float64 = 1.0e-8,
    reactive::Int64 = 0,
    solve::String = "",
    save::String = "",
    pmuset = [],
    pmuvariance = [],
    legacyset = [],
    legacyvariance = [],
)

    path = loadpath(args)
    system, numsys, info = loadsystem(path)
    measurements, num = loadmeasurement(path, system, numsys; pmuvar = pmuvariance, legvar = legacyvariance, runflow = runflow)
    settings = gesettings(num, pmuset, pmuvariance, legacyset, legacyvariance, runflow, save)

    if settings.runflow == 1
        pfsettings = gepfsettings(max, stop, reactive, solve)
        acflow = runacpf(system, numsys, pfsettings, info)
        info = rungenerator(system, measurements, num, numsys, settings, info; flow = acflow)
    else
        info = rungenerator(system, measurements, num, numsys, settings, info)
    end

    return measurements, system, info
end

### State Estimation
function runse(
    args...;
    max::Int64 = 100,
    stop::Float64 = 1.0e-8,
    start = [],
    bad = [],
    lav = [],
    observe = [],
    covariance::Int64 = 0,
    solve::String = "",
    save::String = "",
)

    if loadse(args)
        path = loadpath(args)
        system, numsys, info = loadsystem(path)
        measurements, num = loadmeasurement(path, system, numsys)
    else
        system, numsys, measurements, num, info = loadsedirect(args)
    end
    settings = sesettings(args, system, max, stop, start, bad, lav, observe, covariance, solve, save)

    if settings.algorithm == "dc"
        results = rundcse(system, measurements, num, numsys, settings, info)
    elseif settings.algorithm == "pmu"
        results = runpmuse(system, measurements, num, numsys, settings, info)
    elseif settings.algorithm == "nonlinear"
        results = runacse(system, measurements, num, numsys, settings, info)
    end

    return results, measurements, system, info
end

end # JuliaGrid

max_stars_repo_name: mcosovic/JuliaGrid

max_issues_repo_name: mcosovic/JuliaGrid

max_forks_repo_name: mcosovic/JuliaGrid


text: Dict{String, Any}("group" => Dict{String, Any}("more" => Any[Dict{String, Any}("value" => "42", "type" => "integer"), Dict{String, Any}("value" => "42", "type" => "integer")], "answer" => Dict{String, Any}("value" => "42", "type" => "integer"), "d" => Dict{String, Any}("value" => "1979-05-27T07:32:12Z", "type" => "datetime")))


text:
## to split the cfg files
## Format is :  "keyword : value"
function input_parameters(file)
    let
        param = []
        open(file,"r") do f
            for line in eachline(f)
                l = lstrip(line)
                if length(l) > 0 && l[1] != '#'
                    res = strip.(split(l,":"))
                    push!(param,res)
                end
            end
        end
        return(param)
    end
end

## parse input file
##
function parse_input(input_param)
    let
        fileParam= "default.cfg"
        saveFold=  "Results/"
        verbose= true

        for pair in input_param
            if pair[1] == "File_Parameters"
                fileParam = pair[2]
            elseif pair[1] == "Folder_Results"
                saveFold= pair[2]
            elseif pair[1] == "Display_Screen_Results" && pair[2] == "false"
                verbose= false
            end
        end
        res= inputCfg(fileParam , saveFold , verbose)
        return(res)
    end
end

## parse parameter file for the subarray constraints
## No consistency check yet...
function parse_parameters(input_param)
    let
        arraycfg= 0
        obs= 0
        sub= 0
        wei= 0
        ga= 0

        Array_Configuration_File= "alma.cfg"
        Result_Folder= "."
        Display_Verbose= true
        Observatory_Latitude= -23.0262015
        Source_Declination= -30
        Source_Hour_Angle = -1
        Subarray_Number= 4

        ## array of nsubarray elements
        Pads_Per_Subarray= []
        Weight_Subarray= []
        Subarray_Name= []
        Spatial_Resolution= []
        Maximum_Recoverable_Scale=
        Elongation= []
        Sidelobe_Level= []
        Weight_Spatial_Resolution= []
        Weight_Maximum_Recoverable_Scale= []
        Weight_Elongation= []
        Weight_Sidelobe_Levels= []

        ## GA parameters
        Number_Iterations = 100
        Population_Size= 100
        Termination_Condition= false
        Threshold= -0.05
        Mutation_Rate= 0.05
        Tournament_Size= 5
        Number_Elitism= 5

        for pair in input_param
            if pair[1] == "Array_Configuration_File"
                Array_Configuration_File= pair[2]
            elseif pair[1] == "Result_Folder"
                Result_Folder= pair[2]
            elseif pair[1] == "Display_Verbose"
                Display_Verbose=  pair[2]=="true" ? true : false
            elseif pair[1] == "Observatory_Latitude"
                Observatory_Latitude= parse(Float64, pair[2])
            elseif pair[1] == "Source_Declination"
                Source_Declination= parse(Float64, pair[2])
            elseif pair[1] == "Source_Hour_Angle"
                Source_Hour_Angle= parse(Float64, pair[2])
            elseif pair[1] == "Subarray_Number"
                Subarray_Number= parse(Int, pair[2])
            elseif pair[1] == "Number_Iterations"
                Number_Iterations= parse(Int, pair[2])
            elseif pair[1] == "Population_Size"
                Population_Size= parse(Int, pair[2])
            elseif pair[1] == "Termination_Condition"
                Termination_Condition= pair[2]=="true" ? true : false
            elseif pair[1] == "Threshold"
                Threshold=parse(Float64, pair[2])
            elseif pair[1] == "Mutation_Rate"
                Mutation_Rate=parse(Float64, pair[2])
            elseif pair[1] == "Tournament_Size"
                Tournament_Size= parse(Int, pair[2])
            elseif pair[1] == "Number_Elitism"
                Number_Elitism= parse(Int, pair[2])
            elseif pair[1] == "Pads_Per_Subarray"
                Pads_Per_Subarray= map(x->(v = tryparse(Int,x);
                        v==nothing ? 0.0 : v),split(pair[2],","))
            elseif pair[1] == "Weights_Subarray"
                  Weight_Subarray= map(x->(v = tryparse(Float64,x);
                        v==nothing ? 0.0 : v),split(pair[2],","))
            elseif pair[1] == "Subarray_Name"
                Subarray_Name= split(pair[2],",")
                for i in 1:length(Subarray_Name)
                  Subarray_Name[i]= strip(Subarray_Name[i])
                end
            elseif pair[1] == "Spatial_Resolution"
                Spatial_Resolution= map(x->(v = tryparse(Float64,x);
                        v==nothing ? 0.0 : v),split(pair[2],","))
            elseif pair[1] == "Maximum_Recoverable_Scale"
                Maximum_Recoverable_Scale= map(x->(v = tryparse(Float64,x);
                        v==nothing ? 0.0 : v),split(pair[2],","))
            elseif pair[1] == "Elongation"
                Elongation= map(x->(v = tryparse(Float64,x);
                        v==nothing ? 0.0 : v),split(pair[2],","))
            elseif pair[1] == "Sidelobe_Level"
               Sidelobe_Level= map(x->(v = tryparse(Float64,x);
                        v==nothing ? 0.0 : v),split(pair[2],","))
            elseif pair[1] == "Weight_Spatial_Resolution"
                Weight_Spatial_Resolution= map(x->(v = tryparse(Float64,x);
                        v==nothing ? 0.0 : v),split(pair[2],","))
            elseif pair[1] == "Weight_Maximum_Recoverable_Scale"
                Weight_Maximum_Recoverable_Scale= map(x->(v = tryparse(Float64,x);
                        v==nothing ? 0.0 : v),split(pair[2],","))
            elseif pair[1] == "Weight_Elongation"
                Weight_Elongation= map(x->(v = tryparse(Float64,x);
                        v==nothing ? 0.0 : v),split(pair[2],","))
            elseif pair[1] == "Weight_Sidelobe_Levels"
                Weight_Sidelobe_Levels=map(x->(v = tryparse(Float64,x);
                        v==nothing ? 0.0 : v),split(pair[2],","))
            end
        end

        ## trim leading/trailing spaces...NOT YET

        arrcfg = CSV.read(convert(String,Array_Configuration_File), datarow=4 , header=["X" , "Y", "Z" , "diam" , "name"] ,  ignorerepeated=true, delim= " ")

        subrange = []
        start= 1
            for npad in Pads_Per_Subarray
                push!(subrange,start:start+npad-1)
                start= start+npad
            end

        Antenna_Number= size(arrcfg,1)

        ### setting the struct.
        obs= observation(Array_Configuration_File, Result_Folder, Display_Verbose, Observatory_Latitude , Source_Declination,Source_Hour_Angle,
            Antenna_Number ,Subarray_Number)
        sub= subarrayParameters(Pads_Per_Subarray, Subarray_Name , subrange, Spatial_Resolution,
            Maximum_Recoverable_Scale , Elongation, Sidelobe_Level)
        wei= weight(Weight_Subarray,Weight_Spatial_Resolution,Weight_Maximum_Recoverable_Scale,
            Weight_Elongation,Weight_Sidelobe_Levels)
        ga= GA(Number_Iterations , Population_Size,Termination_Condition ,Threshold ,
            Mutation_Rate , Tournament_Size ,Number_Elitism)

        ### population setting
        res= cfg(arrcfg , obs , sub , wei , ga)

   return(res)
end
end

## Normalize the weights
##
function normalize_weight!(cfg)
   k1= sum(cfg.wei.Weight_Spatial_Resolution) + sum(cfg.wei.Weight_Maximum_Recoverable_Scale) +
     sum(cfg.wei.Weight_Elongation) + sum(cfg.wei.Weight_Sidelobe_Levels)
   cfg.wei.Weight_Spatial_Resolution= cfg.wei.Weight_Spatial_Resolution ./ k1
   cfg.wei.Weight_Maximum_Recoverable_Scale= cfg.wei.Weight_Maximum_Recoverable_Scale ./ k1
   cfg.wei.Weight_Elongation= cfg.wei.Weight_Elongation ./ k1
   cfg.wei.Weight_Sidelobe_Levels= cfg.wei.Weight_Sidelobe_Levels ./ k1

  k2= sum(cfg.wei.Weight_Subarray)
  cfg.wei.Weight_Subarray= cfg.wei.Weight_Subarray ./ k2
end

## main function to read the cfg and check it.
function read_cfg(inpfile, check=true)
    res= input_parameters(inpfile)
    cfg= parse_parameters(res)
    normalize_weight!(cfg)

    if cfg.obs.Display_Verbose
      @printf("## Input Parameters for GASS \n")
      @printf("### Configuration file: %s \n", cfg.obs.Array_Configuration_File)
      @printf("### Result folder: %s \n", cfg.obs.Result_Folder)
      @printf("### Obs. Latitude: %3.3f \n", cfg.obs.Observatory_Latitude)
      @printf("### Source Declination: %3.1f \n", cfg.obs.Source_Declination)
      @printf("### HA: %3.1f \n", cfg.obs.Source_Hour_Angle)
      @printf("### Antenna number: %d \n", cfg.obs.Antenna_Number)
      @printf("### Subarray number: %d \n", cfg.obs.Subarray_Number)
      @printf("##\n")

      @printf("## Subarray Parameters\n")
      println("### Pads per subarray: ",cfg.sub.Pads_Per_Subarray)
      println("### Name: ",cfg.sub.Subarray_Name)
      println("### AR: ",cfg.sub.Spatial_Resolution)
      println("### MRS: ",cfg.sub.Maximum_Recoverable_Scale)
      println("### elongation: ",cfg.sub.Elongation)
      println("### sidelobe: ",cfg.sub.Sidelobe_Level)
      @printf("##\n")

      @printf("## GA parameters\n")
      @printf("### Iterations: %d \n",cfg.ga.Number_Iterations)
      @printf("### Population size: %d \n",cfg.ga.Population_Size)
      @printf("### Mutation rate: %3.3f \n",cfg.ga.Mutation_Rate)
      @printf("### Tournament size: %d \n",cfg.ga.Tournament_Size)
      @printf("### Elitism: %d \n",cfg.ga.Number_Elitism)
      @printf("##\n")

      @printf("## Weights\n")
      println("### Subarray weights: ",cfg.wei.Weight_Subarray)
      println("### AR weights: ",cfg.wei.Weight_Spatial_Resolution)
      println("### MRS weights: ",cfg.wei.Weight_Maximum_Recoverable_Scale)
      println("### elongation weights: ",cfg.wei.Weight_Elongation)
      println("### sidelobe weights: ",cfg.wei.Weight_Sidelobe_Levels)
      @printf("##\n")
    end

    if check
      check_consistency(cfg)
    end

    return(cfg)
end

function check_consistency(cfg::cfg)
    nsubpads= sum(cfg.sub.Pads_Per_Subarray)
    npads= nrow(cfg.arr)
    if nsubpads != npads
        error("##Error: Subarray pads are not equal to the total of pads.")
    end
    return(true)
end

max_stars_repo_name: bosscha/GASS

max_issues_repo_name: bosscha/GASS

max_forks_repo_name: bosscha/GASS


text:
function fib(n::Int)::BigInt
    # https://www.nayuki.io/page/fast-fibonacci-algorithms
    # fast doubling algorithm
    # Returns the tuple (fib(n), fib(n+1)).
    function F(n::Int)::Tuple{BigInt,BigInt}
	      if n == 0
		        (BigInt(0), BigInt(1))
	      else
            k,m = divrem(n,2)
		        Fk,Fkp1 = F(k)
		        c = Fk * (2 * Fkp1 - Fk)
		        d = Fk^2 + Fkp1^2
		        if m == 0
			          (c, d)
		        else
			          (d, c + d)
            end
        end
    end

    if n < 0
        F(n)[1] * ( n%2 == 0 ? -1 : +1 )
    else
        F(n)[1]
    end
end

max_stars_repo_name: UlrichBerntien/Codewars-Katas

max_issues_repo_name: UlrichBerntien/Codewars-Katas

max_forks_repo_name: UlrichBerntien/Codewars-Katas


text:
__precompile__()

module ExpFunctions

struct FunctionSet
   setmodel!::Function
   initstate!::Function
   setcontrol!::Function
   observe!::Function
   reward::Function
   evaluate::Function
end

export FunctionSet

end


text:
push!(LOAD_PATH, "Modules")
using Colors
using JuliRay

φ=MathConstants.φ

Base.@irrational ° 0.0174532925199432957692369076848861271344 (big(pi)/big(180))

function POLYGON2(vertices,O,R;r1=1.0,r2=0.05,r3=0.025)
    n=length(vertices)
    C=O+R*normalize(+(vertices...)/n-O)
    V=rgbColor(csgUnion((p->Sphere(p,r2)).(vertices)),RGB(0.1,0.1,0.1))
    E=rgbColor(csgUnion([Arc(vertices[i],O+R*normalize((vertices[i]+vertices[mod(i,n)+1])/2-O),vertices[mod(i,n)+1],r3) for i in 1:n]),RGB(0.1,0.1,0.1))
    if(R<10^2)
        F=csgClip(
            rgbftColor(
                Sphere(O,R),RGB(1,0.5,1),FT(0.1,0.1)
            ),
            csgIntersection(
                [
                    (v=JuliRay.NormalVector(O,vertices[i],vertices[mod(i,n)+1]);s=sign(dot(v,C-O));Cylinder(O,O+s*R*v,R))
                for i ∈ 1:n]
            )
        )
    else
        F=rgbftColor(Polygon(vertices),RGB(1,0.5,1),FT(0.1,0.1))
    end
    return csgUnion(V,E,F)
end

function Mirror(q,p1,p2,p3)
    e₃=NormalVector(p1,p2,p3)
    e₁=OrthogonalVector(e₃)
    e₂=cross(e₃,e₁)
    A=hcat(e₁,e₂,e₃)*transpose(hcat(e₁,e₂,-e₃))
    fixedpoint=(p1+p2+p3)/3
    b=fixedpoint-A*fixedpoint
    return A*q+b
end

function Put(F,i,O)
    p1=O
    n=length(F)
    p2,p3=F[mod(i-1,n)+1],F[mod(i,n)+1]
    return (p->Mirror(p,p1,p2,p3)).(F)
end

function f(x)
    ifelse(x>0,exp(-1/x),0)
end
function g(x)
    f(2x/√3)/(f(2x/√3)+f(2/√3-2x/√3))
end
function Smooth(a,b,x)
    g((x-a)/(b-a))
end

## Tetrahedron
h=1/3 # 単位球面に内接させた際の内接球半径
n=3 # n角形
m=4 # m面体
M=60 # アニメーションの刻み数

r1=√(1-h^2) # n角形の外接円半径
for i ∈ 1:M
    θ=(1-Smooth(0,1,2(i-1)/M)+Smooth(1,2,2(i-1)/M))*π/2 # パラメータθは0→1→0の順で滑らかに変化
    R=1/cos(θ)
    O=[0,0,tan(θ)]
    N=O+[0,0,R]
    H=O-[0,0,√(R^2-r1^2)]

    F=[[[r1*cos(2π*i/n),r1*sin(2π*i/n),0.0]+H for i in 1:n] for i ∈ 1:m] # 最初のn角形の頂点
    F[2]=Put(F[1],1,O) # 以下はn角形を付け加える操作
    F[3]=Put(F[2],2,O)
    F[4]=Put(F[2],3,O)

    object=csgUnion([POLYGON2(F[i],O,R,r2=0.03,r3=0.01) for i ∈ 1:m]) # 頂点座標をもとにオブジェクトを生成
    render(csgUnion(object),camera=LngLatCamera(lng=30°,lat=30°,pers=0.2,zoom=0.2,width=640,height=360),name="Ns"*string(m),index=i)
end

## Cube
h=1/√3 # 単位球面に内接させた際の内接球半径
n=4 # n角形
m=6 # m面体
M=60 # アニメーションの刻み数

r1=√(1-h^2) # n角形の外接円半径
for i ∈ 1:M
    θ=(1-Smooth(0,1,2(i-1)/M)+Smooth(1,2,2(i-1)/M))*π/2 # パラメータθは0→1→0の順で滑らかに変化
    R=1/cos(θ)
    O=[0,0,tan(θ)]
    N=O+[0,0,R]
    H=O-[0,0,√(R^2-r1^2)]

    F=[[[r1*cos(2π*i/n),r1*sin(2π*i/n),0.0]+H for i in 1:n] for i ∈ 1:m] # 最初のn角形の頂点
    F[2]=Put(F[1],1,O) # 以下はn角形を付け加える操作
    F[3]=Put(F[2],3,O)
    F[4]=Put(F[1],2,O)
    F[5]=Put(F[4],3,O)
    F[6]=Put(F[5],1,O)

    object=csgUnion([POLYGON2(F[i],O,R,r2=0.03,r3=0.01) for i ∈ 1:m]) # 頂点座標をもとにオブジェクトを生成
    render(object,camera=LngLatCamera(lng=30°,lat=30°,pers=0.2,zoom=0.2,width=640,height=360),name="Ns"*string(m),index=i)
end

## Octahedron
h=1/√3 # 単位球面に内接させた際の内接球半径
n=3 # n角形
m=8 # m面体
M=60 # アニメーションの刻み数

r1=√(1-h^2) # n角形の外接円半径
for i ∈ 1:M
    θ=(1-Smooth(0,1,2(i-1)/M)+Smooth(1,2,2(i-1)/M))*π/2 # パラメータθは0→1→0の順で滑らかに変化
    R=1/cos(θ)
    O=[0,0,tan(θ)]
    N=O+[0,0,R]
    H=O-[0,0,√(R^2-r1^2)]

    F=[[[r1*cos(2π*i/n),r1*sin(2π*i/n),0.0]+H for i in 1:n] for i ∈ 1:m] # 最初のn角形の頂点
    F[2]=Put(F[1],1,O) # 以下はn角形を付け加える操作
    F[3]=Put(F[2],3,O)
    F[4]=Put(F[3],2,O)
    F[5]=Put(F[4],3,O)
    F[6]=Put(F[1],2,O)
    F[7]=Put(F[6],3,O)
    F[8]=Put(F[7],2,O)

    object=csgUnion([POLYGON2(F[i],O,R,r2=0.03,r3=0.01) for i ∈ 1:m]) # 頂点座標をもとにオブジェクトを生成
    render(object,camera=LngLatCamera(lng=30°,lat=30°,pers=0.2,zoom=0.2,width=640,height=360),name="Ns"*string(m),index=i)
end

## Dodecahedron
φ=2cos(π/5)
ξ=2sin(π/5)

h=φ/(√3*ξ) # 単位球面に内接させた際の内接球半径
n=5 # n角形
m=12 # m面体
M=60 # アニメーションの刻み数

r1=√(1-h^2) # n角形の外接円半径
for i ∈ 1:M
    θ=(1-Smooth(0,1,2(i-1)/M)+Smooth(1,2,2(i-1)/M))*π/2 # パラメータθは0→1→0の順で滑らかに変化
    R=1/cos(θ)
    O=[0,0,tan(θ)]
    N=O+[0,0,R]
    H=O-[0,0,√(R^2-r1^2)]

    F=[[[r1*cos(2π*i/n),r1*sin(2π*i/n),0.0]+H for i in 1:n] for i ∈ 1:m] # 最初のn角形の頂点
    F[2]=Put(F[1],4,O) # 以下はn角形を付け加える操作
    F[3]=Put(F[2],1,O)
    F[4]=Put(F[2],2,O)
    F[5]=Put(F[2],3,O)
    F[6]=Put(F[2],5,O)
    F[7]=Put(F[1],2,O)
    F[8]=Put(F[7],5,O)
    F[9]=Put(F[8],1,O)
    F[10]=Put(F[8],2,O)
    F[11]=Put(F[8],3,O)
    F[12]=Put(F[8],4,O)

    object=csgUnion([POLYGON2(F[i],O,R,r2=0.03,r3=0.01) for i ∈ 1:m]) # 頂点座標をもとにオブジェクトを生成
    render(object,camera=LngLatCamera(lng=30°,lat=30°,pers=0.2,zoom=0.2,width=640,height=360),name="Ns"*string(m),index=i)
end

## Icosahedron
φ=2cos(π/5)
ξ=2sin(π/5)

h=φ/(√3*ξ) # 単位球面に内接させた際の内接球半径
n=3 # n角形
m=20 # m面体
M=60 # アニメーションの刻み数

r1=√(1-h^2) # n角形の外接円半径
for i ∈ 1:M
    θ=(1-Smooth(0,1,2(i-1)/M)+Smooth(1,2,2(i-1)/M))*π/2 # パラメータθは0→1→0の順で滑らかに変化
    R=1/cos(θ)
    O=[0,0,tan(θ)]
    N=O+[0,0,R]
    H=O-[0,0,√(R^2-r1^2)]

    F=[[[r1*cos(2π*i/n),r1*sin(2π*i/n),0.0]+H for i in 1:n] for i ∈ 1:m] # 最初のn角形の頂点
    F[2]=Put(F[1],1,O) # 以下はn角形を付け加える操作
    F[3]=Put(F[2],3,O)
    F[4]=Put(F[3],2,O)
    F[5]=Put(F[4],1,O)
    F[6]=Put(F[5],3,O)
    F[7]=Put(F[1],2,O)
    F[8]=Put(F[7],3,O)
    F[9]=Put(F[8],1,O)
    F[10]=Put(F[9],2,O)
    F[11]=Put(F[1],3,O)
    F[12]=Put(F[2],2,O)
    F[13]=Put(F[3],1,O)
    F[14]=Put(F[4],3,O)
    F[15]=Put(F[5],2,O)
    F[16]=Put(F[6],1,O)
    F[17]=Put(F[7],1,O)
    F[18]=Put(F[8],2,O)
    F[19]=Put(F[9],3,O)
    F[20]=Put(F[10],1,O)

    object=csgUnion([POLYGON2(F[i],O,R,r2=0.03,r3=0.01) for i ∈ 1:m]) # 頂点座標をもとにオブジェクトを生成
    render(object,camera=LngLatCamera(lng=30°,lat=30°,pers=0.2,zoom=0.2,width=640,height=360),name="Ns"*string(m),index=i)
end

max_stars_repo_name: hyrodium/JuliRay

max_issues_repo_name: hyrodium/JuliRay

max_forks_repo_name: hyrodium/JuliRay


text:
module Ising2D

using Statistics
using StatsBase
using Printf

# abstract types
include("types.jl")

# utility functions
include("utils2D.jl")
export RandomLattice, RandomMLattice, IsingOptions

# abstract sampler
include("sample.jl")
export sample_trajectory!, sample_trajectory, sample_observables

# observables, including energy
include("observables2D.jl")
export Energy, Magnetization, Clusters

# spin flip sampler with Metropolis rates
include("samplers/Metropolis2D.jl")
export Metropolis

# Wolff sampler
include("samplers/Wolff2D.jl")
export Wolff

end # module end


text:
# This file was used to create the list of deprecations when moving from v0.8.1
# to v0.9.0.

using Gurobi

io = open("deprecated_functions.jl", "w")
print(io, """
const _DEPRECATED_ERROR_MESSAGE = \"\"\"
The C API of Gurobi.jl has been rewritten to expose the complete C API, and
all old functions have been removed.

For example:

    model = Gurobi.Optimizer()
    stat = Gurobi.get_status_code(model.inner)

is now:

    model = Gurobi.Optimizer()
    valueP = Ref{Cint}()
    ret = GRBgetintattr(model, "Status", valueP)
    if ret != 0
        # Do something because the call failed
    end
    stat = valueP[]

The new API is more verbose, but the names and function arguments are now
identical to the C API, documentation for which is available at:
https://www.gurobi.com/documentation/9.0/refman/c_api_details.html

To revert to the old API, use:

    import Pkg
    Pkg.add(Pkg.PackageSpec(name = \"Gurobi\", version = v\"0.8.1\"))

Then restart Julia for the change to take effect.
\"\"\"
""")

exported_names = Base.names(Gurobi; all = false)
for name in Base.names(Gurobi; all = true)
    foo = getfield(Gurobi, name)
    if !(foo isa Function)
        continue
    elseif any(startswith.(Ref(string(foo)), ["#", "@", "_"]))
        continue
    end
    println(io, "$(foo)(args...; kwargs...) = error(_DEPRECATED_ERROR_MESSAGE)")
    if name in exported_names
        println(io, "export $(foo)")
    end
    println(io)
end
close(io)


text:
abstract type MTTranslationModel <: TranslationModel end

struct MTTranslationParam <: EoSParam
    acentricfactor::SingleParam{Float64}
end

@newmodelsimple MTTranslation MTTranslationModel MTTranslationParam

export MTTranslation
function MTTranslation(components::Vector{String}; userlocations::Vector{String}=String[], verbose::Bool=false)
    params = getparams(components, ["properties/critical.csv"]; userlocations=userlocations, verbose=verbose)
    acentricfactor = SingleParam(params["w"],"acentric factor")
    packagedparams = MTTranslationParam(acentricfactor)
    model = MTTranslation(packagedparams, verbose=verbose)
    return model
end

function translation(model::CubicModel,V,T,z,translation_model::MTTranslation)
    Tc = model.params.Tc.values
    Pc = model.params.Pc.values
    ω  = translation_model.params.acentricfactor.values
    c = zeros(typeof(T),length(Tc))
    for i ∈ @comps
        ωi = ω[i]
        Zc = evalpoly(ωi,(0.289,-0.0701,-0.0207))
        β  = -10.2447-28.6312*ωi
        Tci = Tc[i]
        RTp = R̄*Tci/Pc[i]
        t0 = RTp*evalpoly(ωi,(-0.014471,0.067498,-0.084852,0.067298,-0.017366))
        tc = RTp*(0.3074-Zc)
        Tr = T/Tci
        c[i] = t0+(tc-t0)*exp(β*abs(1-Tr))
    end
    return c
end