text:
"""
```
plot_prior_posterior(m, param_key; kwargs...)

plot_prior_posterior(m, param_keys = map(x -> x.key, m.parameters),
    include_fixed = false, verbose = :low, kwargs...)
```

Plot prior distribution and histogram of posterior draws for specified
`param_key` or `param_keys`.

### Inputs

- `m::AbstractDSGEModel`
- `param_key::Symbol` or `param_keys::Vector{Symbol}`: e.g. `:α` or `[:α, :ζ_p]`

### Keyword Arguments

- `verbose::Symbol`

See `?priorpost` for additional keyword arguments, all of which can be passed
into `plot_prior_posterior`.

**Method 2 only:**

- `include_fixed::Bool`: whether to plot fixed parameters
"""
function plot_prior_posterior(m::AbstractDSGEModel, param_key::Symbol; kwargs...)

    plots = plot_prior_posterior(m, [param_key]; include_fixed = true, kwargs...)
    return plots[param_key]
end

function plot_prior_posterior(m::AbstractDSGEModel,
                              param_keys::Vector{Symbol} = map(x -> x.key, m.parameters);
                              include_fixed::Bool = false,
                              verbose::Symbol = :low,
                              kwargs...)
    # Load parameter draws from Metropolis-Hastings
    posterior_draws = load_draws(m, :full)

    # Loop through parameters
    plots = OrderedDict{Symbol, Plots.Plot}()
    for key in param_keys
        i = findfirst(x -> x.key == key, m.parameters)
        if i > 0
            param = m.parameters[i]
            posterior = posterior_draws[:, i]

            # Skip fixed parameter if desired
            if !include_fixed && param.fixed
                continue
            end

            # Call recipe
            plots[key] = priorpost(param, posterior; kwargs...)

            # Save plot
            output_file = figurespath(m, "estimate", "prior_posterior_" * detexify(string(param.key)) * ".pdf")
            save_plot(plots[key], output_file, verbose = verbose)

        else
            error("Parameter not found: " * string(param_key))
        end
    end
    return plots
end

@userplot PriorPost

"""
```
priorpost(param, posterior_draws)
```

User recipe called by `plot_prior_posterior`.

### Inputs

- `param::Parameter`
- `posterior_draws::AbstractVector`

### Keyword Arguments

- `prior_color`
- `posterior_color`

Additionally, all Plots attributes (see docs.juliaplots.org/latest/attributes)
are supported as keyword arguments.
"""
priorpost

@recipe function f(pp::PriorPost;
                   prior_color = :red,
                   posterior_color = :blue)
    # Error checking
    if length(pp.args) != 2 || !(typeof(pp.args[1]) <: Parameter) || !(typeof(pp.args[2]) <: AbstractVector)
        error("priorpost must be given a Parameter and an AbstractVector. Got $(typeof(pp.args))")
    end

    param, posterior_draws = pp.args

    title  --> param.tex_label
    legend --> :bottomright

    # Posterior
    @series begin
        seriestype := :histogram
        normalize  := true
        label      := param.fixed ? "Posterior: " * describe_prior(param) : "Posterior"
        color      := posterior_color

        posterior_draws
    end

    # Prior
    @series begin
        label      := "Prior: " * describe_prior(param)
        color      := prior_color
        linewidth --> 4

        if param.fixed
            [minimum(posterior_draws), maximum(posterior_draws)], [1, 1]
        elseif !param.fixed && !isnull(param.prior)
            prior = get(param.prior)
            typeof(prior), prior
        else
            error("Parameter must either be fixed or have a non-null prior")
        end
    end
end


text:
module FunRefsTest

using Comm, GIDs, FunRefs
using Base.Test

# Test basic operations
function test_basic(T::Type)
    i = 1
    ri = FunRef{T}(i)
    @test ri[] == i
    rn = FunRef{Union{}}()
    global DONE = false
    rexec(mod1(2, Comm.nprocs())) do
        remote1(i, ri)
    end
    while !DONE yield() end
end

# Test remote garbage collection
function test_remote(do_gc::Bool)
    # Use a nested function to ensure that ri is garbage collected
    function inner()
        i = 1
        ri = FunRef(i)
        maxcount = 100 # 1000
        global COUNT = 0
        manyrefs(i, ri, maxcount, do_gc)
        while COUNT<maxcount yield() end
        sleep(1)
        @test COUNT==maxcount
    end
    inner()
    global COUNT = 0
    rexec_everywhere(check_refs)
    while COUNT<Comm.nprocs() yield() end
    sleep(1)
end

# Test future behaviour
function test_future()
    i = 1
    ri = FunRef(i)
    @test isready(ri)
    @test islocal(ri)
    @test ri[] == i
    ri = FunRef{Int}()
    @test !isready(ri)
    ri[] = i
    @test isready(ri)
    @test islocal(ri)
    @test ri[] == i
    r0 = FunRef{Void}()
    @test !isready(r0)
    global FUTURE_SET = false
    global FUTURE_DONE = false
    rexec(mod1(2, Comm.nprocs())) do
        future_continued(ri, r0)
    end
    while !FUTURE_SET yield() end
    r0[] = nothing
    while !FUTURE_DONE yield() end
end

function remote1(i::Int, ri::FunRef)
    @test proc(ri) == 1
    @test islocal(ri) == (Comm.myproc()==1)
    rexec(proc(ri)) do
        remote2(i, ri)
    end
end

function remote2(i::Int, ri::FunRef)
    @test islocal(ri)
    @test ri[] == i
    global DONE = true
end

function manyrefs(i::Int, ri::FunRef, count::Int, do_gc::Bool)
    @assert count>0
    if do_gc gc() end
    @test !islocal(ri) || ri[] == i
    if count == 1
        rexec(1) do
            inc()
        end
        return
    end
    proc1 = rand(1:Comm.nprocs())
    proc2 = rand(1:Comm.nprocs())
    proc3 = rand(1:Comm.nprocs())
    counta = rand(0:count)
    countb = rand(0:count)
    count1 = min(counta, countb)
    count2 = max(counta, countb) - count1
    count3 = count - (count1 + count2)
    @assert count1 + count2 + count3 == count
    if count1>0
        rexec(proc1) do
            manyrefs(i, ri, count1, do_gc)
        end
    end
    if count2>0
        rexec(proc2) do
            manyrefs(i, ri, count2, do_gc)
        end
    end
    j = i+1
    rj = FunRef(j)
    if count3>0
        rexec(proc3) do
            manyrefs(j, rj, count3, do_gc)
        end
    end
end

function inc()
    global COUNT += 1
end

function check_refs()
    gc()
    sleep(1)
    @test isempty(GIDs.mgr())
    rexec(1) do
        inc()
    end
end

function future_continued(ri::FunRef, r0::FunRef)
    @test isready(ri)
    @test islocal(ri) == (Comm.nprocs()==1)
    rj = typeof(ri)()
    @test !isready(rj)
    set_from_ref!(rj, ri)
    @test isready(rj)
    @test islocal(rj) == (Comm.nprocs()==1)
    @test !isready(r0)
    rexec(1) do
        future_set()
    end
    wait(r0)
    rexec(1) do
        future_done()
    end
end

function future_set()
    global FUTURE_SET = true
end

function future_done()
    global FUTURE_DONE = true
end

function main()
    test_basic(Any)
    test_basic(Int)
    test_remote(true)
    test_remote(false)
    test_future()
end

end

ext: jl


text:
module ModelCompare

using MathOptFormat
using MathOptInterface

export compare_models, fill_constraint_names, remove_intervals

const MOI = MathOptInterface
const MOIU = MathOptInterface.Utilities

const VI = MOI.VariableIndex
const CI = MOI.ConstraintIndex

include("sorting.jl")
include("modify.jl")
include("utils.jl")

function compare_models(model1::MOI.ModelLike, model2::MOI.ModelLike, tol = 1e-8)
    io = STDIO

    # TODO: give test-friendly feedback instead of errors?
    extra_1_model = MathOptFormat.MOF.Model()
    extra_2_model = MathOptFormat.MOF.Model()

    model1_vnames = [MOI.get(model, MOI.VariableName(), index)
                     for index in MOI.get(model, MOI.ListOfVariableIndices())]
    model1_cnames = String[]
    for (F,S) in MOI.get(model, MOI.ListOfConstraints())
        for index in MOI.get(model, MOI.ListOfConstraintIndices{F,S}())
            push!(model1_cnames, MOI.get(model, MOI.ConstraintName(), index))
        end
    end
    # model1_cnames = vec([MOI.get(model, MOI.ConstraintName(), index)
    #                      for (F,S) in MOI.get(model, MOI.ListOfConstraints()),
    #                          index in MOI.get(model, MOI.ListOfConstraintIndices{F,S}())])

    extra_1_vnames, extra_2_vnames, same_vnames = compare_variablenames(model2, model1_vnames)
    println("Variables in both: $(length(same_vnames))")
    println("Variables only in left: $(length(extra_1_vnames))")
    println("Variables only in right: $(length(extra_2_vnames))")
    extra_1_cnames, extra_2_cnames, same_cnames = compare_constraintnames(model2, model1_cnames)
    println("Constraints in both: $(length(same_cnames))")
    println("Constraints only in left: $(length(extra_1_cnames))")
    println("Constraints only in right: $(length(extra_2_cnames))")

    variablemap_2to1 = Dict{VI,VI}()
    variablemap_1to2 = Dict{VI,VI}()
    for vname in same_vnames
        index1 = MOI.get(model1, VI, vname)
        index2 = MOI.get(model2, VI, vname)
        variablemap_2to1[index2] = index1
        variablemap_1to2[index1] = index2
    end

    variablemap_1_old_to_new = Dict{VI,VI}()
    extra_1_vi = Set{VI}()
    for vname in extra_1_vnames
        index1 = MOI.get(model1, VI, vname)
        # add variables
        x = MOI.add_variable(extra_1_model)
        MOI.set(extra_1_model, MOI.VariableName(), x, vname)
        variablemap_1_old_to_new[index1] = x
        push!(extra_1_vi, index1)
    end

    variablemap_2_old_to_new = Dict{VI,VI}()
    extra_2_vi = Set{VI}()
    for vname in extra_2_vnames
        index2 = MOI.get(model2, VI, vname)
        # add variables
        x = MOI.add_variable(extra_2_model)
        MOI.set(extra_2_model, MOI.VariableName(), x, vname)
        variablemap_2_old_to_new[index2] = x
        push!(extra_2_vi, index2)
    end

    for cname in same_cnames
        index1 = MOI.get(model1, CI, cname)
        index2 = MOI.get(model2, CI, cname)
        f1 = MOI.get(model1, MOI.ConstraintFunction(), index1)
        f2 = MOI.get(model2, MOI.ConstraintFunction(), index2)

        f1_type = typeof(f1)
        f2_type = typeof(f2)
        @assert f1_type == f2_type

        s1 = MOI.get(model1, MOI.ConstraintSet(), index1)
        s2 = MOI.get(model2, MOI.ConstraintSet(), index2)
        s_12 = set_difference(s1,s2)
        s_21 = set_difference(s2,s1)
        is_approx_zero_s_12 = is_approx_zero(s_12, tol)

        if f1_type <: MOI.ScalarAffineFunction#, MOI.ScalarQuadraticFunction]

            f1_extra = MOI.zero(f1_type)
            f2_extra = MOI.zero(f2_type)
            unsafe_move_variables(extra_1_vi, f1, f1_extra)
            unsafe_move_variables(extra_2_vi, f2, f2_extra)

            f_12 = MOIU.canonical(f1 - MOIU.mapvariables(variablemap_2to1, f2))
            f_21 = MOI.zero(f1_type)

            if !MOIU.isapprox_zero(f_12, tol)
                f_21 = MOIU.canonical(f2 - MOIU.mapvariables(variablemap_1to2, f1))
                pass_variables(model1, extra_1_model, f_12, variablemap_1_old_to_new, tol)
                pass_variables(model2, extra_2_model, f_21, variablemap_2_old_to_new, tol)
                f_12 = MOIU.isapprox_zero(f1_extra, tol) ? f_12 : MOIU.canonical(f_12 + f1_extra)
                f_21 = MOIU.isapprox_zero(f2_extra, tol) ? f_21 : MOIU.canonical(f_21 + f2_extra)
                c1 = MOI.add_constraint(extra_1_model, MOIU.mapvariables(variablemap_1_old_to_new, f_12), s_12)
                c2 = MOI.add_constraint(extra_2_model, MOIU.mapvariables(variablemap_2_old_to_new, f_21), s_21)
                MOI.set(extra_1_model, MOI.ConstraintName(), c1, cname)
                MOI.set(extra_2_model, MOI.ConstraintName(), c2, cname)
            else
                if !MOIU.isapprox_zero(f1_extra, tol)
                    c1 = MOI.add_constraint(extra_1_model, MOIU.mapvariables(variablemap_1_old_to_new, f1_extra), s_12)
                    MOI.set(extra_1_model, MOI.ConstraintName(), c1, cname)
                end
                if !MOIU.isapprox_zero(f2_extra, tol)
                    c2 = MOI.add_constraint(extra_2_model, MOIU.mapvariables(variablemap_2_old_to_new, f2_extra), s_21)
                    MOI.set(extra_2_model, MOI.ConstraintName(), c2, cname)
                end
            end
        elseif f1_type <: MOI.SingleVariable
            if f1.variable != MOIU.mapvariables(variablemap_2to1, f2).variable || !is_approx_zero_s_12
                pass_variables(model1, extra_1_model, f1, variablemap_1_old_to_new, tol)
                pass_variables(model2, extra_2_model, f2, variablemap_2_old_to_new, tol)
                c1 = MOI.add_constraint(extra_1_model, f1_extra, s_12)
                c2 = MOI.add_constraint(extra_2_model, f2_extra, s_21)
                MOI.set(extra_1_model, MOI.ConstraintName(), c1, cname)
                MOI.set(extra_2_model, MOI.ConstraintName(), c2, cname)
            end
        elseif f1_type == MOI.VectorOfVariables
            error("Vector os variables comparisson not supported")
        else
            error("Invalid function type")
        end
    end
    for cname in extra_1_cnames
        index = MOI.get(model1, CI, cname)
        s = MOI.get(model1, MOI.ConstraintSet(), index)
        f = MOI.get(model1, MOI.ConstraintFunction(), index)
        pass_variables(model1, extra_1_model, f, variablemap_1_old_to_new, tol)
        c = MOI.add_constraint(extra_1_model, MOIU.mapvariables(variablemap_1_old_to_new, f), s)
        MOI.set(extra_2_model, MOI.ConstraintName(), c, cname)
    end
    for cname in extra_2_cnames
        index = MOI.get(model2, CI, cname)
        s = MOI.get(model2, MOI.ConstraintSet(), index)
        f = MOI.get(model2, MOI.ConstraintFunction(), index)
        pass_variables(model2, extra_2_model, f, variablemap_2_old_to_new, tol)
        c = MOI.add_constraint(extra_2_model, MOIU.mapvariables(variablemap_2_old_to_new, f), s)
        MOI.set(extra_2_model, MOI.ConstraintName(), c, cname)
    end

    sort_model(extra_1_model)
    sort_model(extra_2_model)

    if true
        lp_1 = MathOptFormat.LP.Model()
        MOI.copy_to(lp_1, extra_1_model)
        lp_2 = MathOptFormat.LP.Model()
        MOI.copy_to(lp_2, extra_2_model)
        MOI.write_to_file(lp_1, "A.lp")
        MOI.write_to_file(lp_2, "B.lp")
    end

    if true
        mps_1 = MathOptFormat.MPS.Model()
        MOI.copy_to(mps_1, extra_1_model)
        mps_2 = MathOptFormat.MPS.Model()
        MOI.copy_to(mps_2, extra_2_model)
        MOI.write_to_file(mps_1, "A.mps")
        MOI.write_to_file(mps_2, "B.mps")
    end

    MOI.write_to_file(extra_1_model, "A.mof")
    MOI.write_to_file(extra_2_model, "B.mof")

