text:
@testset "gd/sgd" begin

    # setup
    Random.seed!(123)
    kernel = "../src/pca/pca.jl"
    nworkers = 2
    niterations = 200
    inputdataset = "X"
    outputdataset = "V"
    nsamples, dimension = 20, 10
    ncomponents = div(dimension, 2)

    # generate input dataset
    X = randn(dimension, nsamples)
    inputfile = tempname()
    h5open(inputfile, "w") do file
        file[inputdataset] = X
        flush(file)
    end

    # correct solution (computed via LinearAlgebra.svd)
    V = pca(X, ncomponents)
    ev = explained_variance(X, V)

    ### exact
    outputfile = tempname()
    mpiexec(cmd -> run(```
        $cmd -n $(nworkers+1) julia --project $kernel $inputfile $outputfile
        --niterations $niterations
        --ncomponents $ncomponents
        --saveiterates
        --loadbalance
        ```))
    Vs, _ = test_pca_iterates(;X, niterations, ncomponents, ev, outputfile, outputdataset)

    ### providing an initial iterate
    iteratedataset = "V"
    niterations = 1
    h5open(inputfile, "r+") do fid
        fid[iteratedataset] = Vs[end]
    end
    outputfile = tempname()
    mpiexec(cmd -> run(```
        $cmd -n $(nworkers+1) julia --project $kernel $inputfile $outputfile
        --niterations $niterations
        --ncomponents $ncomponents
        --saveiterates
        --iteratedataset $iteratedataset
        --loadbalance
        ```))
    test_pca_iterates(;X, niterations, ncomponents, ev, outputfile, outputdataset)

    ### exact with replication
    nreplicas = 2
    niterations = 200
    outputfile = tempname()
    mpiexec(cmd -> run(```
        $cmd -n $(nworkers+1) julia --project $kernel $inputfile $outputfile
        --niterations $niterations
        --ncomponents $ncomponents
        --nreplicas $nreplicas
        --saveiterates
        --loadbalance
        ```))
    test_pca_iterates(;X, niterations, ncomponents, ev, outputfile, outputdataset)

    ### replication + ignoring the slowest worker
    nworkers = 2
    outputfile = tempname()
    mpiexec(cmd -> run(```
        $cmd -n $(nworkers+1) julia --project $kernel $inputfile $outputfile
        --ncomponents $ncomponents
        --niterations $niterations
        --nwait $(nworkers-1)
        --nreplicas $nreplicas
        --saveiterates
        --loadbalance
        ```))
    test_pca_iterates(;X, niterations, ncomponents, ev, outputfile, outputdataset)

    ### 12 workers, a factor 3 replication
    nworkers = 12
    nreplicas = 3
    nwait = div(nworkers, nreplicas)
    outputfile = tempname()
    mpiexec(cmd -> run(```
        $cmd -n $(nworkers+1) julia --project $kernel $inputfile $outputfile
        --niterations $niterations
        --ncomponents $ncomponents
        --nreplicas $nreplicas
        --nwait $nwait
        --saveiterates
        --loadbalance
        ```))
    test_pca_iterates(;X, niterations, ncomponents, ev, outputfile, outputdataset)

    ## using mini-batches
    nworkers = 2
    outputfile = tempname()
    stepsize = 1.0
    mpiexec(cmd -> run(```
        $cmd -n $(nworkers+1) julia --project $kernel $inputfile $outputfile
        --ncomponents $ncomponents
        --niterations $niterations
        --nwait $(nworkers-1)
        --stepsize $stepsize
        --saveiterates
        --loadbalance
        ```))
    test_pca_iterates(;X, niterations, ncomponents, ev, outputfile, outputdataset, atol=Inf)

    ### sub-partitioning and stochastic sub-gradients
    nsubpartitions = 2
    stepsize = 1.0 / nsubpartitions
    outputfile = tempname()
    mpiexec(cmd -> run(```
        $cmd -n $(nworkers+1) julia --project $kernel $inputfile $outputfile
        --ncomponents $ncomponents
        --niterations $niterations
        --nwait $(nworkers-1)
        --nsubpartitions $nsubpartitions
        --stepsize $stepsize
        --saveiterates
        --loadbalance
        ```))
    test_pca_iterates(;X, niterations, ncomponents, ev, outputfile, outputdataset, atol=Inf)
end

@time @testset "sparse matrices" begin

    # setup
    Random.seed!(123)
    kernel = "../src/pca/pca.jl"
    nworkers = 2
    niterations = 200
    inputdataset = "X"
    outputdataset = "V"
    nsamples, dimension = 20, 10
    ncomponents = div(dimension, 2)
    p = 0.9 # matrix density

    # generate input dataset
    X = sprand(dimension, nsamples, p)
    inputfile = tempname()
    !isfile(inputfile) || rm(inputfile)
    # fid = h5open(inputfile, "cw")

    # # h5writecsc(inputfile, inputdataset, X)
    A = H5SparseMatrixCSC(inputfile, inputdataset, X)
    flush(A.fid)

    # exact solution (computed via LinearAlgebra.svd)
    V = pca(Matrix(X), ncomponents)
    ev = explained_variance(X, V)

    # exact solution
    outputfile = tempname()
    mpiexec(cmd -> run(```
        $cmd -n $(nworkers+1) julia --project $kernel $inputfile $outputfile
        --niterations $niterations
        --ncomponents $ncomponents
        --saveiterates
        --loadbalance
        ```))
    test_pca_iterates(;X, niterations, ncomponents, ev, outputfile, outputdataset)
end

@time @testset "variance-reduced" begin
    # with variance reduction, the algorithm will always converge to the optimum

    # setup
    Random.seed!(123)
    kernel = "../src/pca/pca.jl"
    nworkers = 2
    niterations = 100
    inputdataset = "X"
    inputdataset_sparse = "S"
    outputdataset = "V"
    nsamples, dimension = 20, 10
    ncomponents = div(dimension, 2)

    # generate input dataset
    X = randn(dimension, nsamples)
    inputfile = tempname()
    h5open(inputfile, "w") do file
        file[inputdataset] = X
        H5SparseMatrixCSC(file, inputdataset_sparse, sparse(X))
    end
    V = pca(X, ncomponents)
    ev = explained_variance(X, V)

    ### partitioning the dataset over the workers
    stepsize = 1
    outputfile = tempname()
    nsubpartitions = 1
    mpiexec(cmd -> run(```$cmd -n $(nworkers+1) julia --project $kernel $inputfile $outputfile
        --ncomponents $ncomponents
        --niterations $niterations
        --stepsize $stepsize
        --nsubpartitions $nsubpartitions
        --nwait $(nworkers-1)
        --variancereduced
        --saveiterates
        ```))
    test_pca_iterates(;X, niterations, ncomponents, ev, outputfile, outputdataset)

    ### with a factor 2 replication
    outputfile = tempname()
    nsubpartitions = 1
    nworkers = 4
    nreplicas = 2
    npartitions = div(nworkers, nreplicas)
    mpiexec(cmd -> run(```$cmd -n $(nworkers+1) julia --project $kernel $inputfile $outputfile
        --ncomponents $ncomponents
        --niterations $niterations
        --stepsize $stepsize
        --nsubpartitions $nsubpartitions
        --nwait $(npartitions-1)
        --variancereduced
        --saveiterates
        ```))
    test_pca_iterates(;X, niterations, ncomponents, ev, outputfile, outputdataset)

    ### same as the previous, but with kickstart enabled
    outputfile = tempname()
    mpiexec(cmd -> run(```$cmd -n $(nworkers+1) julia --project $kernel $inputfile $outputfile
        --ncomponents $ncomponents
        --niterations $niterations
        --stepsize $stepsize
        --nsubpartitions $nsubpartitions
        --nwait $(npartitions-1)
        --variancereduced
        --kickstart
        --saveiterates
        ```))
    test_pca_iterates(;X, niterations, ncomponents, ev, outputfile, outputdataset)

    ### same as the previous, but with nwaitschedule < 1
    nwaitschedule = 0.9
    outputfile = tempname()
    mpiexec(cmd -> run(```$cmd -n $(nworkers+1) julia --project $kernel $inputfile $outputfile
        --ncomponents $ncomponents
        --niterations $niterations
        --stepsize $stepsize
        --nsubpartitions $nsubpartitions
        --nwait $(npartitions-1)
        --variancereduced
        --kickstart
        --saveiterates
        --nwaitschedule $nwaitschedule
        ```))
    test_pca_iterates(;X, niterations, ncomponents, ev, outputfile, outputdataset)

    ### with a factor 3 replication
    outputfile = tempname()
    nsubpartitions = 1
    nworkers = 3
    nreplicas = 3
    npartitions = div(nworkers, nreplicas)
    mpiexec(cmd -> run(```$cmd -n $(nworkers+1) julia --project $kernel $inputfile $outputfile
        --ncomponents $ncomponents
        --niterations $niterations
        --stepsize $stepsize
        --nsubpartitions $nsubpartitions
        --nwait $npartitions
        --variancereduced
        --saveiterates
        ```))
    test_pca_iterates(;X, niterations, ncomponents, ev, outputfile, outputdataset)

    ### sub-partitioning the data stored at each worker
    nworkers = 2
    niterations = 100
    nsubpartitions = 2
    stepsize = 1/2
    outputfile = tempname()
    mpiexec(cmd -> run(```
        $cmd -n $(nworkers+1) julia --project $kernel $inputfile $outputfile
        --ncomponents $ncomponents
        --niterations $niterations
        --stepsize $stepsize
        --nwait $(nworkers-1)
        --nsubpartitions $nsubpartitions
        --variancereduced
        --saveiterates
        ```))
    test_pca_iterates(;X, niterations, ncomponents, ev, outputfile, outputdataset)

    ### same as previous, but for a sparse input matrix
    outputfile = tempname()
    mpiexec(cmd -> run(```
        $cmd -n $(nworkers+1) julia --project $kernel $inputfile $outputfile
        --inputdataset $inputdataset_sparse
        --ncomponents $ncomponents
        --niterations $niterations
        --stepsize $stepsize
        --nwait $(nworkers-1)
        --nsubpartitions $nsubpartitions
        --variancereduced
        --saveiterates
        ```))
    test_pca_iterates(;X, niterations, ncomponents, ev, outputfile, outputdataset)

    # with replication
    nworkers = 4
    nreplicas = 2
    npartitions = div(nworkers, nreplicas)
    niterations = 100
    stepsize = 1/2
    nsubpartitions = 2
    outputfile = tempname()
    mpiexec(cmd -> run(```
        $cmd -n $(nworkers+1) julia --project $kernel $inputfile $outputfile
        --ncomponents $ncomponents
        --niterations $niterations
        --stepsize $stepsize
        --nwait $(npartitions-1)
        --nsubpartitions $nsubpartitions
        --nreplicas $nreplicas
        --variancereduced
        --saveiterates
        ```))
    test_pca_iterates(;X, niterations, ncomponents, ev, outputfile, outputdataset)
end


text:
# Acknowledgement: this adaption settings is mimicing Stan's 3-phase adaptation.

# Ref： https://github.com/stan-dev/math/blob/develop/stan/math/prim/mat/fun/welford_var_estimator.hpp
mutable struct VarEstimator{T<:Real}
    n::Int
    μ::Vector{T}
    M::Vector{T}
end

function reset!(ve::VarEstimator{T}) where T
    ve.n = 0
    ve.μ .= zero(T)
    ve.M .= zero(T)
    return ve
end

function add_sample!(ve::VarEstimator, s::AbstractVector)
    ve.n += 1
    δ = s .- ve.μ
    ve.μ .+= δ ./ ve.n
    ve.M .+= δ .* (s .- ve.μ)
    return ve
end

# https://github.com/stan-dev/stan/blob/develop/src/stan/mcmc/var_adaptation.hpp
function get_var(ve::VarEstimator)
    @assert ve.n >= 2
    return (ve.n / ((ve.n + 5) * (ve.n - 1))) .* ve.M .+ 1e-3 * (5.0 / (ve.n + 5))
end

mutable struct WarmUpManager
    adapt_n::Int
    params::Dict{Symbol, Any}
    ve::VarEstimator
end

getindex(wum::WarmUpManager, param) = wum.params[param]

setindex!(wum::WarmUpManager, value, param) = wum.params[param] = value

function init_warm_up_params(vi::VarInfo, spl::Sampler{<:Hamiltonian})
    D = length(vi[spl])
    ve = VarEstimator{Float64}(0, zeros(D), zeros(D))

    wum = WarmUpManager(1, Dict(), ve)

    # Pre-cond
    wum[:stds] = ones(D)

    # Dual averaging
    wum[:ϵ] = [] # why we are using a vector for ϵ
    restart_da(wum)
    wum[:n_warmup] = spl.alg.n_adapt
    wum[:δ] = spl.alg.delta

    # Initialize by Stan if Stan is installed
    is_init_by_stan = false
    @require Stan="682df890-35be-576f-97d0-3d8c8b33a550" begin
        is_init_by_stan = true

        # Stan.Adapt
        adapt_conf = spl.info[:adapt_conf]

        # Hyper parameters for dual averaging
        wum[:γ] = adapt_conf.gamma
        wum[:t_0] = adapt_conf.t0
        wum[:κ] = adapt_conf.kappa

        # Three phases settings
        wum[:init_buffer] = adapt_conf.init_buffer
        wum[:term_buffer] = adapt_conf.term_buffer
        wum[:window_size] = adapt_conf.window
    end

    # If wum is not initialised by Stan (when Stan is not avaible),
    # initialise wum by common default values.
    if ~is_init_by_stan
        wum[:γ] = 0.05
        wum[:t_0] = 10.0
        wum[:κ] = 0.75
        wum[:init_buffer] = 75
        wum[:term_buffer] = 50
        wum[:window_size] = 25
    end
    wum[:next_window] = wum[:init_buffer] + wum[:window_size] - 1

    @debug wum.params

    spl.info[:wum] = wum
