text:
#=
Copyright (c) 2015, Intel Corporation

Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions are met:

    * Redistributions of source code must retain the above copyright notice,
      this list of conditions and the following disclaimer.
    * Redistributions in binary form must reproduce the above copyright
      notice, this list of conditions and the following disclaimer in the
      documentation and/or other materials provided with the distribution.
    * Neither the name of Intel Corporation nor the names of its contributors
      may be used to endorse or promote products derived from this software
      without specific prior written permission.

THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE
FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
=#

export Neuron, neuron, @neuron
abstract Neuron

function init(neuron::Neuron)
end

"""
Specify a neuron function using the syntax
```
@neuron forward(neuron::NeuronSubType) do
    ...
end
```
"""
macro neuron(expr)
    @assert isa(expr, Expr) "@neuron macro expects an expression"
    if expr.head == :call
        @eval function $(expr.args[1])($(expr.args[3:end]...))
            return $(Expr(:quote, expr.args[2].args[2]))
        end
    elseif expr.head == :type
        type_defn = expr.args[2]
        expr.args[2] = Expr(:(<:), type_defn, :Neuron)
        original = remove_line_nodes(copy(expr.args[3])).args
        T = Float32
        block = expr.args[3].args
        unshift!(block, :(∇inputs :: DenseArray{$T}))
        unshift!(block, :(inputs  :: DenseArray{$T}))
        unshift!(block, :(∇       :: $T))
        unshift!(block, :(value   :: $T))
        args = [a.args[1] for a in original]
        defn = quote
            function $type_defn($(original...))
                $type_defn(zero($T), zero($T), Array($T, 1), Array($T, 1),
                           $(args...))
            end
        end

        result = quote
            $expr
            $defn
        end
        @eval $result
        return Expr(:escape, quote
            typealias $type_defn Latte.$type_defn
        end)
    else
        throw("Not Implemented")
    end
end


text:
## -----------  linear

"""
	$(SIGNATURES)

Admission probability = linear around point (score = 0.5, prob = prob50).

Parameters may vary by college (as indicated by switches). But some colleges are open admission. They have no parameters.

All students get into the lowest college. Parameters only need to be set for colleges 2+.

Change: A lot of code is repeated from the logistic case. +++
    The more general approach is a function with named params (e.g. prob50V, slopeV).
"""
mutable struct AdmProbFctLinearSwitches{F1} <: AbstractAdmProbFctSwitches{F1}
    pvec :: ParamVector
    nColleges :: Int
    nOpenColleges :: Int
end

mutable struct AdmProbFctLinear{F1} <: AbstractAdmProbFct{F1}
    objId :: ObjectId
    switches :: AdmProbFctLinearSwitches{F1}
    prob50V :: Vector{F1}
    slopeV :: Vector{F1}
end

# The order matters here!
param_names(af :: AdmProbFctLinearSwitches{F1}) where F1 =
    param_names(typeof(af));
param_names(::Type{AdmProbFctLinearSwitches{F1}}) where F1 = [:prob50V, :slopeV];

linear_adm_prob(x, p50, slope) = max.(0.0001, min.(0.999, p50 .+ (x .- 0.5) .* slope));

## -------  Constructors

# nc: total no of colleges
function make_test_admprob_fct_linear_switches(nc)
    init_admprob_fct_linear_switches(ObjectId(:test), nc);
end

function make_test_admprob_fct_linear(nc)
    switches = make_test_admprob_fct_linear_switches(nc);
    af = init_admprob_fct(switches);
    @assert validate_admprob_fct(af);
    return af
end

"""
	$(SIGNATURES)

Set up switches for linear admission probability functions.
"""
function init_admprob_fct_linear_switches(
    objId :: ObjectId,  nc :: Integer;
    byCollegeV :: Vector{Symbol} = [:prob50V, :slopeV],
    calibratedV :: Vector{Symbol} = [:prob50V, :slopeV]
    )

    nOpen = 1;
    make_admprob_fct_switches(objId, AdmProbFctLinearSwitches{Float64},
        nc, nOpen, byCollegeV, calibratedV);
end

# notation from symbol table +++++
# Min admission prob for each college.
# Input is no of colleges that need parameters.
function init_prob50(nc :: Integer, isCalibrated :: Bool)
    sz = (nc, );
    return Param(:prob50V, "Adm prob at score = 0.5", "prob50V",
        fill(0.5, sz), fill(0.5, sz),
        fill(0.05, sz), fill(0.85, sz), isCalibrated);
end

function init_slope(nc :: Integer, isCalibrated :: Bool)
    sz = (nc, );
    return Param(:slopeV, "Adm prob function slope", "slopeV",
        fill(0.35, sz), fill(0.35, sz),
        fill(0.0, sz), fill(4.0, sz), isCalibrated);
end

function init_param(::Type{AdmProbFctLinearSwitches{F1}}, pName, nc, isCalibrated) where F1
    if pName == :prob50V
        return init_prob50(nc, isCalibrated);
    elseif pName == :slopeV
        return init_slope(nc, isCalibrated);
    else
        error("Invalid: $pName");
    end
end

# stub ++++++
function validate_admprob_fct_switches(switches :: AdmProbFctLinearSwitches{F1}) where F1
    isValid = true;
    return isValid;
end

"""
	$(SIGNATURES)

Make the object that holds the admission probability functions for a set of colleges.
"""
function init_admprob_fct(switches :: AdmProbFctLinearSwitches{F1}) where F1
    pV = [param_default_value(switches, pName)  for pName in param_names(switches)];
    af = AdmProbFctLinear(get_object_id(switches), switches, pV...);
        # param_default_value(switches, :prob50V),
        # param_default_value(switches, :slopeV));
    @assert validate_admprob_fct(af);
    return af
end

function validate_admprob_fct(switches :: AdmProbFctLinear{F1}) where F1
    isValid = true;
    return isValid;
end

## ---------  Admission prob for one college

"""
	$(SIGNATURES)

Make admission probability function for one college. Maps endowment percentile into [0, 1].
"""
function make_admprob_function(af :: AdmProbFctLinear{F1},
    iCollege :: Integer) where F1

    if iCollege <= n_open_colleges(af)
        fct = x -> 1.0;
    else
        paramV = get_params(af, iCollege);
        fct = x -> linear_adm_prob(x, paramV...);
    end
    return fct
end

"""
	$(SIGNATURES)

Probability of being admitted into a specific college (just based on its admissions prob function).
"""
function prob_admit(af :: AdmProbFctLinear{F1}, iCollege :: Integer, hsGpa) where F1
    if iCollege <= n_open_colleges(af)
        probV = prob_admit_open(hsGpa);
    else
        probV = linear_adm_prob(hsGpa, get_params(af, iCollege)...);
    end
    return probV
end

# --------------------

ext: jl


text:
using Test
#using Documenter
using Distributions, LogNormals

include("fitstats.jl")

include("normal.jl")
include("lognormal.jl")
include("logitnormal.jl")

include("semcor.jl")


text:
function build_stoichiometric_matrix(reactions::Array{String,1};
    expand::Bool=false)::Tuple{Array{Float64,2}, Array{String,1}, Array{String,1}}

	# initialize -
	species_array = Array{String,1}()
	reaction_array = Array{String,1}()
	reaction_dictionary_array = Array{Dict{String,Float64},1}()

    # should we expand the reversible reactions?
    reactions_to_process = reactions;
    if (expand == true)
        reactions_to_process = expand_reversible_reactions(reactions);
    end

	# first: let's discover the species list -
	for reaction_string ∈ reactions_to_process

		# initialize tmp storage -
		tmp_dictionary = Dict{String,Float64}()

		# split the reaction into its components -
		component_array = split(reaction_string,',');

		# reaction name -
		reaction_name = String.(component_array[1]);
		push!(reaction_array, reaction_name);

		# reactant phrase => 2, and product phrase => 3
		reactant_phrase = String.(component_array[2]);
		product_phrase = String.(component_array[3]);

		# generate species lists for the reactants and products, then merge -
		merge!(tmp_dictionary, extract_species_dictionary(reactant_phrase; direction = -1.0))
		merge!(tmp_dictionary, extract_species_dictionary(product_phrase; direction = 1.0))

		# grab the tmp_dictionary for later -
		push!(reaction_dictionary_array, tmp_dictionary)

		# the species that we need to look at are the keys of the tmp_dictionary -
		tmp_species_list = keys(tmp_dictionary)

		# we need a unique species list, so check to see if we have already discovered this species -
		for tmp_species ∈ tmp_species_list

			if (in(tmp_species, species_array) == false)

				# ok, we have *not* seen this species before, let's grab it -
				push!(species_array, tmp_species)
			end
		end
	end

	# sort alphabetically -
	sort!(species_array)

	# we have a *unique* species array, let's initialize some storage for the stoichiometric array
	S = zeros(length(species_array), length(reactions_to_process));

	# last: fill in the values for stoichiometric coefficents -
	for (row_index, species_symbol) ∈ enumerate(species_array)
		for (col_index, reaction_dictionary) ∈ enumerate(reaction_dictionary_array)

			# ok: is this species symbol in my reaction dictionary?
			if (haskey(reaction_dictionary, species_symbol) == true)
				S[row_index,col_index] = reaction_dictionary[species_symbol]
			end
		end
	end

	# return -
	return (S, species_array, reaction_array)
end

function extract_species_dictionary(reaction_phrase::String;
	direction::Float64 = -1.0)::Dict{String,Float64}

	# initialize -
	species_symbol_dictionary = Dict{String,Float64}()

	# ok, do we hve a +?
	component_array = split(reaction_phrase,'+');
	for component ∈ component_array

		if (contains(component,'*') == true)

			tmp_array = split(component,'*')
			st_coeff = direction*parse(Float64,tmp_array[1])
			species_symbol = String(tmp_array[2])

			# don't cache the ∅ -
			if (species_symbol != "∅" && species_symbol != "[]")
				species_symbol_dictionary[species_symbol] = st_coeff
			end
		else

			# strip any spaces -
			species_symbol = component |> lstrip |> rstrip

			# don't cache the ∅ -
			if (species_symbol != "∅" && species_symbol != "[]")
				species_symbol_dictionary[species_symbol] = direction*1.0
			end
		end
	end

	# return -
	return species_symbol_dictionary
end

function binary_stoichiometric_matrix(matrix::Array{Float64,2})::Array{Int64,2}

	# initialize -
	(ℳ,ℛ) = size(matrix)
	B = Array{Int64,2}(undef,ℳ,ℛ)

	for row_index ∈ 1:ℳ
		for col_index ∈ 1:ℛ

			old_value = matrix[row_index,col_index]
			if (old_value == 0.0)
				B[row_index,col_index] = 0
			else
				B[row_index,col_index] = 1
			end
		end
	end

	# return -
	return B
end

function expand_reversible_reactions(reaction_array::Array{String,1})::Array{String,1}

    # initialize -
    processed_reaction_string_array = Array{String,1}()

    # main loop -
    for reaction_string ∈ reaction_array

        # chop up the reaction string -
	    reaction_component_array = split(reaction_string,',');

        # get components -
        rname = reaction_component_array[1]
        forward_phrase = reaction_component_array[2]
        reverse_phrase = reaction_component_array[3]
        is_reversible = parse(Bool, reaction_component_array[4])

        if (is_reversible == true)

            # build new reaction strings (forward, and reverse)
            new_string_forward = "F$(rname),$(forward_phrase),$(reverse_phrase),false"
            new_string_reverse = "R$(rname),$(reverse_phrase),$(forward_phrase),false"
            push!(processed_reaction_string_array,new_string_forward)
            push!(processed_reaction_string_array,new_string_reverse)
        else
            push!(processed_reaction_string_array, reaction_string)
        end
    end

    # return -
    return processed_reaction_string_array
end


