text:
using LinearAlgebra
using Distributed
using SharedArrays

export DSProblem, ProblemSense, SetObjective, SetInitialPoint, SetVariableRange,
       SetVariableRanges, Optimize!, SetIterationLimit, BumpIterationLimit, SetMaxEvals

@enum ProblemSense Min Max
@enum OptimizationStatus Unoptimized PrecisionLimit IterationLimit

"""
	DSProblem{T}(N::Int; poll::AbstractPoll=LTMADS{T}(),
                         search::AbstractSearch=NullSearch(),
                         objective::Union{Function,Nothing}=nothing,
                         initial_point::Vector{T}=zeros(T, N),
                         iteration_limit::Int=1000,
                         sense::ProblemSense=Min)

Return a problem definition for an `N` dimensional problem.

`poll` and `search` specify the poll and search step algorithms to use. The default
choices are (LTMADS)[@ref] and (NullSearch)[@ref] respectively.

Note that if working with `Float64` (normally the case) then the type
parameterisation can be ignored.
"""
mutable struct DSProblem{T} <: AbstractProblem{T}

    #Solver Config
    mesh::AbstractMesh
    poll::AbstractPoll
    search::AbstractSearch

    #Problem Definition
    objective::Function
    constraints::Constraints
    sense::ProblemSense

    #TODO incumbent points should be sets not points, therefore change to vectors of points
    #and remove type unions

    #Feasible Incumbent point
    x::Union{Vector{T},Nothing}
    #Feasible Incumbent point evaluated cost
    x_cost::Union{T,Nothing}

    #Infeasible Incumbent point
    i::Union{Vector{T},Nothing}
    #Infeasible Incumbent point evaluated cost
    i_cost::Union{T,Nothing}

    user_initial_point::Union{Vector{T},Nothing}

    #Problem size
    N::Int

    function_evaluations::Int

    iteration::Int
    iteration_limit::Int
    status::OptimizationStatus

    meshscale::Vector{T}

    cache::AbstractCache

    num_procs::Int
    max_simultanious_evaluations::Int

    #Barrier threshold
    h_max::T

    DSProblem(N::Int;kwargs...) = DSProblem{Float64}(N; kwargs...)

    function DSProblem{T}(N::Int;
                          poll::AbstractPoll=LTMADS{T}(),
                          search::AbstractSearch=NullSearch(),
                          objective::Union{Function,Nothing}=nothing,
                          initial_point::Vector{T}=zeros(T, N),
                          iteration_limit::Int=1000,
                          sense::ProblemSense=Min
                         ) where T

        p = new()

        p.N = N
        p.mesh = Mesh{T}(N)
        p.poll = poll

        p.search = search

        p.iteration = 0
        p.function_evaluations = 0
        p.iteration_limit = iteration_limit
        p.status = Unoptimized
        p.sense = sense

        p.meshscale = ones(p.N)

        p.cache = PointCache{T}()

        p.num_procs = nworkers()
        p.max_simultanious_evaluations = 1

        p.x = nothing
        p.x_cost = nothing
        p.i = nothing
        p.i_cost = nothing

        if objective != nothing
            p.objective = objective
        end

        p.user_initial_point = initial_point

        p.constraints = Constraints{T}()
        return p
    end
end

MeshUpdate!(p::DSProblem, result::IterationOutcome) = MeshUpdate!(p.mesh, p.poll, result)

(GetMeshSize(p::DSProblem{T})::T) where T = p.mesh.Δᵐ

min_mesh_size(::DSProblem{Float64}) = 1.1102230246251565e-16
min_mesh_size(::DSProblem{T}) where T = eps(T)/2

"""
    SetMaxEvals(p::DSProblem, m::Int)

Set the maximum number of simultaneous function evaluations that can be run.
By default this will be set 1.

If (DirectSearch.function_evaluation)[@ref] is not overriden (e.g. for sending
calculation to a cluster) then setting this to a number greater than your PC's
number of threads will result in no improvement.
"""
function SetMaxEvals(p::DSProblem, m::Int)
    if m > p.num_procs
        println("$m is larger than the number of workers that Julia was started with ($(p.num_procs))")
        println("Setting maximum number of workers to $(p.num_procs)")
        println("Start Julia with the option `-p N` where N is the number of additional processes")
        p.max_simultanious_evaluations = p.num_procs
    else
        p.max_simultanious_evaluations = m
    end
end

"""
    SetObjective(p::DSProblem, obj::Function)

Sets the target objective function to solve. `obj` should take a vector and return
a single cost value.
"""
function SetObjective(p::DSProblem, obj::Function)
    if p.sense == Min
        p.objective = obj
    else
        p.objective = x -> -obj(x)
    end
end

"""
    SetIterationLimit(p::DSProblem, i::Int)

Set the maximum number of iterations to `i`.
"""
function SetIterationLimit(p::DSProblem, i::Int)
    if i < p.iteration
        error("Cannot set iteration limit to lower than the number of iterations that have run")
    else
        p.iteration_limit = i
    end
end

"""
    BumpIterationLimit(p::DSProblem, val::Int=100)

Increase the iteration limit by `i`.
"""
function BumpIterationLimit(p::DSProblem; i::Int=100)
    p.iteration_limit += i
end

"""
    SetInitialPoint(p::DSProblem{T}, x::Vector{T}) where T

Set the initial incumbent point to `x`. This must be of the correct dimension. If using
any extreme barrier constraints then it must also satisfy these constraints.
"""
function SetInitialPoint(p::DSProblem{T}, x::Vector{T}) where T
    size(x, 1) == p.N || error("Point dimensions don't match problem definition")
    #TODO Check against constraints
    p.user_initial_point = x
end

function EvaluateInitialPoint(p::DSProblem)
    p.user_initial_point == Nothing() && return
    feasibility = ConstraintEvaluation(p.constraints, p.user_initial_point)
    if feasibility == StrongInfeasible
        error("Initial point must be feasible")
    elseif feasibility == WeakInfeasible
        p.i = p.user_initial_point
        p.i_cost = p.objective(p.user_initial_point)
        CachePush(p, p.i, p.i_cost)
    elseif feasibility == Feasible
        p.x = p.user_initial_point
        p.x_cost = p.objective(p.user_initial_point)
        CachePush(p, p.x, p.x_cost)
    end
    p.user_initial_point = nothing
end

#TODO review if this is the kind of scaling we want to do
"""
	SetVariableRange(p::DSProblem{T}, index::Int, l::T, u::T) where T

Set the expected range of the variable with index `i` to between lower (`l`) and upper (`u`)
values. This **does not** create a constraint, and is only used for scaling when a variable
varies with a significantly different scale to the others.
"""
function SetVariableRange(p::DSProblem{T}, index::Int, l::T, u::T) where T
    1 <= index <= p.N || error("Invalid variable index, should be in range 1 to $(p.N).")
    p.meshscale[index] = 0.1(u-l)
end

"""
	SetVariableRanges(p::DSProblem{T}, l::Vector{T}, u::Vector{T}) where T

Call [`SetVariableRange`](@ref) for each variable. The vectors `l` and `u` should contain
a lower and upper bound for each variable. If it is desired to keep a variable at default
scaling, then give it upper and lower bounds of `-5` and `5` respectively.
"""
function SetVariableRanges(p::DSProblem{T}, l::Vector{T}, u::Vector{T}) where T
    size(l, 1) == p.N || error("Lower bound vector dimensions don't match problem definition")
    size(u, 1) == p.N || error("Upper bound vector dimensions don't match problem definition")
    p.meshscale = @. 0.1(u-l)
end

"""
    Optimize!(p::DSProblem)

Run the direct search algorithm on problem `p`.

`p` must have had its initial point and objective function set. If extreme
barrier constraints have been set then the initial point must be value for
those constraints.
"""
function Optimize!(p::DSProblem)
    #TODO check that problem definition is complete
    EvaluateInitialPoint(p)
    CacheOrderPush(p)

    while p.iteration < p.iteration_limit && GetMeshSize(p) >= min_mesh_size(p)

        result = Search(p)
        #If search fails, run poll
        if result == Unsuccessful
            result = Poll(p)
        end

        result != Unsuccessful && CacheOrderPush(p)

        #pass the result of search/poll to update
        MeshUpdate!(p, result)

        p.iteration += 1
    end

    if p.iteration > p.iteration_limit
        p.status = IterationLimit
    end
    if GetMeshSize(p) <= min_mesh_size(p)
        p.status = PrecisionLimit
    end

    #report_finish(p)
end

"""
    EvaluatePoint!(p::DSProblem{T}, trial_points)::IterationOutcome where T

Determine whether the set of trial points result in a dominating, improving, or unsuccesful
algorithm iteration. Update the feasible and infeasible incumbent points of `p`.
"""
function EvaluatePoint!(p::DSProblem{T}, trial_points)::IterationOutcome where T
    #TODO could split into an evaluation function and an update function
    println("")

    isempty(trial_points) && return Unsuccessful

    #Variables to store the best evaluated points
    feasible_point = isnothing(p.x) ? Inf * ones(p.N) : p.x
    feasible_cost = isnothing(p.x_cost) ? Inf : p.x_cost
    infeasible_point = isnothing(p.i) ? Inf * ones(p.N) : p.i
    infeasible_cost = isnothing(p.i_cost) ? Inf : p.i_cost

    #The current minimum hmax value for all collections
    h_min = GetOldHmaxSum(p.constraints)

    #Iterate over all trial points
    for point in trial_points
        feasibility = ConstraintEvaluation(p.constraints, point)

        # Point violates h_max on at least one constraint collection
        feasibility == StrongInfeasible && continue

        #Point is feasible for relaxed constraints, so evaluate it
        cost = function_evaluation(p, point)

        #To determine if a point is dominating or improving the combined h_max is needed
        h = GetViolationSum(p.constraints, point)

        # Conditions met for a dominant point
        if feasibility == Feasible && cost < feasible_cost
            feasible_point = point
            feasible_cost = cost
        # Conditions met for an improving point (worse cost, but closer to being feasible) or
        # a dominant point (better cost and closer to feasibility)
        # Only record if it offers an improved constraint violation
        elseif feasibility == WeakInfeasible && h < h_min
            infeasible_point = point
            infeasible_cost = cost
            h_min = h
        end
    end

    result = Unsuccessful

    incum_i_cost = isnothing(p.i_cost) ? Inf : p.i_cost
    incum_x_cost = isnothing(p.x_cost) ? Inf : p.x_cost

    # Dominates if there is a feasible improvement, or an infeasible point with
    # reduced violation (determined previously) as well as a cost lower than any
    # yet tested (feasible and infeasible)
    if feasible_cost < incum_x_cost
        result = Dominating
    elseif infeasible_cost < incum_i_cost && h_min < GetOldHmaxSum(p.constraints)
        result = Dominating
    elseif h_min < GetOldHmaxSum(p.constraints)
        result = Improving
    end

    if feasible_cost < incum_x_cost
        p.x = feasible_point
        p.x_cost = feasible_cost
    end

    if infeasible_cost < incum_i_cost || h_min < GetOldHmaxSum(p.constraints)
        p.i = infeasible_point
        p.i_cost = infeasible_cost
    end

    UpdateConstraints(p.constraints, h_min, result, p.x, p.i)

    return result
end

"""
    function_evaluation(p::DSProblem{T}, trial_points::Vector{Vector{T}})::Vector{T} where T

Calculate the cost of the points in `trial_points` and return as a vector.

If the number of available workers is greater than one, and the max_simultanious_evaluations
value of `p` is greater than one then the calculation is distributed across
several cores.

Currently, this has significant overheads and is much slower than evaluating in a single threaded
manner on all testcases. This may give performance benefits when `f` is a heavy, single threaded
operation.

If a specialised way  of calling the function is needed then this function should be overriden, e.g.:

```
function DS.function_evaluation(p::DS.DSProblem{T}, trial_points::Vector{Vector{T}}) where T
	println("I am overriden")
	return map(p.objective, trial_points)
end
```
"""
function function_evaluation(p::DSProblem{T},
                             trial_points::Vector{Vector{T}})::Vector{T} where T
    if p.max_simultanious_evaluations > 1
        costs = SharedArray{T,1}((length(trial_points)))
        #TODO try with threads, might be faster
        @sync @distributed for i in 1:length(trial_points)
            costs[i] = p.objective(trial_points[i])
        end
        return costs
    else
        return map(p.objective, trial_points)
    end
end

"""
	function_evaluation(p::DSProblem{T}, trial_point::Vector{T})::T where T

Evaluate a single trial point with the objective function of `p`.

By default calls the function with the trial point and returns the result. Override to
provide custom evaluation behaviour.
"""
function function_evaluation(p::DSProblem{T},
                             trial_point::Vector{T})::T where T
    CacheQuery(p, trial_point) && return CacheGet(p, trial_point)
    cost = p.objective(trial_point)
    p.function_evaluations += 1
    CachePush(p, trial_point, cost)
	return cost
end


text:
using POMDPs
using POMDPModels
using Test

function test_initial_belief(b0, o0)
    for i=1:5
        if b0[i] != o0
            return false
        end
    end
    return true
end

function test_hist(belief_hist)
    for i=1:length(belief_hist)-1
        if hist.belief_hist[i][2:end] != hist.belief_hist[i+1][1:end-1]
            return false
        end
    end
    return true
end

rng = MersenneTwister(0)
pomdp = RandomPOMDP()

# test constructor
up = KMarkovUpdater(5)

o0 = rand(rng, observations(pomdp))
b0 = initialize_belief(up, fill(o0, up.k))

@test typeof(b0[1]) == typeof(o0)
@test length(b0) == up.k == 5
@test test_initial_belief(b0, o0)

# generate random observation and stack them
o = rand(rng, observations(pomdp))
b = b0
bp = update(up, b, rand(rng, actions(pomdp)), o)
@test bp[end] == o
@test length(bp) == up.k == 5
@test bp[1:end-1] == fill(o0, length(bp)-1)

# check that b is unchanged
@test b == initialize_belief(up, fill(o0, up.k))

b = bp
op = rand(rng, observations(pomdp))
bp = update(up, bp, rand(rng, actions(pomdp)), op)
@test bp[end] == op
@test bp[end-1] == o
@test length(bp) == up.k == 5
@test bp[1:end-2] == fill(o0, length(bp)-2)

# test with history recorder
pomdp = BabyPOMDP()
s0 = initialstate(pomdp, rng)
o0 = generate_o(pomdp, s0, rng)

