text:
include("utils.jl")

""" Nominal facility location model. """
function facility_model(c::Matrix, f::Vector)
    n, m = size(c)
    @assert length(f) == n
    model = Model(Gurobi.Optimizer)
    set_optimizer_attribute(model, "OutputFlag", 0)

    # VARIABLES
    @variable(model, x[1:n], Bin)      # Facility locations
    @variable(model, y[1:n, 1:m] >= 0) # Flow variables
    # CONSTRAINTS
    @constraint(model, [j=1:m], sum(y[i, j] for i = 1:n) >= d[j])      # demand constraint
    @constraint(model, [i=1:n], sum(y[i, j] for j = 1:m) <= s[i]*x[i]) # capacity constraint
    # OBJECTIVE
    @objective(model, Min,
            sum(f[j] * x[j] for j = 1:n) + sum(c[i, j] * y[i, j] for i=1:n, j=1:m))
    return model, x, y
end

""" Robust facility location model, with JuMP using the robust counterpart. """
function robust_facility_model(c::Matrix, f::Vector, rho::Real, Gamm::Real)
    n, m = size(c)
    @assert length(f) == n
    model = Model(Gurobi.Optimizer)
    set_optimizer_attribute(model, "OutputFlag", 0)

    # VARIABLES
    @variable(model, x[1:n], Bin)                 # Facility locations
    @variable(model, u[1:n, 1:m] >= 0)            # Flow variables

    # CONSTRAINTS
    for j=1:m
        infdummy = @variable(model)
        normdummy = @variable(model, [1:m])
        y = @variable(model, [1:m])
        @constraint(model, [i=1:m], normdummy[i] >= y[i])
        @constraint(model, [i=1:m], normdummy[i] >= -y[i])
        @constraint(model, [i=1:m], infdummy >= -y[i] - P[j,i])
        @constraint(model, [i=1:m], infdummy >= y[i] + P[j,i])
        @constraint(model, sum(u[i, j] for i = 1:n) >= d[j] + rho*sum(normdummy) + Gamm * infdummy)
    end
    @constraint(model, [i=1:n], sum(u[i, j] for j = 1:m) <= s[i] * x[i]) # capacity constraint
    # OBJECTIVE
    @objective(model, Min,
            sum(f[j] * x[j] for j = 1:n) + sum(c[i, j] * u[i, j] for i=1:n, j=1:m))
    return model, x, u
end

# Nominal model (3.1)
model, x, y = facility_model(c, f)
optimize!(model)
plt = plot_solution(model, x, y)
@show plt

# Robust model (3.2)
rho = 1
Gamm = 5
model, x, y = robust_facility_model(c, f, rho, Gamm)
optimize!(model)
plt = plot_solution(model, x, y)
@show plt

# Different radii of influence
R_Ds = collect(0:0.05:0.75)
R_Ds[1] = 1e-5

# Unscaled P (3.3)
df = DataFrame("R_D" => [], "sumP" => [], "Gamma" => [], "f+c" => [], "f" => [], "c" => [], "nP" => [], "nx" => [])
for k=1:length(R_Ds)
    global P = [0.2*exp(-1/R_Ds[k] * LinearAlgebra.norm(customers[i, :] .- customers[j, :])[1]) for j=1:m, i=1:m];
    global P = (P .>= 0.2*exp(-1/R_Ds[k] .* R_Ds[k])) .* P
    println("Nonzero Ps: $(sum(P .> 0))")
    model, x, y = robust_facility_model(c, f, rho, Gamm)
    optimize!(model)
    push!(df, Dict("R_D" => R_Ds[k],
                     "sumP" => sum(P),
                     "Gamma" => Gamm,
                     "f+c" => objective_value(model),
                     "f" => value(sum(f[j] * x[j] for j = 1:n)),
                     "c" => value(sum(c[i, j] * y[i, j] for i=1:n, j=1:m)),
                     "nP" => sum(P .> 0),
                     "nx" => sum(value.(x))))
end

# Column-wise scaled P (3.7)
df = DataFrame("R_D" => [], "sumP" => [], "Gamma" => [], "f+c" => [], "f" => [], "c" => [], "nP" => [], "nx" => [])
P_nominal = [0.2*exp(-1/R_D .*LinearAlgebra.norm(customers[i, :] .- customers[j, :])[1]) for j=1:m, i=1:m];
P_nominal = (P_nominal .>= 0.2*exp(-1/R_D .* R_D)) .* P_nominal
for k=1:length(R_Ds)
    global P = [0.2*exp(-1/R_Ds[k] * LinearAlgebra.norm(customers[i, :] .- customers[j, :])[1]) for j=1:m, i=1:m];
    global P = (P .>= 0.2*exp(-1/R_Ds[k] .* R_Ds[k])) .* P
    for i=1:m
        global P[i,:] = P[i,:] .* sum(P_nominal[i,:]) ./ sum(P[i,:])
    end
    println("Nonzero Ps: $(sum(P .> 0))")
    model, x, y = robust_facility_model(c, f, rho, Gamm)
    optimize!(model)
    push!(df, Dict("R_D" => R_Ds[k],
                     "sumP" => sum(P),
                     "Gamma" => Gamm,
                     "f+c" => objective_value(model),
                     "f" => value(sum(f[j] * x[j] for j = 1:n)),
                     "c" => value(sum(c[i, j] * y[i, j] for i=1:n, j=1:m)),
                     "nP" => sum(P .> 0),
                     "nx" => sum(value.(x))))
end


text:
module ODETools

using Test

using ForwardDiff

const ϵ = 100*eps()
export ∂t
export ∂tt
export time_derivative

using Gridap.Fields: VectorValue, TensorValue
using Gridap.Fields: return_type
using Gridap.Arrays: get_array

using Gridap: GridapType
using Gridap.Helpers
using Gridap.Algebra: NonlinearSolver
using Gridap.Algebra: LinearSolver
using Gridap.Algebra: NonlinearOperator
using Gridap.Algebra: AffineOperator

export ODEOperator
export AffineODEOperator
export ConstantODEOperator
export SecondOrderODEOperator
export OperatorType
export Nonlinear
export Affine
export Constant
using Gridap.Algebra: residual
using Gridap.Algebra: jacobian
using Gridap.Algebra: symbolic_setup
using Gridap.Algebra: numerical_setup
using Gridap.Algebra: numerical_setup!
import Gridap.Algebra: residual!
import Gridap.Algebra: jacobian!
import Gridap.Algebra: allocate_residual
import Gridap.Algebra: allocate_jacobian
export allocate_cache
export update_cache!
export jacobian!
export jacobian_t!
export jacobian_and_jacobian_t!
export test_ode_operator

export ODESolver
export solve_step!
export test_ode_solver
import Gridap: solve
import Gridap.Algebra: solve!
import Gridap.Algebra: zero_initial_guess

export BackwardEuler
export ForwardEuler
export MidPoint
export ThetaMethod
export RungeKutta
export Newmark
import Gridap.Algebra: fill_entries!

export ODESolution
export test_ode_solution
import Base: iterate

include("DiffOperators.jl")

include("ODEOperators.jl")

include("ODESolvers.jl")

include("ODESolutions.jl")

end #module

ext: jl


text:
struct PSUserModel

    # data for user -
    alphavantage_api_email::String
    alphavantage_api_key::String

    # constructor -
    function PSUserModel(api_email::String, api_key::String)
        new(api_email,api_key)
    end
end

struct PSDataStoreAPICallModel

    # data -
    ticker::String
    dataType::Symbol
    outputsize::Symbol
    apicall::Function

    function PSDataStoreAPICallModel(apicall::Function, ticker::String; dataType::Symbol=:cvs, output::Symbol = :compact)
        this = new(ticker,dataType,output,apicall)
    end
end


text:
@testset "Testing DescriptorZoo exports" begin
using DescriptorZoo, Test
end


text:
"""
    PolicyGraphGenerator

* `fast = false` use fast (but very picky) alternate parser for .pomdp files
* `graph_max_depth` maximum horizon of the generated policy graph. There is no limit by default.
* `graph_max_branch` maximum number of branches of the generated policy graph. Shown will be top in probability. There is no limit by default.
* `graph_min_prob` minimum probability threshold for a branch to be shown in the policy graph defaults to zero
* `graph_filename` output name for the DOT file to be generated
* `pomdp_filename` pomdp file input
* `policy_filename` policy file input
"""
@with_kw struct PolicyGraphGenerator
    fast::Bool = false
    graph_max_depth::Union{Int,Nothing}=nothing
    graph_max_branch::Union{Int,Nothing}=nothing
    graph_min_prob::Union{Float64,Nothing}=nothing
    graph_filename::AbstractString = "policy_graph.dot"
    pomdp_filename::AbstractString = "model.pomdpx"
    policy_filename::AbstractString = "policy.out"
end

function get_graph_generator_options(graphgen::PolicyGraphGenerator)
    options = AbstractString[]
    push!(options, "--policy-graph")
    push!(options, string(graphgen.graph_filename))
    if graphgen.fast
        push!(options, "--fast")
    end
    if graphgen.graph_max_depth != nothing
        push!(options, "--graph-max-depth")
        push!(options, string(graphgen.graph_max_depth))
    end
    if graphgen.graph_max_branch != nothing
        push!(options, "graph-max-branch")
        push!(options, string(graphgen.graph_max_branch))
    end
    if graphgen.graph_min_prob != nothing
        push!(options, "graph-min-prob")
        push!(options, string(graphgen.graph_min_prob))
    end
    return options
end

"""
    generate_graph(graphgen::PolicyGraphGenerator)

Generate a policy graph, see PolicyGraphGenerator to see the available options.
"""
function generate_graph(graphgen::PolicyGraphGenerator)
    options_list = get_graph_generator_options(graphgen)
    polgraph() do polgraph_path
        run(`$polgraph_path $(graphgen.pomdp_filename) --policy-file $(graphgen.policy_filename) $options_list`)
    end
end


text:
using StatsBase
using Base.Test

# rle & inverse_rle

z = [1, 1, 2, 2, 2, 3, 1, 2, 2, 3, 3, 3, 3]
(vals, lens) = rle(z)
@test vals == [1, 2, 3, 1, 2, 3]
@test lens == [2, 3, 1, 1, 2, 4]
@test inverse_rle(vals, lens) == z

z = [true, true, false, false, true, false, true, true, true]
vals, lens = rle(z)
@test vals == [true, false, true, false, true]
@test lens == [2, 2, 1, 1, 3]
@test inverse_rle(vals, lens) == z

# findat

findat([2, 4, 6, 8, 10], [8, 6, 4]) == [4, 3, 2]
findat([2, 4, 6, 8, 10], [2, 3, 4]) == [1, 0, 2]

# indicatormat

I = bool([0 1 0 0 0;
          1 0 0 0 1;
          0 0 1 1 0])

x = [2, 1, 3, 3, 2]
@test indicatormat(x, 3) == I
@test full(indicatormat(x, 3; sparse=true)) == I

x = ["b", "a", "c", "c", "b"]
@test indicatormat(x) == I
@test full(indicatormat(x; sparse=true)) == I