text:
@testset "Base interop" begin
    @testset "open() and mktemp()" begin
        path = ""
        io1 = nothing
        io2 = nothing
        io3 = nothing
        io4 = nothing
        @context begin
            (path,io1) = @! mktemp()
            msg = "hi from ResourceContexts.jl"
            write(io1, msg)
            flush(io1)
            # Test various forms of open()
            io2 = @! open(path)
            @test read(io2, String) == msg
            io3 = @! open(path, "r")
            @test read(io3, String) == msg
            io4 = @! open(path, read=true)
            @test read(io4, String) == msg
        end
        @test !isempty(path) && !ispath(path)
        @test !isopen(io1)
        @test !isopen(io2)
        @test !isopen(io3)
        @test !isopen(io4)
    end

    @testset "mktemp() and mktempdir()" begin
        parent_dir = nothing
        dir2 = nothing
        filepath = nothing
        @context begin
            parent_dir = @! mktempdir()
            # forms with parent dirs
            dir2 = @! mktempdir(parent_dir)
            (filepath,_) = @! mktemp(parent_dir)
            # All resources are created
            @test isdir(parent_dir)
            @test isdir(dir2)
            @test isfile(filepath)
            @test parent_dir == dirname(dir2)
            @test parent_dir == dirname(filepath)
        end
        # All resources are cleaned up
        @test !ispath(parent_dir)
        @test !ispath(dir2)
        @test !ispath(filepath)
    end

    @testset "cd()" begin
        oldpath = pwd()
        @context begin
            @! cd()
            # Can't just compare paths here due to the fact that some of these
            # could be symlinks on macOS.  Instead use `samefile` with an alias
            # for @test pretty printing
            ≃(a,b) = Base.Filesystem.samefile(a,b)
            @test pwd() ≃ homedir()
            path = @! mktempdir()
            @! cd(path)
            @test pwd() ≃ path
        end
        @test pwd() == oldpath
    end

    @testset "lock()" begin
        lk = ReentrantLock()
        @context begin
            @! lock(lk)
            @test islocked(lk)
        end
        @test !islocked(lk)
    end

    @testset "redirect_stdout / stderr / stdin" begin
        orig_stdin = stdin
        orig_stdout = stdout
        orig_stderr = stderr

        stdout_result = nothing
        @context begin
            (_,io) = @! mktemp()
            @! redirect_stdout(io)
            println("hi")
            flush(io)
            seek(io, 0)
            stdout_result = readline(io)
        end
        @test stdout_result == "hi"

        stderr_result = nothing
        @context begin
            (_,io) = @! mktemp()
            @! redirect_stderr(io)
            println(stderr, "hi")
            flush(io)
            seek(io, 0)
            stderr_result = readline(io)
        end
        @test stderr_result == "hi"

        stdin_result = nothing
        @context begin
            (_,io) = @! mktemp()
            println(io, "hi")
            flush(io)
            seek(io,0)
            @! redirect_stdin(io)
            stdin_result = readline()
        end
        @test stdin_result == "hi"

        @test orig_stdin == stdin
        @test orig_stdout == stdout
        @test orig_stderr == stderr
    end
end


text:
module FamaMacBeth

import StatsBase, StatsBase.coef, StatsBase.stderror
using Statistics: mean, std, cov

export fmb, coef, stderror, tstat, nobs, predict, residuals, r2, adjr2, group

# response object
struct FMBResp
    group
    coef::Array
    nobs::Vector{Int}
end

"""
    FMBModel(X,y,g,resp::FMBResp) <: StatsBase.RegressionModel
"""
mutable struct FMBModel <: StatsBase.RegressionModel
    X
    y
    g
    resp::Union{FMBResp, Nothing}

    function FMBModel(X,y,g,resp)
        # checks on input
        @assert size(X,1)==length(y) "y and X must have the same number of observations"
        @assert length(y)==length(g) "g and y must have the same number of observations"
        new(X,y,g,resp)
    end
end

"""
    FMBModel(X,y,g) -> FMBModel()
"""
FMBModel(X,y,g) = FMBModel(X,y,g,nothing)

"""
    group(m::FMBResp)
    group(m::FMBModel)

Returns the group labels used.
"""
group(m::FMBResp) = m.group
group(m::FMBModel) = group(m.resp)

"""
    coef(m::FMBResp) -> Vector
    coef(m::FMBModel) -> Vector

The estimated coefficients.
"""
coef(m::FMBResp) = mean(m.coef, dims=2)
coef(m::FMBModel) = coef(m.resp)

"""
    stderror(m::FMBResp) -> Vector
    stderror(m::FMBModel) -> Vector

The standard error around the estimated coefficients.
"""
stderror(m::FMBResp) = std(m.coef, dims=2)./sqrt(size(m.coef,2))
stderror(m::FMBModel) = stderror(m.resp)

"""
    vcov(m::FMBResp) -> Matrix
    vcov(m::FMBModel) -> Matrix

Variance-covariance matrix of the estimated coefficients.
"""
vcov(m::FMBResp) = cov(m.coef, dims=2)./size(m.coef,2)
vcov(m::FMBModel) = vcov(m.resp)

"""
    tstat(m::FMBModel) -> Vector
    tstat(m::FMBModel) -> Vector

The t-statistic per coefficient whether the coefficient equals zero or not.
"""
tstat(m::FMBResp) = coef(m)./stderror(m)
tstat(m::FMBModel) = tstat(m.resp)

"""
    nobs(m::FMBModel) -> Vector{Int}

Returns the number of observations, per group as identified by `g`. See also `group(m::FMBResp)`.
"""
nobs(m::FMBModel) = m.resp.nobs

"""
    predict(m::FMBModel) -> Vector
    predict(m::FMBModel, X) -> Vector

Predicts the outcome using the estimated coefficients and provided independent variables.
If independent variables are not provided, uses those from the model.
"""
predict(m::FMBModel, X) = X * coef(m)
predict(m::FMBModel) = predict(m, m.X)

"""
    residuals(m::FMBModel) -> Vector
"""
residuals(m::FMBModel) = m.y .- predict(m)

"""
    r2(m::FMBModel) -> Float64

The un-adjusted R2.
"""
function r2(m::FMBModel)
    SSE = sum(residuals(m).^2)
    SST = sum((m.y .- mean(m.y)).^2)
    1-SSE/SST
end

"""
    adjr2(m::FMBModel) -> Float64

The adjusted R2.
"""
function adjr2(m::FMBModel)
    n = length(m.y)
    p = size(m.X,2)
    1.0-(1.0-r2(m))*(n-1)/(n-p)
end

"""
    fmb(X, y, g) -> FMBModel()

Short-hand for fit(FMBModel(X,y,g))
"""
fmb(X::Array, y, g) = fit(FMBModel(X,y,g))

"""
    fit(m::FMBModel) -> FMBModel()

Runs a Fama-MacBeth (1973) panel regression by first conducting OLS of `y` on `X`
per group indicated by `g`, and subsequently computes the mean and standard error
using the resulting coefficients. An intercept is not added to `X`.
"""
function fit(m::FMBModel)
    # pre-allocate output
    ug = unique(m.g)
    coef = zeros(eltype(m.X), size(m.X,2), length(ug))
    nobs = zeros(Int64, length(ug))

    # per group: OLS
    for (i, gs) in enumerate(ug)
        select = m.g.==gs
        coef[:,i] = m.X[select,:]\m.y[select]
        nobs[i] = sum(select)
    end

    m.resp = FMBResp(ug, coef, nobs)
    return m
end

end


text:
__precompile__(true)

module GCMAES

using Printf, Distributed, LinearAlgebra
using Dates, Random, Statistics, Pkg
using Requires, BSON, FastClosures

include("util.jl")
include("constraint.jl")
include("transform.jl")

function __init__()
    @require MPI="da04e1cc-30fd-572f-bb4f-1f8673147195" include("mpi.jl")
    @require Elemental="902c3f28-d1ec-5e7e-8399-a24c3845ee38" include("elemental.jl")
end

mutable struct CMAESOpt{T, F, G, CO, TR}
    # fixed hyper-parameters
    f::F
    g::G
    N::Int
    σ0::T
    lo::Vector{T}
    hi::Vector{T}
    constr::CO
    trans::TR
    rng::MersenneTwister
    # strategy parameter setting: selection
    λ::Int
    μ::Int
    w::Vector{T}
    μeff::T
    χₙ::T
    # strategy parameter setting: adaptation
    σ::T
    cc::T
    cσ::T
    c1::T
    cμ::T
    dσ::T
    # dynamic (internal) strategy parameters and constants
    x̄::Vector{T}
    pc::Vector{T}
    pσ::Vector{T}
    D::Vector{T}
    B::Matrix{T}
    BD::Matrix{T}
    C::Matrix{T}
    arx::Matrix{T}
    ary::Matrix{T}
    arz::Matrix{T}
    arfitness::Vector{T}
    arindex::Vector{Int}
    # recordings
    xmin::Vector{T}
    fmin::T
    fmins::Vector{T} # history of best fitness
    fmeds::Vector{T} # history of median fitness
    feqls::Vector{T} # history of equal fitness
    # report
    start_time::Float64
    last_report_time::Float64
    pmap_time::Float64
    grad_time::Float64
    ls_time::Float64
    ls_dec::T
    file::String
end

function CMAESOpt(f, g, x0, σ0, lo = -fill(1, size(x0)), hi = fill(1, size(x0));
            λ = 0, constr = NoConstraint(), trans = false, rng = MersenneTwister())
    trans = trans == true ? BoxLinQuadTransform(lo, hi) :
            trans == false ? NoTransform() : trans
    x0′ = inverse(trans, x0)
    f′ = @closure z -> getfitness(f, constr, transform(trans, z))
    g′ = @closure z -> g(transform(trans, z))
    N, x̄, xmin, fmin, σ = length(x0), x0′, x0′, f′(x0′), σ0
    # strategy parameter setting: selection
    λ = λ == 0 ? round(Int, 4 + 3log(N)) : max(4, λ)
    μ = ceil(Int, λ / 2)                   # number of parents/points for recombination
    w = log(μ + 1/2) .- log.(1:μ)          # μXone array for weighted recombination
    normalize!(w, 1)                       # normalize recombination w array
    μeff = 1 / sum(abs2, w)                # variance-effectiveness of sum w_i x_i
    # strategy parameter setting: adaptation
    cc = 4 / (N + 4)                       # time constant for cumulation for C
    cσ = (μeff + 2) / (N + μeff + 3)       # t-const for cumulation for σ control
    c1 = 2 / ((N + 1.3)^2 + μeff)          # learning rate for rank-one update of C
    cμ = min(1 - c1, 2(μeff - 2 + 1 / μeff) / ((N + 2)^2 + μeff) )   # and for rank-μ update
    dσ = 1 + 2 * max(0, sqrt((μeff - 1) / (N + 1)) - 1) + cσ # damping for σ, usually close to 1
    # initialize dynamic (internal) strategy parameters and constants
    pc = zeros(N); pσ = zeros(N)            # evolution paths for C and σ
    D = fill(σ, N)                          # diagonal matrix D defines the scaling
    B = Matrix(I, N, N)                     # B defines the coordinate system
    BD = B .* reshape(D, 1, N)              # B*D for speed up only
    C = diagm(0 => D.^2)                    # covariance matrix == BD*(BD)'
    χₙ = sqrt(N) * (1 -1 / 4N + 1 / 21N^2)   # expectation of  ||N(0,I)|| == norm(randn(N,1))
    # init a few things
    arx, ary, arz = zeros(N, λ), zeros(N, λ), zeros(N, λ)
    arfitness, arindex = zeros(λ), ones(λ)
    @master @printf("%i-%i CMA-ES\n", λ, μ)
    # gradient
    T, F, G = eltype(x̄), typeof(f′), typeof(g′)
    CO, TR = typeof(constr), typeof(trans)
    return CMAESOpt{T, F, G, CO, TR}(
            f′, g′, N, σ0, lo, hi,
            constr, trans, rng,
            λ, μ, w, μeff, χₙ,
            σ, cc, cσ, c1, cμ, dσ,
            x̄, pc, pσ, D, B, BD, C,
            arx, ary, arz, arfitness, arindex,
            xmin, fmin, T[], T[], T[],
            time(), time(), 0, 0, 0, 0, "CMAES.bson")
end

function update_candidates!(opt::CMAESOpt)
    # generate and evaluate λ offspring
    randn!(opt.rng, opt.arz) # resample
    opt.ary = opt.BD * opt.arz
    opt.arx .= opt.x̄ .+ opt.σ .* opt.ary
    arx_cols = [opt.arx[:, k] for k in 1:opt.λ]
    @assert allequal(sum(opt.arx))
    opt.pmap_time = @elapsed opt.arfitness .= pmap(opt.f, arx_cols)
    # sort by fitness and compute weighted mean into x̄
    sortperm!(opt.arindex, opt.arfitness)
    opt.arfitness = opt.arfitness[opt.arindex] # minimization
    # store the best candidate
    if opt.arfitness[1] < opt.fmin
        opt.xmin, opt.fmin = opt.arx[:, opt.arindex[1]], opt.arfitness[1]
    end
    push!(opt.fmins, opt.arfitness[1])
    push!(opt.fmeds, opt.arfitness[ceil(Int, opt.λ / 2)])
    feql = opt.λ > 1 ? opt.arfitness[1] == opt.arfitness[ceil(Int, 1.1 + opt.λ / 4)] : 0
    push!(opt.feqls, feql)
end

function linesearch(f, x0::Array{T}, Δ) where T
    nrm = norm(Δ)
    nrm == 0 && return x0, typemax(T), zero(T)
    rmul!(Δ, 1 / nrm)
    αs = T[0.0; 2.0.^(2 - worldsize():0)]
    xs = [x0 .+ α .* Δ for α in αs]
    fs = pmap(f, xs)
    fx, i = findmin(fs)
    return xs[i], fx, fs[1] - fx
end