initial_obs_vec = fill(o0, up.k)
@test_throws ErrorException initialize_belief(up, o0)
@test_throws ErrorException initialize_belief(up, fill(o0, up.k-1))
@test_throws ErrorException update(up, initial_obs_vec, rand(rng, actions(pomdp)), 1.0)

# solver = RandomSolver(rng=rng)
# policy = solve(solver, pomdp)
# hist = simulate(hr, pomdp, policy, up, initial_obs_vec, s0)
# @test hist.belief_hist[1] == fill(o0, up.k)
# @test test_hist(hist.belief_hist)
# @test_throws ErrorException simulate(hr, pomdp, policy, up)

ext: jl


text:
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md

using JuliaFEM, LinearAlgebra, Test

mesh_file = @__DIR__() * "/testdata/primitives.med"
mesh = aster_read_mesh(mesh_file, "TETRA_TET10_1")

problem = Problem(Heat, "tet", 1)
problem_elements = create_elements(mesh, "TET")
update!(problem_elements, "thermal conductivity", 50.0)
add_elements!(problem, problem_elements)

face = Problem(Heat, "face 4", 1)
face_elements = create_elements(mesh, "FACE4")
update!(face_elements, "external temperature", 20.0)
update!(face_elements, "heat transfer coefficient", 60.0)
add_elements!(face, face_elements)

fixed = Problem(Dirichlet, "fixed face 3", 1, "temperature")
fixed_elements = create_elements(mesh, "FACE2")
update!(fixed_elements, "temperature 1", 0.0)
@debug("number of elements in fixed set", length(fixed_elements))
add_elements!(fixed, fixed_elements)

analysis = Analysis(Linear)
add_problems!(analysis, problem, face, fixed)
run!(analysis)
Temp = problem("temperature", 0.0)
Temp_expected = Dict(
    1 => 1.45606533688540E+01,
    2 => 0.0,
    3 => 0.0,
    4 => 0.0,
    5 => 1.05228712963739E+01,
    6 => 0.0,
    7 => 9.44202309239159E+00,
    8 => 1.05228712963739E+01,
    9 => 0.0,
    10 => 0.0)

@debug("Temperature solution using JuliaFEM", Temp)
@debug("Temperature solution using Code Aster", Temp_expected)
Temp_diff = collect(Temp[i]-Temp_expected[i] for i in 1:10)
@debug("Difference between solutions", Temp_diff)
@test isapprox(Temp_diff, zeros(10); atol=1.0e-9)


text:
using Test
using Distributions
using CalibrateEDMF.DistributionUtils
using EnsembleKalmanProcesses.ParameterDistributions

@testset "Priors" begin
    tmpdir = mktempdir()

    params = Dict("foo" => [bounded(0.1, 1.0)], "bar" => [bounded(0.2, 2.0)])
    priors1 = construct_priors(params, outdir_path = tmpdir, to_file = false, unconstrained_σ = 1.0)

    @test isempty(readdir(tmpdir))
    # Construction from Dict does not retain ordering
    @test priors1.name == ["bar", "foo"]

    # Construction from lists (deserializer) is ordered
    prior_dict = Dict(
        "u_names" => ["foo", "bar"],
        "constraints" => [[bounded(0.1, 1.0)], [bounded(0.2, 2.0)]],
        "distributions" => repeat([Parameterized(Normal(0.0, 1.0))], 2),
    )
    priors2 = deserialize_prior(prior_dict)

    @test priors2.name == ["foo", "bar"]
    @test priors2.name != priors1.name
    @test priors2.constraint != priors1.constraint

    μ_correct = Dict("foo" => [0.4], "bar" => [0.5])
    μ_incorrect1 = Dict("foo" => [0.0], "bar" => [2.5]) # Out of bounds
    μ_incorrect2 = Dict("foo" => [0.1]) # Incorrect shape
    priors3 = construct_priors(params, outdir_path = tmpdir, prior_mean = μ_correct)

    @test !isempty(readdir(tmpdir))
    @test readdir(tmpdir)[1] == "prior.jld2"
    @test_throws AssertionError construct_priors(params, outdir_path = tmpdir, prior_mean = μ_incorrect2)
    @test_throws DomainError construct_priors(params, outdir_path = tmpdir, prior_mean = μ_incorrect1)
end

@testset "Prior Vectors" begin
    tmpdir = mktempdir()
    unc_σ = 0.5

    # vector and float -> correct
    params = Dict("foo" => [bounded(0.1, 1.0)], "bar_vect" => [no_constraint(), bounded(-3.0, 3.0), no_constraint()])
    prior_μ = Dict("foo" => [0.4], "bar_vect" => [1.0, 2.0, 3.0])
    priors = construct_priors(params, outdir_path = tmpdir, unconstrained_σ = unc_σ, prior_mean = prior_μ)

    @test priors.name == ["bar_vect_{1}", "bar_vect_{2}", "bar_vect_{3}", "foo"]
    @test priors.constraint == [no_constraint(), bounded(-3.0, 3.0), no_constraint(), bounded(0.1, 1.0)]
    @test priors.distribution[1].distribution.μ == 1.0
    @test priors.distribution[2].distribution.μ ==
          params["bar_vect"][2].constrained_to_unconstrained(prior_μ["bar_vect"][2])
    @test priors.distribution[4].distribution.μ == params["foo"][1].constrained_to_unconstrained(prior_μ["foo"][1])

    # length mismatch [length(bar_vect mu) != length(bar_vect constraint)]
    params = Dict("foo" => [bounded(0.1, 1.0)], "bar_vect" => [repeat([no_constraint()], 3)...])
    prior_μ = Dict("foo" => [0.4], "bar_vect" => ones(2))
    @test_throws AssertionError construct_priors(
        params,
        outdir_path = tmpdir,
        unconstrained_σ = unc_σ,
        prior_mean = prior_μ,
    )

    # unspecified prior mean should yield μ=0
    params = Dict("foo" => [bounded(0.1, 1.0)], "bar_vect" => [repeat([no_constraint()], 3)...])
    priors = construct_priors(params, outdir_path = tmpdir, unconstrained_σ = unc_σ)
    @test [priors.distribution[i].distribution.μ for i in 1:length(priors.distribution)] == zeros(length(priors.distribution))

end

@testset "Transformations" begin
    μ_log = 1.0
    σ_log = 0.5
    μ, σ = mean_and_std_from_ln(μ_log, σ_log)
    μ_log2, σ_log2 = logmean_and_logstd(μ, σ)

    # Test recovery
    @test σ_log ≈ σ_log2
    @test μ_log ≈ μ_log2

end

@testset "flatten_config_dict" begin
    foo_constraint = [bounded(0.1, 1.0)]
    bar_constraint = [bounded(0.2, 2.0)]
    vect_constraint = [repeat([bounded(-1.0, 1.0)], 2)..., bounded(-0.5, 0.5)]

    params = Dict("foo" => foo_constraint, "bar" => bar_constraint, "vect" => vect_constraint)
    flattened_names, flattened_values = DistributionUtils.flatten_config_dict(params)
    @test flattened_names == ["bar", "vect_{1}", "vect_{2}", "vect_{3}", "foo"]
    @test flattened_values[1] == bar_constraint
    @test flattened_values[2][1] == vect_constraint[1]
    @test flattened_values[4][1] == vect_constraint[3]
    @test flattened_values[5] == foo_constraint
end


text:
using PyPlot

vs = randn(100)
plot(vs)


text:
# ---
# title: 24. Swap Nodes in Pairs
# id: problem24
# author: Qling
# date: 2021-05-17
# difficulty: Medium
# categories: Linked List
# link: <https://leetcode.com/problems/swap-nodes-in-pairs/description/>
# hidden: true
# ---
#
# Given a linked list, swap every two adjacent nodes and return its head.
#
# You may **not** modify the values in the list's nodes. Only nodes itself may
# be changed.
#
#
#
# **Example 1:**
#
# ![](https://assets.leetcode.com/uploads/2020/10/03/swap_ex1.jpg)
#
#
#
#     Input: head = [1,2,3,4]
#     Output: [2,1,4,3]
#
#
# **Example 2:**
#
#
#
#     Input: head = []
#     Output: []
#
#
# **Example 3:**
#
#
#
#     Input: head = [1]
#     Output: [1]
#
#
#
#
# **Constraints:**
#
#   * The number of nodes in the list is in the range `[0, 100]`.
#   * `0 <= Node.val <= 100`
#
#
## @lc code=start
using LeetCode

function swap_pairs(head::ListNode)::ListNode
    pre = dummpy = ListNode()
    pre.next = head

    while !isnothing(next(pre)) && !isnothing(next(next(pre)))
        a = pre.next
        b = a.next

        ## `pre->a->b->b.next` => `pre->b->a->b.next`
        pre.next, b.next, a.next = b, a, b.next

        pre = a
    end

    return dummpy.next
end

## @lc code=end


text:
using Test
using WavePropBase
using WavePropBase.Geometry
using WavePropBase.IO
using WavePropBase.Interpolation
using WavePropBase.Mesh
using StaticArrays

@testset "Sphere geo" begin
    # Test the simple sphere geometry
    fname = joinpath(@__DIR__,"sphere.geo")
    clear_entities!()
    Ω,M = read_geo(fname)
    @test length(Ω) == 1
    @test geometric_dimension(Ω) == 3

    Γi = internal_boundary(Ω)
    Γe = external_boundary(Ω)

    @test length(Γi) == 0
    @test length(Γe) == 1
    @test geometric_dimension(Γe) == 2

    Ce = external_boundary(Γe)
    Ci = internal_boundary(Γe)
    S  = skeleton(Γe)
    @test length(Ce) == 2
    @test geometric_dimension(Ce) == 1
end

@testset "Sphere msh" begin
    # Test the simple sphere geometry
    clear_entities!()
    fname = joinpath(@__DIR__,"sphere.msh")
    Ω, M  = read_msh(fname)
    # Test internal creation of sphere
    Ω, M = WavePropBase.IO.gmsh_sphere(dim=2)
end

@testset "Disk" begin
    # Test internal creation of disk
    clear_entities!()
    Ω, M = WavePropBase.IO.gmsh_disk()
end

@testset "Element iterator" begin
    clear_entities!()
    (lx,ly,lz) = widths = (1.,1.,2.)
    Ω, M  = WavePropBase.IO.gmsh_box(;widths=widths)
    E    = keys(M) |> first
    iter = ElementIterator{E}(M)
    @test eltype(iter) == E
    @test length(iter) == size(M.elements[E],2)
end

@testset "Sub mesh" begin
    clear_entities!()
    (lx,ly,lz) = widths = (1.,1.,2.)
    Ω, M  = WavePropBase.IO.gmsh_box(;widths=widths)
    subM  = SubMesh(M,external_boundary(Ω))
    E    = keys(subM) |> first
    iter = ElementIterator(subM,E)
    @test eltype(iter) == E
    @test length(iter) == size(M.elements[E],2)
end

max_stars_repo_name: WaveProp/WavePropBase

max_issues_repo_name: WaveProp/WavePropBase

max_forks_repo_name: WaveProp/WavePropBase


text:
"""

    @primitive  fx g1 g2...

Define a new primitive operation for AutoGrad and (optionally) specify its gradients.
Non-differentiable functions such as `sign`, and non-numeric functions such as `size` should
be defined using the @zerograd macro instead.

# Examples

    @primitive sin(x::Number)
    @primitive hypot(x1,x2),dy,y

    @primitive sin(x::Number),dy  (dy.*cos(x))
    @primitive hypot(x1,x2),dy,y  (dy.*x1./y)  (dy.*x2./y)

The first example shows that `fx` is a typed method declaration.  Julia supports multiple
dispatch, i.e. a single function can have multiple methods with different arg types.
AutoGrad takes advantage of this and supports multiple dispatch for primitives and
gradients.

The second example specifies variable names for the output gradient `dy` and the output `y`
after the method declaration which can be used in gradient expressions.  Untyped, ellipsis
and keyword arguments are ok as in `f(a::Int,b,c...;d=1)`.  Parametric methods such as
`f(x::T) where {T<:Number}` cannot be used.

The method declaration can optionally be followed by gradient expressions.  The third and
fourth examples show how gradients can be specified.  Note that the parameters, the return
variable and the output gradient of the original function can be used in the gradient
expressions.

# Under the hood

The @primitive macro turns the first example into:

    sin(x::Value{T}) where {T<:Number} = forw(sin, x)

This will cause calls to `sin` with a boxed argument (`Value{T<:Number}`) to be recorded.
The recorded operations are used by AutoGrad to construct a dynamic computational graph.
With multiple arguments things are a bit more complicated.  Here is what happens with the
second example:

    hypot(x1::Value{S}, x2::Value{T}) where {S,T} = forw(hypot, x1, x2)
    hypot(x1::S, x2::Value{T})        where {S,T} = forw(hypot, x1, x2)
    hypot(x1::Value{S}, x2::T)        where {S,T} = forw(hypot, x1, x2)

We want the forw method to be called if any one of the arguments is a boxed `Value`.  There
is no easy way to specify this in Julia, so the macro generates all 2^N-1 boxed/unboxed
argument combinations.

In AutoGrad, gradients are defined using gradient methods that have the following pattern:

    back(f,Arg{i},dy,y,x...) => dx[i]

For the third example here is the generated gradient method:

    back(::typeof(sin), ::Type{Arg{1}}, dy, y, x::Value{T}) where {T<:Number} = dy .* cos(x)

For the last example a different gradient method is generated for each argument:

    back(::typeof(hypot), ::Type{Arg{1}}, dy, y, x1::Value{S}, x2::Value{T}) where {S,T} = (dy .* x1) ./ y
    back(::typeof(hypot), ::Type{Arg{2}}, dy, y, x1::Value{S}, x2::Value{T}) where {S,T} = (dy .* x2) ./ y

In fact @primitive generates four more definitions for the other boxed/unboxed argument
combinations.

# Broadcasting

Broadcasting is handled by extra `forw` and `back` methods. `@primitive` defines the following
so that broadcasting of a primitive function with a boxed value triggers `forw` and `back`.

    broadcasted(::typeof(sin), x::Value{T}) where {T<:Number} = forw(broadcasted,sin,x)
    back(::typeof(broadcasted), ::Type{Arg{2}}, dy, y, ::typeof(sin), x::Value{T}) where {T<:Number} = dy .* cos(x)

If you do not want the broadcasting methods, you can use the `@primitive1` macro. If you
only want the broadcasting methods use `@primitive2`. As a motivating example, here is how
`*` is defined for non-scalars:

    @primitive1 *(x1,x2),dy  (dy*x2')  (x1'*dy)
    @primitive2 *(x1,x2),dy  unbroadcast(x1,dy.*x2)  unbroadcast(x2,x1.*dy)

Regular `*` is matrix multiplication, broadcasted `*` is elementwise multiplication and the
two have different gradients as defined above. `unbroadcast(a,b)` reduces `b` to the same
shape as `a` by performing the necessary summations.
"""
:(@primitive), :(@primitive1), :(@primitive2)

