text:
# extensions to JuMP for multiobjective optimisation

module MultiJuMP

using JuMP
import Ipopt
import JuMP: JuMPTypes, getvalue
import MathProgBase
using RecipesBase, LaTeXStrings
import Combinatorics.combinations
using LinearAlgebra
export MultiModel, SingleObjective, getMultiData

mutable struct SingleObjective
    f # JuMP-expression TODO: use JuMPTypes or something?
    sense::Symbol
    initialvalue::Dict{Symbol,Any} # Variable => Initial value
    # TODO: implement bound in algorithm
    bound::Float64 # Hard lower/upper depending on sense
end

SingleObjective() = SingleObjective(Any, :Min, Dict{Symbol,Any}(), NaN)
SingleObjective(f::JuMPTypes; sense::Symbol = :Min,
                iv::Dict{Symbol,Any} = Dict{Symbol,Any}(),
                bound::Float64 = NaN) = SingleObjective(f, sense, iv, bound)

getvalue(obj::SingleObjective) =  getvalue(obj.f)

function senseValue(obj::SingleObjective)
    # To get Φ and Fstar with right signs
    if obj.sense == :Max
        return -getvalue(obj.f)
    else
        return getvalue(obj.f)
    end
end

getvalue(arr::Array{SingleObjective}) = map(getvalue, arr)
senseValue(arr::Array{SingleObjective}) = map(senseValue, arr)

# stores extension data inside JuMP Model
mutable struct MultiData{Tx,To}
    objectives::Vector{SingleObjective}
    f1::SingleObjective
    f2::SingleObjective
    normalf1
    normalf2

    # Number of points in each direction
    # of Pareto submanifold
    pointsperdim::Int
    #
    # stored values
    utopiavarvalues::Vector{Vector{Tx}}
    utopia::Vector{To}
    nadir::Vector{Tx}
    paretovarvalues::Vector{Vector{To}}
    paretofront

    Phi::Array{To,2}
    Fstar::Vector{To}
end

function MultiModel(;solver=Ipopt.IpoptSolver())
    m = Model(solver=solver)
    m.solvehook = solvehook
    m.ext[:Multi] = MultiData(Array{SingleObjective}(undef, 0),
                              SingleObjective(), SingleObjective(),
                              Any, Any,
                              10,
                              Vector{Float64}[], Float64[],
                              Float64[], Vector{Float64}[], Any[],
                              Array{Float64}(undef, 2,2), Array{Float64}(undef, 2))
    return m
end

function getMultiData(m::Model)
    if haskey(m.ext, :Multi)
        return m.ext[:Multi]::MultiData
    else
        error("This functionality is only available for MultiJuMP models")
    end
end

function betas(levels,parts)
    # Sets up the tree of all possible
    # convex combinations of the objectives
    map(x->([x;levels+parts].-[0;x] .-1)/parts,
        combinations(1:(levels+parts-1),(levels-1)))
end

function _solve_ws(m::Model)
    multim = getMultiData(m)
    objectives = multim.objectives
    sensemap = Dict(:Min => 1.0, :Max => -1.0)

    vararr = [JuMP.Variable(m,i) for i in 1:MathProgBase.numvar(m)]

    numobj = length(objectives)
    Phi = zeros(numobj,numobj)

    # Individual minimisations
    for (i, objective) in enumerate(objectives)
        @NLobjective(m, objective.sense, objective.f)
        for (key, value) in objective.initialvalue
            # TODO: What is the correct way to set these values?
            setvalue(m.objDict[key], value)
        end
        status = solve(m, ignore_solve_hook=true);
        if status != :Optimal
            return status
        end

        push!(multim.utopiavarvalues, [getvalue(var) for var in vararr])
        push!(multim.paretovarvalues, [getvalue(var) for var in vararr])

        Phi[:,i] = senseValue(objectives)

        push!(multim.paretofront, getvalue(objectives))
    end
    Fmax = maximum(Phi, dims=2)
    Fmin = minimum(Phi, dims=2) # == diag(Phi)?
    if Fmax == Fmin
        error("The Nadir and Utopia points are equal") # I think that's what this means?
    end

    multim.Phi = Phi

    beta = zeros(numobj); beta[end] = 1.0
    @NLparameter(m, β[i=1:numobj] == beta[i])

    @NLobjective(m, :Min,
                 sum(β[i]*(sensemap[objectives[i].sense]*objectives[i].f -
                           Fmin[i])/(Fmax[i]-Fmin[i]) for i=1:numobj))

    betatree = betas(numobj, multim.pointsperdim-1)

    for betaval in betatree
        if count(t -> t != 0, betaval) == 1
            # Skip individual optimisations as
            # they are already performed
            continue
        end

        setvalue(β, betaval)

        status = solve(m, ignore_solve_hook=true);
        if status != :Optimal
            return status
        end

        push!(multim.paretofront, getvalue(objectives))
        push!(multim.paretovarvalues, [getvalue(var) for var in vararr])
    end

    return :Optimal
end

function _solve_nbi(m::Model, inequalityconstraint::Bool = false)
    multim = getMultiData(m)
    objectives = multim.objectives
    sensemap = Dict(:Min => 1.0, :Max => -1.0)

    vararr = [JuMP.Variable(m,i) for i in 1:MathProgBase.numvar(m)]

    # Stage 1: Calculate Φ
    numobj = length(objectives)
    Fstar = zeros(numobj)
    Phi = zeros(numobj,numobj)

    # Individual minimisations
    for (i, objective) in enumerate(objectives)
        @NLobjective(m, objective.sense, objective.f)
        for (key, value) in objective.initialvalue
            # TODO: What is the correct way to do this? (anonymous variables)
            setvalue(m.objDict[key], value)
        end
        status = solve(m, ignore_solve_hook=true);
        if status != :Optimal
            return status
        end

        push!(multim.utopiavarvalues, [getvalue(var) for var in vararr])
        push!(multim.paretovarvalues, [getvalue(var) for var in vararr])

        Phi[:,i] = senseValue(objectives)

        push!(multim.paretofront, getvalue(objectives))
    end
    Fstar = diag(Phi)

    for j = 1:numobj
        Phi[:,j] = Phi[:,j] - Fstar
    end
    multim.Phi = Phi
    multim.Fstar = Fstar

    # Stage 2: Create NBI subproblems
    @variable(m, t)
    # TODO: There is a bug in JuMP so it doesn't propagate all
    # the necessary information if we use @objective instead of NLObjective
    # TODO: test this with nlprewrite
    @NLobjective(m, :Max, t)

    beta = zeros(numobj); beta[end] = 1.0
    @NLparameter(m, β[i=1:numobj] == beta[i])

    if inequalityconstraint == false
        # Standard NBI
        for (i, objective) in enumerate(objectives)
            @NLconstraint(m,
                          sum(Phi[i,j]*(β[j]-t)
                              for j = 1:numobj if j != i) ==
                          sensemap[objective.sense]*objective.f-Fstar[i])
        end
    else
        # Pascoletti-Serafini extension
        for (i, objective) in enumerate(objectives)
            @NLconstraint(m,
                          sum(Phi[i,j]*(β[j]-t)
                              for j=1:numobj if j != i) >=
                          sensemap[objective.sense]*objective.f-Fstar[i])
        end
    end

    # Stage 3: Solve NBI subproblems
    betatree = betas(numobj, multim.pointsperdim-1)

    for betaval in betatree
        if count(t -> t != 0, betaval) == 1
            # Skip individual optimisations as
            # they are already performed
            continue
        end
        setvalue(β, betaval)
        status = solve(m, ignore_solve_hook=true);
        if status != :Optimal
            return status
        end

        push!(multim.paretofront, getvalue(objectives))
        push!(multim.paretovarvalues, [getvalue(var) for var in vararr])
    end

    return :Optimal
end

function _solve_eps(m::Model)
    multim = getMultiData(m)
    objectives = multim.objectives
    sensemap = Dict(:Min => 1.0, :Max => -1.0)

    vararr = [JuMP.Variable(m,i) for i in 1:MathProgBase.numvar(m)]

    numobj = length(objectives)
    Phi = zeros(numobj,numobj)

    if numobj > 2
        # TODO:
        # The logic  here becomes difficult, as the feasible region will require
        # a dependency between the constraints
        Base.error(":EPS is thought through for > 2 objectives yet")
    end

    # Individual minimisations
    for (i, objective) in enumerate(objectives)
        @NLobjective(m, objective.sense, objective.f)
        for (key, value) in objective.initialvalue
            # TODO: What is the correct way to do this? (anonymous variables)
            setvalue(m.objDict[key], value)
        end
        status = solve(m, ignore_solve_hook=true);
        if status != :Optimal
            return status
        end

        push!(multim.utopiavarvalues, [getvalue(var) for var in vararr])
        push!(multim.paretovarvalues, [getvalue(var) for var in vararr])

        Phi[:,i] = senseValue(objectives)

        push!(multim.paretofront, getvalue(objectives))
    end
    Fmax = maximum(Phi, dims=2)
    Fmin = minimum(Phi, dims=2) # == diag(Phi)?

    multim.Phi = Phi

    @NLobjective(m, objectives[end].sense, objectives[end].f)

    beta = zeros(numobj); beta[end] = 1.0
    @NLparameter(m, β[i=1:numobj] == beta[i])

    @NLconstraint(m, objconstr[i=1:numobj-1],
                  sensemap[objectives[i].sense]*objectives[i].f
                  <= β[i]*Fmin[i]+(1-β[i])*Fmax[i])

    betatree = betas(numobj, multim.pointsperdim-1)

    for betaval in betatree
        if count(t -> t != 0, betaval) == 1
            # Skip individual optimisations as
            # they are already performed
            continue
        end

        setvalue(β, betaval)

        status = solve(m, ignore_solve_hook=true);
        if status != :Optimal
            return status
        end

        push!(multim.paretofront, getvalue(objectives))
        push!(multim.paretovarvalues, [getvalue(var) for var in vararr])
    end

    return :Optimal
end

function solvehook(m::Model; method::Symbol = :NBI,
                   inequalityconstraint::Bool = false, kwargs...)
    if method == :WS
        # Weighted sums
        status = _solve_ws(m)
    elseif method == :EPS
        # Epsilon constraint method
        status = _solve_eps(m)
    elseif method == :NBI
        # Normal boundary intersection
        status = _solve_nbi(m, inequalityconstraint)
    else
        Base.error("Multiobjective method not recognized.")
    end

    return status
end

@recipe function f{Tx,To}(md::MultiData{Tx,To})
    seriestype --> :scatter
    xlabel --> L"$f_1$"
    ylabel --> L"$f_2$"
    zlabel --> L"$f_3$"
    label --> ""
    numobjectives = length(md.objectives)
    if numobjectives > 3 || numobjectives < 2
        Base.error("Only plotting 2d and 3d fronts")
    end

    f1arr = convert(Array{Float64},
                    [val[1] for val in md.paretofront])
    f2arr = convert(Array{Float64},
                    [val[2] for val in md.paretofront])

    title --> "Pareto front with $(length(f1arr)) points"

    if numobjectives == 2
        xyz = (f1arr, f2arr)
    else
        f3arr = convert(Array{Float64},
                        [val[3] for val in md.paretofront])
        xyz = (f1arr, f2arr, f3arr)
    end
    (xyz...,)
end

@recipe function f(m::JuMP.Model)
    plot(getMultiData(m))
end

end


text:
@testset "7.1 Inverse hyperbolic sine" begin
    include("7.1.2 (d x)^m (a+b arcsinh(c x))^n.jl")
    include("7.1.4 (f x)^m (d+e x^2)^p (a+b arcsinh(c x))^n.jl")
    include("7.1.5 Inverse hyperbolic sine functions.jl")
end

ext: jl

max_stars_repo_path: test/7 Inverse hyperbolic functions/7.1 Inverse hyperbolic sine/.jl


text:
using Peacock
using PyPlot, SpecialFunctions, LinearAlgebra
r = 0.2
epf(x, y) = x^2 + y^2 < r^2 ? 9 : 1
muf(x,y) = 1
a1 = [1, 0]; a2 = [0, 1]
d1 = 0.001; d2 = 0.001
geometry = Geometry(epf, muf, a1, a2, d1, d2)
solver = Solver(geometry, 7, 7)
m = 3

kx = [i for i in -m:m]
kx = repeat(kx,1, 2*m+1)
ky = copy(kx)
kx = kx'
K = hcat(kx[:],ky[:])#为何没有*2pi
Kp = K

Kxp = [i for i in Kp[:,1]]
Kxp = repeat(Kxp,1,length(Kp[:,1]))
Kx = copy(Kxp)
Kxp = Kxp'

Kyp = [i for i in Kp[:,2]]
Kyp = repeat(Kyp,1,length(Kp[:,1]))
Ky = copy(Kyp)
Kyp = Kyp'

function geo(kx,ky)
    k=sqrt.(kx.^2 .+ ky.^2)
    eb = 1
    eab=8
    r=0.2
    ep=eab*r*besselj.(1,2*pi*k*r)./k.*cos.(pi*kx)
    ep[k.==0] .= (eb+eab*pi*r^2)
    return ep
end

ep1 = solver.epc
#ep1 = geo(Kx-Kxp, Ky-Kyp)