end

function compare_variablenames(model, variablenames)
    missing_names = String[] # TODO use sets?
    extra_names = String[]
    same_names = String[]
    seen_name = Dict(name => false for name in variablenames)
    for index in MOI.get(model, MOI.ListOfVariableIndices())
        vname = MOI.get(model, MOI.VariableName(), index)
        if !haskey(seen_name, vname)
            push!(extra_names, vname)
            # error("Variable with name $vname present in model but not expected list of variable names.")
        elseif seen_name[vname]
            println("Variable with name $vname present twice in model (shouldn't happen!)")
        else
            seen_name[vname] = true
            push!(same_names, vname)
        end
    end
    for (vname,seen) in seen_name
        if !seen
            push!(missing_names, vname)
            # error("Did not find variable with name $vname in intance.")
        end
    end
    return missing_names, extra_names, same_names
end
function compare_constraintnames(model, constraintnames)
    missing_names = String[]
    extra_names = String[]
    same_names = String[]
    seen_name = Dict(name => false for name in constraintnames)
    for (F,S) in MOI.get(model, MOI.ListOfConstraints())
        for index in MOI.get(model, MOI.ListOfConstraintIndices{F,S}())
            cname = MOI.get(model, MOI.ConstraintName(), index)
            if !haskey(seen_name, cname)
                push!(extra_names, cname)
                # error("Constraint with name $cname present in model but not expected list of constraint names.")
            elseif seen_name[cname]
                println("Constraint with name $cname present twice in model (shouldn't happen!)")
            else
                seen_name[cname] = true
                push!(same_names, cname)
            end
        end
    end
    for (cname,seen) in seen_name
        if !seen
            push!(missing_names, cname)
            # error("Did not find constraint with name $cname in intance.")
        end
    end
    return missing_names, extra_names, same_names
end

end


text:
function nlsolve!(F::OnceDiffed, x, approach::Tuple=(Newton(), Static(1)); maxiter=200, f_abstol=1e-8, f_reltol=1e-12)
    scheme, linesearch = approach
    xp, d, Fx, Jx = copy(x), copy(x), copy(x), x*x'

    Fx, Jx = F(x, Fx, Jx)
    normFx = norm(Fx, 2)
    T = typeof(normFx)
    stoptol = T(f_reltol)*normFx + T(f_abstol)
    if normFx < stoptol
        return x, Fx, 0
    end
    iter = 1
    while iter ≤ maxiter
        # Update the point of evaluation
        d .= -(Jx\Fx)

        merit = MeritObjective(F, Fx, Jx, d)
        # φ = LineObjective!(F, ∇fz, z, x, d, fx, dot(∇fx, d))
        φ = LineObjective(merit, x, d, (normFx^2)/2, -normFx^2)
        # Perform line search along d
        α, f_α, ls_success = find_steplength(linesearch, φ, T(1.0))

        x .= x .+ α.*d
        Fx, Jx = F(x, Fx, Jx)
        normFx = norm(Fx, 2)
        if normFx < stoptol
            break
        end
        iter += 1
    end
    return x, Fx, iter
end


text:
using PolynomialZonotopes, Base.Test

# =====================================
# Testing the polynomial zonotope type
# =====================================
@time @testset "LazySets.PolynomialZonotope" begin include("unit_PolynomialZonotope.jl") end


text:
using ApproxFun, Test

d=Interval(-1000.,5.)
x=Fun(identity,d)

A=Derivative(d)^2-x

@time A[1:12500,1:12505]  # 0.002191
@time A.ops[1][1:12500,1:12505]  # 0.000253
@time A.ops[2][1:12500,1:12505]  # 0.002507

# T=TimesOperator([A.ops[2].ops[1].op.ops...,A.ops[2].ops[2]])
# @time T[1:12500,1:12505]
# bandwidth(A.ops[2].ops[2],2)
@time a=A.ops[2].ops[1][1:12500,1:12504]  # 0.000770
@time b=A.ops[2].ops[2][1:12504,1:12505]  # 0.000117

@time A.ops[2].ops[1].op.ops[1][1:12500,1:12504]  # 0.000177
@time A.ops[2].ops[1].op.ops[2][1:12500,1:12504]  # 0.000167

@time A.ops[2].ops[2][1:12500,1:12504]  # 0.000091

#
# @profile for k=1:100
#     A.ops[2].ops[2][1:12500,1:12504]
# end
#
# Profile.print()
@time a*b  #0.0014


text:
### A Pluto.jl notebook ###
# v0.12.21

using Markdown
using InteractiveUtils

# ╔═╡ 62405c38-8dfb-11eb-0dfe-b7545eb40047
using ImageDraw, ColorTypes, BenchmarkTools,ColorVectorSpace

# ╔═╡ 0a8bbcd6-8df2-11eb-1010-9fa5c200e91b
 begin
vert=CartesianIndex{2}[]
	push!(vert, CartesianIndex(2,2))
	push!(vert, CartesianIndex(3,2))
    push!(vert, CartesianIndex(4,2))
	push!(vert, CartesianIndex(5,2))
    push!(vert, CartesianIndex(6,3))
	push!(vert, CartesianIndex(7,4))
	push!(vert, CartesianIndex(8,4))
	push!(vert, CartesianIndex(9,3))
	push!(vert, CartesianIndex(10,2))
	push!(vert, CartesianIndex(11,3))
	push!(vert, CartesianIndex(12,4))
	push!(vert, CartesianIndex(13,5))
	push!(vert, CartesianIndex(12,6))
	push!(vert, CartesianIndex(11,7))
	push!(vert, CartesianIndex(10,6))#u
	push!(vert, CartesianIndex(9,7))
	push!(vert, CartesianIndex(8,8))
	push!(vert, CartesianIndex(8,9))
	push!(vert, CartesianIndex(8,10))
	push!(vert, CartesianIndex(8,11))
	push!(vert, CartesianIndex(9,11))
	push!(vert, CartesianIndex(10,11))
	push!(vert, CartesianIndex(11,11))
	push!(vert, CartesianIndex(12,11))
	push!(vert, CartesianIndex(13,11))
	push!(vert, CartesianIndex(13,12))
	push!(vert, CartesianIndex(13,13))
	push!(vert, CartesianIndex(12,13))
	push!(vert, CartesianIndex(11,13))
	push!(vert, CartesianIndex(10,13))
	push!(vert, CartesianIndex(9,13))
	push!(vert, CartesianIndex(8,13))
	push!(vert, CartesianIndex(7,13))
	push!(vert, CartesianIndex(6,13))
	push!(vert, CartesianIndex(6,12))
	push!(vert, CartesianIndex(6,11))#u
	push!(vert, CartesianIndex(5,11))
	push!(vert, CartesianIndex(4,11))
	push!(vert, CartesianIndex(3,11))
	push!(vert, CartesianIndex(2,11))#u
	push!(vert, CartesianIndex(3,10))
	push!(vert, CartesianIndex(4,9))
	push!(vert, CartesianIndex(5,8))
	push!(vert, CartesianIndex(4,7))
	push!(vert, CartesianIndex(4,6))
	push!(vert, CartesianIndex(5,5))
	push!(vert, CartesianIndex(4,4))
    push!(vert, CartesianIndex(3,3))
	push!(vert, CartesianIndex(2,2))
end

# ╔═╡ 76b68af4-8dfb-11eb-11df-3f5ecff30aef
img = draw(zeros(Gray{Bool},14,14), Polygon(vert));

# ╔═╡ c037b994-8e39-11eb-308c-db84d51ff35c
struct edgetabletuple
	initial::CartesianIndex
	final::CartesianIndex
end

# ╔═╡ e728998c-8e03-11eb-383c-63adbbb0a9a7
begin
function createedgetable()
edgetable= []
push!(edgetable, edgetabletuple(CartesianIndex(2,2),CartesianIndex(5,2)))
push!(edgetable, edgetabletuple(CartesianIndex(5,2),CartesianIndex(7,4)))
push!(edgetable, edgetabletuple(CartesianIndex(7,4),CartesianIndex(8,4)))
push!(edgetable, edgetabletuple(CartesianIndex(8,4),CartesianIndex(10,2)))
push!(edgetable, edgetabletuple(CartesianIndex(10,2),CartesianIndex(13,5)))
push!(edgetable, edgetabletuple(CartesianIndex(13,5),CartesianIndex(11,7)))
push!(edgetable, edgetabletuple(CartesianIndex(11,7),CartesianIndex(10,6)))
push!(edgetable, edgetabletuple(CartesianIndex(10,6),CartesianIndex(8,8)))
push!(edgetable, edgetabletuple(CartesianIndex(8,8),CartesianIndex(8,11)))
push!(edgetable, edgetabletuple(CartesianIndex(8,11),CartesianIndex(13,11)))
push!(edgetable, edgetabletuple(CartesianIndex(13,11),CartesianIndex(13,13)))
push!(edgetable, edgetabletuple(CartesianIndex(13,13),CartesianIndex(6,13)))
push!(edgetable, edgetabletuple(CartesianIndex(6,13),CartesianIndex(6,11)))
push!(edgetable, edgetabletuple(CartesianIndex(6,11),CartesianIndex(2,11)))
push!(edgetable, edgetabletuple(CartesianIndex(2,11),CartesianIndex(5,8)))
push!(edgetable, edgetabletuple(CartesianIndex(5,8),CartesianIndex(4,8)))
push!(edgetable, edgetabletuple(CartesianIndex(4,8),CartesianIndex(4,7)))
push!(edgetable, edgetabletuple(CartesianIndex(4,7),CartesianIndex(4,6)))
push!(edgetable, edgetabletuple(CartesianIndex(4,6),CartesianIndex(5,5)))
push!(edgetable, edgetabletuple(CartesianIndex(5,5),CartesianIndex(2,2)))
return edgetable
end
end

# ╔═╡ cce6a1e6-8e39-11eb-024f-7bd13f877a7a
#find ymin and ymax
function yminmax(vert)
	ymax = vert[1][1];
	ymin = vert[1][1];
	for i in 1:length(vert)
	   if(vert[i][2] > ymax)
		 ymax = vert[i][2]
	   elseif(vert[i][2] < ymin)
		 ymax = vert[i][2]
	   end
	end
	return ymin, ymax
end

# ╔═╡ f81db174-8e39-11eb-03d0-addb3b5b039c
function findintersections(edgetable, yvalue)
	points =[]
    for i in 1:length(edgetable)
		if(edgetable[i].final[2] > yvalue && edgetable[i].initial[2] > yvalue)
			continue
		end
		if(edgetable[i].final[2] <= yvalue && edgetable[i].initial[2] <= yvalue)
			continue
		end
		x=1
		deltay = edgetable[i].final[1]-edgetable[i].initial[1]
		deltax = edgetable[i].final[2]-edgetable[i].initial[2]
		alpha = deltay/deltax
		constant = edgetable[i].initial[1]
		x = alpha * (yvalue - edgetable[i].initial[2]) + constant
		if (x == -Inf || isnan(x) || x == Inf) continue end
		push!(points,CartesianIndex(ceil(Int,x),yvalue))
	end
	return points
end

# ╔═╡ fc3e886a-8dfc-11eb-28e2-a56245654c9a
function timetocolor(points, fill_color::T)where T<:Colorant
   for i in 1:2:length(points)-1
	  draw!(img, LineSegment(points[i], points[i+1]),fill_color)
   end
end

# ╔═╡ da3d93c2-8e39-11eb-0850-93179bd9e20c
function scanline(vert, fill_color::T) where T<:Colorant
	ymin, ymax = yminmax(vert)
	edgetable = createedgetable()
	pyramid = []
	for i in ymin:ymax
		points = findintersections(edgetable, i)
		points = sort!(points, by = x -> x[1])
		timetocolor(points, fill_color)
		push!(pyramid,img[:,:])
	end
	pyramid
end

# ╔═╡ 7f369560-8e01-11eb-214b-83d8210ca559
#with_terminal() do
begin
	fill_color = Gray{Bool}(1.0)
    pyramid = scanline(vert, fill_color)
end

# ╔═╡ Cell order:
# ╠═62405c38-8dfb-11eb-0dfe-b7545eb40047
# ╠═0a8bbcd6-8df2-11eb-1010-9fa5c200e91b
# ╠═76b68af4-8dfb-11eb-11df-3f5ecff30aef
# ╠═e728998c-8e03-11eb-383c-63adbbb0a9a7
# ╠═c037b994-8e39-11eb-308c-db84d51ff35c
# ╠═cce6a1e6-8e39-11eb-024f-7bd13f877a7a
# ╠═f81db174-8e39-11eb-03d0-addb3b5b039c
# ╠═fc3e886a-8dfc-11eb-28e2-a56245654c9a
# ╠═da3d93c2-8e39-11eb-0850-93179bd9e20c
# ╠═7f369560-8e01-11eb-214b-83d8210ca559