macro primitive(f,g...)                         # @primitive sin(x::Number),dy,y  (dy.*cos.(x))
    (f,dy,y) = fparse(f)
    b = Expr(:block)
    forwcall = fcall(f)                	    	# forw(sin,x)
    forwcast = fcall(f,broadcasted=true)        # forw(broadcasted,sin,x)
    for fx in fsigs(f)                          # sin(x::Value{T}) where {T<:Number}
        push!(b.args, :($fx = $forwcall))
        bfx = f2b(fx)                           # broadcasted(::typeof(sin), x::Value{T}) where {T<:Number}
        push!(b.args, :($bfx = $forwcast))
        for i=1:length(g)
            gx = gsig(fx,dy,y,i)                # back(::typeof(sin), ::Type{Arg{1}}, dy, y, x::Value{T}) where {T<:Number}
            push!(b.args, :($gx = $(g[i])))     # '' = (dy.*cos.(x))
            bgx = bsig(fx,dy,y,i)               # back(::typeof(broadcasted), ::Type{Arg{2}}, dy, y, ::typeof(sin), x::Value{T}) where {T<:Number}
            push!(b.args, :($bgx = $(g[i])))
        end
    end
    return esc(b)
end

macro primitive1(f,g...)        # non-broadcasting version
    (f,dy,y) = fparse(f)
    b = Expr(:block)
    forwcall = fcall(f)
    for fx in fsigs(f)
        push!(b.args, :($fx = $forwcall))
        for i=1:length(g)
            gx = gsig(fx,dy,y,i)
            push!(b.args, :($gx = $(g[i])))
        end
    end
    return esc(b)
end

macro primitive2(f,g...)        # broadcasting-only version
    (f,dy,y) = fparse(f)
    b = Expr(:block)
    forwcast = fcall(f,broadcasted=true)
    for fx in fsigs(f)
        bfx = f2b(fx)
        push!(b.args, :($bfx = $forwcast))
        for i=1:length(g)
            bgx = bsig(fx,dy,y,i)
            push!(b.args, :($bgx = $(g[i])))
        end
    end
    return esc(b)
end

"""

    @zerograd f(args...; kwargs...)

Define `f` as an AutoGrad primitive operation with zero gradient.

# Example:

    @zerograd  floor(x::Float32)

`@zerograd` allows `f` to handle boxed `Value` inputs by unboxing them like a `@primitive`,
but unlike `@primitive` it does not record its actions or return a boxed `Value` result.
Some functions, like `sign()`, have zero gradient.  Others, like `length()` have discrete or
constant outputs.  These need to handle `Value` inputs, but do not need to record anything
and can return regular values.  Their output can be treated like a constant in the program.
Use the `@zerograd` macro for those.  Use the `@zerograd1` variant if you don't want to
define the broadcasting version and `@zerograd2` if you only want to define the broadcasting
version. Note that `kwargs` are NOT unboxed.

"""
:(@zerograd), :(@zerograd1), :(@zerograd2)

macro zerograd(f)                               # @zerograd sign(x::Number)
    (f,dy,y) = fparse(f)
    b = Expr(:block)
    for fx in fsigs(f)                          # sign(x::Value{T}) where {T<:Number}
        zx = zcall(fx)                          # sign(value(x))
        push!(b.args, esc(:($fx = $zx)))
        (bfx,bzx) = bzcall(fx,zx)
        push!(b.args, esc(:($bfx = $bzx)))      # broadcasted(::typeof(sign), x::Value{T}) where {T <: Number} = broadcasted(sign, value(x))
        bx = v2b(fx)                            # sign(x::Bcasted{T}) where {T<:Number}
        push!(b.args, esc(:($bx = AutoGrad.Bcasted($bzx)))) # ... = Bcasted(broadcasted(sign, value(x)))
    end
    return b
end

macro zerograd1(f)   # non-broadcasting version
    (f,dy,y) = fparse(f)
    b = Expr(:block)
    for fx in fsigs(f)
        zx = zcall(fx)
        push!(b.args, esc(:($fx = $zx)))
    end
    return b
end

macro zerograd2(f)   # broadcasting-only version
    (f,dy,y) = fparse(f)
    b = Expr(:block)
    for fx in fsigs(f)
        zx = zcall(fx)
        #push!(b.args, esc(:($fx = $zx)))
        (bfx,bzx) = bzcall(fx,zx)
        push!(b.args, esc(:($bfx = $bzx)))
        bx = v2b(fx)
        push!(b.args, esc(:($bx = AutoGrad.Bcasted($bzx))))
    end
    return b
end

function fparse(f)
    isa(f,Expr) || error("'$f' not a method signature")
    if f.head == :tuple # Using f(x),dy,y to indicate return variable for gradients
        if length(f.args) == 3
            (f,dy,y) = f.args
        elseif length(f.args) == 2
            (f,dy) = f.args; y = gensym()
        else
            error("The first arg '$f' should have the format f(x),dy,y")
        end
    else
        dy = gensym(); y = gensym()
    end
    f.head == :call || error("'$f' not a method signature")
    isa(dy,Symbol) || error("Output gradient '$dy' not a symbol")
    isa(y,Symbol) || error("Return variable '$y' not a symbol")
    return (f,dy,y)
end

function fcall(f; broadcasted = false)
    rx = notypes(f)
    fn = rx.args[1]
    rx.args[1]=:(AutoGrad.forw)
    # Need to fix kwargs
    r2 = rx.args[2]
    if isa(r2,Expr) && r2.head == :parameters
        for i in 1:length(r2.args)
            k = r2.args[i]
            if isa(k,Symbol); r2.args[i] = Expr(:kw,k,k)
            elseif !isa(k,Expr); error("Bad kwarg '$k'")
            elseif k.head == :(...); continue
            elseif k.head != :kw; error("Bad kwarg '$k'")
            elseif !isa(k.args[1],Symbol); error("Bad kwarg '$k'")
            else; k.args[2]=k.args[1]; end
        end
        a = 3
    else
        a = 2
    end
    insert!(rx.args,a,fn)
    if broadcasted; insert!(rx.args,a,:(Base.Broadcast.broadcasted)); end
    return rx
end

# Input is of the form: (where (call f (:: x (curly (. AutoGrad Value) T))) (<: T Int))
function zcall(f)
    z = copy(f.args[1])
    z1 = z.args[1]
    isa(z1,Expr) && z1.head==:curly && (z.args[1]=z1.args[1]) # This should not be needed in 0.7
    for i=2:length(z.args)
        zi = z.args[i]
        if isa(zi,Symbol)
            # all done
        elseif !isa(zi,Expr)
            error("Unrecognized argtype '$zi'")
        elseif zi.head==:(::)
            (v,t) = zi.args
            if t==:(AutoGrad.Value) || (isa(t,Expr) && t.head==:curly && t.args[1]==:(AutoGrad.Value))
                z.args[i] = :(value($v))
            else
                z.args[i] = v
            end
        elseif zi.head==:(...)  # done
        elseif zi.head==:parameters # done
        else
            error("Unrecognized argtype '$zi'")
        end
    end
    return notypes(z)
end

function bzcall(fx,zx)
    bfx = copy(fx)
    g = bfx.args[1]
    fname = g.args[1]
    g.args[1] = :(Base.Broadcast.broadcasted)
    if g.args[2].head == :parameters; a = 3; else; a = 2; end
    insert!(g.args, a, :(::typeof($fname)))
    bzx = copy(zx)
    bzx.args[1] = :(Base.Broadcast.broadcasted)
    insert!(bzx.args, a, fname)
    return (bfx,bzx)
end

# eliminate type declarations from a function call
function notypes(ex)
    if isa(ex, Expr)
        if (ex.head == :(::) || ex.head == :curly)
            return notypes(ex.args[1])
        else
            return Expr(ex.head, map(notypes, ex.args)...)
        end
    else
        return ex
    end
end

# input: (where (call f (:: x (curly Value T))) (<: T Int))
# output: (where (call broadcasted (:: (call typeof f)) (:: x (curly Value T))) (<: T Int))
function f2b(fx)
    bx = copy(fx)               # where...
    cx = bx.args[1]             # call...
    f = cx.args[1]              # func
    cx.args[1] = :(Base.Broadcast.broadcasted)
    if cx.args[2].head == :parameters; a = 3; else; a = 2; end
    insert!(cx.args, a, :(::typeof($f)))
    return bx
end

# change AutoGrad.Value -> AutoGrad.Bcasted
function v2b(fx)
    if fx == :(AutoGrad.Value)
        :(AutoGrad.Bcasted)
    elseif isa(fx, Expr)
        Expr(fx.head, v2b.(fx.args)...)
    else
        fx
    end
end

# create type signatures for f where one or more args are Value's.
# With multiple args add Value to each subset combinatorially.
# The input has the form (call f (:: x Int))
# The 0.6 output was     (call (curly f (<: T Int)) (:: x (curly Value T)))
# The 0.7 output is      (where (call f (:: x (curly Value T))) (<: T Int))
function fsigs(f)               # sin(x::Number) => sin(x::Value{T}) where {T <: Number}
    f1 = copy(f)                # sin(x::Number)
    fname = f1.args[1]
    if isa(fname, Symbol) && isdefined(@__MODULE__, fname) # @__MODULE__ here resolves to AutoGrad when compiling Knet
        m = which(@__MODULE__, fname)                      # which does not work for symbols undefined in AutoGrad
        f1.args[1] = :(($m).$fname) # Base.sin
    end
    a1 = Expr(:where,f1)        # sin(x::Number) where {}
    nargs = 0
    for i=2:length(f1.args)
        ai = f1.args[i]         # x::Number
        if isa(ai,Symbol)
            nargs+=1
            ti = gensym()
            push!(a1.args, Expr(:<:, ti, Any))
            f1.args[i] = Expr(:(::),ai,ti)
        elseif !isa(ai,Expr)
            error("Neither Symbol nor Expr: $ai")
        elseif in(ai.head, (:parameters, :(...)))
            continue
        elseif ai.head == :(::)
            nargs+=1
            ti = gensym()
            push!(a1.args, Expr(:<:,ti,ai.args[2]))
            ai.args[2] = ti
        else
            error("Argtype not supported: '$ai'")
        end
    end
    flist = []
    for nodes=0:(1<<nargs-2)
        fn = copy(a1)           # sin(x::T) where {T <: Number}
        f1 = fn.args[1]         # sin(x::T)
        iargs = 0
        for i=2:length(f1.args)
            ai = f1.args[i]     # (x::T)
            in(ai.head, (:parameters, :(...))) && continue
            ai.head == :(::) || error("Bad arg '$ai'")
            if nodes & (1<<iargs) == 0
                ai.args[2] = :(AutoGrad.Value{$(ai.args[2])}) #Expr(:curly,:Value,ai.args[2])
            end
            iargs += 1
        end
        push!(flist, fn)        # sin(x::Value{T}) where {T <: Number}
    end
    return flist
end

# The first input to gsig is an output of fsigs, e.g.
# (where (call f (:: x (curly Value T))) (<: T Int))
function gsig(f,dy,y,i)
    fcopy = copy(f)
    g = fcopy.args[1]
    fname = g.args[1]
    g.args[1] = :(AutoGrad.back)
    if g.args[2].head == :parameters; a = 3; else; a = 2; end
    insert!(g.args, a, :(::typeof($fname)))
    insert!(g.args, a+1, :(::Type{AutoGrad.Arg{$i}}))
    insert!(g.args, a+2, dy)
    insert!(g.args, a+3, y)
    return fcopy
end

# This is for the broadcast version
# Input: (where (call f (:: x (curly Value T))) (<: T Int))
# Output: (where (call broadcasted :(::Type{Grad{2}}) dy y :(::typeof(f)) :(x::Value{T})) (<: T Int))
function bsig(f,dy,y,i)
    fcopy = copy(f)
    g = fcopy.args[1]
    fname = g.args[1]
    g.args[1] = :(AutoGrad.back)
    if g.args[2].head == :parameters; a = 3; else; a = 2; end
    insert!(g.args, a, :(::typeof(Base.Broadcast.broadcasted)))
    insert!(g.args, a+1, :(::Type{AutoGrad.Arg{$(i+1)}}))
    insert!(g.args, a+2, dy)
    insert!(g.args, a+3, y)
    insert!(g.args, a+4, :(::typeof($fname)))
    return fcopy
end


text:
using JGUI
using Base.Test

## connect values
w = window()
f = vbox(w)
push!(f)

ed = lineedit(f, "")
lab = label(f, "")
append!(f, [ed, lab])

lab[:value] = ed
ed[:value] = "new value"
@assert lab[:value] == ed[:value]

raise(w)

if !JGUI.isqt()

    ## cairographics
    ## manipulate like
    w = window()
    f=vbox(w)
    push!(f)

    n = slider(f, 1:10)
    g = cairographic(f)

    append!(f, [n, g])

    @wlift display(g, plot(sin, 0, n*pi))

end


text:
module Variational

using ..Core, ..Core.RandomVariables, ..Utilities
using Distributions, Bijectors
using ProgressMeter, LinearAlgebra
using ..Turing: PROGRESS
using ..Turing: Model, SampleFromPrior, SampleFromUniform
using ..Turing: Turing
using Random: AbstractRNG

using ForwardDiff
using Tracker

import ..Core: getchunksize, getADtype

using Requires
function __init__()
    @require Flux="587475ba-b771-5e3f-ad9e-33799f191a9c" apply!(o, x, Δ) = Flux.Optimise.apply!(o, x, Δ)
end

export
    vi,
    ADVI,
    ELBO,
    TruncatedADAGrad

abstract type VariationalInference{AD} end