end

function restart_da(wum::WarmUpManager)
    wum[:m] = 0
    wum[:x_bar] = 0.0
    wum[:H_bar] = 0.0
end

# See NUTS paper sec 3.2.1
function update_da_μ(wum::WarmUpManager, ϵ::Float64)
    wum[:μ] = log(10 * ϵ)
end

# Ref: https://github.com/stan-dev/stan/blob/develop/src/stan/mcmc/stepsize_adaptation.hpp
function adapt_step_size!(wum::WarmUpManager, stats::Real)

    @debug "adapting step size ϵ..."
    @debug "current α = $(stats)"
    wum[:m] = wum[:m] + 1
    m = wum[:m]

    # Clip average MH acceptance probability.
    stats = stats > 1 ? 1 : stats

    γ = wum[:γ]; t_0 = wum[:t_0]; κ = wum[:κ]; δ = wum[:δ]
    μ = wum[:μ]; x_bar = wum[:x_bar]; H_bar = wum[:H_bar]

    η_H = 1.0 / (m + t_0)
    H_bar = (1.0 - η_H) * H_bar + η_H * (δ - stats)

    x = μ - H_bar * sqrt(m) / γ            # x ≡ logϵ
    η_x = m^(-κ)
    x_bar = (1.0 - η_x) * x_bar + η_x * x

    ϵ = exp(x)
    @debug "new ϵ = $(ϵ), old ϵ = $(wum[:ϵ][end])"

    if isnan(ϵ) || isinf(ϵ) || ϵ <= 1e-3
        @warn "Incorrect ϵ = $ϵ; ϵ_previous = $(wum[:ϵ][end]) is used instead."
    else
        push!(wum[:ϵ], ϵ)
        wum[:x_bar], wum[:H_bar] = x_bar, H_bar
    end

    if m == wum[:n_warmup]
        @debug " Adapted ϵ = $ϵ, $m HMC iterations is used for adaption."
    end

end

# Ref: https://github.com/stan-dev/stan/blob/develop/src/stan/mcmc/windowed_adaptation.hpp
function in_adaptation(wum::WarmUpManager)
    return (wum.adapt_n >= wum[:init_buffer]) &&
        (wum.adapt_n < wum[:n_warmup] - wum[:term_buffer]) &&
        (wum.adapt_n != wum[:n_warmup])
end

function is_window_end(wum::WarmUpManager)
    return (wum.adapt_n == wum[:next_window]) && (wum.adapt_n != wum[:n_warmup])
end

function compute_next_window(wum::WarmUpManager)

    if ~(wum[:next_window] == wum[:n_warmup] - wum[:term_buffer] - 1)

        wum[:window_size] *= 2
        wum[:next_window] = wum.adapt_n + wum[:window_size]

        if ~(wum[:next_window] == wum[:n_warmup] - wum[:term_buffer] - 1)
            next_window_boundary = wum[:next_window] + 2 * wum[:window_size]

            if (next_window_boundary >= wum[:n_warmup] - wum[:term_buffer])
                wum[:next_window] = wum[:n_warmup] - wum[:term_buffer] - 1
            end
        end
    end
end

# Ref: https://github.com/stan-dev/stan/blob/develop/src/stan/mcmc/hmc/nuts/adapt_diag_e_nuts.hpp
function update_pre_cond!(wum::WarmUpManager, θ_new::AbstractVector{<:Real})

  if in_adaptation(wum)
    add_sample!(wum.ve, θ_new)
  end
    if is_window_end(wum)
        compute_next_window(wum)
        var = get_var(wum.ve)
        wum[:stds] = sqrt.(var)
        reset!(wum.ve)
        return true
    end
    return false
end

function adapt!(wum::WarmUpManager, stats::Real, θ_new; adapt_ϵ=false, adapt_M=false)

    if wum.adapt_n < wum[:n_warmup]

        if adapt_ϵ
            adapt_step_size!(wum, stats)
            if is_window_end(wum)
                ϵ = exp(wum[:x_bar])
                push!(wum[:ϵ], ϵ)
                update_da_μ(wum, ϵ)
                restart_da(wum)
            end
        end

        if adapt_M
            update_pre_cond!(wum, θ_new)  # window is updated implicitly.
        else   # update window explicitly.
            is_window_end(wum) && compute_next_window(wum)
        end

        wum.adapt_n += 1

    elseif wum.adapt_n == wum[:n_warmup]

        if adapt_ϵ
            ϵ = exp(wum[:x_bar])
            push!(wum[:ϵ], ϵ)
        end
    end
end

ext: jl


text:
@testset "DisplayArray" begin

    array = DS.DisplayArray(10, 40, background='.')

    @test repr(array) == "DisplayArray(" *
        "size=(10, 40), " *
        "background char=Char(46)" *
    ")"

    DS.render(array)
    @test T.read_stream(T.out_stream)  ==
        "\e[s........................................\e[u\e[1B" *
        "\e[s........................................\e[u\e[1B" *
        "\e[s........................................\e[u\e[1B" *
        "\e[s........................................\e[u\e[1B" *
        "\e[s........................................\e[u\e[1B" *
        "\e[s........................................\e[u\e[1B" *
        "\e[s........................................\e[u\e[1B" *
        "\e[s........................................\e[u\e[1B" *
        "\e[s........................................\e[u\e[1B" *
        "\e[s........................................\e[u\e[1B"

    array[6, 8] = '噢'
    DS.render(array)
    @test T.read_stream(T.out_stream)  ==
        "\e[s........................................\e[u\e[1B" *
        "\e[s........................................\e[u\e[1B" *
        "\e[s........................................\e[u\e[1B" *
        "\e[s........................................\e[u\e[1B" *
        "\e[s........................................\e[u\e[1B" *
        "\e[s.......噢...............................\e[u\e[1B" *
        "\e[s........................................\e[u\e[1B" *
        "\e[s........................................\e[u\e[1B" *
        "\e[s........................................\e[u\e[1B" *
        "\e[s........................................\e[u\e[1B"

    array[4:6, 8:19] = "這是一段\n非常感性的\n字串"
    DS.render(array)
    @test T.read_stream(T.out_stream)  ==
        "\e[s........................................\e[u\e[1B" *
        "\e[s........................................\e[u\e[1B" *
        "\e[s........................................\e[u\e[1B" *
        "\e[s.......這是一段.........................\e[u\e[1B" *
        "\e[s.......非常感性的.......................\e[u\e[1B" *
        "\e[s.......字串.............................\e[u\e[1B" *
        "\e[s........................................\e[u\e[1B" *
        "\e[s........................................\e[u\e[1B" *
        "\e[s........................................\e[u\e[1B" *
        "\e[s........................................\e[u\e[1B"

    # pos
    DS.render(array, pos=(5, 6))
    @test T.read_stream(T.out_stream)  ==
        "\e[5;6H" *
        "\e[s........................................\e[u\e[1B" *
        "\e[s........................................\e[u\e[1B" *
        "\e[s........................................\e[u\e[1B" *
        "\e[s.......這是一段.........................\e[u\e[1B" *
        "\e[s.......非常感性的.......................\e[u\e[1B" *
        "\e[s.......字串.............................\e[u\e[1B" *
        "\e[s........................................\e[u\e[1B" *
        "\e[s........................................\e[u\e[1B" *
        "\e[s........................................\e[u\e[1B" *
        "\e[s........................................\e[u\e[1B"

    # last index
    array[2, end-1] = '哈'
    DS.render(array)
    @test T.read_stream(T.out_stream)  ==
        "\e[s........................................\e[u\e[1B" *
        "\e[s......................................哈\e[u\e[1B" *
        "\e[s........................................\e[u\e[1B" *
        "\e[s.......這是一段.........................\e[u\e[1B" *
        "\e[s.......非常感性的.......................\e[u\e[1B" *
        "\e[s.......字串.............................\e[u\e[1B" *
        "\e[s........................................\e[u\e[1B" *
        "\e[s........................................\e[u\e[1B" *
        "\e[s........................................\e[u\e[1B" *
        "\e[s........................................\e[u\e[1B"

    array[3, (end-4):end] = "哈哈"
    DS.render(array)
    @test T.read_stream(T.out_stream)  ==
        "\e[s........................................\e[u\e[1B" *
        "\e[s......................................哈\e[u\e[1B" *
        "\e[s...................................哈哈.\e[u\e[1B" *
        "\e[s.......這是一段.........................\e[u\e[1B" *
        "\e[s.......非常感性的.......................\e[u\e[1B" *
        "\e[s.......字串.............................\e[u\e[1B" *
        "\e[s........................................\e[u\e[1B" *
        "\e[s........................................\e[u\e[1B" *
        "\e[s........................................\e[u\e[1B" *
        "\e[s........................................\e[u\e[1B"

    array[2:end-1, 3] = "直\n的\n文\n字"
    DS.render(array)
    @test T.read_stream(T.out_stream)  ==
        "\e[s........................................\e[u\e[1B" *
        "\e[s..直..................................哈\e[u\e[1B" *
        "\e[s..的...............................哈哈.\e[u\e[1B" *
        "\e[s..文...這是一段.........................\e[u\e[1B" *
        "\e[s..字...非常感性的.......................\e[u\e[1B" *
        "\e[s.......字串.............................\e[u\e[1B" *
        "\e[s........................................\e[u\e[1B" *
        "\e[s........................................\e[u\e[1B" *
        "\e[s........................................\e[u\e[1B" *
        "\e[s........................................\e[u\e[1B"

end


text:
__precompile__()
module AmplNLWriter

using MathProgBase
importall MathProgBase.SolverInterface

debug = false
setdebug(b::Bool) = global debug = b

solverdata_dir = joinpath(Pkg.dir("AmplNLWriter"), ".solverdata")

include("nl_linearity.jl")
include("nl_params.jl")
include("nl_convert.jl")

export AmplNLSolver, BonminNLSolver, CouenneNLSolver, IpoptNLSolver,
       getsolvername, getsolveresult, getsolveresultnum, getsolvemessage,
       getsolveexitcode

immutable AmplNLSolver <: AbstractMathProgSolver
    solver_command::AbstractString
    options::Vector{ASCIIString}
    filename::AbstractString
end

osl = isdir(Pkg.dir("CoinOptServices"))
ipt = isdir(Pkg.dir("Ipopt"))

if osl; import CoinOptServices; end
if ipt; import Ipopt; end

function AmplNLSolver(solver_command::AbstractString;
                      filename::AbstractString="")
  AmplNLSolver(solver_command,
               ASCIIString[],
               filename=filename)
end

function AmplNLSolver(solver_command::AbstractString,
                      options::Vector{ASCIIString}=ASCIIString[];
                      filename::AbstractString="")
    AmplNLSolver(solver_command, options, filename)
end

function BonminNLSolver(options=ASCIIString[]; filename::AbstractString="")
    osl || error("CoinOptServices not installed. Please run\n",
                 "Pkg.add(\"CoinOptServices\")")
    AmplNLSolver(CoinOptServices.bonmin, options; filename=filename)
end
function CouenneNLSolver(options=ASCIIString[]; filename::AbstractString="")
    osl || error("CoinOptServices not installed. Please run\n",
                 "Pkg.add(\"CoinOptServices\")")
    AmplNLSolver(CoinOptServices.couenne, options; filename=filename)
end
function IpoptNLSolver(options=ASCIIString[]; filename::AbstractString="")
    ipt || error("Ipopt not installed. Please run\nPkg.add(\"Ipopt\")")
    AmplNLSolver(Ipopt.amplexe, options; filename=filename)
end

getsolvername(s::AmplNLSolver) = basename(s.solver_command)

type AmplNLMathProgModel <: AbstractMathProgModel
    options::Vector{ASCIIString}

    solver_command::AbstractString

    x_l::Vector{Float64}
    x_u::Vector{Float64}
    g_l::Vector{Float64}
    g_u::Vector{Float64}

    nvar::Int
    ncon::Int

    obj
    constrs::Vector{Any}

    lin_constrs::Vector{Dict{Int, Float64}}
    lin_obj::Dict{Int, Float64}

    r_codes::Vector{Int}
    j_counts::Vector{Int}

    vartypes::Vector{Symbol}
    varlinearities_con::Vector{Symbol}
    varlinearities_obj::Vector{Symbol}
    conlinearities::Vector{Symbol}
    objlinearity::Symbol

    v_index_map::Dict{Int, Int}
    v_index_map_rev::Dict{Int, Int}
    c_index_map::Dict{Int, Int}
    c_index_map_rev::Dict{Int, Int}

    sense::Symbol

    x_0::Vector{Float64}

    file_basename::AbstractString
    probfile::AbstractString
    solfile::AbstractString

    objval::Float64
    solution::Vector{Float64}

    status::Symbol
    solve_exitcode::Int
    solve_result_num::Int
    solve_result::AbstractString
    solve_message::AbstractString
    solve_time::Float64

    d::AbstractNLPEvaluator

    function AmplNLMathProgModel(solver_command::AbstractString,
                                 options::Vector{ASCIIString},
                                 filename::AbstractString)
        new(options,
            solver_command,
            zeros(0),
            zeros(0),
            zeros(0),
            zeros(0),
            0,
            0,
            :(0),
            [],
            Dict{Int, Float64}[],
            Dict{Int, Float64}(),
            Int[],
            Int[],
            Symbol[],
            Symbol[],
            Symbol[],
            Symbol[],
            :Lin,
            Dict{Int, Int}(),
            Dict{Int, Int}(),
            Dict{Int, Int}(),
            Dict{Int, Int}(),
            :Min,
            zeros(0),
            filename,
            "",
            "",
            NaN,
            zeros(0),
            :NotSolved,
            -1,
            -1,
            "?",
            "",
            NaN)
    end