max_stars_repo_name: ashwani-rathee/Demos

max_issues_repo_name: ashwani-rathee/Demos

max_forks_repo_name: ashwani-rathee/Demos


text:
using Test, ASE, NeighbourLists, JuLIP
using LinearAlgebra: norm

h1("Misc ASE Tests")

# ======================================================================

h3("Check that the cubic unit cell for Al is reproduced correctly")
a0 = 2.025
p = [ [0.0;0.0;0.0] [0.0;a0;a0] [a0;0.0;a0] [a0;a0;0.0] ]
at = bulk("Al", cubic=true)
println(@test (positions(at) |> mat) == p)

# ======================================================================

h2("Check neighbourlist without periodicity")
# TODO: implement a test with periodicity as well???
at = bulk("Al", cubic=true) * 3
set_pbc!(at, (false,false,false))
rcut = 1.7 * a0
nlist = neighbourlist(at, rcut)

@info "   ... assemble neighbour list ..."
# create a neighbourlist via a naive double-loop
simple = zeros(length(at), length(at))
X = positions(at)
for n = 2:length(at), m = 1:n-1
   if norm(X[m]-X[n]) <= rcut
      simple[n,m] = simple[m,n] = 1
   end
end

@info "   ... check the bond-iterator ... "
pass_bonds_test = true
for (i,j,r,R) in pairs(nlist)
   if !( (simple[i,j] == 1) && (abs(norm(X[i]-X[j]) - r) < 1e-12) &&
         (norm(X[j]-X[i] - R) < 1e-12)  )
      global pass_bonds_test = false
      break
   end
   # switch the flag
   simple[i,j] = -1
end

h3("check that no false bonds have been found")
println(@test pass_bonds_test)
h3("check that all bonds have been found")
println(@test maximum(simple) == 0)

# revert to original
simple *= -1
h3("   ... check the site iterator ... ")
pass_site_test = true
for (i,j,r,R) in sites(nlist)
   for n = 1:length(j)
      if simple[i,j[n]] != 1
         global pass_site_test = false
         break
      end
      simple[i,j[n]] = -1
   end
   if maximum(simple[i,:]) != 0
      global pass_site_test = false
      break
   end
end
println(@test pass_site_test)
println(@test maximum(simple) == 0)

# ======================================================================

# h3("Checking `***_data`, `***_array`, `***_info`, `***_transient`")
#
# at = bulk("Cu") * 3
# N = length(at)
# # set an array and test that it is read back correctly
# z = rand(N)
# set_array!(at, "z", z)
# @test get_array(at, "z") == z
# # set some info and test reading
# i = "some info"
# set_info!(at, "i", i)
# @test get_info(at, "i") == i
# # ***_transient should be tested automatically via the calculators.
# # test the has_***
# @test has_array(at, "z")
# @test has_info(at, "i")
# @test !has_array(at, "i")
# @test !has_info(at, "z")

# println("Checking momenta and velocites")
# p = rand(3,N) |> vecs
# set_momenta!(at, p)
# @test isapprox(momenta(at) |> mat, p |> mat)
# v = p ./ masses(at)
# @test isapprox(velocities(at) |> mat, v |> mat)
# set_velocities!(at, v/2.0)
# @test isapprox(velocities(at) |> mat, v/2.0 |> mat)
# @test isapprox(momenta(at) |> mat, p/2.0 |> mat)
#
# println("Checking chemical_symbols")
# @test chemical_symbols(at) == ["Cu" for i in 1:length(at)]

# ======================================================================

# h3("testing static_neighbourlist")
# nlist = static_neighbourlist(at, rnn("Cu") * 1.1)
# rattle!(at, 0.05)
# nlist1 = static_neighbourlist(at, rnn("Cu") * 1.1)
# @test nlist === nlist1

# ======================================================================

h3("testing read and write")
fname = "test.xyz"
at = bulk("Cu") * 2
write_xyz(fname, at)
at2 = read_xyz(fname)
println(@test positions(at) == positions(at2))
rm(fname)

# ---------
h3("Conversion between JuLIP.Atoms and ASEAtoms")
at1 = bulk(:Si)          # JuLIP.Atoms
at2 = ASEAtoms(at1)      # ASE.ASEAtoms
@show cell(at2)
at3 = Atoms(at2)         # JuLIP.Atoms
println(@test at1 == at3)

# ---------
h3("Conversion between JuLIP.Atoms and ASEAtoms, Multi-species")
at1 = bulk(:Al) * 3         # JuLIP.Atoms
at1.Z[1:(length(at1)÷2)] .= atomic_number(:Ti)
at1.M[:] .= 1.0 .+ rand(length(at1))
at2 = ASEAtoms(at1)      # ASE.ASEAtoms
@show cell(at2)
at3 = Atoms(at2)         # JuLIP.Atoms
println(@test at1 == at3)


text:
module CorrelationCurveFit

using LsqFit

export HyperModel, GeomModel, ExpModel, FitHyperModel, FitGeomModel, FitExpModel

HyperModel(xdata, p) = 1 ./ (p[1] + p[2] .* xdata)
function FitHyperModel(xData, yData)
	p0 = linreg(xData, 1 ./ yData)
	fit = curve_fit(HyperModel, xData, yData, p0)
	return fit
end

GeomModel(xdata, p) = 1 ./ (p[1] + p[2] * (1- (1 - p[3]).^xdata))
function FitGeomModel(xData, yData)
	hyperFit = FitHyperModel(xData, yData)
	b0 = hyperFit.param[1] + hyperFit.param[2] * 100
	b1 = hyperFit.param[2] * 100
	b2 = 0.01
	p0 = [b0, b1, b2]
	fit = curve_fit(GeomModel, xData, yData, p0)
end

ExpModel(xdata, p) = 1 ./ (p[1] + p[2] - p[2] .* exp(-xdata ./ p[3]))
function FitExpModel(xData, yData)
	hyperFit = FitHyperModel(xData, yData)
	b0 = hyperFit.param[1] + hyperFit.param[2] * 100
	b1 = hyperFit.param[2] * 100
	b2 = 100
	p0 = [b0, b1, b2]
	fit = curve_fit(ExpModel, xData, yData, p0)
	return fit
end

end # module


text:
# verbose flag
const VERBOSE = Ref(false)

"""
    set_verbose(::Bool)

For debugging - set flag to print macro expansions
"""
function set_verbose(verbose::Bool)
    VERBOSE[] = verbose
end


text:
using Images, AxisArrays
using Base: OneTo
using Unitful: Time, Length, s

# TODO: Make this be generic over type of SIM acquisition like so
# SIMAcquisition{T} where T <: AbstractSIMType <19-05-21, kunzaatko> #
struct SIMAcquisition <: AbstractAcquisition # {{{
    "Acquired reconstructed images"
    reconstructedimages::AxisArray{<:AbstractGray,3}
    "Acquired unreconstructed images stored with there acquisition times and orientations"
    unreconstructedimages::AxisArray{<:AbstractGray,3}
    "Numerical aperature"
    NA::Real
    "Magnification"
    M::Unsigned
    "Fluorophore marker"
    marker::FluorophoreMarker
    "Pixel size"
    pixelsize::Length
end # }}}

function SIMAcquisition( # {{{
    reconstructedimages::AbstractArray{<:Color,3},
    unreconstructedimages::AbstractArray{<:Color,3},
    acquisitiontimes::AbstractVector{<:Time},
    NA,
    M,
    marker::FluorophoreMarker,
    pixelsize::Length,
)
    if !(eltype(reconstructedimages) <: AbstractGray)
        reconstructedimages = Gray.(reconstructedimages)
    end

    if !(eltype(unreconstructedimages) <: AbstractGray)
        unreconstructedimages = Gray.(unreconstructedimages)
    end

    # FIX: Add check for acquisition times having the same size <19-05-21, kunzaatko> #
    # TODO: This should be generic over AbstractSIMType <19-05-21, kunzaatko> #
    # TODO: x and y axes should be indexed by length instead of index <19-05-21, kunzaatko> #
    unreconstructedimages_axis = AxisArray(
        unreconstructedimages,
        Axis{:x}(OneTo(size(unreconstructedimages, 1))),
        Axis{:y}(OneTo(size(unreconstructedimages, 2))),
        Axis{:time}(acquisitiontimes),
    )

    # TODO: This should make an average for the times of acquisition for the unreconstructed sequence that soruced this <19-05-21, kunzaatko> #
    # TODO: x and y axes should be indexed by length instead of index <19-05-21, kunzaatko> #
    reconstructedimages_axis = AxisArray(reconstructedimages, :x, :y, :time)

    return SIMAcquisition(
        reconstructedimages_axis, unreconstructedimages_axis, NA, M, marker, pixelsize
    )
end # }}}

function SIMAcquisition( # {{{
    reconstructedimages::AbstractArray{<:Color,3},
    unreconstructedimages::AbstractArray{<:Color,3},
    acquisitiontimestep::Time,
    NA,
    M,
    marker::FluorophoreMarker,
    pixelsize::Length,
)
    return SIMAcquisition(
        reconstructedimages,
        unreconstructedimages,
        range(0s; length=size(unreconstructedimages, 3), step=acquisitiontimestep),
        NA,
        M,
        marker,
        pixelsize,
    )
end # }}}

function times(acquisition::SIMAcquisition; images=:unreconstructedimages) # {{{
    if images == :unreconstructedimages
        return acquisition.unreconstructedimages[Axis{:time}].val
    elseif images == :reconstructedimages
        return acquisition.reconstructedimages[Axis{:time}].val
    else
        # TODO: error handling when not reconstructed or unreconstructed <20-05-21, kunzaatko> #
        return error("Must be reconstructed or unreconstructed")
    end
end # }}}


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

_mapreduce(f, op, ::IndexCartesian, A::AbstractQuasiArray) = mapfoldl(f, op, A)

## Functions to compute the reduced shape

# This differs from Arrays in that we always have Base.OneTo
reduced_index(i::Inclusion) = Base.OneTo(1)
reduced_indices(a::AbstractQuasiArray, region) = reduced_indices(axes(a), region)

function reduced_indices(inds::NTuple{N,Any}, d::Int) where N
    d < 1 && throw(ArgumentError("dimension must be ≥ 1, got $d"))
    if d == 1
        return (reduced_index(inds[1]), tail(inds)...)
    elseif 1 < d <= N
        return tuple(inds[1:d-1]..., reduced_index(inds[d]), inds[d+1:N]...)
    else
        return inds
    end
end

###### Generic reduction functions #####

# reducedim_initarray is called by
reducedim_initarray(A::AbstractQuasiArray, region, init, ::Type{R}) where {R} = fill!(similar(A,R,reduced_indices(A,region)), init)
reducedim_initarray(A::AbstractQuasiArray, region, init::T) where {T} = reducedim_initarray(A, region, init, T)

# function reducedim_init(f, op::Union{typeof(+),typeof(add_sum)}, A::AbstractQuasiArray, region)
#     _reducedim_init(f, op, zero, sum, A, region)
# end
# function reducedim_init(f, op::Union{typeof(*),typeof(mul_prod)}, A::AbstractQuasiArray, region)
#     _reducedim_init(f, op, one, prod, A, region)
# end

# # initialization when computing minima and maxima requires a little care
# for (f1, f2, initval, typeextreme) in ((:min, :max, :Inf, :typemax), (:max, :min, :(-Inf), :typemin))
#     @eval function reducedim_init(f, op::typeof($f1), A::AbstractQuasiArray, region)
#         # First compute the reduce indices. This will throw an ArgumentError
#         # if any region is invalid
#         ri = reduced_indices(A, region)

#         # Next, throw if reduction is over a region with length zero
#         any(i -> isempty(axes(A, i)), region) && _empty_reduce_error()

#         # Make a view of the first slice of the region
#         A1 = view(A, ri...)

#         if isempty(A1)
#             # If the slice is empty just return non-view version as the initial array
#             return copy(A1)
#         else
#             # otherwise use the min/max of the first slice as initial value
#             v0 = mapreduce(f, $f2, A1)

#             T = _realtype(f, promote_union(eltype(A)))
#             Tr = v0 isa T ? T : typeof(v0)

#             # but NaNs and missing need to be avoided as initial values
#             if (v0 == v0) === false
#                 # v0 is NaN
#                 v0 = $initval
#             elseif isunordered(v0)
#                 # v0 is missing or a third-party unordered value
#                 Tnm = nonmissingtype(Tr)
#                 # TODO: Some types, like BigInt, don't support typemin/typemax.
#                 # So a Matrix{Union{BigInt, Missing}} can still error here.
#                 v0 = $typeextreme(Tnm)
#             end
#             # v0 may have changed type.
#             Tr = v0 isa T ? T : typeof(v0)

#             return reducedim_initarray(A, region, v0, Tr)
#         end
#     end
# end
# reducedim_init(f::Union{typeof(abs),typeof(abs2)}, op::typeof(max), A::AbstractQuasiArray{T}, region) where {T} =
#     reducedim_initarray(A, region, zero(f(zero(T))), _realtype(f, T))

reducedim_init(f, op::typeof(&), A::AbstractQuasiArray, region) = reducedim_initarray(A, region, true)
reducedim_init(f, op::typeof(|), A::AbstractQuasiArray, region) = reducedim_initarray(A, region, false)

# specialize to make initialization more efficient for common cases

let
    BitIntFloat = Union{BitInteger, IEEEFloat}
    T = Union{
        [AbstractQuasiArray{t} for t in uniontypes(BitIntFloat)]...,
        [AbstractQuasiArray{Complex{t}} for t in uniontypes(BitIntFloat)]...}

    global function reducedim_init(f, op::Union{typeof(+),typeof(add_sum)}, A::T, region)
        z = zero(f(zero(eltype(A))))
        reducedim_initarray(A, region, op(z, z))
    end
    global function reducedim_init(f, op::Union{typeof(*),typeof(mul_prod)}, A::T, region)
        u = one(f(one(eltype(A))))
        reducedim_initarray(A, region, op(u, u))
    end