getchunksize(::T) where {T <: VariationalInference} = getchunksize(T)
getchunksize(::Type{<:VariationalInference{AD}}) where AD = getchunksize(AD)
getADtype(alg::VariationalInference) = getADtype(typeof(alg))
getADtype(::Type{<: VariationalInference{AD}}) where {AD} = AD

abstract type VariationalObjective end

const VariationalPosterior = Distribution{Multivariate, Continuous}

"""
    rand(vi::VariationalInference, num_samples)

Produces `num_samples` samples for the given VI method using number of samples equal
to `num_samples`.
"""
function rand(vi::VariationalPosterior, num_samples) end

"""
    objective(vi::VariationalInference, q::VariationalPosterior, model::Model, args...)

Computes the variational objective to be optimized for a given VI method.
"""
function objective(
    vi::VariationalInference,
    q::VariationalPosterior,
    model::Model,
    num_samples)
end

"""
    grad!(vo, vi::VariationalInference, q::VariationalPosterior, model::Model, θ, out, args...)

Computes the gradients used in `optimize!`. Default implementation is provided for
`VariationalInference{AD}` where `AD` is either `ForwardDiffAD` or `TrackerAD`.
This implicitly also gives a default implementation of `optimize!`.

Variance reduction techniques, e.g. control variates, should be implemented in this function.
"""
function grad!(
    vo, vi::VariationalInference{AD},
    q::VariationalPosterior,
    model::Model,
    θ,
    out,
    args...
) where AD
    error("Turing.Variational.grad!: unmanaged variational inference algorithm: "
          * "$(typeof(alg))")
end

"""
    vi(model::Model, alg::VariationalInference)
    vi(model::Model, alg::VariationalInference, q::VariationalPosterior)

Constructs the variational posterior from the `model` and performs the optimization
following the configuration of the given `VariationalInference` instance.
"""
function vi(model::Model, alg::VariationalInference)
    error("Turing.Variational.vi: variational inference algorithm $(typeof(alg)) "
          * "is not implemented")
end
function vi(model::Model, alg::VariationalInference, q::VariationalPosterior)
    error("Turing.Variational.vi: variational inference algorithm $(typeof(alg)) "
          * "is not implemented")
end

# default implementations
function grad!(
    vo,
    alg::VariationalInference{AD},
    q::VariationalPosterior,
    model::Model,
    θ::AbstractVector{T},
    out::DiffResults.MutableDiffResult,
    args...
) where {T <: Real, AD <: ForwardDiffAD}
    # TODO: this probably slows down executation quite a bit; exists a better way
    # of doing this?
    f(θ_) = - vo(alg, q, model, θ_, args...)

    chunk_size = getchunksize(alg)
    # Set chunk size and do ForwardMode.
    chunk = ForwardDiff.Chunk(min(length(θ), chunk_size))
    config = ForwardDiff.GradientConfig(f, θ, chunk)
    ForwardDiff.gradient!(out, f, θ, config)
end

function grad!(
    vo,
    alg::VariationalInference{AD},
    q::VariationalPosterior,
    model::Model,
    θ::AbstractVector{T},
    out::DiffResults.MutableDiffResult,
    args...
) where {T <: Real, AD <: TrackerAD}
    θ_tracked = [Tracker.param(θ[i]) for i ∈ eachindex(θ)]
    y = - vo(alg, q, model, θ_tracked, args...)
    Tracker.back!(y, 1.0)

    DiffResults.value!(out, Tracker.data(y))
    DiffResults.gradient!(out, [Tracker.grad(θ_tracked[i]) for i ∈ eachindex(θ_tracked)])
end

"""
    optimize!(vo, alg::VariationalInference{AD}, q::VariationalPosterior, model::Model, θ; optimizer = TruncatedADAGrad())

Iteratively updates parameters by calling `grad!` and using the given `optimizer` to compute
the steps.
"""
function optimize!(
    vo,
    alg::VariationalInference{AD},
    q::VariationalPosterior,
    model::Model,
    θ;
    optimizer = TruncatedADAGrad()
) where AD
    # TODO: should we always assume `samples_per_step` and `max_iters` for all algos?
    alg_name = alg_str(alg)
    samples_per_step = alg.samples_per_step
    max_iters = alg.max_iters

    # number of previous gradients to use to compute `s` in adaGrad
    stepsize_num_prev = 10

    num_params = length(q)

    # TODO: really need a better way to warn the user about potentially
    # not using the correct accumulator
    if (optimizer isa TruncatedADAGrad) && (θ ∉ keys(optimizer.acc))
        # this message should only occurr once in the optimization process
        @info "[$alg_name] Should only be seen once: optimizer created for θ" objectid(θ)
    end

    diff_result = DiffResults.GradientResult(θ)

    i = 0
    prog = if PROGRESS[]
        ProgressMeter.Progress(max_iters, 1, "[$alg_name] Optimizing...", 0)
    else
        0
    end

    # add criterion? A running mean maybe?
    time_elapsed = @elapsed while (i < max_iters) # & converged
        grad!(vo, alg, q, model, θ, diff_result, samples_per_step)

        # apply update rule
        Δ = DiffResults.gradient(diff_result)
        Δ = apply!(optimizer, θ, Δ)
        @. θ = θ - Δ

        Turing.DEBUG && @debug "Step $i" Δ DiffResults.value(diff_result)
        PROGRESS[] && (ProgressMeter.next!(prog))

        i += 1
    end

    return θ
end

# distributions
include("distributions.jl")

# objectives
include("objectives.jl")

# optimisers
include("optimisers.jl")

# VI algorithms
include("advi.jl")

end


text:
using ITensorNetworkAD
using Test

@testset "ITensorNetworkAD.jl" begin
  for filename in [
    "ITensorAutoHOOT/runtests.jl",
    "ITensorChainRules/runtests.jl",
    "ITensorNetworks/runtests.jl",
    "Optimizations/runtests.jl",
  ]
    println("Running $filename")
    include(filename)
  end
end


text:
sprand_reward{R}(::Type{R}, states, actions) =
    sprand(states, actions, 1/3, x -> rand_reward(R, x))

if VERSION < v"0.4-"
    # Adapted from
    # https://github.com/JuliaLang/julia/blob/v0.3.7/base/sparse/sparsematrix.jl
    function sprand_IJ(m::Integer, n::Integer, density::FloatingPoint)
        ((m < 0) || (n < 0)) && throw(ArgumentError("invalid Array dimensions"))
        0 <= density <= 1 || throw(ArgumentError("$density not in [0,1]"))
        N = n*m

        I, J = Array(Int, 0), Array(Int, 0) # indices of nonzero elements
        sizehint(I, int(N*density))
        sizehint(J, int(N*density))

        # density of nonzero columns:
        L = log1p(-density)
        coldensity = -expm1(m*L) # = 1 - (1-density)^m
        colsparsity = exp(m*L) # = 1 - coldensity
        L = 1/L

        rows = Array(Int, 0)
        for j in randsubseq(1:n, coldensity)
            # To get the right statistics, we *must* have a nonempty column j
            # even if p*m << 1.   To do this, we use an approach similar to
            # the one in randsubseq to compute the expected first nonzero row k,
            # except given that at least one is nonzero (via Bayes' rule);
            # carefully rearranged to avoid excessive roundoff errors.
            k = ceil(log(colsparsity + rand()*coldensity) * L)
            ik = k < 1 ? 1 : k > m ? m : int(k) # roundoff-error/underflow paranoia
            randsubseq!(rows, 1:m-ik, density)
            push!(rows, m-ik+1)
            append!(I, rows)
            nrows = length(rows)
            Jlen = length(J)
            resize!(J, Jlen+nrows)
            @inbounds for i = Jlen+1:length(J)
                J[i] = j
            end
        end
        I, J
    end
else
    import Base.SparseMatrix.sprand_IJ
end

function sprand_transition{T}(::Type{T}, states)
    I, J = VERSION < v"0.4-" ?
        sprand_IJ(states, states, 1/3) :
        sprand_IJ(Base.Random.GLOBAL_RNG, states, states, 1/3)
    V = rand(T, length(I))
    for j = 1:states
        idx = find(J .== j)
        if length(idx) == 0
            i = rand(1:states)
            push!(I, i)
            push!(J, j)
            push!(V, one(T))
        else
            V[idx] = V[idx] ./ sum(V[idx])
        end
    end
    sparse(I, J, V, states, states)
end

function sprandom{P<:FloatingPoint,R<:FloatingPoint}(
    ::Type{P},
    ::Type{R},
    states,
    actions,
)
    transition = SparseMatrixCSC{P,Int}[
        sprand_transition(P, states) for a = 1:actions
    ]
    reward = sprand_reward(R, states, actions)
    (transition, reward)
end

sprandom(states, actions) = sprandom(Float64, Float64, states, actions)


text:
const menubarChoices = Menubar.AbstractMenuItem[
    Menubar.MenuChoices(
        "File",
        Menubar.AbstractMenuItem[
            Menubar.MenuChoice("New", createNewMap, Image(icon_name="document-new-symbolic", size=:MENU)),
            Menubar.MenuChoice("Open", showFileOpenDialog, Image(icon_name="document-open-symbolic")),
            Menubar.MenuSeparator(),
            Menubar.MenuChoice("Save", menuFileSave, Image(icon_name="document-save-symbolic")),
            Menubar.MenuChoice("Save as", showFileSaveDialog, Image(icon_name="document-save-as-symbolic")),
            Menubar.MenuSeparator(),
            Menubar.MenuChoice("Exit", ExitWindow.exitAhorn, Image(icon_name="application-exit-symbolic")),
        ]
    ),
    Menubar.MenuChoices(
        "Edit",
        Menubar.AbstractMenuItem[
            Menubar.MenuChoice("Undo", History.undo!, Image(icon_name="edit-undo-symbolic")),
            Menubar.MenuChoice("Redo", History.redo!, Image(icon_name="edit-redo-symbolic")),
            Menubar.MenuSeparator(),
            Menubar.MenuChoice("Settings", SettingsWindow.showSettings, Image(icon_name="preferences-system-symbolic")),
        ]
    ),
    Menubar.MenuChoices(
        "View",
        Menubar.AbstractMenuItem[
            Menubar.MenuChoices(
                "RoomList",
                Menubar.AbstractMenuItem[
                    Menubar.MenuRadioGroup(
                        [
                            ("All", showAllRoomListColumns),
                            ("Room Name", showRoomListNameColumn),
                            ("None", hideAllRoomListColumns),
                        ],
                        getRoomListVisibilityRadioIndex
                    )
                ]
            ),
            Menubar.MenuSeparator(),
            Menubar.MenuCheck("Foreground Tiles", (w) -> setGlobalLayerVisiblity("fgTiles", w), true),
            Menubar.MenuCheck("Background Tiles", (w) -> setGlobalLayerVisiblity("bgTiles", w), true),
            Menubar.MenuCheck("Entities", (w) -> setGlobalLayerVisiblity("entities", w), true),
            Menubar.MenuCheck("Triggers", (w) -> setGlobalLayerVisiblity("triggers", w), true),
            Menubar.MenuCheck("Foreground Decals", (w) -> setGlobalLayerVisiblity("fgDecals", w), true),
            Menubar.MenuCheck("Background Decals", (w) -> setGlobalLayerVisiblity("bgDecals", w), true),
        ]
    ),
    Menubar.MenuChoices(
        "Map",
        Menubar.AbstractMenuItem[
            Menubar.MenuChoice("Stylegrounds", StylegroundWindow.createStylegroundWindow, Image(icon_name="document-page-setup-symbolic")),
            Menubar.MenuChoice("Metadata", MetadataWindow.configureMetadata, Image(icon_name="document-properties-symbolic")),
            Menubar.MenuSeparator(),
            Menubar.MenuChoice("Save Map Image", MapImageDumper.dumpMapImageDialog, Image(icon_name="image-x-generic-symbolic"))
        ]
    ),
    Menubar.MenuChoices(
        "Room",
        Menubar.AbstractMenuItem[
            Menubar.MenuChoice("Add", RoomWindow.createRoom, Image(icon_name="star-new-symbolic")),
            Menubar.MenuChoice("Configure", RoomWindow.configureRoom, Image(icon_name="system-run-symbolic")),
            Menubar.MenuSeparator(),
            Menubar.MenuChoice("Delete", deleteRoomCallback, Image(icon_name="edit-delete-symbolic")),
        ]
    ),
]

const menubarDebugChoices = Menubar.AbstractMenuItem[
    Menubar.MenuChoices(
        "Debug",
        Menubar.AbstractMenuItem[
            Menubar.MenuChoices(
                "Reload",
                Menubar.AbstractMenuItem[
                    Menubar.MenuChoice("Tools", (w) -> debug.reloadTools!()),
                    Menubar.MenuChoice("Entities", (w) -> debug.reloadEntities!()),
                    Menubar.MenuChoice("Triggers", (w) -> debug.reloadTriggers!()),
                    Menubar.MenuChoice("Effects", (w) -> debug.reloadEffects!()),
                    Menubar.MenuChoice("Libraries", (w) -> debug.reloadLibraries!()),
                    Menubar.MenuChoice("External Sprites", (w) -> loadAllExternalSprites!(force=true)),
                    Menubar.MenuChoice("Language Files", (w) -> debug.reloadLangdata()),
                ]
            ),
            Menubar.MenuChoice("Clear Map Render Cache", (w) -> debug.clearMapDrawingCache!()),
            Menubar.MenuChoice("Force Draw All Rooms", (w) -> debug.forceDrawWholeMap!()),
            Menubar.MenuChoice("Check Tooltip Coverage", (w) -> debug.checkTooltipCoverage()),
            Menubar.MenuChoice("Dump Sprites", SpriteDumperWindow.dumpSpritesDialog),
            Menubar.MenuChoice("Test Console", TestEvalWindow.showEvalWindow),
        ]
    )
]

const menubarEndChoices = Menubar.AbstractMenuItem[
    Menubar.MenuChoices(
        "Help",
        Menubar.AbstractMenuItem[
            Menubar.MenuChoice("Check for Updates", UpdateWindow.updateAhorn, Image(icon_name="software-update-available-symbolic")),
            Menubar.MenuChoice("About", AboutWindow.showAboutWindow, Image(icon_name="help-about-symbolic")),
        ]
    ),
]

menubar = Gtk.MenuBar()

max_stars_repo_name: CelestialCartographers/Ahorn