text:
#=
given data (xᵢ, yᵢ), find a polynomial p to solve
minimize    ∑ᵢℓ(p(xᵢ), yᵢ)
subject to  ρⱼ × dᵏp/dtⱼᵏ ≥ 0 ∀ t ∈ D
where
- dᵏp/dtⱼᵏ is the kᵗʰ derivative of p in direction j,
- ρⱼ determines the desired sign of the derivative,
- D is a domain such as a box or an ellipsoid,
- ℓ is a convex loss function.
see e.g. Chapter 8 of thesis by G. Hall (2018)
=#

import DelimitedFiles
import Distributions

struct ShapeConRegrJuMP{T <: Real} <: ExampleInstanceJuMP{T}
    X::Matrix{T}
    y::Vector{T}
    deg::Int
    use_wsos::Bool # use WSOS cone formulation, else SDP formulation
    use_L1_obj::Bool # in objective function use L1 norm, else L2 norm
    use_monotonicity::Bool # if true add monotonicity constraints, else don't
    use_convexity::Bool # if true add convexity constraints, else don't
    is_fit_exact::Bool
end

function ShapeConRegrJuMP{Float64}(
    data_name::Symbol,
    args...)
    Xy = DelimitedFiles.readdlm(joinpath(@__DIR__, "data", "$data_name.txt"))
    (X, y) = (Xy[:, 1:(end - 1)], Xy[:, end])
    # TODO assert X data is on the domain [-1, 1]^n
    return ShapeConRegrJuMP{Float64}(X, y, args...)
end

function ShapeConRegrJuMP{Float64}(
    n::Int,
    num_points::Int,
    func::Symbol,
    signal_ratio::Real,
    args...)
    X = rand(Distributions.Uniform(-1, 1), num_points, n)
    f = shapeconregr_data[func]
    y = Float64[f(X[p, :]) for p in 1:num_points]
    if !iszero(signal_ratio)
        noise = randn(num_points)
        noise .*= norm(y) / sqrt(signal_ratio) / norm(noise)
        y .+= noise
    end
    return ShapeConRegrJuMP{Float64}(X, y, args...)
end