end

## generic (map)reduction

function _mapreducedim!(f, op, R::AbstractQuasiArray, A::AbstractQuasiArray)
    lsiz = check_reducedims(R,A)
    isempty(A) && return R

    indsAt, indsRt = safe_tail(axes(A)), safe_tail(axes(R)) # handle d=1 manually
    keep, Idefault = Broadcast.shapeindexer(indsRt)
    if reducedim1(R, A)
        # keep the accumulator as a local variable when reducing along the first dimension
        i1 = first(axes1(R))
        @inbounds for IA in QuasiCartesianIndices(indsAt)
            IR = Broadcast.newindex(IA, keep, Idefault)
            r = R[i1,IR]
            for i in axes(A, 1)
                r = op(r, f(A[i, IA]))
            end
            R[i1,IR] = r
        end
    else
        @inbounds for IA in QuasiCartesianIndices(indsAt)
            IR = Broadcast.newindex(IA, keep, Idefault)
            for i in axes(A, 1)
                R[i,IR] = op(R[i,IR], f(A[i,IA]))
            end
        end
    end
    return R
end

mapreducedim!(f, op, R::AbstractQuasiArray, A::AbstractQuasiArray) =
    (_mapreducedim!(f, op, R, A); R)

reducedim!(op, R::AbstractQuasiArray{RT}, A::AbstractQuasiArray) where {RT} =
    mapreducedim!(identity, op, R, A)

mapreduce(f, op, A::AbstractQuasiArray; dims=:, init=_InitialValue()) =
    _mapreduce_dim(f, op, init, A, dims)
mapreduce(f, op, A::AbstractQuasiArray...; kw...) =
    reduce(op, map(f, A...); kw...)

_mapreduce_dim(f, op, nt, A::AbstractQuasiArray, ::Colon) =
    mapfoldl_impl(f, op, nt, A)

_mapreduce_dim(f, op, ::_InitialValue, A::AbstractQuasiArray, ::Colon) =
    _mapreduce(f, op, IndexStyle(A), A)

_mapreduce_dim(f, op, nt, A::AbstractQuasiArray, dims) =
    mapreducedim!(f, op, reducedim_initarray(A, dims, nt), A)

_mapreduce_dim(f, op, ::_InitialValue, A::AbstractQuasiArray, dims) =
    mapreducedim!(f, op, reducedim_init(f, op, A, dims), A)

count(A::AbstractQuasiArray; dims=:, init=0) = count(identity, A; dims, init)
count(f, A::AbstractQuasiArray; dims=:, init=0) = _count(f, A, dims, init)

_count(f, A::AbstractQuasiArray, dims::Colon, init) = _simple_count(f, A, init)
_count(f, A::AbstractQuasiArray, dims, init) = mapreduce(_bool(f), add_sum, A; dims, init)

count!(r::AbstractQuasiArray, A::AbstractQuasiArray; init::Bool=true) = count!(identity, r, A; init=init)
count!(f, r::AbstractQuasiArray, A::AbstractQuasiArray; init::Bool=true) =
    mapreducedim!(_bool(f), add_sum, initarray!(r, add_sum, init, A), A)

for (fname, _fname, op) in [(:sum,     :_sum,     :add_sum), (:prod,    :_prod,    :mul_prod),
                            (:maximum, :_maximum, :max),     (:minimum, :_minimum, :min)]
    @eval begin
        # User-facing methods with keyword arguments
        @inline ($fname)(a::AbstractQuasiArray; dims=:, kw...) = ($_fname)(a, dims; kw...)
        @inline ($fname)(f, a::AbstractQuasiArray; dims=:, kw...) = ($_fname)(f, a, dims; kw...)
    end
end

any(a::AbstractQuasiArray; dims=:)              = _any(a, dims)
any(f::Function, a::AbstractQuasiArray; dims=:) = _any(f, a, dims)
all(a::AbstractQuasiArray; dims=:)              = _all(a, dims)
all(f::Function, a::AbstractQuasiArray; dims=:) = _all(f, a, dims)

for (fname, op) in [(:sum, :add_sum), (:prod, :mul_prod),
                    (:maximum, :max), (:minimum, :min),
                    (:all, :&),       (:any, :|)]
    fname! = Symbol(fname, '!')
    _fname = Symbol('_', fname)
    @eval begin
        $(fname!)(f::Function, r::AbstractQuasiArray, A::AbstractQuasiArray; init::Bool=true) =
            mapreducedim!(f, $(op), initarray!(r, $(op), init, A), A)
        $(fname!)(r::AbstractQuasiArray, A::AbstractQuasiArray; init::Bool=true) = $(fname!)(identity, r, A; init=init)
    end
end

# ##### findmin & findmax #####
# # The initial values of Rval are not used if the corresponding indices in Rind are 0.
# #
# function findminmax!(f, Rval, Rind, A::AbstractQuasiArray{T,N}) where {T,N}
#     (isempty(Rval) || isempty(A)) && return Rval, Rind
#     lsiz = check_reducedims(Rval, A)
#     for i = 1:N
#         axes(Rval, i) == axes(Rind, i) || throw(DimensionMismatch("Find-reduction: outputs must have the same indices"))
#     end
#     # If we're reducing along dimension 1, for efficiency we can make use of a temporary.
#     # Otherwise, keep the result in Rval/Rind so that we traverse A in storage order.
#     indsAt, indsRt = safe_tail(axes(A)), safe_tail(axes(Rval))
#     keep, Idefault = Broadcast.shapeindexer(indsRt)
#     ks = keys(A)
#     y = iterate(ks)
#     zi = zero(eltype(ks))
#     if reducedim1(Rval, A)
#         i1 = first(axes1(Rval))
#         @inbounds for IA in CartesianIndices(indsAt)
#             IR = Broadcast.newindex(IA, keep, Idefault)
#             tmpRv = Rval[i1,IR]
#             tmpRi = Rind[i1,IR]
#             for i in axes(A,1)
#                 k, kss = y::Tuple
#                 tmpAv = A[i,IA]
#                 if tmpRi == zi || f(tmpRv, tmpAv)
#                     tmpRv = tmpAv
#                     tmpRi = k
#                 end
#                 y = iterate(ks, kss)
#             end
#             Rval[i1,IR] = tmpRv
#             Rind[i1,IR] = tmpRi
#         end
#     else
#         @inbounds for IA in CartesianIndices(indsAt)
#             IR = Broadcast.newindex(IA, keep, Idefault)
#             for i in axes(A, 1)
#                 k, kss = y::Tuple
#                 tmpAv = A[i,IA]
#                 tmpRv = Rval[i,IR]
#                 tmpRi = Rind[i,IR]
#                 if tmpRi == zi || f(tmpRv, tmpAv)
#                     Rval[i,IR] = tmpAv
#                     Rind[i,IR] = k
#                 end
#                 y = iterate(ks, kss)
#             end
#         end
#     end
#     Rval, Rind
# end

# function findmin!(rval::AbstractQuasiArray, rind::AbstractQuasiArray, A::AbstractQuasiArray;
#                   init::Bool=true)
#     findminmax!(isgreater, init && !isempty(A) ? fill!(rval, first(A)) : rval, fill!(rind,zero(eltype(keys(A)))), A)
# end

# findmin(A::AbstractQuasiArray; dims=:) = _findmin(A, dims)

# function findmax!(rval::AbstractQuasiArray, rind::AbstractQuasiArray, A::AbstractQuasiArray;
#                   init::Bool=true)
#     findminmax!(isless, init && !isempty(A) ? fill!(rval, first(A)) : rval, fill!(rind,zero(eltype(keys(A)))), A)
# end
# findmax(A::AbstractQuasiArray; dims=:) = _findmax(A, dims)
# argmin(A::AbstractQuasiArray; dims=:) = findmin(A; dims=dims)[2]
# argmax(A::AbstractQuasiArray; dims=:) = findmax(A; dims=dims)[2]

# support overloading sum by MemoryLayout
_sum(V::AbstractQuasiArray, dims) = __sum(MemoryLayout(V), V, dims)
_sum(V::AbstractQuasiArray, ::Colon) = __sum(MemoryLayout(V), V, :)

_cumsum(A, dims) = __cumsum(MemoryLayout(A), A, dims)
cumsum(A::AbstractQuasiArray; dims::Integer) = _cumsum(A, dims)
cumsum(x::AbstractQuasiVector) = cumsum(x, dims=1)

# sum is equivalent to hitting by ones(n) on the left or right

__cumsum(::QuasiArrayLayout, A, ::Colon) = QuasiArray(cumsum(parent(A)), axes(A))
__cumsum(::QuasiArrayLayout, A, d::Int) = QuasiArray(cumsum(parent(A),dims=d), axes(A))

for Sum in (:sum, :cumsum)
    __Sum = Symbol("__", Sum)
    @eval function $__Sum(LAY::ApplyLayout{typeof(*)}, V::AbstractQuasiMatrix, d::Int)
        a = arguments(LAY, V)
        if d == 1
            *($Sum(first(a); dims=1), tail(a)...)
        else
            @assert d == 2
            *(most(a)..., $Sum(last(a); dims=2))
        end
    end
end

function __cumsum(LAY::ApplyLayout{typeof(*)}, V::AbstractQuasiVector, dims)
    a = arguments(LAY, V)
    apply(*, cumsum(a[1]; dims=dims), tail(a)...)
end

function __sum(LAY::ApplyLayout{typeof(*)}, V::AbstractQuasiVector, ::Colon)
    a = arguments(LAY, V)
    first(apply(*, sum(a[1]; dims=1), tail(a)...))
end

__sum(_, A, dims) = _sum(identity, A, dims)


text:
using DataDeps

register(DataDep(
  "thyroid-disease-allhypo",
  "http://archive.ics.uci.edu/ml/datasets/Thyroid+Disease",
  "http://archive.ics.uci.edu/ml/machine-learning-databases/thyroid-disease/allhypo.data",
  "273bfa43c6ad226204cb42d17365143723786db447b074abb71a23a632c1ff53",
  post_fetch_method=(path -> begin
    # The ids are appended to each row after `.|`. Replace with comma
    write(path, replace(read(path, String), ".|" => ","))

    UCIData.process_dataset(path,
        target_index=30,
        feature_indices=1:29,
        id_indices=31,
        categoric_indices=[2:17; 19; 21; 23; 25; 27; 29],
    )
  end),
))

max_stars_repo_name: JackDunnNZ/uci-data

max_issues_repo_name: JackDunnNZ/uci-data

max_forks_repo_name: JackDunnNZ/uci-data


text:
"""
    msmplot(T; pi_i=nothing, x=nothing, y=nothing, filename=nothing,
            edgewidth_scale=10.0, arrow_scale=0.1, nodesize=0.5, fontsize=10, names=[], dpi=100)

Visualize the graphical structure of the given Markov state model parameters.
`T` is a transition probability matrix whose elements Tij represents the probablity of transition from state i to j.
`T` should satisfy the detailed balance condition. `pi_i` is a vector whose elements are equilibrium probabilities of states.
`x` and `y` are X and Y coordinates of states, respectively.

# Examples
```julia-repl
julia> nstate = 5
julia> T, pi_i = msmtransitionmatrix(rand(nstate, nstate))
julia> x = rand(nstate); y = rand(nstate)
julia> msmplot(T, pi_i=pi_i, x=x, y=y)
```
"""
function msmplot(T; pi_i=nothing, x=nothing, y=nothing, filename=nothing,
                 edgewidth_scale=3.0, arrow_scale=0.0001, nodesize=0.5, fontsize=10, names=[], dpi=100)
  n = size(T, 1)
  g = MetaDiGraph(n, 0.0)
  U = typeof(T[1, 1])

  for i = 1:n
    for j = 1:n
      if !iszero(T[i, j])
        Graphs.SimpleGraphs.add_edge!(g, i, j)
      end
     end
  end

  for i = 1:n
    for j = 1:n
      if !iszero(T[i, j])
        set_prop!(g, Edge(i, j), :weight, T[i, j])
      end
    end
  end

  if isnothing(pi_i)
    pi_i = ones(1, n)./n
    for i = 1:(n*100)
      pi_i .= pi_i * T
    end
    pi_i = pi_i./sum(pi_i)
    pi_i = pi_i[:]
  end

  if !isnothing(x)
    std_x = std(x)
    std_y = std(y)
    max_xy = maximum([maximum(x),maximum(y)])
    xx = (x[:].-mean(x))./std_x .- std_x
    yy = (y[:].-mean(y))./std_y .- std_y
  else
    xx = nothing
    yy = nothing
  end

  gr()
#  pyplot()
  p = graphplot(g, edge_width=(s,d,w) -> get_prop(g, Edge(s,d), :weight)*edgewidth_scale,
            #arrow=arrow(:closed, :head, 0.1, 0.1),
            arrow=arrow(:closed, :head, arrow_scale, arrow_scale),
            #markersize = markersize,
            node_weights = pi_i[:],
            #node_weights = 1:n,
            #markercolor = range(colorant"yellow", stop=colorant"red", length=n),
            markercolor = :white,
            names = names,
            nodesize = nodesize,
            fontsize = fontsize,
            fontcolor = :white,
            nodeshape = :circle,
            linecolor = :darkgrey,
            #shorten=0.95,
            x = xx,
            y = yy,
            curves = true,
            curvature_scalar = 0.8,
            dpi = dpi
            )

  if !isnothing(filename)
    savefig(p, filename)
  end

  return p
end

function msmsample(p)
  p_cum = cumsum(p) ./ sum(p)
  r = rand()
  for i = 1:length(p_cum)
    if r <= p_cum[i]
      return i
    end
  end
end

"""
    msmgenerate(nframe::Int, T, pi_i)
                states = zeros(typeof(nframe), nframe)

Randomly samples a state trajectory from the the given transition matrix, and equilibrium probabilities of states.

# Examples
```julia-repl
julia> T, pi_i = msmtransitionmatrix(C)
julia> states = msmgenerate(1000, T, pi_i)
```
"""
function msmgenerate(nframe::Int, T, pi_i)
  states = zeros(typeof(nframe), nframe)

  states[1] = msmsample(pi_i)
  for iframe = 2:nframe
    states[iframe] = msmsample(T[states[iframe-1], :])
  end
  return states
end