max_issues_repo_name: CelestialCartographers/Ahorn

max_forks_repo_name: CelestialCartographers/Ahorn


text:
var2(x) = mean(x.^2) .- mean(x)^2

function waic( ll::AbstractArray; pointwise=false , log_lik="log_lik" , kwargs... )

    n_samples, n_obs = size(ll)
    pD = zeros(n_obs);

    lpd = reshape(logsumexp(ll .- log(n_samples); dims=1), n_obs);
    for i in 1:n_obs
        pD[i] = var2(ll[:,i])
    end

    waic_vec = (-2) .* ( lpd - pD );
    if pointwise == false
        waics = sum(waic_vec)
        lpd = sum(lpd)
        pD = sum(pD)
    else
        waics = waic_vec
    end

    local se
    try
        se = sqrt( n_obs*var2(waic_vec) )
    catch e
        println(e)
        se = nothing
    end

    (WAIC=waics, lppd=lpd, penalty=pD, std_err=se)
end

export
    var2,
    waic


text:
module MountainCarStepsExperiment

using DeepRL
using Flux
using Random
using ProgressMeter
using Plots
using TensorBoardLogger
using Logging
using LinearAlgebra
using BSON

glorot_uniform(rng::Random.AbstractRNG, dims...) = (rand(rng, Float32, dims...) .- 0.5f0) .* sqrt(24.0f0/sum(dims))
glorot_normal(rng::Random.AbstractRNG, dims...) = randn(rng, Float32, dims...) .* sqrt(2.0f0/sum(dims))

function construct_agent(s, num_actions, tn_counter_init = 50)

    ϵ=0.1
    γ=1.0f0
    batch_size=32
    # tn_counter_init = 50

    model = Chain(Dense(length(s), 128, Flux.relu; initW=Flux.glorot_uniform),
                  Dense(128, 128, Flux.relu; initW=Flux.glorot_uniform),
                  Dense(128, 32, Flux.relu; initW=Flux.glorot_uniform),
                  Dense(32, num_actions; initW=Flux.glorot_uniform))

    target_network = mapleaves(Flux.Tracker.data, deepcopy(model)::typeof(model))

    return DQNAgent(model,
                    target_network,
                    ADAM(0.001),
                    QLearning(γ),
                    ϵGreedy(ϵ),
                    1000000,
                    γ,
                    batch_size,
                    tn_counter_init,
                    s)
end

function plot_layers(agent::DQNAgent, data_range)

    y = maximum.(Flux.data.(agent.model.(data_range)))

end

function get_heatmap(agent::DQNAgent, data_range; kwargs...)
    data = collect.(collect(Iterators.product(data_range...)))
    ret = Flux.data(agent.model(hcat(reshape(data, length(data), 1)...)))
    ret_max_action = reshape(getindex.(findmax(ret; dims=1)[2], 1), size(data)...)

    return heatmap(
        collect(data_range[1]),
        collect(data_range[2]),
        ret_max_action; kwargs...)
end

mutable struct HeatmapAnimation
    anim::Animation
    counter_init::Int64
    counter::Int64
    display::Bool
    num_frames::Int64
    HeatmapAnimation(anim, counter_init, display) = new(anim, counter_init, counter_init, display, 0)
end

HeatmapAnimation(counter_init::Int64, display=false) = HeatmapAnimation(Animation(), counter_init, display)

function add_heatmap!(hmanim::HeatmapAnimation, agent::DQNAgent, data_range)
    if hmanim.counter == 1
        plt = get_heatmap(agent, data_range; title="$(hmanim.num_frames*hmanim.counter_init)", clim=(1.0,3.0))
        hmanim.num_frames += 1
        if hmanim.display
            gui(plt)
        else
            frame(hmanim.anim, plt)
        end
        hmanim.counter = hmanim.counter_init
    else
        hmanim.counter -= 1
    end
end

mutable struct ModelSnapshot
    loc::String
    base_file::String
    counter_init::Int64
    counter::Int64
    num_snapshots::Int64
    function ModelSnapshot(loc::AbstractString, snapshot_time::Int64, base_file::AbstractString="model")
        if !isdir(loc)
            mkdir(loc)
        else
            throw("Directory $(loc) already exists")
        end
        new(loc, base_file, snapshot_time, 1, 0)
    end
end

function take_snapshot!(mss::ModelSnapshot, agent::DQNAgent)
    if mss.counter == 1
        bson(joinpath(mss.loc, "$(mss.base_file)_$(mss.num_snapshots*mss.counter_init).bson"), model=agent.model)
        mss.counter = mss.counter_init
        mss.num_snapshots += 1
    else
        mss.counter -= 1
    end
end

function episode!(env, agent, rng, max_steps, anim, data_range, mss)
    terminal = false
    s_t = start!(env, rng)
    action = start!(agent, s_t, rng)
    add_heatmap!(anim, agent, data_range)
    take_snapshot!(mss, agent)

    total_rew = 0
    steps = 0

    while !terminal

        s_tp1, rew, terminal = step!(env, action)
        if steps == max_steps
            terminal = true
        end
        action = step!(agent, s_tp1, rew, terminal, rng)
        total_rew += rew
        steps += 1
        add_heatmap!(anim, agent, data_range)
        take_snapshot!(mss, agent)
    end
    return total_rew, steps
end

function main_experiment(seed, num_episodes, tn_counter=50; gif_save_file="my_gif.gif", model_snapshot_loc="dqn_snapshots")

    # lg=TBLogger("tensorboard_logs/run", min_level=Logging.Info)

    mc = MountainCar(0.0, 0.0, true)
    Random.seed!(Random.GLOBAL_RNG, seed)

    s = JuliaRL.get_state(mc)

    agent = construct_agent(s, length(JuliaRL.get_actions(mc)), tn_counter)::DQNAgent

    total_rews = zeros(num_episodes)
    steps = zeros(Int64, num_episodes)

    front = ['▁' ,'▂' ,'▃' ,'▄' ,'▅' ,'▆', '▇']
    p = ProgressMeter.Progress(
        num_episodes;
        dt=0.01,
        desc="Episode:",
        barglyphs=ProgressMeter.BarGlyphs('|','█',front,' ','|'),
        barlen=Int64(floor(300/length(front))))

    gr()
    if !isdir("plts")
        mkdir("plts")
    end

    data_ranges = (0.0:0.005:1.0, 0.0:0.005:1.0)
    anim = HeatmapAnimation(50, false)
    mss = ModelSnapshot(model_snapshot_loc, 50)

    # with_logger(lg) do
    for e in 1:num_episodes
        total_rews[e], steps[e] = episode!(mc, agent, Random.GLOBAL_RNG, 50000, anim, data_ranges, mss)
        next!(p)
    end

    if !anim.display
        gif(anim.anim, gif_save_file)
    end
    # end

    return total_rews, steps

end

end


text:
#=
    Sorting networks that process 1,2,..16 values are given.
    Sorting networks that process 17..25 values are given in "swapsort_17to25.jl".

    These sorting networks should minimize conditional swaps.
    The first sixteen are known to minimize conditional swaps.

    Functions are written with groups of internally parallelizable
    statements given together and sequential steps are separated.

    `swapsort(x1::T, .., xn::T)::NTuple{n,T}`
    `swapsort(tup::NTuple{n,T})::NTuple{n,T}`
    `swapsort(vec::Vector{T})::NTuple{N,T}`

    networks were selected using software by John Gamble
        http://pages.ripco.net/~jgamble/nw.html
        http://search.cpan.org/dist/Algorithm-Networksort/
=#

const VALS = ([Val{i} for i=1:ITEMS_MAX]...,)

@inline function swapsort(vec::Vector{T}) where {T}
    n = length(vec)
    n > ITEMS_MAX && throw(ErrorException("swapsort(vector) not defined where length(vector) > $ITEMS_MAX"))
    return swapsort(VALS[length(vec)], vec)
end

#    sort 1 value with 0 minmaxs in 1 stage

function swapsort(a::T) where {T}
    return (a,)
end

@inline swapsort(x::NTuple{1,T}) where {T} = swapsort(x[1])
@inline swapsort(::Type{Val{1}}, x::Vector{T}) where {T} = swapsort(x[1])

#    sort 2 values with 1 minmaxs in 1 stage

function swapsort(a::T, b::T) where {T}

    a, b = min_max(a, b)

    return a,b
end

@inline swapsort(x::NTuple{2,T}) where {T} = swapsort(x[1], x[2])
@inline swapsort(::Type{Val{2}}, x::Vector{T}) where {T} = swapsort(x[1], x[2])

#    sort 3 values with 3 minmaxs in 3 parallel stages

function swapsort(a::T, b::T, c::T) where {T}

    b, c = min_max(b, c)

    a, c = min_max(a, c)

    a, b = min_max(a, b)

    return a, b, c
end

@inline swapsort(x::NTuple{3,T}) where {T} = swapsort(x[1], x[2], x[3])
@inline swapsort(::Type{Val{3}}, x::Vector{T}) where {T} = swapsort(x[1], x[2], x[3])

#    sort 4 values with 5 minmaxs in 3 parallel stages

function swapsort(a::T, b::T, c::T, d::T) where {T}

    a, b = min_max(a, b)
    c, d = min_max(c, d)

    a, c = min_max(a, c)
    b, d = min_max(b, d)

    b, c = min_max(b, c)

    return a, b, c, d
end

@inline swapsort(x::NTuple{4,T}) where {T} = swapsort(x[1], x[2], x[3], x[4])
@inline swapsort(::Type{Val{4}}, x::Vector{T}) where {T} = swapsort(x[1], x[2], x[3], x[4])

#    sort 5 values with 9 minmaxs in 5 parallel stages

function swapsort(a::T, b::T, c::T, d::T, e::T) where {T}

    b, c = min_max(b, c)
    d, e = min_max(d, e)

    b, d = min_max(b, d)
    a, c = min_max(a, c)

    c, e = min_max(c, e)
    a, d = min_max(a, d)

    a, b = min_max(a, b)
    c, d = min_max(c, d)

    b, c = min_max(b, c)

    return a, b, c, d, e
end

@inline swapsort(x::NTuple{5,T}) where {T} = swapsort(x[1], x[2], x[3], x[4], x[5])
@inline swapsort(::Type{Val{5}}, x::Vector{T}) where {T} = swapsort(x[1], x[2], x[3], x[4], x[5])

#    sort 6 values with 12 minmaxs in 6 parallel stages

function swapsort(a::T, b::T, c::T, d::T, e::T,
              f::T) where {T}

    b, c = min_max(b, c)
    e, f = min_max(e, f)

    a, c = min_max(a, c)
    d, f = min_max(d, f)

    a, b = min_max(a, b)
    d, e = min_max(d, e)
    c, f = min_max(c, f)

    a, d = min_max(a, d)
    b, e = min_max(b, e)

    c, e = min_max(c, e)
    b, d = min_max(b, d)

    c, d = min_max(c, d)

    return a, b, c, d, e, f
end

@inline swapsort(x::NTuple{6,T}) where {T} = swapsort(x[1], x[2], x[3], x[4], x[5], x[6])
@inline swapsort(::Type{Val{6}}, x::Vector{T}) where {T} =
    swapsort(x[1], x[2], x[3], x[4], x[5], x[6])

#    sort 7 values with 16 minmaxs in 6 parallel stages

function swapsort(a::T, b::T, c::T, d::T, e::T,
              f::T, g::T) where {T}

    a, e = min_max(a, e)
    b, f = min_max(b, f)
    c, g = min_max(c, g)

    a, c = min_max(a, c)
    b, d = min_max(b, d)
    e, g = min_max(e, g)

    c, e = min_max(c, e)
    d, f = min_max(d, f)
    a, b = min_max(a, b)

    c, d = min_max(c, d)
    e, f = min_max(e, f)

    b, e = min_max(b, e)
    d, g = min_max(d, g)

    b, c = min_max(b, c)
    d, e = min_max(d, e)
    f, g = min_max(f, g)

    return a, b, c, d, e, f, g
end

@inline swapsort(x::NTuple{7,T}) where {T} = swapsort(x[1], x[2], x[3], x[4], x[5], x[6], x[7])
@inline swapsort(::Type{Val{7}}, x::Vector{T}) where {T} =
    swapsort(x[1], x[2], x[3], x[4], x[5], x[6], x[7])

#    sort 8 values with 19 minmaxs in 7 parallel stages

function swapsort(a::T, b::T, c::T, d::T, e::T,
              f::T, g::T, h::T) where {T}

    a, b = min_max(a, b)
    c, d = min_max(c, d)
    e, f = min_max(e, f)
    g, h = min_max(g, h)

    a, c = min_max(a, c)
    b, d = min_max(b, d)
    e, g = min_max(e, g)
    f, h = min_max(f, h)

    b, c = min_max(b, c)
    f, g = min_max(f, g)
    a, e = min_max(a, e)
    d, h = min_max(d, h)

    b, f = min_max(b, f)
    c, g = min_max(c, g)

    b, e = min_max(b, e)
    d, g = min_max(d, g)

    c, e = min_max(c, e)
    d, f = min_max(d, f)

    d, e = min_max(d, e)

    return a, b, c, d, e, f, g, h
end

@inline swapsort(x::NTuple{8,T}) where {T} = swapsort(x[1], x[2], x[3], x[4], x[5], x[6], x[7], x[8])
@inline swapsort(::Type{Val{8}}, x::Vector{T}) where {T} =
    swapsort(x[1], x[2], x[3], x[4], x[5], x[6], x[7], x[8])

#    sort 9 values with 25 minmaxs in 9 parallel stages

function swapsort(a::T, b::T, c::T, d::T, e::T,
              f::T, g::T, h::T, i::T) where {T}

    a, b = min_max(a, b)
    d, e = min_max(d, e)
    g, h = min_max(g, h)

    b, c = min_max(b, c)
    e, f = min_max(e, f)
    h, i = min_max(h, i)

    a, b = min_max(a, b)
    d, e = min_max(d, e)
    g, h = min_max(g, h)
    c, f = min_max(c, f)

    a, d = min_max(a, d)
    b, e = min_max(b, e)
    f, i = min_max(f, i)

    d, g = min_max(d, g)
    e, h = min_max(e, h)
    c, f = min_max(c, f)

    a, d = min_max(a, d)
    b, e = min_max(b, e)
    f, h = min_max(f, h)
    c, g = min_max(c, g)

    b, d = min_max(b, d)
    e, g = min_max(e, g)

    c, e = min_max(c, e)
    f, g = min_max(f, g)

    c, d = min_max(c, d)

    return a, b, c, d, e, f, g, h, i