end
type AmplNLLinearQuadraticModel <: AbstractLinearQuadraticModel
    inner::AmplNLMathProgModel
end
type AmplNLNonlinearModel <: AbstractNonlinearModel
    inner::AmplNLMathProgModel
end

include("nl_write.jl")

NonlinearModel(s::AmplNLSolver) = AmplNLNonlinearModel(
    AmplNLMathProgModel(s.solver_command, s.options, s.filename)
)
LinearQuadraticModel(s::AmplNLSolver) = AmplNLLinearQuadraticModel(
    AmplNLMathProgModel(s.solver_command, s.options, s.filename)
)

function loadproblem!(outer::AmplNLNonlinearModel, nvar::Integer, ncon::Integer,
                      x_l, x_u, g_l, g_u, sense::Symbol,
                      d::AbstractNLPEvaluator)
    m = outer.inner

    m.nvar, m.ncon = nvar, ncon
    loadcommon!(m, x_l, x_u, g_l, g_u, sense)

    m.d = d
    initialize(m.d, [:ExprGraph])

    # Process constraints
    m.constrs = map(1:m.ncon) do i
        c = constr_expr(m.d, i)

        # Remove relations and bounds from constraint expressions
        if length(c.args) == 3
            if VERSION < v"0.5-"
                expected_head = :comparison
                expr_index = 1
                rel_index = 2
            else
                expected_head = :call
                expr_index = 2
                rel_index = 1
            end

            @assert c.head == expected_head
            # Single relation constraint: expr rel bound
            rel = c.args[rel_index]
            m.r_codes[i] = relation_to_nl[rel]
            if rel == [:<=, :(==)]
                m.g_u[i] = c.args[3]
            end
            if rel in [:>=, :(==)]
                m.g_l[i] = c.args[3]
            end
            c = c.args[expr_index]
        else
            # Double relation constraint: bound <= expr <= bound
            @assert c.head == :comparison
            m.r_codes[i] = relation_to_nl[:multiple]
            m.g_u[i] = c.args[5]
            m.g_l[i] = c.args[1]
            c = c.args[3]
        end

        # Convert non-linear expression to non-linear, linear and constant
        c, constant, m.conlinearities[i] = process_expression!(
            c, m.lin_constrs[i], m.varlinearities_con)

        # Update bounds on constraint
        m.g_l[i] -= constant
        m.g_u[i] -= constant

        # Update jacobian counts using the linear constraint variables
        for j in keys(m.lin_constrs[i])
            m.j_counts[j] += 1
        end
        c
    end

    # Process objective
    m.obj = obj_expr(m.d)
    if length(m.obj.args) < 2
        m.obj = nothing
    else
        # Convert non-linear expression to non-linear, linear and constant
        m.obj, constant, m.objlinearity = process_expression!(
            m.obj, m.lin_obj, m.varlinearities_obj)

        # Add constant back into non-linear expression
        if constant != 0
            m.obj = add_constant(m.obj, constant)
        end
    end
    m
end

function loadproblem!(outer::AmplNLLinearQuadraticModel, A::AbstractMatrix,
                      x_l, x_u, c, g_l, g_u, sense)
    m = outer.inner
    m.ncon, m.nvar = size(A)

    loadcommon!(m, x_l, x_u, g_l, g_u, sense)

    # Load A into the linear constraints
    @assert (m.ncon, m.nvar) == size(A)
    load_A!(m, A)
    m.constrs = zeros(m.ncon)  # Dummy constraint expression trees

    # Load c
    for (index, val) in enumerate(c)
        m.lin_obj[index] = val
    end
    m.obj = 0  # Dummy objective expression tree

    # Process variables bounds
    for j = 1:m.ncon
        lower = m.g_l[j]
        upper = m.g_u[j]
        if lower == -Inf
            if upper == Inf
                error("Neither lower nor upper bound on constraint $j")
            else
                m.r_codes[j] = 1
            end
        else
            if lower == upper
                m.r_codes[j] = 4
            elseif upper == Inf
                m.r_codes[j] = 2
            else
                m.r_codes[j] = 0
            end
        end
    end
    m
end

function load_A!(m::AmplNLMathProgModel, A::SparseMatrixCSC{Float64})
    for var = 1:A.n, k = A.colptr[var] : (A.colptr[var + 1] - 1)
        m.lin_constrs[A.rowval[k]][var] = A.nzval[k]
        m.j_counts[var] += 1
    end
end

function load_A!(m::AmplNLMathProgModel, A::Matrix{Float64})
    for con = 1:m.ncon, var = 1:m.nvar
        val = A[con, var]
        if val != 0
            m.lin_constrs[A.rowval[k]][var] = A.nzval[k]
            m.j_counts[var] += 1
        end
    end
end

function loadcommon!(m::AmplNLMathProgModel, x_l, x_u, g_l, g_u, sense)
    @assert m.nvar == length(x_l) == length(x_u)
    @assert m.ncon == length(g_l) == length(g_u)

    m.x_l, m.x_u = x_l, x_u
    m.g_l, m.g_u = g_l, g_u
    setsense!(m, sense)

    m.lin_constrs = [Dict{Int, Float64}() for _ in 1:m.ncon]
    m.j_counts = zeros(Int, m.nvar)

    m.r_codes = Array(Int, m.ncon)

    m.varlinearities_con = fill(:Lin, m.nvar)
    m.varlinearities_obj = fill(:Lin, m.nvar)
    m.conlinearities = fill(:Lin, m.ncon)
    m.objlinearity = :Lin

    m.vartypes = fill(:Cont, m.nvar)
    m.x_0 = zeros(m.nvar)
end

getvartype(m::AmplNLMathProgModel) = copy(m.vartypes)
function setvartype!(m::AmplNLMathProgModel, cat::Vector{Symbol})
    @assert all(x-> (x in [:Cont,:Bin,:Int]), cat)
    m.vartypes = copy(cat)
end

getsense(m::AmplNLMathProgModel) = m.sense
function setsense!(m::AmplNLMathProgModel, sense::Symbol)
    @assert sense == :Min || sense == :Max
    m.sense = sense
end

setwarmstart!(m::AmplNLMathProgModel, v::Vector{Float64}) = m.x_0 = v

function optimize!(m::AmplNLMathProgModel)
    m.status = :NotSolved
    m.solve_exitcode = -1
    m.solve_result_num = -1
    m.solve_result = "?"
    m.solve_message = ""

    # There is no non-linear binary type, only non-linear discrete, so make
    # sure binary vars have bounds in [0, 1]
    for i in 1:m.nvar
        if m.vartypes[i] == :Bin
            if m.x_l[i] < 0
                m.x_l[i] = 0
            end
            if m.x_u[i] > 1
                m.x_u[i] = 1
            end
        end
    end

    make_var_index!(m)
    make_con_index!(m)

    if length(m.file_basename) == 0
        # No filename specified - write to a randomly named file
        file_basepath, f_prob = mktemp(solverdata_dir)
    else
        file_basepath = joinpath(solverdata_dir, "$(m.file_basename)")
        f_prob = open(file_basepath, "w")
    end
    m.probfile = "$file_basepath.nl"
    m.solfile = "$file_basepath.sol"

    write_nl_file(f_prob, m)
    close(f_prob)

    # Rename file to have .nl extension (this is required by solvers)
    mv(file_basepath, m.probfile)

    # Run solver and save exitcode
    t = time()
    proc = spawn(pipeline(
        `$(m.solver_command) $(m.probfile) -AMPL $(m.options)`, stdout=STDOUT))
    wait(proc)
    kill(proc)
    m.solve_exitcode = proc.exitcode
    m.solve_time = time() - t

    if m.solve_exitcode == 0
        read_results(m)
    else
        m.status = :Error
        m.solution = fill(NaN,m.nvar)
        m.solve_result = "failure"
        m.solve_result_num = 999
    end

    # Clean up temp files
    if !debug
        for temp_file in [m.probfile; m.solfile]
            isfile(temp_file) && rm(temp_file)
        end
    end
end

function process_expression!(nonlin_expr::Expr, lin_expr::Dict{Int, Float64},
                             varlinearities::Vector{Symbol})
    # Get list of all variables in the expression
    extract_variables!(lin_expr, nonlin_expr)
    # Extract linear and constant terms from non-linear expression
    tree = LinearityExpr(nonlin_expr)
    tree = pull_up_constants(tree)
    _, tree, constant = prune_linear_terms!(tree, lin_expr)
    # Make sure all terms remaining in the tree are .nl-compatible
    nonlin_expr = convert_formula(tree)

    # Track which variables appear nonlinearly
    nonlin_vars = Dict{Int, Float64}()
    extract_variables!(nonlin_vars, nonlin_expr)
    for j in keys(nonlin_vars)
        varlinearities[j] = :Nonlin
    end

    # Remove variables at coeff 0 that aren't also in the nonlinear tree
    for (j, coeff) in lin_expr
        if coeff == 0 && !(j in keys(nonlin_vars))
            delete!(lin_expr, j)
        end
    end

    # Mark constraint as nonlinear if anything is left in the tree
    linearity = nonlin_expr != 0 ? :Nonlin : :Lin

    return nonlin_expr, constant, linearity
end
function process_expression!(nonlin_expr::Real, lin_expr, varlinearities)
    # Special case where body of constraint is constant
    # Return empty nonlinear and linear parts, and use the body as the constant
    0, nonlin_expr, :Lin
end

status(m::AmplNLMathProgModel) = m.status
getsolution(m::AmplNLMathProgModel) = copy(m.solution)
getobjval(m::AmplNLMathProgModel) = m.objval
numvar(m::AmplNLMathProgModel) = m.nvar
numconstr(m::AmplNLMathProgModel) = m.ncon
getsolvetime(m::AmplNLMathProgModel) = m.solve_time

# Access to AMPL solve result items
get_solve_result(m::AmplNLMathProgModel) = m.solve_result
get_solve_result_num(m::AmplNLMathProgModel) = m.solve_result_num
get_solve_message(m::AmplNLMathProgModel) = m.solve_message
get_solve_exitcode(m::AmplNLMathProgModel) = m.solve_exitcode

# We need to track linear coeffs of all variables present in the expression tree
extract_variables!(lin_constr::Dict{Int, Float64}, c) = c
extract_variables!(lin_constr::Dict{Int, Float64}, c::LinearityExpr) =
    extract_variables!(lin_constr, c.c)
function extract_variables!(lin_constr::Dict{Int, Float64}, c::Expr)
    if c.head == :ref
        if c.args[1] == :x
            @assert isa(c.args[2], Int)
            lin_constr[c.args[2]] = 0
        else
            error("Unrecognized reference expression $c")
        end
    else
        map(arg -> extract_variables!(lin_constr, arg), c.args)
    end
end

add_constant(c, constant::Real) = c + constant
add_constant(c::Expr, constant::Real) = Expr(:call, :+, c, constant)

function make_var_index!(m::AmplNLMathProgModel)
    nonlin_cont = Int[]
    nonlin_int = Int[]
    lin_cont = Int[]
    lin_int = Int[]
    lin_bin = Int[]

    for i in 1:m.nvar
        if m.varlinearities_obj[i] == :Nonlin ||
           m.varlinearities_con[i] == :Nonlin
            if m.vartypes[i] == :Cont
                push!(nonlin_cont, i)
            else
                push!(nonlin_int, i)
            end
        else
            if m.vartypes[i] == :Cont
                push!(lin_cont, i)
            elseif m.vartypes[i] == :Int
                push!(lin_int, i)
            else
                push!(lin_bin, i)
            end
        end
    end

    # Index variables in required order
    for var_list in (nonlin_cont, nonlin_int, lin_cont, lin_bin, lin_int)
        add_to_index_maps!(m.v_index_map, m.v_index_map_rev, var_list)
    end
end

function make_con_index!(m::AmplNLMathProgModel)
    nonlin_cons = Int[]
    lin_cons = Int[]

    for i in 1:m.ncon
        if m.conlinearities[i] == :Nonlin
            push!(nonlin_cons, i)
        else
            push!(lin_cons, i)
        end
    end
    for con_list in (nonlin_cons, lin_cons)
        add_to_index_maps!(m.c_index_map, m.c_index_map_rev, con_list)
    end
end

function add_to_index_maps!(forward_map::Dict{Int, Int},
                            backward_map::Dict{Int, Int},
                            inds::Array{Int})
    for i in inds
        # Indices are 0-prefixed so the next index is the current dict length
        index = length(forward_map)
        forward_map[i] = index
        backward_map[index] = i
    end
end