function update_mean!(opt::CMAESOpt)
    get(ENV, "GCMAES", "true") == "false" && return
    opt.grad_time = @elapsed Δ = -opt.g(opt.x̄)
    opt.ls_time = @elapsed opt.x̄, fx, opt.ls_dec = linesearch(opt.f, opt.x̄, Δ)
    if fx < opt.fmin copyto!(opt.xmin, opt.x̄); opt.fmin = fx end
end

function update_parameters!(opt::CMAESOpt, iter, lazydecomp)
    indμ = opt.arindex[1:opt.μ]
    # calculate new x̄, this is selection and recombination
    x̄old = copy(opt.x̄)                                # for speed up of Eq. (2) and (3)
    mul!(opt.x̄, opt.arx[:, indμ], opt.w)
    # use parallel line search to determin the best α s.t. x += α * Δ achieves minimum
    update_mean!(opt)
    # calculate new z̄
    z̄ = opt.arz[:, indμ] * opt.w                       # == D^-1 * B' * (x̄ - x̄old) / σ
    # cumulation: update evolution paths
    BLAS.gemv!('N', sqrt(opt.cσ * (2 - opt.cσ) * opt.μeff), opt.B, z̄, 1 - opt.cσ, opt.pσ)  # i.e. pσ = (1 - cσ) * pσ + sqrt(cσ * (2 - cσ) * μeff) * (B * z̄) Eq.(4)
    hsig = norm(opt.pσ) / sqrt(1 - (1 - opt.cσ)^2iter) / opt.χₙ < 1.4 + 2 / (opt.N + 1)
    rmul!(opt.pc, 1 - opt.cc)
    BLAS.axpy!(hsig * sqrt(opt.cc * (2 - opt.cc) * opt.μeff) / opt.σ, opt.x̄ - x̄old, opt.pc) # i.e. pc = (1 - cc) * pc + (hsig * sqrt(cc * (2 - cc) * μeff) / σ) * (x̄ - x̄old)
    # adapt covariance matrix C
    rmul!(opt.C, (1 - opt.c1 - opt.cμ + (1 - hsig) * opt.c1 * opt.cc * (2 - opt.cc)))      # discard old C
    BLAS.syr!('U', opt.c1, opt.pc, opt.C)               # rank 1 update C += c1 * pc * pc'
    artmp = opt.ary[:, indμ]                            # μ difference vectors
    artmp = (artmp .* reshape(opt.w, 1, opt.μ)) * transpose(artmp)
    BLAS.axpy!(opt.cμ, artmp, opt.C)
    # adapt step size σ
    opt.σ *= exp((norm(opt.pσ) / opt.χₙ - 1) * opt.cσ / opt.dσ)  #Eq. (5)
    # update B and D from C
    if !lazydecomp || mod(iter, 1 / (opt.c1 + opt.cμ) / opt.N / 10) < 1 # if counteval - eigeneval > λ / (c1 + cμ) / N / 10  # to achieve O(N^2)
        (opt.D, opt.B) = deigen(opt.C)                  # eigen decomposition, B == normalized eigenvectors
        opt.D .= sqrt.(opt.D)                           # D contains standard deviations now
        opt.BD .= opt.B .* reshape(opt.D, 1, opt.N)     # O(n^2)
    end
end

function terminate(opt::CMAESOpt, equal_best = Inf)
    histiter = 10 + round(Int, 30opt.N / opt.λ)
    stagiter = round(Int, 0.2 * length(opt.fmins) + 120 + 30opt.N / opt.λ)
    stagiter = min(stagiter, 20000)
    hist = opt.fmins[max(1, end - histiter):end]
    condition = Dict{String, Bool}()
    # FlatFitness: warn if 70% candidates' fitnesses are identical
    if opt.arfitness[1] == opt.arfitness[ceil(Int, 0.7opt.λ)]
        opt.σ *= exp(0.2 + opt.cσ / opt.dσ)
        @master println("warning: flat fitness, consider reformulating the objective")
    end
    # Stop conditions:
    # MaxIter: the maximal number of iterations in each run of CMA-ES
    length(opt.fmins) >= 100 + 50(opt.N + 3)^2 / sqrt(opt.λ) && get!(condition, "MaxIter", true)
    # TolHistFun: the range of the best function values during the last
    #             10 + ⌈30N/λ⌉ iterations is smaller than TolHistFun = 1e-13
    length(hist) > histiter && ptp(hist) < 1e-13 && get!(condition, "TolHistFun", true)
    # EqualFunVals: in more than 1/3rd of the last D iterations the objective function value of the best and
    #               the k-th best solution are identical, that is f(x1:λ) = f(xk:λ), where k = 1 + ⌈0.1 + λ/4⌉
    length(opt.fmins) > opt.N && mean(opt.feqls[end - opt.N:end]) > 0.33 && get!(condition, "EqualFunVals", true)
    # TolX: all components of pc and all square roots of diagonal components of C,
    #       multiplied by σ/σ0, are smaller than TolX = 1e-12.
    sqrt(maximum(abs, diag(opt.C))) * opt.σ / opt.σ0 < 1e-11 &&
    maximum(abs, opt.pc) * opt.σ / opt.σ0 < 1e-11 &&
    get!(condition, "TolX", true)
    # TolUpSigma: σ/σ0 >= maximum(D) * TolUpSigma = 1e20
    opt.σ / opt.σ0 > 1e20 * maximum(opt.D) && get!(condition, "TolUpSigma", true)
    # Stagnation: we track a history of the best and the median fitness in each iteration
    #             over the last 20% but at least 120 + 30n/λ and no more than 20000 iterations.
    #             We stop, if in both histories the median of the last (most recent) 30% values
    #             is not better than the median of the first 30%
    length(opt.fmins) > stagiter &&
    median(opt.fmins[end - 20:end]) >= median(opt.fmins[end - stagiter:end - stagiter + 20]) &&
    median(opt.fmeds[end - 20:end]) >= median(opt.fmeds[end - stagiter:end - stagiter + 20]) &&
    get!(condition, "Stagnation", true)
    # ConditionCov: the condition number of C exceeds 1e14
    maximum(opt.D) > 1e7 * minimum(opt.D) && get!(condition, "ConditionCov", true)
    # TolFun: stop if the range of the best objective function values of the last 10+⌈30n/λ⌉
    #         generations and all function values of the recent generation is below TolFun = 1e-12.
    length(hist) > histiter && ptp(vcat(hist, opt.arfitness)) < 1e-11 && get!(condition, "TolFun", true)
    # EqualBest: terminate is the range of the fmins of the last EqualBest iterations is smaller than 1e-12
    length(opt.fmins) > equal_best &&
    ptp(accumulate(min, opt.fmins)[end - equal_best:end]) < 1e-12 &&
    get!(condition, "EqualBest", true)
    # NoEffectAxis: stop if adding a 0.1-standard deviation vector in any principal axis
    #               direction of C does not change m.
    # NoEffectCoor: stop if adding 0.2-standard deviations in any single coordinate does not
    #               change m (i.e. mᵢ equals mᵢ + 0.2σcᵢᵢ, for any i)
    # TolFacUpX: terminate when step-size increases by TolFacUpx = 1e3 (diverges). That is,
    #            the initial step-size was chosen far too small and better solutions were found
    #            far away from the initial solution x0.
    # https://github.com/DEAP/deap/blob/8b8c82186083af4b78bf67d9bb3508cf1489a4be/examples/es/cma_bipop.py
    # http://cma.gforge.inria.fr/cmaes.m
    # Benchmarking a BI-Population CMA-ES on the BBOB-2009 Function Testbed
    termination = false
    for (k, v) in condition
        @master v && printstyled("Termination Condition Satisfied: ", k, '\n', color = :red)
        termination = termination | v
    end
    return termination
end

function restart(opt::CMAESOpt)
    @master @printf("restarting...\n")
    opt′ = CMAESOpt(opt.f, sample(opt.lo, opt.hi), opt.σ0, opt.lo, opt.hi;
                    λ = 2opt.λ, constr = opt.constr, trans = opt.constr)
    opt′.xmin, opt′.fmin = opt.xmin, opt.fmin
    return opt′
end

function trace_state(opt::CMAESOpt, iter, fcount)
    elapsed_time = time() - opt.last_report_time
    # display some information every iteration
    @master @printf("time:%s  iter:%d  elapsed-time:%.2f  ", Time(now()), iter, elapsed_time)
    @master @printf("pmap-time:%.2f  grad-time:%.2f  ls-time:%.2f ls-dec:%2.2e\n",
            opt.pmap_time, opt.grad_time, opt.ls_time, opt.ls_dec)
    @master @printf("fcount:%d  fval:%2.2e  fmin:%2.2e  ", fcount, opt.arfitness[1], opt.fmin)
    @master @printf("norm:%2.2e  axis-ratio:%2.2e  free-mem:%.2fGB\n",
            norm(opt.arx[:, opt.arindex[1]]), maximum(opt.D) / minimum(opt.D), Sys.free_memory() / 1024^3)
    opt.last_report_time = time()
    return nothing
end

function load!(opt::CMAESOpt, resume)
    (resume == false || !isfile(opt.file)) && return
    data = BSON.load(opt.file)
    data[:N] != opt.N && return
    loadvars = [:σ, :cc, :cσ, :c1, :cμ, :dσ, :x̄, :pc, :pσ, :D, :B, :BD, :C]
    resume == :full && append!(loadvars, [:xmin, :fmin, :fmins, :fmeds, :feqls])
    for s in loadvars
        haskey(data, s) && setfield!(opt, s, data[s])
    end
end

function save(opt::CMAESOpt, saveall = false)
    data = Dict{Symbol, Any}(:N => opt.N)
    savevars = [:σ, :cc, :cσ, :c1, :cμ, :dσ, :x̄, :pc, :pσ, :D, :B, :BD, :C, :xmin, :fmin, :fmins, :fmeds, :feqls]
    if !saveall && Base.summarysize(opt) >= 1024^2 * 20
        savevars = setdiff(savevars, [:D, :B, :BD, :C])
    end
    for s in savevars
        data[s] = copy(getfield(opt, s))
    end
    @master BSON.bson(opt.file, data)
end