end

@inline swapsort(x::NTuple{9,T}) where {T} = swapsort(x[1], x[2], x[3], x[4], x[5], x[6], x[7], x[8], x[9])
@inline swapsort(::Type{Val{9}}, x::Vector{T}) where {T} =
    swapsort(x[1], x[2], x[3], x[4], x[5], x[6], x[7], x[8], x[9])

#    sort 10 values with 29 minmaxs in 9 parallel stages

function swapsort(a::T,b::T,c::T,d::T,e::T,
              f::T,g::T,h::T,i::T,j::T) where {T}

    e, j = min_max(e, j)
    d, i = min_max(d, i)
    c, h = min_max(c, h)
    b, g = min_max(b, g)
    a, f = min_max(a, f)

    b, e = min_max(b, e)
    g, j = min_max(g, j)
    a, d = min_max(a, d)
    f, i = min_max(f, i)

    a, c = min_max(a, c)
    d, g = min_max(d, g)
    h, j = min_max(h, j)

    a, b = min_max(a, b)
    c, e = min_max(c, e)
    f, h = min_max(f, h)
    i, j = min_max(i, j)

    b, c = min_max(b, c)
    e, g = min_max(e, g)
    h, i = min_max(h, i)
    d, f = min_max(d, f)

    c, f = min_max(c, f)
    g, i = min_max(g, i)
    b, d = min_max(b, d)
    e, h = min_max(e, h)

    c, d = min_max(c, d)
    g, h = min_max(g, h)

    d, e = min_max(d, e)
    f, g = min_max(f, g)

    e, f = min_max(e, f)

    return a,b,c,d,e,f,g,h,i,j
end

@inline swapsort(x::NTuple{10,T}) where {T} =
    swapsort(x[1], x[2], x[3], x[4], x[5], x[6], x[7], x[8], x[9], x[10])
@inline swapsort(::Type{Val{10}}, x::Vector{T}) where {T} =
    swapsort(x[1], x[2], x[3], x[4], x[5], x[6], x[7], x[8], x[9], x[10])

#    sort 11 values with 35 minmaxs in 9 parallel stages

function swapsort(a::T, b::T, c::T, d::T, e::T,
              f::T, g::T, h::T, i::T, j::T,
              k::T) where {T}

    a, b = min_max(a, b)
    c, d = min_max(c, d)
    e, f = min_max(e, f)
    g, h = min_max(g, h)
    i, j = min_max(i, j)

    b, d = min_max(b, d)
    f, h = min_max(f, h)
    a, c = min_max(a, c)
    e, g = min_max(e, g)
    i, k = min_max(i, k)

    b, c = min_max(b, c)
    f, g = min_max(f, g)
    j, k = min_max(j, k)
    a, e = min_max(a, e)
    d, h = min_max(d, h)

    b, f = min_max(b, f)
    g, k = min_max(g, k)
    e, i = min_max(e, i)

    f, j = min_max(f, j)
    c, g = min_max(c, g)
    a, e = min_max(a, e)
    d, i = min_max(d, i)

    b, f = min_max(b, f)
    g, k = min_max(g, k)
    c, d = min_max(c, d)
    i, j = min_max(i, j)

    b, e = min_max(b, e)
    h, k = min_max(h, k)
    d, f = min_max(d, f)
    g, i = min_max(g, i)

    c, e = min_max(c, e)
    h, j = min_max(h, j)
    f, g = min_max(f, g)

    d, e = min_max(d, e)
    h, i = min_max(h, i)

    return a,b,c,d,e,f,g,h,i,j,k
end

@inline swapsort(x::NTuple{11,T}) where {T} =
    swapsort(x[1], x[2], x[3], x[4], x[5], x[6], x[7], x[8], x[9], x[10],
             x[11])
@inline swapsort(::Type{Val{11}}, x::Vector{T}) where {T} =
    swapsort(x[1], x[2], x[3], x[4], x[5], x[6], x[7], x[8], x[9], x[10],
             x[11])

#    sort 12 values with 39 minmaxs in 9 parallel stages

function swapsort(a::T,b::T,c::T,d::T,e::T,
              f::T,g::T,h::T,i::T,j::T,
              k::T,l::T) where {T}

    a, b = min_max(a, b)
    c, d = min_max(c, d)
    e, f = min_max(e, f)
    g, h = min_max(g, h)
    i, j = min_max(i, j)
    k, l = min_max(k, l)

    b, d = min_max(b, d)
    f, h = min_max(f, h)
    j, l = min_max(j, l)
    a, c = min_max(a, c)
    e, g = min_max(e, g)
    i, k = min_max(i, k)

    b, c = min_max(b, c)
    f, g = min_max(f, g)
    j, k = min_max(j, k)
    a, e = min_max(a, e)
    h, l = min_max(h, l)

    b, f = min_max(b, f)
    g, k = min_max(g, k)
    d, h = min_max(d, h)
    e, i = min_max(e, i)

    f, j = min_max(f, j)
    c, g = min_max(c, g)
    a, e = min_max(a, e)
    h, l = min_max(h, l)
    d, i = min_max(d, i)

    b, f = min_max(b, f)
    g, k = min_max(g, k)
    c, d = min_max(c, d)
    i, j = min_max(i, j)

    b, e = min_max(b, e)
    h, k = min_max(h, k)
    d, f = min_max(d, f)
    g, i = min_max(g, i)

    c, e = min_max(c, e)
    h, j = min_max(h, j)
    f, g = min_max(f, g)

    d, e = min_max(d, e)
    h, i = min_max(h, i)

    return a,b,c,d,e,f,g,h,i,j,k,l
end

@inline swapsort(x::NTuple{12,T}) where {T} =
    swapsort(x[1], x[2], x[3], x[4], x[5], x[6], x[7], x[8], x[9], x[10],
             x[11], x[12])
@inline swapsort(::Type{Val{12}}, x::Vector{T}) where {T} =
    swapsort(x[1], x[2], x[3], x[4], x[5], x[6], x[7], x[8], x[9], x[10],
             x[11], x[12])

#    sort 13 values with 45 minmaxs in 10 parallel stages

function swapsort(a::T, b::T, c::T, d::T, e::T,
              f::T, g::T, h::T, i::T, j::T,
              k::T, l::T, m::T) where {T}

    b, h = min_max(b, h)
    j, l = min_max(j, l)
    d, e = min_max(d, e)
    f, i = min_max(f, i)
    a, m = min_max(a, m)
    c, g = min_max(c, g)

    a, b = min_max(a, b)
    c, d = min_max(c, d)
    e, g = min_max(e, g)
    i, l = min_max(i, l)
    h, m = min_max(h, m)
    f, j = min_max(f, j)

    a, c = min_max(a, c)
    d, h = min_max(d, h)
    k, l = min_max(k, l)
    b, e = min_max(b, e)
    g, m = min_max(g, m)

    h, i = min_max(h, i)
    l, m = min_max(l, m)
    e, j = min_max(e, j)
    g, k = min_max(g, k)

    d, e = min_max(d, e)
    f, g = min_max(f, g)
    i, j = min_max(i, j)
    k, l = min_max(k, l)
    b, h = min_max(b, h)

    c, g = min_max(c, g)
    j, l = min_max(j, l)
    b, d = min_max(b, d)
    e, h = min_max(e, h)
    i, k = min_max(i, k)
    a, f = min_max(a, f)

    c, f = min_max(c, f)
    g, i = min_max(g, i)
    j, k = min_max(j, k)

    b, c = min_max(b, c)
    d, f = min_max(d, f)
    h, i = min_max(h, i)
    e, g = min_max(e, g)

    c, d = min_max(c, d)
    e, f = min_max(e, f)
    g, h = min_max(g, h)
    i, j = min_max(i, j)

    d, e = min_max(d, e)
    f, g = min_max(f, g)

    return a,b,c,d,e,f,g,h,i,j,k,l,m
end

@inline swapsort(x::NTuple{13,T}) where {T} =
    swapsort(x[1], x[2], x[3], x[4], x[5], x[6], x[7], x[8], x[9], x[10],
             x[11], x[12], x[13])
@inline swapsort(::Type{Val{13}}, x::Vector{T}) where {T} =
    swapsort(x[1], x[2], x[3], x[4], x[5], x[6], x[7], x[8], x[9], x[10],
             x[11], x[12], x[13])

#    sort 14 values with 51 minmaxs in 10 parallel stages

function swapsort(a::T, b::T, c::T, d::T, e::T,
              f::T, g::T, h::T, i::T, j::T,
              k::T, l::T, m::T, n::T) where {T}

    a, b = min_max(a, b)
    c, d = min_max(c, d)
    e, f = min_max(e, f)
    g, h = min_max(g, h)
    i, j = min_max(i, j)
    k, l = min_max(k, l)
    m, n = min_max(m, n)

    a, c = min_max(a, c)
    e, g = min_max(e, g)
    i, k = min_max(i, k)
    b, d = min_max(b, d)
    f, h = min_max(f, h)
    j, l = min_max(j, l)

    a, e = min_max(a, e)
    i, m = min_max(i, m)
    b, f = min_max(b, f)
    j, n = min_max(j, n)
    c, g = min_max(c, g)
    d, h = min_max(d, h)

    a, i = min_max(a, i)
    b, j = min_max(b, j)
    c, k = min_max(c, k)
    d, l = min_max(d, l)
    e, m = min_max(e, m)
    f, n = min_max(f, n)

    f, k = min_max(f, k)
    g, j = min_max(g, j)
    d, m = min_max(d, m)
    h, l = min_max(h, l)
    b, c = min_max(b, c)
    e, i = min_max(e, i)

    b, e = min_max(b, e)
    h, n = min_max(h, n)
    c, i = min_max(c, i)
    f, g = min_max(f, g)
    j, k = min_max(j, k)

    c, e = min_max(c, e)
    l, n = min_max(l, n)
    d, i = min_max(d, i)
    h, m = min_max(h, m)

    g, i = min_max(g, i)
    k, m = min_max(k, m)
    d, f = min_max(d, f)
    h, j = min_max(h, j)

    d, e = min_max(d, e)
    f, g = min_max(f, g)
    h, i = min_max(h, i)
    j, k = min_max(j, k)
    l, m = min_max(l, m)

    g, h = min_max(g, h)
    i, j = min_max(i, j)

    return a,b,c,d,e,f,g,h,i,j,k,l,m,n
end

@inline swapsort(x::NTuple{14,T}) where {T} =
    swapsort(x[1], x[2], x[3], x[4], x[5], x[6], x[7], x[8], x[9], x[10],
             x[11], x[12], x[13], x[14])
@inline swapsort(::Type{Val{14}}, x::Vector{T}) where {T} =
    swapsort(x[1], x[2], x[3], x[4], x[5], x[6], x[7], x[8], x[9], x[10],
             x[11], x[12], x[13], x[14])

#    sort 15 values with 56 minmaxs in 10 parallel stages

function swapsort(a::T, b::T, c::T, d::T, e::T,
              f::T, g::T, h::T, i::T, j::T,
              k::T, l::T, m::T, n::T, o::T) where {T}

    a, b = min_max(a, b)
    c, d = min_max(c, d)
    e, f = min_max(e, f)
    g, h = min_max(g, h)
    i, j = min_max(i, j)
    k, l = min_max(k, l)
    m, n = min_max(m, n)

    a, c = min_max(a, c)
    e, g = min_max(e, g)
    i, k = min_max(i, k)
    m, o = min_max(m, o)
    b, d = min_max(b, d)
    f, h = min_max(f, h)
    j, l = min_max(j, l)

    a, e = min_max(a, e)
    i, m = min_max(i, m)
    b, f = min_max(b, f)
    j, n = min_max(j, n)
    c, g = min_max(c, g)
    k, o = min_max(k, o)
    d, h = min_max(d, h)

    a, i = min_max(a, i)
    b, j = min_max(b, j)
    c, k = min_max(c, k)
    d, l = min_max(d, l)
    e, m = min_max(e, m)
    f, n = min_max(f, n)
    g, o = min_max(g, o)

    f, k = min_max(f, k)
    g, j = min_max(g, j)
    d, m = min_max(d, m)
    n, o = min_max(n, o)
    h, l = min_max(h, l)
    b, c = min_max(b, c)
    e, i = min_max(e, i)

    b, e = min_max(b, e)
    h, n = min_max(h, n)
    c, i = min_max(c, i)
    l, o = min_max(l, o)
    f, g = min_max(f, g)
    j, k = min_max(j, k)

    c, e = min_max(c, e)
    l, n = min_max(l, n)
    d, i = min_max(d, i)
    h, m = min_max(h, m)

    g, i = min_max(g, i)
    k, m = min_max(k, m)
    d, f = min_max(d, f)
    h, j = min_max(h, j)

    d, e = min_max(d, e)
    f, g = min_max(f, g)
    h, i = min_max(h, i)
    j, k = min_max(j, k)
    l, m = min_max(l, m)

    g, h = min_max(g, h)
    i, j = min_max(i, j)

    return a,b,c,d,e,f,g,h,i,j,k,l,m,n,o
end

@inline swapsort(x::NTuple{15,T}) where {T} =
    swapsort(x[1], x[2], x[3], x[4], x[5], x[6], x[7], x[8], x[9], x[10],
             x[11], x[12], x[13], x[14], x[15])
@inline swapsort(::Type{Val{15}}, x::Vector{T}) where {T} =
    swapsort(x[1], x[2], x[3], x[4], x[5], x[6], x[7], x[8], x[9], x[10],
             x[11], x[12], x[13], x[14], x[15])

#    sort 16 values with 60 minmaxs in 10 parallel stages