function read_results(m::AmplNLMathProgModel)
    did_read_solution = read_sol(m)

    # Convert solve_result
    if 0 <= m.solve_result_num < 100
        m.status = :Optimal
        m.solve_result = "solved"
    elseif 100 <= m.solve_result_num < 200
        # Used to indicate solution present but likely incorrect.
        m.status = :Optimal
        m.solve_result = "solved?"
        warn("The solver has returned the status :Optimal, but indicated that there might be an error in the solution. The status code returned by the solver was $(m.solve_result_num). Check the solver documentation for more info.""")
    elseif 200 <= m.solve_result_num < 300
        m.status = :Infeasible
        m.solve_result = "infeasible"
    elseif 300 <= m.solve_result_num < 400
        m.status = :Unbounded
        m.solve_result = "unbounded"
    elseif 400 <= m.solve_result_num < 500
        m.status = :UserLimit
        m.solve_result = "limit"
    elseif 500 <= m.solve_result_num < 600
        m.status = :Error
        m.solve_result = "failure"
    end

    # If we didn't get a valid solve_result_num, try to get the status from the
    # solve_message string.
    # Some solvers (e.g. SCIP) don't ever print the suffixes so we need this.
    if m.status == :NotSolved
        message = lowercase(m.solve_message)
        if contains(message, "optimal")
            m.status = :Optimal
        elseif contains(message, "infeasible")
            m.status = :Infeasible
        elseif contains(message, "unbounded")
            m.status = :Unbounded
        elseif contains(message, "limit")
            m.status = :UserLimit
        elseif contains(message, "error")
            m.status = :Error
        end
    end

    if did_read_solution
        # Calculate objective value from nonlinear and linear parts
        obj_nonlin = eval(substitute_vars!(deepcopy(m.obj), m.solution))
        obj_lin = evaluate_linear(m.lin_obj, m.solution)
        m.objval = obj_nonlin + obj_lin
    end
end

function read_sol(m::AmplNLMathProgModel)
    # Reference implementation:
    # https://github.com/ampl/mp/tree/master/src/asl/solvers/readsol.c

    f = open(m.solfile, "r")
    stat = :Undefined
    line = ""

    # Keep building solver message by reading until first truly empty line
    while true
        line = readline(f)
        isempty(chomp(line)) && break
        m.solve_message *= line
    end

    # Skip over empty lines
    while true
        line = readline(f)
        !isempty(chomp(line)) && break
    end

    # Read through all the options. Direct copy of reference implementation.
    @assert line[1:7] == "Options"
    options = [parse(Int, chomp(readline(f))) for _ in 1:3]
    num_options = options[1]
    3 <= num_options <= 9 || error("expected num_options between 3 and 9; " *
                                   "got $num_options")
    need_vbtol = false
    if options[3] == 3
        num_options -= 2
        need_vbtol = true
    end
    for j = 3:num_options
        eof(f) && error()
        push!(options, parse(Int, chomp(readline(f))))
    end

    # Read number of constraints
    num_cons = parse(Int, chomp(readline(f)))
    @assert(num_cons == m.ncon)

    # Read number of duals to read in
    num_duals_to_read = parse(Int, chomp(readline(f)))
    @assert(num_duals_to_read in [0; m.ncon])

    # Read number of variables
    num_vars = parse(Int, chomp(readline(f)))
    @assert(num_vars == m.nvar)

    # Read number of variables to read in
    num_vars_to_read = parse(Int, chomp(readline(f)))
    @assert(num_vars_to_read in [0; m.nvar])

    # Skip over vbtol line if present
    need_vbtol && readline(f)

    # Skip over duals
    # TODO do something with these?
    for index in 0:(num_duals_to_read - 1)
        eof(f) && error("End of file while reading duals.")
        line = readline(f)
    end

    # Next, read for the variable values
    x = fill(NaN, m.nvar)
    m.objval = NaN
    for index in 0:(num_vars_to_read - 1)
        eof(f) && error("End of file while reading variables.")
        line = readline(f)

        i = m.v_index_map_rev[index]
        x[i] = float(chomp(line))
    end
    m.solution = x

    # Check for status code
    while !eof(f)
        line = readline(f)
        linevals = split(chomp(line), " ")
        num_vals = length(linevals)
        if num_vals > 0 && linevals[1] == "objno"
            # Check for objno == 0
            @assert parse(Int, linevals[2]) == 0
            # Get solve_result
            m.solve_result_num = parse(Int, linevals[3])

            # We can stop looking for the 'objno' line
            break
        end
    end

    close(f)
    return num_vars_to_read > 0
end

substitute_vars!(c, x::Array{Float64}) = c
function substitute_vars!(c::Expr, x::Array{Float64})
    if c.head == :ref
        if c.args[1] == :x
            index = c.args[2]
            @assert isa(index, Int)
            c = x[index]
        else
            error("Unrecognized reference expression $c")
        end
    else
        if c.head == :call
            # Convert .nl unary minus (:neg) back to :-
            if c.args[1] == :neg
                c.args[1] = :-
            # Convert .nl :sum back to :+
            elseif c.args[1] == :sum
                c.args[1] = :+
            end
        end
        map!(arg -> substitute_vars!(arg, x), c.args)
    end
    c
end

function evaluate_linear(linear_coeffs::Dict{Int, Float64}, x::Array{Float64})
    total = 0.0
    for (i, coeff) in linear_coeffs
        total += coeff * x[i]
    end
    total
end

# Wrapper functions
for f in [:getvartype,:getsense,:optimize!,:status,:getsolution,:getobjval,:numvar,:numconstr,:get_solve_result,:get_solve_result_num,:get_solve_message,:get_solve_exitcode,:getsolvetime]
    @eval $f(m::AmplNLNonlinearModel) = $f(m.inner)
    @eval $f(m::AmplNLLinearQuadraticModel) = $f(m.inner)
end
for f in [:setvartype!,:setsense!,:setwarmstart!]
    @eval $f(m::AmplNLNonlinearModel, x) = $f(m.inner, x)
    @eval $f(m::AmplNLLinearQuadraticModel, x) = $f(m.inner, x)
end

# Utility method for deleting any leftover debug files
function clean_solverdata()
    for file in readdir(solverdata_dir)
        ext = splitext(file)[2]
        (ext == ".nl" || ext == ".sol") && rm(joinpath(solverdata_dir, file))
    end
end

end


text:
"""
Ensures that the LP model in `state.gurobi_model` is up-to-date for the current
`node`. Depending on how things are configured, this could require us doing
some or all of: rebuilding the model from scratch, applying new branchings,
and reformulating disjunctions.
"""
function populate_lp_model!(
    state::CurrentState,
    form::DMIPFormulation,
    node::Node,
    node_result::NodeResult,
    config::AlgorithmConfig,
)
    if state.rebuild_model
        create_base_model!(state, form, node, config)
        apply_branchings!(state, node)
        formulate_disjunctions!(state, form, node, node_result, config)
    else
        # If this case we would like to reuse the same model and not rebuild
        # from scratch...
        if !state.on_a_dive
            # ...however, upon backtracking we need to reset bounds and reapply
            # (or reset) disjunctive formulations.
            remove_all_branchings!(state, form)
        end
        apply_branchings!(state, node)
        if config.formulation_tightening_strategy == TIGHTEN_AT_EACH_NODE
            error(
                "Cerberus does not currently support tightening formulations via problem modification; you must rebuild from scratch.",
            )
        end
        return nothing
    end
    return nothing
end

# NOTE: At all times, the first `num_variables(form)` variables in the model
# will correspond to those registered with `form`. The disjunctive formulaters
# may add additional continuous variables, but they must come after this chunk.
# This means that these variables will present, in the same order, in every
# node LP.
function create_base_model!(
    state::CurrentState,
    form::DMIPFormulation,
    node::Node,
    config::AlgorithmConfig,
)
    reset_formulation_state!(state)
    model = config.lp_solver_factory(state, config)::Gurobi.Optimizer
    state.gurobi_model = model
    for cvi in all_variables(form)
        l, u = get_bounds(form, cvi)
        # Cache the above updates in formulation. Even better,
        # batch add variables.
        vi, ci = MOI.add_constrained_variable(model, IN(l, u))
        attach_index!(state, vi)
        push!(state.constraint_state.base_state.var_constrs, ci)
    end
    for lt_constr in get_constraints(form, CCI{LT})
        # Invariant: constraint was normalized via MOIU.normalize_constant.
        ci = MOI.add_constraint(
            model,
            instantiate(lt_constr.f, state),
            lt_constr.s,
        )
        push!(state.constraint_state.base_state.lt_constrs, ci)
    end
    for gt_constr in get_constraints(form, CCI{GT})
        # Invariant: constraint was normalized via MOIU.normalize_constant.
        ci = MOI.add_constraint(
            model,
            instantiate(gt_constr.f, state),
            gt_constr.s,
        )
        push!(state.constraint_state.base_state.gt_constrs, ci)
    end
    for et_constr in get_constraints(form, CCI{ET})
        # Invariant: constraint was normalized via MOIU.normalize_constant.
        ci = MOI.add_constraint(
            model,
            instantiate(et_constr.f, state),
            et_constr.s,
        )
        push!(state.constraint_state.base_state.et_constrs, ci)
    end
    MOI.set(model, MOI.ObjectiveFunction{SAF}(), instantiate(form.obj, state))
    MOI.set(model, MOI.ObjectiveSense(), MOI.MIN_SENSE)
    state.rebuild_model = false
    state.total_model_builds += 1
    return nothing
end

function formulate_disjunctions!(
    state::CurrentState,
    form::DMIPFormulation,
    node::Node,
    node_result::NodeResult,
    config::AlgorithmConfig,
)
    for (formulater, cvis) in form.disjunction_formulaters
        formulate!(
            state,
            form,
            formulater,
            # If we use a static formulation, it must be valid at the root. In
            # a bit of a hack, in this case we will just pass in an empty Node.
            if config.formulation_tightening_strategy == STATIC_FORMULATION
                Node()
            else
                node
            end,
            node_result,
            config,
        )
    end
    return nothing
end

function _unattached_bounds(cs::ConstraintState, node::Node, ::Type{LT})
    return view(
        node.lt_bounds,
        Colon()(cs.branch_state.num_lt_branches + 1, length(node.lt_bounds)),
    )
end
function _unattached_bounds(cs::ConstraintState, node::Node, ::Type{GT})
    return view(
        node.gt_bounds,
        Colon()(cs.branch_state.num_gt_branches + 1, length(node.gt_bounds)),
    )
end
function _unattached_constraints(cs::ConstraintState, node::Node, ::Type{LT})
    return view(
        node.lt_general_constrs,
        Colon()(
            length(cs.branch_state.lt_general_constrs) + 1,
            length(node.lt_general_constrs),
        ),
    )
end
function _unattached_constraints(cs::ConstraintState, node::Node, ::Type{GT})
    return view(
        node.gt_general_constrs,
        Colon()(
            length(cs.branch_state.gt_general_constrs) + 1,
            length(node.gt_general_constrs),
        ),
    )
end

function remove_all_branchings!(state::CurrentState, form::DMIPFormulation)
    model = state.gurobi_model
    cs = state.constraint_state
    for cvi in all_variables(form)
        l, u = get_bounds(form, cvi)
        ci = cs.base_state.var_constrs[index(cvi)]
        MOI.set(model, MOI.ConstraintSet(), ci, IN(l, u))
    end
    # TODO: Can potentially be smarter about not deleting all of these
    # constraints on a backtrack.
    MOI.delete(model, cs.branch_state.lt_general_constrs)
    MOI.delete(model, cs.branch_state.gt_general_constrs)
    empty!(cs.branch_state)
    return nothing
end

# We skip any bounds or constraints that are already added to the model. Note
# that we do this a bit dangerously--we merely count the number of
# bounds/constraints attached (as recorded in constraint_state). Therefore,
# this will only work on dives.
function apply_branchings!(state::CurrentState, node::Node)
    model = state.gurobi_model
    cs = state.constraint_state
    for lt_bound in _unattached_bounds(cs, node, LT)
        cvi = lt_bound.cvi
        vi = instantiate(cvi, state)
        ci = cs.base_state.var_constrs[index(cvi)]
        interval = MOI.get(model, MOI.ConstraintSet(), ci)
        new_interval = IN(interval.lower, min(lt_bound.s.upper, interval.upper))
        MOI.set(model, MOI.ConstraintSet(), ci, new_interval)
        cs.branch_state.num_lt_branches += 1
    end
    for gt_bound in _unattached_bounds(cs, node, GT)
        cvi = gt_bound.cvi
        vi = instantiate(cvi, state)
        ci = cs.base_state.var_constrs[index(cvi)]
        interval = MOI.get(model, MOI.ConstraintSet(), ci)
        new_interval = IN(max(gt_bound.s.lower, interval.lower), interval.upper)
        MOI.set(model, MOI.ConstraintSet(), ci, new_interval)
        cs.branch_state.num_gt_branches += 1
    end
    for ac in _unattached_constraints(cs, node, LT)
        # Invariant: constraint was normalized via MOIU.normalize_constant.
        ci = MOI.add_constraint(model, instantiate(ac.f, state), ac.s)
        push!(cs.branch_state.lt_general_constrs, ci)
    end
    for ac in _unattached_constraints(cs, node, GT)
        # Invariant: constraint was normalized via MOIU.normalize_constant.
        ci = MOI.add_constraint(model, instantiate(ac.f, state), ac.s)
        push!(cs.branch_state.gt_general_constrs, ci)
    end
    return nothing
end

# TODO: Unit test
function _update_lp_solution!(state::CurrentState, form::DMIPFormulation)
    # NOTE: I believe that the following is costless if the size is correct.
    resize!(state.current_solution, num_variables(form))
    for cvi in all_variables(form)
        state.current_solution[index(cvi)] = MOI.get(
            state.gurobi_model,
            MOI.VariablePrimal(),
            instantiate(cvi, state),
        )
    end
    return nothing
end