omegaTM = []
omegaTE = []
for l=0:0.05:1.5;
    #0:0.25:1.5;
    if l<=0.5
        k=[1 0]*l;
    elseif l<=1;
        k=[0.5 l-0.5];
    else
        k=[0.5 0.5]-[1 1]*(l-1);
    end

    #find \omega^2
    kk=ones(length(K[:,1]))*k
    #和lecture是等价的，相当于：kx*kx+ky*ky
    omegaTM2, HkTM=eigen(4*pi^2*(kk+K)*((kk+Kp)').*inv(ep1),sortby=real)
    temp = (kk+K)*(kk+Kp)'
    temp = temp[diagind(temp)]
    temp = diagm(temp)
    omegaTE2, HkTE=eigen(temp,ep1,sortby=real)

#    omegaTM = [omegaTM, sqrt.(omegaTM2)]
    push!(omegaTM,sqrt.(Complex.(omegaTM2)))
    push!(omegaTE,sqrt.(Complex.(omegaTE2)))
end
figure()
plot(real.(omegaTM)/2/pi) #这里也没有÷2pi
ylim([0,0.8])


text:
using AccuracyAtTopKernels

using Test, Distributions, LinearAlgebra, Random
import AccuracyAtTopKernels.MLKernels

include("./tests_surrogates.jl")
include("./tests_datasets.jl")
include("./tests_projections.jl")
include("./tests_primal_problems.jl")
include("./tests_dual_problems.jl")
include("./tests_kernels.jl")
include("./tests_scores.jl")

@time @testset "all tests" begin
    @testset "tests surrogates" begin
        test_hinge()
        test_quadratic()
    end

    @testset "tests datasets" begin
        test_primal()
        test_dual()
    end

    @testset "tests projections" begin
        test_projections()
    end

    @testset "tests primal problems" begin
        test_primal_problems()
    end

    @testset "tests dual problems" begin
        test_dual_problems()
    end

    @testset "tests kernels" begin
        test_kernels()
    end

    @testset "tests scores" begin
        test_scores()
    end
end;


text:
using DataFrames
using CSV
using CSVFiles # for providing consistent reading syntax with ExcelFiles
using ExcelFiles
using DelimitedFiles
using XLSX
using BenchmarkTools
using Plots
using StatsPlots

##
## TODO: BUILD AN AFTER ANALYSIS
# [X] BUILD a DATAFRAME for result
# [X]BUILD a DATAFRAME for models parameters to test the running time
    # [X] SET UP EXCEL FiLE FOR DATA INPUT
    # [X] WRITE CODE THAT INPUT THE SETS AND PARAM DIRECTLY FROM THE EXCEL FILE
# [X] CODE FOR RUNNING MULTPLE ROUNDS
    # [X] MULTPLE INPUT
    # [X] MULTPLE OUTPUT
# [X]BUILD Plot to show the sensitivity of the running time regardings the adjusting parameters
    # [X] Create DataFrames - ALL VARS
    # [X] Create DataFrames - ALL STATS
    # [X] Plot Line chart for modeling time and solving time

###################### BUILDING DATAFRAMES #################################################
model_first = 0;
model_last = 3;
#### Variable DataFrames ####
var_Dict = Dict();
var_df = DataFrame(load("output_var0.csv"));
for i in (model_first:model_last)
    var_Dict["scenario$i"] = DataFrame(load("output_var$i.csv"))
end;

for i in (model_first+1:model_last)
    var_df = outerjoin(var_df, var_Dict["scenario$i"],on= [:variable_type,:index],makeunique=true)
end;

var_df = rename(var_df, :value => :scenario_0);
for i in (model_first+1:model_last)
    var_df = rename(var_df, "value_$i" => "scenario_$i")
end;

println(var_df);
# Export the merged dataframe
CSV.write("ScenarioCompare_vars.csv",var_df);

#### Stats DataFrames ####
stat_Dict = Dict();
stat_df = DataFrame(load("output_stat0.csv"));
for i in (model_first:model_last)
    stat_Dict["scenario$i"] = DataFrame(load("output_stat$i.csv"))

end;

for i in (model_first+1:model_last)
    stat_df = append!(stat_df, stat_Dict["scenario$i"])
end;

stat_df[!,:scenario] = convert.(Int,stat_df[!,:scenario,]);

println(stat_df);
# Export the merged stat dataframe
CSV.write("ScenarioCompare_stats.csv",stat_df);

########################## Plotting #################################################
#### DOES NOT WORK ANYMORE!!!! REVISION NEEDED###############
## Interactive Data Exploration (Uncomment to use!)
# stat_df |> Voyager()
# var_df  |> Voyager()
# stat_df
# p1_sce_sol = stat_df |> @vlplot( :point, x=:scenario, y=:solve_time)
# p2_sce_mod = stat_df |> @vlplot( :point, x=:scenario, y=:modeling_time) # not really correct as zero always have a higher model building time
# p3_sce_obj = stat_df |> @vlplot( :point, x=:scenario, y=:objective_time) # something is wrong with y axis
# p4_sce_sol = stat_df |> @vlplot( :line, x=:scenario, y=:solve_time) # maybe will be more useful when have more scenarios
stat_df
@df stat_df plot(:scenario,:objective_value)
@df stat_df plot(:scenario,:solve_time)
@df stat_df plot(:scenario,:modeling_time)


text:
module StatsPlots

using Reexport
import RecipesBase: recipetype
import Tables: istable, columntable, select, schema, rows
@reexport using Plots
import Plots: _cycle
using Plots.PlotMeasures
using StatsBase
using Distributions
import DataValues: DataValue
using Widgets, Observables
import Observables: AbstractObservable, @map, observe
import Widgets: @nodeps
import DataStructures: OrderedDict
import Clustering: Hclust, nnodes
using Interpolations
import MultivariateStats: MDS, eigvals, projection, principalvars,
                            principalratio, transform

import KernelDensity
@recipe f(k::KernelDensity.UnivariateKDE) = k.x, k.density
@recipe f(k::KernelDensity.BivariateKDE) = k.x, k.y, permutedims(k.density)

@shorthands cdensity

export @df, dataviewer

include("df.jl")
include("interact.jl")
include("corrplot.jl")
include("cornerplot.jl")
include("distributions.jl")
include("boxplot.jl")
include("dotplot.jl")
include("violin.jl")
include("hist.jl")
include("marginalhist.jl")
include("bar.jl")
include("dendrogram.jl")
include("andrews.jl")
include("ordinations.jl")

end # module


text:
using ExcelFiles
using Mimi
using Random
using Statistics
using Test

include("../../wip/load_empirical_dist.jl")

# function load_vector(path, range, header=false)
#     tups = collect(load(path, range, header=header))
#     name = fieldnames(tups[1])[1]   # field name of first item in NamedTuple
#     map(obj -> getfield(obj, name), tups)
# end

# function load_empirical_dist(path::AbstractString,
#                              values_range::AbstractString,
#                              probs_range::AbstractString="")
#     println("Reading from '$path'")
#     values = load_vector(path, values_range)
#     probs = probs_range == "" ? nothing : load_vector(path, probs_range)
#     d = Mimi.EmpiricalDistribution(values, probs)
#     println("returning distribution $(typeof(d))")
#     return d
# end

filename = joinpath(@__DIR__, "RB-ECS-distribution.xls")
d = load_empirical_dist(filename, "Sheet1!A2:A1001", "Sheet1!B2:B1001")

# Set the seed to get repeatable results (with some caveats...)
Random.seed!(1234567)

# Get the statistical outline of the distribution
q = quantile.(Ref(d), [0.01, 0.10, 0.25, 0.50, 0.75, 0.90, 0.99])

expected = [1.47, 1.92, 2.33, 3.01, 4.16, 5.85, 9.16]
# @info "quantiles: $q\n expected: $expected"

@test isapprox(q, expected, atol=0.01)

# Test that the EmpiricalDistribution gets saved as a SampleStore and values
# get re-used across multiple scenarios

num_scenarios = 4
num_trials = 5
output_dir = "./out/"

results = zeros(num_scenarios, num_trials)

_values = collect(1:10)
_probs = 0.1 * ones(10)

sd = @defsim begin
    sampling(LHSData)
    p = EmpiricalDistribution(_values, _probs)
end

@defcomp test begin
    p = Parameter(default = 5)
    function run_timestep(p, v, d, t) end
end

scenario_args = [
    :num => collect(1:num_scenarios)
]

function scenario_func(sim_inst::SimulationInstance{T}, tup::Tuple) where T <: AbstractSimulationData
    nothing
end

function post_trial_func(sim_inst::SimulationInstance{T}, trialnum::Int, ntimesteps::Int,
                         tup::Tuple)  where T <: AbstractSimulationData
    m, = sim_inst.models
    scenario_num, = tup
    results[scenario_num, trialnum] = m[:test, :p]
end

m = Model()
set_dimension!(m, :time, 2000:2001)
add_comp!(m, test)

si = run(sd, m, num_trials;
    results_output_dir = output_dir,
    scenario_args = scenario_args,
    scenario_func = scenario_func,
    post_trial_func = post_trial_func
    )

for rv in values(si.sim_def.rvdict)
    @test rv.dist isa Mimi.SampleStore
end

for i = 1:num_scenarios
    @test results[i, :] == results[1, :]
end

rm(output_dir, recursive = true)

si1 = run(sd, m, num_trials)
trial1 = copy(collect(values(si1.sim_def.rvdict))[1].dist.values)

for rv in values(si1.sim_def.rvdict)
    @test rv.dist isa Mimi.SampleStore
end

si2 = run(sd, m, num_trials)
trial2 = copy(collect(values(si2.sim_def.rvdict))[1].dist.values)

@test trial1!=trial2


text:
# Author: Youngjun Kim, youngjun@stanford.edu
# Date: 06/25/2014

if !isdefined(:PACKAGE_PATH)
    PACKAGE_PATH = Pkg.dir("SISLES", "src")
end

push!(LOAD_PATH, "$PACKAGE_PATH/include")
push!(LOAD_PATH, "$PACKAGE_PATH/common")
push!(LOAD_PATH, "$PACKAGE_PATH/Encounter")
push!(LOAD_PATH, "$PACKAGE_PATH/Encounter/CorrAEMImpl")
push!(LOAD_PATH, "$PACKAGE_PATH/Encounter/LLAEMImpl")
push!(LOAD_PATH, "$PACKAGE_PATH/PilotResponse")
push!(LOAD_PATH, "$PACKAGE_PATH/PilotResponse/SimplePilotResponseImpl")
push!(LOAD_PATH, "$PACKAGE_PATH/DynamicModel")
push!(LOAD_PATH, "$PACKAGE_PATH/DynamicModel/SimpleADMImpl")
push!(LOAD_PATH, "$PACKAGE_PATH/WorldModel")
push!(LOAD_PATH, "$PACKAGE_PATH/WorldModel/AirspaceImpl")
push!(LOAD_PATH, "$PACKAGE_PATH/Sensor")
push!(LOAD_PATH, "$PACKAGE_PATH/Sensor/SimpleTCASSensorImpl")
push!(LOAD_PATH, "$PACKAGE_PATH/CollisionAvoidanceSystem")
push!(LOAD_PATH, "$PACKAGE_PATH/CollisionAvoidanceSystem/SimpleTCASImpl")
push!(LOAD_PATH, "$PACKAGE_PATH/Simulator")
push!(LOAD_PATH, "$PACKAGE_PATH/Simulator/TCASSimulatorImpl")

using Util
using Encounter
using PilotResponse
using DynamicModel
using WorldModel
using Sensor
using CollisionAvoidanceSystem
using Simulator


text:
function extract_kwargs(kwargs, keys)
    return (k => v for (k, v) in pairs(kwargs) if k in keys)
end


text:
using MXNet

#--------------------------------------------------------------------------------
# define lenet

# input
data = mx.Variable(:data)

# first conv
conv1 = @mx.chain mx.Convolution(data=data, kernel=(5,5), num_filter=20)  =>
                  mx.Activation(act_type=:tanh) =>
                  mx.Pooling(pool_type=:max, kernel=(2,2), stride=(2,2))

# second conv
conv2 = @mx.chain mx.Convolution(data=conv1, kernel=(5,5), num_filter=50) =>
                  mx.Activation(act_type=:tanh) =>
                  mx.Pooling(pool_type=:max, kernel=(2,2), stride=(2,2))

# first fully-connected
fc1   = @mx.chain mx.Flatten(data=conv2) =>
                  mx.FullyConnected(num_hidden=500) =>
                  mx.Activation(act_type=:tanh)

# second fully-connected
fc2   = mx.FullyConnected(data=fc1, num_hidden=10)

# softmax loss
lenet = mx.SoftmaxOutput(data=fc2, name=:softmax)

#--------------------------------------------------------------------------------
# load data
batch_size = 100
include("mnist-data.jl")
train_provider, eval_provider = get_mnist_providers(batch_size; flat=false)

#--------------------------------------------------------------------------------
# fit model
model = mx.FeedForward(lenet, context=mx.gpu())

# optimizer
optimizer = mx.SGD(lr=0.05, momentum=0.9, weight_decay=0.00001)

# fit parameters
mx.fit(model, optimizer, train_provider, n_epoch=20, eval_data=eval_provider)


text:
export UniformMassDistribution

"Initial uniform mass distribution."
struct UniformMassDistribution{N} <: InitialMassDistribution{N}
    shape::SVector{N,Int}
    lower_bounds::SVectorF{N}
    cell_extent::Float
    ρ::Float

    """
        UniformMassDistribution(
            bounds::Bounds{N,<:Number},
            ρ::Number,
            number_of_particles::Integer,
        )

    Create a new uniform mass distribution with the given mass density `ρ`
    between the given bounds, with the given number of particles.
    """
    function UniformMassDistribution(
        bounds::Bounds{N,<:Number},
        ρ::Number,
        number_of_particles::Integer,
    ) where {N}
        @check ρ > 0
        @check number_of_particles > 0
        grid_extents = extents(bounds)
        volume = prod(grid_extents)
        cell_extent = (volume / number_of_particles)^(1 / N)
        shape =
            max.(
                ones(SVector{N,Int}),
                SVector{N,Int}(round.(grid_extents / cell_extent)),
            )
        new_extents = shape * cell_extent
        extensions = new_extents - grid_extents
        new_bounds = extended(bounds, extensions)
        new_lower_bounds = lower(new_bounds)
        new_cell_extent = (prod(new_extents) / prod(shape))^(1 / N)
        new{N}(shape, new_lower_bounds, new_cell_extent, ρ)
    end
end

Base.show(io::IO, ::Type{UniformMassDistribution{N}}) where {N} =
    print(io, "UniformMassDistribution")

function (distribution::UniformMassDistribution{N})() where {N}
    number_of_particles = prod(distribution.shape)
    m = distribution.ρ * distribution.cell_extent

    coordinates = [
        range(
            distribution.lower_bounds[dim] + 0.5 * distribution.cell_extent,
            step = distribution.cell_extent,
            length = distribution.shape[dim],
        ) for dim = 1:N
    ]
    indices = CartesianIndices(Tuple(distribution.shape))
    positions = [
        SVectorF{N}([coordinates[dim][indices[i][dim]] for dim = 1:N]) for i = 1:number_of_particles
    ]

    mass_densities = fill(distribution.ρ, number_of_particles)

    MassDistributionData(m, positions, mass_densities)
end

max_stars_repo_name: lars-frogner/Whirl

max_issues_repo_name: lars-frogner/Whirl

max_forks_repo_name: lars-frogner/Whirl


text:
using BinaryProvider

dependencies = Dict(
    "https://github.com/JuliaMath/GMPBuilder/releases/download/v6.1.2-2/build_GMP.v6.1.2.jl" =>
        "91b2197c4d2209985b2c156e04f3fc11e4beff633cd86bb0053bfefe71cdfba6",

    "https://github.com/JuliaMath/MPFRBuilder/releases/download/v4.0.1-3/build_MPFR.v4.0.1.jl" =>
        "40f3d06abf45b571b23135a079d66827282e587be15d540874bb45dac163e2c7",

    "https://github.com/isuruf/MPCBuilder/releases/download/v1.1.0-2/build_MPC.v1.1.0.jl" =>
        "85cac0057832da9c9d965531e9a1bada7150032aea4dbead59ff76e95bbdc47f",

    "https://github.com/symengine/SymEngineBuilder/releases/download/v0.3.0-2/build_SymEngine.v0.3.0.jl" =>
        "664b7df2b2e173625fa5742aa194e63392692489b73e6ba4005dcbf661093c9d",
)

prefix = joinpath(@__DIR__, "symengine-0.3")
downloads_dir = joinpath(prefix, "downloads")
all_products = LibraryProduct[]

for (url, hash) in dependencies
    tmp_file = joinpath(downloads_dir, "build_$hash.jl")
    download_verify(url, hash, tmp_file)
    contents = read(tmp_file, String)
    m = Module(:__anon__)
    Core.eval(m, quote
        using BinaryProvider
        function write_deps_file(path, products) end
        ARGS = [$prefix]
    end)
    Base.include_string(m, contents)
    products = Core.eval(m, :(products))
    append!(all_products, products)
end

for product in all_products
    locate(product, verbose=true)
end

# Write out a deps.jl file that will contain mappings for our products
write_deps_file(joinpath(@__DIR__, "deps.jl"), all_products)


text:
using DataFrames, CSV, Distributions

# Test Data
    # source_data = rand(Normal(), 1000, 5)
    # FileName = "sip_export"

"""
    dfnames(headings)

Adds ":" in front of all items in an array. This format is required by the DataFrame packages when giving names to columns.
"""
function dfnames(headings)
h_list = []
    for i = 1:size([headings...],1)
        name = ":" * headings[i]
        println(name)
        push!(h_list, name)

    end
    return h_list
end

"""
    sip2csv(FileName, source_df, s_name="SLURP", s_origin="Julia Language")

- In order to export a SIP Library from Julia, you simply need to have a DataFrame.
- Only specify header fields if your data does not contain any.
- Also make to specify the full filename, including extension
"""
function sip2csv(FileName, source_df, s_name="SLURP", s_origin="Julia Language")

    csvfilename = "$FileName.csv"
    csvfile = open(csvfilename, "w") #Create CSV csvfile
    numSamples = size(source_df, 1)
    colCount = size(source_df, 2)
    df_names = SIPNames(source_df)

    #colwise methods
    meta_mean = lineprep([[mean(col) for col = eachcol(source_df)]])
    meta_median = lineprep([[median(col) for col = eachcol(source_df)]])
    meta_min = lineprep([[minimum(col) for col = eachcol(source_df)]])
    meta_max = lineprep([[maximum(col) for col = eachcol(source_df)]])

    #Load MetaData Table
    metafilename = "$s_name"*"_meta.csv"

    if isfile(metafilename) == true
        meta_df = CSV.read(metafilename, datarow=4)
        AttributeCount = size(meta_df, 1)
        SipAttrsEnd = 23 + AttributeCount
        SipTlc = SipAttrsEnd + 4

        #metadata = convert(Matrix, meta_df)
        # Setup SIP Header
        write(csvfile,"CSV,,,,,,Always have something in A1,,,,,,,,,,,,,,,,,,,,,,,,,,,
        ,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,
        Control,,SheetName,N/A,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,
        ,,FilePath,$FileName.csv,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,
        ,,SlurpAttrs,C14:C16,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,
        ,,SipAttrs,C19:C$SipAttrsEnd,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,
        ,,SipTlc,A$SipTlc,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,
        ,,ClearFirst,TRUE,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,
        ,,numSamples,$numSamples,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,
        ,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,
        ,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,
        ,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,
        ,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,
        SLURP,,name,$s_name,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,
        ,,count,$colCount,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,
        ,,origin,$s_origin,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,
        ,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,
        ,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,
        SIP,,name,$df_names
        ,,mean,$meta_mean
        ,,median,$meta_median
        ,,min,$meta_min
        ,,max,$meta_max,\n")
        close(csvfile)

        CSV.write("$FileName.csv", meta_df, append=true)
        csvfile = open(csvfilename, "a")
        write(csvfile,",,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,
        ,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,\n")

        close(csvfile)

        CSV.write("$FileName.csv", source_df, append=true, writeheader=true)

    elseif isfile(metafilename) == false

        # Setup SIP Header
        write(csvfile,"CSV,,,,,,Always have something in A1,,,,,,,,,,,,,,,,,,,,,,,,,,,
        ,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,
        Control,,SheetName,N/A,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,
        ,,FilePath,$FileName.csv,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,
        ,,SlurpAttrs,C14:C16,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,
        ,,SipAttrs,C19:C23,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,
        ,,SipTlc,A27,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,
        ,,ClearFirst,TRUE,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,
        ,,numSamples,$numSamples,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,
        ,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,
        ,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,
        ,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,
        ,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,
        SLURP,,name,$s_name,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,
        ,,count,$colCount,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,
        ,,origin,$s_origin,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,
        ,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,
        ,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,
        SIP,,name,$df_names
        ,,mean,$meta_mean
        ,,median,$meta_median
        ,,min,$meta_min
        ,,max,$meta_max
        ,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,
        ,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,\n")
        close(csvfile)
        CSV.write("$FileName.csv", source_df, append=true, writeheader=true)
    end

end

"""
    SIPNames(source_df)

Preps a vector of names for the CSV export of SIPs by removing unnecessary characters in front of all items in an array.
"""
function SIPNames(source_df)
    df_names=names(source_df)
    df_names=string([df_names])
    df_names = replace(df_names, "Array{Symbol,1}[[:" =>"")
    df_names = replace(df_names, "]]" =>"")
    df_names = replace(df_names, ":" =>"")
    return df_names
end
"""
    lineprep(line_array)

When you want to strip "[]" from Arrays for export to CSV

"""
function lineprep(line_array)
    val = string([line_array])
    val = replace(val, "Array{Array{Float64,1},1}" =>"")
    val = replace(val, "[" =>"")
    val = replace(val, "]" =>"")
    return val
end

#Loading a SIPloadsip
"""
    importxlsip(FileName, source="")

This function allows to import SIPs in CSV format from Excel using the SIP 2.0 Standard (Stochastic Information Packets, Savage[2009]) to build unified simulation models.
Based on the CSV package, the function cleans up meta data and cleans out redundant columns. You can also set a *source string* that will be included in your SIP DataFrame as a dimension.
"""
function importxlsip(FileName, source="")
println("Processing SIP Data...\n")
    open(FileName) do file
        linecounter = 1
            for ln in eachline(file)
                if occursin("SIP,,name,", ln) == true
                    println("header_row: $linecounter")
                    global header_row = linecounter
                    linecounter +=1
                elseif occursin("SIPs,,,", ln) == true
                    println("data_row: $linecounter")
                    global start_row = linecounter
                else linecounter +=1
                end
            end
    end
    println("")
    println("SIPs are acceccisble via Dataframe. \n")
    println("header=$header_row, datarow=$start_row")
    SLURP = CSV.read(FileName, header=header_row, datarow=start_row)
    select!(SLURP, Not([:SIP, :Column2, :name])) #replaces deletecols!
    #SLURP[!, :source] .= source # new syntax. now uses .functions to broadcast
    println("")
    print(describe(SLURP))
    println("\n")
    return SLURP
end

"""
    importsip(FileName, source="Not specified by user")

This function allows to import SIPs in CSV format from Julia using the SIP 2.0 Standard (Stochastic Information Packets, Savage[2009]) to build unified simulation models.
Based on the CSV package, the function cleans up meta data and cleans out redundant columns. You can also set a *source string* that will be included in your SIP DataFrame as a dimension.
"""
function importsip(FileName, source="Not specified by user")
    open(FileName) do file
        global linecounter = 1
            for ln in eachline(file)
            linecounter +=1
            end

            if linecounter < 1000
                rounding_digit = -2
            else rounding_digit = -3
            end

        Trial_Count = round(linecounter, digits=rounding_digit)
        global start_row = Int64(linecounter - Trial_Count)
        global header_row = start_row-1
        end
    println("Header Row: $header_row and Data_Row: $start_row")
    SLURP = CSV.read(FileName, header=header_row, datarow=start_row)
    #SLURP[!, :source] .= source # new syntax. now uses .functions to broadcast
    println("")
    print(describe(SLURP))
    println("\n")
return SLURP
end

"""
    genmeta(source_df, s_name="SLURPName")

SIPs using the SIP 2.0 Standard (Stochastic Information Packets, Savage[2009]) require that Meta Data be available and maintained. This function creates a file template to add MetaData to your SIP Library. This is a seperate file must accompany the DataFrame to generate the SIP Library correctly. If ommited, the file will export the SIP Library without any metadata.
"""
function genmeta(source_df, s_name="SLURPName")
    #File Operations
    csvfilename = string("$s_name","_meta.csv")
    #csvfilename = "$FileName.csv"
    csvfile = open(csvfilename, "w") #Create CSV csvfile
    ColHeadings = SIPNames(source_df)

    #CSV Body
    write(csvfile,"This file is used to store SIP level MetaData and dimensions for SIP Library : $s_name,,
    ,,,
    ,,SIP Attribute, $ColHeadings,,
    ,,,\n")

#Close File
    close(csvfile)
    CSV.write("$csvfilename.csv")
end


text:
module TestCreate

using Test
using IonBase
using IonBase.Templates
using Configurations

include(joinpath(pkgdir(IonBase), "test", "utils.jl"))

with_test_ion() do
    IonBase.create(joinpath(test_dir, "Basic"); user="abc", force=true)
    IonBase.create(joinpath(test_dir, "Academic"); user= "abc", template = "academic", force=true)
    IonBase.create(joinpath(test_dir, "Package"); user= "abc", template = "package", force=true)
    IonBase.create(joinpath(test_dir, "Comonicon"); user= "abc", template = "comonicon", force=true)
    IonBase.create(joinpath(test_dir, "ComoniconSysImg"); user= "abc", template = "comonicon-sysimg", force=true)
end

end


text:
function irfvar(y,p,H)

# Luca Brugnolini 2014
# Based on Kilian code
# Point estimates and stdev of structural IRFs identified by recursive structure
# Input : data y, lag order p, horizon for impulse response H
# Output: IRF, Std error of IRF, A0inv for structural identification,
#         slope parameter estimates B, cov matrix of estimated residual SIG
#         estimated residual U, reduced form IRF

(T,K) = size(y);

# OLS
Y = (y[p+1:T,:])';                           # LHS variable for OLS
X = ones(1,T-p);                             # constant term

for j=1:p
    X = [X; y[p+1-j:T-j,:]'];                # RHS variable for OLS

end

B   = (Y*X')/(X*X');                         # OLS estimates **
A   = B[:,2:end];                            # OLS estimates w/o intercept
U   = Y-B*X;                                 # estimated residual
SIG = (Y-B*X)*(Y-B*X)'/(T-p-p*K-1);          # estimated sigma(U)
SIGb = kron(inv(X*X'/(T-p)),SIG);            # varcov of vec(B):(K^2*p+K)^2
SIGa = SIGb[K+1:end,K+1:end];                # varcov of vec(A):(K^2*p)^2

# IRF and stdev(IRF)
#  Calculation of stdev follows Lutkepohl(2005) p.111,93
A = [A; eye(K*(p-1)) zeros(K*(p-1),K)];         # Companion form
A0inv = chol(SIG)';                              # recursive str assumption
rfIRF = zeros(K^2,H+1);                          # reduced form IRF
IRF   = zeros(K^2,H+1);                          # structural IRF
STD   = zeros(K^2,H+1);                          # stdev of structural IRF
COV2   = zeros(K^2,H+1);                         # COV of structural IRF due to SIG
J = [eye(K) zeros(K,K*(p-1))];                   # J matrix
L = elimat(K);                                   # elimination matrix
Kk = commutation(K,K);                           # commutation matrix
Hk = L'/(L*(eye(K^2)+Kk)*kron(A0inv,eye(K))*L');  # H matrix **
D = duplication(K);                              # duplication matrix
Dplus = (D'*D)\D';
SIGsig = 2*Dplus*kron(SIG,SIG)*Dplus';           # sigma(vech(sigma))
Cbar0 = kron(eye(K),J*eye(K*p)*J')*Hk;           # Cbar
rfIRF[:,1] = vec((J*eye(K*p)*J')');
IRF[:,1] = vec((J*eye(K*p)*J'*A0inv)');
STD[:,1] = vec((reshape(diag(real(sqrt(complex(Cbar0*SIGsig*Cbar0'/(T-p))))),K,K))')
COV2[:,1] = vec((reshape(diag((Cbar0*SIGsig*Cbar0'/(T-p))),K,K))');

for h=1:H
    Gi = zeros(K^2,K^2*p);

    for m=0:(h-1)
        Gi += kron(J*(A')^(h-1-m),J*(A^m)*J');
    end;

    C = kron(A0inv',eye(K))*Gi;
    Cbar = kron(eye(K),J*A^h*J')*Hk;

    rfIRF[:,h+1] = vec((J*A^h*J')');
    IRF[:,h+1] = vec((J*A^h*J'*A0inv)');
    STD[:,h+1] = vec((reshape(diag(real(sqrt(complex(C*SIGa*C'+Cbar*SIGsig*Cbar')/(T-p)))),K,K))')
    COV2[:,h+1] = vec((reshape(diag(((Cbar*SIGsig*Cbar')/(T-p))),K,K))');

end

    return IRF,STD,A0inv,B,SIG,U,rfIRF,COV2;

end


text:
using MLStyle

function parse_line(l)
    mv, i = split(l)
    mv, parse(Int, i)
end

const data = parse_line.(readlines())

function solve_1()
    pos, depth = 0, 0
    for (mv, i) in data
        @match (mv, i) begin
            ("forward", i) => (pos += i)
            ("up", i) => (depth -= i)
            ("down", i) => (depth += i)
        end
    end
    pos * depth
end

function solve_2()
    pos, depth, aim = 0, 0, 0
    for (mv, i) in data
        @match (mv, i) begin
            ("forward", i) => (pos += i; depth += aim * i)
            ("up", i) => (aim -= i)
            ("down", i) => (aim += i)
        end
    end
    pos * depth
end

solve_1() |> println
solve_2() |> println

max_stars_repo_name: pLOPeGG/JuliaOfCode

max_issues_repo_name: pLOPeGG/JuliaOfCode

max_forks_repo_name: pLOPeGG/JuliaOfCode


text:
function compute_log_joint(observations, partition, tau0, tau1, sigma, theta)
  n = length(observations)
  k = length(partition)
  prob = k*log(sigma) + lgamma(theta) + lgamma(theta/sigma + k) - lgamma(theta/sigma) - lgamma(theta + n)
  for cluster in partition
    prob += lgamma(length(cluster) - sigma) - lgamma(1 - sigma)
    prob += compute_log_conditonal_observations(observations, cluster, tau0, tau1)
  end
  prob
end

function compute_log_conditonal_observations(observations, cluster, tau0, tau1)
  nl = length(cluster)
  prob = (nl/2)*log(tau1) - (nl/2)*log(2*pi) + 0.5*log(tau0) + 0.5*log(tau0+nl)
  prob += -tau1/2*(sum(observations)) + 0.5*(tau0*mu_0+tau1*sum(observations[cluster]))^2/(tau0+nl*tau1)
  prob
end

# Test of similarity between distributions
function correct_posterior(empirical_probs, data, partitions, τ0, τ1, σ, θ)
    true_log_probs = map(p -> compute_log_joint(data, p, τ0, τ1, σ, θ), partitions)
    true_probs = exp.(true_log_probs)
    true_probs /= sum(true_probs)

    empirical_probs /= sum(empirical_probs)

    # compare distribitions
    # L2
    L2 = sum((empirical_probs - true_probs).^2)

    # Discrepancy
    discr = maximum(abs.(empirical_probs - true_probs))
    return L2, discr
end


text:
@testset "ModelKit" begin
    @testset "System (no parameters)" begin
        @var x y z
        function F_test(x)
            let (x, y, z) = (x[1], x[2], x[3])
                [
                    (x^2 + y + z + 2)^2 - 3,
                    4 * x^2 * z^2 * y + 4z - 6x * y * z^7,
                    (-z) / (x + y),
                ]
            end
        end

        f = System(F_test([x, y, z]), [x, y, z])

        @test sprint(show, f) == """
        System
         variables: x, y, z

         (x ^ 2 + y + z + 2) ^ 2 - 3
         (4 * x ^ 2 * z ^ 2 * y + 4z) - 6 * x * y * z ^ 7
         -z / (x + y)"""

        F = compile(f)

        @test size(F) == (3, 3)
        @test length(F) == 3
        @test size(F, 1) == 3
        @test size(F, 2) == 3
        u = [0.0, 0.0, 0]
        U = zeros(3, 3)
        v = [0.4814, -0.433, -0.82709]

        fdm = FD.central_fdm(5, 1)
        û = F_test(v)
        @test evaluate(F, v) ≈ û
        @test F(v) ≈ û
        @test ModelKit.interpreted(F)(v) ≈ û
        @test evaluate(ModelKit.interpreted(F), v) ≈ û
        @test evaluate!(u, F, v) ≈ û
        @test u ≈ û
        Û = FD.jacobian(fdm, F_test, v)[1]
        @test jacobian(F, v) ≈ Û
        @test jacobian(ModelKit.interpreted(F), v) ≈ Û
        @test jacobian!(U, F, v) ≈ Û
        @test U ≈ Û rtol = 1e-8
        u .= 0
        U .= 0
        @test evaluate_and_jacobian!(u, U, F, v) === nothing
        @test u ≈ û
        @test U ≈ Û
    end

    @testset "System (parameters)" begin
        @var x y z a b c
        function F_test(x, p)
            let (x, y, z, b, c, a) = (x[1], x[2], x[3], p[1], p[2], p[3])
                [
                    (x^2 + y + z * a^2 + 2)^2 - 3,
                    4 * x^2 * z^2 * y + b + 4z - 6x * y * z^2,
                    a * c * (-z) / (x + y),
                ]
            end
        end

        f = System(
            [
                (x^2 + y + z * a^2 + 2)^2 - 3,
                4 * x^2 * z^2 * y + b + 4z - 6x * y * z^2,
                a * c * (-z) / (x + y),
            ],
            [x, y, z],
            [b, c, a],
        )
        @test sprint(show, f) == """
        System
         variables: x, y, z
         parameters: b, c, a

         (x ^ 2 + y + z * a ^ 2 + 2) ^ 2 - 3
         (4 * x ^ 2 * z ^ 2 * y + b + 4z) - 6 * x * y * z ^ 2
         (a * c * -z) / (x + y)"""

        F = compile(f)
        @test startswith(sprint(show, F), "Compiled{TSystem{3,3,3,#")

        @test size(F) == (3, 3)
        @test length(F) == 3
        @test size(F, 1) == 3
        @test size(F, 2) == 3
        u = [0.0, 0.0, 0]
        U = zeros(3, 3)
        v = [0.4814, -0.433, -0.82709]
        p = [0.1312, 2.1232, 1.2411]

        fdm = FD.central_fdm(5, 1)
        û = F_test(v, p)
        @test evaluate(F, v, p) ≈ û
        @test F(v, p) ≈ û
        @test evaluate(ModelKit.interpreted(F), v, p) ≈ û
        @test ModelKit.interpreted(F)(v, p) ≈ û
        @test evaluate!(u, F, v, p) ≈ û
        @test u ≈ û
        Û = FD.jacobian(fdm, v -> F_test(v, p), v)[1]
        @test jacobian(F, v, p) ≈ Û
        @test jacobian(ModelKit.interpreted(F), v, p) ≈ Û
        @test jacobian!(U, F, v, p) ≈ Û
        @test U ≈ Û rtol = 1e-8
        u .= 0
        U .= 0
        @test evaluate_and_jacobian!(u, U, F, v, p) === nothing
        @test u ≈ û
        @test U ≈ Û
        v, V = evaluate_and_jacobian(F, v, p)
        @test v ≈ û
        @test V ≈ Û
    end

    @testset "Homotopy" begin
        @var x y z t

        function H_test(x, s)
            let (x, y, z, t) = (x[1], x[2], x[3], s)
                [x^2 + y + z + 2t, 4 * x^2 * z^2 * y + 4z - 6x * y * z^2]
            end
        end
        h = Homotopy(H_test([x, y, z], t), [x, y, z], t)

        @test sprint(show, h) == """
        Homotopy in t
         variables: x, y, z

         x ^ 2 + y + z + 2t
         (4 * x ^ 2 * z ^ 2 * y + 4z) - 6 * x * y * z ^ 2"""

        @test h == Homotopy(H_test([x, y, z], t), [x, y, z], t)

        H = compile(h)
        @test typeof(H) == typeof(compile(h))
        @test size(H) == (2, 3)
        @test size(H, 1) == 2
        @test size(H, 2) == 3
        @test length(H) == 2

        u = [0.0, 0.0]
        U = zeros(2, 3)
        v = [0.4814, -0.433, -0.82709]
        s = 0.359

        fdm = FD.central_fdm(5, 1)
        û = H_test(v, s)
        @test evaluate(H, v, s) ≈ û
        @test H(v, s) ≈ û
        @test evaluate(ModelKit.interpreted(H), v, s) ≈ û
        @test ModelKit.interpreted(H)(v, s) ≈ û
        @test evaluate!(u, H, v, s) ≈ û
        @test u ≈ û
        Û = FD.jacobian(fdm, v -> H_test(v, s), v)[1]
        @test jacobian(H, v, s) ≈ Û
        @test jacobian(ModelKit.interpreted(H), v, s) ≈ Û
        @test jacobian!(U, H, v, s) ≈ Û
        @test U ≈ Û rtol = 1e-8
        u .= 0
        U .= 0
        @test evaluate_and_jacobian!(u, U, H, v, s) === nothing
        @test u ≈ û
        @test U ≈ Û
        #
        # # û = FD.grad(fdm, s -> H_test(v, s), s)[1]
        # @test dt(H, v, s) ≈ û
        # @test dt(ModelKit.interpreted(H), v, s) ≈ û
        # @test dt!(u, H, v, s) ≈ û
        # u .= 0
        # U .= 0
        # @test jacobian_and_dt!(U, u, H, v, s) === nothing
        # @test u ≈ û
        # @test U ≈ Û
    end

    @testset "Homotopy (parameters)" begin
        @var x y z t a b c

        function H_test(x, s, p)
            let (x, y, z, t, a, b, c) = (x[1], x[2], x[3], s, p[1], p[2], p[3])
                [
                    x^2 + y + z + a^2 - b + 2t,
                    4 * x^2 * z^2 * a * b * y + c * 4z - 6x * y * z^2,
                ]
            end
        end

        h = Homotopy(H_test([x, y, z], t, [a, b, c]), [x, y, z], t, [a, b, c])
        @test sprint(show, h) == """
        Homotopy in t
         variables: x, y, z
         parameters: a, b, c

         ((x ^ 2 + y + z + a ^ 2) - b) + 2t
         (4 * x ^ 2 * z ^ 2 * a * b * y + c * 4 * z) - 6 * x * y * z ^ 2"""

        H = compile(h)
        @test startswith(sprint(show, H), "Compiled{THomotopy{2,3,3,#")

        @test size(H) == (2, 3)
        @test size(H, 1) == 2
        @test size(H, 2) == 3
        @test length(H) == 2

        u = [0.0, 0.0]
        U = zeros(2, 3)
        v = [0.4814, -0.433, -0.82709]
        s = 0.359
        p = [-2.123, 0.412, 1.21321]

        fdm = FD.central_fdm(5, 1)
        û = H_test(v, s, p)
        @test evaluate(H, v, s, p) ≈ û
        @test H(v, s, p) ≈ û
        @test evaluate(ModelKit.interpreted(H), v, s, p) ≈ û
        @test ModelKit.interpreted(H)(v, s, p) ≈ û
        @test evaluate!(u, H, v, s, p) ≈ û
        @test u ≈ û
        Û = FD.jacobian(fdm, v -> H_test(v, s, p), v)[1]
        @test jacobian(H, v, s, p) ≈ Û
        @test jacobian(ModelKit.interpreted(H), v, s, p) ≈ Û
        @test jacobian!(U, H, v, s, p) ≈ Û
        @test U ≈ Û rtol = 1e-8
        u .= 0
        U .= 0
        @test evaluate_and_jacobian!(u, U, H, v, s, p) === nothing
        @test u ≈ û
        @test U ≈ Û

        # û = FD.grad(fdm, s -> H_test(v, s, p), s)[1]
        # @test dt(H, v, s, p) ≈ û
        # @test dt(ModelKit.interpreted(H), v, s, p) ≈ û
        # @test dt!(u, H, v, s, p) ≈ û
        # u .= 0
        # U .= 0
        # @test jacobian_and_dt!(U, u, H, v, s, p) === nothing
        # @test u ≈ û
        # @test U ≈ Û
        #
        # v, V = jacobian_and_dt(H, v, s, p)
        # @test v ≈ û
        # @test V ≈ Û
    end

    @testset "Modeling" begin
        @testset "Bottleneck" begin
            @var x y z
            f = [
                (0.3 * x^2 + 0.5z + 0.3x + 1.2 * y^2 - 1.1)^2 +
                    (0.7 * (y - 0.5x)^2 + y + 1.2 * z^2 - 1)^2 - 0.3,
            ]

            I = let
                x = variables(f)
                n, m = length(x), length(f)
                @unique_var y[1:n] v[1:m] w[1:m]
                J = [differentiate(fᵢ, xᵢ) for fᵢ in f, xᵢ in x]
                f′ = [subs(fᵢ, x => y) for fᵢ in f]
                J′ = [subs(gᵢ, x => y) for gᵢ in J]
                Nx = (x - y) - J' * v
                Ny = (x - y) - J′' * w
                System([f; f′; Nx; Ny], [x; y; v; w])
            end
            @test I isa System
            @test size(I) == (8, 8)
        end

        @testset "Steiner" begin
            @var x[1:2] a[1:5] c[1:6] y[1:2, 1:5]

            #tangential conics
            f = sum([a; 1] .* monomials(x, 2))
            ∇ = differentiate(f, x)
            #5 conics
            g = sum(c .* monomials(x, 2))
            ∇_2 = differentiate(g, x)
            #the general system
            #f_a_0 is tangent to g_b₀ at x₀
            function Incidence(f, a₀, g, b₀, x₀)
                fᵢ = f(x => x₀, a => a₀)
                ∇ᵢ = [∇ᵢ(x => x₀, a => a₀) for ∇ᵢ in ∇]
                Cᵢ = g(x => x₀, c => b₀)
                ∇_Cᵢ = [∇ⱼ(x => x₀, c => b₀) for ∇ⱼ in ∇_2]

                [fᵢ; Cᵢ; det([∇ᵢ ∇_Cᵢ])]
            end
            @var v[1:6, 1:5]
            I = vcat(map(i -> Incidence(f, a, g, v[:, i], y[:, i]), 1:5)...)
            F = System(I, [a; vec(y)], vec(v))
            @test size(F) == (15, 15)
        end

        @testset "Reach plane curve" begin
            @var x y
            f = (x^3 - x * y^2 + y + 1)^2 * (x^2 + y^2 - 1) + y^2 - 5
            ∇ = differentiate(f, [x; y]) # the gradient
            H = differentiate(∇, [x; y]) # the Hessian

            g = ∇ ⋅ ∇
            v = [-∇[2]; ∇[1]]
            h = v' * H * v
            dg = differentiate(g, [x; y])
            dh = differentiate(h, [x; y])

            ∇σ = g .* dh - ((3 / 2) * h) .* dg

            F = System([v ⋅ ∇σ; f], [x, y])
            @test size(F) == (2, 2)
        end
    end
end


text:
# TODO: add documentation and example! (MPC)

# terminal set around origin => how to shift and still satisfy constraints?
# just shift \bbX ???? and apply algo?
# tube-based MPC => X,U are X_bar, U_bar!!!
function terminal_set(A, B, 𝕏::LazySet, 𝕌::LazySet, K; max_iter::Integer=100)
    return maximum_invariant_set(A, B, 𝕏, 𝕌, K, max_iter)
end

# How do we compute P?
function terminal_set(center, P, α, 𝕏::LazySet)
    # TODO: extend ELlipsoid to have biggest ELlipsoid in X
    E =  Ellipsoid(center, inv(P)/α)
    if E ⊆ 𝕏
        return E
    end
    error("E not contained in 𝕏")
end

function feasible_set(A, B, 𝕏::LazySet, 𝕌::LazySet, 𝕏f::LazySet, N::Integer;
                       preset=_preset_vrep)
    𝕏_feas = preset(A, B, 𝕌, 𝕏f) ∩ᶜ 𝕏
    for i=2:N
        𝕏_feas = preset(A, B, 𝕌, 𝕏_feas) ∩ᶜ 𝕏
    end
    return 𝕏_feas
end

# TODO: Add once the paper is published and write Tests!
"""
    tightened_constraints(A, B, W, X, U, K)

Source Code not yet published due to impending paper.
"""
function tightened_constraints() end
function tightened_input_constraint() end
function tightened_state_constraint() end

max_stars_repo_name: mforets/InvariantSets

max_issues_repo_name: mforets/InvariantSets


text:
# isstruct(::Type{T}) where T = !isbitstype(T) && !(T <: AbstractArray)
# isstruct(obj) = !isbits(obj) && !isa(obj, AbstractArray)
"Check if an object is of a struct type, i.e. not a number or array"
isstruct(::Type{T}) where T = !isempty(fieldnames(T))
isstruct(obj) = isstruct(typeof(obj))

if !isdefined(@__MODULE__, :__EXPRESSION_HASHES__)
    __EXPRESSION_HASHES__ = Set{AbstractString}()
end

"""
If loaded twice without changes, evaluate expression only for the first time.
"""
macro runonce(expr)
    h = string(expr)
    return esc(quote
        if !in($h, __EXPRESSION_HASHES__)
            push!(__EXPRESSION_HASHES__, $h)
            $expr
        end
    end)
end


text:
mutable struct MiddlewareStore{T} <: AbstractStore{T}
    store::Store{T}
    middlewares::Vector
end
MiddlewareStore(store::Store) = MiddlewareStore(store, [])

function identity_middleware(dispatch)
    new_dispatch = (s::Store, a::AbstractAction) -> begin
        # prologue
        result = dispatch(s, a)
        # epilogue
        return result
    end
    return new_dispatch
end

function dispatch!(store::MiddlewareStore, action::AbstractAction)
    isempty(store.middlewares) && dispatch!(store.store, action)
    dispatch = identity_middleware(dispatch!)
    for middleware in store.middlewares
        dispatch = middleware(dispatch)
    end
    return dispatch(store.store, action)
end

get_state(store::MiddlewareStore) = get_state(store.store)
subscribe!(store::MiddlewareStore, listener::Base.Callable) = subscribe!(store.store, listener)

function create_store(reducer::Base.Callable, preloaded_state, middlewares)
    store = create_store(reducer, preloaded_state)
    return MiddlewareStore(store, middlewares)
end


text:
using packageTrial
using Test

@testset "myExample.jl" begin
    # 2x+3y
    @test my_f(2,1) == 7
    @test my_f(2,3) == 13
end

#@testset "Derivative" begin
#    @test ** == **
#end


text:
include("../benchmark.jl")
using OnePhase, advanced_timer

function run_cutest_problems_on_solver(problems::Array{String,1}, test_name::String, solver)
    summary = Dict{String, problem_summary2}()

    if_mkdir("../results/$test_name")
    if_mkdir("../results/$test_name/log")

    for problem_name in problems
          println("RUNNING $problem_name")

          ORG_STDOUT = STDOUT
          file = open("../results/$(test_name)/log/$(problem_name).txt", "w")
          redirect_stdout(file)
          nlp_raw = CUTEstModel(problem_name)
          summary[problem_name] = problem_summary2()
          tic()

          try
            t = 0
            function intermediate_ipopt(alg_mod::Int,iter_count::Int,obj_value::Float64,inf_pr::Float64, inf_du::Float64,  mu::Float64, d_norm::Float64, regularization_size::Float64, alpha_du::Float64, alpha_pr::Float64, ls_trials::Int)
              t = iter_count
              return true  # Keep going
            end

            mp = NLPModels.NLPtoMPB(nlp_raw, solver)
            setIntermediateCallback(mp.inner, intermediate_ipopt)
            MathProgBase.optimize!(mp)
            #x = MathProgBase.getsolution(mp)
            #solver = MathProgBase.getrawsolver(mp)

            status = MathProgBase.status(mp)
            summary[problem_name].status = status;
            summary[problem_name].it_count = t;

            x = MathProgBase.getsolution(mp)
            set_cutest_info_ipopt!(summary[problem_name], mp.inner, nlp_raw, x)
            println("here4")
          catch(e)
            println("Uncaught error in algorithm!!!")
            @show e;
            summary[problem_name].status = :ERR
            summary[problem_name].it_count = -1;
          end

          summary[problem_name].total_time = toc();

          redirect_stdout(ORG_STDOUT)
          finalize(nlp_raw)
          close(file)

          println("it count = ", summary[problem_name].it_count)
          println("status = ", summary[problem_name].status)

          save("../results/$(test_name)/summary.jld","summary",summary)

          summary_file = open("../results/$(test_name)/summary.txt", "w")
          write_summary(summary_file, summary)
          close(summary_file)
    end
end

function run_cutest_problems_using_our_solver(problems::Array{String,1}, test_name::String, par::Class_parameters)

    if_mkdir("../results/$test_name")
    if_mkdir("../results/$test_name/log")
    if_mkdir("../results/$test_name/jld")

    if isfile("../results/$(test_name)/summary.jld")
      summary = load("../results/$(test_name)/summary.jld","summary")
    else
      summary = Dict{String, problem_summary2}()
    end

    already_solved_problems = keys(summary)

    par_file = open("../results/$(test_name)/par.txt", "w")
    write_pars(par_file, par)
    close(par_file)

    master_timer = class_advanced_timer()

    for problem_name in problems
      if problem_name in already_solved_problems
          println("$problem_name already solved")
      else
          println("RUNNING $problem_name")
          ORG_STDOUT = STDOUT
          file = open("../results/$(test_name)/log/$(problem_name).txt", "w")
          redirect_stdout(file)
          summary[problem_name] = problem_summary2()

          nlp_raw = false
          nlp = false

          try
              try
                nlp_raw = CUTEstModel(problem_name)
                #nlp = Class_CUTEst(nlp_raw)
              catch(e)
                println("$problem_name failed to load")
                throw(e)
              end

              tic()

              #=timer = class_advanced_timer()
              start_advanced_timer(timer)
              #include("include.jl")
              #intial_it = initial_point_satisfy_bounds(nlp, par)
              start_advanced_timer(timer, "INIT")
              intial_it = OnePhase.init(nlp, par, timer)
              pause_advanced_timer(timer, "INIT")

              #intial_it = initial_point_generic(nlp, par, nlp_raw.meta.x0)

              @assert(OnePhase.is_feasible(intial_it, par.ls.comp_feas))
              iter, status, history, t, err = OnePhase.one_phase_IPM(intial_it, par, timer);

              pause_advanced_timer(timer)

              print_timer_stats(timer)
              =#
              iter, status, history, t, err, timer = OnePhase.one_phase_solve(nlp_raw, par)

              master_timer = merge_timers(timer, master_timer)

              save("../results/$(test_name)/jld/$(problem_name).jld","history",history, "timer", timer)

              set_info_me!(summary[problem_name], history, status)
              #.it_count = t;

              summary[problem_name].total_time = toc();
          catch(e)
              println("Uncaught error in algorithm!!!")
              @show e;

              if isa(e, Eval_NaN_error)
                summary[problem_name].status = :NaN_ERR
                summary[problem_name].it_count = -1;
              else
                summary[problem_name].status = :ERR
                summary[problem_name].it_count = -1;
              end
          end

          redirect_stdout(ORG_STDOUT)
          finalize(nlp_raw)
          close(file)

          println("it count = ", summary[problem_name].it_count)
          println("status = ", summary[problem_name].status)

          save("../results/$(test_name)/summary.jld","summary",summary, "pars", par)

          summary_file = open("../results/$(test_name)/summary.txt", "w")
          write_summary(summary_file, summary)
          close(summary_file)

          timer_file = open("../results/$(test_name)/timer.txt", "w")
          print_timer_stats(timer_file, master_timer)
          close(timer_file)
      end
    end
end

#function run_cutest_problems_using_IPOPT(problems::Array{String,1}, test_name::String)
#
#end

#run_cutest_problems(["DISCS"], "test")

function test_problems(problem_list::Array{String,1}, start::Int64)
    for i = start:length(problem_list)
      problem_name = problem_list[i]
      println(problem_name)
      nlp_raw = CUTEstModel(problem_name)
      println(i)
      finalize(nlp_raw)
    end
end

function filter_string_array(problem_list::Array{String,1}, remove_list::Array{String,1})
    new_problem_list = Array{String,1}()
    for problem_name in problem_list
      if !(problem_name in remove_list)
        push!(new_problem_list, problem_name)
      end
    end
    return new_problem_list
end

function get_problem_list(filter::Function)
  #problem_list = CUTEst.select(max_var=100, max_var=1000, min_con=100, max_con=3000)
  problem_list = CUTEst.select(custom_filter=filter);
  problem_list = convert(Array{String,1},problem_list);

  remove_list = ["MISRA1D","OSBORNE1","LANCZOS2","MEYER3NE","ROSZMAN1","INTEQNE","MGH10",
  "ECKERLE4","RAT43","MISRA1A","KOWOSBNE","MOREBVNE",
  "JENSMPNE","MISRA1C","VARDIMNE","SANTA","BARDNE","MGH10S",
  "THURBER","MGH17S","BOX3NE","PENLT1NE","GAUSS2","GULFNE","BA-L1SP","ENSO","CHWIRUT1", "CHWIRUT2","NELSON","HAHN1","KIRBY2","MUONSINE","OSBORNE2","BENNETT5","BA-L1","GAUSS3",
  "GAUSS1","LANCZOS3","BIGGS6NE","LANCZOS1","VANDANIUMS","MISRA1B","MGH09","MGH17","WATSONNE",
  "DMN37142LS","VESUVIOU","PENLT2NE", "DMN37143LS","INTEQNELS","BROYDN3DLS","DMN15332LS",
  "DMN15102LS","BROYDNBDLS","VESUVIA","BA-L1SPLS","BA-L1SPLS","SANTALS","VESUVIO","DMN15333LS",
  "ARGTRIGLS","RAT42","SSINE","LSC1","LSC2","BOXBOD","POWELLSE","FREURONE","DANWOOD","HELIXNE",
  "DMN37142","DMN37143","DMN15102","DMN15333","DMN15332","DMN15103",
  "VANDERM4"]; # manual remove because of NaNs
  problem_list = filter_string_array(problem_list, remove_list);
  #test_problems(problem_list, 1)
  #CUTEstModel("GAUSS2")

  return problem_list
end

function get_problem_list(min_size::Int64, max_size::Int64)
  function filter_cutest(problem)
      min_size_ok = problem["constraints"]["number"] >= min_size && problem["variables"]["number"] >= min_size
      #max_size_ok = problem["constraints"]["number"] <= max_size && problem["variables"]["number"] <= max_size
      max_size_ok = problem["constraints"]["number"] + problem["variables"]["number"] <= max_size

      correct_size = min_size_ok && max_size_ok

      regular = problem["derivative_order"] >= 2 && problem["regular"] == true
      if correct_size && regular
          return true
      else
        return false
      end
  end

  return get_problem_list(filter_cutest)
end

function get_problem_list(min_ncon::Int64, max_ncon::Int64, min_nvar::Int64, max_nvar::Int64)
  function filter_cutest(problem)
      min_size_ok = problem["constraints"]["number"] >= min_ncon && problem["variables"]["number"] >= min_nvar
      max_size_ok = problem["constraints"]["number"] <= max_ncon && problem["variables"]["number"] <= max_nvar
      #max_size_ok = problem["constraints"]["number"] + problem["variables"]["number"] <= max_size

      correct_size = min_size_ok && max_size_ok

      regular = problem["derivative_order"] >= 2 && problem["regular"] == true
      if correct_size && regular
          @show problem["constraints"]["number"], problem["variables"]["number"]
          return true
      else
        return false
      end
  end

  return get_problem_list(filter_cutest)
end

function select_problem_with_sparse_rows(problem_list::Array{String,1},max_density::Int64)
    sparse_names = Array{String,1}()
    dense_names = Array{String,1}()
    for problem_name in problem_list
        nlp_raw = CUTEstModel(problem_name)

        J = jac(nlp_raw, nlp_raw.meta.x0)
        if densest_row(J) <= max_density
          push!(sparse_names, problem_name)
        else
          push!(dense_names, problem_name)
        end

        finalize(nlp_raw)
    end

    return sparse_names, dense_names
end

function default_list(sparsify=false)
    problem_list = get_problem_list(10,10000)
    # only run problems with max row density 1000.
    if sparsify
      sparse_names, dense_names = select_problem_with_sparse_rows(problem_list, 1000)
      problem_list = sparse_names
      println("these problems are included ...")
      @show sparse_names
      println("these problems are excluded ...")
      @show dense_names
    end

    return problem_list
end


text:
export MultiThreadEnv

using Base.Threads: @spawn

"""
    MultiThreadEnv(envs::Vector{<:AbstractEnv})

Wrap multiple instances of the same environment type into one environment.
Each environment will run in parallel by leveraging `Threads.@spawn`.
So remember to set the environment variable `JULIA_NUM_THREADS`!
"""
struct MultiThreadEnv{E,S,R,AS,SS,L} <: AbstractEnv
    envs::Vector{E}
    states::S
    rewards::R
    terminals::BitArray{1}
    action_space::AS
    state_space::SS
    legal_action_space_mask::L
end

function Base.show(io::IO, t::MIME"text/markdown", env::MultiThreadEnv)
    print(io, "MultiThreadEnv($(length(env)) x $(nameof(env[1])))")
end

"""
    MultiThreadEnv(f, n::Int)

`f` is a lambda function which creates an `AbstractEnv` by calling `f()`.
"""
MultiThreadEnv(f, n::Int) = MultiThreadEnv([f() for _ in 1:n])

function MultiThreadEnv(envs::Vector{<:AbstractEnv})
    n = length(envs)
    S = state_space(envs[1])
    s = state(envs[1])
    if S isa Space
        S_batch = similar(S, size(S)..., n)
        s_batch = similar(s, size(s)..., n)
        for j in 1:n
            Sₙ = state_space(envs[j])
            sₙ = state(envs[j])
            for i in CartesianIndices(size(S))
                S_batch[i, j] = Sₙ[i]
                s_batch[i, j] = sₙ[i]
            end
        end
    else
        S_batch = Space(state_space.(envs))
        s_batch = state.(envs)
    end

    A = action_space(envs[1])
    if A isa Space
        A_batch = similar(A, size(A)..., n)
        for j in 1:n
            Aⱼ = action_space(envs[j])
            for i in CartesianIndices(size(A))
                A_batch[i, j] = Aⱼ[i]
            end
        end
    else
        A_batch = Space(action_space.(envs))
    end

    r_batch = reward.(envs)
    t_batch = is_terminated.(envs)
    if ActionStyle(envs[1]) === FULL_ACTION_SET
        m_batch = BitArray(undef, size(A)..., n)
        for j in 1:n
            L = legal_action_space_mask(envs[j])
            for i in CartesianIndices(size(A))
                m_batch[i, j] = L[i]
            end
        end
    else
        m_batch = nothing
    end
    MultiThreadEnv(envs, s_batch, r_batch, t_batch, A_batch, S_batch, m_batch)
end

MacroTools.@forward MultiThreadEnv.envs Base.getindex, Base.length, Base.iterate

function (env::MultiThreadEnv)(actions)
    N = ndims(actions)
    @sync for i in 1:length(env)
        @spawn begin
            if N == 1
                env[i](actions[i])
            else
                env[i](selectdim(actions, N, i))
            end
        end
    end
end

function RLBase.reset!(env::MultiThreadEnv; is_force = false)
    if is_force
        for i in 1:length(env)
            reset!(env[i])
        end
    else
        @sync for i in 1:length(env)
            if is_terminated(env[i])
                @spawn begin
                    reset!(env[i])
                end
            end
        end
    end
end

const MULTI_THREAD_ENV_CACHE = IdDict{AbstractEnv,Dict{Symbol,Array}}()

function RLBase.state(env::MultiThreadEnv)
    N = ndims(env.states)
    @sync for i in 1:length(env)
        @spawn selectdim(env.states, N, i) .= state(env[i])
    end
    env.states
end

function RLBase.reward(env::MultiThreadEnv)
    env.rewards .= reward.(env.envs)
    env.rewards
end

function RLBase.is_terminated(env::MultiThreadEnv)
    env.terminals .= is_terminated.(env.envs)
    env.terminals
end

function RLBase.legal_action_space_mask(env::MultiThreadEnv)
    N = ndims(env.states)
    @sync for i in 1:length(env)
        @spawn selectdim(env.legal_action_space_mask, N, i) .=
            legal_action_space_mask(env[i])
    end
    env.legal_action_space_mask
end

RLBase.action_space(env::MultiThreadEnv) = env.action_space
RLBase.state_space(env::MultiThreadEnv) = env.state_space
RLBase.legal_action_space(env::MultiThreadEnv) = Space(legal_action_space.(env.envs))
# RLBase.current_player(env::MultiThreadEnv) = current_player.(env.envs)

for f in RLBase.ENV_API
    if endswith(String(f), "Style")
        @eval RLBase.$f(x::MultiThreadEnv) = $f(x[1])
    end
end

#####
# Patches
#####

(env::MultiThreadEnv)(action::EnrichedAction) = env(action.action)

function (π::QBasedPolicy)(env::MultiThreadEnv, ::MinimalActionSet, A)
    [A[i][a] for (i, a) in enumerate(π.explorer(π.learner(env)))]
end

function (π::QBasedPolicy)(env::MultiThreadEnv, ::FullActionSet, A)
    [
        A[i][a] for
        (i, a) in enumerate(π.explorer(π.learner(env), legal_action_space_mask(env)))
    ]
end

function (π::QBasedPolicy)(
    env::MultiThreadEnv,
    ::MinimalActionSet,
    ::Space{<:Vector{<:Base.OneTo{<:Integer}}},
)
    π.explorer(π.learner(env))
end

function (π::QBasedPolicy)(
    env::MultiThreadEnv,
    ::FullActionSet,
    ::Space{<:Vector{<:Base.OneTo{<:Integer}}},
)
    π.explorer(π.learner(env), legal_action_space_mask(env))
end


text:
using LiveServer

export
    build_templates,
    serve_templates

"""
    newsite(topdir; template="basic", cd=true)

Generate a new folder (an error is thrown if it already exists) that contains
the skeleton of a website that can be processed by Franklin. The user can
specify a `template` out of the list of available templates.
If `topdir` is specified as `"."` then the current directory is used.

* `template="basic"`: name of the template to use,
* `cd=true`:          whether to change the current directory to the newly
                      created folder or not.
* `verbose=true`:     whether to display information or not.

### Example

```julia
newsite("MyNewWebsite", template="pure-sm")
```
"""
function newsite(topdir::String="TestWebsite";
                 template::String="basic", changedir::Bool=true,
                 verbose::Bool=true)

    template = lowercase(template)
    template ∈ LIST_OF_TEMPLATES ||
        throw(ArgumentError("Template $template doesn't exist."))

    # create the top-directory
    if topdir == "."
        topdir    = pwd()
        changedir = true
    else
        topdir = mkdir(topdir)
    end

    common_dir   = joinpath(TEMPL_PATH, "common")
    template_dir = joinpath(TEMPL_PATH, template)

    mergefolders(common_dir,   topdir)
    mergefolders(template_dir, topdir)

    # for sandbox, remove aux pages menu1,menu2,menu3
    if template in ("sandbox", "sandbox-extended")
        for i in 1:3
            rm(joinpath(topdir, "menu$i.md"))
        end
    end

    # For 'added' packages, Pkg.jl makes some files read-only, so here we
    # restore to 644 to guarantee that they are r/w
    for (root, _, files) ∈ walkdir(topdir)
        for file in files
            chmod(joinpath(root, file), 0o644)
        end
    end

    # check if the template has a pre-specified index.html, if so rename
    # index.md to avoid confusion
    if isfile(joinpath(topdir, "index.html"))
        mv(joinpath(topdir, "index.md"), joinpath(topdir, "franklin.md"))
    end

    # move to the directory if relevant
    changedir && cd(topdir)

    # display information as adequate
    verbose && begin
        print("✓ Website folder generated at \"$(topdir)\"")
        println(ifelse(changedir, " (now the current directory)." , "."))
        println("→ Use serve() from Franklin to see the website in your browser.")
    end

    return nothing
end

"""
    modify(template_name)

Convenience function to quickly instantiate a directory with a given template
once satisfied with the changes, they can be added to the repo by doing
`addtemplate(path)`.
"""
function modify(name)
    @assert name in LIST_OF_TEMPLATES "unknown template"
    newsite(name, template=name)
    return nothing
end

"""
    addtemplate(path)

Assuming you're working with FranklinTemplates on development mode (i.e.
`] dev FranklinTemplates` this function takes a path where you would have
worked on a template and incorporates it into FranklinTemplates, putting files
at the right place; this makes iterations and PRs easier.

Note: if you use this without FranklinTemplates on dev mode, bad things will
likely occur.
"""
function addtemplate(path::String)
    path = abspath(path)
    @assert isdir(path) "'$path' does not seem to be a valid directory."

    # check whether it looks like a sensible dir
    rd = readdir(path)
    if !all(p->p in rd, ("_css", "_layout", "_libs"))
        error("The directory $path does not contain a `_css`, `_layout` and " *
              "a `_libs` folder and so will be ignored.")
    end

    # Add the folder at the right place in FranklinTemplates
    #
    # does it exist?
    d, tname = splitdir(path)
    isempty(tname) && (tname = splitdir(d)[2]) # it ended with separator

    dest = joinpath(TEMPL_PATH, tname, "")
    common = joinpath(TEMPL_PATH, "common", "")

    if tname ∉ readdir(TEMPL_PATH)
        cp(path, dest)
    else
        # it exists already --> merge the folders
        mergefolders(path, dest)
    end

    # postprocess
    if isdir(joinpath(dest, "__site"))
        rm(joinpath(dest, "__site"), recursive=true, force=true)
    end
    # remove duplicates
    deduplicate(common, dest)
    return nothing
end

"""
    mergefolders(src, dst)

Internal function to looks at what's inside `src/` and put it in `dst/`. If
there are paths that match, the files are merged. It is assumed that files will
not clash, if they clash files in `dst` are replaced.
See also [`newsite`](@ref).
"""
function mergefolders(src, dst)
    for (root, _, files) ∈ walkdir(src)
        for file ∈ files
            newpath = replace(root, Regex("^$(escape_string(src))")=>"$dst")
            isdir(newpath) || mkpath(newpath)
            newpathfile = joinpath(newpath, file)
            cp(joinpath(root, file), newpathfile; force=true)
        end
    end
end

"""
    deduplicate(src, dest)

Consider the files provided by `src` and remove all identical files from `dest`
to avoid duplication.
"""
function deduplicate(src, dest)
    # Post processs by removing files which are already provided by src
    dest = joinpath(dest, "") # ensure it has a separator at the end
    for (root, _, files) in walkdir(dest)
        for file in files
            fpath = joinpath(root, file)
            rpath = replace(fpath, dest => "")
            cpath = joinpath(src, rpath)
            isfile(cpath) || continue
            if filecmp(fpath, cpath)
                rm(fpath)
            end
        end
    end
    # there may be empty folders (these would be ignored by GitHub anyway)
    emptydirs = []
    for (root, dirs, _) in walkdir(dest)
        for dir in dirs
            dpath = joinpath(root, dir)
            if isdir(dpath) && isempty(readdir(dpath))
                push!(emptydirs, dpath)
            end
        end
    end
    for dpath in emptydirs
        try; rm(dpath, recursive=true); catch; end
    end
    return nothing
end

"""
    filecmp(path1, path2)

Take 2 absolute paths and check if the files are different (return false if
different and true otherwise).
This code was suggested by Steven J. Johnson on discourse:
https://discourse.julialang.org/t/how-to-obtain-the-result-of-a-diff-between-2-files-in-a-loop/23784/4
"""
function filecmp(path1::AbstractString, path2::AbstractString)
    stat1, stat2 = stat(path1), stat(path2)
    if !(isfile(stat1) && isfile(stat2)) || filesize(stat1) != filesize(stat2)
        return false
    end
    stat1 == stat2 && return true # same file
    open(path1, "r") do file1
        open(path2, "r") do file2
            buf1 = Vector{UInt8}(undef, 32768)
            buf2 = similar(buf1)
            while !eof(file1) && !eof(file2)
                n1 = readbytes!(file1, buf1)
                n2 = readbytes!(file2, buf2)
                n1 != n2 && return false
                0 != Base._memcmp(buf1, buf2, n1) && return false
            end
            return eof(file1) == eof(file2)
        end
    end
end

"""
    build_templates()

Build all the templates.
"""
function build_templates()
    root = dirname(dirname(@__FILE__))
    make_path = joinpath(root, "docs", "make.jl")
    include(make_path)
    for (root, _, files) ∈ walkdir(joinpath(root, "docs", "build"))
        for file ∈ files
            endswith(file, ".html") || continue;
            path  = joinpath(root, file);
            htmls = read(path, String);
            if get(ENV, "CI", nothing) == "true"
                htmls = replace(htmls, "href=\"/" => "href=\"/FranklinTemplates.jl/");
                htmls = replace(htmls, "src=\"/" => "src=\"/FranklinTemplates.jl/");
            end
            write(path, htmls);
        end
    end
end

"""
    serve_templates()

Rebuild all templates upon detecting changes.
"""
function serve_templates()
    build_templates()

    root = dirname(dirname(@__FILE__))
    function custom_callback(fp::AbstractString)
        if contains(fp, "src/")
            build_templates()
            cd(joinpath(root, "docs", "build"))
            LiveServer.file_changed_callback(fp)
        end
    end

    sw = LiveServer.SimpleWatcher(custom_callback)

    # Source files to watch for changes.
    for (root, dirs, files) in walkdir(joinpath(root, "src"))
        for file in files
            path = joinpath(root, file)
            LiveServer.watch_file!(sw, path)
        end
    end

    cd(joinpath(root, "docs", "build"))
    LiveServer.serve(sw; verbose=true)
end


text:
module NURBS
	export bezier, dbezier
	export knot, knotu, knotc
	export basis, dbasis, dbasisu
	export bspline, bsplineu, dbspline, dbsplineu, matpbspl, nmatrix
	export rbasis, rbspline, rbsplinu
	export bezsurf
	export bsplsurf, bspsurfu, dbsurf
	export rbspsurf
	export raise23, raise45

	export bezsurfj, bezsurf1, bezsurf2, bezsurf3, runn, runner

	include("bezier.jl")
	include("knot.jl")
	include("basis.jl")
	include("bspline.jl")
	include("rbspline.jl")
	include("beziersurfaces.jl")
	include("bsurfaces.jl")
	include("rbsurfaces.jl")
	include("raise.jl")

	include("bezsurftest.jl")
end


text:
export
    mPulseAPIException,
    mPulseAPIAuthException,
    mPulseAPIRequestException,
    mPulseAPIResultFormatException,
    mPulseAPIBugException

import HttpCommon.Response

"""
Thrown when the REST API has a problem and returns something other than a 2xx response.

#### Fields
`msg::AbstractString`
:    The error message

`response::Response`
:    The response object from the REST API call.  You can inspect headers, data, cookies, redirects, and the initiating request.

`responseBody::AbstractString`
:    The body of the HTTP response from the server

"""
immutable mPulseAPIException <: Exception
    msg::AbstractString
    response::Response
    responseBody::AbstractString

    mPulseAPIException(msg::AbstractString, response::Response) = new(msg, response, join(map(Char, response.data)))
end

"""
Thrown when the token used to authenticate with the REST API is invalid or has expired

#### Fields
`msg::AbstractString`
:    This message is always set to "Error Authenticating with REST API"

`response::Response`
:    The response object from the REST API call.  You can inspect headers, data, cookies, redirects, and the initiating request.

`responseBody::AbstractString`
:    The body of the HTTP response from the server

"""
immutable mPulseAPIAuthException <: Exception
    msg::AbstractString
    response::Response
    responseBody::AbstractString

    mPulseAPIAuthException(response::Response) = new("Error Authenticating with REST API", response, join(map(Char, response.data)))
end

"""
Thrown when a request parameter is invalid

#### Fields
`msg::AbstractString`
:    The error message sent from the mPulse server

`code::AbstractString`
:    The error code sent from the mPulse server

`parameter::AbstractString`
:    The parameter that the mPulse server had a problem with

`value::AbstractString`
:    The value of the parameter that the mPulse server had a problem with

`response::Response`
:    The response object from the REST API call.  You can inspect headers, data, cookies, redirects, and the initiating request.

"""
immutable mPulseAPIRequestException <: Exception
    msg::AbstractString
    code::AbstractString
    parameter::AbstractString
    value::AbstractString
    response::Union{Response, Void}
end

"""
Thrown when the result returned by an API call was not in the expected format

#### Fields
`msg::AbstractString`
:    The error message

`data::Any`
:    The actual data returned

"""
immutable mPulseAPIResultFormatException <: Exception
    msg::AbstractString
    data::Any
end

"""
Thrown when the REST API has an internal server error and returns a `500 Internal Server Error`

#### Fields
`msg::AbstractString`
:    The string "Internal Server Error, please report this. Timestamp: <current unix timestamp in seconds since the epoch>"

`response::Response`
:    The response object from the REST API call.  You can inspect headers, data, cookies, redirects, and the initiating request.

`responseBody::AbstractString`
:    The body of the HTTP response from the server

"""
immutable mPulseAPIBugException <: Exception
    msg::AbstractString
    response::Response
    responseBody::Response

    mPulseAPIBugException(resp::Response) = new("Internal Server Error, please report this. Timestamp: $(round(Int, datetime2unix(now())))", resp, join(map(Char, resp.data)))
end


text:
using JuMP, EAGO

                     m = Model()

                     EAGO.register_eago_operators!(m)

                     @variable(m, -1 <= x[i=1:3] <= 1)

                     @variable(m, -5.992552298143351 <= q <= 4.750344613921084)

                     add_NL_constraint(m, :((0.47930083070023954 + -0.2811056889538768*(-0.30287094621097754 + 0.8508304098923269*$(x[1]) + 0.544182872577263*$(x[2]) + -0.2945971031198069*$(x[3]))*(1 + erf((-0.30287094621097754 + 0.8508304098923269*$(x[1]) + 0.544182872577263*$(x[2]) + -0.2945971031198069*$(x[3]))/sqrt(2)))/2 + -0.8918376604733185*(0.6507439169566887 + 0.9833163152014732*$(x[1]) + -0.36891084266703134*$(x[2]) + -0.5752518327899794*$(x[3]))*(1 + erf((0.6507439169566887 + 0.9833163152014732*$(x[1]) + -0.36891084266703134*$(x[2]) + -0.5752518327899794*$(x[3]))/sqrt(2)))/2 + -0.8083246377220759*(-0.1240101222843979 + -0.304115262711357*$(x[1]) + -0.19393989104697695*$(x[2]) + -0.50508829642924*$(x[3]))*(1 + erf((-0.1240101222843979 + -0.304115262711357*$(x[1]) + -0.19393989104697695*$(x[2]) + -0.50508829642924*$(x[3]))/sqrt(2)))/2)*(1 + erf((0.47930083070023954 + -0.2811056889538768*(-0.30287094621097754 + 0.8508304098923269*$(x[1]) + 0.544182872577263*$(x[2]) + -0.2945971031198069*$(x[3]))*(1 + erf((-0.30287094621097754 + 0.8508304098923269*$(x[1]) + 0.544182872577263*$(x[2]) + -0.2945971031198069*$(x[3]))/sqrt(2)))/2 + -0.8918376604733185*(0.6507439169566887 + 0.9833163152014732*$(x[1]) + -0.36891084266703134*$(x[2]) + -0.5752518327899794*$(x[3]))*(1 + erf((0.6507439169566887 + 0.9833163152014732*$(x[1]) + -0.36891084266703134*$(x[2]) + -0.5752518327899794*$(x[3]))/sqrt(2)))/2 + -0.8083246377220759*(-0.1240101222843979 + -0.304115262711357*$(x[1]) + -0.19393989104697695*$(x[2]) + -0.50508829642924*$(x[3]))*(1 + erf((-0.1240101222843979 + -0.304115262711357*$(x[1]) + -0.19393989104697695*$(x[2]) + -0.50508829642924*$(x[3]))/sqrt(2)))/2)/sqrt(2)))/2 + (-0.740550943643274 + -0.007926005977438866*(-0.30287094621097754 + 0.8508304098923269*$(x[1]) + 0.544182872577263*$(x[2]) + -0.2945971031198069*$(x[3]))*(1 + erf((-0.30287094621097754 + 0.8508304098923269*$(x[1]) + 0.544182872577263*$(x[2]) + -0.2945971031198069*$(x[3]))/sqrt(2)))/2 + 0.328386397895855*(0.6507439169566887 + 0.9833163152014732*$(x[1]) + -0.36891084266703134*$(x[2]) + -0.5752518327899794*$(x[3]))*(1 + erf((0.6507439169566887 + 0.9833163152014732*$(x[1]) + -0.36891084266703134*$(x[2]) + -0.5752518327899794*$(x[3]))/sqrt(2)))/2 + 0.6266373732109729*(-0.1240101222843979 + -0.304115262711357*$(x[1]) + -0.19393989104697695*$(x[2]) + -0.50508829642924*$(x[3]))*(1 + erf((-0.1240101222843979 + -0.304115262711357*$(x[1]) + -0.19393989104697695*$(x[2]) + -0.50508829642924*$(x[3]))/sqrt(2)))/2)*(1 + erf((-0.740550943643274 + -0.007926005977438866*(-0.30287094621097754 + 0.8508304098923269*$(x[1]) + 0.544182872577263*$(x[2]) + -0.2945971031198069*$(x[3]))*(1 + erf((-0.30287094621097754 + 0.8508304098923269*$(x[1]) + 0.544182872577263*$(x[2]) + -0.2945971031198069*$(x[3]))/sqrt(2)))/2 + 0.328386397895855*(0.6507439169566887 + 0.9833163152014732*$(x[1]) + -0.36891084266703134*$(x[2]) + -0.5752518327899794*$(x[3]))*(1 + erf((0.6507439169566887 + 0.9833163152014732*$(x[1]) + -0.36891084266703134*$(x[2]) + -0.5752518327899794*$(x[3]))/sqrt(2)))/2 + 0.6266373732109729*(-0.1240101222843979 + -0.304115262711357*$(x[1]) + -0.19393989104697695*$(x[2]) + -0.50508829642924*$(x[3]))*(1 + erf((-0.1240101222843979 + -0.304115262711357*$(x[1]) + -0.19393989104697695*$(x[2]) + -0.50508829642924*$(x[3]))/sqrt(2)))/2)/sqrt(2)))/2 + (-0.3322291763851517 + -0.3879865017411772*(-0.30287094621097754 + 0.8508304098923269*$(x[1]) + 0.544182872577263*$(x[2]) + -0.2945971031198069*$(x[3]))*(1 + erf((-0.30287094621097754 + 0.8508304098923269*$(x[1]) + 0.544182872577263*$(x[2]) + -0.2945971031198069*$(x[3]))/sqrt(2)))/2 + 0.13951112377587105*(0.6507439169566887 + 0.9833163152014732*$(x[1]) + -0.36891084266703134*$(x[2]) + -0.5752518327899794*$(x[3]))*(1 + erf((0.6507439169566887 + 0.9833163152014732*$(x[1]) + -0.36891084266703134*$(x[2]) + -0.5752518327899794*$(x[3]))/sqrt(2)))/2 + -0.16669374723112762*(-0.1240101222843979 + -0.304115262711357*$(x[1]) + -0.19393989104697695*$(x[2]) + -0.50508829642924*$(x[3]))*(1 + erf((-0.1240101222843979 + -0.304115262711357*$(x[1]) + -0.19393989104697695*$(x[2]) + -0.50508829642924*$(x[3]))/sqrt(2)))/2)*(1 + erf((-0.3322291763851517 + -0.3879865017411772*(-0.30287094621097754 + 0.8508304098923269*$(x[1]) + 0.544182872577263*$(x[2]) + -0.2945971031198069*$(x[3]))*(1 + erf((-0.30287094621097754 + 0.8508304098923269*$(x[1]) + 0.544182872577263*$(x[2]) + -0.2945971031198069*$(x[3]))/sqrt(2)))/2 + 0.13951112377587105*(0.6507439169566887 + 0.9833163152014732*$(x[1]) + -0.36891084266703134*$(x[2]) + -0.5752518327899794*$(x[3]))*(1 + erf((0.6507439169566887 + 0.9833163152014732*$(x[1]) + -0.36891084266703134*$(x[2]) + -0.5752518327899794*$(x[3]))/sqrt(2)))/2 + -0.16669374723112762*(-0.1240101222843979 + -0.304115262711357*$(x[1]) + -0.19393989104697695*$(x[2]) + -0.50508829642924*$(x[3]))*(1 + erf((-0.1240101222843979 + -0.304115262711357*$(x[1]) + -0.19393989104697695*$(x[2]) + -0.50508829642924*$(x[3]))/sqrt(2)))/2)/sqrt(2)))/2 - $q <= 0.0))

                     @objective(m, Min, q)

                     return m

max_stars_repo_name: PSORLab/RSActivationFunctions

max_issues_repo_name: PSORLab/RSActivationFunctions

max_forks_repo_name: PSORLab/RSActivationFunctions


text:
using FiniteVolumeMethod
using Test

@testset "FiniteVolumeMethod.jl" begin

    init_steps = 3
    ref_steps = 2

    for dim=1:3
        uniform = TreeGrid(zeros(dim))
        vec = 1 : dim
        vec2 = dim : -1 : 1

        # A uniform grid
        for step=1:init_steps
            FiniteVolumeMethod.refine!(uniform, [uniform.tree], recurse=true)
        end

        # Locally refined grid
        refined = TreeGrid(zeros(dim))
        child = refined.tree
        for step=1:init_steps
            FiniteVolumeMethod.refine!(refined, [child])
            child = child.children[1+mod(step-1, 1<<dim)]
        end

        grids = [uniform, refined]
        expected_orders = [2, 1]

        div_lInf_constant = zeros(ref_steps)
        div_lInf_linear = zeros(ref_steps)
        div_lInf_nonlinear = zeros(ref_steps)
        grad_lInf_constant = zeros(ref_steps)
        grad_lInf_linear = zeros(ref_steps)
        grad_lInf_nonlinear = zeros(ref_steps)
        for (i, grid) ∈ enumerate(grids)
            for step=1:ref_steps
                FiniteVolumeMethod.refine!(grid, [grid.tree], recurse=true)

                # Gradient:
                p = CellVar(x -> pi, grid)
                g = gradient(p, grid)
                grad_lInf_constant[step] = maximum(abs.(g.data))

                p = CellVar(x -> pi + sum(vec .* x), grid)
                gExact = FaceVar(x -> vec, grid, filter = active)
                g = gradient(p, grid)
                grad_lInf_linear[step] = maximum(abs.(g.data - gExact.data))

                p = CellVar(x -> pi + sum(vec .* cos.(x)), grid)
                gExact = FaceVar(x -> -vec .* sin.(x), grid, filter = active)
                g = gradient(p, grid)
                grad_lInf_nonlinear[step] = maximum(abs.(g.data - gExact.data))

                # divergence:
                u = FaceVar(x -> vec, grid)
                d = divergence(u, grid)
                div_lInf_constant[step] = maximum(abs.(d.data))

                u = FaceVar(x -> vec .* x + vec2, grid)
                dExact = CellVar(x -> sum(vec), grid, filter = active)
                d = divergence(u, grid)
                div_lInf_linear[step] = maximum(abs.(d.data - dExact.data))

                u = FaceVar(x -> sum(x.^2) .+ vec .* x + cos.(x), grid)
                dExact = CellVar(x -> 2sum(x) + sum(vec) - sum(sin.(x)), grid, filter = active)
                d = divergence(u, grid)
                div_lInf_nonlinear[step] = maximum(abs.(d.data - dExact.data))
            end

            @test all(grad_lInf_constant .== 0.0)
            @test all(div_lInf_constant .== 0.0)

            # Approximate derivative of a linear function should be exact
            @test all(grad_lInf_linear .== 0.0)
            @test all(div_lInf_linear .== 0.0)

            # Approximate derivative of a nonlinear function should be second order accurate (uniform)
            # and first order accurate at refinement interfaces
            grad_order = log2.(grad_lInf_nonlinear[1:ref_steps-1] ./ grad_lInf_nonlinear[2:ref_steps])
            div_order = log2.(div_lInf_nonlinear[1:ref_steps-1] ./ div_lInf_nonlinear[2:ref_steps])
            @test all(grad_order .>= expected_orders[i] - 1E-1)
            @test all(div_order .>= expected_orders[i] - 1E-1)
        end
    end
end


text:
function hypot(x::DoubleFloat{T}, y::DoubleFloat{T}) where {T<:IEEEFloat}
    ax = abs(x)
    ay = abs(y)
    ay, ax = minmax(ax, ay)

    r = ay
    if !iszero(ax)
        r /= ax
    end
    r *= r

    rr = ax * sqrt(one(DoubleFloat{T}) + r)

    # from Base
    # use type of rr to make sure that return type
    #    is the same for all branches
    if isnan(ay)
        isinf(ax) && return oftype(rr, Inf)
        isinf(ay) && return oftype(rr, Inf)
        return oftype(rr, r)
    end
    return rr
 end

"""
    normalize(x,y)

    x_normalized^2 + y_normalized^2 == one(promote_type(typeof(x),typeof(y)))
"""
function normalize(x::DoubleFloat{T}, y::DoubleFloat{T}) where {T<:IEEEFloat}
    hypotenuse = hypot(x, y)
    if !iszero(hypotenuse) && isfinite(hypotenuse)
        xnorm = x / hypotenuse
        ynorm = y / hypotenuse
    else
        # try for xnorm, ynorm such that xnorm + ynorm == 1
        xnorm = x/(x+y)
        ynorm = y/(x+y)
    end
    return xnorm, ynorm
end


text:
@testset "COMPAS Dataset" begin
	X, y = @load_compas
	fpath = joinpath(Fairness.DATA_DIR, "compas-scores-two-years.csv")
	@test isfile(fpath)
	@test scitype(y) == AbstractArray{Multiclass{2}, 1}
	cols = ["sex", "age", "age_cat", "race", "c_charge_degree", "priors_count", "days_b_screening_arrest", "decile_score"]
	@test string.(names(X)) == cols
end

@testset "Adult Dataset" begin
	X, y = @load_adult
	fpath = joinpath(Fairness.DATA_DIR, "adult.data")
	@test isfile(fpath)
	@test scitype(y) == AbstractArray{Multiclass{2}, 1}
	cols = ["age", "workclass", "fnlwgt", "education",
		"education_num", "marital_status", "occupation",
		"relationship", "race", "sex", "capital_gain",
		"capital_loss", "hours_per_week", "native_country",
	]
	@test string.(names(X)) == cols
end

@testset "German Dataset" begin
	X, y = @load_german
	fpath = joinpath(Fairness.DATA_DIR, "german.data")
	@test isfile(fpath)
	@test scitype(y) == AbstractArray{Multiclass{2},1}
	cols = [
		"status", "duration", "credit_history",
		"purpose", "credit_amount", "savings", "employment",
		"installment_rate", "gender_status",
		"other_debtors", "residence_since", "property", "age",
		"installment_plans", "housing", "existing_credits",
		"skill_level", "people_liable", "telephone",
		"foreign_worker"
	]
	@test string.(names(X)) == cols
end

@testset "Bank Marketing Dataset" begin
	X, y = @load_bank_marketing
	fpath = joinpath(Fairness.DATA_DIR, "bank-marketing.csv")
	@test isfile(fpath)
	@test scitype(y) == AbstractArray{Multiclass{2},1};
end

@testset "Communities Crime Dataset" begin
	X, y = @load_communities_crime
	fpath = joinpath(Fairness.DATA_DIR, "communities_crime.data")
	@test isfile(fpath)
	@test scitype(y) == AbstractArray{Multiclass{2},1};
end

@testset "Student Performance" begin
	X, y = @load_student_performance
	fpath = joinpath(Fairness.DATA_DIR, "student-performance.csv")
	@test isfile(fpath)
	@test scitype(y) == AbstractArray{Multiclass{2},1};
end


text:
module TypeTransform

export @transform, transform

include("fexpr.jl")
include("subtypes.jl")
include("inverse_hasmethod.jl")

"""
    @transform

Transform the given type to another type during defining a method.

Use `@transform` and the function that transforms the type to another type. The function should return an `Array` of types that you want the method to be defined for.

For example, we use `allsubtypes()` type transform function to define specific methods for all of subtypes of a given type (fix ambiguity error!).
```julia
using TypeTransform
abstract type A end
abstract type B <:A end
abstract type C <:B end

@transform function foo(a, b::allsubtypes(A))
    println("a new method")
end
```
Since `allsubtypes(A)` returns the array of types `[A, B, C]`, three methods are defined
```julia
julia> methods(foo)
# 3 methods for generic function "foo":
[1] foo(a, b::C) in Main at none:2
[2] foo(a, b::B) in Main at none:2
[3] foo(a, b::A) in Main at none:2
```
Note that you could use `subtypes()` instead of `allsubtypes()`, which defines methods only for the direct subtypes (`[B]` in this case).

If you want that only specific functions to be considered in transformation by `@transform`, give an `Array` of `Symbol`s that contains the function names you want to be transformed.

```julia
@transform [:subtypes, :allsubtypes], function foo_array(a, b::allsubtypes(A))
    println("a new method")
end
```

It is possible to use the function names inside curly expressions like `Union{A, subtypes{B}}` or `Type{allsubtypes{A}}` or use arguments without a name:
```julia
@transform function foo_curly(a, ::Union{T,allsubtypes(A)}, c::T) where {T<:Int64}
    println("a new method")
end
```

# Motivation
The first motivation for this package was to fix ambiguity error by defining specific methods.

If you run the following program
```julia
abstract type A end
abstract type B <:A end

# my general vector method
foo(a::Vector, b::Type{<:A}) = print("vector method")

# my special B mwthod
foo(a, b::Type{B}) = print("B method")
```
`foo([1,2], B)` will give an ambiguity error, while if you use `allsubtypes`, you can fix the issue.

```julia
# my general vector method
@transform foo(a::Vector, b::allsubtypes(A)) = print("vector method")
```
"""
macro transform(expr::Expr)
    #TODO: support for multiple function transforms in a @transform
    modul = __module__
    macroexpand(modul, expr)
    out = transform(modul, expr)

    return out
end

function transform(modul::Module, expr)
    if expr.head == :block
        expr = expr.args[2]
    end

    if expr.head == :tuple # some functions are specified
        funclist = eval(expr.args[1])
        isfunclist = true

        fexpr = expr.args[2]
    else
        fexpr = expr
        isfunclist = false
    end

    f, args, wherestack, body = unwrap_fun(fexpr, true, true)

    fmethods = Expr[]
    for (iArg, arg) in enumerate(args)
        if arg isa Expr && arg.head == :(::) &&
            arg.args[end] isa Expr

            if arg.args[end].head == :call

                # skip this function if it is not in funclist
                if isfunclist && !(arg.args[end].args[1] in funclist)
                    continue
                end

                funcname = arg.args[end].args[1]
                intype = arg.args[end].args[2]

                isCurly =false

            elseif arg.args[end].head == :curly

                # string match is faster
                strarg = string(arg.args[end])

                # match any function name
                m = match(r"([a-zA-Z\_][a-zA-Z0-9\_]*)\((.*)\)", strarg)
                if m === nothing
                    continue
                end

                funcname = Meta.parse(m.captures[1])

                 # skip this function if it is not in funclist
                if isfunclist && !(funcname in funclist)
                    continue
                end

                intype = Meta.parse(m.captures[2])

                isCurly =true

            else
                continue
            end

            outtypes =Core.eval(modul,
                quote
                    $funcname($intype)
                end)

            outtypes_len = length(outtypes)
            fmethod = Vector{Expr}(undef, outtypes_len)
            for (iouttype, outtypeI) in enumerate(outtypes)
                # replacing with actual trasformed type
                if !isCurly
                    args[iArg].args[end] = outtypeI
                else
                    args[iArg].args[end] = Meta.parse(replace(strarg, m.match=>string(outtypeI)))
                end
                fmethod[iouttype] = copy(wrap_fun(f, args, wherestack, body))
            end
            append!(fmethods, fmethod)
        end
    end
    if isempty(fmethods)
        error("No method defined")
    end
    # print(fmethods)

    out = quote
        Base.@__doc__(function $f end) # supports docuementation
        $(esc.(fmethods)...)
    end
    return out
end

end