function minimize(fg, x0, a...; maxfevals = 0, gcitr = false, maxiter = 0,
                resume = false, cb = [], seed = nothing, autodiff = false,
                saveall = false, lazydecomp = false, equal_best = Inf, print_interval, maxtime = 300,ka...)
    rng = bcast(MersenneTwister(seed))
    f, g = fg isa Tuple ? fg : autodiff ? (fg, fg') : (fg, zero)
    opt = CMAESOpt(f, g, bcast(x0), a...; rng = rng, ka...)
    cb = runall([throttle(x -> save(opt, saveall), 300); cb])
    maxfevals = (maxfevals == 0) ? 1e3 * length(x0)^2 : maxfevals
    maxfevals = maxiter != 0 ? maxiter * opt.λ : maxfevals
    load!(opt, resume)
    iter = length(opt.fmins)
    fcount = iter * opt.λ
    status = 0
    while fcount < maxfevals && ((time() - opt.start_time) < maxtime)
        iter += 1; fcount += opt.λ
        if opt.λ == 1
            update_mean!(opt)
            trace_state(opt, iter, fcount)
            continue
        end
        update_candidates!(opt)
        update_parameters!(opt, iter, lazydecomp)
        if (iter % print_interval == 0)
            trace_state(opt, iter, fcount)
        end
        gcitr && @everywhere GC.gc(true)
        cb(opt.xmin) == :stop && break
        terminate(opt, equal_best) && (status = 1; break)
        # if terminate(opt, equal_best) opt, iter = restart(opt), 0 end
    end
    xmin = inverse(opt.trans, opt.xmin)
    return xmin, opt.fmin, status, opt
end

function maximize(fg, args...; kws...)
    f, g = fg isa Tuple ? fg : (fg, zero)
    fg′ = (x -> -f(x), x -> -g(x))
    xmin, fmin, status = minimize(fg′, args...; kws...)
    return xmin, -fmin, status
end

end


text:
using Statistics

function get_timings_for_file(input_file_string, mean_timings,
  median_timings)
  #== read in input file ==#
  input_file::Array{String,1} = []
  open(input_file_string) do file
    input_file = readlines(file)
  end

  for time_to_parse in collect(keys(mean_timings))
    #== extract timings from input file ==#
    timings_string::Array{String,1} = []
    for line in input_file
      if occursin("$time_to_parse TIME",line)
        push!(timings_string, line)
      end
    end
    #display(timings_string)

    timings::Array{Float64,1} = []
    for line in timings_string
      timing_match = eachmatch(r"[0-9]+.[0-9]+", line)
      timing = parse(Float64, collect(timing_match)[1].match)
      push!(timings, timing)
    end
    #display(timings)
    #println()

    avg = mean(timings)
    med = median(timings)
    stdev = std(timings)
    reldev = 100*stdev/avg

    push!(mean_timings["$time_to_parse"], avg)
    push!(median_timings["$time_to_parse"], med)
  end
end

#=
function analyze_timing(input_file_string::String)
  #== read in input file ==#
  input_file::Array{String,1} = []
  open(input_file_string) do file
    input_file = readlines(file)
  end

  for time_to_parse in ["DECOMPOSE", "ASSIGN", "INTEGRALS", "FOCK",
    "KERNEL SUM", "TWOEI"]
    #== extract timings from input file ==#
    timings_string::Array{String,1} = []
    for line in input_file
      if occursin("$time_to_parse TIME",line)
        push!(timings_string, line)
      end
    end
    #display(timings_string)

    timings::Array{Float64,1} = []
    for line in timings_string
      timing_match = eachmatch(r"[0-9]+.[0-9]+", line)
      timing = parse(Float64, collect(timing_match)[1].match)
      push!(timings, timing)
    end
    #display(timings)
    #println()

    avg = mean(timings)
    med = median(timings)
    stdev = std(timings)
    reldev = 100*stdev/avg

    println("Mean $time_to_parse time: $avg")
    println("Median $time_to_parse time: $med")
    println("Standard Deviation $time_to_parse time: $stdev")
    println("Relative Deviation $time_to_parse time: $reldev")
    println()
  end
end
=#

files = [ "01_MP2.log",
          "02_MP2.log",
          "03_MP2.log",
          "04_MP2.log",
          "05_MP2.log",
          "06_MP2.log",
          "07_MP2.log",
          "08_MP2.log",
          "09_MP2.log",
          "10_MP2.log",
          "11_MP2.log",
          "12_MP2.log",
          "13_MP2.log",
          "14_MP2.log",
          "15_MP2.log",
          "16_MP2.log",
          "17_MP2.log",
          "18_MP2.log",
          "19_MP2.log",
          "20_MP2.log",
          "21_MP2.log",
          "22_MP2.log"
        ]

mean_timings = Dict("DECOMPOSE" => Vector{Float64}([]),
                    "ASSIGN" => Vector{Float64}([]),
                    "INTEGRALS" => Vector{Float64}([]),
                    "FOCK" => Vector{Float64}([]),
                    "KERNEL SUM" => Vector{Float64}([]),
                    "TWOEI" => Vector{Float64}([])
                   )

median_timings = Dict("DECOMPOSE" => Vector{Float64}([]),
                      "ASSIGN" => Vector{Float64}([]),
                      "INTEGRALS" => Vector{Float64}([]),
                      "FOCK" => Vector{Float64}([]),
                      "KERNEL SUM" => Vector{Float64}([]),
                      "TWOEI" =>Vector{Float64}([])
                     )

for file in files
  file_location = joinpath(@__DIR__, "example_inputs/S22-profile/$file")
  get_timings_for_file(file_location, mean_timings, median_timings)
end

for timing_key in collect(keys(mean_timings))
  println("MEAN $timing_key TIMINGS:")
  for timing in mean_timings[timing_key]
    println(timing)
  end
    println()
end


text:
using Distributions
import Distributions
using Random: AbstractRNG

export invert, primdist, distapply
export StdNormal, StdUniform
export PrimDist, ExoRandVar

"Primitive Distribution"
abstract type PrimDist end

"An exogenous random variable"
ExoRandVar{F <: PrimDist, ID} = Member{F, ID}

struct StdNormal{T<:Real} <: PrimDist end
Base.eltype(::Type{StdNormal{T}}) where T = T
Distributions.logpdf(::StdNormal{T}, x) where T =
  Distributions.logpdf(Normal(zero(T), one(T)), x)
Base.rand(rng::AbstractRNG, ::StdNormal{T}) where {T} = rand(rng, Normal(zero(T), one(T)))

# This is called from dispatch
@inline (d::Normal{T})(id, ω) where T =
  Member(id, StdNormal{T}())(ω) * d.σ + d.μ

struct StdUniform{T} <: PrimDist end
Base.eltype(::Type{StdUniform{T}}) where T = T
Base.rand(rng::AbstractRNG, ::StdUniform{T}) where T =
  rand(rng, Uniform(zero(T), one(T)))
# @inline Space.recurse(d::Distribution, id, ω) =
#   quantile(d, resolve(StdUniform(), id, ω))

@inline (d::Bernoulli)(id, ω) =
  Member(id, StdUniform{Float64}())(ω) < d.p

@inline (d::Distribution)(id, ω) =
  quantile(d, Member(id, StdUniform{Float64}())(ω))

"""`primdist(d::Distribution)``
Primitive (parameterless) distribution that `d` is defined in terms of"""
function primdist end

primdist(d::Distribution) = StdUniform()
primdist(d::Normal) = StdNormal()

"""`invert(d::Distribution, val)`
If output of `val` is `val` what must its primitives have been?`"""
function invert end

invert(o::Normal, val) = (val / o.σ) - o.μ
invert(d::Distribution, val) = cdf(d, val)


text:
module SankeyPlots

using LayeredLayouts
using LightGraphs
using Plots
using RecipesBase
using SimpleWeightedGraphs
using SparseArrays

export sankey, sankey!

"""
    sankey(src, dst, weights; kwargs..., plotattributes...)
    sankey(g::SimpleWeightedDiGraph; kwargs..., plotattributes...)

Plot a sankey diagram.

In addition to [Plots.jl attributes](http://docs.juliaplots.org/latest/attributes/) the following keyword arguments are supported.

| Keyword argument | Default value | Options |
|---|---|----|
| `node_labels` | `nothing` | `AbstractVector{<:String}` |
| `node_colors` | `nothing` | Vector of [color specifications supported by Plots.jl](http://docs.juliaplots.org/latest/colors/) or [color palette](http://docs.juliaplots.org/latest/generated/colorschemes/#ColorPalette) |
| `edge_color` | `:gray` | Plots.jl supported [color](http://docs.juliaplots.org/latest/colors/) or color selection from connected nodes with `:src`, `:dst` or `:gradient` |
| `label_position` | `:inside` | `:legend`, `:node`, `:left`, `:right`, `:top` or `:bottom` |
| `label_size` | `8` | `Int` |
| `compact` | `false` | `Bool` |
"""
@userplot Sankey

@recipe function sankey(
    s::Sankey;
    node_labels=nothing,
    node_colors=nothing,
    edge_color=:gray,
    label_position=:inside,
    label_size=8,
    compact=false,
)
    g = sankey_graph(s.args...)
    names = sankey_names(g, node_labels)
    if node_colors === nothing
        node_colors = palette(get(plotattributes, :color_palette, :default))
    end

    x, y, mask = sankey_layout!(g)
    perm = sortperm(y, rev=true)

    vw = vertex_weight.(Ref(g), vertices(g))
    m = maximum(vw)

    if compact == true
        y = make_compact(x, y, vw / m)
    end

    src_offsets = get_src_offsets(g, perm) ./ m
    dst_offsets = get_dst_offsets(g, perm) ./ m

    if label_position ∉ (:inside, :left, :right, :top, :bottom, :node, :legend)
        error("label_position :$label_position not supported")
    elseif label_position !== :legend
        xlabs = Float64[]
        ylabs = Float64[]
        lab_orientations = Symbol[]
    end

    for (i, v) in enumerate(vertices(g))
        h = vw[i] / (2m)

        if !(mask[i])
            @series begin
                seriestype := :shape
                label := label_position === :legend ? names[i] : ""
                fillcolor --> node_colors[mod1(i, end)]
                [x[i]-0.1, x[i]+0.1, x[i]+0.1, x[i]-0.1], [y[i]-h, y[i]-h, y[i]+h, y[i]+h]
            end

            for (j, w) in enumerate(vertices(g))
                if has_edge(g, v, w)
                    y_src = y[i] + h - src_offsets[i, j]
                    h_edge = g.weights[j, i] / (2m)

                    sankey_y = Float64[]
                    x_start = x[i] + 0.1
                    k = j
                    l = i
                    while mask[k]
                        y_dst = y[k] + vw[k] / (2m) - dst_offsets[k, l]
                        x_coords = range(0, 1, length=length(x_start:0.01:x[k]))
                        y_coords =
                            remap(1 ./ (1 .+ exp.(6 .* (1 .- 2 .* x_coords))), y_src, y_dst)

                        append!(sankey_y, y_coords)

                        x_start = x[k] + 0.01
                        y_src = y_dst
                        l = k
                        k = findfirst(==(first(outneighbors(g, k))), vertices(g))
                    end

                    y_dst = y[k] + vw[k] / (2m) - dst_offsets[k, l]
                    x_coords = range(0, 1, length=length(x_start:0.01:x[k]-0.1))
                    y_coords = remap(1 ./ (1 .+ exp.(6 .* (1 .- 2 .* x_coords))), y_src, y_dst)
                    append!(sankey_y, y_coords)
                    sankey_x = range(x[i]+0.1, x[k]-0.1, step=0.01)

                    @series begin
                        seriestype := :path
                        primary := false
                        linecolor := nothing
                        linewidth := false
                        fillrange := sankey_y .- 2h_edge
                        fillalpha --> 0.5
                        if edge_color === :gradient
                            fillcolor := getindex(
                                cgrad(node_colors[[i, k]]),
                                range(0, 1, length=length(sankey_x)),
                            )
                        elseif edge_color === :src
                            fillcolor := node_colors[i]
                        elseif edge_color === :dst
                            fillcolor := node_colors[k]
                        else
                            fillcolor := edge_color
                        end
                        sankey_x, sankey_y
                    end
                end
            end

            if label_position !== :legend
                xlab, orientation = if label_position in (:node, :top, :bottom)
                    olab = if label_position === :top
                        :bottom
                    elseif label_position === :bottom
                        :top
                    else
                        :center
                    end
                    x[i], olab
                elseif label_position === :right ||
                        (label_position === :inside && x[i] < maximum(x))
                    x[i] + 0.15, :left
                elseif label_position === :left ||
                        (label_position === :inside && x[i] == maximum(x))
                    x[i] - 0.15, :right
                else
                    error("label_position :$label_position not supported")
                end
                ylab = if label_position === :top
                    y[i] + h + 0.02
                elseif label_position === :bottom
                    y[i] - h - 0.0025 * label_size
                else
                    y[i]
                end
                push!(xlabs, xlab)
                push!(ylabs, ylab)
                push!(lab_orientations, orientation)
            end
        end
    end

    if label_position !== :legend
        @series begin
            primary := :false
            seriestype := :scatter
            markersize := 0
            series_annotations := text.(names, lab_orientations, label_size)
            xlabs, ylabs
        end

        # extend axes for labels
        if label_position in (:left, :right)
            x_extra = label_position === :left ? minimum(xlabs) - 0.4 : maximum(xlabs) + 0.5
            @series begin
                primary := false
                seriestype := :scatter
                markersize := 0
                [x_extra], [ylabs[1]]
            end
        end
    end

    primary := false
    framestyle --> :none
    legend --> label_position === :legend ? :outertopright : false
    ()
end

sankey_graph(src, dst, w) = SimpleWeightedDiGraph(src, dst, w)
sankey_graph(g::SimpleWeightedDiGraph) = copy(g)
sankey_graph(args...) = error("Check `?sankey` for supported signatures.")

sankey_names(g, names) = names
sankey_names(g, ::Nothing) = string.("Node", eachindex(vertices(g)))

function sankey_layout!(g)
    xs, ys, paths = solve_positions(Zarate(), g)
    mask = falses(length(xs))
    for (edge, path) in paths
        s = edge.src
        px, py = path
        if length(px) > 2
            for i in 2:length(px)-1
                add_vertex!(g)
                v = last(vertices(g))
                add_edge!(g, s, v, edge.weight)
                push!(xs, px[i])
                push!(ys, py[i])
                push!(mask, true)
                s = v
            end
            add_edge!(g, s, edge.dst, edge.weight)
            rem_edge!(g, edge)
        end
    end
    return xs, ys, mask
end

function vertex_weight(g, v)
    max(
        sum0(weight, Iterators.filter(e -> src(e) == v, edges(g))),
        sum0(weight, Iterators.filter(e -> dst(e) == v, edges(g))),
    )
end
sum0(f, x) = isempty(x) ? 0.0 : sum(f, x)

in_edges(g, v) = Iterators.filter(e -> dst(e) == v, edges(g))
out_edges(g, v) = Iterators.filter(e -> src(e) == v, edges(g))

function get_src_offsets(g, perm)
    verts = vertices(g)
    n = nv(g)
    p = spzeros(n, n)
    for (i, v) in enumerate(verts)
        offset = 0.0
        for j in perm
            if has_edge(g, v, verts[j])
                if offset > 0
                    p[i, j] = offset
                end
                offset += g.weights[j, i]
            end
        end
    end
    return p
end

function get_dst_offsets(g, perm)
    verts = vertices(g)
    n = nv(g)
    p = spzeros(n, n)
    for (i, v) in enumerate(verts)
        offset = 0.0
        for j in perm
            if has_edge(g, verts[j], i)
                if offset > 0
                    p[i, j] = offset
                end
                offset += g.weights[i, j]
            end
        end
    end
    return p
end

function remap(x, lo, hi)
    xlo, xhi = extrema(x)
    lo .+ (x .- xlo) / (xhi - xlo) * (hi - lo)
end

function make_compact(x, y, w)
    x = round.(Int, x)
    ux = unique(x)
    heights = zeros(length(ux))
    uinds = [findall(==(uxi), x) for uxi in ux]
    for (i, inds) in enumerate(uinds)
        perm = sortperm(view(y, inds))
        start = 0
        for j in inds[perm]
            y[j] = start + w[j] / 2
            start += 0.1 + w[j]
        end
        heights[i] = start - 0.1
    end
    maxh = maximum(heights)
    for (i, inds) in enumerate(uinds)
        y[inds] .+= (maxh - heights[i]) / 2
    end
    return y
end

end


text:
#!/usr/bin/env julia --project

## @article{10.1175/1520-0469(2000)057<1052:ALESSO>2.0.CO;2,
##     author = {Kosović, Branko and Curry, Judith A.},
##     title = "{A Large Eddy Simulation Study of a Quasi-Steady,
##               Stably Stratified Atmospheric Boundary Layer}",
##     journal = {Journal of the Atmospheric Sciences},
##     volume = {57},
##     number = {8},
##     pages = {1052-1068},
##     year = {2000},
##     month = {04},
##     issn = {0022-4928},
##     doi = {10.1175/1520-0469(2000)057<1052:ALESSO>2.0.CO;2},
##     url = {https://doi.org/10.1175/1520-0469(2000)057<1052:ALESSO>2.0.CO;2},
## }

## @article{doi:10.1029/2018MS001534,
## author = {Nishizawa, S. and Kitamura, Y.},
## title = {A Surface Flux Scheme Based on the Monin-Obukhov Similarity for Finite Volume Models},
## journal = {Journal of Advances in Modeling Earth Systems},
## volume = {10},
## number = {12},
## pages = {3159-3175},
## year = {2018}
## doi = {10.1029/2018MS001534},
## url = {https://agupubs.onlinelibrary.wiley.com/doi/abs/10.1029/2018MS001534},
## }

using ClimateMachine
ClimateMachine.init(parse_clargs = true)

using ClimateMachine.Atmos
using ClimateMachine.Orientations
using ClimateMachine.ConfigTypes
using ClimateMachine.DGMethods.NumericalFluxes
using ClimateMachine.Diagnostics
using ClimateMachine.GenericCallbacks
using ClimateMachine.Mesh.Filters
using ClimateMachine.Mesh.Grids
using ClimateMachine.ODESolvers
using ClimateMachine.Thermodynamics
using ClimateMachine.TurbulenceClosures
using ClimateMachine.VariableTemplates

using Distributions
using Random
using StaticArrays
using Test
using DocStringExtensions
using LinearAlgebra

using CLIMAParameters
using CLIMAParameters.Planet: R_d, cp_d, cv_d, MSLP, grav, day
struct EarthParameterSet <: AbstractEarthParameterSet end
const param_set = EarthParameterSet()

import ClimateMachine.Atmos: atmos_source!, flux_second_order!
using ClimateMachine.Atmos: altitude, recover_thermo_state

"""
  StableBL Geostrophic Forcing (Source)
"""
struct StableBLGeostrophic{FT} <: Source
    "Coriolis parameter [s⁻¹]"
    f_coriolis::FT
    "Eastward geostrophic velocity `[m/s]` (Base)"
    u_geostrophic::FT
    "Eastward geostrophic velocity `[m/s]` (Slope)"
    u_slope::FT
    "Northward geostrophic velocity `[m/s]`"
    v_geostrophic::FT
end
function atmos_source!(
    s::StableBLGeostrophic,
    atmos::AtmosModel,
    source::Vars,
    state::Vars,
    diffusive::Vars,
    aux::Vars,
    t::Real,
    direction,
)

    f_coriolis = s.f_coriolis
    u_geostrophic = s.u_geostrophic
    u_slope = s.u_slope
    v_geostrophic = s.v_geostrophic

    z = altitude(atmos, aux)
    # Note z dependence of eastward geostrophic velocity
    u_geo = SVector(u_geostrophic + u_slope * z, v_geostrophic, 0)
    ẑ = vertical_unit_vector(atmos, aux)
    fkvector = f_coriolis * ẑ
    # Accumulate sources
    source.ρu -= fkvector × (state.ρu .- state.ρ * u_geo)
    return nothing
end

"""
  StableBL Sponge (Source)
"""
struct StableBLSponge{FT} <: Source
    "Maximum domain altitude (m)"
    z_max::FT
    "Altitude at with sponge starts (m)"
    z_sponge::FT
    "Sponge Strength 0 ⩽ α_max ⩽ 1"
    α_max::FT
    "Sponge exponent"
    γ::FT
    "Eastward geostrophic velocity `[m/s]` (Base)"
    u_geostrophic::FT
    "Eastward geostrophic velocity `[m/s]` (Slope)"
    u_slope::FT
    "Northward geostrophic velocity `[m/s]`"
    v_geostrophic::FT
end
function atmos_source!(
    s::StableBLSponge,
    atmos::AtmosModel,
    source::Vars,
    state::Vars,
    diffusive::Vars,
    aux::Vars,
    t::Real,
    direction,
)

    z_max = s.z_max
    z_sponge = s.z_sponge
    α_max = s.α_max
    γ = s.γ
    u_geostrophic = s.u_geostrophic
    u_slope = s.u_slope
    v_geostrophic = s.v_geostrophic

    z = altitude(atmos, aux)
    u_geo = SVector(u_geostrophic + u_slope * z, v_geostrophic, 0)
    ẑ = vertical_unit_vector(atmos, aux)
    # Accumulate sources
    if z_sponge <= z
        r = (z - z_sponge) / (z_max - z_sponge)
        β_sponge = α_max * sinpi(r / 2)^s.γ
        source.ρu -= β_sponge * (state.ρu .- state.ρ * u_geo)
    end
    return nothing
end

"""
  Initial Condition for StableBoundaryLayer LES
"""
function init_problem!(problem, bl, state, aux, (x, y, z), t)
    # Problem floating point precision
    FT = eltype(state)
    R_gas::FT = R_d(bl.param_set)
    c_p::FT = cp_d(bl.param_set)
    c_v::FT = cv_d(bl.param_set)
    p0::FT = MSLP(bl.param_set)
    _grav::FT = grav(bl.param_set)
    γ::FT = c_p / c_v
    # Initialise speeds [u = Eastward, v = Northward, w = Vertical]
    u::FT = 8
    v::FT = 0
    w::FT = 0
    # Assign piecewise quantities to θ_liq and q_tot
    θ_liq::FT = 0
    q_tot::FT = 0
    # Piecewise functions for potential temperature and total moisture
    z1 = FT(100)
    if z <= z1
        θ_liq = FT(265)
    else
        θ_liq = FT(265) + FT(0.01) * (z - z1)
    end
    θ = θ_liq
    π_exner = FT(1) - _grav / (c_p * θ) * z # exner pressure
    ρ = p0 / (R_gas * θ) * (π_exner)^(c_v / R_gas) # density
    # Establish thermodynamic state and moist phase partitioning
    TS = PhaseEquil_ρθq(bl.param_set, ρ, θ_liq, q_tot)

    # Compute momentum contributions
    ρu = ρ * u
    ρv = ρ * v
    ρw = ρ * w

    # Compute energy contributions
    e_kin = FT(1 // 2) * (u^2 + v^2 + w^2)
    e_pot = _grav * z
    ρe_tot = ρ * total_energy(e_kin, e_pot, TS)

    # Assign initial conditions for prognostic state variables
    state.ρ = ρ
    state.ρu = SVector(ρu, ρv, ρw)
    state.ρe = ρe_tot
    state.moisture.ρq_tot = ρ * q_tot

    if z <= FT(50) # Add random perturbations to bottom 50m of model
        state.ρe += rand() * ρe_tot / 100
    end
end

function surface_temperature_variation(state, aux, t)
    FT = eltype(state)
    ρ = state.ρ
    q_tot = state.moisture.ρq_tot / ρ
    θ_liq_sfc = FT(265) - FT(1 / 4) * (t / 3600)
    TS = PhaseEquil_ρθq(param_set, ρ, θ_liq_sfc, q_tot)
    return air_temperature(TS)
end

function config_problem(::Type{FT}, N, resolution, xmax, ymax, zmax) where {FT}

    ics = init_problem!     # Initial conditions

    C_smag = FT(0.23)     # Smagorinsky coefficient
    C_drag = FT(0.001)    # Momentum exchange coefficient
    z_sponge = FT(300)     # Start of sponge layer
    α_max = FT(0.75)       # Strength of sponge layer (timescale)
    γ = 2                  # Strength of sponge layer (exponent)
    u_geostrophic = FT(8)        # Eastward relaxation speed
    u_slope = FT(0)              # Slope of altitude-dependent relaxation speed
    v_geostrophic = FT(0)        # Northward relaxation speed
    f_coriolis = FT(1.39e-4) # Coriolis parameter
    q_sfc = FT(0)
    u_star = FT(0.30)

    # Assemble source components
    source = (
        Gravity(),
        StableBLSponge{FT}(
            zmax,
            z_sponge,
            α_max,
            γ,
            u_geostrophic,
            u_slope,
            v_geostrophic,
        ),
        StableBLGeostrophic{FT}(
            f_coriolis,
            u_geostrophic,
            u_slope,
            v_geostrophic,
        ),
    )

    # Choose default IMEX solver
    ode_solver_type = ClimateMachine.ExplicitSolverType()

    # Set up problem initial and boundary conditions
    moisture_flux = FT(0)
    problem = AtmosProblem(
        init_state_prognostic = ics,
        boundarycondition = (
            AtmosBC(
                momentum = Impenetrable(DragLaw(
                    # normPu_int is the internal horizontal speed
                    # P represents the projection onto the horizontal
                    (state, aux, t, normPu_int) -> (u_star / normPu_int)^2,
                )),
                energy = BulkFormulaEnergy(
                    (state, aux, t, normPu_int) -> C_drag,
                    (state, aux, t) -> (
                        surface_temperature_variation(state, aux, t),
                        q_sfc,
                    ),
                ),
                moisture = BulkFormulaMoisture(
                    (state, aux, t, normPu_int) -> C_drag,
                    (state, aux, t) -> q_sfc,
                ),
            ),
            AtmosBC(),
        ),
    )

    # Assemble model components
    model = AtmosModel{FT}(
        AtmosLESConfigType,
        param_set;
        problem = problem,
        turbulence = SmagorinskyLilly{FT}(C_smag),
        moisture = EquilMoist{FT}(; maxiter = 5, tolerance = FT(0.1)),
        source = source,
    )

    # Assemble configuration
    config = ClimateMachine.AtmosLESConfiguration(
        "StableBoundaryLayer",
        N,
        resolution,
        xmax,
        ymax,
        zmax,
        param_set,
        init_problem!,
        solver_type = ode_solver_type,
        model = model,
    )
    return config
end

function config_diagnostics(driver_config)
    default_dgngrp = setup_atmos_default_diagnostics(
        AtmosLESConfigType(),
        "2500steps",
        driver_config.name,
    )
    core_dgngrp = setup_atmos_core_diagnostics(
        AtmosLESConfigType(),
        "2500steps",
        driver_config.name,
    )
    return ClimateMachine.DiagnosticsConfiguration([
        default_dgngrp,
        core_dgngrp,
    ])
end

function main()
    FT = Float64

    # DG polynomial order
    N = 4
    # Domain resolution and size
    Δh = FT(20)
    Δv = FT(20)

    resolution = (Δh, Δh, Δv)

    # Prescribe domain parameters
    xmax = FT(100)
    ymax = FT(100)
    zmax = FT(400)

    t0 = FT(0)

    # Required simulation time == 9hours
    timeend = FT(3600 * 0.1)
    CFLmax = FT(0.4)

    driver_config = config_problem(FT, N, resolution, xmax, ymax, zmax)
    solver_config = ClimateMachine.SolverConfiguration(
        t0,
        timeend,
        driver_config,
        init_on_cpu = true,
        Courant_number = CFLmax,
    )
    dgn_config = config_diagnostics(driver_config)

    cbtmarfilter = GenericCallbacks.EveryXSimulationSteps(1) do
        Filters.apply!(
            solver_config.Q,
            ("moisture.ρq_tot",),
            solver_config.dg.grid,
            TMARFilter(),
        )
        nothing
    end

    # State variable
    Q = solver_config.Q
    # Volume geometry information
    vgeo = driver_config.grid.vgeo
    M = vgeo[:, Grids._M, :]
    # Unpack prognostic vars
    ρ₀ = Q.ρ
    ρe₀ = Q.ρe
    # DG variable sums
    Σρ₀ = sum(ρ₀ .* M)
    Σρe₀ = sum(ρe₀ .* M)
    cb_check_cons = GenericCallbacks.EveryXSimulationSteps(3000) do
        Q = solver_config.Q
        δρ = (sum(Q.ρ .* M) - Σρ₀) / Σρ₀
        δρe = (sum(Q.ρe .* M) .- Σρe₀) ./ Σρe₀
        @show (abs(δρ))
        @show (abs(δρe))
        nothing
    end

    result = ClimateMachine.invoke!(
        solver_config;
        diagnostics_config = dgn_config,
        user_callbacks = (cbtmarfilter, cb_check_cons),
        check_euclidean_distance = true,
    )
end

main()

max_stars_repo_name: mwarusz/CLIMA

max_issues_repo_name: mwarusz/CLIMA

max_forks_repo_name: mwarusz/CLIMA


text:
function create_sfactor(w::Direction, s::DerivativeDirection, g::Grid, ω::Float ; m=3.5, lnR=-12)
    (ϵ₀, μ₀, c₀) = normalize_parameters(g);
    dw = dh(g, w);
    Tw = g.Npml[Int(w)]*dw;
    Nw = size(g, Int(w));
    Nw_pml = g.Npml[Int(w)];

    σmax = -(m+1)*lnR/(2*η₀*Tw);
    σw = l -> σmax*(l/Tw)^m;
    S = l -> 1-1im*σw(l)/(ω*ϵ₀);

    sfactor_array = ones(ComplexF64, Nw);

    for i in 1:Nw
        if s == Forward
            if i <= Nw_pml
                sfactor_array[i] = S(dw*(Nw_pml-i+0.5));
            elseif i > Nw - Nw_pml
                sfactor_array[i] = S(dw*(i-(Nw-Nw_pml)-0.5));
            end
        elseif s == Backward
            if i <= Nw_pml
                sfactor_array[i] = S(dw*(Nw_pml-i+1));
            elseif i > Nw - Nw_pml
                sfactor_array[i] = S(dw*(i-(Nw-Nw_pml)-1));
            end
        end
    end

    return sfactor_array
end

function S_create(g::Grid, ω::Float)
    # Create the sfactor in each direction and for forward and backward
    s_vector_x_f = create_sfactor(x̂, Forward,  g, ω);
    s_vector_x_b = create_sfactor(x̂, Backward, g, ω);
    s_vector_y_f = create_sfactor(ŷ, Forward,  g, ω);
    s_vector_y_b = create_sfactor(ŷ, Backward, g, ω);

    # Fill the 2D space with layers of appropriate s-factors
    Sx_f_2D = zeros(ComplexF64, size(g));
    Sx_b_2D = zeros(ComplexF64, size(g));
    Sy_f_2D = zeros(ComplexF64, size(g));
    Sy_b_2D = zeros(ComplexF64, size(g));

    for i = 1:size(g, 1)
        Sy_f_2D[i, :] = s_vector_y_f.^-1;
        Sy_b_2D[i, :] = s_vector_y_b.^-1;
    end

    for i = 1:size(g, 2)
        Sx_f_2D[:, i] = s_vector_x_f.^-1;
        Sx_b_2D[:, i] = s_vector_x_b.^-1;
    end

    # Construct the 1D total s-array into a diagonal matrix
    Sx_f = sparse(Diagonal(Sx_f_2D[:]));
    Sx_b = sparse(Diagonal(Sx_b_2D[:]));
    Sy_f = sparse(Diagonal(Sy_f_2D[:]));
    Sy_b = sparse(Diagonal(Sy_b_2D[:]));

    return (Sx_f, Sx_b, Sy_f, Sy_b)
end


text:
using DelimitedFiles

const program = readdlm("day08.txt")

function process(inst::Array, acc::Int, idx::Int)
    if inst[1] == "nop"
        idx += 1
    elseif inst[1] == "acc"
        acc += inst[2]
        idx += 1
    else
        idx += inst[2]
    end
    return acc, idx
end

function loop(program::Array)
    acc = 0
    idx = 1
    seen = Set{Int}()
    while !(idx in seen)
        push!(seen, idx)
        acc, idx = process(program[idx, :], acc, idx)
    end
    return acc
end

println("Part 1: $(loop(program))")

function alterprogram(program::Array{Any,2}, start::Int)
    new_program = copy(program)
    # find next instruction to modify
    # We need to provide a CartesianIndex but we are only interested in the first axis
    new_start = findnext(x -> x != "acc", program, CartesianIndex(start + 1, 1))[1]
    changes = Dict("nop" => "jmp", "jmp" => "nop")
    new_program[new_start] = changes[new_program[new_start]]
    return new_program, new_start
end

"""
return the acc value before exiting or looping and a flag False if looped or True if
exited normally.
"""
function loop_or_exit(program::Array)
    acc = 0
    idx = 1
    seen = Set{Int}()
    exited = false
    while !(idx in seen)
        push!(seen, idx)
        if idx > size(program)[1]
            exited = true
            break
        end
        acc, idx = process(program[idx, :], acc, idx)
    end
    return acc, exited
end

function testallprograms(program::Array)
    start = 1
    exited = false
    acc = 0
    while !exited
        new_program, start = alterprogram(program, start)
        acc, exited = loop_or_exit(new_program)
    end
    return acc
end

println("Part 2: $(testallprograms(program))")


text:
using Documenter
using ErdosExtras

makedocs(
    modules = [ErdosExtras]
    )

deploydocs(
    deps = Deps.pip("pygments", "mkdocs", "mkdocs-material", "python-markdown-math"),
    repo   = "github.com/CarloLucibello/ErdosExtras.jl.git",
    julia  = "1.0"
    )


text:
#
# Estimate-pi.jl
#
sumsq(x,y) = x*x + y*y;

N=1000000;
x = 0;
for i = 1:N
  if sumsq(rand(), rand()) < 1.0
    x += 1
  end
end
@printf "Estimate of PI for %d trials is %8.5f\n" N 4.0*(x / N);

max_stars_repo_path: Module 3/Chapter01/estimate-pi.jl

max_issues_repo_path: Module 3/Chapter01/estimate-pi.jl

max_forks_repo_path: Module 3/Chapter01/estimate-pi.jl


text:
using Documenter
using Literate
using Vulkan

push!(LOAD_PATH, joinpath(dirname(@__DIR__), "generator"))

using VulkanGen

function julia_files(dir)
    files = reduce(vcat, [joinpath(root, file) for (root, dirs, files) in walkdir(dir) for file in files])
    sort(filter(endswith(".jl"), files))
end

function replace_edit(content)
    haskey(ENV, "JULIA_GITHUB_ACTIONS_CI") && return content
    # Linking does not work locally, but we can make
    # the warning go away with a hard link to the repo.
    replace(
        content,
        r"EditURL = \".*<unknown>/(.*)\"" => s"EditURL = \"https://github.com/JuliaGPU/Vulkan.jl/tree/master/\1\"",
    )
end

function generate_markdowns()
    dir = joinpath(@__DIR__, "src")
    Threads.@threads for file in julia_files(dir)
        Literate.markdown(
            file,
            dirname(file);
            postprocess = replace_edit,
            documenter = true,
        )
    end
end

generate_markdowns()

makedocs(;
    modules=[Vulkan, VulkanGen],
    format=Documenter.HTML(prettyurls = true),
    pages=[
        "Home" => "index.md",
        "Introduction" => "intro.md",
        "Glossary" => "glossary.md",
        "Tutorial" => [
            "Getting started" => "tutorial/getting_started.md",
            "Error handling" => "tutorial/error_handling.md",
            "Resource management" => "tutorial/resource_management.md",
            "In-depth tutorial" => "tutorial/indepth.md",
            "Running compute shaders" => "tutorial/minimal_working_compute.md",
        ],
        "How to" => [
            "Specify package options" => "howto/preferences.md",
            "Debug an application" => "howto/debugging.md",
            "Manipulate handles" => "howto/handles.md",
            "Compile a SPIR-V shader from Julia" => "howto/shaders.md",
        ],
        "Reference" => [
            "Wrapper types" => "reference/wrapper_types.md",
            "Wrapper functions" => "reference/wrapper_functions.md",
            "API function dispatch" => "reference/dispatch.md",
            "Package options" => "reference/options.md",
        ],
        "Explanations" => [
            "Motivations" => "about/motivations.md",
            "Extension mechanism" => "about/extension_mechanism.md",
        ],
        "API" => "api.md",
        "Utility" => "utility.md",
        "Troubleshooting" => "troubleshooting.md",
        "Developer documentation" => [
            "Overview" => "dev/overview.md",
            "Vulkan specification" => "dev/spec.md",
            "Generator" => "dev/gen.md",
            "Next chains" => "dev/next_chains.md",
        ],
    ],
    repo="https://github.com/JuliaGPU/Vulkan.jl/blob/{commit}{path}#L{line}",
    sitename="Vulkan.jl",
    authors="serenity4 <cedric.bel@hotmail.fr>",
    strict=false,
    doctest=false,
    checkdocs=:exports,
    linkcheck=:false,
)

deploydocs(
    repo = "github.com/JuliaGPU/Vulkan.jl.git",
    push_preview = true,
)

max_stars_repo_name: JuliaGPU/VulkanAbstraction

max_issues_repo_name: JuliaGPU/VulkanAbstraction

max_forks_repo_name: JuliaGPU/VulkanAbstraction


text:
using Random

abstract type Modifier end

struct Modifiers <: Modifier
    modifiers::Array{Modifier}
end

Modifiers(modifiers::Modifier ...) = Modifiers([modifiers...])

function apply!(modifiers::Modifiers, data::DataFrame, event::Event)
    for modifier in modifiers.modifiers
        apply!(modifier, data, event)
    end

    return nothing
end

struct NoiseModifier <: Modifier
    seed::Random.MersenneTwister
    attribute::Float64 # The fraction of attribute values to disturb
end

NoiseModifier(attribute::Float64, seed::Int) = NoiseModifier(Random.seed!(seed), attribute)
NoiseModifier(attribute::Float64) = NoiseModifier(Random.default_rng(), attribute)

function apply!(modifier::NoiseModifier, data::DataFrame, event::Event)
    return nothing
end

struct FilterModifier <: Modifier
    columns::Array{Symbol}

    function FilterModifier(columns)
        _columns = unique(columns)
        if length(_columns) != length(columns) @warn "There are duplicate columns." end
        return new(_columns)
    end
end

FilterModifier(columns::Symbol...) = FilterModifier([columns...])

function apply!(modifier::FilterModifier, data::DataFrame, event::Event)
    columns = Symbol[]
    for col in modifier.columns
        if !(col in propertynames(data))
            throw(ArgumentError("stream doesn't have the column $col"))
        else
            push!(columns, col)
        end
    end

    select!(data, columns)
    return nothing
end

struct AlterDataModifier <: Modifier
    alter!::Function
end

function apply!(modifier::AlterDataModifier, data::DataFrame, event::Event)
    modifier.alter!(data, event)

    return nothing
end


text:
using LinearAlgebra, Statistics

import MPI, LinearAlgebra

function montecarlo(mc_eval::Function, mc_monitor::Function,
                    comm::MPI.Comm, # MPI communicator
                    n_evals::Integer, # total number of evaluations
                    n_returns::Integer, # number of return values per evaluation
                    batchsize::Integer=1) # transmission size
    rank = MPI.Comm_rank(comm)
    commsize = MPI.Comm_size(comm)

    @assert commsize > 1

    # bookkeeping
    if mod(n_evals, commsize-1) != 0
        if rank == 0
            println("Error (montecarlo):")
            println("    n_evals=$n_evals, commsize=$commsize")
            println("Choose commsize so that n_evals/(commsize-1) is integer")
        end
        MPI.Finalize()
        exit(1)
    end
    n_pernode = div(n_evals, commsize-1)

    # number of evaluations per worker
    if mod(n_pernode, batchsize) != 0
        batchsize = div(n_pernode, commsize-1)
    end
    if mod(n_pernode, batchsize) != 0
        if rank == 0
            println("Error (montecarlo):")
            println("    n_pernode=$n_pernode, batchsize=$batchsize")
            println("Choose commsize so that n_pernode/batchsize is integer")
        end
        MPI.Finalize()
        exit(1)
    end

    if rank > 0
        # workers
        contrib = zeros(batchsize, n_returns)
        @inbounds for i = 1:div(n_pernode, batchsize)
            # do work
            for j = 1:batchsize
                contrib[j,:] .= mc_eval()
            end
            MPI.Send(contrib, 0, 0, comm)
        end
    else
        # manager
        results = zeros(n_evals, n_returns)
        contrib = zeros(batchsize, n_returns)
        for nresults = 1:div(n_evals, batchsize)
            MPI.Recv!(contrib, MPI.ANY_SOURCE, 0, comm)
            results[nresults*batchsize-batchsize+1:nresults*batchsize,:] = contrib
            mc_monitor(nresults*batchsize, results)
        end
    end
end

function pi_wrapper()
    4.0 * (norm(rand(2)) < 1)
end

# this function reports intermediate results during MC runs
function pi_monitor(sofar, results)
    # Examine results every 12.5*10^5 draws
    if mod(sofar, 10^6) == 0
        m = mean(results[1:sofar,:],dims=1)
        println("reps: $sofar, pihat: $m")
    end
end

# do the monte carlo: 10^6 reps of single draws
function main()
    MPI.Init()
    n_evals = 10^7 # desired number of MC reps
    n_returns = 1
    batchsize = 10^5
    montecarlo(pi_wrapper, pi_monitor,
               MPI.COMM_WORLD, n_evals, n_returns, batchsize)
    MPI.Finalize()
end

main()

max_stars_repo_name: terasakisatoshi/juliaExer

max_issues_repo_name: terasakisatoshi/juliaExer

max_forks_repo_name: terasakisatoshi/juliaExer


text:
module ManifestTests

using  Test, UUIDs, Dates, TOML
import ..Pkg, LibGit2
using  ..Utils

# used with the reference manifests in `test/manifest/formats`
# ensures the manifests are valid and restored after test
function reference_manifest_isolated_test(f, dir::String; v1::Bool=false)
    env_dir = joinpath(@__DIR__, "manifest", "formats", dir)
    env_manifest = joinpath(env_dir, "Manifest.toml")
    env_project = joinpath(env_dir, "Project.toml")
    cp(env_manifest, string(env_manifest, "_backup"))
    cp(env_project, string(env_project, "_backup"))
    try
        isfile(env_manifest) || error("Reference manifest is missing")
        if Base.is_v1_format_manifest(Base.parsed_toml(env_manifest)) == !v1
            error("Reference manifest file at $(env_manifest) is invalid")
        end
        isolate(loaded_depot=true) do
            f(env_dir, env_manifest)
        end
    finally
        cp(string(env_manifest, "_backup"), env_manifest, force = true)
        rm(string(env_manifest, "_backup"))
        cp(string(env_project, "_backup"), env_project, force = true)
        rm(string(env_project, "_backup"))
    end
end

##

@testset "Manifest.toml formats" begin
    @testset "Default manifest format is v2" begin
        isolate(loaded_depot=true) do
            io = IOBuffer()
            Pkg.activate(; io=io, temp=true)
            output = String(take!(io))
            @test occursin(r"Activating.*project at.*", output)
            Pkg.add("Profile")
            env_manifest = Pkg.Types.Context().env.manifest_file
            @test Base.is_v1_format_manifest(Base.parsed_toml(env_manifest)) == false
            @test Pkg.Types.Context().env.manifest.manifest_format == v"2.0.0"
        end
    end

    @testset "v1.0: activate, change, maintain manifest format" begin
        reference_manifest_isolated_test("v1.0", v1 = true) do env_dir, env_manifest
            io = IOBuffer()
            Pkg.activate(env_dir; io=io)
            output = String(take!(io))
            @test occursin(r"Activating.*project at.*`.*v1.0`", output)
            @test Base.is_v1_format_manifest(Base.parsed_toml(env_manifest))

            Pkg.add("Profile")
            @test Base.is_v1_format_manifest(Base.parsed_toml(env_manifest))

            Pkg.rm("Profile")
            @test Base.is_v1_format_manifest(Base.parsed_toml(env_manifest))
        end
    end

    @testset "v2.0: activate, change, maintain manifest format" begin
        reference_manifest_isolated_test("v2.0") do env_dir, env_manifest
            io = IOBuffer()
            Pkg.activate(env_dir; io=io)
            output = String(take!(io))
            @test occursin(r"Activating.*project at.*`.*v2.0`", output)
            @test Base.is_v1_format_manifest(Base.parsed_toml(env_manifest)) == false

            Pkg.add("Profile")
            @test Base.is_v1_format_manifest(Base.parsed_toml(env_manifest)) == false

            Pkg.rm("Profile")
            @test Base.is_v1_format_manifest(Base.parsed_toml(env_manifest)) == false

            m = Pkg.Types.read_manifest(env_manifest)
            @test m.other["some_other_field"] == "other"
            @test m.other["some_other_data"] == [1,2,3,4]

            mktemp() do path, io
                Pkg.Types.write_manifest(io, m)
                m2 = Pkg.Types.read_manifest(env_manifest)
                @test m.deps == m2.deps
                @test m.julia_version == m2.julia_version
                @test m.manifest_format == m2.manifest_format
                @test m.other == m2.other
            end
        end
    end

    @testset "v3.0: unknown format, warn" begin
        # the reference file here is not actually v3.0. It just represents an unknown manifest format
        reference_manifest_isolated_test("v3.0_unknown") do env_dir, env_manifest
            io = IOBuffer()
            @test_logs (:warn,) Pkg.activate(env_dir; io=io)
        end
    end

    @testset "Pkg.upgrade_manifest()" begin
        reference_manifest_isolated_test("v1.0", v1 = true) do env_dir, env_manifest
            io = IOBuffer()
            Pkg.activate(env_dir; io=io)
            output = String(take!(io))
            @test occursin(r"Activating.*project at.*`.*v1.0`", output)
            @test Base.is_v1_format_manifest(Base.parsed_toml(env_manifest))

            Pkg.upgrade_manifest()
            @test Base.is_v1_format_manifest(Base.parsed_toml(env_manifest)) == false
            Pkg.activate(env_dir; io=io)
            output = String(take!(io))
            @test occursin(r"Activating.*project at.*`.*v1.0`", output)
            @test Pkg.Types.Context().env.manifest.manifest_format == v"2.0.0"
        end
    end
end

@testset "Manifest metadata" begin
    @testset "julia_version" begin
        @testset "new environment: value is `nothing`, then `VERSION` after resolve" begin
            isolate(loaded_depot=true) do
                Pkg.activate(; temp=true)
                @test Pkg.Types.Context().env.manifest.julia_version == nothing
                Pkg.add("Profile")
                @test Pkg.Types.Context().env.manifest.julia_version == VERSION
            end
        end
        @testset "activating old environment: maintains old version, then `VERSION` after resolve" begin
            reference_manifest_isolated_test("v2.0") do env_dir, env_manifest
                Pkg.activate(env_dir)
                @test Pkg.Types.Context().env.manifest.julia_version == v"1.7.0-DEV.1199"

                Pkg.add("Profile")
                @test Pkg.Types.Context().env.manifest.julia_version == VERSION
            end
        end
        @testset "instantiate manifest from different julia_version" begin
            reference_manifest_isolated_test("v1.0", v1 = true) do env_dir, env_manifest
                Pkg.activate(env_dir)
                @test_logs (:warn, r"The active manifest file") Pkg.instantiate()
                @test Pkg.Types.Context().env.manifest.julia_version == nothing
            end
            if VERSION >= v"1.8"
                reference_manifest_isolated_test("v2.0") do env_dir, env_manifest
                    Pkg.activate(env_dir)
                    @test_logs (:warn, r"The active manifest file") Pkg.instantiate()
                    @test Pkg.Types.Context().env.manifest.julia_version == v"1.7.0-DEV.1199"
                end
            end
        end
    end
end

end # module


text:
using Millboard
using Base.Test

macro test_colored(expr)
  #println()
  #println(expr.args[3])
  #println(expr.args[1])
  :( @test $expr )
end

cyan = colored(:cyan)
@test_colored "\e[1m\e[36m1\e[0m\e[1m" == 1 |> cyan |> string
@test_colored "\e[1m\e[36mπ = 3.1415926535897...\e[0m\e[1m" == pi |> cyan |> string

green = colored(:green)
@test_colored "Millboard.Coat(:green,sin)" == sin |> green |> string

magenta = colored(:magenta)
@test_colored "\e[1m\e[35mABC\e[0m\e[1m" == "ABC" |> magenta |> string
@test_colored "\e[1m\e[35mABC\e[0m\e[1m" == "ABC" |> magenta(s->s) |> string
@test_colored "\e[1m\e[1m\e[35mABC\e[0m\e[1m\e[1m" == "ABC" |> magenta(s->s[1:end]) |> string

@test_colored "\e[1m\e[1m\e[35mA\e[0m\e[1m\e[1mBC" == "ABC" |> magenta(s->s[]) |> string
@test_colored "\e[1m\e[1m\e[35mA\e[0m\e[1m\e[1mBC" == "ABC" |> magenta(s->s[1]) |> string

@test_colored "\e[1mA\e[1m\e[35mB\e[0m\e[1m\e[1mC" == "ABC" |> magenta(s->s[2]) |> string
@test_colored "\e[1mAB\e[1m\e[35mC\e[0m\e[1m\e[1m" == "ABC" |> magenta(s->s[3]) |> string
@test_colored "\e[1mAB\e[1m\e[35mC\e[0m\e[1m\e[1m" == "ABC" |> magenta(s->s[end]) |> string

str = "Lorem ipsum dolor sit amet"
@test_colored "\e[1m\e[33mLorem ipsum dolor sit amet\e[0m\e[1m" == str |> colored(:yellow) |> string
@test_colored "\e[1m\e[33mLorem ipsum dolor sit amet\e[0m\e[1m" == str |> colored(:yellow, s->s) |> string
@test_colored "\e[1m\e[33mLorem ipsum dolor sit amet\e[0m\e[1m" == colored(:yellow, s->s, str) |> string
@test_colored "\e[1mLorem \e[1m\e[33mipsum dol\e[0m\e[1m\e[1mor sit amet" == colored(:yellow, s->s[7:15], str) |> string

# test colored tables

colnames = [:blue :green :red :magenta]
board = map(colnames) do name
  colored(name, s->s[1:end], string(name))
end
@test_colored "+===+=======+========+======+==========+\n|   | :blue | :green | :red | :magenta |\n+===+=======+========+======+==========+\n| 1 |  \e[1m\e[1m\e[34mblue\e[0m\e[1m\e[1m |  \e[1m\e[1m\e[32mgreen\e[0m\e[1m\e[1m |  \e[1m\e[1m\e[31mred\e[0m\e[1m\e[1m |  \e[1m\e[1m\e[35mmagenta\e[0m\e[1m\e[1m |\n+---+-------+--------+------+----------+" == table(board, :colnames=>colnames) |> string

colored_colnames = [colored(:blue, "FIRST")]
@test_colored "+===+=======+=======+=====+=========+\n|   | \e[1m\e[34mFIRST\e[0m\e[1m |     2 |   3 |       4 |\n+===+=======+=======+=====+=========+\n| 1 |  \e[1m\e[1m\e[34mblue\e[0m\e[1m\e[1m | \e[1m\e[1m\e[32mgreen\e[0m\e[1m\e[1m | \e[1m\e[1m\e[31mred\e[0m\e[1m\e[1m | \e[1m\e[1m\e[35mmagenta\e[0m\e[1m\e[1m |\n+---+-------+-------+-----+---------+" == table(board, :colnames=>colored_colnames) |>string

board = [5 6; 7 8]
@test_colored "+===+===+\n|   | 1 |\n+===+===+\n| 1 | \e[1m\e[36m5\e[0m\e[1m |\n+---+---+\n| 2 | \e[1m\e[36m7\e[0m\e[1m |\n+---+---+\n| 3 | \e[1m\e[36m6\e[0m\e[1m |\n+---+---+\n| 4 | \e[1m\e[36m8\e[0m\e[1m |\n+---+---+" == board |> cyan(b->b[:]) |> table |> string
@test_colored "+===+===+\n|   | 1 |\n+===+===+\n| 1 | \e[1m\e[36m5\e[0m\e[1m |\n+---+---+\n| 2 | \e[1m\e[36m7\e[0m\e[1m |\n+---+---+\n| 3 | \e[1m\e[36m6\e[0m\e[1m |\n+---+---+\n| 4 | \e[1m\e[36m8\e[0m\e[1m |\n+---+---+" == board |> cyan(b->b[:,]) |> table |> string

@test_colored "+===+===+===+\n|   | 1 | 2 |\n+===+===+===+\n| 1 | \e[1m\e[36m5\e[0m\e[1m | \e[1m\e[36m6\e[0m\e[1m |\n+---+---+---+\n| 2 | \e[1m\e[36m7\e[0m\e[1m | \e[1m\e[36m8\e[0m\e[1m |\n+---+---+---+" == board |> cyan(b->b[:,:]) |> table |> string
@test_colored "+===+===+===+\n|   | 1 | 2 |\n+===+===+===+\n| 1 | \e[1m\e[36m5\e[0m\e[1m | 6 |\n+---+---+---+\n| 2 | \e[1m\e[36m7\e[0m\e[1m | 8 |\n+---+---+---+" == board |> cyan(b->b[1:2]) |> table |> string
@test_colored "+===+===+===+\n|   | 1 | 2 |\n+===+===+===+\n| 1 | \e[1m\e[36m5\e[0m\e[1m | \e[1m\e[36m6\e[0m\e[1m |\n+---+---+---+\n| 2 | \e[1m\e[36m7\e[0m\e[1m | \e[1m\e[36m8\e[0m\e[1m |\n+---+---+---+" == board |> cyan(b->b[1:end]) |> table |> string
@test_colored "+===+===+===+\n|   | 1 | 2 |\n+===+===+===+\n| 1 | \e[1m\e[36m5\e[0m\e[1m | 6 |\n+---+---+---+\n| 2 | 7 | 8 |\n+---+---+---+" == board |> cyan(b->b[1:1]) |> table |> string
@test_colored "+===+===+===+\n|   | 1 | 2 |\n+===+===+===+\n| 1 | 5 | 6 |\n+---+---+---+\n| 2 | 7 | 8 |\n+---+---+---+" == board |> cyan(b->b[2:1]) |> table |> string
@test_colored "+===+===+===+\n|   | 1 | 2 |\n+===+===+===+\n| 1 | 5 | 6 |\n+---+---+---+\n| 2 | \e[1m\e[36m7\e[0m\e[1m | 8 |\n+---+---+---+" == board |> cyan(b->b[2,1]) |> table |> string
@test_colored "+===+===+===+\n|   | 1 | 2 |\n+===+===+===+\n| 1 | 5 | \e[1m\e[36m6\e[0m\e[1m |\n+---+---+---+\n| 2 | \e[1m\e[36m7\e[0m\e[1m | \e[1m\e[36m8\e[0m\e[1m |\n+---+---+---+" == board |> cyan(b->b[2:end]) |> table |> string
@test_colored "+===+===+===+\n|   | 1 | 2 |\n+===+===+===+\n| 1 | 5 | 6 |\n+---+---+---+\n| 2 | \e[1m\e[36m7\e[0m\e[1m | \e[1m\e[36m8\e[0m\e[1m |\n+---+---+---+" == board |> cyan(b->b[2,:]) |> table |> string

board = ([1 2], [5 6;7 8], [9 10; 11 12])
@test_colored "+===+=====+=====+=======+\n|   |   1 |   2 |     3 |\n+===+=====+=====+=======+\n|   | \e[1m\e[36m1 2\e[0m\e[1m |\e[1m\e[36m 5 6\e[0m\e[1m |\e[1m\e[36m  9 10\e[0m\e[1m |\n| 1 |     |\e[1m\e[36m 7 8\e[0m\e[1m |\e[1m\e[36m 11 12\e[0m\e[1m |\n+---+-----+-----+-------+" == board |> cyan |> table |> string
@test_colored "+===+=====+=====+=======+\n|   |   1 |   2 |     3 |\n+===+=====+=====+=======+\n|   | \e[1m\e[36m1 2\e[0m\e[1m |\e[1m\e[36m 5 6\e[0m\e[1m |\e[1m\e[36m  9 10\e[0m\e[1m |\n| 1 |     |\e[1m\e[36m 7 8\e[0m\e[1m |\e[1m\e[36m 11 12\e[0m\e[1m |\n+---+-----+-----+-------+" == board |> cyan(b->b) |> table |> string
@test_colored "+===+=====+=====+=======+\n|   |   1 |   2 |     3 |\n+===+=====+=====+=======+\n|   | \e[1m\e[36m1 2\e[0m\e[1m |\e[1m\e[36m 5 6\e[0m\e[1m |\e[1m\e[36m  9 10\e[0m\e[1m |\n| 1 |     |\e[1m\e[36m 7 8\e[0m\e[1m |\e[1m\e[36m 11 12\e[0m\e[1m |\n+---+-----+-----+-------+" == board |> cyan(b->b[:,:]) |> table |> string
@test_colored "+===+=====+=====+=======+\n|   |   1 |   2 |     3 |\n+===+=====+=====+=======+\n|   | \e[1m\e[36m1 2\e[0m\e[1m | 5 6 |  9 10 |\n| 1 |     | 7 8 | 11 12 |\n+---+-----+-----+-------+" == board |> cyan(b->b[]) |> table |> string

# coverage
import Millboard: Coating, coating
A = Coating(:green, [1,2,3])
buf = IOBuffer()
show(buf, A)
coating(A, AbstractArray)
coating(A, String)

import Millboard: ANSIEscaped, ANSIEscapedString
E = ANSIEscaped(:green, [1,2,3])
endof(E)
next(E, true)

s = ANSIEscapedString("hello", 5, 5)
endof(s)
next(s, 1)


text:
# https://github.com/alan-turing-institute/MLJ.jl/issues/540
@testset "Levels" begin
    X, y = MLJBase.@load_iris
    df = DataFrame(hcat(DataFrame(X), y))
    X = filter(r -> r.x1 != "virginica", df)
    y = X.x1
    X = DataFrames.select(X, Not(:x1))
    mdl = LogisticClassifier()
    mach = MLJBase.machine(mdl, X, y)
    MLJBase.fit!(mach)
    p = MLJBase.predict(mach, rows=1:2)[1]
    @test typeof(p) <: MLJBase.UnivariateFinite

    X, y = MLJBase.@load_iris
    mach = MLJBase.machine(mdl, X, y)
    MLJBase.fit!(mach)
    p = MLJBase.predict(mach, rows=1:2)[1]
    @test typeof(p) <: MLJBase.UnivariateFinite
end


text:
__precompile__(true)
module WorldMapData

using Printf
using Shapefile

const VALID_RESOLUTIONS = [
    :FULL,
    :HIGH,
    :INTERMEDIATE,
    :LOW,
    :CRUDE
]

export valid_resolution
"""
Return whether the input symbol is a valid resolution for GSHHS data.

Arguments:
    - `res::Symbol` Symbol to validate

Returns:
    - `valid::Bool` Boolean of resolution validity.
"""
function valid_resolution(res::Symbol)
    return res in VALID_RESOLUTIONS ? true : false
end

const RESOLUTION_LOOKUP = Dict(
    :FULL => "f",
    :HIGH => "h",
    :INTERMEDIATE => "i",
    :LOW => "l",
    :CRUDE => "c",
)

"""
Ensure that the shapefile exists at the specified path otherwise unzip the
base file.

Arguments:
    - `path` Path to datafile
"""
function load_data(path)

    shp_path = path * ".shp"

    # if !isfile(shp_path)
    #     unzip(zip_path)
    # end

    # @debug "Loading datafile: $shp_path"
    handle = open(shp_path, "r") do io
        read(io, Shapefile.Handle)
    end

    return handle
end

export load_GSHHS_data
"""
Load GSHHS Shapefile data and return all polygons in the data set.

All data sets come in 5 different resolutions:
	- `:FULL`: Full resolution.  These contain the maximum resolution
	    of this data and has not been decimated.
	- `:HIGH`: High resolution.  The Douglas-Peucker line reduction was
	    used to reduce data size by ~80% relative to full.
	- `:INTERMEDIATE`: Intermediate resolution.  The Douglas-Peucker line reduction was
	    used to reduce data size by ~80% relative to high.
	- `:LOW`: Low resolution.  The Douglas-Peucker line reduction was
	    used to reduce data size by ~80% relative to intermediate.
	- `:CRUDE`: Crude resolution.  The Douglas-Peucker line reduction was
	    used to reduce data size by ~80% relative to low.

For each resolution there are several levels; these depends on the data type.

The shoreline data are distributed in 6 levels:
    - `1`: Continental land masses and ocean islands, except Antarctica.
    - `2`: Lakes
    - `3`: Islands in lakes
    - `4`: Ponds in islands within lakes
    - `5`: Antarctica based on ice front boundary.
    - `6`: Antarctica based on grounding line boundary.

Note that because GIS software confusingly seem to assume a Cartesian geometry,
any polygon straddling the Dateline is broken into an east and west component.
The most obvious example is Antarctica.
"""
function load_GSHHS_data(;resolution::Symbol=:CRUDE, level::Integer=1)
    # Check for valid inputs
    if !valid_resolution(resolution)
        throw(ArgumentError("Invalid GSHHS product resolution: $(string(resolution))"))
    end

    if !(1 <= level <= 6)
        throw(ArgumentError("Invalid data level $level - Must be between 1 and 6."))
    end

    if (resolution == :CRUDE) && (level == 4)
        throw(ArgumentError("Invalid combination resolution $resolution - $level."))
    end

    path = joinpath(@__DIR__,
                    "..",
                    "data",
                    "GSHHS_shp",
                    RESOLUTION_LOOKUP[resolution],
                    "GSHHS_$(RESOLUTION_LOOKUP[resolution])_L$level")

    return load_data(path)
end

export load_WDBII_border_data
"""
Load WDBI Shapefile data and return all polygons in the data set.

All data sets come in 5 different resolutions:
    - `:FULL`: Full resolution.  These contain the maximum resolution
	    of this data and has not been decimated.
    - `:HIGH`: High resolution.  The Douglas-Peucker line reduction was
	    used to reduce data size by ~80% relative to full.
    - `:INTERMEDIATE`: Intermediate resolution.  The Douglas-Peucker line reduction was
	    used to reduce data size by ~80% relative to high.
    - `:LOW`: Low resolution.  The Douglas-Peucker line reduction was
	    used to reduce data size by ~80% relative to intermediate.
    - `:CRUDE`: Crude resolution.  The Douglas-Peucker line reduction was
	    used to reduce data size by ~80% relative to low.

For each resolution there are several levels; these depends on the data type.

The political boundary data come in 3 levels:
    - `1`: National boundaries.
    - `2`: Internal (state) boundaries for the 8 largest countries only.
    - `3`: Maritime boundaries.
"""
function load_WDBII_border_data(;resolution::Symbol=:CRUDE, level::Integer=1)
    # Check for valid inputs
    if !valid_resolution(resolution)
        throw(ArgumentError("Invalid WDBII product resolution: $(string(resolution))"))
    end

    if !(1 <= level <= 3)
        throw(ArgumentError("Invalid data level $level - Must be between 1 and 3."))
    end

    path = joinpath(@__DIR__,
                    "..",
                    "data",
                    "WDBII_shp",
                    RESOLUTION_LOOKUP[resolution],
                    "WDBII_border_$(RESOLUTION_LOOKUP[resolution])_L$level")

    return load_data(path)
end

export load_WDBII_river_data
"""
Load WDBI Shapefile data and return all polygons in the data set.

All data sets come in 5 different resolutions:
    - `:FULL`: Full resolution.  These contain the maximum resolution
	    of this data and has not been decimated.
    - `:HIGH`: High resolution.  The Douglas-Peucker line reduction was
	    used to reduce data size by ~80% relative to full.
    - `:INTERMEDIATE`: Intermediate resolution.  The Douglas-Peucker line reduction was
	    used to reduce data size by ~80% relative to high.
    - `:LOW`: Low resolution.  The Douglas-Peucker line reduction was
	    used to reduce data size by ~80% relative to intermediate.
    - `:CRUDE`: Crude resolution.  The Douglas-Peucker line reduction was
	    used to reduce data size by ~80% relative to low.

For each resolution there are several levels; these depends on the data type.

The river database come with 11 levels:
    - ` `1: Double-lined rivers (river-lakes).
    - ` `2: Permanent major rivers.
    - ` `3: Additional major rivers.
    - ` `4: Additional rivers.
    - ` `5: Minor rivers.
    - ` `6: Intermittent rivers - major.
    - ` `7: Intermittent rivers - additional.
    - ` `8: Intermittent rivers - minor.
    - ` `9: Major canals.
    - `1`0: Minor canals.
    - `1`1: Irrigation canals.
"""
function load_WDBII_river_data(;resolution::Symbol=:CRUDE, level::Integer=1)
    # Check for valid inputs
    if !valid_resolution(resolution)
        throw(ArgumentError("Invalid WDBII product resolution: $(string(resolution))"))
    end

    if !(1 <= level <= 11)
        throw(ArgumentError("Invalid data level $level - Must be between 1 and 11."))
    end

    path = joinpath(@__DIR__,
                    "..",
                    "data",
                    "WDBII_shp",
                    RESOLUTION_LOOKUP[resolution],
                    "WDBII_river_$(RESOLUTION_LOOKUP[resolution])_L$(@sprintf("%02d", level))")

    return load_data(path)
end

end # module


text:
export ft_setusbparameters

function ft_setusbparameters(ft_handle::Culong,
                             intransfersize::Integer = 4096,
                             outtransfersize::Integer = 4096)
  @assert intransfersize > 63
  @assert outtransfersize > 63
  @assert intransfersize < 2^16+1
  @assert outtransfersize < 2^16+1
  @assert (intransfersize>>4)*16 == intransfersize
  @assert (outtransfersize>>4)*16 == outtransfersize
  ft_status = ccall((:FT_SetUSBParameters, d2xx),
                     Cuint,
                     (Culong, Cuint, Cuint),
                     ft_handle, intransfersize, outtransfersize)
  checkstatus(ft_status)
  return nothing
end