function get_basis(state::CurrentState)::Basis
    basis = Basis()
    cs = state.constraint_state
    for ci in cs.base_state.var_constrs
        push!(
            basis.base_var_constrs,
            MOI.get(state.gurobi_model, MOI.ConstraintBasisStatus(), ci),
        )
    end
    for ci in cs.base_state.lt_constrs
        push!(
            basis.base_lt_constrs,
            MOI.get(state.gurobi_model, MOI.ConstraintBasisStatus(), ci),
        )
    end
    for ci in cs.base_state.gt_constrs
        push!(
            basis.base_gt_constrs,
            MOI.get(state.gurobi_model, MOI.ConstraintBasisStatus(), ci),
        )
    end
    for ci in cs.base_state.et_constrs
        push!(
            basis.base_et_constrs,
            MOI.get(state.gurobi_model, MOI.ConstraintBasisStatus(), ci),
        )
    end
    for ci in cs.branch_state.lt_general_constrs
        push!(
            basis.branch_lt_constrs,
            MOI.get(state.gurobi_model, MOI.ConstraintBasisStatus(), ci),
        )
    end
    for ci in cs.branch_state.gt_general_constrs
        push!(
            basis.branch_gt_constrs,
            MOI.get(state.gurobi_model, MOI.ConstraintBasisStatus(), ci),
        )
    end
    return basis
end

function set_basis_if_available!(state::CurrentState, node::Node)
    if haskey(state.warm_starts, node)
        _set_basis!(
            state.gurobi_model,
            state.constraint_state,
            state.warm_starts[node],
        )
        state.total_warm_starts += 1
    end
    return nothing
end

function _set_basis!(
    model::MOI.AbstractOptimizer,
    cs::ConstraintState,
    basis::Basis,
)::Nothing
    # TODO: Check that basis is, in fact, a basis after modification
    @debug "Basis is being set ($(length(basis)) elements)"
    if isempty(basis.base_var_constrs) &&
       isempty(basis.base_lt_constrs) &&
       isempty(basis.base_gt_constrs) &&
       isempty(basis.base_et_constrs) &&
       isempty(basis.branch_lt_constrs) &&
       isempty(basis.branch_gt_constrs)
        throw(ArgumentError("You are attempting to set an empty basis."))
    end
    for (i, bs) in enumerate(basis.base_var_constrs)
        ci = cs.base_state.var_constrs[i]
        MOI.set(model, MOI.ConstraintBasisStatus(), ci, bs)
    end
    for (i, bs) in enumerate(basis.base_lt_constrs)
        ci = cs.base_state.lt_constrs[i]
        MOI.set(model, MOI.ConstraintBasisStatus(), ci, bs)
    end
    for (i, bs) in enumerate(basis.base_gt_constrs)
        ci = cs.base_state.gt_constrs[i]
        MOI.set(model, MOI.ConstraintBasisStatus(), ci, bs)
    end
    for (i, bs) in enumerate(basis.base_et_constrs)
        ci = cs.base_state.et_constrs[i]
        MOI.set(model, MOI.ConstraintBasisStatus(), ci, bs)
    end
    for (i, bs) in enumerate(basis.branch_lt_constrs)
        ci = cs.branch_state.lt_general_constrs[i]
        MOI.set(model, MOI.ConstraintBasisStatus(), ci, bs)
    end
    for (i, bs) in enumerate(basis.branch_gt_constrs)
        ci = cs.branch_state.gt_general_constrs[i]
        MOI.set(model, MOI.ConstraintBasisStatus(), ci, bs)
    end
    return nothing
end


text:
using Plots
using ProgressMeter

function eeg_animate(eeg, n, idx, skip=1; kwargs...)
    pr = ProgressMeter.Progress(Int64(floor((length(eeg.activations)-n)/skip)), 0.01, "Anim: ")
    anim = @animate for i in 1:skip:(length(eeg.activations)-n)
        l = @layout grid(length(eeg.activations[1]), 1)
        plt_a = plot(
            hcat(getindex.(eeg.activations, 1)...)',
            xlims=(i, n+i), title="Inputs", label="")
        plt_acts = [plot(
            hcat(getindex.(eeg.activations, idx)...)',
            xlims=(i,n+i), label="", title="Layer $(idx)") for idx in 2:length(eeg.activations[1])]
        plot(plt_a, plt_acts...; layout=l, kwargs...)
        next!(pr)
        # gui()
    end
    anim
end

function plot_eeg(eeg)

end


text:
using ProtoBuf
using Test
import ProtoBuf.meta

function testmetalock()
    current_async_safety = ProtoBuf.MetaLock.lck !== nothing
    for async_safety in (true, false)
        ProtoBuf.enable_async_safety(async_safety)
        str = ProtoBuf.metalock() do
            "test with lock enabled $async_safety"
        end
        @test endswith(str, string(async_safety))
    end
    ProtoBuf.enable_async_safety(current_async_safety)
end

@testset "Metadata locking" begin
    testmetalock()
end


text:
push!(LOAD_PATH,"../src/")
using Documenter, gmsh

include("generate.jl")

GENERATEDEXAMPLES = [joinpath("examples", f) for f in (
    "t1.md",
    "t2.md",
    "t3.md",
    "t4.md",
    "t5.md",
    "t10.md",
    #"t16.md"
    )]

makedocs(sitename="Gmsh API for Julia",
         modules=[gmsh],
         authors="Maximilian Köhler",
         pages=["Home"=> "index.md",
                "Tutorial" => GENERATEDEXAMPLES,
                "Built-in Geometry Kernel" => Any["Geometry" => "geo/geometry.md",
                                                  "Mesh Constraints" => "geo/mesh.md"],
                "Model" => Any["Model" => "model/model.md",
                               "Meshing" => "model/mesh.md"],
                "Options" => "options/options.md",
                "OpenCASCADE Kernel" => "occ/occ.md"])
deploydocs(
    repo = "github.com/koehlerson/gmsh.jl.git",
    push_preview=true,
)


text:
# SimLynx/src/process.jl
# Licensed under the MIT License. See LICENSE.md file in the project root for
# full license information.

"""
    Acceptor

Represents the function to accept a message to a process.

# Fields
- `name::Symbol`: the name of the message to be accepted
- `proc::Union{Function, Nothing}`: the function to be executed to accept a
      message from a process or nothing if there is no critical section
"""
struct Acceptor
    name::Symbol
    proc::Union{Function, Nothing}
end

"""
    Process

Represents a process in a simulation.

# Fields
- `name::String`: the name of the process (used in tracing)
- `state::Symbol`: the process state
- `task::Task`: the task executing the process
- `notice::Notice`: the notice used to schedule the process
- `queue::Array{AbstractMessage, 1}`: the queue of messages to the process
- `acceptors::Array{Acceptor, 1}`: the list of current acceptors for the process
- `response::Any`: the response to a message send to a process
- `storage::IdDict{Symbol, Any}`: an IdDict to store the process properties

# States
- `:created`: the process has been created but is not yet running
- `:active`: the process is running
- `:working`: the process is working or waiting
- `:delayed`: the process is delayed (waiting for some resource)
- `:interrupted`: the process has been interrupted by another process
- `:suspended`: the process has suspended itself
- `:terminated`: the process has terminated
"""
mutable struct Process <: ActiveElement
    name::String
    state::Symbol
    task::Task
    notice::Notice
    queue::Array{AbstractMessage, 1}
    acceptors::Array{Acceptor, 1}
    response::Any
    storage::IdDict{Symbol, Any}
    function Process(form::String, task::Task)
        x = new(form,
                :created,
                task)
        x.notice = Notice(0.0, x)
        x.queue = Array{AbstractMessage, 1}(undef, 0)
        x.acceptors = Array{Acceptor, 1}(undef, 0)
        x.response = nothing
        x.storage = IdDict{Symbol, Any}()
        return x
    end
end

"Return the value of the property key for the specified process."
process_store(process::Process, key::Symbol) = process.storage[key]

"Sets the value of the property key for the specified process."
function process_store!(process::Process, key::Symbol, value)
    process.storage[key] = value
end

"Prints the process name and state to the IO stream."
Base.show(io::IO, process::Process) =
    print(io, "Process $(process.name) [$(process.state)]")

"Return the process state for the specified process."
process_state(process::Process) = process.state

"Set the process state for a the specified process."
function process_state!(process::Process, state::Symbol)
    process.state = state
    if current_trace()
        @printf("%9.3f: %s\n", current_time(), process)
    end
end

"""
    @process <sig> begin
        <body>
    end

Define a simulation process with the specified signature and implemented by the
given body.
"""
macro process(sig, body)
    @capture(sig, f_Symbol(xs__)) ||
        throw(ArgumentError("the first argument must be a signature, " *
                            "given $sig"))
    @capture(body, begin exprs__ end) ||
        throw(ArgumentError("the second argument must be a body, " *
                            "given $body"))
    # Extract the argument identifiers
    args = [isa(arg, Symbol) ? arg : arg.args[1] for arg in xs]
    quote
        $(esc(sig)) =
            let _argvals = [$(esc.(args)...)] # Build array of argument values
                Process(string($(esc(f))) * "(" * join(string.(_argvals), ", ") * ")",
                        @task begin
                            try
                                $(esc(body))
                            catch e
                                showerror(stdout, e, catch_backtrace())
                            finally
                                process_state!(current_process(), :terminated)
                                yieldto(current_simulation.control_task)
                            end
                        end)
            end
    end
end


text:
module Optimise

export update!, params, train!,
  SGD, ADAM, Momentum, Nesterov, RMSProp, ADAGrad, ADADelta

struct Param{T}
  x::T
  Δ::T
end

Base.convert(::Type{Param}, x::AbstractArray) = Param(x, zeros(x))

include("optimisers.jl")
include("interface.jl")
include("train.jl")

using Flux.Tracker: TrackedArray

Base.convert(::Type{Param}, x::TrackedArray) = Param(x.data, x.grad)

end


text:
@testset "Test count word occurrence in file" begin
    msg = "one two one two three four"
    (fname, testio) = mktemp()
    println(testio, msg)
    close(testio)
    @test count_word_occurrence_in_file(fname, "one") == 2
    @test count_word_occurrence_in_file(fname, "three") == 1
    @test count_word_occurrence_in_file(fname, "six") == 0
end


text: Dict{String,Any}("ints" => Dict{String,Any}("value" => Any[Dict{String,Any}("value" => "1","type" => "integer"), Dict{String,Any}("value" => "2","type" => "integer"), Dict{String,Any}("value" => "3","type" => "integer")],"type" => "array"))

max_stars_repo_name: SmiteWindows/julia

max_issues_repo_name: SmiteWindows/julia

max_forks_repo_name: SmiteWindows/julia


text:
#################################################
#
#   Abstract types
#
#################################################

abstract Functor{N}  # N is the number of arguments

typealias UnaryFunctor Functor{1}
typealias BinaryFunctor Functor{2}
typealias TernaryFunctor Functor{3}

## macros for defining functors

# On 0.4 we can do `function evaluate end` (althought it's probably better)
# to use the call overload when we require julia 0.4
function evaluate() end

macro functor1(F, fun, T)
    eF = esc(F)
    efun = esc(fun)
    eT = esc(T)
    quote
        type $eF <: Functor{1} end
        NumericFuns.evaluate(::$eF, x::$eT) = $efun(x)
    end
end

macro functor2(F, fun, T)
    eF = esc(F)
    efun = esc(fun)
    eT = esc(T)
    quote
        type $eF <: Functor{2} end
        NumericFuns.evaluate(::$eF, x::$eT, y::$eT) = $efun(x, y)
    end
end

default_functorsym(f::Symbol) =
    (fstr = string(f); symbol(string(uppercase(fstr[1]), fstr[2:end], "Fun")))

macro functor1a(fun, T)
    F = default_functorsym(fun)
    eF = esc(F)
    efun = esc(fun)
    eT = esc(T)
    quote
        type $eF <: Functor{1} end
        NumericFuns.evaluate(::$eF, x::$eT) = $efun(x)
    end
end

macro functor2a(fun, T)
    F = default_functorsym(fun)
    eF = esc(F)
    efun = esc(fun)
    eT = esc(T)
    quote
        type $eF <: Functor{2} end
        NumericFuns.evaluate(::$eF, x::$eT, y::$eT) = $efun(x, y)
    end
end

macro functor1a_ord(fun, T)
    F = default_functorsym(fun)
    eF = esc(F)
    efun = esc(fun)
    eT = esc(T)
    quote
        immutable $eF{OT<:Real} <: Functor{1}
            order::OT
        end
        $eF{OT<:Real}(ord::OT) = $eF{OT}(ord)
        NumericFuns.evaluate(f::$eF, x::$eT) = $efun(f.order, x)
    end
end

#################################################
#
#   Functors
#
#################################################

## arithmetic operators

export Negate, Add, Subtract, Multiply, Divide, RDivide, Pow
export FixPow, FixAbsPow

@functor1(Negate,   -, Number)
@functor2(Add,      +, Number)
@functor2(Subtract, -, Number)
@functor2(Multiply, *, Number)
@functor2(Divide,   /, Number)
@functor2(RDivide,  \, Number)
@functor2(Pow,      ^, Number)

immutable FixPow{T<:Real} <: Functor{1}
    p::T
end
evaluate(f::FixPow, x::Real) = x ^ f.p

immutable FixAbsPow{T<:Real} <: Functor{1}
    p::T
end
evaluate(f::FixAbsPow, x::Real) = abs(x) ^ f.p

## comparison operators

export LT, GT, LE, GE, EQ, NE

@functor2(LT, <,  Real)
@functor2(GT, >,  Real)
@functor2(LE, <=, Real)
@functor2(GE, >=, Real)
@functor2(EQ, ==, Number)
@functor2(NE, !=, Number)