function swapsort(a::T, b::T, c::T, d::T, e::T,
              f::T, g::T, h::T, i::T, j::T,
              k::T, l::T, m::T, n::T, o::T,
              p::T) where {T}

    a, b = min_max(a, b)
    c, d = min_max(c, d)
    e, f = min_max(e, f)
    g, h = min_max(g, h)
    i, j = min_max(i, j)
    k, l = min_max(k, l)
    m, n = min_max(m, n)
    o, p = min_max(o, p)

    a, c = min_max(a, c)
    e, g = min_max(e, g)
    i, k = min_max(i, k)
    m, o = min_max(m, o)
    b, d = min_max(b, d)
    f, h = min_max(f, h)
    j, l = min_max(j, l)
    n, p = min_max(n, p)

    a, e = min_max(a, e)
    i, m = min_max(i, m)
    b, f = min_max(b, f)
    j, n = min_max(j, n)
    c, g = min_max(c, g)
    k, o = min_max(k, o)
    d, h = min_max(d, h)
    l, p = min_max(l, p)

    a, i = min_max(a, i)
    b, j = min_max(b, j)
    c, k = min_max(c, k)
    d, l = min_max(d, l)
    e, m = min_max(e, m)
    f, n = min_max(f, n)
    g, o = min_max(g, o)
    h, p = min_max(h, p)

    f, k = min_max(f, k)
    g, j = min_max(g, j)
    d, m = min_max(d, m)
    n, o = min_max(n, o)
    h, l = min_max(h, l)
    b, c = min_max(b, c)
    e, i = min_max(e, i)

    b, e = min_max(b, e)
    h, n = min_max(h, n)
    c, i = min_max(c, i)
    l, o = min_max(l, o)
    f, g = min_max(f, g)
    j, k = min_max(j, k)

    c, e = min_max(c, e)
    l, n = min_max(l, n)
    d, i = min_max(d, i)
    h, m = min_max(h, m)

    g, i = min_max(g, i)
    k, m = min_max(k, m)
    d, f = min_max(d, f)
    h, j = min_max(h, j)

    d, e = min_max(d, e)
    f, g = min_max(f, g)
    h, i = min_max(h, i)
    j, k = min_max(j, k)
    l, m = min_max(l, m)

    g, h = min_max(g, h)
    i, j = min_max(i, j)

    return a,b,c,d,e,f,g,h,i,j,k,l,m,n,o,p
end

@inline swapsort(x::NTuple{16,T}) where {T} =
    swapsort(x[1], x[2], x[3], x[4], x[5], x[6], x[7], x[8], x[9], x[10],
             x[11], x[12], x[13], x[14], x[15], x[16])
@inline swapsort(::Type{Val{16}}, x::Vector{T}) where {T} =
    swapsort(x[1], x[2], x[3], x[4], x[5], x[6], x[7], x[8], x[9], x[10],
             x[11], x[12], x[13], x[14], x[15], x[16])


text:
@testset "Plotting" begin
  Random.seed!(123)
  z₁ = randn(10000)
  z₂ = z₁ + randn(10000)
  d = georef(DataFrame(z₁=z₁,z₂=z₂), RegularGrid(100,100))

  @testset "DistPlot1D" begin
    if visualtests
      @plottest distplot1d(d,:z₁) joinpath(datadir,"distplot1D.png") !istravis
    end
  end

  @testset "DistPlot2D" begin
    if visualtests
      @plottest distplot2d(d,:z₁,:z₂) joinpath(datadir,"distplot2D.png") !istravis
    end
  end

  @testset "CornerPlot" begin
    if visualtests
      @plottest cornerplot(d) joinpath(datadir,"cornerplot.png") !istravis
    end
  end
end


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

using ScoreDrivenERGM
using LightGraphs

using Libdl
using RCall
R""" library("ergm") """

ergmC = Libdl.dlopen("C:\\Users\\digan\\Documents\\R\\win-library\\4.0\\ergm\\libs\\x64\\ergm")

# define test matrix
N = 100
A = rand(Bool,N,N)

g = SimpleDiGraph(A)

edgeList = collect(edges(g))
tails = getproperty.(edgeList, :src)
heads = getproperty.(edgeList, :dst)
dnedges = ne(g)
dflag = Int(is_directed(g))

# try to understand what is wl input
@rput A
R" net <- network(A)"
R"to_ergm_Cdouble(net)"

# store the sufficient statistics and change statistics in R
R"""
    net <- network(A)
    chStat_t <- ergmMPLE(formula_ergm)
    changeStats_t_R <- cbind(chStat_t$response, chStat_t$predictor,     chStat_t$weights)
        """

# create julia variables (Refs ) that can be passed to MPLE_wrapper fun

# test MPLE_wrapper function
tailsIn = Ref(tails)
headsIn = Ref(heads)
ccall(:GetRNGstate, Int32)

# MPLE_wrapper(int *tails, int *heads, int *dnedges,
# 		  double *wl,
# 		  int *dn, int *dflag, int *bipartite, int *nterms,
# 		  char **funnames, char **sonames, double *inputs,
# 		  int *responsevec, double *covmat,
# 		  int *weightsvector,
# 		  int *maxDyads, int *maxDyadTypes)

# ABANDONING THE IDEA BECAUSE IT WOULD REQUIRE TOO MUCH WORK FOR THE MOMENT


text:
using Altro
using TrajectoryOptimization
using Test
using StaticArrays
using ForwardDiff
using RobotDynamics
using LinearAlgebra
const RD = RobotDynamics
const TO = TrajectoryOptimization

function alcost(::TO.Equality, con, λ, μ, z)
    c = RD.evaluate(con, z)
    Iμ = Diagonal(μ)
    λ'c + 0.5*c'Iμ*c
end

function alcost(::TO.Inequality, con, λ, μ, z)
    c = RD.evaluate(con, z)
    a = (λ .> 0) .| (c .> 0)
    Iμ = Diagonal(a .* μ)
    λ'c + 0.5*c'Iμ*c
end

alcost(::TO.SecondOrderCone, con, λ, μ, z) = alcost(con, λ, μ, z)

function alcost(con, λ, μ, z)
    cone = TO.dualcone(TO.sense(con))
    c = RD.evaluate(con, z)
    λbar = λ .- μ .* c
    λp = TO.projection(cone, λbar)
    Iμ = Diagonal(inv.(μ))
    return 0.5*(λp'Iμ*λp - λ'Iμ*λ)
end

function algrad(::TO.Equality, con, λ, μ, ∇c, z)
    c = RD.evaluate(con, z)
    λbar = λ .+ μ .* c
    ∇c'λbar
end

function algrad(::TO.Inequality, con, λ, μ, ∇c, z)
    c = RD.evaluate(con, z)
    a = (c .>= 0) .| (λ .> 0)
    λbar = λ .+ (a .* μ) .* c
    ∇c'λbar
end

algrad(::TO.SecondOrderCone, con, λ, μ, ∇c, z) = algrad(con, λ, μ, ∇c, z)

function algrad(con, λ, μ, ∇c, z)
    cone = TO.dualcone(TO.sense(con))
    c = RD.evaluate(con, z)
    λbar = λ .-  μ .* c
    λp = TO.projection(cone, λbar)
    Iμ = Diagonal(μ)
    p = length(λ)
    ∇proj = zeros(eltype(z), p,p)
    TO.∇projection!(cone, ∇proj, λbar)
    -∇c'Iμ*∇proj'*(Iμ\λp)
end

##
T = Float64
n, m, N = 3, 2, 11
tf = 2.0
xf = [1,2,3.]

con_eq = TO.GoalConstraint(xf, SA[1,3])
inds_eq = 1:10
confun_eq(x) = (x - xf)[[1,3]]

xc, yc, rad = [4,2,2,0.], [0,1,2,3.], [0.2, 0.1, 0.3, 0.4]
con_in = TO.CircleConstraint(n, xc, yc, rad)
inds_in = 2:6
confun_in(x) = rad.^2 .- (x[1] .- xc).^2 .- (x[2] .- yc).^2

con_so = TO.NormConstraint(n, m, 3, TO.SecondOrderCone(), :state)
inds_so = N:N
confun_so(x) = [x; 3]

# Other data needed for constructor
Z = SampledTrajectory{n,m}(randn(n,N), randn(m,N-1), dt=0.1)
opts = SolverOptions()
costs = zeros(N)

##
@testset "$(TO.conename(TO.sense(con))) ALConstraint" for (con,inds,confun) in [
    (con_eq, inds_eq, confun_eq),
    (con_in, inds_in, confun_in),
    (con_so, inds_so, confun_so),
]
# con, inds, confun =
#     (con_so, inds_so, confun_so)
alcon = Altro.ALConstraint{T}(Z, con, inds, costs, solveropts=opts)

# Test that it got filled in correctly
p = TO.output_dim(alcon.con)
P = length(inds)
@test alcon.nx[1] == n
@test alcon.nu[1] == m
@test alcon.inds == inds
for field in (:vals, :λ, :μ, :μinv, :λbar, :λproj, :λscaled, :viol)
    @test all(x->length(x)==p, getfield(alcon, field))
    @test length(getfield(alcon, field)) == P
end
w = RD.input_dim(con)
@test size(alcon.c_max) == (N,)
@test length(alcon.jac) == P
@test all(x->size(x)==(p,w), alcon.jac)
@test length(alcon.∇proj) == P
@test all(x->size(x)==(p,p), alcon.∇proj)
@test length(alcon.∇²proj) == P
@test all(x->size(x)==(p,p), alcon.∇²proj)
@test length(alcon.grad) == P
@test all(x->size(x)==(n+m,), alcon.grad)
@test length(alcon.hess) == P
@test all(x->size(x)==(n+m,n+m), alcon.hess)
@test size(alcon.tmp_jac) == (p,w)
@test all(x->x == ones(p), alcon.μ)

# Test evaluation methods
Altro.settraj!(alcon, Z)
Altro.evaluate_constraints!(alcon)
c = map(Z[inds]) do z
    confun(TO.state(z))
    # (TO.state(z) - xf)[[1,3]]
end
@test c ≈ alcon.vals

Altro.constraint_jacobians!(alcon)
C = map(Z[inds]) do z
    J = ForwardDiff.jacobian(confun, TO.state(z))
    [J zeros(p,w-n)]
end
@test C ≈ alcon.jac

# Set non-zero duals
for i = 1:P
    alcon.λ[i] .= randn(p)
end
λ = alcon.λ
μ = alcon.μ

# Test alcost
alcon.opts.use_conic_cost = false
cone = TO.sense(con)

Altro.evaluate_constraints!(alcon)
for i = 1:P
    Altro.alcost(alcon, i)
    Altro.algrad!(cone, alcon, i)
    z = Z[inds[i]]
    _z = RD.StaticKnotPoint{n,m}(MVector{n+m}(z.z), z.t, z.dt)
    grad = ForwardDiff.gradient(
        x->alcost(cone, alcon.con, SVector{p}(λ[i]), SVector{p}(μ[i]), x),
        _z
    )
    @test grad ≈ alcon.grad[i]
    @test alcon.grad[i][1:n] ≈ algrad(
        cone, alcon.con, SVector{p}(λ[i]), SVector{p}(μ[i]), alcon.jac[i], _z
    )[1:n]
end

for i = 1:P
    cone = TO.sense(alcon.con)
    Altro.alhess!(cone, alcon, i)
    hess = zeros(n+m, n+m)
    # Get Gauss-Newton approximation
    z = Z[inds[i]]
    _z = RD.StaticKnotPoint{n,m}(MVector{n+m}(z.z), z.t, z.dt)
    ix = Altro.getinputinds(alcon, i)
    jac = ForwardDiff.jacobian(
        x->algrad(cone, alcon.con, SVector{p}(λ[i]), SVector{p}(μ[i]), alcon.jac[i], x),
        _z
    )
    hess[ix,:] .= jac
    @test hess[ix,ix] ≈ alcon.hess[i][ix,ix]
end

# Use Conic cost
Altro.evaluate_constraints!(alcon)
alcon.opts.use_conic_cost = true
for i = 1:P
    z = Z[inds[i]]
    _z = RD.StaticKnotPoint{n,m}(MVector{n+m}(z.z), z.t, z.dt)
    @test Altro.alcost(alcon, i) ≈ alcost(alcon.con, SVector{p}(λ[i]), SVector{p}(μ[i]), _z)
end

Altro.algrad!(alcon, 1)
for i = 1:P
    z = Z[inds[i]]
    _z = RD.StaticKnotPoint{n,m}(MVector{n+m}(z.z), z.t, z.dt)
    grad = ForwardDiff.gradient(
        x->alcost(alcon.con, SVector{p}(λ[i]), SVector{p}(μ[i]), x), _z
    )
    Altro.algrad!(alcon, i)
    @test alcon.grad[i] ≈ grad
    grad = algrad(alcon.con, SVector{p}(λ[i]), SVector{p}(μ[i]), alcon.jac[i], _z)
    @test alcon.grad[i][1:n] ≈ grad[1:n]
end

let i = P
    z = Z[inds[i]]
    _z = RD.StaticKnotPoint{n,m}(MVector{n+m}(z.z), z.t, z.dt)
    Altro.alhess!(alcon, i)
    hess = zeros(n+m, n+m)
    jac = ForwardDiff.jacobian(x->algrad(alcon.con, SVector{p}(λ[i]), SVector{p}(μ[i]),
        alcon.jac[i], x), _z
    )
    iz = Altro.getinputinds(alcon, i)
    hess[iz,:] .= jac
    @test alcon.hess[i] ≈ hess
end

## Dual and penalty updates
alcon.opts.use_conic_cost = false
λ0 = deepcopy(alcon.λ)
Altro.alcost(alcon)
Altro.dualupdate!(alcon)
for i = 1:P
    if cone == TO.Equality()
        @test λ[i] ≈ λ0[i] + μ[i] .* c[i]
    elseif cone == TO.Inequality()
        @test λ[i] ≈ clamp.(λ0[i] + μ[i] .* c[i], 0, Inf)
    else
        @test λ[i] ≈ alcon.λproj[i]
    end
    λ[i] .= λ0[i]  # reset them to what they were before
end

alcon.opts.use_conic_cost = true
Altro.dualupdate!(alcon)
for i = 1:P
    @test λ[i] ≈ alcon.λproj[i]
end

Altro.reset_penalties!(alcon)
Altro.penaltyupdate!(alcon)
@test Altro.max_penalty(alcon) == 10.0

Altro.reset_penalties!(alcon)
@test Altro.max_penalty(alcon) == 1.0

Altro.reset_duals!(alcon)
@test norm(λ) ≈ 0

end


text:
export genepop