"""
    msmgenerate(nframe, T, pi_i, emission) -> states, observations

Randomly samples a state trajectory and observations from the given transition matrix, equilibrium probabilities of states, and emissions.

# Examples
```julia-repl
julia> T, pi_i = msmtransitionmatrix(C)
julia> states, observations = msmgenerate(1000, T, pi_i, emission)
```
"""
function msmgenerate(nframe::Int, T, pi_i, emission)
  states = zeros(typeof(nframe), nframe)
  observations = zeros(typeof(nframe), nframe)

  states[1] = msmsample(pi_i)
  observations[1] = msmsample(emission[states[1], :])
  for iframe = 2:nframe
    states[iframe] = msmsample(T[states[iframe-1], :])
    observations[iframe] = msmsample(emission[states[iframe], :])
  end
  return states, observations
end

"""
    msmcountmatrix(indexOfCluster; tau=1) -> C::Matrix

Transition count matrix from a sinlgle binned trajectory or a set of binned trajectories.
`indexOfCluster` is a vector or a set of vectors containing binned trajectorie(s).
Returns count matrix for transitions from state i to state j in a time step of tau.

# Examples
```julia-repl
julia> ta = mdload("ak.dcd")
julia> X = compute_distancemap(ta["atomname CA"])
julia> F = clusterkcenters(X)
julia> c = msmcountmatrix(F.indexOfCluster, tau=10)
```
"""
function msmcountmatrix(indexOfCluster; tau=1)
  if eltype(indexOfCluster) <: Number
    indexOfCluster2 = [indexOfCluster]
  else
    indexOfCluster2 = indexOfCluster
  end

  nstate = maximum(maximum(indexOfCluster2))
  ntraj = length(indexOfCluster2)
  U = typeof(indexOfCluster2[1][1])
  c = spzeros(U, nstate, nstate)
  for itraj = 1:ntraj
    nframe = length(indexOfCluster2[itraj])

    index_from = 1:(nframe-tau)
    index_to   = (1+tau):nframe
    indexOfCluster_from = indexOfCluster2[itraj][index_from]
    indexOfCluster_to   = indexOfCluster2[itraj][index_to]

    s = ones(U, length(indexOfCluster_from))
    c_1 = sparse(indexOfCluster_from, indexOfCluster_to, s, nstate, nstate)
    c .= c .+ c_1
  end

  return Array(c)
end

"""
    msmtransitionmatrix(C; TOLERANCE=10^(-4), verbose=true) -> T::Matrix, p::Vector

Estimate the transition probability matrix from count matrix `C`.
Detailed balance is implicitly imposed in the estimation.

Returns the transition probability matrix and the equilibrium probabilities of states.

# Examples
```julia-repl
julia> ta = mdload("ak.dcd")
julia> X = compute_distancemap(ta["atomname CA"])
julia> F = clusterkcenters(X)
julia> C = msmcountmatrix(F.indexOfCluster, tau=10)
julia> T, p = msmtransitionmatrix(C)
```
# References
```
This routines uses the reversible maximum likelihood estimator described in
K. A. Beauchamp, G. R. Bowman, T. J. Lane, L. Maibaum, I. S. Haque, and V. S. Pande,
MSMBuilder2: Modeling Conformational Dynamics on the Picosecond to Millisecond Scale,
J. Chem. Theory Comput. 7, 3412 (2011).
```
"""
function msmtransitionmatrix(C; TOLERANCE=10^(-4), verbose=true)
  c = Matrix{Float64}(C)
  nstate = size(c, 1)

  c_sym = c + c'
  x = c_sym

  c_rs = sum(c, dims=2)
  x_rs = sum(x, dims=2)

  logL_old = 2.0*TOLERANCE
  logL = 0.0
  count_iteration = 0

  x = zeros(Float64, nstate, nstate)
  x_new = zeros(Float64, nstate, nstate)

  while abs(logL_old - logL) > TOLERANCE
    count_iteration = count_iteration + 1
    logL_old = logL

    # fixed-point method
    for i = 1:nstate
      for j = 1:nstate
        denom = (c_rs[i]/x_rs[i]) + (c_rs[j]/x_rs[j])
        x_new[i, j] = (c[i, j] + c[j, i]) / denom
      end
    end

    # update
    x_rs .= sum(x_new, dims=2)
    x .= x_new
    logL = 0.0
    for i = 1:nstate
      for j = 1:nstate
        if !iszero(x[i, j]) & !iszero(x_rs[i])
          logL = logL + c[i, j] * log(x[i, j] / x_rs[i]);
        end
      end
    end

    if verbose & (mod(count_iteration, 10) == 0)
      Printf.@printf("%d iteration  LogLikelihood = %8.5e  delta = %8.5e  tolerance = %8.5e\n", count_iteration, logL, abs(logL_old-logL), TOLERANCE)
    end
  end

  pi_i = x_rs./sum(x_rs)
  pi_i = pi_i[:]
  t = x ./ x_rs

  return t, pi_i
end

function msmforward(data_list, T, pi_i, emission)
    ndata = length(data_list)
    nstate = length(T[1, :])
    logL = zeros(Float64, ndata)
    alpha_list = []
    factor_list = []
    for idata = 1:ndata
        data = data_list[idata]
        nframe = length(data)
        alpha  = zeros(Float64, (nframe, nstate))
        factor = zeros(Float64, nframe)
        alpha[1, :] = pi_i.*emission[:, data[1]]
        factor[1] = sum(alpha[1, :])
        alpha[1, :] = alpha[1, :]./factor[1]
        for iframe = 2:nframe
            alpha[iframe, :] = sum(alpha[iframe-1, :] .* T, dims=1)' .* emission[:, data[iframe]]
            factor[iframe] = sum(alpha[iframe, :])
            alpha[iframe, :] = alpha[iframe, :]./factor[iframe]
        end
        logL[idata] = sum(log.(factor))
        push!(alpha_list, alpha)
        push!(factor_list, factor)
    end
    logL, alpha_list, factor_list
end

function msmforward_missing(data_list, T, pi_i, emission)
  ndata = length(data_list)
  nstate = length(T[1, :])
  logL = zeros(Float64, ndata)
  alpha_list = []
  factor_list = []
  for idata = 1:ndata
      data = data_list[idata]
      nframe = length(data)
      alpha  = zeros(Float64, (nframe, nstate))
      factor = zeros(Float64, nframe)
      if data[1] === missing
        alpha[1, :] = pi_i
      else
        alpha[1, :] = pi_i.*emission[:, data[1]]
      end
      factor[1] = sum(alpha[1, :])
      alpha[1, :] = alpha[1, :]./factor[1]
      for iframe = 2:nframe
          if data[iframe] === missing
            alpha[iframe, :] = sum(alpha[iframe-1, :] .* T, dims=1)'
          else
            alpha[iframe, :] = sum(alpha[iframe-1, :] .* T, dims=1)' .* emission[:, data[iframe]]
          end
          factor[iframe] = sum(alpha[iframe, :])
          alpha[iframe, :] = alpha[iframe, :]./factor[iframe]
      end
      logL[idata] = sum(log.(factor))
      push!(alpha_list, alpha)
      push!(factor_list, factor)
  end
  logL, alpha_list, factor_list
end