## logical & bitwise operators

export Not, And, Or
export BitwiseNot, BitwiseAnd, BitwiseOr, BitwiseXor

@functor1(Not, !, Bool)
@functor2(And, &, Bool)
@functor2(Or,  |, Bool)

@functor1(BitwiseNot, ~, Integer)
@functor2(BitwiseAnd, &, Integer)
@functor2(BitwiseOr,  |, Integer)
@functor2(BitwiseXor, $, Integer)

## arithmetic functions

export DivFun, FldFun, RemFun, ModFun
export AbsFun, Abs2Fun, RealFun, ImagFun, MaxFun, MinFun
export SqrFun, RcpFun, SignFun, SignbitFun

@functor2a div Real
@functor2a fld Real
@functor2a rem Real
@functor2a mod Real

@functor1a abs  Number
@functor1a abs2 Number
@functor1a real Number
@functor1a imag Number
@functor2a max Real
@functor2a min Real

@functor1a sqr Number
@functor1a rcp Number

@functor1a sign     Real
@functor1a signbit  Real

## rounding functions

export FloorFun, CeilFun, TruncFun, RoundFun
export IfloorFun, IceilFun, ItruncFun, IroundFun

@functor1a floor Real
@functor1a ceil  Real
@functor1a trunc Real
@functor1a round Real

_ifloor(x) = floor(Integer, x)
_iceil(x) = ceil(Integer, x)
_itrunc(x) = trunc(Integer, x)
_iround(x) = round(Integer, x)

@functor1 IfloorFun _ifloor Real
@functor1 IceilFun  _iceil  Real
@functor1 ItruncFun _itrunc Real
@functor1 IroundFun _iround Real

## number classification

export IsnanFun, IsinfFun, IsfiniteFun

@functor1a isnan Real
@functor1a isinf Real
@functor1a isfinite Real

## algebraic functions

export SqrtFun, CbrtFun, HypotFun
export RsqrtFun, RcbrtFun

@functor1a sqrt  Number
@functor1a cbrt  Real
@functor2a hypot Real

@functor1a rsqrt Number
@functor1a rcbrt Real

## exponential & logarithm

export ExpFun, Exp2Fun, Exp10Fun, Expm1Fun
export LogFun, Log2Fun, Log10Fun, Log1pFun
export XlogxFun, XlogyFun, LogisticFun, LogitFun
export SoftplusFun, InvsoftplusFun, LogsumexpFun

@functor1a exp   Number
@functor1a exp2  Number
@functor1a exp10 Number
@functor1a expm1 Real

@functor1a log   Number
@functor1a log2  Number
@functor1a log10 Number
@functor1a log1p Real

@functor1a xlogx Real
@functor2a xlogy Real
@functor1a logistic Real
@functor1a logit Real
@functor1a softplus Real
@functor1a invsoftplus Real
@functor2a logsumexp Real

## trigonometric functions

export SinFun, CosFun, TanFun, CotFun, SecFun, CscFun
export AsinFun, AcosFun, AtanFun, AcotFun, AsecFun, AcscFun, Atan2Fun
export SincFun, CoscFun, SinpiFun, CospiFun
export SindFun, CosdFun, TandFun, CotdFun, SecdFun, CscdFun
export AsindFun, AcosdFun, AtandFun, AcotdFun, AsecdFun, AcscdFun

@functor1a sin Number
@functor1a cos Number
@functor1a tan Number
@functor1a cot Number
@functor1a sec Number
@functor1a csc Number

@functor1a asin Number
@functor1a acos Number
@functor1a atan Number
@functor1a acot Number
@functor1a asec Number
@functor1a acsc Number
@functor2a atan2 Real

@functor1a sinc Number
@functor1a cosc Number
@functor1a sinpi Number
@functor1a cospi Number

@functor1a sind Real
@functor1a cosd Real
@functor1a tand Real
@functor1a cotd Real
@functor1a secd Real
@functor1a cscd Real

@functor1a asind Real
@functor1a acosd Real
@functor1a atand Real
@functor1a acotd Real
@functor1a asecd Real
@functor1a acscd Real

## hyperbolic functions

export SinhFun, CoshFun, TanhFun, CothFun, SechFun, CschFun
export AsinhFun, AcoshFun, AtanhFun, AcothFun, AsechFun, AcschFun

@functor1a sinh Number
@functor1a cosh Number
@functor1a tanh Number
@functor1a coth Number
@functor1a sech Number
@functor1a csch Number

@functor1a asinh Number
@functor1a acosh Number
@functor1a atanh Number
@functor1a acoth Number
@functor1a asech Number
@functor1a acsch Number

## special functions

export ErfFun, ErfcFun, ErfinvFun, ErfcinvFun, ErfiFun, ErfcxFun
export GammaFun, LgammaFun, DigammaFun, EtaFun, ZetaFun, BetaFun, LbetaFun
export AiryFun, AiryprimeFun, AiryaiFun, AiryaiprimeFun, AirybiFun, AirybiprimeFun

export Besselj0Fun, Besselj1Fun, Bessely0Fun, Bessely1Fun
export BesseliFun, BesseljFun, BesselkFun, BesselyFun
export Hankelh1Fun, Hankelh2Fun

@functor1a erf Number
@functor1a erfc Number
@functor1a erfinv Real
@functor1a erfcinv Real

@functor1a erfi Number
@functor1a erfcx Number

@functor1a gamma Number
@functor1a lgamma Number
@functor1a digamma Real

@functor1a eta Number
@functor1a zeta Number

@functor2a beta Real
@functor2a lbeta Real

@functor1a airy Number
@functor1a airyprime Number
@functor1a airyai Number
@functor1a airyaiprime Number
@functor1a airybi Number
@functor1a airybiprime Number

@functor1a besselj0 Number
@functor1a besselj1 Number
@functor1a bessely0 Number
@functor1a bessely1 Number

@functor1a_ord besseli Number
@functor1a_ord besselj Number
@functor1a_ord besselk Number
@functor1a_ord bessely Number

@functor1a_ord hankelh1 Number
@functor1a_ord hankelh2 Number

#######################################
#
#  Ternary functors
#
#######################################

export FMA, IfelseFun

type FMA <: Functor{3} end
evaluate(::FMA, x::Number, y::Number, z::Number) = (x + y * z)

type IfelseFun <: Functor{3} end
evaluate{T<:Number}(::IfelseFun, c::Bool, x::T, y::T) = ifelse(c, x, y)


text:
using CellListMap
using StaticArrays
import Random

#
# In this example we compute the "gravitational potential", assigning to each particle
# has a different mass. In this case, the closure is used to pass the masses to the
# function that computes the potential.
#
function gravitational_potential(;N=100_000,parallel=true,x=nothing)

    # Number of particles, sides and cutoff
    sides = @SVector [250,250,250]
    cutoff = 10.
    box = Box(sides, cutoff)

    # Particle positions
    Random.seed!(321)
    if x === nothing
        x = [ sides .* rand(SVector{3,Float64}) for i in 1:N ]
    end

    # Initialize auxiliary linked lists
    cl = CellList(x, box, parallel=parallel)

    # masses
    mass = [ 5 * x[i][1] for i in 1:N ]

    # Function to be evalulated for each pair: build distance histogram
    function potential(i, j, d2, u, mass)
        d = sqrt(d2)
        u = u - 9.8 * mass[i] * mass[j] / d
        return u
    end

    # Run pairwise computation
    u = map_pairwise(
        (x, y, i, j, d2, u) -> potential(i, j, d2, u, mass),
        0.0,box,cl,
        parallel=parallel
    )

    return u

end


text:
struct Interaction{F}
    func::F
    func_type::String
end

function func_type end
func_type(I::Interaction{F}) where {F} = I.func_type
Interaction(func) = Interaction{typeof(func)}(func, func_type(func))
@inline (I::Interaction)(v::Vector{T}) where T<:Real = I.func(v)

function nearest_neighbours(v::Vector{T}) where T<:Real
    margin = 1e-9
    if abs(norm(v)-1.) < margin
        return 1.
    else
        return 0.
    end
end

struct MDI{D}
    dir::NTuple{D, Int64}
    ndir::Vector{Float64}
    function MDI(dir::NTuple{D, Int64}) where {D}
        ndir = collect(dir)/norm(collect(dir))
        new{D}(dir, ndir)
    end
end
@inline MDI(dir...) = MDI(dir)

function (mdi::MDI)(v::Vector{T}) where T<:Real
    cos_theta = dot(mdi.ndir, v)/norm(v)
    return (1-3*cos_theta^2)/norm(v)^3
end

function func_type(Func)
    if isequal(Func, nearest_neighbours)
        return "(nn)"
    end
    if typeof(Func)<:MDI
        return "(md)$(Func.dir)"
    else
        error("no string for this function\n")
    end
end


text:
using Clang
using Clang.Generators
using SDL2_jll
using SDL2_mixer_jll
using SDL2_image_jll
using SDL2_ttf_jll

cd(@__DIR__)

sdl2_include_dir = joinpath(SDL2_jll.artifact_dir, "include") |> normpath
sdl_mixer_h = joinpath(SDL2_mixer_jll.artifact_dir, "include", "SDL2", "SDL_mixer.h") |> normpath
sdl_image_h = joinpath(SDL2_image_jll.artifact_dir, "include", "SDL2", "SDL_image.h") |> normpath
sdl_ttf_h = joinpath(SDL2_ttf_jll.artifact_dir, "include", "SDL2", "SDL_ttf.h") |> normpath
sdl_h = joinpath(sdl2_include_dir, "SDL2", "SDL.h")

local_include_dir = joinpath(@__DIR__, "include")
isdir(local_include_dir) && rm(local_include_dir, recursive=true)
cp(sdl2_include_dir, local_include_dir, force=true)
cp(sdl_mixer_h, joinpath(@__DIR__, "include", "SDL2", basename(sdl_mixer_h)))
cp(sdl_image_h, joinpath(@__DIR__, "include", "SDL2", basename(sdl_image_h)))
cp(sdl_ttf_h, joinpath(@__DIR__, "include", "SDL2", basename(sdl_ttf_h)))

sdl_mixer_h = joinpath(local_include_dir, "SDL2", "SDL_mixer.h")
sdl_image_h = joinpath(local_include_dir, "SDL2", "SDL_image.h")
sdl_ttf_h = joinpath(local_include_dir, "SDL2", "SDL_ttf.h")
sdl_h = joinpath(local_include_dir, "SDL2", "SDL.h")

options = load_options(joinpath(@__DIR__, "generator.toml"))

args = get_default_args()
push!(args, "-I$local_include_dir")

ctx = create_context([sdl_h, sdl_mixer_h, sdl_image_h, sdl_ttf_h], args, options)

build!(ctx)


text:
module FusionSystems

  using Roots
  using SymEngine

  using SimplePlots

  using ImportMacros
  @import SpecialFunctions.gamma as Gamma

  const A = 2.5 # D-T Fuel [AMU's]

  const H_max = 4
  const I_M_max = 25
  const n_20_max = 4
  const nu_J_max = 2
  const T_k_max = 40
  const B_C_max = 40

  const checked_min = 0.1
  const nu_J_min = -1/2

  const R_0_min = 1.0
  const R_0_max = 10.0

  const mu_0 = 4 * pi * 1e-7

  include("coeff.jl")
  include("bisection_indices.jl")

  include("structs/index.jl")
  include("decks/index.jl")
  include("coeffs/index.jl")
  include("solvers/index.jl")

end


text:
import Base.+

function +(buffer::Array{String,1},content::String)
    push!(buffer,content)
end

function +(buffer::Array{String,1}, content_array::Array{String,1}; padding::Union{String,Nothing}=nothing)
    for line in content_array

        if (padding === nothing)
            push!(buffer, line)
        else
            push!(buffer,"$(padding)$(line)")
        end
    end
end

max_stars_repo_name: varnerlab/CFMG

max_issues_repo_name: varnerlab/CFMG

max_forks_repo_name: varnerlab/CFMG


text:
module RemoteREPL

export connect_repl, serve_repl

# Technically, server and client could be completely separate packages, but
# having them together seems simplest for now.

const protocol_magic = "RemoteREPL"
const protocol_version = UInt32(1);

include("repl_server.jl")
include("repl_client.jl")

end


text:
## sigmoidal functions

export logistic, logistic_prime, weibull, weibull_prime

logistic(x) = 1/(1 + exp(-x))
logistic_prime(x) = exp(-x)/(1 + exp(-x))^2
weibull(x, coef, scale, shape) = (1 - exp(-((x/scale)^shape)))*coef
weibull_prime(x, coef, scale, shape) = coef*(shape/scale)*((x/scale)^(shape-1))*exp(-((x/scale)^shape))


text:
CURRENT_YEAR=2017

# Define shapes
abstract Shape
type Position
     x::Real
     y::Real
end
type Circ <: Shape
     center::Position
     radius::Real
end
type Square <: Shape
     upper_left::Position
     length::Real
end
type Rect <: Shape
     upper_left::Position
     width::Real
     height::Real
end

# Define area function
function area(shape::Rect)
	 shape.width * shape.height
end
function area(shape::Circ)
	 return (pi * (shape.radius * shape.radius))
end
function area(shape::Square)
	 return (shape.length * shape.length)
end

# Define in_shape function
function in_shape(shape::Square, position::Position)
	 !(position.x > shape.upper_left.x ||
	 	      position.x < shape.upper_left.x + shape.length ||
		      position.y > shape.upper_left.y ||
		      position.y < shape.upper_left.y - shape.length)