function build(inst::ShapeConRegrJuMP{T}) where {T <: Float64}
    (X, y, deg) = (inst.X, inst.y, inst.deg)
    n = size(X, 2)
    num_points = size(X, 1)
    mono_dom = PolyUtils.BoxDomain{T}(-ones(size(X, 2)), ones(size(X, 2)))
    conv_dom = mono_dom
    mono_profile = ones(Int, size(X, 2))
    conv_profile = 1

    # setup interpolation (not actually using FreeDomain, just need points here)
    halfdeg = div(deg + 1, 2)
    free_dom = PolyUtils.FreeDomain{T}(n)
    (U, points, Ps, V) = PolyUtils.interpolate(free_dom,
        halfdeg, calc_V = true)
    F = qr!(Array(V'), ColumnNorm()) # TODO reuse QR parts
    V_X = PolyUtils.make_chebyshev_vandermonde(X, 2halfdeg)
    X_points_polys = F \ V_X'

    model = JuMP.Model()
    JuMP.@variable(model, regressor[1:U])
    JuMP.@variable(model, z)
    JuMP.@objective(model, Min, z)

    # objective epigraph
    norm_vec = y - X_points_polys' * regressor
    if inst.use_L1_obj || (num_points <= U)
        obj_K = (inst.use_L1_obj ? MOI.NormOneCone : MOI.SecondOrderCone)(
            1 + num_points)
        JuMP.@constraint(model, vcat(z, norm_vec) in obj_K)
    else
        # using L2 norm objective and number of samples exceeds variables,
        # so use qr trick to reduce dimension
        coef_mat = zeros(num_points, U + 1)
        for (i, expr_i) in enumerate(norm_vec)
            for (c, v) in JuMP.linear_terms(expr_i)
                coef_mat[i, JuMP.index(v).value] = c
            end
            coef_mat[i, end] = JuMP.constant(expr_i)
        end
        coef_R = qr(coef_mat).R
        JuMP.@constraint(model, vcat(z, coef_R * vcat(regressor, 1)) in
            MOI.SecondOrderCone(2 + U))
    end

    # monotonicity
    if inst.use_monotonicity
        gradient_halfdeg = div(deg, 2)
        (mono_U, mono_points, mono_Ps) =
            PolyUtils.interpolate(mono_dom, gradient_halfdeg)
        univ_chebs_der = [PolyUtils.calc_univariate_chebyshev(
            mono_points[:, i], 2halfdeg, calc_gradient = true) for i in 1:n]

        for j in 1:n
            iszero(mono_profile[j]) && continue

            univ_chebs_g = [univ_chebs_der[i][(i == j) ? 2 : 1] for i in 1:n]
            V_g = PolyUtils.make_product_vandermonde(univ_chebs_g,
                PolyUtils.n_deg_exponents(n, 2halfdeg))
            scal = inv(maximum(abs, V_g) / 10)
            if scal < 1e-7
                @warn("model is numerically challenging to set up", maxlog = 1)
            end
            lmul!(scal, V_g)
            g_points_polys = F \ V_g'

            gradient_interp = mono_profile[j] * g_points_polys' * regressor

            if inst.use_wsos
                JuMP.@constraint(model, gradient_interp in
                    Hypatia.WSOSInterpNonnegativeCone{T, T}(mono_U, mono_Ps))
            else
                psd_vars = []
                for (r, Pr) in enumerate(mono_Ps)
                    Lr = size(Pr, 2)
                    psd_r = JuMP.@variable(model, [1:Lr, 1:Lr], Symmetric)
                    push!(psd_vars, psd_r)
                    JuMP.@constraint(model, psd_r in JuMP.PSDCone())
                end
                coeffs_lhs = JuMP.@expression(model, [u in 1:mono_U],
                    sum(sum(Pr[u, k] * Pr[u, l] * psd_r[k, l] * (k == l ? 1 : 2)
                    for k in 1:size(Pr, 2) for l in 1:k)
                    for (Pr, psd_r) in zip(mono_Ps, psd_vars)))
                JuMP.@constraint(model, coeffs_lhs .== gradient_interp)
            end
        end
    end

    # convexity
    if inst.use_convexity && !iszero(conv_profile)
        hessian_halfdeg = div(deg - 1, 2)
        (conv_U, conv_points, conv_Ps) =
            PolyUtils.interpolate(conv_dom, hessian_halfdeg)
        univ_chebs_der = [PolyUtils.calc_univariate_chebyshev(
            conv_points[:, i], 2halfdeg, calc_gradient = true,
            calc_hessian = true) for i in 1:n]

        deriv_num(i, j, k) = (k != i && k != j && return 1;
            k == i && k == j && return 3; return 2)

        V_Hs = Matrix{T}[]
        for i in 1:n, j in 1:i
            univ_chebs_H = [univ_chebs_der[k][deriv_num(i, j, k)] for k in 1:n]
            V_H = PolyUtils.make_product_vandermonde(univ_chebs_H,
                PolyUtils.n_deg_exponents(n, 2halfdeg))
            push!(V_Hs, V_H)
        end
        scal = inv(maximum(maximum(abs, V_H) for V_H in V_Hs))
        if scal < 1e-7
            @warn("model is numerically challenging to set up", maxlog = 1)
        end
        lmul!.(scal, V_Hs)

        H_interp = conv_profile * vcat([(F \ V_H')' *
            regressor for V_H in V_Hs]...)

        if inst.use_wsos
            if n == 1
                conv_K = Hypatia.WSOSInterpNonnegativeCone{T, T}(conv_U, conv_Ps)
            else
                Cones.scale_svec!(H_interp, sqrt(T(2)), incr = conv_U)
                conv_K = Hypatia.WSOSInterpPosSemidefTriCone{T}(n, conv_U, conv_Ps)
            end
            JuMP.@constraint(model, H_interp in conv_K)
        else
            psd_vars = []
            for (r, Pr) in enumerate(conv_Ps)
                Lr = size(Pr, 2)
                psd_r = JuMP.@variable(model, [1:(Lr * n), 1:(Lr * n)], Symmetric)
                push!(psd_vars, psd_r)
                JuMP.@constraint(model, psd_r in JuMP.PSDCone())
            end
            # for readability
            Ls = [size(Pr, 2) for Pr in conv_Ps]
            offset = 0
            for x1 in 1:n, x2 in 1:x1
                offset += 1
                coeffs_lhs = JuMP.@expression(model, [u in 1:conv_U],
                    sum(sum(conv_Ps[r][u, k] * conv_Ps[r][u, l] *
                    psd_vars[r][(x1 - 1) * Ls[r] + k, (x2 - 1) * Ls[r] + l]
                    for k in 1:Ls[r] for l in 1:Ls[r]) for r in eachindex(Ls)))
                JuMP.@constraint(model, coeffs_lhs .==
                    H_interp[conv_U .* (offset - 1) .+ (1:conv_U)])
            end
        end
    end

    return model
end

function test_extra(inst::ShapeConRegrJuMP{T}, model::JuMP.Model) where T
    stat = JuMP.termination_status(model)
    @test stat == MOI.OPTIMAL
    if (stat == MOI.OPTIMAL) && inst.is_fit_exact
        # check objective value is correct
        tol = eps(T)^0.25
        @test JuMP.objective_value(model) ≈ 0 atol=tol rtol=tol
    end
end

shapeconregr_data = Dict(
    :func1 => (x -> sum(x .^ 2)),
    :func2 => (x -> sum(x .^ 3)),
    :func3 => (x -> sum(x .^ 4)),
    :func4 => (x -> exp(norm(x)^2 / length(x)) - 1),
    :func5 => (x -> -inv(1 + exp(-10 * norm(x)))),
    :func6 => (x -> sum((x .+ 1) .^ 4)),
    :func7 => (x -> sum((x / 2 .+ 1) .^ 3)),
    :func8 => (x -> sum((x .+ 1) .^ 5 .- 2)),
    :func9 => (x -> (5x[1] + x[2] / 2 + x[3])^2 + sqrt(x[4]^2 + x[5]^2)),
    :func10 => (x -> sum(exp.(x))),
    )


text:
struct Shortcut
    s
end
Flux.@functor Shortcut
(s::Shortcut)(mx, x) = mx + s.s(x)

struct ResidualBlock
    block
end
Flux.@functor ResidualBlock
(b::ResidualBlock)(x) = x |> b.block .|> relu

function BasicBlock(
    channels::Pair{Int64, Int64}, connection; stride::Int64 = 1,
)
    layer = Chain(
        Conv((3, 3), channels; stride, pad=1, bias=false),
        BatchNorm(channels[2], relu),
        Conv((3, 3), channels[2]=>channels[2]; pad=1, bias=false),
        BatchNorm(channels[2]),
    )
    ResidualBlock(SkipConnection(layer, connection))
end

function Bottleneck(
    channels::Pair{Int64, Int64}, connection;
    stride::Int = 1, expansion::Int = 4,
)
    layer = Chain(
        Conv((1, 1), channels, bias=false),
        BatchNorm(channels[2], relu),
        Conv((3, 3), channels[2]=>channels[2]; stride, pad=1, bias=false),
        BatchNorm(channels[2], relu),
        Conv((1, 1), channels[2]=>(channels[2] * expansion); bias=false),
        BatchNorm(channels[2] * expansion),
    )
    ResidualBlock(SkipConnection(layer, connection))
end

function make_layer(
    block, channels::Pair{Int64, Int64}, repeat::Int64,
    expansion::Int64, stride::Int64 = 1,
)
    expanded_channels = channels[2] * expansion
    if stride == 1 && channels[1] == channels[2]
        connection = +
    else
        connection = Shortcut(Chain(
            Conv((1, 1), channels; stride, bias=false),
            BatchNorm(channels[2]),
        ))
    end

    layer = ResidualBlock[]
    push!(layer, block(channels, connection; stride))
    for i in 2:repeat
        push!(layer, block(expanded_channels=>channels[2], +))
    end
    Chain(layer...)
end


text:
using Oceananigans.Grids: domain_string, short_show
using Oceananigans.Fields: show_location

for op_string in ("UnaryOperation", "BinaryOperation", "MultiaryOperation", "Derivative")
    op = eval(Symbol(op_string))
    @eval begin
        operation_name(::$op) = $op_string
    end
end

Base.show(io::IO, operation::AbstractOperation) =
    print(io,
          operation_name(operation), " at ", show_location(operation), '\n',
          "├── grid: ", typeof(operation.grid), '\n',
          "│   ├── size: ", size(operation.grid), '\n',
          "│   └── domain: ", domain_string(operation.grid), '\n',
          "└── tree: ", "\n", "    ", tree_show(operation, 1, 0))

"Return a representaion of number or function leaf within a tree visualization of an `AbstractOperation`."
tree_show(a::Union{Number, Function}, depth, nesting) = string(a)

"Fallback for displaying a leaf within a tree visualization of an `AbstractOperation`."
tree_show(a, depth, nesting) = short_show(a) # fallback

"Returns a string corresponding to padding characters for a tree visualization of an `AbstractOperation`."
get_tree_padding(depth, nesting) = "    "^(depth-nesting) * "│   "^nesting

"Return a string representaion of a `UnaryOperation` leaf within a tree visualization of an `AbstractOperation`."
function tree_show(unary::UnaryOperation{X, Y, Z}, depth, nesting)  where {X, Y, Z}
    padding = get_tree_padding(depth, nesting)

    return string(unary.op, " at ", show_location(X, Y, Z), " via ", unary.▶, '\n',
                  padding, "└── ", tree_show(unary.arg, depth+1, nesting))
end

"Return a string representaion of a `BinaryOperation` leaf within a tree visualization of an `AbstractOperation`."
function tree_show(binary::BinaryOperation{X, Y, Z}, depth, nesting) where {X, Y, Z}
    padding = get_tree_padding(depth, nesting)

    return string(binary.op, " at ", show_location(X, Y, Z), " via ", binary.▶op, '\n',
                  padding, "├── ", tree_show(binary.a, depth+1, nesting+1), '\n',
                  padding, "└── ", tree_show(binary.b, depth+1, nesting))
end

"Return a string representaion of a `MultiaryOperation` leaf within a tree visualization of an `AbstractOperation`."
function tree_show(multiary::MultiaryOperation{X, Y, Z, N}, depth, nesting) where {X, Y, Z, N}
    padding = get_tree_padding(depth, nesting)

    out = string(multiary.op, " at ", show_location(X, Y, Z), '\n',
        ntuple(i -> padding * "├── " * tree_show(multiary.args[i], depth+1, nesting+1) * '\n', Val(N-1))...,
                    padding * "└── " * tree_show(multiary.args[N], depth+1, nesting)
                )
    return out
end

"Return a string representaion of a `Derivative` leaf within a tree visualization of an `AbstractOperation`."
function tree_show(deriv::Derivative{X, Y, Z}, depth, nesting)  where {X, Y, Z}
    padding = get_tree_padding(depth, nesting)

    return string(deriv.∂, " at ", show_location(X, Y, Z), " via ", deriv.▶, '\n',
                  padding, "└── ", tree_show(deriv.arg, depth+1, nesting))
end


text:
using LightXML
using Iterators
import HDF5

function parameters2xml(params::Dict, prefix::ASCIIString)
  for (i, param_set) in enumerate(product(values(params)...))
    xdoc = XMLDocument()
    xroot = create_root(xdoc, "SIMULATION")
    pn = new_child(xroot, "PARAMETERS")
    for (k, p) in zip(keys(params), param_set)
      pc = new_child(pn, "PARAMETER")
      add_text(pc, string(p))
      set_attribute(pc, "name", string(k))
    end
    save_file(xdoc, prefix * ".task" * string(i) * ".in.xml")
  end
end

function xml2parameters(prefix::ASCIIString, idx::Int)
  println("Prefix is ", prefix, " and idx is ", idx)
  params = Dict{Any, Any}()
  xdoc = parse_file(prefix * ".task" * string(idx) * ".in.xml")
  xroot = root(xdoc)
  for c in child_nodes(xroot)
    if is_elementnode(c)
      for p in child_nodes(c)
        if is_elementnode(p)
          e = XMLElement(p)
          if name(e) == "PARAMETER"
            println(attribute(e, "name"), " = ", content(e))
            params[attribute(e, "name")] = content(e)
          end
        end
      end
    end
  end
  return params
end

function parameters2hdf5(params::Dict, filename::ASCIIString)
  f = HDF5.h5open(filename, "r+")

  for (k, v) in params
    try
      if HDF5.exists(f, "parameters/" * k)
        if read(f["parameters/"*k]) != v
          error(k, " exists but differs from current ")
        end
      else
        f["parameters/" * k] = v
      end
    catch e
    end

  end

  close(f)
end

max_stars_repo_name: pebroecker/CythonHardcoreBosons

max_issues_repo_name: pebroecker/CythonHardcoreBosons

max_forks_repo_name: pebroecker/CythonHardcoreBosons


text:
export CentralFlux
export RusanovFlux
export EntropyConservativeFlux
export KennedyGruberFlux
export RoeFlux
export MatrixFlux

avg(s⁻, s⁺) = (s⁻ + s⁺) / 2
function logavg(a, b)
    ζ = a / b
    f = (ζ - 1) / (ζ + 1)
    u = f^2
    ϵ = eps(eltype(u))

    if u < ϵ
        F = @evalpoly(u, one(u), one(u) / 3, one(u) / 5, one(u) / 7, one(u) / 9)
    else
        F = log(ζ) / 2f
    end

    (a + b) / 2F
end
roe_avg(ρ⁻, ρ⁺, s⁻, s⁺) = (sqrt(ρ⁻) * s⁻ + sqrt(ρ⁺) * s⁺) / (sqrt(ρ⁻) + sqrt(ρ⁺))

abstract type AbstractNumericalFlux end
function twopointflux end
function surfaceflux(flux::AbstractNumericalFlux, law::AbstractBalanceLaw, n⃗, q⁻, aux⁻, q⁺, aux⁺)
  twopointflux(flux, law, q⁻, aux⁻, q⁺, aux⁺)' * n⃗
end

struct CentralFlux <: AbstractNumericalFlux end
function twopointflux(::CentralFlux, law::AbstractBalanceLaw, q₁, aux₁, q₂, aux₂)
  (flux(law, q₁, aux₁) + flux(law, q₂, aux₂)) / 2
end

struct RusanovFlux <: AbstractNumericalFlux end
function surfaceflux(::RusanovFlux, law::AbstractBalanceLaw, n⃗, q⁻, aux⁻, q⁺, aux⁺)
  fc = surfaceflux(CentralFlux(), law, n⃗, q⁻, aux⁻, q⁺, aux⁺)
  ws⁻ = wavespeed(law, n⃗, q⁻, aux⁻)
  ws⁺ = wavespeed(law, n⃗, q⁺, aux⁺)
  fc - max(ws⁻, ws⁺) * (q⁺ - q⁻) / 2
end

struct EntropyConservativeFlux <: AbstractNumericalFlux end
struct KennedyGruberFlux <: AbstractNumericalFlux end
struct RoeFlux <: AbstractNumericalFlux end
struct MatrixFlux <: AbstractNumericalFlux end


text:
function reptag!(atag::tag, rtag::tag, repcount)
  tag(atag, fill(rtag, repcount));
  end;

function reptag!(taglist::Array{tag}, replist::Array{Int})
  n = size(replist,1);
  n2 = size(taglist,1)
  if(n!=n2 || replist[1]!= 1)
    println("$(n) != $(n2) || replist[1] != 1");
    return [];
    end;

  for iter = n2:-1:2
    reptag!(taglist[iter-1], taglist[iter], replist[iter]);
    ##tag(taglist[iter-1], fill(taglist[iter], replist[iter]));
    end;
  end;

function clonetag(atag::tag)
  return tag(atag.name, atag.pairs);
  end;

##function cyclic_check(atag::tag, xtaglist::Array{tag})
##  return any(x->x===atag, xtaglist)
##  end;

function cc_nodes(xnodelist, xnodeptr)
  znode = [];
  i = xnodeptr;
  while(i>0)
    push!(znode, i);
    i = xnodelist[i];
    end;
  return znode;
  end;

function cc_mode(xptr, xtag, xlist, ccstack, xnodes)
  if( !cyclic_check(xtag, ccstack))
    push!(ccstack, xtag);
    return false;
    end;
  xid = findfirst(x->x===xtag, ccstack);
  return cyclic_check(xtag, xlist[cc_nodes(xnodes, xptr)]);
  end;

function dereptags(atag::tag)

  taglist = tag[atag];
  clonelist = tag[];
  nodelist = [0];
  cctagstack = [atag];
  tagptr = 1;
  itercount = 1;

  while tagptr <= itercount ##tp_ptr>0
    local ctag = nodelist[tagptr]; ##tagptr[tp_ptr];
    local xtag;

    if(ctag == 0)
      xtag = taglist[tagptr];
      push!(clonelist, clonetag(xtag));
      m = size(xtag.attach,1);
      append!(taglist, xtag.attach);
      append!(nodelist, fill(tagptr, m));
      itercount += m;
      tagptr += 1;
      tagptr > itercount && return clonelist[1];
      ctag = nodelist[tagptr];
      end;

    if(!cc_mode(ctag, taglist[tagptr], taglist, cctagstack, nodelist))
      xtag = clonetag(taglist[tagptr]); ##.attach[ctag];
      tag(clonelist[nodelist[tagptr]], xtag);
      push!(clonelist, xtag);
      m = size(taglist[tagptr].attach,1);
      append!(taglist, taglist[tagptr].attach);
      append!(nodelist, fill(tagptr, m));
      itercount += m;
      end;
    tagptr += 1;
    end;

  return clonelist[1];
  end;


text:
using Pkg
Pkg.activate(".")
Pkg.instantiate()

using Weave

const content_folder = joinpath(pwd(), "content")

scss_files = ["configuration", "index"]
for scss_file in scss_files
    run(`sass themes/shebang/static/css/$(scss_file).scss themes/shebang/static/css/$(scss_file).css`)
end

for content_type in ["lessons", "primers", "capstones"]
    this_content_folder = joinpath(content_folder, content_type)
    raw_files = filter(x -> endswith(x, ".Jmd"), readdir(this_content_folder))
    for this_file in raw_files
        target_file = replace(this_file, ".Jmd" => ".md")
        file_time = mtime(joinpath(this_content_folder, this_file))
        trgt_time = mtime(joinpath(this_content_folder, target_file))

        if file_time > trgt_time
            @info "Updating:   $(joinpath(content_type, target_file))"
            weave(joinpath(this_content_folder, this_file), doctype="hugo")
        else
            @info "Up-to-date: $(joinpath(content_type, target_file))"
        end
    end
end


text:
# test Base Extensions
@testset "Base Extensions" begin
    # Setup data
    m = InfiniteModel();
    gvref = GeneralVariableRef(m, 1, PointVariableIndex)
    gvref2 = GeneralVariableRef(m, 2, PointVariableIndex)
    ivref = InfiniteVariableRef(m, InfiniteVariableIndex(1))
    hvref = HoldVariableRef(m, HoldVariableIndex(1))
    # test Base.copy
    @testset "Base.copy" begin
        @test copy(gvref) == gvref
        @test copy(ivref) == ivref
    end
    # test Base.:(==)
    @testset "Base.:(==)" begin
        @test gvref == gvref
        @test ivref == ivref
        @test !(gvref == ivref)
        @test !(hvref == ivref)
        @test !(gvref == gvref2)
    end
    # test Base.broadcastable
    @testset "Base.broadcastable" begin
        @test Base.broadcastable(gvref) isa Base.RefValue{GeneralVariableRef}
        @test Base.broadcastable(ivref) isa Base.RefValue{InfiniteVariableRef}
    end
    # test Base.length
    @testset "Base.length" begin
        @test length(gvref) == 1
        @test length(ivref) == 1
    end
    # test JuMP.isequal_canonical
    @testset "JuMP.isequal_canonical" begin
        @test JuMP.isequal_canonical(gvref, gvref)
        @test JuMP.isequal_canonical(ivref, ivref)
    end
    # test JuMP.variable_type
    @testset "JuMP.variable_type" begin
        @test JuMP.variable_type(m) == GeneralVariableRef
    end
    # test _remove_name_index
    @testset "_remove_name_index" begin
        pref = @independent_parameter(m, test, set = IntervalSet(0, 1))
        pref2 = @independent_parameter(m, θ, set = IntervalSet(0, 1))
        prefs = @independent_parameter(m, x[1:2], set = IntervalSet(0, 1))
        @test InfiniteOpt._remove_name_index(pref) == "test"
        @test InfiniteOpt._remove_name_index(prefs[1]) == "x"
        @test InfiniteOpt._remove_name_index(pref2) == "θ"
    end
    # test Basic getters
    @testset "Attribute Accessors" begin
        @test InfiniteOpt._index_type(gvref) == PointVariableIndex
        @test InfiniteOpt._raw_index(gvref) == 1
        @test InfiniteOpt._param_index(gvref) == -1
    end
end

# test Reference Accessers
@testset "Reference Accessers" begin
    # Setup data
    m = InfiniteModel();
    gvref1 = GeneralVariableRef(m, 1, PointVariableIndex)
    gvref2 = GeneralVariableRef(m, 1, DependentParameterIndex, 2)
    ivref = InfiniteVariableRef(m, InfiniteVariableIndex(1))
    # test JuMP.index (GeneralVariableRef)
    @testset "JuMP.index (GeneralVariableRef)" begin
        @test index(gvref1) == PointVariableIndex(1)
        @test index(gvref2) == DependentParameterIndex(DependentParametersIndex(1), 2)
    end
    # test JuMP.index (DispatchVariableRef)
    @testset "JuMP.index (DispatchVariableRef)" begin
        @test index(ivref) == InfiniteVariableIndex(1)
    end
    # test JuMP.owner_model (GeneralVariableRef)
    @testset "JuMP.owner_model (GeneralVariableRef)" begin
        @test owner_model(gvref1) == m
        @test owner_model(gvref2) == m
    end
    # test JuMP.owner_model (DispatchVariableRef)
    @testset "JuMP.owner_model (DispatchVariableRef)" begin
        @test owner_model(ivref) == m
    end
end

# test Dispatch Variabel Methods
@testset "dispatch_variable_ref" begin
    # Setup data
    m = InfiniteModel();
    idx = TestIndex(1)
    dvref = TestVariableRef(m, idx)
    gvref = GeneralVariableRef(m, 1, TestIndex)
    # test the test type
    @testset "TestVariableRef" begin
        @test dispatch_variable_ref(m, idx) == dvref
    end
    # test GeneralVariableRef input
    @testset "GeneralVariableRef" begin
        @test dispatch_variable_ref(gvref) == dvref
    end
end

# test Naming Methods
@testset "Naming Methods" begin
    # Setup data
    m = InfiniteModel();
    idx = TestIndex(1)
    dvref = TestVariableRef(m, idx)
    gvref = GeneralVariableRef(m, 1, TestIndex)
    # test JuMP.name (Fallback)
    @testset "JuMP.name (Fallback)" begin
        @test_throws ArgumentError name(dvref)
    end
    # test JuMP.name (GeneralVariableRef)
    @testset "JuMP.name (GeneralVariableRef)" begin
        @test_throws ArgumentError name(gvref)
    end
    # test JuMP.set_name (Fallback)
    @testset "JuMP.set_name (Fallback)" begin
        @test_throws ArgumentError set_name(dvref, "new")
    end
    # test JuMP.set_name (GeneralVariableRef)
    @testset "JuMP.set_name (GeneralVariableRef)" begin
        @test_throws ArgumentError set_name(gvref, "new")
    end
end

# test Core data methods
@testset "Core Data Methods" begin
    # Setup data
    m = InfiniteModel();
    gvref = GeneralVariableRef(m, 1, TestIndex)
    idx = TestIndex(1)
    dvref = TestVariableRef(m, idx)
    # test _data_dictionary with GeneralVariableRef
    @testset "_data_dictionary (GeneralVariableRef)" begin
        @test_throws MethodError InfiniteOpt._data_dictionary(gvref)
    end
    # test _data_object with GeneralVariableRef
    @testset "_data_object (GeneralVariableRef)" begin
        @test_throws MethodError InfiniteOpt._data_object(gvref)
    end
    # test _delete_data_object with DispatchVariableRef
    @testset "_delete_data_object (DispatchVariableRef)" begin
        @test_throws MethodError InfiniteOpt._delete_data_object(dvref)
    end
end

# test Validity Methods
@testset "Validity Methods" begin
    # Setup data
    m = InfiniteModel();
    idx = TestIndex(1)
    dvref = TestVariableRef(m, idx)
    gvref = GeneralVariableRef(m, 1, TestIndex)
    # test JuMP.is_valid (Fallback)
    @testset "JuMP.is_valid (Fallback)" begin
        @test_throws MethodError is_valid(m, dvref)
    end
    # test JuMP.is_valid (GeneralVariableRef)
    @testset "JuMP.is_valid (GeneralVariableRef)" begin
        @test_throws MethodError is_valid(m, gvref)
    end
end

# test Core Object Methods
@testset "Core Object Methods" begin
    # Setup data
    m = InfiniteModel();
    gvref = GeneralVariableRef(m, 1, TestIndex)
    # test _core_variable_object (GeneralVariableRef)
    @testset "_core_variable_object (GeneralVariableRef)" begin
        @test_throws MethodError InfiniteOpt._core_variable_object(gvref)
    end
end

# test Dependency Methods
@testset "Dependency Methods" begin
    # Setup data
    m = InfiniteModel();
    gvref = GeneralVariableRef(m, 1, TestIndex)
    # Loop over the various functions
    for f in (:_infinite_variable_dependencies, :_reduced_variable_dependencies,
              :_point_variable_dependencies, :_measure_dependencies,
              :_constraint_dependencies, :_derivative_dependencies,
              :_derivative_constraint_dependencies,
              :_parameter_function_dependencies, :_generative_measures)
        @test_throws MethodError InfiniteOpt.eval(f)(gvref)
    end
end

# test Usage Methods
@testset "Usage Methods" begin
    # Setup data
    m = InfiniteModel();
    idx = TestIndex(1)
    dvref = TestVariableRef(m, idx)
    gvref = GeneralVariableRef(m, 1, TestIndex)
    for f in (:used_by_infinite_variable, :used_by_reduced_variable,
              :used_by_point_variable, :used_by_measure, :used_by_constraint,
              :used_by_objective, :used_by_derivative, :is_used,
              :has_derivative_constraints, :used_by_parameter_function)
        @test_throws ArgumentError eval(f)(dvref)
        @test_throws ArgumentError eval(f)(gvref)
    end
end

# test Delete Methods
@testset "Delete Methods" begin
    # Setup data
    m = InfiniteModel();
    idx = TestIndex(1)
    dvref = TestVariableRef(m, idx)
    gvref = GeneralVariableRef(m, 1, TestIndex)
    # test JuMP.delete (Fallback)
    @testset "JuMP.delete (Fallback)" begin
        @test_throws ArgumentError JuMP.delete(m, dvref)
    end
    # test JuMP.delete (GeneralVariableRef)
    @testset "JuMP.delete (GeneralVariableRef)" begin
        @test_throws ArgumentError JuMP.delete(m, gvref)
    end
    # test JuMP.delete (GeneralVariableRefs)
    @testset "JuMP.delete (GeneralVariableRefs)" begin
        @test_throws ArgumentError JuMP.delete(m, [gvref])
    end
end

# test Parameter Methods
@testset "Parameter Methods" begin
    # Setup data
    m = InfiniteModel();
    idx = TestIndex(1)
    dvref = TestVariableRef(m, idx)
    gvref = GeneralVariableRef(m, 1, TestIndex)
    gvrefs = [GeneralVariableRef(m, 1, TestIndex),
              GeneralVariableRef(m, 2, TestIndex)]
    sparse_gvrefs = convert(JuMPC.SparseAxisArray, gvrefs)
    dense_gvrefs = JuMPC.DenseAxisArray(gvrefs, 3:4)
    sset = IntervalSet(0, 1)
    aset = CollectionSet([sset, sset])
    # test _parameter_number
    @testset "_parameter_number" begin
        @test_throws MethodError InfiniteOpt._parameter_number(gvref)
    end
    # test _object_number
    @testset "_object_number" begin
        @test_throws MethodError InfiniteOpt._object_number(gvref)
    end
    # test 1 argument methods and fallbacks
    for f in (:significant_digits, :parameter_value, :derivative_method,
              :has_generative_supports, :has_internal_supports,
              :add_generative_supports, :raw_function, :generative_support_info)
        @test_throws ArgumentError eval(f)(dvref)
        @test_throws ArgumentError eval(f)(gvref)
    end
    # test infinite_set (Fallback)
    @testset "infinite_set (Fallback)" begin
        @test_throws ArgumentError infinite_set(dvref)
    end
    # test infinite_set (GeneralVariableRef)
    @testset "infinite_set (GeneralVariableRef)" begin
        @test_throws ArgumentError infinite_set(gvref)
    end
    # test infinite_set (GeneralVariableRef Array)
    @testset "infinite_set (GeneralVariableRef Array)" begin
        @test_throws ArgumentError infinite_set(gvrefs)
        @test_throws ArgumentError infinite_set(dense_gvrefs)
        @test_throws ArgumentError infinite_set(sparse_gvrefs)
    end
    # test set_infinite_set (Fallback)
    @testset "set_infinite_set (Fallback)" begin
        @test_throws ArgumentError set_infinite_set(dvref, sset)
    end
    # test set_infinite_set (GeneralVariableRef)
    @testset "set_infinite_set (GeneralVariableRef)" begin
        @test_throws ArgumentError set_infinite_set(gvref, sset)
    end
    # test set_infinite_set (GeneralVariableRef Array)
    @testset "set_infinite_set (GeneralVariableRef Array)" begin
        @test_throws ArgumentError set_infinite_set(gvrefs, aset)
        @test_throws ArgumentError set_infinite_set(dense_gvrefs, aset)
        @test_throws ArgumentError set_infinite_set(sparse_gvrefs, aset)
    end
    # test num_supports (Fallback)
    @testset "num_supports (Fallback)" begin
        @test_throws ArgumentError num_supports(dvref, label = All)
    end
    # test num_supports (GeneralVariableRef)
    @testset "num_supports (GeneralVariableRef)" begin
        @test_throws ArgumentError num_supports(gvref, label = All)
    end
    # test num_supports (GeneralVariableRef Array)
    @testset "num_supports (GeneralVariableRef Array)" begin
        @test_throws ArgumentError num_supports(gvrefs, label = All)
        @test_throws ArgumentError num_supports(dense_gvrefs, label = All)
        @test_throws ArgumentError num_supports(sparse_gvrefs, label = All)
    end
    # test has_supports (Fallback)
    @testset "has_supports (Fallback)" begin
        @test_throws ArgumentError has_supports(dvref)
    end
    # test has_supports (GeneralVariableRef)
    @testset "has_supports (GeneralVariableRef)" begin
        @test_throws ArgumentError has_supports(gvref)
    end
    # test has_supports (GeneralVariableRef Array)
    @testset "has_supports (GeneralVariableRef Array)" begin
        @test_throws ArgumentError has_supports(gvrefs)
        @test_throws ArgumentError has_supports(dense_gvrefs)
        @test_throws ArgumentError has_supports(sparse_gvrefs)
    end
    # test supports (Fallback)
    @testset "supports (Fallback)" begin
        @test_throws ArgumentError supports(dvref, label = All)
    end
    # test supports (GeneralVariableRef)
    @testset "supports (GeneralVariableRef)" begin
        @test_throws ArgumentError supports(gvref, label = All)
    end
    # test supports (GeneralVariableRef Array)
    @testset "supports (GeneralVariableRef Array)" begin
        @test_throws ArgumentError supports(gvrefs, label = All)
        @test_throws ArgumentError supports(dense_gvrefs, label = All)
        @test_throws ArgumentError supports(sparse_gvrefs, label = All)
    end
    # test set_supports (Fallback)
    @testset "set_supports (Fallback)" begin
        @test_throws ArgumentError set_supports(dvref, 2, label = All, force = true)
    end
    # test set_supports (GeneralVariableRef)
    @testset "set_supports (GeneralVariableRef)" begin
        @test_throws ArgumentError set_supports(gvref, 2, label = All, force = true)
    end
    # test set_supports (GeneralVariableRef Array)
    @testset "set_supports (GeneralVariableRef Array)" begin
        @test_throws ArgumentError set_supports(gvrefs, ones(1, 1), label = All,
                                                force = true)
        @test_throws ArgumentError set_supports(dense_gvrefs, ones(1, 1),
                                                label = All, force = true)
        @test_throws ArgumentError set_supports(sparse_gvrefs, ones(1, 1),
                                                label = All, force = true)
    end
    # test add_supports (Fallback)
    @testset "add_supports (Fallback)" begin
        @test_throws ArgumentError add_supports(dvref, 2, label = All, check = true)
    end
    # test add_supports (GeneralVariableRef)
    @testset "add_supports (GeneralVariableRef)" begin
        @test_throws ArgumentError add_supports(gvref, 2, label = All, check = true)
    end
    # test add_supports (GeneralVariableRef Array)
    @testset "add_supports (GeneralVariableRef Array)" begin
        @test_throws ArgumentError add_supports(gvrefs, ones(1, 1), label = All,
                                                check = true)
        @test_throws ArgumentError add_supports(dense_gvrefs, ones(1, 1),
                                                label = All, check = true)
        @test_throws ArgumentError add_supports(sparse_gvrefs, ones(1, 1),
                                                label = All, check = true)
    end
    # test delete_supports (Fallback)
    @testset "delete_supports (Fallback)" begin
        @test_throws ArgumentError delete_supports(dvref)
    end
    # test delete_supports (GeneralVariableRef)
    @testset "delete_supports (GeneralVariableRef)" begin
        @test_throws ArgumentError delete_supports(gvref)
    end
    # test delete_supports (GeneralVariableRef Array)
    @testset "delete_supports (GeneralVariableRef Array)" begin
        @test_throws ArgumentError delete_supports(gvrefs)
        @test_throws ArgumentError delete_supports(dense_gvrefs)
        @test_throws ArgumentError delete_supports(sparse_gvrefs)
    end
    # test fill_in_supports! (Fallback)
    @testset "fill_in_supports! (Fallback)" begin
        @test_throws ArgumentError fill_in_supports!(dvref, num_supports = 1)
    end
    # test fill_in_supports! (GeneralVariableRef)
    @testset "fill_in_supports! (GeneralVariableRef)" begin
        @test_throws ArgumentError fill_in_supports!(gvref, num_supports = 1)
    end
    # test fill_in_supports! (GeneralVariableRef Array)
    @testset "fill_in_supports! (GeneralVariableRef Array)" begin
        @test_throws ArgumentError fill_in_supports!(gvrefs, num_supports = 1)
        @test_throws ArgumentError fill_in_supports!(dense_gvrefs, num_supports = 1)
        @test_throws ArgumentError fill_in_supports!(sparse_gvrefs,num_supports = 1)
    end
    # test JuMP.set_value (Fallback)
    @testset "JuMP.set_value (Fallback)" begin
        @test_throws ArgumentError set_value(dvref, 2)
    end
    # test JuMP.set_value (GeneralVariableRef)
    @testset "JuMP.set_value (GeneralVariableRef)" begin
        @test_throws ArgumentError set_value(gvref, 2)
    end
    # test set_derivative_method (Fallback)
    @testset "set_derivative_method (Fallback)" begin
        @test_throws ArgumentError set_derivative_method(dvref, FiniteDifference())
    end
    # test set_derivative_method (GeneralVariableRef)
    @testset "set_derivative_method (GeneralVariableRef)" begin
        @test_throws ArgumentError set_derivative_method(gvref, FiniteDifference())
    end
    # test parameter status setters
    for f in (:_set_has_generative_supports, :_set_has_internal_supports,
              :_set_has_derivative_constraints)
        @test_throws ArgumentError InfiniteOpt.eval(f)(dvref, true)
        @test_throws ArgumentError InfiniteOpt.eval(f)(gvref, true)
    end
end

# test Variable Methods
@testset "Variable Methods" begin
    # Setup data
    m = InfiniteModel();
    idx = TestIndex(1)
    dvref = TestVariableRef(m, idx)
    gvref = GeneralVariableRef(m, 1, TestIndex)
    # test 1 argument methods
    for f in (:raw_parameter_refs, :parameter_refs, :parameter_list,
              :start_value_function, :reset_start_value_function,
              :infinite_variable_ref, :eval_supports, :raw_parameter_values,
              :parameter_values, :parameter_bounds, :delete_parameter_bounds)
        @test_throws ArgumentError eval(f)(dvref)
        @test_throws ArgumentError eval(f)(gvref)
    end
    # test has_parameter_vbounds (GeneralVariableRef)
    @testset "has_parameter_bounds (GeneralVariableRef)" begin
        @test !has_parameter_bounds(gvref)
    end
    # test set_parameter_bounds (Fallback)
    @testset "set_parameter_bounds (Fallback)" begin
        @test_throws ArgumentError set_parameter_bounds(dvref, ParameterBounds())
    end
    # test set_parameter_vbounds (GeneralVariableRef)
    @testset "set_parameter_bounds (GeneralVariableRef)" begin
        @test_throws ArgumentError set_parameter_bounds(gvref, ParameterBounds(),
                                                        force = true)
    end
    # test add_parameter_bounds (Fallback)
    @testset "add_parameter_bounds (Fallback)" begin
        @test_throws ArgumentError add_parameter_bounds(dvref, ParameterBounds())
    end
    # test add_parameter_vbounds (GeneralVariableRef)
    @testset "add_parameter_bounds (GeneralVariableRef)" begin
        @test_throws ArgumentError add_parameter_bounds(gvref, ParameterBounds())
    end
end

# test Measure Methods
@testset "Measure Methods" begin
    # Setup data
    m = InfiniteModel();
    idx = TestIndex(1)
    dvref = TestVariableRef(m, idx)
    gvref = GeneralVariableRef(m, 1, TestIndex)
    # loop of the methods
    for f in (:measure_function, :measure_data, :is_analytic, :expand)
        @test_throws ArgumentError eval(f)(dvref)
        @test_throws ArgumentError eval(f)(gvref)
    end
end

# test Derivative Methods
@testset "Derivative Methods" begin
    # Setup data
    m = InfiniteModel();
    idx = TestIndex(1)
    dvref = TestVariableRef(m, idx)
    gvref = GeneralVariableRef(m, 1, TestIndex)
    # loop through methods
    for f in (:derivative_argument, :operator_parameter, :evaluate,
              :derivative_constraints, :delete_derivative_constraints)
        @test_throws ArgumentError eval(f)(dvref)
        @test_throws ArgumentError eval(f)(gvref)
    end
end

# test variable info methods
@testset "Variable Info Methods" begin
    # Setup data
    m = InfiniteModel();
    idx = TestIndex(1)
    dvref = TestVariableRef(m, idx)
    gvref = GeneralVariableRef(m, 1, TestIndex)
    # test 1 argument methods
    for f in (:has_lower_bound, :has_upper_bound, :is_fixed, :is_binary,
              :is_integer, :lower_bound, :upper_bound, :fix_value, :start_value,
              :set_binary, :set_integer, :LowerBoundRef, :UpperBoundRef, :FixRef,
              :BinaryRef, :IntegerRef, :delete_lower_bound, :delete_upper_bound,
              :unfix, :unset_binary, :unset_integer)
        @test_throws ArgumentError eval(f)(dvref)
        @test_throws ArgumentError eval(f)(gvref)
    end
    # test setting methods
    for f in (:set_lower_bound, :set_upper_bound, :set_start_value)
        @test_throws ArgumentError eval(f)(dvref, 42)
        @test_throws ArgumentError eval(f)(gvref, 42)
    end
    # test JuMP.fix (Fallback)
    @testset "JuMP.fix (Fallback)" begin
        @test_throws ArgumentError fix(dvref, 42, force = true)
    end
    # test JuMP.fix (GeneralVariableRef)
    @testset "JuMP.fix (GeneralVariableRef)" begin
        @test_throws ArgumentError fix(gvref, 42, force = true)
    end
    # test set_start_value_function (Fallback)
    @testset "set_start_value_function (Fallback)" begin
        @test_throws ArgumentError set_start_value_function(dvref, sin)
    end
    # test set_start_value_function (GeneralVariableRef)
    @testset "set_start_value_function (GeneralVariableRef)" begin
        @test_throws ArgumentError set_start_value_function(gvref, sin)
    end
end


text:
########## Routines related to setting up interpolators.

##### Routines related to Barycentric interpolation.

# Barycentric Interpolations
function Barycentric𝓵(x::Real, 𝓧)
    return prod( x-𝓧[i] for i = 1:length(𝓧))
end

# This is l' in thesis (see Eq. 3.15).
function Barycentric𝓵prime(k::Int, 𝓧)
    return prod( 𝓧[k]-𝓧[i] for i = 1:length(𝓧) if i != k)
end

function get𝑤(k::Int, 𝓧::Vector{T})::T where T
    return one(T)/Barycentric𝓵prime(k, 𝓧)
end

function get𝑤(𝓧::Vector{T})::Vector{T} where T
    return collect( get𝑤(k,𝓧) for k = 1:length(𝓧) )
end

function get𝑤scaled!(w::Vector{T}, 𝓧::Vector{T})::Nothing where T
    L = length(𝓧) - 2
    factor::T = L*log(one(T)+one(T)) # C is 0.5 for the interval [-1,1]

    resize!(w,length(𝓧))
    for k = 0:length(𝓧)-1

        # convert to Int to speed up the multiplication used in prod().
        numerator = prod( convert(Int,sign(𝓧[k+1]-𝓧[i])) for i = 1:length(𝓧) if i != k+1)

        # log-sum-exp trick for adding ill-conditioned numbers.
        denominator = exp( factor + sum( log(abs(𝓧[k+1]-𝓧[i])) for i = 1:length(𝓧) if i != k+1 ) )

        w[k+1] = numerator/denominator
        if !isfinite(w[k+1])
            println("𝓧 is ", 𝓧)
            println("k is ", k)
            println("𝓧[k] is ", 𝓧[k])
            println("numerator is ", numerator)
            println("denominator is ", denominator)
        end
        @assert isfinite(w[k+1])
    end

    return nothing
end

function get𝑤scaled(𝓧::Vector{T})::Vector{T} where T

    w = Vector{T}(undef,length(𝓧))
    get𝑤scaled!(w,𝓧)

    return w
end

# uses SharedArray{Float64} for internal computation.
function get𝑤scaledfast(𝓧::Vector{T})::Vector{T} where T

    L = length(𝓧) - 2
    factor::T = L*log(one(T)+one(T)) # C is 0.5 for the interval [-1,1]

    w = SharedArray{Float64}(length(𝓧))
    @sync @distributed for k = 0:length(𝓧)-1

        numerator::T = prod( sign(𝓧[k+1]-𝓧[i]) for i = 1:length(𝓧) if i != k+1)
        denominator::T = exp( factor + sum( log(abs(𝓧[k+1]-𝓧[i])) for i = 1:length(𝓧) if i != k+1 ) )

        w[k+1] = convert(Float64,numerator/denominator)
    end

    # sanity check.
    for k = 0:length(𝓧)-1
        if !isfinite(w[k+1])
            println("𝓧 is ", 𝓧)
            println("k is ", k)
            println("𝓧[k] is ", 𝓧[k])
            println("numerator is ", numerator)
            println("denominator is ", denominator)
        end
        @assert isfinite(w[k+1])
    end

    return convert(Vector{T},w)
end

##### Routines related to Chebyshev interpolation at nodes.

function getChebyshevinterpolator(  f::Function,
                                    L::Int,
                                    νf::Vector{T},
                                    ν::Vector{T})::Tuple{Vector{T},Vector{T}} where T
    Lp1 = L + 1

    a = Vector{BigFloat}(undef,Lp1)
    for k = 0:L

        term1 = 1/L * f(νf[1]) * Chebyshev1st(ν[1],k)

        term2 = BigFloat("0")
        for l = 1:L-1
            term2 += f(νf[l+1]) * Chebyshev1st(ν[l+1],k)
        end
        term2 *= 2/L

        term3 = 1/L * f(νf[Lp1]) * Chebyshev1st(ν[Lp1],k)

        a[k+1] = term1 + term2 + term3
    end
    a[1] = a[1]/BigFloat(2)

    return a, ν
end

# efficient version. Less accurate.
function getChebyshevinterpolatorfast(  f::Function,
                                        L::Int,
                                        νf::SharedArray{Float64,1},
                                        ν::SharedArray{Float64,1})::Tuple{SharedArray{Float64,1},SharedArray{Float64,1}}
    Lp1 = L + 1

    ### pre-compute.

    # f(νf), length νf.
    f_νf = SharedArray{Float64}(length(νf))
    @sync @distributed for i = 1:length(f_νf)
        f_νf[i] = convert(Float64, f( convert(BigFloat,νf[i]) ))
    end

    ###

    a = SharedArray{Float64}(Lp1)
    for k = 0:L

        term1 = 1/L * f_νf[1] * Chebyshev1st(ν[1],k)

        term2 = zero(Float64)
        for l = 1:L-1
            term2 += f_νf[l+1] * Chebyshev1st(ν[l+1],k)
        end
        term2 *= 2/L

        term3 = 1/L * f_νf[Lp1] * Chebyshev1st(ν[Lp1],k)

        a[k+1] = term1 + term2 + term3
    end
    a[1] = a[1]/2

    return a, ν
end

# Barycentric interpolation coefficients to Chebyshev interpolation coefficients.
# Used in the minimax solution's conversion to impulse response in packagefiltersolution().
function BarycentrictoChebyshev(𝑋, wfunc, dfunc, type_num)
    𝑤 = get𝑤scaled(𝑋)
    h = geth(𝑤, wfunc, dfunc, 𝑋)
    𝑓 = collect( dfunc(𝑋[i]) - (-1)^(i-1)*h/wfunc(𝑋[i]) for i = 1:length(𝑋))

    f_itp = xx->getp(xx, 𝑓, 𝑤, h, wfunc, dfunc, 𝑋)
    if type_num == 2
        f_itp = xx->getp(xx, 𝑓, 𝑤, h, wfunc, dfunc, 𝑋)/type2weightfunc(xx)

    elseif type_num == 3
        f_itp = xx->getp(xx, 𝑓, 𝑤, h, wfunc, dfunc, 𝑋)/type3weightfunc(xx)

    elseif type_num == 4
        f_itp = xx->getp(xx, 𝑓, 𝑤, h, wfunc, dfunc, 𝑋)/type4weightfunc(xx)

    end

    # get Chebyshev coefficients.
    c, ν = getChebyshevinterpolator(f_itp,length(𝑋)-1)

    # sanity check.
    f_itp_Cheby = xx->evalunivariateChebyshevpolynomial(c,xx)

    xq = convert(Vector{BigFloat}, collect(LinRange(-1,1,200)) )
    yq = f_itp.(xq)
    yq2 = f_itp_Cheby.(xq)
    discrepancy = LinearAlgebra.norm(yq-yq2)
    return c, ν, discrepancy, xq, yq, yq2
end

# does not use BigFloats.
function BarycentrictoChebyshevfast(𝑋, wfunc, dfunc, type_num)

    𝑤 = get𝑤scaled(𝑋)
    h = geth(𝑤, wfunc, dfunc, 𝑋)

    𝑓 = collect( dfunc(𝑋[i]) - (-1)^(i-1)*h/wfunc(𝑋[i]) for i = 1:length(𝑋))

    f_itp = xx->getp(xx, 𝑓, 𝑤, h, wfunc, dfunc, 𝑋)
    if type_num == 2
        f_itp = xx->getp(xx, 𝑓, 𝑤, h, wfunc, dfunc, 𝑋)/type2weightfunc(xx)

    elseif type_num == 3
        f_itp = xx->getp(xx, 𝑓, 𝑤, h, wfunc, dfunc, 𝑋)/type3weightfunc(xx)

    elseif type_num == 4
        f_itp = xx->getp(xx, 𝑓, 𝑤, h, wfunc, dfunc, 𝑋)/type4weightfunc(xx)

    end

    # get Chebyshev coefficients.
    c_SharedArray, ν_SharedArray = getChebyshevinterpolatorfast(f_itp,length(𝑋)-1)

    c = convert(Vector{BigFloat},c_SharedArray)
    ν = convert(Vector{BigFloat},ν_SharedArray)

    # sanity check.
    f_itp_Cheby = xx->clenshaweval(c,xx)

    xq = convert( Vector{BigFloat}, collect(LinRange(-1,1,200)) )

    yq = f_itp.(xq)
    yq2 = f_itp_Cheby.(xq)

    discrepancy = LinearAlgebra.norm(yq-yq2)
    return c, ν, discrepancy, xq, yq, yq2
end

max_stars_repo_name: RoyCCWang/FIRRemez

max_issues_repo_name: RoyCCWang/FIRRemez

max_forks_repo_name: RoyCCWang/FIRRemez


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

# The Computer Language Benchmarks Game
# http://shootout.alioth.debian.org/
#
# Contributed by David Campbell
# Based on the Go version

import Base.isless

function count(data::AbstractString, n::Int)
    counts = Dict{AbstractString, Int}()
    top = length(data) - n + 1
    for i = 1:top
        s = data[i : i+n-1]
        if haskey(counts, s)
            counts[s] += 1
        else
            counts[s] = 1
        end
    end
    counts
end

function count_one(data::AbstractString, s::AbstractString)
    count(data, length(s))[s]
end

type KNuc
    name::AbstractString
    count::Int
end

# sort down
function isless(x::KNuc, y::KNuc)
    if x.count == y.count
        return x.name > y.name
    end
    x.count > y.count
end

function sorted_array(m::Dict{AbstractString, Int})
    kn = Array(KNuc, length(m))
    i = 1
    for elem in m
        kn[i] = KNuc(elem...)
        i += 1
    end
    sort(kn)
end

function print_knucs(a::Array{KNuc, 1})
    sum = 0
    for kn in a
        sum += kn.count
    end
    for kn in a
        @printf("%s %.3f\n", kn.name, 100.0kn.count/sum)
    end
    println()
end

function k_nucleotide(infile="knucleotide-input.txt")
    input = open(infile, "r")
    three = ">THREE "
    while true
        line = readline(input)
        if length(line) >= length(three) && line[1:length(three)] == three
            break
        end
    end
    data = collect(readall(input))
    # delete the newlines and convert to upper case
    i, j = 1, 1
    while i <= length(data)
        if data[i] != '\n'
            data[j] = uppercase(data[i])
            j += 1
        end
        i += 1
    end
    str = join(data[1:j-1], "")

    arr1 = sorted_array(count(str, 1))
    arr2 = sorted_array(count(str, 2))
    close(input)
end


text:
"""
    InvWarpedView(img, tinv, [indices]) -> wv

Create a view of `img` that lazily transforms any given index `I`
passed to `wv[I]` so that `wv[I] == img[inv(tinv)(I)]`.

The conceptual difference to [`WarpedView`](@ref) is that
`InvWarpedView` is intended to be used when reasoning about the
image is more convenient that reasoning about the indices.
Furthermore, `InvWarpedView` allows simple nesting of
transformations, in which case the transformations will be
composed into a single one.

See [`invwarpedview`](@ref) for a convenient constructor of `InvWarpedView`.

For detailed explaination of warp, associated arguments and parameters,
please refer to [`warp`](@ref).
"""
struct InvWarpedView{T,N,A,F,I,FI<:Transformation,E} <: AbstractArray{T,N}
    inner::WarpedView{T,N,A,F,I,E}
    inverse::FI
end

function InvWarpedView(inner::WarpedView{T,N,TA,F,I,E}) where {T,N,TA,F,I,E}
    tinv = _round(inv(inner.transform))
    InvWarpedView{T,N,TA,F,I,typeof(tinv),E}(inner, tinv)
end

function InvWarpedView(A::AbstractArray, tinv::Transformation, inds::Tuple = autorange(A, tinv))
    InvWarpedView(WarpedView(A, inv(tinv), inds), tinv)
end

function InvWarpedView(inner::InvWarpedView, outer_tinv::Transformation)
    tinv = compose(outer_tinv, inner.inverse)
    InvWarpedView(parent(inner), tinv)
end

function InvWarpedView(inner::InvWarpedView, outer_tinv::Transformation, inds::Tuple)
    tinv = compose(outer_tinv, inner.inverse)
    InvWarpedView(parent(inner), tinv, inds)
end

Base.parent(A::InvWarpedView) = parent(A.inner)
@inline Base.axes(A::InvWarpedView) = axes(A.inner)

IndexStyle(::Type{T}) where {T<:InvWarpedView} = IndexCartesian()
@inline Base.getindex(A::InvWarpedView{T,N}, I::Vararg{Int,N}) where {T,N} = A.inner[I...]

Base.size(A::InvWarpedView)    = size(A.inner)
Base.size(A::InvWarpedView, d) = size(A.inner, d)

function Base.showarg(io::IO, A::InvWarpedView, toplevel)
    print(io, "InvWarpedView(")
    Base.showarg(io, parent(A), false)
    print(io, ", ")
    print(io, A.inverse)
    if toplevel
        print(io, ") with eltype ", eltype(parent(A)))
    else
        print(io, ')')
    end
end

"""
    invwarpedview(img, tinv, [indices]; kwargs...) -> wv

Create a view of `img` that lazily transforms any given index `I`
passed to `wv[I]` so that `wv[I] == img[inv(tinv)(I)]`.

Except for the lazy evaluation, the following two lines are equivalent:

```julia
warp(img, inv(tform), [indices]; kwargs...)
invwarpedview(img, tform, [indices]; kwargs...)
```

For detailed explaination of warp, associated arguments and parameters,
please refer to [`warp`](@ref).
"""
function invwarpedview(A::AbstractArray, tinv::Transformation, indices::Tuple=autorange(A, tinv); kwargs...)
    InvWarpedView(box_extrapolation(A; kwargs...), tinv, indices)
end

# For SubArray:
# 1. We can exceed the boundary of SubArray by using its parent and thus trick Interpolations in
#    order to get better extrapolation result around the border. Otherwise it will just fill it.
# 2. For default indices, we use `IdentityUnitRange`, which guarantees `r[i] == i`, to preserve the view indices.
function invwarpedview(A::SubArray, tinv::Transformation; kwargs...)
    default_indices = map(IdentityUnitRange, autorange(CartesianIndices(A.indices), tinv))
    invwarpedview(A, tinv, default_indices; kwargs...)
end
function invwarpedview(A::SubArray, tinv::Transformation, indices::Tuple; kwargs...)
    inner = parent(A)
    new_inner = InvWarpedView(inner, tinv, autorange(inner, tinv))
    view(new_inner, indices...)
end


text:
module JuliaTex

using Formatting
using StatsFuns

export texifycols, textable

const siglev_def = [0.05, 0.01, 0.001]
const dsgn = "\$"
const RealOrVec = Union{T,AbstractVector{T}} where {T<:Real}

function nstars(p::Real, siglev::Vector{<:Real}=siglev_def)
    0.0 <= p <= 1.0 || throw(DomainError("p must be in [0,1]"))
    return "^{" * (rpad("*"^sum(p .< siglev), length(siglev))) * "}"
end

fmt_xstar(x::Real, p::Real, d::Integer, siglev::Vector{<:Real}=siglev_def) = format("{:> 9.$(d)f}", x) * dsgn * nstars(p, siglev) * dsgn
fmt_xstar(x::AbstractVector, p::AbstractVector, d::Integer, siglev::Vector{<:Real}=siglev_def) = [fmt_xstar(xx, pp, d, siglev) for (xx,pp) in zip(x,p)]
padstrvec(x::Vector{String}) = rpad.(x, maximum(length.(x)))

function texifycols(x::RealOrVec, se::RealOrVec, t::RealOrVec, p::RealOrVec; dx::Integer=3, dse::Integer=3, dt::Integer=2, dp::Integer=2, siglev::Vector{<:Real}=siglev_def)
    length(t) == length(p) == length(x) == length(se) || throw(DimensionMismatch())

    strx = padstrvec( fmt_xstar(x, p, dx, siglev) )
    strse = padstrvec( format.("({:.$(dse)f})", se) )
    strt  = padstrvec( format.("{: .$(dt)f}", t)     )
    strp  = padstrvec( format.("{:.$(dp)f}", p)     )

    return (strx, strse, strt, strp)
end

function texifycols(x::RealOrVec, se::RealOrVec; kwargs...)
    length(x) == length(se) || throw(DimensionMismatch())
    t = x ./ se
    p = 2*normcdf.(-abs.(t))
    return texifycols(x, se, t, p; kwargs...)
end

function textable(strcols::NTuple{4,Vector{String}}, nms::Vector{String}=[""]; caption="Caption", label="label", note="", fileheader="", kwargs...)
    strx, strse, strt, strp = strcols
    if length(nms) != length(strx)
        strnms = padstrvec(["\$ \\theta_{$i} \$" for i in 1:length(strx)])
    elseif length(nms) == length(strx)
        strnms = padstrvec(nms)
    else
        throw(DimensionMismatch())
    end
    tbl = strnms .* " & " .* strx .* " & " .* strse .* " & " .* strt .* " & " .* strp .* " \\\\"

    str1 =[
        fileheader,
        "\\begin{table}[ht]",
            "\t\\centering",
            "\t\\caption{$caption}",
            "\t\\label{$label}",
            "\\begin{adjustbox}{max size={.95\\textwidth}{.45\\textheight}}",
            "\t\\begin{tabular}{lcccc}",
            "\t\tParameter & Estimate & SE & \$t\$-statistic & \$p\$-value \\\\",
            "\t\t\\midrule",
        ]
    str2 = "\t\t" .* tbl
    str3 = [
            "\t\t\\midrule",
            "\t\t\\multicolumn{5}{l}{\\scriptsize{\$^{***}p<0.001\$, \$^{**}p<0.01\$, \$^*p<0.05\$}} \\\\",
            length(note) > 0 ? "\t\t\\multicolumn{5}{l}{\\scriptsize $note} \\\\" : "",
            "\t\\end{tabular}",
            "\t\\end{adjustbox}",
        "\\end{table}",
    ]

    strout = join(vcat(str1,str2,str3), "\n")
end

# module end
end


text:
# One step of the heat bath algorithm on a Ising's spin grid
function stepheatbath!(grid::Array{Int, 2};              # Spin grid
                       h::Float64                = 0.0,  # External field
                       temp::Float64             = 1.0,  # Temperature
                       verbose::Bool             = true) # Verbose flag

    @assert temp != 0 "ArgumentError: Temperature can't be zero"

    # Randomly generate a new position within the grid and flip it with the right probability
    x   = rand(1:size(grid, 1))
    y   = rand(1:size(grid, 2))

    # Get the probability of a up spin on (x,y) spot
    m   = nspins(grid, x, y) |> sum
    eplus  = -m - h
    probup = exp(-eplus/temp) / (exp(-eplus/temp) + exp(eplus/temp))

    # Change spin accordingly
    grid[x,y] = rand() < probup ? 1: -1

    if verbose println("Changed $x,$y vertice spin to $(grid[x,y])") end
    return grid
end

# Several steps of the heat bath algorithm on a Ising's spin grid
function heatbath!(grid::Array{Int, 2};               # Spin grid
                   h::Float64                = 0.0,   # External field
                   temp::Float64             = 1.0,   # Temperature
                   iters::Integer            = 50000, # Number of iterations
                   plot::Bool                = true,  # Plot flag
                   verbose::Bool             = true)  # Verbose flag

    @assert iters > 100 "ArgumentError: \"iters\" must be higher than 100"
    m = [stepheatbath!(grid, h=h, temp=temp, verbose=false) |> magnetization for i in 1:iters]

    if verbose println("(T=$temp) Heat bath ended with magnetization $(m[end]) after $(length(m)) iterations") end
    if plot
        PyPlot.plot(1:length(m), m, "o", color="blue")
        PyPlot.plot(1:length(m), m, "-", color="blue")
        PyPlot.title("Heat bath for T=$temp for size $size ($(length(m)) iterations and h=$h)")
        PyPlot.xlabel("Number of Iterations")
        PyPlot.ylabel("Magnetization")
        PyPlot.ylim(0,1.1)
        PyPlot.savefig("heatbath_T=$temp.png")
        PyPlot.close()
    end

    return m
end


text:
@test begin
    acft = Aircraft()

    i1 = add_point!(acft, zeros(3))
    i2 = add_point!(acft, [1.0, 0.0, 0.0])
    i3 = add_point!(acft, [1.0, 1e-2, 0.0])
    i4 = add_point!(acft, [1.0, -1e-2, 0.0])

    add_beam!(
        acft, i1, i2;
        m = 1.0,
        Ixx = 1.0,
        EA = 10.0
    )

    add_mass!(acft, i3, 1.0)
    add_mass!(acft, i4, 1.0)

    add_link!(acft, i1)

    u = get_state(acft, 0.0)
    set_u(u, i2, 1.5)
    qd, _ = state_space(acft, u, 0.0)

    acc_free = get_̇u(qd)

    add_link!(acft, i2, i3)
    add_link!(acft, i2, i4)

    u = get_state(acft, 0.0)
    set_u(u, i2, 1.5)
    qd, _ = state_space(acft, u, 0.0)

    acc_const = get_̇u(qd)

    ms = [
        m.m for m in acft.masses
    ]

    @assert ms ⋅ acc_const ≈ ms ⋅ acc_free

    ϕ, ω = get_eigenmodes(acft, 1)

    @assert isapprox(ω[1], 2.0; rtol = 1e-2)

    true
end

max_stars_repo_name: Equipe-AeroDesign-ITA/WingBiology

max_issues_repo_name: Equipe-AeroDesign-ITA/WingBiology

max_forks_repo_name: Equipe-AeroDesign-ITA/WingBiology


text:
using Test
using XLSX: readxlsx
using DataFrames
using Mimi
using MimiDICE2016
using CSVFiles

using MimiDICE2016: getparams

@testset "MimiDICE2016" begin

    #------------------------------------------------------------------------------
    #   1. Run tests on the whole model
    #------------------------------------------------------------------------------

    @testset "MimiDICE2016-model" begin

        m = MimiDICE2016.get_model();
        run(m)

        f = readxlsx(joinpath(@__DIR__, "../data/DICE2016R-090916ap-v2.xlsm"))

        #Test Precision
        Precision = 1.0e-10

        #Time Periods
        T=100

        #TATM Test (temperature increase)
        True_TATM = getparams(f, "B99:CW99", :all, "Base", T);
        @test maximum(abs, m[:climatedynamics, :TATM] .- True_TATM) ≈ 0. atol = Precision

        #MAT Test (carbon concentration atmosphere)
        True_MAT = getparams(f, "B87:CW87", :all, "Base", T);
        @test maximum(abs, m[:co2cycle, :MAT] .- True_MAT) ≈ 0. atol = Precision

        #DAMFRAC Test (damages fraction)
        True_DAMFRAC = getparams(f, "B105:CW105", :all, "Base", T);
        @test maximum(abs, m[:damages, :DAMFRAC] .- True_DAMFRAC) ≈ 0. atol = Precision

        #DAMAGES Test (damages $)
        True_DAMAGES = getparams(f, "B106:CW106", :all, "Base", T);
        @test maximum(abs, m[:damages, :DAMAGES] .- True_DAMAGES) ≈ 0. atol = Precision

        #E Test (emissions)
        True_E = getparams(f, "B112:CW112", :all, "Base", T);
        @test maximum(abs, m[:emissions, :E] .- True_E) ≈ 0. atol = Precision

        #YGROSS Test (gross output)
        True_YGROSS = getparams(f, "B104:CW104", :all, "Base", T);
        @test maximum(abs, m[:grosseconomy, :YGROSS] .- True_YGROSS) ≈ 0. atol = Precision

        #AL test (total factor productivity)
        True_AL = getparams(f, "B21:CW21", :all, "Base", T);
        @test maximum(abs, m[:totalfactorproductivity, :AL] .- True_AL) ≈ 0. atol = Precision

        #CPC Test (per capita consumption)
        True_CPC = getparams(f, "B126:CW126", :all, "Base", T);
        @test maximum(abs, m[:neteconomy, :CPC] .- True_CPC) ≈ 0. atol = Precision

        #FORCOTH Test (exogenous forcing)
        True_FORCOTH = getparams(f, "B73:CW73", :all, "Base", T);
        @test maximum(abs, m[:radiativeforcing, :FORCOTH] .- True_FORCOTH) ≈ 0. atol = Precision

        #FORC Test (radiative forcing)
        True_FORC = getparams(f, "B100:CW100", :all, "Base", T);
        @test maximum(abs, m[:radiativeforcing, :FORC] .- True_FORC) ≈ 0. atol = Precision

        #Utility Test
        True_UTILITY = getparams(f, "B129:B129", :single, "Base", T);
        @test maximum(abs, m[:welfare, :UTILITY] .- True_UTILITY) ≈ 0. atol = Precision

    end #MimiDICE2016-model testset

    #------------------------------------------------------------------------------
    #   2. Run tests on SCC
    #------------------------------------------------------------------------------

    @testset "Standard API" begin

        m = MimiDICE2016.get_model()
        run(m)

        # Test the errors
        @test_throws ErrorException MimiDICE2016.compute_scc()  # test that it errors if you don't specify a year
        @test_throws ErrorException MimiDICE2016.compute_scc(year=2021)  # test that it errors if the year isn't in the time index
        @test_throws ErrorException MimiDICE2016.compute_scc(last_year=2299)  # test that it errors if the last_year isn't in the time index
        @test_throws ErrorException MimiDICE2016.compute_scc(year=2105, last_year=2100)  # test that it errors if the year is after last_year

        # Test the SCC
        scc1 = MimiDICE2016.compute_scc(year=2020)
        @test scc1 isa Float64

        # Test that it's smaller with a shorter horizon
        scc2 = MimiDICE2016.compute_scc(year=2020, last_year=2200)
        @test scc2 < scc1

        # Test that it's smaller with a larger prtp
        scc3 = MimiDICE2016.compute_scc(year=2020, last_year=2200, prtp=0.02)
        @test scc3 < scc2

        # Test with a modified model
        m = MimiDICE2016.get_model()
        update_param!(m, :t2xco2, 5)
        scc4 = MimiDICE2016.compute_scc(m, year=2020)
        @test scc4 > scc1   # Test that a higher value of climate sensitivty makes the SCC bigger

        # Test compute_scc_mm
        result = MimiDICE2016.compute_scc_mm(year=2030)
        @test result.scc isa Float64
        @test result.mm isa Mimi.MarginalModel
        marginal_temp = result.mm[:climatedynamics, :TATM]
        @test all(marginal_temp[1:findfirst(isequal(2030), MimiDICE2016.model_years)] .== 0.)
        @test all(marginal_temp[findfirst(isequal(2035), MimiDICE2016.model_years):end] .!= 0.)

    end

    # ------------------------------------------------------------------------------
    #   3. Test Deterministic SCC values
    # ------------------------------------------------------------------------------
    @testset "SCC values" begin

        atol = 1e-6 # TODO what is a reasonable tolerance given we test on a few different machines etc.

        # Test several validation configurations against the pre-saved values from MimiDICE2016
        specs = Dict([
            :year => [2020, 2050],
            :eta => [0, 1.5],
            :prtp => [0.015, 0.03],
            :last_year => [2200, 2305],
        ])

        results = DataFrame(year = [], eta = [], prtp = [], last_year = [], SC = [])

        for year in specs[:year]
            for eta in specs[:eta]
                for prtp in specs[:prtp]
                    for last_year in specs[:last_year]
                        sc = MimiDICE2016.compute_scc(year=Int(year), eta=eta, prtp=prtp, last_year=Int(last_year))
                        push!(results, (year, eta, prtp, last_year, sc))
                    end
                end
            end
        end

        validation_results = load(joinpath(@__DIR__, "..", "data", "SC validation data", "deterministic_sc_values_v0-2-0.csv")) |> DataFrame
        # println("MAXIMUM DIFF IS $(maximum(results[!, :SC] - validation_results[!, :SC]))")
        @test all(isapprox.(results[!, :SC], validation_results[!, :SC], atol = atol))

    end # SCC values testset

end #MimiDICE2016 testset

nothing


text:
using Test
function is_ci()
    get(ENV, "TRAVIS", "") == "true" ||
    get(ENV, "APPVEYOR", "") in ("true", "True") ||
    get(ENV, "CI", "") in ("true", "True")
end

if is_ci()
    @info "CI testing is disabled due to the need for a Digilent WaveForms-compatible camera during tests."
else
    using DigilentWaveForms

    devlist = devicelist()
    if length(devlist) < 1
        error("""DigilentWaveForms could not find a compatible device.
        Tests require that a Digilent WaveForms-compatible camera is connected and discoverable.
        If such a device is connected, check that it is discoverable in the WaveForms app""")
    else
        @testset "Device Interaction" begin
            dev = devlist[1]
            @test typeof(dev) == DigilentWaveForms.DWFDevice
        end
    end
end


text:
using PSSFSS
using PSSFSS.Constants: c₀
using PSSFSS.GSMs: GSM
using LinearAlgebra: norm
using Test
using Logging: Error, ConsoleLogger, default_metafmt, global_logger

testlogger = ConsoleLogger(stderr, Error,
                       meta_formatter=default_metafmt, show_limited=true,
                       right_justify=0)
oldlogger = global_logger(testlogger)

"""
    grating(kP, nterms=30) -> (Γ, T)

Compute the normal incidence reflecton and transmission coefficients of a symmetric grid of
zero-thickness conducting strips.  The product of the period of the strips and the incident
electric field wavenumber is `kP` (dimensionless).  The incident electric field is perpendicular
to the direction along the axis of the strips.  The formula is from Problem 10.6 of Collin,
*Field Theory of Guided Waves, Second Ed.*, IEEE Press, 1991.  The series have been accelerated
by appying a Kummer transformation, using the first two terms in the Maclaurin series for the
inverse sin function.  `kP` must be in the half-open interval [0,1). The default number of summed
terms `nterms` yeilds better than 10 digits of accuracy over the interval [0.01,0.99].
"""
function grating(kP;nterms=30)
    sum1 = 1.3862943611198906 # \sum_{n=1}^{\infty} 1/(n-1/2) - 1/n = log(4)
    sum3 = 7.2123414189575710 # \sum_{n=1}^{\infty} (n-1/2)^{-3} - n^{-3} = 6 * \zeta(3)
    x = kP/(4π)
    θ = x*sum1 + x^3/6 * sum3
    for n = 1:nterms
        xonmhalf = x/(n - 0.5)
        xon = x/n
        term = asin(xonmhalf) - (xonmhalf + (xonmhalf)^3/6) -
              (asin(xon) - (xon + xon^3/6))
        θ += term
    end
    Γ = sin(θ) * cis(-π/2 - θ)
    T = 1+Γ
    return (Γ, T)
end

FGHz = 0.8
steering = (θ=0, ϕ=0)
## Compute analytic answer
(r,t) = grating(2π*FGHz)
gsmexact = GSM(2,2)
gsmexact.s11 = [r 0.0; 0.0 -t]
gsmexact.s22 = gsmexact.s11
gsmexact.s21 = [t 0; 0 -r]
gsmexact.s12 = gsmexact.s21

c_inchghz = c₀ * 100/2.54 * 1e-9
period = c_inchghz # so the period/wavelength = freq in GHz
Py = period
Ly = period/2
Px = Lx = Ly/10
Ny = 40
Nx = round(Int, Ny*Lx/Ly)

@testset "JtypeSymmetricStrip" begin

    sheet = rectstrip(;Px, Py, Lx, Ly, Nx, Ny, units=inch)
    strata = [Layer()
              sheet
              Layer()]

    results = analyze(strata, FGHz, steering, logfile=devnull,
                      resultfile=devnull, showprogress=false)
    gsmj = results[1].gsm
    for m in 1:2, n in 1:2
        @test norm(gsmj[m,n] - gsmexact[m,n], Inf) < 0.01
    end
end

@testset "MtypeSymmetricStrip" begin

    sheet = rectstrip(;Px, Py, Lx, Ly, Nx, Ny, units=inch, class='M')
    strata = [Layer()
              sheet
              Layer() ]

    results = analyze(strata, FGHz, steering, logfile=devnull,
                      resultfile=devnull, showprogress=false)
    gsmm = results[1].gsm
    for m in 1:2, n in 1:2
        @test  norm(gsmm[m,n] - gsmexact[m,n], Inf) < 0.01
    end
end

global_logger(oldlogger)
nothing