function msmbackward(data_list, factor_list, T, pi_i, emission)
    ndata = length(data_list)
    nstate = length(T[1, :])
    logL = zeros(Float64, ndata)
    beta_list = []
    for idata = 1:ndata
        data   = data_list[idata]
        factor = factor_list[idata]
        nframe = length(data)
        beta   = zeros(Float64, (nframe, nstate))
        beta[nframe, :] .= 1.0
        for iframe = (nframe-1):-1:1
            beta[iframe, :] = sum((T .* (emission[:, data[iframe+1]] .* beta[iframe+1, :])'), dims=2) ./ factor[iframe+1]
        end
        logL[idata] = sum(log.(factor))
        push!(beta_list, beta)
    end
    logL, beta_list
end

function msmbackward_missing(data_list, factor_list, T, pi_i, emission)
    ndata = length(data_list)
    nstate = length(T[1, :])
    logL = zeros(Float64, ndata)
    beta_list = []
    for idata = 1:ndata
        data   = data_list[idata]
        factor = factor_list[idata]
        nframe = length(data)
        beta   = zeros(Float64, (nframe, nstate))
        beta[nframe, :] .= 1.0
        for iframe = (nframe-1):-1:1
            if data[iframe+1] === missing
                beta[iframe, :] = sum((T .* beta[iframe+1, :]'), dims=2) ./ factor[iframe+1]
            else
                beta[iframe, :] = sum((T .* (emission[:, data[iframe+1]] .* beta[iframe+1, :])'), dims=2) ./ factor[iframe+1]
            end
        end
        logL[idata] = sum(log.(factor))
        push!(beta_list, beta)
    end
    logL, beta_list
end

"""
    msmbaumwelch(observations, T_init, p_init, emission_init; TOLERANCE=10.0^(-4), MAXITERATION=Inf64) -> T::Matrix, p::Vector, emission::Matrix

Baum-Welch algorithm estimates the most probable transition probabilities from the given observation data.
In this function, detailed balance is implicitly imposed in the estimation,
so the equilibrium probabilities can be determined from the estimated transition probabilities.
Also, unlike the original Baum-Welch algorithm, the emission is NOT estimated in this function,
because the emission probabilites are usually known a priori in cases of molecular experiments.
`observations` is a set of observed vectors. `T_init`, `p_init` are initial transition probabilities and equilibrium probabilities, respectively.
`emission_init` is a matrix whose rows correspond to states, and columns correspond to observations.

Returns the transition probability matrix, the equilibrium probabilities of states,
and the emission probabilities (though the emission does not change).

# Examples
```julia-repl
julia> ta = mdload("ak.dcd")
julia> X = compute_distancemap(ta["atomname CA"])
julia> F = clusterkcenters(X)
julia> C = msmcountmatrix(F.indexOfCluster, tau=10)
julia> T, p = msmtransitionmatrix(C)
```
# References
```
The algorithm of this routines is based on the descriptions in PRML book by C. Bishop.
```
"""
function msmbaumwelch(data_list, T0, pi_i0, emission0; TOLERANCE=10.0^(-4), MAXITERATION=Inf64)
    ## setup
    check_convergence = Inf64
    count_iteration = 0
    logL_old = 1.0
    #if not isinstance(data_list, list):
    #    data_list = [data_list]
    ndata = length(data_list)
    nobs = length(emission0[1, :])
    nstate = length(T0[1, :])
    T = similar(T0)
    emission = similar(emission0)
    pi_i = similar(pi_i0)
    while (check_convergence > TOLERANCE) & (count_iteration <= MAXITERATION)
        ## E-step
        logL, alpha_list, factor_list = msmforward(data_list, T0, pi_i0, emission0)
        #print("1"); println(logL)
        logL2, beta_list = msmbackward(data_list, factor_list, T0, pi_i0, emission0)
        #print("2"); println(logL2)
        log_alpha_list = []
        for a in alpha_list
            push!(log_alpha_list, log.(a))
        end
        log_beta_list = []
        for b in beta_list
            push!(log_beta_list, log.(b))
        end
        log_T0 = log.(T0)
        log_emission0 = log.(emission0)

        ## M-step
        # pi
        # pi = np.zeros(nstate, dtype=np.float64)
        # log_gamma_list = []
        # for idata in range(ndata):
        #     log_gamma_list.append(log_alpha_list[idata] + log_beta_list[idata])
        #     pi = pi + np.exp(log_gamma_list[idata][0, :])
        # pi = pi/np.sum(pi)
        # pi_i = pi_i0
        # emission
        # emission = np.zeros((nstate, nobs), dtype=np.float64)
        # for idata in range(ndata):
        #     data = data_list[idata]
        #     for istate in range(nstate):
        #         for iobs in range(nobs):
        #             id = (data == iobs)
        #             if np.any(id):
        #                 emission[istate, iobs] = emission[istate, iobs] + np.sum(np.exp(log_gamma_list[idata][id, istate]))
        # emission[np.isnan(emission)] = 0.0
        # emission = emission / np.sum(emission, axis=1)[:, None]
        emission = emission0
        # T
        T = zeros(Float64, (nstate, nstate))
        for idata = 1:ndata
          data = data_list[idata]
          nframe = length(data)
          for iframe = 2:nframe
            #log_xi = bsxfun(@plus, log_alpha{idata}(iframe-1, :)', log_beta{idata}(iframe, :));
            log_xi = log_alpha_list[idata][iframe-1, :] .+ log_beta_list[idata][iframe, :]'
            #T = T .+ exp(bsxfun(@plus, log_xi, log_emission0(:, data(iframe))') + log_T0)./factor{idata}(iframe);
            T = T .+ exp.((log_xi .+ log_emission0[:, data[iframe]]') .+ log_T0) ./ factor_list[idata][iframe]
          end
        end
        #T[np.isnan(T)] = 0.0
        T = T ./ sum(T, dims=2)

        ## reversible T
        T, pi_i = msmtransitionmatrix(T, TOLERANCE=10^(-8), verbose=false)

        ## Check convergence
        count_iteration += 1
        logL = sum(logL)
        check_convergence = abs(logL_old - logL)
        if mod(count_iteration, 100) == 0
            Printf.@printf("%d iteration LogLikelihood = %e  delta = %e  tolerance = %e\n" , count_iteration, logL, check_convergence, TOLERANCE)
        end
        logL_old = logL
        pi_i0 = pi_i
        emission0 = emission
        T0 = T
    end
    T, pi_i, emission
end

function msmbaumwelch_missing(data_list, T0, pi_i0, emission0; TOLERANCE=10.0^(-4), MAXITERATION=Inf64)
  ## setup
  check_convergence = Inf64
  count_iteration = 0
  logL_old = 1.0
  #if not isinstance(data_list, list):
  #    data_list = [data_list]
  ndata = length(data_list)
  nobs = length(emission0[1, :])
  nstate = length(T0[1, :])
  T = similar(T0)
  emission = similar(emission0)
  pi_i = similar(pi_i0)
  while (check_convergence > TOLERANCE) & (count_iteration <= MAXITERATION)
      ## E-step
      logL, alpha_list, factor_list = msmforward_missing(data_list, T0, pi_i0, emission0)
      #print("1"); println(logL)
      logL2, beta_list = msmbackward_missing(data_list, factor_list, T0, pi_i0, emission0)
      #print("2"); println(logL2)
      log_alpha_list = []
      for a in alpha_list
          push!(log_alpha_list, log.(a))
      end
      log_beta_list = []
      for b in beta_list
          push!(log_beta_list, log.(b))
      end
      log_T0 = log.(T0)
      log_emission0 = log.(emission0)

      ## M-step
      # pi
      # pi = np.zeros(nstate, dtype=np.float64)
      # log_gamma_list = []
      # for idata in range(ndata):
      #     log_gamma_list.append(log_alpha_list[idata] + log_beta_list[idata])
      #     pi = pi + np.exp(log_gamma_list[idata][0, :])
      # pi = pi/np.sum(pi)
      # pi_i = pi_i0
      # emission
      # emission = np.zeros((nstate, nobs), dtype=np.float64)
      # for idata in range(ndata):
      #     data = data_list[idata]
      #     for istate in range(nstate):
      #         for iobs in range(nobs):
      #             id = (data == iobs)
      #             if np.any(id):
      #                 emission[istate, iobs] = emission[istate, iobs] + np.sum(np.exp(log_gamma_list[idata][id, istate]))
      # emission[np.isnan(emission)] = 0.0
      # emission = emission / np.sum(emission, axis=1)[:, None]
      emission = emission0
      # T
      T = zeros(Float64, (nstate, nstate))
      for idata = 1:ndata
        data = data_list[idata]
        nframe = length(data)
        for iframe = 2:nframe
          log_xi = log_alpha_list[idata][iframe-1, :] .+ log_beta_list[idata][iframe, :]'
          if data[iframe] === missing
            T = T .+ exp.(log_xi .+ log_T0) ./ factor_list[idata][iframe]
          else
            T = T .+ exp.((log_xi .+ log_emission0[:, data[iframe]]') .+ log_T0) ./ factor_list[idata][iframe]
          end
        end
      end
      #T[np.isnan(T)] = 0.0
      T = T ./ sum(T, dims=2)

      ## reversible T
      T, pi_i = msmtransitionmatrix(T, TOLERANCE=10^(-8), verbose=false)

      ## Check convergence
      count_iteration += 1
      logL = sum(logL)
      check_convergence = abs(logL_old - logL)
      if mod(count_iteration, 100) == 0
          Printf.@printf("%d iteration LogLikelihood = %e  delta = %e  tolerance = %e\n" , count_iteration, logL, check_convergence, TOLERANCE)
      end
      logL_old = logL
      pi_i0 = pi_i
      emission0 = emission
      T0 = T
  end
  T, pi_i, emission
end

"""
    msmviterbi(observation, T, p, emission) -> states::Vector

Viterbi algorithm estimates the most probable hidden state sequence from the observation data.
`observations` is a set of observed vectors. `T`, `p` are the transition probabilities and equilibrium probabilities, respectively.
`emission` is a matrix whose rows correspond to states, and columns correspond to observations.

Returns the transition probability matrix, the equilibrium probabilities of states,
and the emission probabilities (though the emission does not change).

# Examples
```julia-repl
julia> nframe = 1000
julia> states, observations = msmgenerate(nframe, T, pi_i, emission)
julia> states_estimated = msmviterbi(T, pi_i, emission, observation)
```
# References
```
The algorithm of this routines is based on the descriptions in PRML book by C. Bishop.
```
"""
function msmviterbi(observation, T, pi_i, emission)
    T2 = T .+ 10^(-100)
    T2 .= T2 ./ sum(T2, dims=2)
    pi_i2 = pi_i .+ 10^(-100)
    pi_i2 = pi_i2 ./sum(pi_i2)
    emission2 = emission .+ 10^(-100)
    emission2 .= emission2 ./ sum(emission2, dims=1)

    nframe = size(observation, 1)
    nstate = size(T2, 1)
    P = zeros(eltype(T), nstate, nframe)
    I = zeros(eltype(T), nstate, nframe)
    state_estimated = zeros(eltype(observation), nframe)

    # initialization
    if observation[1] === missing
      P[:, 1] .= log.(pi_i2)
    else
      P[:, 1] .= log.(pi_i2) .+ log.(emission2[:, observation[1]])
    end
    I[:, 1] .= zeros(eltype(T2), nstate)

    # argmax forward
    Z = zeros(eltype(T2), nstate, nstate)
    for t = 2:nframe
        Z .= P[:, (t-1):(t-1)] .+ log.(T2)
        #if observation[t] === missing
        #  Z .= P[:, (t-1):(t-1)] .+ log.(T2)
        #else
        #  Z .= P[:, (t-1):(t-1)] .+ log.(T2) .+ log.(emission2[:, observation[t]:observation[t]]')
        #end
        I[:, t] .= getindex.(argmax(Z, dims=1), 1)[:]
        if observation[t] === missing
          P[:, t] .= maximum(Z, dims=1)[:]
        else
          P[:, t] .= maximum(Z, dims=1)[:] .+ log.(emission[:, observation[t]])
        end
    end

    # termination
    P_star = maximum(P[:, nframe])
    state_estimated[nframe] = argmax(P[:, nframe])
    #@show P[:, nframe]

    # decoding
    for t = (nframe-1):-1:1
        state_estimated[t] = I[state_estimated[t+1], t+1]
    end

    return state_estimated
end

function msmimpliedtime(indexOfCluster, tau)
  n = length(tau)
  k = min(10, n)
  nstate = Int(maximum(indexOfCluster))
  #implied_timescale = zeros(Float64, n, k)
  implied_timescale = Matrix{Union{Float64, Missing}}(undef, n, nstate-1)
  for i in 1:n
    t = tau[i]
    #@printf "calculating tau = %d\n" t
    C = msmcountmatrix(indexOfCluster, tau=t)
    # TODO fixme; tarjan algorithm to remove non-ergodic graphs
    T, p = msmtransitionmatrix(C, TOLERANCE=10^(-8), verbose=false)
    F = eigen(T, sortby = x -> -real(x))
    #@show F.values
    for istate = 2:nstate
      r = real(F.values[istate])
      if r >  0.0
        implied_timescale[i, istate-1] = - t / log(r)
      else
        implied_timescale[i, istate-1] = missing
      end
    end
  end
  return implied_timescale
end


text:
using WaterLily
using GLMakie
using FileIO
using MeshIO

# Error ArgumentError: invalid base 10 digit '\\' in "\\"
path = abspath(joinpath(@__DIR__, "assets", "Santasleigh.obj"))
sleigh = load(path)


text:
import Interpolations: interpolate, scale

mutable struct InterpKDE{K,I} <: AbstractKDE
    kde::K
    itp::I
    InterpKDE{K,I}(kde::K, itp::I) where {K,I} = new{K,I}(kde, itp)
end

function InterpKDE(kde::UnivariateKDE, opts...)
    itp_u = interpolate(kde.density, opts...)
    itp = scale(itp_u, kde.x)
    etp = extrapolate(itp, zero(eltype(kde.density)))
    InterpKDE{typeof(kde),typeof(etp)}(kde, etp)
end
InterpKDE(kde::UnivariateKDE) = InterpKDE(kde, BSpline(Quadratic(Line(OnGrid()))))

function InterpKDE(kde::BivariateKDE, opts...)
    itp_u = interpolate(kde.density,opts...)
    itp = scale(itp_u, kde.x, kde.y)
    etp = extrapolate(itp, zero(eltype(kde.density)))
    InterpKDE{typeof(kde),typeof(etp)}(kde, etp)
end
InterpKDE(kde::BivariateKDE) = InterpKDE(kde::BivariateKDE, BSpline(Quadratic(Line(OnGrid()))))

pdf(ik::InterpKDE,x::Real...) = ik.itp(x...)
pdf(ik::InterpKDE,xs::AbstractVector) = [ik.itp(x) for x in xs]
pdf(ik::InterpKDE,xs::AbstractVector,ys::AbstractVector) = [ik.itp(x,y) for x in xs, y in ys]

pdf(k::UnivariateKDE,x) = pdf(InterpKDE(k),x)
pdf(k::BivariateKDE,x,y) = pdf(InterpKDE(k),x,y)


text:
using Rebugger, Test

@info "These tests manipulate the console. Wait until you see \"Done\""
include("edit.jl")
include("interpret.jl")
if Sys.isunix() && VERSION >= v"1.1.0"
    include("interpret_ui.jl")
else
    @warn "Skipping UI tests"
end
println("Done")


text:
import Base: size

export MultiBasis, getbasis

"""
$(TYPEDEF)

An immutable structure to hold the Nx-cartesian product of the 1D basis of functions B
MultiB = B × B× ... × B (Nx elements). The basis B contains m elements.

## Fields

$(TYPEDFIELDS)

"""
struct MultiBasis
    B::Basis
    Nx::Int64
end

"""
$(TYPEDSIGNATURES)

Returns the size of the `MultiBasis` `B`.
"""
size(B::MultiBasis) = (B.B.m, B.Nx)

"""
$(TYPEDSIGNATURES)

Returns the kind of the underlying basis of the `MultiBasis` `B`.
"""
getbasis(B::MultiBasis) = string(typeof(B.B))


text:
"""
Authors : Adrien BLASSIAU Corentin JUVIGNY
TP RODD - 2019-2020

Programme pour lancer une batterie de tests.

Utilisation :
Supposons pour l'exemple que vous avez un algorithme dans le fichier "pl.jl"
et un dossier "instances" contenant des fichier ".input".
Placez vous à la racine du dossier du projet (dans le dossier contenant pl.jl et helpers)
et tapez la commande suivante.

julia helpers/tester.jl algorithm_test.jl instances/

Cette commande exécutera l'algorithme pl.jl sur tous les fichiers ".input" du dossier.

"""

println("Chargement des bibliothèques.")

using Distributions
using JuMP
using CPLEX

m = Model(with_optimizer(CPLEX.Optimizer))

println("Chargement terminé.")

# Ne pas enlever
if length(ARGS) == 2
  algorithm = ARGS[1]
  test_folder = ARGS[2]
  if test_folder[length(test_folder)] != '/'
    test_folder *= "/"
  end

  while length(ARGS) > 0
    pop!(ARGS)
  end
  include("../" * algorithm)
  for input_file in readdir(test_folder)
    if(endswith(input_file, ".dat"))
      println(input_file)
      main(test_folder * input_file)
    end
  end
end

max_stars_repo_name: AdrienBlassiau/RODD

max_issues_repo_name: AdrienBlassiau/RODD

max_forks_repo_name: AdrienBlassiau/RODD


text:
@testset "VariableSelection" begin
    include("mindomain.jl")
    include("random.jl")
    include("failureBased.jl")
end


text:
module Types

export jint, jlong, jbyte, jboolean, jchar, jshort, jfloat, jdouble, jsize,
    jvoid, jobject, jclass, jthrowable, jweak, jmethodID, jfieldID, jstring, jarray,
    JNINativeMethod, jobjectArray, jbooleanArray, jbyteArray, jshortArray, jintArray,
    jlongArray, jfloatArray, jdoubleArray, jcharArray, jvalue, jobjectRefType

# jni_md.h
const jint = Cint
#ifdef _LP64 /* 64-bit Solaris */
# typedef long jlong;
const jlong = Clonglong
const jbyte = Cchar

# jni.h

const jboolean = Cuchar         # typedef unsigned char   jboolean;
const jchar = Cushort           # typedef unsigned short  jchar;
const jshort = Cshort           # typedef short           jshort;
const jfloat = Cfloat           # typedef float           jfloat;
const jdouble = Cdouble         # typedef double          jdouble;

const jsize = jint              # typedef jint            jsize;

const jvoid = Cvoid

struct _jobject end             # struct _jobject;

const jobject = Ptr{_jobject}   # typedef struct _jobject *jobject;
const jclass = jobject          # typedef jobject jclass;
const jthrowable = jobject      # typedef jobject jthrowable;
const jstring = jobject         # typedef jobject jstring;
const jarray = jobject          # typedef jobject jarray;

const jbooleanArray = jarray    # typedef jarray jbooleanArray;
const jbyteArray = jarray       # typedef jarray jbyteArray;
const jcharArray = jarray       # typedef jarray jcharArray;
const jshortArray = jarray      # typedef jarray jshortArray;
const jintArray = jarray        # typedef jarray jintArray;
const jlongArray = jarray       # typedef jarray jlongArray;
const jfloatArray = jarray      # typedef jarray jfloatArray;
const jdoubleArray = jarray     # typedef jarray jdoubleArray;
const jobjectArray = jarray     # typedef jarray jobjectArray;

const jweak = jobject           # typedef jobject jweak;

const jvalue = Union{           # typedef union jvalue {
    jboolean,                   #   jboolean z;
    jbyte,                      #   jbyte    b;
    jchar,                      #   jchar    c;
    jshort,                     #   jshort   s;
    jint,                       #   jint     i;
    jlong,                      #   jlong    j;
    jfloat,                     #   jfloat   f;
    jdouble,                    #   jdouble  d;
    jobject                     #   jobject  l;
}                               # } jvalue;

struct _jfieldID end                # struct _jfieldID;
const jfieldID = Ptr{_jfieldID}     # typedef struct _jfieldID *jfieldID;

struct _jmethodID end               # struct _jmethodID;
const jmethodID = Ptr{_jmethodID}   # typedef struct _jmethodID *jmethodID;

@enum jobjectRefType begin          # typedef enum _jobjectType {
    JNIInvalidRefType    = 0        #   JNIInvalidRefType    = 0,
    JNILocalRefType      = 1        #   JNILocalRefType      = 1,
    JNIGlobalRefType     = 2        #   JNIGlobalRefType     = 2,
    JNIWeakGlobalRefType = 3        #   JNIWeakGlobalRefType = 3
end                                 # } jobjectRefType;

struct JNINativeMethod      # typedef struct {
    name::Cstring           #   char *name;
    signature::Cstring      #   char *name;
    fnPtr::Ptr{Cvoid}       #   void *fnPtr;
end                         # } JNINativeMethod;

end


text:
module NumericalDifferentiation

export differentiate, integrationoperator, TotalVariation, Tikhonov

using LinearAlgebra, SparseArrays
using LinearMaps, IterativeSolvers

abstract type DifferentiationMethod end

struct TotalVariation{T<:AbstractFloat} <: DifferentiationMethod
    ϵ::T
end
TotalVariation() = TotalVariation(1e-12)

struct Tikhonov <: DifferentiationMethod end

"""
	diffmatrix(x::AbstractVector)

Generate differentiation matrix from node positions `x`.
"""
function diffmatrix(x::AbstractVector)
    D = zeros(eltype(x), length(x)-1, length(x))
    for i in 1:size(D,1)
        D[i,i]   = -1/(x[i+1]-x[i])
        D[i,i+1] =  1/(x[i+1]-x[i])
    end
    return D
end

"""
	spdiffmatrix(x::AbstractVector)

Generate sparse differentiation matrix from node positions `x`.
"""
function spdiffmatrix(x::AbstractVector)
    el = 1 ./diff(x)
    I = [1:length(x)-1;1:length(x)-1]
    J = [1:length(x)-1;2:length(x)]
    K = [-el;el]

    return sparse(I,J,K)
end

"""
	intgmatrix(x::AbstractVector)

Generate integration matrix using Trapezoidal rule from node positions `x`.
"""
function intgmatrix(x::AbstractVector)
    K = zeros(eltype(x), length(x), length(x))
    for  j in 2:size(K,2), i in j:size(K,1)
        K[i,j-1] += .5*(x[j] - x[j-1])
        K[i,j]   += .5*(x[j] - x[j-1])
    end
    return K
end

"""
	integrationoperator(x,y)

Integration operator using trapezoidal rule. Return the integral ``\\int_0^x y(u)du`` for all `x`. This code is reproduced from NumericalIntegration.jl (copied to avoid unnecessary dependencies and compatibility issues).
"""
function integrationoperator(x,y)
    retarr = zeros(eltype(x),size(x)) + zeros(eltype(y),size(y))
    for i in 2 : length(y)
        retarr[i] = retarr[i-1] + (x[i] - x[i-1]) * (y[i-1] + y[i])
    end
    return 0.5 * retarr
end

"""
	integrationadjoint(x,y)

Adjoint of integration operator. Constructed such that the corresponding matrix is truly the transpose of the matrix of the integration operator.
"""
function integrationadjoint(x,y)
    z = zeros(eltype(x),size(x)) + zeros(eltype(y),size(y))
    N = length(z)

    z[N] = y[N]*(x[N]-x[N-1])

    sumy = y[N]

    for i in N-1:-1:2
        z[i] = (x[i]-x[i-1])*y[i] + (x[i+1]-x[i-1])*sumy
        sumy += y[i]
    end

    z[1] = (x[2]-x[1])*sumy

    return 0.5*z
end

"""
	diffoperator(du,Δx,D)

Construct differentiation operator for TV method.
"""
function diffoperator(du,Δx,D,epsilon)
    psi = Diagonal(1 ./sqrt.(du.^2 .+ epsilon))
    δx  = Diagonal(Δx)
    return (D')*(δx*psi)*D
end

"""
	differentiate(x::AbstractVector, f::AbstractVector, tv::TotalVariation, α, tol; maxit=2000, pbsize=:auto, verbose=false, precond=false, u0=[0;diff(f)])

Compute the derivative of vector `f` defined at positions `x`, using Total Variation regularisation with parameter `α`. The method is described in Chartrand, R. (2011), "Numerical Differentiation of Noisy, Nonsmooth Data", ISRN Appl. Math., doi:10.5402/2011/164564.

# Arguments
- `x::AbstractVector`: x coordinates where vector f is defined.
- `f::AbstractVector`: input vector to take derivative of.
- `tv::TotalVariation`: instance of TotalVariation type to specify method.
- 'α': regularisation parameter.
- 'tol': tolerance on magnitude of gradient to stop iterations.

# Optional keyword arguments
- `maxit=2000`: maximum number of iterations.
- 'pbsize=:auto': if `:auto` (default), selects the matrix-based computation method only for problems with less than 1001 points. Enforce use of matrix-based method if set to `:small`, and enforce use of matrix-free method is set to `:large`.
- verbose=false. if true, shows convergence info.
- precond=false. if true, use simple preconditioner. Not used for small systems.
- 'u0=[0;diff(f)]': Can be used to set initial guess for derivative.

"""
function differentiate(x::AbstractVector, f::AbstractVector, tv::TotalVariation, α, tol; maxit=2000, pbsize=:auto, verbose=false, precond=false, u0=[0;diff(f)])

    (length(x) == length(f)) || throw(DimensionMismatch("Inputs `x` and `f` should have same dimension"))

    c = 0
    u = u0
    Δx = diff(x)

    if pbsize==:small || (pbsize==:auto && length(f)<1001)
        D = diffmatrix(x)
        A = intgmatrix(x)
        df = A'*(f .- f[1])

        keep = true

        while keep && c<maxit
            c += 1

            # make diff operator
            du = D*u
            L = diffoperator(du,Δx,D,tv.ϵ)
            # gradient
            g = A'*A*u - df + α*L*u
            # approx. hessian
            H = A'*A + α*L
            # QN step
            s = -H\g

            u += s

            keep = (norm(s)/norm(u)>tol)

            !verbose || print("iteration $c | change $(norm(s)/norm(u)) | grad $(norm(g))\r")

        end
    elseif pbsize==:large || (pbsize==:auto && length(f)>=1001)
        D = spdiffmatrix(x)
        df = integrationadjoint(x, f .- f[1])

        keep = true

        while keep && c<maxit
            c += 1

            # make diff operator
            du = D*u
            L = diffoperator(du,Δx,D,tv.ϵ)
            # gradient
            g = integrationadjoint(x, integrationoperator(x,u)) - df + α*L*u

            # approx. hessian
            H = LinearMap(v -> integrationadjoint(x, integrationoperator(x,v)) + α*L*v, length(x))
            if precond
                P = Diagonal(α*diag(L) .+1)
            else
                P = I
            end
            # QN step
            s = cg(H, -g, Pl=P)

            u += s

            keep = (norm(s)/norm(u)>tol)

            !verbose || print("iteration $c | change $(norm(s)/norm(u)) | grad $(norm(g))\r")

        end

    end

    !verbose || println("\nMax. iterations $c")

    return u
end

"""
	differentiate(x::AbstractVector, f::AbstractVector, ::Tikhonov, α; pbsize=:auto)

Compute the derivative of vector `f` defined at positions `x`, using Tikhonov regularisation with parameter `α`. Optional keyword argument `pbsize` can be set to `:small` (uses matrix-based method), `:large` (uses matrix-free method), or `:auto` (default; chooses matrix-based method if `f` has less than 1001 points).

"""
function differentiate(x::AbstractVector, f::AbstractVector, ::Tikhonov, α; pbsize=:auto, maxit=length(x))
    if pbsize==:small || (pbsize==:auto && length(f)<1001)
        D = diffmatrix(x)
        A = intgmatrix(x)
        r = A'*(f .- f[1])
        δ = tr(D'D)/length(x)^3
        H = (α/δ)*D'D + A'A
        u = H\r
    elseif pbsize==:large || (pbsize==:auto && length(f)>=1001)
        D = spdiffmatrix(x)
        δ = tr(D'D)/length(x)^3
        L = (α/δ)*D'D
        r = integrationadjoint(x, f .- f[1])
        H = LinearMap(v -> integrationadjoint(x, integrationoperator(x,v)) + L*v, length(x))
        u, ch = cg(H, r; log=true, maxiter=maxit)

        ch.isconverged || @warn "Solver did not converge after $(ch.iters) iterations. Result might be inaccurate."
    end

    return u
end

end


text:
###########################################
# defines methods to add to point types
# Note the "geotransform" function handles all transformations
# while the "convert" function converts types leaving values unchanged
###########################################

# TODO: This whole thing should probably be switched to a codegen format

# Convenient union types
Geodesy_fam = Union{WorldPosition, WorldSurfacePosition, LocalPosition}  # try to ctach everything
World_fam = Union{WorldPosition, WorldSurfacePosition}
Local_fam = Union{LocalPosition}
LL_fam = Union{LLA, LL}

Proj4_fam = Union{WorldPosition, Proj4Type}      # acceptable types to give to Proj4

Vec3_fam = Union{WorldPosition, LocalPosition}   # for three element point types
Vec2_fam = Union{LL}                             # for three element point types

ELL_param_fam = Union{LLA, LL, ECEF}             # things where an ellipse / psuedo datum are the template param
LL_param_fam = Union{ENU}                        # things where an LLA points is the template param

#########################################
# Accessors
#########################################

# Point translators (lgacy support)
getX{T <: Geodesy_fam}(X::T) = X[1]
getY{T <: Geodesy_fam}(X::T) = X[2]
getZ{T <: Vec3_fam}(X::T) = X[3]

get_lon{T <: Geodesy_fam}(X::T) = X[1]
get_lat{T <: Geodesy_fam}(X::T) = X[2]
get_alt{T <: Vec3_fam}(X::T) = X[3]

get_east{T <: Geodesy_fam}(X::T) = X[1]
get_north{T <: Geodesy_fam}(X::T) = X[2]
get_up{T <: Vec3_fam}(X::T) = X[3]

#####################################################
# Template parameter manipultation
#####################################################

# replace a TypeVar parameter with a DataType parameter where needed

#TODO: Should these be promote rules or something?

#add_param{T <: ELL_param_fam}(::Type{T}) = typeof(T.parameters[1]) == DataType ? T : T{UnknownDatum}  # not type safe :-(
@generated add_param{T <: ELL_param_fam}(::Type{T}) = (typeof(T.parameters[1]) == DataType) ? :(T) :  :(T{UnknownDatum})

#add_param{T <: LL_param_fam}(::Type{T}) = typeof(T.parameters[1]) == DataType ? T : T{UnknownRef}      # not type safe :-(
@generated add_param{T <: LL_param_fam}(::Type{T}) = (typeof(T.parameters[1]) == DataType) ? :(T) :  :(T{UnknownRef})

add_param{T <: CRS}(::Type{T}) = typeof(T.parameters[1]) == DataType ? T : error("always specify an SRID when using the CRS type")

# retrieve datums and ellipsoids
@generated ellipsoid{T <: ELL_param_fam}(::Type{T}) = :($(ellipsoid(T.parameters[1])))           # reference ellipsoid for the position

# methods to pull the reference point from the template
@generated ELL_type{T <: ELL_param_fam}(::Type{T}) = :($(T.parameters[1]))
@generated ELL_type{T <: ELL_param_fam}(::T) = :($(T.parameters[1]))
@generated ELL_type{T <: LL_param_fam}(::Type{T}) = :($(ELL_type(T.parameters[1])))
@generated LL_ref{T <: LL_param_fam}(::Type{T}) = :($(T.parameters[1]))
@generated LL_ref{T <: LL_param_fam}(::T) = :($(T.parameters[1]))

# add the ENU ref or chnage it to LL if its LLA
@generated add_LL_ref{T <: ENU}(::Type{T}) = typeof(T.parameters[1]) == TypeVar ? :(ENU_NULL) : ((typeof(T.parameters[1]) <: LLA) ? :($(ENU{LL(T.parameters[1])})) :  :($(ENU{T.parameters[1]})))

#####################################################
# Force the template parameter when constructing
#####################################################

call(::Type{LL}, x::Real, y::Real) = add_param(LL)(x,y)
call(::Type{LLA}, x::Real, y::Real, z::Real) = add_param(LLA)(x,y,z)
call(::Type{ECEF}, x::Real, y::Real, z::Real) = add_param(ECEF)(x,y,z)
call(::Type{ENU}, x::Real, y::Real, z::Real) = add_param(ENU)(x,y,z)

#
# allow crazy parameter construction for ALL types (legacy support)
#
call{T <: Vec3_fam}(::Type{T}, x::Real, y::Real) = add_param(T)(x,y,0.0)
call{T <: Vec2_fam}(::Type{T}, x::Real, y::Real, z::Real) = add_param(T)(x,y)

####################################
### Proj4 projections
### Try to get the SRID for the point type
####################################

# build methods to get a proj 4 projection (a coordinate reference system) for stuff we know about
get_projection{T <: Union{WorldPosition, WorldSurfacePosition}}(X::T) = get_projection(T)
get_projection{T <: Union{WorldPosition, WorldSurfacePosition}}(::Type{T}) = get_projection(SRID(T))

#####################################################
# Build abstract calls for the transformations
#####################################################

#
# Allow LL <-> LLA
#
geotransform{T}(::Type{LL{T}}, lla::LLA{T}) = LL{T}(lla.lon, lla.lat)
geotransform{T}(::Type{LL}, lla::LLA{T}) = LL{T}(lla.lon, lla.lat)

geotransform{T}(::Type{LLA{T}}, ll::LL{T}) = LLA{T}(lla.lon, lla.lat, 0.0)
geotransform{T}(::Type{LLA}, ll::LL{T}) = LLA{T}(lla.lon, lla.lat, 0.0)

# No conversions, stripping or adding the template parameters
macro default_convs(type_name, known_subtype, null_subtype)
    eval(quote

        # geotransform it to itself (i.e. do nothing)
        geotransform{T}(::Type{$(type_name){T}}, X::$(type_name){T}) = X

        # allow stripping the template parameter from the template
        geotransform{T <: $(known_subtype)}(::Type{$(type_name){$(null_subtype)}}, X::$(type_name){T}) = $(type_name){$(null_subtype)}(X...)

        # allow inserting a template parameter into the template
        geotransform{T <: $(known_subtype)}(::Type{$(type_name){T}}, X::$(type_name){$(null_subtype)}) = $(type_name){T}(X...)

    end)
end

# generate default conversions for everything
@default_convs(LLA, KnownDatum, UnknownDatum)
@default_convs(LL, KnownDatum, UnknownDatum)
@default_convs(ECEF, KnownDatum, UnknownDatum)
@default_convs(ENU, LLA, UnknownRef)

# dont macro srid type's, we don't want to allow stripping of the SRID
geotransform{T <: SRID}(::Type{CRS}, X::CRS{T}) = X
geotransform{T <: SRID}(::Type{CRS{T}}, X::CRS{T}) = X

#
# templates for constructing one point type from another (I want to leave convert free to do value preserving conversions)
#
macro add_cross_const(Type1, Type2)
    eval(quote

        # copy constructor for the first one
        call{T <: $(Type1), U <: $(Type1)}(::Type{T}, X::U) = geotransform(T, X)

        # and the cross versions
        call{T <: $(Type1), U <: $(Type2)}(::Type{T}, X::U) = geotransform(T, X)
        call{T <: $(Type2), U <: $(Type1)}(::Type{T}, X::U) = geotransform(T, X)
    end)
end
@add_cross_const(WorldPosition, LocalPosition)
@add_cross_const(LocalPosition, WorldSurfacePosition)
@add_cross_const(WorldSurfacePosition, WorldPosition)

#
# constructing one point type from another and a reference
#
call{T <: Local_fam}(::Type{T}, X::World_fam, ll_ref::LL_fam) = geotransform(T, X, ll_ref)
call{T <: World_fam}(::Type{T}, X::Local_fam, ll_ref::LL_fam) = geotransform(T, X, ll_ref)

#
# allow construction from a matrix via convert as its value preserving
#
function convert{T <: Geodesy_fam}(::Type{Vector{T}}, X::AbstractMatrix; row::Bool=true)
    oT = add_param(T)
    n = (row) ? size(X,1) : size(X,2)
    Xout = Vector{oT}(n)  # cant make list comprehesion get the output type right
    if (T <: Vec2_fam) && row
        for i = 1:n; Xout[i] = oT(X[i,1], X[i,2]); end
    elseif (row)
        for i = 1:n; Xout[i] = oT(X[i,1], X[i,2], X[i,3]); end
    elseif (T <: Vec2_fam)
        for i = 1:n; Xout[i] = oT(X[1,i], X[2,i]); end
    else
        for i = 1:n; Xout[i] = oT(X[1,i], X[2,i], X[3,i]); end
    end
    return Xout
end

###################################################################
### Methods to add template parameters based on other arguments ###
###################################################################

# we want to make sure any created Vector have the template parameter in them
# add_param{T <: ELL_param_fam, U <: ELL_param_fam}(::Type{T}, ::Type{U}) = T ==  add_param(T) ? T : T{U.parameters[1]}                              # not type safe :-(
@generated add_param{T <: ELL_param_fam, U <: ELL_param_fam}(::Type{T}, ::Type{U}) = (T == add_param(T)) ? :(T) : :(T{$(ELL_type(U))})

# add_param{T <: ELL_param_fam, U <: LL_param_fam}(::Type{T}, ::Type{U}) = T ==  add_param(T) ? T : T{typeof(U.parameters[1]).parameters[1]}        # not type safe :-(
@generated add_param{T <: ELL_param_fam, U <: LL_param_fam}(::Type{T}, ::Type{U}) = T ==  add_param(T) ? :(T) : :(T{$(typeof(U.parameters[1]).parameters[1])})

# default to not include reference position in local types
add_param{T <: LL_param_fam, U <: Geodesy_fam}(::Type{T}, ::Type{U}) = add_param(T)

# add parameters when it might be an SRID
@generated add_param{T <: Geodesy_fam, U <: Geodesy_fam}(::Type{T}, ::Type{U}) = (T <: CRS) ? :(T) : ((U <: CRS) ? :(add_param(T)) : :(add_param(T, U)))

####################################
### Transformations of Vectors   ###
####################################

# a vectorized way to perform transformations
function geotransform{T <: Geodesy_fam, U <: Geodesy_fam}(::Type{T}, X::Vector{U})

    # make sure the output parameter is filled
    oT =  add_param(T, U)

    # is Proj4 involved?
    if (T <: CRS) || (U <: CRS)
        Xout = proj4_vectorized(oT, X)
    else
        Xout = Vector{oT}(length(X))
        @inbounds for i = 1:length(X)
            Xout[i] = geotransform(oT, X[i])
        end
    end
    return Xout
end

# a vectorized way to perform transformations with reference points
function geotransform{T <: Geodesy_fam, U <: Geodesy_fam, V <: LL_fam}(::Type{T}, X::Vector{U}, ll_ref::V)

    # get inputs and desired output types
    oT = add_param(T, V)

    Xout = Vector{oT}(length(X))
    @inbounds for i = 1:length(X)
        Xout[i] = geotransform(oT, X[i], ll_ref)
    end
    return Xout
end

# worker code for proj4 vectorization (looping prj4 is slow)
function proj4_vectorized{T <: Proj4_fam, U <: Proj4_fam}(::Type{T}, X::Vector{U})

    if !((T <: CRS) || (U <: CRS))
        warn("Unexpected: using Proj4 to geotransform between Geodesy point types.  How'd this happen")
    end

    # convert to a matrix
    mat = Matrix{Float64}(length(X), 3)
    @inbounds for i = 1:length(X)
        mat[i, 1] = X[i][1]
        mat[i, 2] = X[i][2]
        mat[i, 3] = U <: Vec2_fam ? 0.0 : X[i][3]
    end

    # perform it
    Proj4.transform!(Geodesy.get_projection(U), Geodesy.get_projection(T), mat, false)

    # and assign the output
    X = Vector{T}(length(X))
    @inbounds for i = 1:length(X)
        X[i] = T(mat[i,1], mat[i,2], mat[i,3])
    end

    return X
end


text:
module CrossModule

using BinaryTraits
using Test, Logging

module X
    using BinaryTraits
    using BinaryTraits.Prefix: Is
    @trait RowTable
    @implement Is{RowTable} by row(_, ::Integer)
    __init__() = init_traits(@__MODULE__)
end

module Y
    using Test
    using BinaryTraits
    using BinaryTraits.Prefix: Is
    using ..X
    struct AwesomeTable end
    @assign AwesomeTable with Is{X.RowTable}
    r = @check(AwesomeTable)
    @test r.implemented |> length == 0
    X.row(::AwesomeTable, ::Number) = 1
    __init__() = init_traits(@__MODULE__)
end

function test()
        r = @check(Y.AwesomeTable)
        @test r.implemented |> length == 1
end

end # module

using Test
@testset "Cross Module" begin
    import .CrossModule
    CrossModule.test()
end


text:
include("layers.jl")
include("initializers.jl")

"""
A custom constructor for Conv2D + BatchNorm + Activation layer
    w1          : First dimension of the kernel
    w2          : Second dimension of the kernel
    input_dim   : Channel dimension of the input
    output_dim  : Filter size
    f           : For leaky relu pass "leaky_relu"
    alpha       : Activation constant for leaky relu, please leave as default or set to 0 for relu
    pdrop       : Dropout ratio
    padding     : Padding size
    stride      : Stride size
    dilation    : Dilation size
"""
mutable struct ConvBn conv; bn; f; alpha; end

function ConvBn(w1::Int, w2::Int, input_dim::Int, output_dim::Int; init=kaiming_uniform, f=nothing, alpha=0,
    pdrop=0, bias=true, padding=0, stride=1, dilation=1, dtype=Array{Float32}, momentum=0.1)
    return ConvBn(
        Conv2D(w1, w2, input_dim, output_dim, pdrop=pdrop, dtype=dtype, init=init,
                padding=padding, stride=stride, dilation=dilation, bias=bias),
        BatchNorm(momentum=momentum, channels=output_dim, dtype=dtype),
        f, alpha
    )
end

function (cbr::ConvBn)(x; train=true)
    x_val = cbr.conv(x, train=train)
    # print("Sum Conv: ", sum(Array(value(abs.(x_val)))) / prod(size(x_val)), "\n")
    x_val = cbr.bn(x_val, train=train)
    # print("Sum BN: ", sum(Array(value(abs.(x_val)))) / prod(size(x_val)), "\n")
    if cbr.f === nothing return x_val
    elseif cbr.f == "leaky_relu" return max.(cbr.alpha .* x_val, x_val)
    else return cbr.f.(x_val)
    end
end

"""
Residual Layer structure with kernel sizes 1x1, 3x3, 1x1 in order. Used mostly in ResNet
Network with more than 50 layers.
"""
mutable struct Residual_1x3x1  downsample; conv_bn1; conv_bn2; conv_bn3; dtype; end

function Residual_1x3x1(input_dim, filter_sizes; downsample=false, ds_3x3_stride=1, init=kaiming_uniform,
                        bias=false, pdrop=0, dtype=Array{Float32})
    ds_layer = nothing
    if downsample
        ds_layer = ConvBn(1, 1, input_dim, filter_sizes[3], bias=bias, init=init, dtype=dtype, stride=ds_3x3_stride)
    end

    return Residual_1x3x1(
        ds_layer,
        ConvBn(1, 1, input_dim,       filter_sizes[1], init=init, bias=bias, dtype=dtype, pdrop=pdrop, f=relu),
        ConvBn(3, 3, filter_sizes[1], filter_sizes[2], init=init, bias=bias, dtype=dtype, pdrop=pdrop, f=relu, padding=1, stride=ds_3x3_stride),
        ConvBn(1, 1, filter_sizes[2], filter_sizes[3], init=init, bias=bias, dtype=dtype, pdrop=pdrop),
        dtype
    )
end

function (r::Residual_1x3x1)(x; train=true)
    x_val = r.conv_bn1(x, train=train)
    x_val = r.conv_bn2(x_val, train=train)
    x_val = r.conv_bn3(x_val, train=train)
    if r.downsample === nothing x_val += x
    else x_val += r.downsample(x, train=train)
    end
    return relu.(x_val)
end

function load_mat_weights(block::Residual_1x3x1, conv_w, bn_mom, bn_b, bn_mult)
    idx = 1; mom_idx = 1;

    if block.downsample !== nothing
        block.downsample.conv.w = Param(conv_w[idx])
        block.downsample.bn.bn_moments = bnmoments(mean=bn_mom[mom_idx], var=bn_mom[mom_idx+1])
        block.downsample.bn.bn_params = Param(vcat(vec(bn_mult[idx]), vec(bn_b[idx])))
        idx += 1; mom_idx += 2;
    end

    block.conv_bn1.conv.w = Param(conv_w[idx])
    block.conv_bn1.bn.bn_moments = bnmoments(mean=bn_mom[mom_idx], var=bn_mom[mom_idx+1])
    block.conv_bn1.bn.bn_params = Param(vcat(vec(bn_mult[idx]), vec(bn_b[idx])))
    idx += 1; mom_idx += 2;

    block.conv_bn2.conv.w = Param(conv_w[idx])
    block.conv_bn2.bn.bn_moments = bnmoments(mean=bn_mom[mom_idx], var=bn_mom[mom_idx+1])
    block.conv_bn2.bn.bn_params = Param(vcat(vec(bn_mult[idx]), vec(bn_b[idx])))
    idx += 1; mom_idx += 2;

    block.conv_bn3.conv.w = Param(conv_w[idx])
    block.conv_bn3.bn.bn_moments = bnmoments(mean=bn_mom[mom_idx], var=bn_mom[mom_idx+1])
    block.conv_bn3.bn.bn_params = Param(vcat(vec(bn_mult[idx]), vec(bn_b[idx])))

    return block
end

"""
To add layers on top of each other.
"""
mutable struct Chain; layers; end
(c::Chain)(x; train=true) = (for l in c.layers; x = l(x, train=train); end; x)

max_stars_repo_name: barisbatuhan/retinaface

max_issues_repo_name: barisbatuhan/retinaface

max_forks_repo_name: barisbatuhan/retinaface


text:
using QuantumStateDB
using LibPQ
using QuantumStateBase
using QuantumStateDistributions
using DataFrames
using Test

load_config!()
const PG = current_dbconfig()
@show PG[:user], PG[:password], PG[:dbname]

function drop_db(dbname::String)
    connection = LibPQ.Connection(QuantumStateDB.to_config_string(PG))
        @info "Drop database $dbname if exists!"
        execute(connection, "DROP DATABASE IF EXISTS $dbname;")
    close(connection)

    return 0
end

const DBNAME = "test_db"

@testset "QuantumStateDB.jl" begin
    include("utils.jl")
    include("db.jl")
    include("table.jl")
    include("api.jl")
end


text:
using SIMD

abstract type AbstractDataStream end

mutable struct DataBuffer{N,T}
    head ::Union{Nothing, Vec{N, T}}
    tail ::AbstractArray{T}
    DataBuffer(block_size::Val{N}, data::AbstractArray{T}) where {N,T} = new{N,T}(Vec(tuple(data[1:N]...)), (@view data[N+1:end]))
end

mutable struct MergeNode{N,T}
    head ::Union{Nothing, Vec{N, T}}
    left ::Union{MergeNode{N,T}, DataBuffer{N,T}}
    right ::Union{MergeNode{N,T}, DataBuffer{N,T}}
    MergeNode(left ::Union{MergeNode{N,T}, DataBuffer{N,T}}, right ::Union{MergeNode{N,T}, DataBuffer{N,T}}) where {N,T} =
        new{N,T}(pop!(first_stream(left, right)), left, right)
end

function build_multi_merger(block_size, data...)
    k = length(data)
    if k == 1
        DataBuffer(block_size, data[1])
    else
        MergeNode(build_multi_merger(block_size, data[1:div(k,2)]...),
                  build_multi_merger(block_size, data[1+div(k,2):end]...))
    end
end

"""
Returns the stream with the smallest first element. If the two streams are empty, returns `nothing`.
"""
function first_stream(left ::Union{MergeNode{N,T}, DataBuffer{N,T}}, right ::Union{MergeNode{N,T}, DataBuffer{N,T}}) where {N,T}
    if left.head == nothing && right.head == nothing
        nothing
    elseif left.head == nothing && right.head != nothing
        right
    elseif left.head != nothing && right.head == nothing
        left
    elseif left.head[1] < right.head[1]
        left
    else
        right
    end
end

function pop!(dbuf::DataBuffer{N}) where N
    output = dbuf.head
    new_head = @view dbuf.tail[1:min(N, length(dbuf.tail))]
    dbuf.head = if length(new_head)>0 vload(Vec{N, eltype(dbuf.tail)}, dbuf.tail, 1) else nothing end
    dbuf.tail = @view dbuf.tail[(N+1):end]
    output
end

function pop!(node::MergeNode)
    smallest_tail = first_stream(node.left, node.right)

    if smallest_tail == nothing
        output = node.head
        node.head = nothing
        output
    else
        input = pop!(smallest_tail)
        output, node.head = bitonic_merge(input, node.head)
        output
    end
end