end
function in_shape(shape::Circ, position::Position)
	 ((position.x - shape.center.x)^2 + (position.y - shape.center.y)^2) <= shape.radius^2
end
function in_shape(shape::Rect, position::Position)
	 !(position.x < shape.upper_left.x ||
	 	      position.x > shape.upper_left.x + shape.width ||
		      position.y > shape.upper_left.y ||
		      position.y < shape.upper_left.y - shape.height)
end

# Define 'Pixel'
type Pixel
     r::Real
     b::Real
     g::Real
end

# Define greyscale function
function greyscale(picture::Array{Pixel,2})
	 map(x -> Pixel((x.r+x.b+x.g)/3, (x.r+x.b+x.g)/3, (x.r+x.b+x.g)/3), picture)
end

# Define invert function
function invert(picture::Array{Pixel,2})
	 map(x -> Pixel(255 - x.r, 255 - x.b, 255 - x.g), picture)
end

# Define data trees
abstract TreeItem
type Person <: TreeItem
     name::AbstractString
     birthyear::Integer
     eyecolor::Symbol
     father::TreeItem
     mother::TreeItem
end
type Unknown <: TreeItem
end

# Define count_persons function
function count_persons(tree)
	 isa(tree, Person) ? 1 + count_persons(tree.mother) + count_persons(tree.father) : 0
end

# Define average_age function
function total_age(tree)
	 isa(tree, Person) ? (CURRENT_YEAR - tree.birthyear) + total_age(tree.mother) + total_age(tree.father) : 0
end
function average_age(tree)
	 isa(tree, Person) ? total_age(tree) / count_persons(tree) : 0
end

# Define tree_map function
function tree_map(f, tree)
	 if isa(tree, Person)
	    treeitem=Person(tree.name, tree.birthyear, tree.eyecolor, tree.father, tree.mother)
	    treeitem.mother = tree_map(f,treeitem.mother)
	    treeitem.father = tree_map(f,treeitem.father)
	    afterfunc=f(treeitem)
	    afterfunc == nothing ? treeitem : afterfunc
	 else
	    Unknown()
	 end
end

# Define add_last_name function
function add_last_name(name::AbstractString, tree)
	 tree_map(x -> Person(x.name * name, x.birthyear, x.eyecolor, x.father, x.mother), tree)
end

# Define eye_colors function
function eye_colors(tree)
	 isa(tree, Person) ? [tree.eyecolor; eye_colors(tree.father); eye_colors(tree.mother)] : []
end

#=
# Begin test cases

# Test shape constructors
rect = Rect(Position(1, 1), 3, 4)
circ = Circ(Position(1, 1), 2)
square = Square(Position(1, 1), 2)

# Test area function
println(area(rect))
println(area(circ))
println(area(square))
println(in_shape(square, Position(0,2)))

# Test pixel constructors
px1=Pixel(128,128,128)
px2=Pixel(0,255,0)
px3=Pixel(255,0,255)
px4=Pixel(0,128,0)
picture=[px1 px2 ; px3 px4]
println(picture)
# Test greyscale function
println(greyscale(picture))
# Test invert function
println(invert(picture))
# Test tree constructrs
person1=Person("name1",1992,:C1,Unknown(),Unknown())
person2=Person("name2",1993,:C2,Unknown(),Unknown())
person3=Person("name3",1994,:C3,person1,person2)
person4=Person("name4",1995,:C4,person3,Unknown())
# Test count_persons function
println(count_persons(Unknown()))
println(count_persons(person1))
println(count_persons(person2))
println(count_persons(person3))
println(count_persons(person4))
# Test average_age function
println(average_age(Unknown()))
println(average_age(person1))
println(average_age(person2))
println(average_age(person3))
println(average_age(person4))
# Test tree_map function
println("unkown:")
tree_map(println, Unknown())
println("person1:")
tree_map(println, person1)
println("person2:")
tree_map(println, person2)
println("person3:")
tree_map(println, person3)
println("person4:")
tree_map(println, person4)
# Test add_last_name function
println(add_last_name(" last",Unknown()))
println(add_last_name(" last",person1))
println(add_last_name(" last",person2))
println(add_last_name(" last",person3))
println(add_last_name(" last",person4))
# Test eye_colors function
println(eye_colors(Unknown()))
println(eye_colors(person1))
println(eye_colors(person2))
println(eye_colors(person3))
println(eye_colors(person4))
# Finish test cases
println("Finished test cases\n\n\n\n\n")
=#

max_stars_repo_name: scottcrossen/JuliaDevelopment

max_issues_repo_name: scottcrossen/JuliaDevelopment

max_forks_repo_name: scottcrossen/JuliaDevelopment


text:
#The utility functions
#Lanfeng Pan
#Oct 29, 2014

function stopRule(pa::Vector, pa_old::Vector; tol=.005)
    maximum(abs.(pa .- pa_old)./(abs.(pa).+.001)) < tol
end

function marginallikelihood(X::Matrix{Float64}, Y::AbstractArray{Bool, 1}, groupindex::IntegerVector, wi::Vector{Float64}, mu::Vector{Float64}, sigmas::Vector{Float64}, β::Array{Float64,1}; ngh::Int=100)
    ghx, ghw = gausshermite(ngh)
    N, J = size(X)
    M = ngh * length(wi)
    n = maximum(groupindex)
    marginallikelihood(β, X, Y, groupindex, n, wi, mu, sigmas, ghx, ghw, zeros(N), zeros(n), zeros(N), zeros(n, M))
end

function marginallikelihood(beta_new::Array{Float64,1}, X::Matrix{Float64}, Y::AbstractArray{Bool, 1}, groupindex::IntegerVector, nF::Int64, wi::Vector{Float64}, mu::Vector{Float64}, sigmas::Vector{Float64}, ghx::Vector{Float64}, ghw::Vector{Float64}, llvec::Vector, ll_nF::Vector, xb::Vector, sumlogmat::Matrix)

    N,J = size(X)
    M = length(ghx)
    C = length(wi)
    # xb = X*beta_new
    mul!(xb, X, beta_new)
    ll = 0.0
    #ll_nF=zeros(nF)
    # sumlog_nF = zeros(nF)
    #llvec = zeros(N)
    # sumlogmat = zeros(nF, M*C)
    #fill!(sumlogmat, 0.0)
    for jcom in 1:C
        for ix in 1:M
            # fill!(llvec, ghx[ix]*sigmas[jcom]*sqrt(2)+mu[jcom])
            # add!(llvec, llvec, xb)
            # negateiftrue!(llvec, Y)
            #partially devecterize, speed up by 5%; But exp and log from Yeppp are much faster than devecterized ones.
            xtmp = ghx[ix]*sigmas[jcom]*sqrt(2)+mu[jcom]
            wtmp = log(ghw[ix])+log(wi[jcom])
            for i in 1:N
                @inbounds llvec[i] = ifelse(Y[i], -xtmp - xb[i], xtmp + xb[i])
            end
            # Yeppp.exp!(llvec, llvec)
            # log1p!(llvec)
            log1pexp!(llvec, llvec, N)
            ixM = ix+M*(jcom-1)
            for i in 1:nF
                sumlogmat[i, ixM] = wtmp
            end
            for i in 1:N
                @inbounds sumlogmat[groupindex[i], ixM] -= llvec[i]
            end
        end
    end
    for i in 1:nF
        ll += logsumexp(sumlogmat[i,:])
    end

    ll - nF*log(pi)/2
end

function predictgamma(X::Matrix, Y::Vector{Bool}, groupindex::IntegerVector, wi::Vector, mu::Vector, sigmas::Vector, β::Vector; ngh::Int=100)

    ncomponent = length(wi)
    n = maximum(groupindex)
    M = ngh*ncomponent
    N, J = size(X)

    xb = zeros(N)
    llN = zeros(N)
    llN2 = zeros(N)
    gammaM = zeros(M)
    gammahat = zeros(n)

    ghx, ghw = gausshermite(ngh)
    Wim = zeros(n, M)
    for ix in 1:ngh, jcom in 1:ncomponent
        ixM = ix+ngh*(jcom-1)
        gammaM[ixM] = ghx[ix]*sigmas[jcom]*sqrt(2)+mu[jcom]
    end

    mul!(xb, X, β)
    integralweight!(Wim, X, Y, groupindex, gammaM, wi, ghw, llN, llN2, xb, N, J, n, ncomponent, ngh)
    for i in 1:n
        for j in 1:M
            gammahat[i] += gammaM[j] * Wim[i,j]
        end
    end

    return gammahat
end
function FDR(X::Matrix, Y::Vector{Bool}, groupindex::IntegerVector, wi::Vector, mu::Vector, sigmas::Vector, β::Vector, C0::IntegerVector; ngh::Int=100, alphalevel::Real=.05)

    ncomponent = length(wi)
    n = maximum(groupindex)
    M = ngh*ncomponent
    N, J = size(X)
    if maximum(C0) > ncomponent || minimum(C0) < 1
        error("C0 must be within 1 to number of components")
    end

    xb = zeros(N)
    llN = zeros(N)
    llN2 = zeros(N)
    gammaM = zeros(M)
    piposterior = zeros(n, ncomponent)
    clFDR = zeros(n)

    ghx, ghw = gausshermite(ngh)
    Wim = zeros(n, M)
    for ix in 1:ngh, jcom in 1:ncomponent
        ixM = ix+ngh*(jcom-1)
        gammaM[ixM] = ghx[ix]*sigmas[jcom]*sqrt(2)+mu[jcom]
    end

    mul!(xb, X, β)
    integralweight!(Wim, X, Y, groupindex, gammaM, wi, ghw, llN, llN2, xb, N, J, n, ncomponent, ngh)
    for i in 1:n
        for jcom in 1:ncomponent
            for ix in 1:ngh
                ixM = ix+ngh*(jcom-1)
                piposterior[i, jcom] += Wim[i,ixM]
            end
        end
    end
    for i in 1:n
        for jcom in 1:ncomponent
            clFDR[i] = sum(piposterior[i, C0])
        end
    end
    order = sortperm(clFDR)
    n0 = 0
    for i in 1:n
        if mean(clFDR[order[1:i]]) < alphalevel
            n0 = i
        end
    end
    if n0 == 0
        println("Nothing is rejected")
    else
        println("The rejected groups are: ", order[1:n0])
        println("with FDR values: ", round.(clFDR[order[1:n0]], 4))
    end
    return clFDR, order[1:n0]