"""
    genepop(infile::String; kwargs...)
Load a Genepop format file into memory as a PopData object.
- `infile::String` : path to Genepop file

### Keyword Arguments
- `digits::Integer`: number of digits denoting each allele (default: `3`)
- `popsep::String` : word that separates populations in `infile` (default: "POP")
- `silent::Bool`   : whether to print file information during import (default: `false`)
- `allow_monomorphic::Bool` : whether to keep monomorphic loci in the dataset (default: `false`)

### File must follow standard Genepop formatting:
- First line is a comment (and skipped)
- Loci are listed after first line as one-per-line without commas or in single comma-separated row
- A line with a particular keyword (default `POP`) must delimit populations
- Sample name is immediately followed by a *comma*
- File is *tab or space delimted* (but not both!)

### Genepop file example:
```
wasp_hive.gen: Wasp populations in New York
Locus1
Locus2
Locus3
pop
Oneida_01,  250230  564568  110100
Oneida_02,  252238  568558  100120
Oneida_03,  254230  564558  090100
pop
Newcomb_01, 254230  564558  080100
Newcomb_02, 000230  564558  090080
Newcomb_03, 254230  000000  090100
Newcomb_04, 254230  564000  090120
```
## Example
```
waspsNY = genepop("wasp_hive.gen", digits = 3, popsep = "pop")
```
"""
function genepop(
    infile::String;
    digits::Int = 3,
    popsep::String = "POP",
    silent::Bool = false,
    allow_monomorphic::Bool = false
)
    isfile(infile) || throw(ArgumentError("$infile not found."))
    # open the file as lines of strings to suss out loci names, pop idx, and popcounts
    gpop_readlines = readlines(infile)

    # find the #samples per population
    pop_idx = Vector{Int}()
    for (line, gtext) in enumerate(gpop_readlines)
        gtext == popsep && push!(pop_idx, line)
    end
    length(pop_idx) == 0 &&
    error("No populations found in $infile using separator \"$popsep\". Please check the spelling and try again.")

    # create a theoretical place were the last popsep would be (preserve last population for counting)
    push!(pop_idx, countlines(infile) + 1)
    popcounts = (pop_idx[2:end] .- pop_idx[1:end-1]) .- 1

    # check for the delimiter type in the first sample record
    firstrecord = gpop_readlines[pop_idx[1]+1]
    if occursin("\t", firstrecord) & occursin(" ", firstrecord)
        error("$infile contains both tab and space delimiters. Please format the file so it uses either one or the other.")
    elseif occursin("\t", firstrecord)
        delim = "\t"
        delim_txt = "tab"
    elseif occursin(" ", firstrecord)
        delim = " "
        delim_txt = "space"
    else
        error("Please format $infile to be either tab or space delimited")
    end

    # check for horizontal formatting, where popsep would appear on the third line
    if pop_idx[1] <= 3
        # second line should have all the loci
        locus_name_raw = replace(gpop_readlines[2], "." => "_")
        locinames = strip.(split(locus_name_raw |> join, ","))
        format = "horizontal"
    else
        #standard  vertical formatting
        locinames = gpop_readlines[2:pop_idx[1]-1]
        format = "vertical"
    end

    if !silent
        @info "\n $(truncatepath(abspath(infile)))\n formatting: delimiter = $(delim_txt), loci = $(format)\n data: loci = $(length(locinames)), samples = $(sum(popcounts)), populations = $(length(popcounts))"
        println()
    end

    # load in samples and genotypes
    pushfirst!(locinames, "name")

    geno_parse = CSV.read(
        infile,
        DataFrame,
        delim = delim,
        header = locinames,
        skipto = pop_idx[1] + 1,
        comment = popsep,
        missingstring = ["-9", ""],
        normalizenames = true,
        ignorerepeated = true
    )

    popnames = string.(1:length(popcounts))
    popnames = fill.(popnames,popcounts) |> Base.Iterators.flatten |> collect
    insertcols!(geno_parse, 2, :population => popnames)
    geno_parse.name .= strip.(geno_parse.name, ',')
    # wide to long format
    geno_parse = DataFrames.stack(geno_parse, DataFrames.Not([1,2]))
    rename!(geno_parse, [:name, :population, :locus, :genotype])
    select!(
        geno_parse,
        :name => (i -> PooledArray(i, compress = true)) => :name,
        :population => (i -> PooledArray(i, compress = true)) => :population,
        :locus => (i -> PooledArray(Array(i), compress = true)) => :locus,
        :genotype
    )

    # try to compress the alleles into Int8 (snps) or In16 (msat)
    try
        geno_parse.genotype = PooledArray(phase.(geno_parse.genotype, Int8, digits), compress = true)
    catch
        geno_parse.genotype = phase.(geno_parse.genotype, Int16, digits)
    end
    pd_out = PopData(geno_parse)
    !allow_monomorphic && dropmonomorphic!(pd_out, silent = silent)
    return pd_out
end

"""
    genepop(data::PopData; filename::String = "output.gen", digits::Int = 3, format::String = "vertical", miss::Int = 0)
Writes a `PopData` object to a Genepop-formatted file.
- `data`: the `PopData` object you wish to convert to a Genepop file
### keyword arguments
- `filename`: a `String` of the output filename
- `digits` : an `Integer` indicating how many digits to format each allele as (e.g. `(1, 2)` => `001002` for `digits = 3`)
- `format` : a `String` indicating whether loci should be formatted
    - vertically (`"v"` or `"vertical"`)
    - hortizontally (`"h"`, or `"horizontal"`)
    - Genepop Isolation-By-Distance (`"ibd"`) where each sample is a population with long/lat data prepended
- `miss` : an `Integer` for how you would like missing values written
    - `0` : As a genotype represented as a number of zeroes equal to `digits × ploidy` like `000000` (default)
    - `-9` : As a single value `-9`

## Example
```julia
cats = @nancycats;
fewer_cats = omit(cats, name = samplenames(cats)[1:10]);
genepop(fewer_cats, filename = "filtered_nancycats.gen", digits = 3, format = "h")
```
"""
function genepop(data::PopData; filename::String = "output.gen", digits::Int = 3, format::String = "vertical", miss::Int = 0)
    outfile = open(filename, "w")
    println(outfile, "generated from PopData by PopGen.jl")
    if format in ["h", "horizontal"]
        join(outfile, loci(data), ",")
    else
        join(outfile, loci(data), "\n")
    end
    print(outfile, "\n")
    if lowercase(format) != "ibd"
        pops = Vector{String}()
        for (keys, sample) in pairs(groupby(data.genodata, [:name, :population]))
            if keys.population ∉ pops
                push!(pops, keys.population)
                println(outfile, "POP")
            end
            samplename = keys.name
            sample_ploidy = convert(Int, data.sampleinfo.ploidy[data.sampleinfo.name .== samplename][1])
            format_geno = unphase.(sample.genotype, digits = digits, ploidy = sample_ploidy, miss = miss)
            join(outfile, vcat(samplename * ",", format_geno), "\t")
            print(outfile, "\n")
        end
    else
        for (keys, sample) in pairs(groupby(data.genodata, :name))
            samplename = keys.name
            sample_ploidy = convert(Int, data.sampleinfo.ploidy[data.sampleinfo.name .== samplename][1])
            println(outfile, "POP")
            long = data.sampleinfo[data.sampleinfo.name .== keys.name, :longitude][1]
            lat = data.sampleinfo[data.sampleinfo.name .== keys.name, :latitude][1]
            coords = join(string.([long, lat]), "\t") * ","
            format_geno = unphase.(sample.genotype, digits = digits, ploidy = sample_ploidy, miss = miss)
            join(outfile, vcat(coords, format_geno), "\t")
            print(outfile, "\n")
        end
    end
    close(outfile)
end


text:
###########
## asArr ##
###########

## asArr: try conversion to Array, implicitly dealing with NAs
##
## using DataFrames, columns may be of different type and hence are
## possibly not compatible with each other as columns of an array

function asArr{T}(da::DataArray{T, 1}, typ::Type=Any,
                  replaceNA=NA)
    ## convert DataArray to typ and replace NAs with replaceNA
    vals = convert(Array{typ, 1}, da.data)
    vals[da.na] = replaceNA
    return vals
end

function asArr(df::DataFrame, typ::Type=Any,
               replaceNA=NA)
    nObs, nAss = size(df)
    vals = Array(typ, nObs, nAss)
    ## copying into this array implicitly forces conversion to this
    ## type afterwards
    for ii=1:nAss
        if isa(df.columns[ii], DataArray) # deal with NAs
            if eltype(df.columns[ii]) == NAtype # column of NAs only
                vals[:, ii] = replaceNA
            else
                vals[:, ii] = df.columns[ii].data
                vals[df.columns[ii].na, ii] = replaceNA
            end
        else # no NAs in simple Array
            vals[:, ii] = df.columns[ii]
        end
    end
    return vals
end

function asArr(dfr::DataFrameRow, typ::Type=Any,
               replaceNA=NA)
    nObs = length(dfr)
    res = Array(typ, 1, nObs)
    for ii=1:nObs
        if isna(dfr[ii])
            res[1, ii] = replaceNA
        else
            res[1, ii] = dfr[ii]
        end
    end
    return res
end

function asArr(tn::AbstractTimedata, typ::Type=Any,
               replaceNA=NA)
    return asArr(tn.vals, typ, replaceNA)
end

## problem: columns with NAs only
## problem: NaN in boolean context will be true
## problem: DataArray{Any, 1} also exist: da.data then contains NAs

#############
## Timenum ##
#############

import Base.convert
## to Timedata: conversion upwards - always works
function convert(::Type{Timedata}, tn::Timenum)
    return Timedata(tn.vals, tn.idx)
end

## to Timematr: conversion downwards - fails for NAs
function convert(::Type{Timematr}, tn::Timenum, rmNA = false)
    if rmNA
        tm = narm(tn)
        tm = Timematr(tm.vals, tm.idx)
    else
        tm = Timematr(tn.vals, tn.idx)
    end
    return tm
end

##############
## Timematr ##
##############

## to Timedata: conversion upwards - always works
function convert(::Type{Timedata}, tm::Timematr)
    Timedata(tm.vals, tm.idx)
end

## to Timenum: conversion upwards - always works
function convert(::Type{Timenum}, tm::Timematr)
    Timenum(tm.vals, tm.idx)
end

##############
## Timedata ##
##############

## to Timenum: conversion downwards - fails for non-numeric values
function convert(::Type{Timenum}, td::Timedata)
    Timenum(td.vals, td.idx)
end

## to Timematr: conversion downwards - fails for NAs
function convert(::Type{Timematr}, td::Timedata)
    Timematr(td.vals, td.idx)
end

################
## DataFrames ##
################

## convert DataFrame with dates column (as String) to TimeData object

function convert(::Type{AbstractTimedata}, df::DataFrame)
    ## test if some column already is of type Date
    if any(eltypes(df) .== Date)
        ## take first occuring column
        dateCol = find(eltypes(df) .== Date)[1]
        idx = convert(Array, df[dateCol])
        delete!(df, [dateCol])

    else ## find column that contain dates as String

        # find columns that have been parsed as Strings by readtable
        col_to_test = Array(Symbol, 0)

        nCols = size(df, 2)
        for ii=1:nCols
            isa(df[1, ii], String)?
            push!(col_to_test, names(df)[ii]):
            nothing
        end

        # test each column's data to see if Datetime will parse it
        col_that_pass = Array(Symbol, 0)

        for colname in col_to_test
            d = match(r"[-|\s|\/|.]", df[1, colname])
            d !== nothing? (bar = split(df[1, colname], d.match)): (bar = [])
            if length(bar) == 3
                push!(col_that_pass, colname)
            end
        end

        # parse first column that passes the Datetime regex test
        idx = Date[Date(d) for d in df[col_that_pass[1]]] # without
        # Date it would fail chkIdx
        # in constructor

        delete!(df, [col_that_pass[1]])
    end

    ## try whether DataFrame fits subtypes
    try
        td = Timematr(df, idx)
        return td
    catch
        try
            td = Timenum(df, idx)
            return td
        catch
            td = Timedata(df, idx)
            return td
        end
    end

    return td
end

################
## TimeArrays ##
################

function convert(::Type{AbstractTimedata}, ta::TimeArray)

    namesAsSymbols = [DataFrames.identifier(nam) for nam in ta.colnames]

    df = composeDataFrame(ta.values, namesAsSymbols)
    idx = ta.timestamp

    try
        td = Timematr(df, idx)
        return td
    catch
        try
            td = Timenum(df, idx)
            return td
        catch
            td = Timedata(df, idx)
            return td
        end
    end
    td
end

function convert(::Type{TimeArray}, tn::AbstractTimenum)
    dats = idx(tn)
    nams = UTF8String[string(names(tn)[ii]) for ii = 1:size(tn, 2)]
    return TimeArray(dats, asArr(tn, Float64, NaN), nams)
end


text:
using Weave

let start_dir = pwd()
    cd(@__DIR__)
    outpath = normpath(joinpath(@__DIR__, "..", "docs", "src"))
    @show outpath
    if !isdirpath(outpath)
        mkpath(outpath)
    end

    weaveit(name) = weave(
        joinpath(@__DIR__, "$name.jmd"),
        out_path = joinpath(outpath, "$name.md"),
        doctype = "github",
    )

    weaveit.((
        "gettingstarted",
        "OByStoic",
        "fluoryield",
        "meanionizationpotential",
        "Au60Ag40unc",
        "eta",
        "bremsstrahlung",
    ))

    cd(start_dir)
end


text:
using GLVisualize, GeometryTypes, GLAbstraction, Colors

if !isdefined(:runtests)
    window = glscreen()
    timesignal = bounce(linspace(0.0, 1.0, 360))
end

_view(visualize(AABB{Float32}(Vec3f0(-2), Vec3f0(4)), :grid))

if !isdefined(:runtests)
    renderloop(window)
end


text:
using PTMCMC
using Test

@testset "PTMCMC.jl" begin
    # Write your tests here.
end


text:
type DataFrameView <: TreeView
  handle :: Ptr{GObject}
  store :: DataFrameStore
  function DataFrameView(df :: AbstractDataFrame)
    store = DataFrameStore(df)
    view = @TreeView(TreeModel(store))
    renderer = @CellRendererText()
    for i in eachindex(names(df))
      push!(view, @TreeViewColumn(string(names(df)[i]), renderer, Dict("text"=>i-1)))
    end
    Gtk.gobject_move_ref(new(view.handle, store), view)
  end
end

max_stars_repo_name: barcharcraz/eve-julia-gui

max_issues_repo_name: barcharcraz/eve-julia-gui

max_forks_repo_name: barcharcraz/eve-julia-gui