end
function asymptoticdistribution(X::Matrix{Float64}, Y::AbstractArray{Bool, 1}, groupindex::IntegerVector, wi::Vector{Float64}, mu::Vector{Float64}, sigmas::Vector{Float64}, betas::Array{Float64,1}; ngh::Int=100, nrep::Int=10000, debuginfo::Bool=false)

    N,J = size(X)
    nF = maximum(groupindex)
    M = ngh
    C = length(wi)
    if C == 1
        return rand(Chisq(2), nrep)
    end
    ghx, ghw = gausshermite(ngh)
    xb = zeros(N)
    mul!(xb, X, betas)
    llvec = zeros(N)
    llN2 = zeros(N)
    ll_nF = zeros(nF, C)
    sumlogmat = zeros(nF, ngh*C)
    summat_beta = zeros(nF, ngh*C, J)
    S_β = zeros(nF, J)
    S_π = zeros(nF, C-1)
    S_μσ = zeros(nF, 2*C)
    S_λ = zeros(nF, 2*C)
    ml = zeros(nF)
    xtmp = zeros(C*M)
    for jcom in 1:C
        for ix in 1:M
            ixM = ix+M*(jcom-1)
            xtmp[ixM] = ghx[ix]*sigmas[jcom]*sqrt(2)+mu[jcom]
            for i in 1:N
                @inbounds llvec[i] = ifelse(Y[i], -xtmp[ixM] - xb[i], xtmp[ixM] + xb[i])
            end

            # Yeppp.exp!(llvec, llvec)
            copy!(llN2, llvec)
            logistic!(llN2)
            negateiffalse!(llN2, Y)
            for i in 1:N
                for j in 1:J
                    summat_beta[groupindex[i], ixM, j] += llN2[i] * X[i,j]
                end
            end
            log1pexp!(llvec)

            for i in 1:N
                @inbounds sumlogmat[groupindex[i], ixM] -= llvec[i]
            end
            for i in 1:nF
                sumlogmat[i, ixM] +=  log(ghw[ix])
            end
        end
    end
    for i in 1:nF
        u = maximum(sumlogmat[i, :])
        for jcol in 1:C*M
            @inbounds sumlogmat[i, jcol] = sumlogmat[i, jcol] - u
        end
    end
    for i in 1:nF
        for kcom in 1:C
            ll_nF[i, kcom] = sumexp(sumlogmat[i,(1+M*(kcom-1)):M*kcom])
        end
    end
    for i in 1:nF
        for jcom in 1:C
            for ix in 1:M
                sumlogmat[i, ix+M*(jcom-1)] += log(wi[jcom])
            end
        end
        ml[i]=sumexp(sumlogmat[i, :])
    end
    for kcom in 1:(C-1)
        S_π[:, kcom] = (ll_nF[:, kcom] .- ll_nF[:, C]) ./ ml
    end
    for i in 1:nF
        for kcom in 1:C
            ind = (1+M*(kcom-1)):M*kcom
            # w[i, kcom] = ll_nF[i, kcom] * wi[kcom] / ml[i]

            S_μσ[i, 2*kcom-1] = sumexp(sumlogmat[i, ind], H1(xtmp[ind], mu[kcom], sigmas[kcom])) / ml[i]
            S_μσ[i, 2*kcom] = sumexp(sumlogmat[i, ind], H2(xtmp[ind], mu[kcom], sigmas[kcom]))/ml[i]
            S_λ[i, 2*kcom-1] = sumexp(sumlogmat[i, ind], H3(xtmp[ind], mu[kcom], sigmas[kcom]))/ml[i]
            S_λ[i, 2*kcom] = sumexp(sumlogmat[i, ind], H4(xtmp[ind], mu[kcom], sigmas[kcom]))/ml[i]
        end
        for j in 1:J
            S_β[i, j] = sumexp(sumlogmat[i,:], summat_beta[i, :, j])/ml[i]
        end
    end
    S_η = hcat(S_β, S_π, S_μσ)
    debuginfo && println(round.(sum(S_η, 1)./sqrt(nF), 6))
    I_η = S_η'*S_η./nF
    I_λη = S_λ'*S_η./nF
    I_λ = S_λ'*S_λ./nF
    I_all = vcat(hcat(I_η, I_λη'), hcat(I_λη, I_λ))
    if 1/cond(I_all) < eps(Float64)
        @warn("Information Matrix is singular!")
        D, V = eigen(I_all)
        debuginfo && println(D)
        tol2 = maximum(abs.(D)) * 1e-14
        D[D.<tol2] .= tol2
        I_all = V*Diagonal(D)*V'
    end
    debuginfo && println(round.(I_all, 6))
    I_λ_η = I_all[(J+3*C):(J+5*C-1), (J+3*C):(J+5*C-1)] - I_all[(J+3*C):(J+5*C-1), 1:(J+3*C-1)] * inv(I_all[1:(J+3*C-1), 1:(J+3*C-1)]) * I_all[1:(J+3*C-1),(J+3*C):(J+5*C-1) ]
    debuginfo && println(round.(I_λ_η, 6))
    D, V = eigen(I_λ_η)
    D[D.<0.] .= 0.
    I_λ_η2 = V * Diagonal(sqrt.(D)) * V'
    u = randn(nrep, 2*C) * I_λ_η2
    EM = zeros(nrep, C)
    T = zeros(nrep)
    for kcom in 1:C
        EM[:, kcom] = sum(u[:, (2*kcom-1):(2*kcom)] * inv(I_λ_η[(2*kcom-1):(2*kcom), (2*kcom-1):(2*kcom)]) .* u[:, (2*kcom-1):(2*kcom)], dims=2)
    end
    for i in 1:nrep
        T[i] = maximum(EM[i, :])
    end
    T
    #sum(ll_nF) - nF*log(pi)/2
end
function vcov(X::Matrix{Float64}, Y::AbstractArray{Bool, 1}, groupindex::IntegerVector, wi::Vector{Float64}, mu::Vector{Float64}, sigmas::Vector{Float64}, betas::Array{Float64,1}; ngh::Int=100, debuginfo::Bool=false)

    N,J = size(X)
    nF = maximum(groupindex)
    M = ngh
    C = length(wi)
    ghx, ghw = gausshermite(ngh)
    xb = zeros(N)
    mul!(xb, X, betas)
    llvec = zeros(N)
    llN2 = zeros(N)
    ll_nF = zeros(nF, C)
    sumlogmat = zeros(nF, ngh*C)
    summat_beta = zeros(nF, ngh*C, J)
    S_β = zeros(nF, J)
    S_π = zeros(nF, C-1)
    S_μσ = zeros(nF, 2*C)

    ml = zeros(nF)
    xtmp = zeros(C*M)
    for jcom in 1:C
        for ix in 1:M
            ixM = ix+M*(jcom-1)
            xtmp[ixM] = ghx[ix]*sigmas[jcom]*sqrt(2)+mu[jcom]
            for i in 1:N
                @inbounds llvec[i] = ifelse(Y[i], -xtmp[ixM] - xb[i], xtmp[ixM] + xb[i])
            end

            # Yeppp.exp!(llvec, llvec)
            copy!(llN2, llvec)
            logistic!(llN2)
            negateiffalse!(llN2, Y)
            for i in 1:N
                for j in 1:J
                    summat_beta[groupindex[i], ixM, j] += llN2[i] * X[i,j]
                end
            end
            log1pexp!(llvec)

            for i in 1:N
                @inbounds sumlogmat[groupindex[i], ixM] -= llvec[i]
            end
            for i in 1:nF
                sumlogmat[i, ixM] +=  log(ghw[ix])
            end
        end
    end
    for i in 1:nF
        u = maximum(sumlogmat[i, :])
        for jcol in 1:C*M
            @inbounds sumlogmat[i, jcol] = sumlogmat[i, jcol] - u
        end
    end
    for i in 1:nF
        for kcom in 1:C
            ll_nF[i, kcom] = sumexp(sumlogmat[i,(1+M*(kcom-1)):M*kcom])
        end
    end
    for i in 1:nF
        for jcom in 1:C
            for ix in 1:M
                sumlogmat[i, ix+M*(jcom-1)] += log(wi[jcom])
            end
        end
        ml[i]=sumexp(sumlogmat[i, :])
    end
    for kcom in 1:(C-1)
        S_π[:, kcom] = (ll_nF[:, kcom] .- ll_nF[:, C]) ./ ml
    end
    for i in 1:nF
        for kcom in 1:C
            ind = (1+M*(kcom-1)):M*kcom
            S_μσ[i, 2*kcom-1] = sumexp(sumlogmat[i, ind], H1(xtmp[ind], mu[kcom], sigmas[kcom])) / ml[i]
            S_μσ[i, 2*kcom] = sumexp(sumlogmat[i, ind], H2(xtmp[ind], mu[kcom], sigmas[kcom]))/ml[i]
        end
        for j in 1:J
            S_β[i, j] = sumexp(sumlogmat[i,:], summat_beta[i, :, j])/ml[i]
        end
    end
    S_η = hcat(S_β, S_π, S_μσ)
    debuginfo && println(round.(S_η[1:5,:], 5))
    debuginfo && println(round.(sum(S_η, 1)./sqrt(nF), 6))
    I_η = S_η'*S_η./nF
    if 1/cond(I_η) < eps(Float64)
        @warn("Information Matrix is singular!")
        D, V = eigen(I_η)
        debuginfo && println(D)
        tol2 = maximum(abs.(D)) * 1e-14
        D[D.<tol2] = tol2
        I_η = V*diagm(D)*V'
    end
    inv(I_η)[1:J,1:J]./nF
end
function coefpvalue(X::Matrix{Float64}, Y::AbstractArray{Bool, 1}, groupindex::IntegerVector, wi::Vector{Float64}, mu::Vector{Float64}, sigmas::Vector{Float64}, betas::Array{Float64,1}; ngh::Int=100, debuginfo::Bool=false)
    m = vcov(X, Y, groupindex, wi, mu, sigmas, betas, ngh=ngh, debuginfo=debuginfo)
    shat = sqrt.(diag(m))
    p=cdf.(Normal(), betas./shat[1:length(betas)])
    2*min.(p, 1 .- p)
end

####End of utility functions


text:
# top-level functions to run neutral.jl
# Example run:  julia run.jl examples/example1
# Example run:  julia run.jl examples/example1 <seed>    # run with specified random number seed
# Example run:  julia -p 4 run.jl examples/example1 <seed>  # run with 4 processes (cores)
try
  using Distributed
  using Random
  using Dates
catch
end
@everywhere include("NeutralEvolution.jl")

function run_trials( simname::AbstractString )
  global num_trials
  global use_population
  global save_populations
  global log_error    # use log normal error:  work in log space
  global wright_fisher_copy,  wright_fisher_copy_list
  global conformist_probability, conformist_probability_list
  global fit_slope
  global N, mutstddev, initial_value
  if !(simtype == 3 || simtype == 4)
    println("Warning! running neutral.jl on simtype == 2 file.")
  end
  # set defaults for parameters that may not be included in the config file.
  if !isdefined(:use_population)
    use_population = true
  end
  if !isdefined(:save_populations)
    save_populations = false
  end
  if !isdefined(:log_error)
    #println("log_error not defined")
    log_error = false    # default is to not use log error
  end
  if isdefined(:wright_fisher_copy_list)
    wright_fisher_copy = wright_fisher_copy_list[1]
  end
  if !isdefined(:wright_fisher_copy)
    #println("wright_fisher_copy not defined")
    wright_fisher_copy = true     # default is to use wright-fisher copy
  end
  if !isdefined(:wright_fisher_copy_list)
    wright_fisher_copy_list = [wright_fisher_copy]
  end
  if isdefined(:conformist_probability_list)
    conformist_probability = conformist_probability_list[1]
  end
  if !isdefined(:conformist_probability)
    #println("conformist_probability not defined")
    conformist_probability = 0.0     # default is no conformist copy
  end
  if !isdefined(:conformist_probability_list)
    conformist_probability_list = [conformist_probability]
  end
  if !isdefined(:fit_slope)
    fit_slope = 0.0
  end
  if !isdefined(:initial_value)
    initial_value = ideal
  end
  if !isdefined(:N) && isdefined(:N_list)
    N = N_list[1]
    println("Warning! Using N = N_list[1] ")
  end
  if !isdefined(:mutstddev) && isdefined(:mutation_stddev_list)
    mutstddev = mutation_stddev_list[1]
    println("Warning! Using mutstddev = mutation_stddev_list[1]")
  end
  sn = simple_neutral_init( simtype, N, mutstddev, ngens, initial_value, num_trials, record_interval, use_population,
      save_populations, log_error, wright_fisher_copy, conformist_probability, neutral, fit_slope )
  #print_simple_neutral_params( sn )
  csn = cummulative_neutral_init( sn )
  if !csn.save_populations
    writeheader( STDOUT, csn )
    open("$(simname).csv","w") do stream
      writeheader( stream, csn )
    end
  else
    writeheader_populations( STDOUT, csn )
    open("$(simname)_pops.csv","w") do stream
      writeheader_populations( stream, csn )
    end
  end
  for cf_prob in conformist_probability_list
    for wf_copy in wright_fisher_copy_list
      csn = cummulative_neutral_init( sn )
      sn.wright_fisher_copy = wf_copy
      sn.conformist_probability = cf_prob
      if nworkers() > 1
        sn_list_run = simple_neutral_type[]
        run_sim = (sn.N > 1 || sn.use_population) ? simple_neutral_simulation : ces
        for t = 1:num_trials
          Base.push!(sn_list_run,deepcopy(sn))
        end
        sn_list_result = pmap( run_sim, sn_list_run )
        #sn_list_result = map( run_sim, sn_list_run )  # debuging is simpler using map() instead of pmap()
        for rsn in sn_list_result
          accumulate_results( rsn, csn )
        end
      else
        for t = 1:num_trials
          if sn.N > 1 || use_population
            sn = simple_neutral_simulation( sn )
          elseif sn.N == 1
            #sn = ces_log( sn )
            sn = ces( sn )
          else
            error(" sn.N must be positive ")
          end
          accumulate_results( sn, csn )
        end
      end
      if !csn.save_populations
        writerows( STDOUT, csn )
        open("$(simname).csv","a") do stream
          writerows( stream, csn )
        end
      else
        writerows_populations( STDOUT, csn )
        open("$(simname)_pops.csv","a") do stream
          writerows_populations( stream, csn )
        end
      end
    end
  end
end

function coef_var( lst )
  return std(lst)/mean(lst)
end

if length(ARGS) == 0
  simname = "examples/sn_example1"
else
  simname = ARGS[1]
  if length(ARGS) >= 2   # second command-line argument is random number seed
    seed = parse(Int,ARGS[2])
    println("seed: ",seed)
    srand(seed)
  end
end
println("simname: ",simname)
include("$(simname).jl")
println("simtype: ",simtype)
run_trials( simname )


text:
using GLMNet

const GLMNetDefaultParams = Dict(
    :alpha       => 1.0,  # 1.0 means lasso regression, 0.0 means Ridge regression
    :tol         => 1e-7, # convergence criterion
    :standardize => true, # standardize predictors? (beta values always returned on orig scale regardless of this.)
)

abstract GLMNetRegression <: RegressionLearner

# Result is saved
type GLMNetRegressionResult <: RegressionResult
    cv
    best_model::Int
    GLMNetRegressionResult(cv) = new(cv, indmin(cv.meanloss))
end

function coefficients(res::GLMNetRegressionResult)
    res.cv.path.betas[:, res.best_model]
end

function learn(r::GLMNetRegression, dataset::DataSet)
    X = indepvars_as_row_major_matrix(dataset)
    y = depvar_as_vector(dataset)

    # Create params for this specific call based on default ones plus the ones for this
    # regression datum.
    default_params = copy(GLMNetDefaultParams)
    set_parameters!(r, default_params, dataset)
    params = merge(default_params, r.params)

    # Now call GLMNet
    cv = glmnetcv(X, y)
    GLMNetRegressionResult(cv)
end

function coefficients(r::GLMNetRegression, dataset::DataSet)
    res = learn(r, dataset)
    coefficients(res)
end

type LassoRegression <: GLMNetRegression
    params::Dict{Symbol, Any}
    LassoRegression(; params...) = new(Dict(params))
end

function set_parameters!(s::LassoRegression, params, ds::DataSet)
    params[:alpha] = 1.0
    params
end

type RidgeRegression <: GLMNetRegression
    params::Dict{Symbol, Any}
    RidgeRegression(; params...) = new(Dict(params))
end

function set_parameters!(s::RidgeRegression, params, ds::DataSet)
    params[:alpha] = 0.0
    params
end