text:
# Autogenerated wrapper script for libCEED_jll for x86_64-apple-darwin
export libceed

JLLWrappers.@generate_wrapper_header("libCEED")
JLLWrappers.@declare_library_product(libceed, "@rpath/libceed.dylib")
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_library_product(
        libceed,
        "lib/libceed.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
function NodeDetector(arg0::Orbit, arg1::Frame)
    return NodeDetector((Orbit, Frame), arg0, arg1)
end

function NodeDetector(arg0::jdouble, arg1::Orbit, arg2::Frame)
    return NodeDetector((jdouble, Orbit, Frame), arg0, arg1, arg2)
end

function g(obj::NodeDetector, arg0::SpacecraftState)
    return jcall(obj, "g", jdouble, (SpacecraftState,), arg0)
end

function get_frame(obj::NodeDetector)
    return jcall(obj, "getFrame", Frame, ())
end

ext: jl


text:
using SpinGlassNetworks
using LabelledGraphs
using LightGraphs
using MetaGraphs
using Logging
using Test

Base.:(==)(e1::LabelledEdge, e2::LabelledEdge) = src(e1) == src(e2) && dst(e1) == dst(e2)

my_tests = ["ising.jl", "factor.jl"]

for my_test ∈ my_tests
    include(my_test)
end


text:
"""
    ClimaParams

This is a module for staging experimental
clima parameters-- that we are unsure will
be used in the final model. Once they are
no longer experimental, they should be moved
to CLIMAParameters.jl.

TODO: move non-experimental parameters to CLIMAParameters
"""
module ClimaParams

import CLIMAParameters
const CP = CLIMAParameters
const APS = CP.AbstractEarthParameterSet

""" divergence factor for bubble case (zero otherwise) """
entrainment_massflux_div_factor(ps::APS) = ps.nt.c_div
entrainment_sigma(ps::APS) = ps.nt.μ
area_limiter_scale(ps::APS) = ps.nt.γ_lim
area_limiter_power(ps::APS) = ps.nt.β_lim
static_stab_coeff(ps::APS) = ps.nt.c_b
l_max(ps::APS) = ps.nt.l_max

"""Additional entrainment/detrainment parameters."""
c_gen(ps::APS) = ps.nt.c_gen

""" stochastic parameters """
stoch_ε_lognormal_var(ps::APS) = ps.nt.stoch_ε_lognormal_var
stoch_δ_lognormal_var(ps::APS) = ps.nt.stoch_δ_lognormal_var

sde_ϵ_θ(ps::APS) = ps.nt.sde_ϵ_θ
sde_ϵ_σ(ps::APS) = ps.nt.sde_ϵ_σ
sde_δ_θ(ps::APS) = ps.nt.sde_δ_θ
sde_δ_σ(ps::APS) = ps.nt.sde_δ_σ

end


text:
push!(LOAD_PATH, "/Users/megorov/Desktop/projects/stanford/sniper/egorov/sniperRepo/new-version")

using Collada_
using UrbanMap_
using SniperPOMDP_
using SniperMDP_
using ValueIterations
using Helpers
using RayCaster_
using Simulator
using SimServers
using PolicyAssignment

mapPath = "/Users/megorov/Desktop/projects/stanford/sniper/egorov/sniperRepo/new-version/maps/demo_map_border2D.dae"
coll = ColladaObjects(mapPath)
map = UrbanMap(coll, 15, 15) # 15x15 map

t   = [(7.0,7.0), (2.0,14.0), (8.0,4.0), (11.0,13.0)]
mdp = SniperMDP(map, nSnipers = 1, nMonitors = 4, targets=t)

p1 = load_policy("policies/demo-mdp-1.policy")
p2 = load_policy("policies/demo-mdp-2.policy")
p3 = load_policy("policies/demo-mdp-3.policy")
p4 = load_policy("policies/demo-mdp-4.policy")

policies = merge(p1.policy,p2.policy,p3.policy,p4.policy)
values   = merge(p1.V,p2.V,p3.V,p4.V)

policy = AssignmentPolicy(4, policies, values, mdp)

# create the simulation type
#init  = 8_888_888_803 # sniper at (3,1), all monitors at (9,9) for 10x10 4 agent world
#init = 97 # sniper at (7,7), all monitors at (1,1) for 15x15 4 agent world
init = 20_583_164_347 # sniper at (7,7), all monitors at (1,i) for 15x15 4 agent world
sim  = MDPSimulation(mdp, init, 100)

# set up the server
socket = 4444
server = SniperServer(socket)

# begin simulation and send results to client
simAndSend(server, sim, mdp, policy)


text:
using NetworkEpidemics
using Test

@testset "NetworkEpidemics.jl" begin
    include("categorical_tree.jl")
end


text:
fname = :ARCLDLt_abs
c = Combi(hessian_dense,PDataLDLt,solve_modelARCDiagAbs,preprocessLDLt,decreaseFact,Tparam())
include("Template.jl")


text:
struct ReflectedShader
    ir::IR
end

reflect(shader::ShaderFile{FormatSPIRV}) = ReflectedShader(IR(SPIRV.Module(shader.file)))


text:
# Simulations are for running

function stop(sim)
    time_before = time()

    for sc in sim.stop_criteria
        if sc(sim)
            time_after = time()
            sim.run_time += time_after - time_before
            return true
        end
    end

    time_after = time()
    sim.run_time += time_after - time_before

    return false
end

function iteration_limit_exceeded(sim)
    if sim.model.clock.iteration >= sim.stop_iteration
          @info "Simulation is stopping. Model iteration $(sim.model.clock.iteration) " *
                "has hit or exceeded simulation stop iteration $(sim.stop_iteration)."
          return true
    end
    return false
end

function stop_time_exceeded(sim)
    if sim.model.clock.time >= sim.stop_time
          @info "Simulation is stopping. Model time $(prettytime(sim.model.clock.time)) " *
                "has hit or exceeded simulation stop time $(prettytime(sim.stop_time))."
          return true
    end
    return false
end

function wall_time_limit_exceeded(sim)
    if sim.run_time >= sim.wall_time_limit
          @info "Simulation is stopping. Simulation run time $(prettytime(sim.run_time)) " *
                "has hit or exceeded simulation wall time limit $(prettytime(sim.wall_time_limit))."
          return true
    end
    return false
end

get_Δt(Δt) = Δt
get_Δt(wizard::TimeStepWizard) = wizard.Δt
get_Δt(simulation::Simulation) = get_Δt(simulation.Δt)

"""
    run!(simulation)

Run a `simulation` until one of the stop criteria evaluates to true. The simulation
will then stop.
"""
function run!(sim)
    model = sim.model
    clock = model.clock

    [open(out) for out in values(sim.output_writers)]

    while !stop(sim)
        time_before = time()

        if clock.iteration == 0
            [run_diagnostic(sim.model, diag) for diag in values(sim.diagnostics)]
            [write_output(sim.model, out)    for out  in values(sim.output_writers)]
        end

        for n in 1:sim.progress_frequency
            euler = clock.iteration == 0 || (sim.Δt isa TimeStepWizard && n == 1)
            time_step!(model, get_Δt(sim.Δt), euler=euler)

            [time_to_run(clock, diag) && run_diagnostic(sim.model, diag) for diag in values(sim.diagnostics)]
            [time_to_run(clock, out)  && write_output(sim.model, out)    for out  in values(sim.output_writers)]
        end

        sim.progress isa Function && sim.progress(sim)
        sim.Δt isa TimeStepWizard && update_Δt!(sim.Δt, model)

        time_after = time()
        sim.run_time += time_after - time_before
    end

    [close(out) for out in values(sim.output_writers)]

    return nothing
end


text:
@time using DamagedShearBand
const DSB = DamagedShearBand
@time using Plots
using Printf
##
# imposed values
ϵ̇₁₁ = -1e-5
ϵ̇ⁱξη = 1e-5
σ₃ = -1e6 # negative in compression

# initial higher compressive stress multiplier (S = σ₁₁/σ₂₂)
S_i = 1.0
θ = 60.0

# Create a Rheology type instance containing elastic moduli and damage parameters. Change default values by  supplying keywords arguments.
r = DSB.Rheology(D₀=0.2,n=5.0)

# initial damage
#D_i = r.D₀
Dᵒ, _ = DSB.get_KI_mininizer_D_on_S_range(r,S_i,100,σ₃)
Dⁱ = min(Dᵒ+0.15, 0.9)

# time integration parameters :
ϵⁱξη_goal = 0.01
tspan = (0.0,ϵⁱξη_goal/ϵ̇ⁱξη) # initial and final simulation time
#tspan = (0.0,500)
Δt_i = 0.01 # initial timestep

# Build stress and strain tensors
# σᵢⱼ_i = DSB.build_principal_stress_tensor(r,S_i,σ₃,D_i ; abstol=1e-15) # takes care of the plane strain constraint by solving non linear out of plane strain wrt σ₃₃ using Newton algorithm
# ϵᵢⱼ_i = DSB.compute_ϵij(r,D_i,σᵢⱼ_i)

n_iter_saved = 20000
n_iter_printed = 10
op = DSB.OutputParams(save_frequency = n_iter_saved/(tspan[2]-tspan[1]),
                      save_period = 1000,
                      print_frequency = n_iter_printed/(tspan[2]-tspan[1]),
                      print_period = 1000)

sp = DSB.SolverParams(newton_abstol = 1e-12,
                      newton_maxiter = 100,
                      time_maxiter = nothing,
                      e₀ = (D=1e-2, S=1e-2, σ=100.0, ϵ=1e-2))
# sp = DSB.SolverParams(newton_abstol = 1e-12,
#                       newton_maxiter = 100,
#                       time_maxiter = nothing,
#                       e₀ = (D=1e-3, S=1e-3, σ=1.0, ϵ=1e-5))

p = DSB.Params(op,sp)
# Integrate over time, by solving for the unknowns [σ₁₁next, σ₃₃next, ϵ₂₂next] at each timestep
# -- e₀ is the absolute tolerance used to adapt time stepping, may be a scalar or a NamedTuple with keys (D,σ,ϵ)
# -- abstol is the absolute tolerance on the unknowns in the Newton solver

t_vec, S_vec, σⁱᵢⱼ_vec, σᵒᵢⱼ_vec, ϵⁱᵢⱼ_vec, Dⁱ_vec, Dᵒ_vec = DSB.adaptative_time_integration_2_points(r,p,S_i,σ₃,Dⁱ,Dᵒ,ϵ̇₁₁,ϵ̇ⁱξη,Δt_i,θ,tspan);

#
σᵒ_p_vec = DSB.principal_coords.(σᵒᵢⱼ_vec,Ref(θ))
σⁱ_p_vec = DSB.principal_coords.(σⁱᵢⱼ_vec,Ref(θ))
σᵒ₁₁_p_vec = [σᵢⱼ[1,1] for σᵢⱼ in σᵒ_p_vec]
σᵒ₂₂_p_vec = [σᵢⱼ[2,2] for σᵢⱼ in σᵒ_p_vec]
σᵒ₃₃_p_vec = [σᵢⱼ[3,3] for σᵢⱼ in σᵒ_p_vec]

σᵒ₁₁_vec = [σᵢⱼ[1,1] for σᵢⱼ in σᵒᵢⱼ_vec]
σⁱ₁₁_vec = [σᵢⱼ[1,1] for σᵢⱼ in σⁱᵢⱼ_vec]
σᵒ₁₂_vec = [σᵢⱼ[1,2] for σᵢⱼ in σᵒᵢⱼ_vec]
σⁱ₁₂_vec = [σᵢⱼ[1,2] for σᵢⱼ in σⁱᵢⱼ_vec]
σᵒ₂₁_vec = [σᵢⱼ[2,1] for σᵢⱼ in σᵒᵢⱼ_vec]
σⁱ₂₁_vec = [σᵢⱼ[2,1] for σᵢⱼ in σⁱᵢⱼ_vec]
σᵒ₂₂_vec = [σᵢⱼ[2,2] for σᵢⱼ in σᵒᵢⱼ_vec]
σⁱ₂₂_vec = [σᵢⱼ[2,2] for σᵢⱼ in σⁱᵢⱼ_vec]
σᵒ₃₃_vec = [σᵢⱼ[3,3] for σᵢⱼ in σᵒᵢⱼ_vec]
σⁱ₃₃_vec = [σᵢⱼ[3,3] for σᵢⱼ in σⁱᵢⱼ_vec]

ϵⁱ₁₁_vec = [ϵᵢⱼ[1,1] for ϵᵢⱼ in ϵⁱᵢⱼ_vec]
ϵⁱ₁₂_vec = [ϵᵢⱼ[1,2] for ϵᵢⱼ in ϵⁱᵢⱼ_vec]
ϵⁱ₂₁_vec = [ϵᵢⱼ[2,1] for ϵᵢⱼ in ϵⁱᵢⱼ_vec]
ϵⁱ₂₂_vec = [ϵᵢⱼ[2,2] for ϵᵢⱼ in ϵⁱᵢⱼ_vec]
ϵⁱ₃₃_vec = [ϵᵢⱼ[3,3] for ϵᵢⱼ in ϵⁱᵢⱼ_vec]

τⁱ_vec =  DSB.get_τ.(σⁱᵢⱼ_vec; input=:σ)
τᵒ_vec =  DSB.get_τ.(σᵒᵢⱼ_vec; input=:σ)
pⁱ_vec = -1/3 * tr.(σⁱᵢⱼ_vec)
pᵒ_vec = -1/3 * tr.(σⁱᵢⱼ_vec)
KIⁱ_vec = [DSB.compute_KI(r,σⁱᵢⱼ_vec[i], Dⁱ_vec[i]) for i in eachindex(σⁱᵢⱼ_vec)]
KIᵒ_vec = [DSB.compute_KI(r,σᵒᵢⱼ_vec[i], Dᵒ_vec[i]) for i in eachindex(σᵒᵢⱼ_vec)]
γⁱ_vec = [ sqrt(2*ϵᵢⱼ⊡ϵᵢⱼ) for ϵᵢⱼ in ϵⁱᵢⱼ_vec]
ϵⁱkk_vec = [tr(ϵᵢⱼ) for ϵᵢⱼ in ϵⁱᵢⱼ_vec]
deviation_out_vec = DSB.get_stress_deviation_from_far_field.(σᵒ_p_vec ; offdiagtol=1e-5)
deviation_in_vec = DSB.get_stress_deviation_from_far_field.(σⁱ_p_vec ; offdiagtol=1e-5)

σᵒᵢⱼ_i_p, σⁱᵢⱼ_i_p, ϵᵒᵢⱼ_i_p, ϵⁱᵢⱼ_i_p = DSB.initialize_state_var_D(r,p,S_i,σ₃,Dⁱ,Dᵒ,θ ; coords=:principal)
σᵒᵢⱼ_i, σⁱᵢⱼ_i, ϵᵒᵢⱼ_i, ϵⁱᵢⱼ_i = DSB.initialize_state_var_D(r,p,S_i,σ₃,Dⁱ,Dᵒ,θ ; coords=:band)
deviation_out_check = DSB.get_stress_deviation_from_far_field(σᵒᵢⱼ_i_p ; offdiagtol=1e-12)
deviation_in_check = DSB.get_stress_deviation_from_far_field(σⁱᵢⱼ_i_p ; offdiagtol=1e-12)
deviation_out_eps = DSB.get_stress_deviation_from_far_field(ϵᵒᵢⱼ_i_p ; offdiagtol=1e-12)
deviation_in_eps = DSB.get_stress_deviation_from_far_field(ϵⁱᵢⱼ_i_p ; offdiagtol=1e-12)

deviation_out_check = DSB.get_stress_deviation_from_far_field(DSB.principal_coords(σᵒᵢⱼ_i,θ))
deviation_in_check = DSB.get_stress_deviation_from_far_field(DSB.principal_coords(σⁱᵢⱼ_i,θ))
deviation_out_eps = DSB.get_stress_deviation_from_far_field(DSB.principal_coords(ϵᵒᵢⱼ_i,θ))
deviation_in_eps = DSB.get_stress_deviation_from_far_field(DSB.principal_coords(ϵⁱᵢⱼ_i,θ))
## Plots ##
# figure size and backend
plotlyjs(size=(1000,300))

# S vs time
plot(t_vec,γⁱ_vec)
plot(t_vec,ϵⁱ₁₂_vec)
plot(t_vec,S_vec)
plot(ϵⁱ₁₂_vec,S_vec)
# principal stress orientation deviation from far field vs time
plot(t_vec,deviation_in_vec)
plot!(t_vec,deviation_out_vec)
#plot!(t_vec,-(Dⁱ_vec .- Dⁱ).*60)

plot(deviation_in_vec)
# Damage vs time
plot(t_vec,Dⁱ_vec)
plot!(t_vec,Dᵒ_vec)
# KI/KIC vs time
plot(t_vec,KIⁱ_vec./r.K₁c)
plot!(t_vec,KIᵒ_vec./r.K₁c)
# shear stress on band vs time
plot(t_vec,σⁱ₁₂_vec)
plot!(t_vec,σᵒ₁₂_vec)
# vol/shear strain invariants vs time
plot(t_vec, ϵⁱkk_vec./γⁱ_vec)
plot(t_vec[1:end-1], diff(ϵⁱkk_vec)./diff(γⁱ_vec))

# plot τ in and out
plot(t_vec, τⁱ_vec)
plot!(t_vec, τᵒ_vec)
## composite plot :
#lotlyjs(size=(1000,600))
pyplot(size=(1000,600))
layout = @layout [a ; b ; c]

p1 = plot(ϵⁱ₁₂_vec.*100,S_vec,label="S")
title!(p1,"D₀ = $(r.D₀) & D_in/out = ("*@sprintf("%.2f",Dⁱ)*", "*@sprintf("%.2f",Dᵒ)*") & gammadot = $(2*ϵ̇ⁱξη) & p=$(-σ₃)")
ylabel!("σᵒᵘᵗ₁ multiplier (S)")

p2 = plot(ϵⁱ₁₂_vec.*100,Dⁱ_vec,label="inside")
plot!(p2,γⁱ_vec.*100,Dᵒ_vec,label="outside")
ylabel!("Damage")

p3 = plot(ϵⁱ₁₂_vec.*100,deviation_in_vec,label="inside")
plot!(p3,ϵⁱ₁₂_vec.*100,deviation_out_vec,label="outside")
xlabel!("ϵⁱ₁₂ (%)")
ylabel!("stress orientation in band")

pl = plot(p1,p2,p3,layout=layout)

#savefig(pl,"./examples/figures/bifurcation_D0=$(r.D₀)_p=$(-σ₃)_ϵ̇ⁱξη=$(ϵ̇ⁱξη).pdf")


text:
###############################################################################
#
# Xylem K ratio functions
#
###############################################################################
"""
    xylem_k_ratio(
                vc::AbstractXylemVC{FT},
                p_25::FT,
                vis::FT
    ) where {FT<:AbstractFloat}

Returns the relative hydraulic conductance, given
- `vc` Xylem vulnerability curve
- `p` Xylem pressure at 298.15 K in `[MPa]`
- `p_25` Equivalent xylem pressure at 298.15 K in `[MPa]`
- `vis` Relative viscosity. If missing, vis = 1.
"""
function xylem_k_ratio(
            vc::WeibullSingle{FT},
            p_25::FT,
            vis::FT = FT(1)
) where {FT<:AbstractFloat}
    @unpack b,c = vc;

    if p_25<0
        kr = max( FT(1e-4), exp( -1 * (-p_25 / vc.b) ^ vc.c ) ) / vis;
    else
        kr = 1 / vis;
    end

    return kr
end

function xylem_k_ratio(
            vc::WeibullDual{FT},
            p_25::FT,
            vis::FT = FT(1)
) where {FT<:AbstractFloat}
    @unpack b1,c1,f1,b2,c2,f2 = vc;

    if p_25<0
        k1 = exp( -1 * (-p_25 / b1) ^ c1 ) * f1;
        k2 = exp( -1 * (-p_25 / b2) ^ c2 ) * f2;
        kr = max( FT(1e-4), (k1+k2) ) / vis;
    else
        kr = 1 / vis;
    end

    return kr
end

###############################################################################
#
# Xylem P_crit functions
#
###############################################################################
"""
    xylem_p_crit(vc::AbstractXylemVC{FT}, f_st::FT) where {FT<:AbstractFloat}

Returns the relative hydraulic conductance, given
- `vc` Xylem vulnerability curve
- `st` Relative surface tension. If missing, vis = 1.
"""
function xylem_p_crit(
            vc::WeibullSingle{FT},
            f_st::FT = FT(1)
) where {FT<:AbstractFloat}
    @unpack b,c = vc;

    return -b * log( FT(1000) ) ^ (1 / c) * f_st
end

function xylem_p_crit(
            vc::WeibullDual{FT},
            f_st::FT = FT(1)
) where {FT<:AbstractFloat}
    @unpack b1,c1,b2,c2 = vc;

    p1 = -b1 * log( FT(1000) ) ^ (1 / c1) * f_st
    p2 = -b2 * log( FT(1000) ) ^ (1 / c2) * f_st

    return min(p1, p2)
end

###############################################################################
#
# Fit xylem vulnerability curve
#
###############################################################################

function fit_xylem_VC(xs::Array{FT,1}, ys::Array{FT,1}; label="TPLC") where {FT<:AbstractFloat}
    _ps = xs;
    _ks = ys;

    # if x is tension and y is plc
    if label=="TPLC"
        _ps = -xs;
        _ks = (100 .- ys) ./ 100;
    end

    @inline f(x) = (
        @show x;
        _vc = WeibullSingle{FT}(b=x[1], c=x[2]);
        _kp = xylem_k_ratio.([_vc], _ps);
        return -sum((_kp .- _ks) .^ 2);
    );

    _st = SolutionToleranceND{FT}([1e-3, 1e-3], 30);
    _ms = ReduceStepMethodND{FT}(x_mins = FT[1e-3, 1e-3],
                                 x_maxs = FT[ 100, 100],
                                 x_inis = [1, 1],
                                 Δ_inis = FT[1, 1]);
    _bc = find_peak(f, _ms, _st);

    @show _bc;

    return _bc
end

max_stars_repo_name: bbuman/Land

max_issues_repo_name: bbuman/Land

max_forks_repo_name: bbuman/Land


text:
using Bonobo
using Test

using JuMP
using HiGHS

const BB = Bonobo

@testset "End2End tests" begin
    include("end2end/dummy.jl")
    include("end2end/mip.jl")
end

@testset "Unit tests" begin
    include("unit/add_node.jl")
end


text:
using MatrixNetworks
using DelimitedFiles
A = MatrixNetworks.readSMAT("commutes-all.smat")
xy = readdlm("commutes-all.xy")
##
using SparseArrays
i = rand(1:size(A,1))
scatter(xy[:,1],xy[:,2],markersize=1, markerstrokewidth=0)
neighs,vals = findnz(A[:,i])
println(A[:,i].nzval)
scatter!([xy[i,1]], [xy[i,2]], markersize=9)
scatter!(xy[neighs,1],xy[neighs,2], markerstrokewidth=0, markersize=5, marker_z = vals)

##
##
using Plots
using GraphRecipes
using SparseArrays
src,dst = findnz(A)[1:2]
graphplot(src,dst,x = xy[:,1], y = xy[:,2],
  markersize=0, linewidth=0.5, linealpha=0.1)

max_stars_repo_name: empet/more-graphs

max_issues_repo_name: empet/more-graphs

max_forks_repo_name: empet/more-graphs


text:
module MeanFilters

export meanfilter

const Window{N} = NTuple{N,UnitRange{Int}}

function alloc_output(x::AbstractArray)
    T = float(eltype(x))
    return similar(x, T)
end

function meanfilter1d(v::AbstractVector, window::UnitRange{Int})
    out = alloc_output(v)
    meanfilter1d!(out, v, window)
end

function meanfilter1d!(out, v::AbstractVector, window::UnitRange{Int})
    # out and v may not alias!
    if isempty(v)
        return out
    end
    firstindex_v = firstindex(v)
    lastindex_v = lastindex(v)
    eachindex_v = eachindex(v)

    max_window_len = length(window)
    i_pix   = firstindex_v
    i_left  = i_pix + first(window)
    i_right = i_pix + last(window)
    inds_first_window = clamp(i_left, firstindex_v, lastindex_v):clamp(i_right, firstindex_v, lastindex_v)
    window_len = length(inds_first_window)
    T = eltype(out)
    sum = T(Base.sum(view(v, inds_first_window)))
    out[i_pix] = sum/window_len
    # could factor this out into a fringe_loop function
    # and implement meanfilter by
    # fringe_loop(left)
    # branchless_main_loop(middle)
    # fringe_loop(right)
    while true
        i_pix += 1
        i_right += 1
        i_left += 1
        if i_pix > lastindex_v
            return out
        end
        if i_right in eachindex_v
            window_len += 1
            sum += v[i_right]
        end
        if i_left-1 in eachindex_v
            window_len -= 1
            sum -= v[i_left-1]
        end
        out[i_pix] = sum / window_len
    end
    return out
end

function resolve_range(r::AbstractRange)
    resolve_range(UnitRange{Int}(r))
end
function resolve_range(r::UnitRange{Int})
    if first(r) <= 0 <= last(r)
        return r
    else
        msg = """
        Range must contain 0. Got $r
        """
        throw(ArgumentError(msg))
    end
end
function resolve_range(x::Number)
    i = Int(x)
    if i < 1
        msg = """
        Window size must be at least 1. Got $x instead.
        """
        throw(ArgumentError(msg))
    elseif isodd(i)
        h = i÷2
        return (-h:h)
    else
        msg = """
        Window size must be odd. Got $x instead.
        """
        throw(ArgumentError(msg))
    end
end

function resolve_window(n::Integer)
    r = resolve_range(n)
    return (r,)
end
function resolve_window(r::AbstractRange)
    return resolve_window((r,))
end
function resolve_window(window)
    map(resolve_range, window)
end

function _meanfilter(x::AbstractArray{T,N}, window::Window{N}) where {T,N}
    out = alloc_output(x)
    return _meanfilter!(out, x, window)
end

"""
    out, buf2 = _meanfilter_along_axis!!!(buf, inp, axis)

`buf` and `inp` may not alias. `inp` contains the input and `buf` may contain
arbitrary junk data. After calling this function `buf, buf2, inp` contain arbitrary junk.
It `out === buf && buf2 === inp` or `out === inp && buf2 === buf` will hold.
"""
function _meanfilter_along_axis!!!(out, inp, window::Window, axis)
    N = ndims(out)
    @assert axis in 1:N
    # return out, buf such that
    r = window[axis]
    if window === 0:0
        return (inp, out)
    end
    if N === 1
        meanfilter1d!(out, inp, r)
    elseif N === 2
        if axis === 1
            Threads.@threads for iy in axes(out, 2)
                meanfilter1d!(view(out,:,iy), view(inp,:,iy), r)
            end
        else
            @assert axis === 2
            Threads.@threads for ix in axes(out, 1)
                meanfilter1d!(view(out,ix,:), view(inp,ix,:), r)
            end
        end
    else
        grid = Base.setindex(axes(out), Base.OneTo(1), axis)
        Threads.@threads for I in CartesianIndices(grid)
            inds = Base.setindex(Tuple(I), (:), axis)
            meanfilter1d!(view(out, inds...), view(inp,inds...), r)
        end
    end
    return (out, inp)
end

function _meanfilter!(out::AbstractVector, arr::AbstractVector, window::Window{1})
    inp = similar(out)
    copy!(inp, arr)
    out, _ = _meanfilter_along_axis!!!(out, inp, window, 1)
    return out
end

function _meanfilter!(out::AbstractMatrix, arr::AbstractMatrix, window::Window{2})
    inp = similar(out)
    copy!(inp, arr)
    inp, out = _meanfilter_along_axis!!!(out, inp, window, 1)
    inp, out = _meanfilter_along_axis!!!(out, inp, window, 2)
    return inp
end
function _meanfilter!(out::AbstractArray{T2,N}, arr::AbstractArray{T1,N}, window::Window{N}) where {T1,T2,N}
    inp = similar(out)
    copy!(inp, arr)
    for axis in 1:ndims(out)
        inp, out = _meanfilter_along_axis!!!(out, inp, window, axis)
    end
    return inp
end

function meanfilter!(out, arr, window)
    w::Window = resolve_window(window)
    if ndims(arr) !== length(w)
        msg = """
        Number of array and window dimensions must coincide. Got:
        ndims(arr) = $(ndims(arr))
        length(window) = $(length(w))
        """
        throw(ArgumentError(msg))
    end
    out2 = _meanfilter!(out, arr, w)
    if out2 !== out
        copy!(out, out2)
    end
    return out
end
function meanfilter(arr, window)
    out = alloc_output(arr)
    return meanfilter!(out, arr, window)
end

end


text:
# Test code for identifiying inside/outside points

using Base.Test

push!(LOAD_PATH, string(pwd(),"/src"))
using PyPlot
using ModifiedStokesSolver
import AnalyticDomains

@testset "InsideOutside" begin
    # Discretize
    numpanels = 10
    panelorder = 16
    curve = AnalyticDomains.starfish(n_arms = 7, amplitude=0.3)
    dcurve = discretize(curve, numpanels, panelorder)

    # Easy test
    rgrid = linspace(0.01, 0.7-1e-8, 50)
    tgrid = linspace(0, 2*pi, 500)
    R, T = ndgrid(rgrid, tgrid)
    X = R.*cos.(T)
    Y = R.*sin.(T)
    zt = [X[:] Y[:]]'
    interior = interior_points(dcurve, zt)
    exterior = .!interior
    @test all(interior)

    # clf()
    # plot(zt[1,interior], zt[2,interior],".b")
    # plot(zt[1,exterior], zt[2,exterior],".r")
    # plot(dcurve.points[1,:], dcurve.points[2,:], "-k")
    # axis("image")

    # Hard test
    rgrid = linspace(0, 2*pi, 1000)
    igrid = linspace(-5e-2, 5e-2, 100)
    R, I = ndgrid(rgrid, igrid)
    Z = @. curve.tau(R + 1im*I)
    X = real(Z)
    Y = imag(Z)

    zt = [X[:] Y[:]]'
    interior = interior_points(dcurve, zt)
    exterior = .!interior

    ref_interior = I[:] .> 0
    ref_exterior = I[:] .< 0

    # clf()
    # plot(zt[1,interior], zt[2,interior],".b")
    # plot(zt[1,exterior], zt[2,exterior],".r")
    # plot(dcurve.points[1,:], dcurve.points[2,:], "-k")
    # axis("image")

    @test !any(ref_interior .& exterior)
    @test !any(ref_exterior .& interior)
end


text:
@testset "Pseudodivision" begin
    R, (x, y, z) = Nemo.PolynomialRing(Nemo.QQ, ["x", "y", "z"])

    @test iszero(pseudodivision(x^2 - y^2, x - y, x))
    @test pseudodivision(y * x^4 + x^3 + x, z * x^2, x) == x
end


text:
@testset "test_simplification" begin
## SMINREAL ##
G = ss([-5 0 0 0; 0 -1 -2.5 0; 0 4 0 0; 0 0 0 -6], [2 0; 0 1; 0 0; 0 2],
       [0 3 0 0; -2 0 0 1], [0 0; 1 0])
@test sminreal(G) == G
@test sminreal(G[1, 1]) == ss(0)
@test sminreal(G[1, 2]) == ss([-1 -2.5; 4 0], [1; 0], [3 0], [0])
@test sminreal(G[2, 1]) == ss([-5], [2], [-2], [1])
@test sminreal(G[2, 2]) == ss([-6], [2], [1], [0])

using ControlSystems.DemoSystems: resonant
R = resonant()*resonant()
@test sminreal(R) == R # https://github.com/JuliaControl/ControlSystems.jl/issues/409

# https://github.com/JuliaControl/ControlSystems.jl/issues/475
A = [
 0.0  0.0   0.0;
 1.0  0.0  -0.00578297;
 0.0  1.0   0.0;
]
B = [1.0; 0.0; 0.0]
C = [0.0 0.0 1.0]
D = [0]

ss_sys = ss(A,B,C,D)
@test sminreal(ss_sys) == ss_sys

## MINREAL ##

s = tf("s")
P = [1/(s+1) 2/(s+3); 1/(s+1) 1/(s+1)]
sys = ss(P)
sysmin = minreal(sys)

@test size(sysmin.A,1) == 3 # Test that the reduction of sys worked

@test hinfnorm(sys - sysmin)[1] < 1e-15 # And that the answer is correct

@test_broken balreal(sys-sysmin)

@test all(sigma(sys-sysmin, [0.0, 1.0, 2.0])[1] .< 1e-15)  # Previously crashed because of zero dimensions in tzeros

t = 0:0.1:10
y1,x1 = step(sys,t)[[1,3]]
y2,x2 = step(sysmin,t)[[1,3]]
@test sum(abs2,y1.-y2) < 1e-6 # Test that the output from the two systems are the same

end


text:
module Accessors
using MacroTools
using MacroTools: isstructdef, splitstructdef, postwalk
using Requires: @require

include("setindex.jl")
include("optics.jl")
include("sugar.jl")
include("functionlenses.jl")
include("testing.jl")

if !isdefined(Base, :only)
    # Julia pre-1.4
    function only(x)
        length(x) == 1 || throw(ArgumentError("Collection contains $(length(x)) elements, must contain exactly 1 element"))
        first(x)
    end
end

function __init__()
    @require StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" include("staticarrays.jl")
end

end


text:
# ---
# title: 950. Reveal Cards In Increasing Order
# id: problem950
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Array
# link: <https://leetcode.com/problems/reveal-cards-in-increasing-order/description/>
# hidden: true
# ---
#
# In a deck of cards, every card has a unique integer.  You can order the deck
# in any order you want.
#
# Initially, all the cards start face down (unrevealed) in one deck.
#
# Now, you do the following steps repeatedly, until all cards are revealed:
#
#   1. Take the top card of the deck, reveal it, and take it out of the deck.
#   2. If there are still cards in the deck, put the next top card of the deck at the bottom of the deck.
#   3. If there are still unrevealed cards, go back to step 1.  Otherwise, stop.
#
# Return an ordering of the deck that would reveal the cards in **increasing
# order.**
#
# The first entry in the answer is considered to be the top of the deck.
#
#
#
# **Example 1:**
#
#
#
#     Input: [17,13,11,2,3,5,7]
#     Output: [2,13,3,11,5,17,7]
#     Explanation:
#     We get the deck in the order [17,13,11,2,3,5,7] (this order doesn't matter), and reorder it.
#     After reordering, the deck starts as [2,13,3,11,5,17,7], where 2 is the top of the deck.
#     We reveal 2, and move 13 to the bottom.  The deck is now [3,11,5,17,7,13].
#     We reveal 3, and move 11 to the bottom.  The deck is now [5,17,7,13,11].
#     We reveal 5, and move 17 to the bottom.  The deck is now [7,13,11,17].
#     We reveal 7, and move 13 to the bottom.  The deck is now [11,17,13].
#     We reveal 11, and move 17 to the bottom.  The deck is now [13,17].
#     We reveal 13, and move 17 to the bottom.  The deck is now [17].
#     We reveal 17.
#     Since all the cards revealed are in increasing order, the answer is correct.
#
#
#
#
# **Note:**
#
#   1. `1 <= A.length <= 1000`
#   2. `1 <= A[i] <= 10^6`
#   3. `A[i] != A[j]` for all `i != j`
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
using Printf
using Distances

function check_for_random_data()
    nbspectra = 20000
    nbval = 1000
    @debug "check_for_random_data: dispatching nbspecta : " nbspectra
    spectra = Array{Array{Float64,1},1}(undef, nbspectra)
    for i in 1:nbspectra
        spectra[i] = rand(nbval)
    end
    D=Jaccard()
    depth = 8
    argument = RPTreeArg(D, 8, 1.5)
    rptree=RPTree(argument, spectra)
    @time leafCenters = randomProjection(rptree)
    leaves, diameters = analyzeSplittingInfo(rptree)
    #
    if length(leaves) != 2^depth
        @printf stdout "\n bad number of leaves , shoud be %d, got %d " 2^8  length(leaves)
        return false
    end
    #
    return true
end


text:
import YAML
@info "Starting diagnostic script with" ARGS
config_file = ARGS[1]
cfg = YAML.load_file(config_file)
out_file = cfg["setting_name"]
@info "Copying file to" out_file
Base.Filesystem.cp(config_file, out_file)
@info "Done"

max_stars_repo_name: cffbots/ESMValTool

max_issues_repo_name: cffbots/ESMValTool

max_forks_repo_name: cffbots/ESMValTool


text:
# LLVM Headers

addHeaderDir(joinpath(basepath,"usr/include"))
addHeaderDir(joinpath(basepath,"deps/llvm-svn/tools/clang/"))
addHeaderDir(joinpath(basepath,"deps/llvm-svn/tools/clang/include/"))

# Load LLVM and clang libraries

defineMacro("__STDC_LIMIT_MACROS")
defineMacro("__STDC_CONSTANT_MACROS")
cxxinclude("llvm/Support/MemoryBuffer.h")
cxxinclude("llvm/Bitcode/ReaderWriter.h")
cxxinclude("llvm/IR/Module.h")
cxxinclude("llvm/IR/IRBuilder.h")
cxxinclude("llvm/IR/Constants.h")
cxxinclude("llvm/IR/CFG.h")
cxxinclude("llvm/Support/GraphWriter.h")
cxxinclude("llvm/ExecutionEngine/ExecutionEngine.h")
cxxinclude("lib/CodeGen/CGValue.h")
cxxinclude("lib/CodeGen/CodeGenTypes.h")
cxxinclude("clang/AST/DeclBase.h")
cxxinclude("clang/AST/Type.h")
cxxinclude("clang/Basic/SourceLocation.h")
cxxinclude("clang/Frontend/CompilerInstance.h")
cxxinclude("clang/AST/DeclTemplate.h")
cxxinclude("llvm/ADT/ArrayRef.h")
cxxinclude("llvm/Analysis/CallGraph.h")


text:
"""
    SymMatrixSpace()

Use in the [`@variable`](@ref) macro to constrain a matrix of variables to be
symmetric.

## Examples

```jldoctest; setup=:(model = Model())
julia> @variable(model, Q[1:2, 1:2] in SymMatrixSpace())
2×2 LinearAlgebra.Symmetric{VariableRef,Array{VariableRef,2}}:
 Q[1,1]  Q[1,2]
 Q[1,2]  Q[2,2]
```
"""
struct SymMatrixSpace end

"""
    SkewSymmetricMatrixSpace()

Use in the [`@variable`](@ref) macro to constrain a matrix of variables to be
skew-symmetric.

## Examples

```jldoctest; setup=:(model = Model())
@variable(model, Q[1:2, 1:2] in SkewSymmetricMatrixSpace())
```
"""
struct SkewSymmetricMatrixSpace end

"""
    PSDCone

Positive semidefinite cone object that can be used to constrain a square matrix
to be positive semidefinite in the [`@constraint`](@ref) macro. If the matrix
has type `Symmetric` then the columns vectorization (the vector obtained by
concatenating the columns) of its upper triangular part is constrained to belong
to the `MOI.PositiveSemidefiniteConeTriangle` set, otherwise its column
vectorization is constrained to belong to the
`MOI.PositiveSemidefiniteConeSquare` set.

## Examples

Consider the following example:
```jldoctest PSDCone; setup = :(using JuMP)
julia> model = Model();

julia> @variable(model, x)
x

julia> a = [ x 2x
            2x  x];

julia> b = [1 2
            2 4];

julia> cref = @SDconstraint(model, a ⪰ b)
[x - 1    2 x - 2;
 2 x - 2  x - 4  ] ∈ PSDCone()

julia> jump_function(constraint_object(cref))
4-element Array{GenericAffExpr{Float64,VariableRef},1}:
 x - 1
 2 x - 2
 2 x - 2
 x - 4

julia> moi_set(constraint_object(cref))
MathOptInterface.PositiveSemidefiniteConeSquare(2)
```
We see in the output of the last command that the matrix the vectorization of the
matrix is constrained to belong to the `PositiveSemidefiniteConeSquare`.

```jldoctest PSDCone
julia> using LinearAlgebra # For Symmetric

julia> cref = @constraint(model, Symmetric(a - b) in PSDCone())
[x - 1    2 x - 2;
 2 x - 2  x - 4  ] ∈ PSDCone()

julia> jump_function(constraint_object(cref))
3-element Array{GenericAffExpr{Float64,VariableRef},1}:
 x - 1
 2 x - 2
 x - 4

julia> moi_set(constraint_object(cref))
MathOptInterface.PositiveSemidefiniteConeTriangle(2)
```
As we see in the output of the last command, the vectorization of only the upper
triangular part of the matrix is constrained to belong to the
`PositiveSemidefiniteConeSquare`.
"""
struct PSDCone end

"""
    SymmetricMatrixShape

Shape object for a symmetric square matrix of `side_dimension` rows and columns.
The vectorized form contains the entries of the upper-right triangular part of
the matrix given column by column (or equivalently, the entries of the
lower-left triangular part given row by row).
"""
struct SymmetricMatrixShape <: AbstractShape
    side_dimension::Int
end
function reshape_vector(
    vectorized_form::Vector{T},
    shape::SymmetricMatrixShape,
) where {T}
    matrix = Matrix{T}(undef, shape.side_dimension, shape.side_dimension)
    k = 0
    for j in 1:shape.side_dimension
        for i in 1:j
            k += 1
            matrix[j, i] = matrix[i, j] = vectorized_form[k]
        end
    end
    return Symmetric(matrix)
end
function reshape_set(
    ::MOI.PositiveSemidefiniteConeTriangle,
    ::SymmetricMatrixShape,
)
    return PSDCone()
end
function vectorize(matrix::Matrix, ::SymmetricMatrixShape)
    n = LinearAlgebra.checksquare(matrix)
    return [matrix[i, j] for j in 1:n for i in 1:j]
end

"""
    SkewSymmetricMatrixShape

Shape object for a skew symmetric square matrix of `side_dimension` rows and
columns. The vectorized form contains the entries of the upper-right triangular
part of the matrix (without the diagonal) given column by column (or
equivalently, the entries of the lower-left triangular part given row by row).
The diagonal is zero.
"""
struct SkewSymmetricMatrixShape <: AbstractShape
    side_dimension::Int
end
function reshape_vector(
    vectorized_form::Vector{T},
    shape::SkewSymmetricMatrixShape,
) where {T}
    NewType = Base.promote_type(
        T,
        _MA.promote_operation(-, T),
        _MA.promote_operation(zero, T),
    )
    matrix = Matrix{NewType}(undef, shape.side_dimension, shape.side_dimension)
    k = 0
    for j in 1:shape.side_dimension
        for i in 1:(j-1)
            k += 1
            matrix[i, j] = vectorized_form[k]
            matrix[j, i] = -vectorized_form[k]
        end
        matrix[j, j] = zero(NewType)
    end
    return matrix
end

function vectorize(matrix::Matrix, ::SkewSymmetricMatrixShape)
    n = LinearAlgebra.checksquare(matrix)
    return [matrix[i, j] for j in 1:n for i in 1:j-1]
end

"""
    SquareMatrixShape

Shape object for a square matrix of `side_dimension` rows and columns. The
vectorized form contains the entries of the the matrix given column by column
(or equivalently, the entries of the lower-left triangular part given row by
row).
"""
struct SquareMatrixShape <: AbstractShape
    side_dimension::Int
end
function reshape_vector(
    vectorized_form::Vector{T},
    shape::SquareMatrixShape,
) where {T}
    return reshape(vectorized_form, shape.side_dimension, shape.side_dimension)
end
function reshape_set(::MOI.PositiveSemidefiniteConeSquare, ::SquareMatrixShape)
    return PSDCone()
end
vectorize(matrix::Matrix, ::SquareMatrixShape) = vec(matrix)

function vectorize(matrix, shape::Union{SymmetricMatrixShape,SquareMatrixShape})
    return vectorize(Matrix(matrix), shape)
end

function _square_side(_error::Function, ::Containers.SparseAxisArray)
    return _error("Cannot have index dependencies in symmetric variables.")
end
function _square_side(_error::Function, ::Containers.DenseAxisArray)
    return _error(
        "Index sets for symmetric variables must be ranges of the form 1:N.",
    )
end
function _square_side(_error::Function, ::Array)
    return _error("Symmetric variables must be 2-dimensional.")
end
function _square_side(_error::Function, variables::Matrix)
    n, m = size(variables)
    if n != m
        _error("Symmetric variables must be square. Got size ($n, $m).")
    end
    return n
end

function _vectorize_variables(_error::Function, matrix::Matrix)
    n = LinearAlgebra.checksquare(matrix)
    for j in 1:n
        for i in 1:j
            if matrix[i, j] != matrix[j, i]
                _error(
                    "Non-symmetric bounds, integrality or starting values for symmetric variable.",
                )
            end
        end
    end
    return vectorize(matrix, SymmetricMatrixShape(n))
end

"""
    build_variable(_error::Function, variables, ::SymMatrixSpace)

Return a `VariablesConstrainedOnCreation` of shape [`SymmetricMatrixShape`](@ref)
creating variables in `MOI.Reals`, i.e. "free" variables unless they are
constrained after their creation.

This function is used by the [`@variable`](@ref) macro as follows:
```julia
@variable(model, Q[1:2, 1:2], Symmetric)
```
"""
function build_variable(
    _error::Function,
    variables::Matrix{<:ScalarVariable},
    ::SymMatrixSpace,
)
    n = _square_side(_error, variables)
    set = MOI.Reals(MOI.dimension(MOI.PositiveSemidefiniteConeTriangle(n)))
    shape = SymmetricMatrixShape(n)
    return VariablesConstrainedOnCreation(
        _vectorize_variables(_error, variables),
        set,
        shape,
    )
end

"""
    build_variable(_error::Function, variables, ::SkewSymmetricMatrixSpace)

Return a `VariablesConstrainedOnCreation` of shape [`SkewSymmetricMatrixShape`](@ref)
creating variables in `MOI.Reals`, i.e. "free" variables unless they are
constrained after their creation.

This function is used by the [`@variable`](@ref) macro as follows:
```julia
@variable(model, Q[1:2, 1:2] in SkewSymmetricMatrixSpace())
```
"""
function build_variable(
    _error::Function,
    variables::Matrix{<:ScalarVariable},
    ::SkewSymmetricMatrixSpace,
)
    n = _square_side(_error, variables)
    set = MOI.Reals(div(n^2 - n, 2))
    shape = SkewSymmetricMatrixShape(n)
    return VariablesConstrainedOnCreation(
        vectorize(variables, SkewSymmetricMatrixShape(n)),
        set,
        shape,
    )
end

"""
    build_variable(_error::Function, variables, ::PSDCone)

Return a `VariablesConstrainedOnCreation` of shape [`SymmetricMatrixShape`](@ref)
constraining the variables to be positive semidefinite.

This function is used by the [`@variable`](@ref) macro as follows:
```julia
@variable(model, Q[1:2, 1:2], PSD)
```
"""
function build_variable(
    _error::Function,
    variables::Matrix{<:ScalarVariable},
    ::PSDCone,
)
    n = _square_side(_error, variables)
    set = MOI.PositiveSemidefiniteConeTriangle(n)
    shape = SymmetricMatrixShape(n)
    return VariablesConstrainedOnCreation(
        _vectorize_variables(_error, variables),
        set,
        shape,
    )
end

"""
    build_constraint(_error::Function, Q::Symmetric{V, M},
                     ::PSDCone) where {V <: AbstractJuMPScalar,
                                       M <: AbstractMatrix{V}}

Return a `VectorConstraint` of shape [`SymmetricMatrixShape`](@ref) constraining
the matrix `Q` to be positive semidefinite.

This function is used by the [`@constraint`](@ref) macros as follows:
```julia
@constraint(model, Symmetric(Q) in PSDCone())
```
The form above is usually used when the entries of `Q` are affine or quadratic
expressions, but it can also be used when the entries are variables to get the
reference of the semidefinite constraint, e.g.,
```julia
@variable model Q[1:2,1:2] Symmetric
# The type of `Q` is `Symmetric{VariableRef, Matrix{VariableRef}}`
var_psd = @constraint model Q in PSDCone()
# The `var_psd` variable contains a reference to the constraint
```
"""
function build_constraint(
    _error::Function,
    Q::Symmetric{V,M},
    ::PSDCone,
) where {V<:AbstractJuMPScalar,M<:AbstractMatrix{V}}
    n = LinearAlgebra.checksquare(Q)
    shape = SymmetricMatrixShape(n)
    return VectorConstraint(
        vectorize(Q, shape),
        MOI.PositiveSemidefiniteConeTriangle(n),
        shape,
    )
end

"""
    build_constraint(_error::Function,
                     Q::AbstractMatrix{<:AbstractJuMPScalar},
                     ::PSDCone)

Return a `VectorConstraint` of shape [`SquareMatrixShape`](@ref) constraining
the matrix `Q` to be symmetric and positive semidefinite.

This function is used by the [`@constraint`](@ref) and [`@SDconstraint`](@ref)
macros as follows:
```julia
@constraint(model, Q in PSDCone())
@SDconstraint(model, P ⪰ Q)
```
The [`@constraint`](@ref) call above is usually used when the entries of `Q` are
affine or quadratic expressions, but it can also be used when the entries are
variables to get the reference of the semidefinite constraint, e.g.,
```julia
@variable model Q[1:2,1:2]
# The type of `Q` is `Matrix{VariableRef}`
var_psd = @constraint model Q in PSDCone()
# The `var_psd` variable contains a reference to the constraint
```
"""
function build_constraint(
    _error::Function,
    Q::AbstractMatrix{<:AbstractJuMPScalar},
    ::PSDCone,
)
    n = LinearAlgebra.checksquare(Q)
    shape = SquareMatrixShape(n)
    return VectorConstraint(
        vectorize(Q, shape),
        MOI.PositiveSemidefiniteConeSquare(n),
        shape,
    )
end

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
# ellipse to multiple forward simulated shapes
# only final landmark positions observed

using BridgeLandmarks
using Random
using Distributions
using StaticArrays
using LinearAlgebra
using JLD2

workdir = @__DIR__
println(workdir)
cd(workdir)

n = 18  # in the final experiment we'll take n = 18
nshapes = 1
T = 1.0

a = 2.0     # Hamiltonian kernel parameter (the larger, the stronger landmarks behave similarly)
c = 0.1     # multiplicative factor in kernel
γ = 0.7     # Noise level

xobs0 = [PointF(2.0cos(t), sin(t))  for t in collect(0:(2pi/n):2pi)[2:end]]

θ, ψ =  π/4, 0.2
rot =  SMatrix{2,2}(cos(θ), sin(θ), -sin(θ), cos(θ))
stretch = SMatrix{2,2}(1.0 + ψ, 0.0, 0.0, 1.0 - ψ)
shift = [0.5, 0.0]
xobsT = [rot * stretch * xobs0[i] + shift  for i in 1:n ]

JLD2.@save "data_exp1.jld2" xobs0 xobsT n nshapes


text:
using CrystalBalls
using StaticArrays
using BenchmarkTools

lat = LatticeToroidal(2, 1.0, 3)

x = SVector(0, 0)

# U = GaugeFieldSU2(lat)
U = rand(GaugeFieldSU2, lat)


text:
"""
	msfm()

	Calculate the shortest distance from a list of
	points to all other pixels in an image, using the
	Multistencil Fast Marching Method (MSFM).

	This method gives more accurate
	distances by using second order derivatives and cross neighbours.

	T=msfm2d(F, SourcePoints, UseSecond, UseCross)

	# inputs,
	  F: The speed image. The speed function must always be larger
				than zero (min value 1e-8), otherwise some regions will
				never be reached because the time will go to infinity.
	  SourcePoints : A list of starting points [2 x N] (distance zero)
	  UseSecond : Boolean Set to true if not only first but also second
	               order derivatives are used (default)
	  UseCross : Boolean Set to true if also cross neighbours
	               are used (default)
	# outputs
	  T : Image with distance from SourcePoints to all pixels
"""
function msfm(speedimage::AbstractArray{T,2}, SourcePoints::AbstractArray{T}, usesecond::Bool=true, usecross::Bool=true,Ed=false::Bool) where T

	distanceimage = fill(T(-1.0),size(speedimage))
	# Augmented Fast Marching (For skeletonize)
	#Ed=false # Was: nargout>1;

	# Euclidian distance image
	if(Ed)
		Y = zeros(T,size(speedimage))
	end

	# Pixels which are processed and have a final distance are frozen
	Frozen = falses(size(speedimage));

	# Free memory to store neighbours of the (segmented) region
	neg_free = 100000
	neg_pos=0
	if Ed
		neg_list = zeros(4,neg_free)
	else
		neg_list = zeros(3,neg_free)
	end

	# (There are 3 pixel classes:
	#   - frozen (processed)
	#   - narrow band (boundary) (in list to check for the next pixel with smallest distance)
	#   - far (not yet used)

 	# Neighbours
	ne = [-1  0;
		   1  0;
		   0 -1;
		   0  1]

	# SourcePoints=round.(Int,floor.(SourcePoints))

	# set all starting points to distance zero and frozen
	for z=1:size(SourcePoints,2)
		# starting point
		x= SourcePoints[1,z]
		y= SourcePoints[2,z]
		# Set starting point to frozen and distance to zero
		Frozen[x,y]=true
		distanceimage[x,y]=zero(T)
	end

	# Add all neighbours of the starting points to narrow list
	for z=1:size(SourcePoints,2)
		# starting point
		x=SourcePoints[1,z]
		y=SourcePoints[2,z]
		for k=1:4
			# Location of neighbour
			i=x+ne[k,1]
			j=y+ne[k,2]
			# Check if current neighbour is not yet frozen and inside the
			# picture
			if((i>0)&&(j>0)&&(i<=size(speedimage,1))&&(j<=size(speedimage,2))&&(~Frozen[i,j]))
				Tt=1/max(speedimage[i,j],eps(T))
				Ty=1;
				# Update distance in neigbour list or add to neigbour list
				if distanceimage[i,j]>zero(T)
					if neg_list[1,distanceimage[i,j]]>Tt
						neg_list[1,distanceimage[i,j]]=Tt
					end
					if Ed
						neg_list[4,distanceimage[i,j]]=min(Ty,neg_list[4,distanceimage[i,j]]);
					end
				else
					neg_pos=neg_pos+1
					# If running out of memory at a new block
					if neg_pos>neg_free
						neg_free = neg_free + 100000
						if Ed
							neg_list = hcat(neg_list,zeros(4,neg_free))# Was: neg_list(1,neg_free)=0;
						else
							neg_list = hcat(neg_list,zeros(3,neg_free))# Was: neg_list(1,neg_free)=0;
						end
					end
					if Ed
						@views neg_list[:,neg_pos]=[Tt;i;j;Ty]
					else
						@views neg_list[:,neg_pos]=[Tt;i;j]
					end
					distanceimage[i,j]=neg_pos
				end
			end
		end
	end

	# Reuse information in loop
	si_size = size(speedimage)
	si_size_1 = size(speedimage,1)
	si_size_2 = size(speedimage,2)
	Tpatch = @MMatrix fill(T(Inf),5,5)
	Order = @MArray zeros(Int8, 1, 4);
	Tm = @MArray zeros(T,1,4)
	Tm2 = @MArray zeros(T,1,4)
	Coeff = @MArray zeros(T,3)
	Tt = zero(T) #@MArray zeros(T,2)
	Tt2 = zero(T)
	# Loop through all pixels of the image
	for itt=1:length(speedimage)
		# Get the pixel from narrow list (boundary list) with smallest
		# distance value and set it to current pixel location

		if neg_pos==0
			break
		end
		@views (t,index) = findmin(neg_list[1,1:neg_pos])
		x = round(Int, neg_list[2,index])
		y = round(Int, neg_list[3,index])
		Frozen[x,y] = true
		distanceimage[x, y] = neg_list[1,index]

		if Ed
			Y[x,y] = neg_list[4, index]
		end

		# Remove min value by replacing it with the last value in the array
		if index<neg_pos
			@views neg_list[:,index] = neg_list[:, neg_pos]
			x2 = round(Int, neg_list[2,index])
			y2 = round(Int, neg_list[3,index])
			distanceimage[x2, y2] = index
		end
		neg_pos = neg_pos - 1

		# Loop through all 4 neighbours of current pixel
		for k = 1:4
			# Location of neighbour
			i = x + ne[k,1]
			j = y + ne[k,2]

			# Check if current neighbour is not yet frozen and inside the
			# picture
			if (i>0) && (j>0) && (i<= si_size_1) && (j<= si_size_2) && (~Frozen[i,j])

				Tt = calculatedistance!(distanceimage,Tpatch,Order,Coeff,Tm,Tm2,Tt,Tt2,speedimage[i,j],si_size,i,j,usesecond,usecross,Frozen)
				if Ed
					Ty=calculatedistance!(Y,Tpatch,Order,Coeff,Tm,Tm2,Tt,Tt2,1,si_size,i,j,usesecond,usecross,Frozen)
				end

				# Update distance in neigbour list or add to neigbour list
				if distanceimage[i,j]>0
					neg_list[1,round(Int,distanceimage[i,j])] = minimum((Tt,neg_list[1,round(Int,distanceimage[i,j])]))
					if Ed
						neg_list[4,round(Int,distanceimage[i,j])] = minimum((Ty,neg_list[4,round(Int,distanceimage[i,j])]))
					end
				else
					neg_pos=neg_pos+1
					# If running out of memory at a new block
					if neg_pos>neg_free
						neg_free = neg_free + 100_000
						#Was: neg_list(1,neg_free)=0;
						if Ed
							neg_list = hcat(neg_list,zeros(4,neg_free))
						else
							neg_list = hcat(neg_list,zeros(3,neg_free))
						end
					end
					if Ed
						@views neg_list[:,neg_pos]=[Tt;i;j;Ty]
					else
						@views neg_list[:,neg_pos]=[Tt;i;j]
					end
					distanceimage[i,j]=neg_pos
				end
			end
		end
	end
	return distanceimage
end # End function msfm

function calculatedistance(TI::AbstractArray{T}, Tpatch, Fij, sizeF, i, j, usesecond, usecross, Frozen) where T
	Order = @MArray zeros(Int8, 1, 4)
	Tm = @MArray zeros(T,1,4)
	Tm2 = @MArray zeros(T,1,4)
	Coeff = @MArray zeros(T,3)
	Tt = zero(T)#@MArray zeros(T,2)
	Tt2 = zero(T)
	calculatedistance!(TI, Tpatch, Order, Coeff, Tm, Tm2, Tt, Tt2, Fij, sizeF, i, j, usesecond, usecross, Frozen)
end

function calculatedistance!(TI::AbstractArray{T}, Tpatch, Order, Coeff, Tm, Tm2, Tt, Tt2, Fij, sizeF, i, j, usesecond, usecross, Frozen) where T

	fill!(Tpatch,T(Inf))
	for nx = -2:2
		for ny = -2:2
			i_n = i + nx
			j_n = j + ny
			if i_n > 0 && j_n > 0 && i_n <= sizeF[1] && j_n <= sizeF[2] && Frozen[i_n, j_n]
				Tpatch[nx+3, ny+3] = TI[i_n, j_n]
			end
		end
	end

	# The values in order is 0 if no neighbours in that direction
	# 1 if 1e order derivatives is used and 2 if second order
	# derivatives are used
	fill!(Tm,T(0))

	# Make 1e order derivatives in x and y direction
	Tm[1] = min(Tpatch[2,3], Tpatch[4,3] )
	if isfinite(Tm[1])
		Order[1] = 1
	end

	Tm[2] = min(Tpatch[3,2], Tpatch[3,4])
	if isfinite(Tm[2])
		Order[2] = 1
	end

	# Make 1e order derivatives in cross directions
	if usecross
		Tm[3] = minimum((Tpatch[2,2], Tpatch[4,4]))
		if isfinite(Tm[3])
			Order[3] = 1
		end
		Tm[4] = min(Tpatch[2,4], Tpatch[4,2])
		if isfinite(Tm[4])
			Order[4] = 1
		end
	end

	# Make 2e order derivatives
	if usesecond
		fill!(Tm2,T(0))
		# pixels with a pixeldistance 2 from the center must be
		# lower in value otherwise use other side or first order
		ch1 = (Tpatch[1,3] < Tpatch[2,3]) && isfinite(Tpatch[2,3])
		ch2 = (Tpatch[5,3] < Tpatch[4,3]) && isfinite(Tpatch[4,3])

		if ch1 && ch2
			Tm2[1] = minimum(( (4*Tpatch[2,3] - Tpatch[1,3])/3 , (4*Tpatch[4,3] - Tpatch[5,3])/3 ))
			Order[1] = 2
		elseif ch1
			Tm2[1] = (4*Tpatch[2,3] - Tpatch[1,3])/3
			Order[1] = 2
		elseif ch2
			Tm2[1] =(4*Tpatch[4,3] - Tpatch[5,3])/3
			Order[1] = 2
		end

		ch1 = (Tpatch[3,1] < Tpatch[3,2]) && isfinite(Tpatch[3,2])
		ch2 = (Tpatch[3,5] < Tpatch[3,4]) && isfinite(Tpatch[3,4])

		if ch1 && ch2
			Tm2[2] = min((4*Tpatch[3,2] - Tpatch[3,1])/3 , (4*Tpatch[3,4] - Tpatch[3,5])/3)
			Order[2] = 2
		elseif ch1
			Tm2[2] = (4*Tpatch[3,2] - Tpatch[3,1])/3
			Order[2] = 2
		elseif ch2
			Tm2[2] = (4*Tpatch[3,4] - Tpatch[3,5])/3
			Order[2] = 2
		end

		if usecross
			ch1 = (Tpatch[1,1] < Tpatch[2,2]) && isfinite(Tpatch[2,2])
			ch2 = (Tpatch[5,5] < Tpatch[4,4]) && isfinite(Tpatch[4,4])
			if ch1 && ch2
				Tm2[3] = min( (4*Tpatch[2,2] - Tpatch[1,1])/3 , (4*Tpatch[4,4] - Tpatch[5,5])/3)
				Order[3] = 2
			elseif ch1
				Tm2[3] = (4*Tpatch[2,2] - Tpatch[1,1])/3
				Order[3] = 2
			elseif ch2
				Tm2[3] = (4*Tpatch[4,4] - Tpatch[5,5])/3
				Order[3] = 2
			end

			ch1 = (Tpatch[1,5] < Tpatch[2,4]) && isfinite(Tpatch[2,4])
			ch2 = (Tpatch[5,1] < Tpatch[4,2]) && isfinite(Tpatch[4,2])
			if ch1 && ch2
				Tm2[4] = minimum(( (4*Tpatch[2,4] - Tpatch[1,5])/3 , (4*Tpatch[4,2] - Tpatch[5,1])/3))
				Order[4] = 2
			elseif ch1
				Tm2[4] = (4*Tpatch[2,4] - Tpatch[1,5])/3
				Order[4] = 2
			elseif ch2
				Tm2[4] = (4*Tpatch[4,2] - Tpatch[5,1])/3
				Order[4] = 2
			end
		end
	else
		fill!(Tm2,T(0))
	end

	# Calculate the distance using x and y direction
	# Coeff = [0 0 -1/(max(Fij^2,eps()))]

	Coeff .= ((0, 0, -1/(max(Fij^2,eps(T)))))
	for t = 1:2
		if Order[t] == 1
			Coeff .+= (1, -2*Tm[t], Tm[t]^2)
		elseif Order[t] == 2
			Coeff .+=  (1, -2*Tm2[t], Tm2[t]^2) .* 2.25
		end
	end

	Tt = roots(Coeff)
	Tt = maximum(Tt)
	# Calculate the distance using the cross directions
	if usecross
		Coeff  .+= ((0, 0, -1/(max(Fij^2,eps(T)))))
		for t = 3:4
			if Order[t] == 1
				Coeff .+= 0.5 .* (1, -2*Tm[t], Tm[t]^2)
			elseif Order[t] == 2
				Coeff .+= 0.5 .* (1, -2*Tm2[t], Tm2[t]^2) .* 2.25
			end
		end
		T2 = roots(Coeff)
		Tt2=maximum(Tt2)

		# Select minimum distance value of both stensils
		if ~isempty(Tt2)
			Tt = minimum((Tt,Tt2))
		end
	end

	# Upwind condition check, current distance must be larger
	# then direct neighbours used in solution
	#DirectNeigbInSol=Tm(isfinite(Tm));
	#if(nnz(DirectNeigbInSol>=Tt)>0) # Will this ever happen?
	#    Tt=min(DirectNeigbInSol)+(1/(max(Fij,eps)));
	#end

end

function roots!(z, Coeff::AbstractArray{T}) where T
	a = Coeff[1]
	b = Coeff[2]
	c = Coeff[3]
	d = max((b*b)-4.0*a*c,zero(T))
	if a != 0
		z[1] = (-b - sqrt(d)) / (T(2.0)*a)
		z[2] = (-b + sqrt(d)) / (T(2.0)*a)
	else
		z[1] = (T(2.0)*c)/(-b - sqrt(d))
		z[2] = (T(2.0)*c)/(-b + sqrt(d))
	end
end

function roots_exp!(z, Coeff::AbstractArray{T}) where T
	a = Coeff[1]
	b = Coeff[2]
	c = Coeff[3]
	d = max((b*b)-4.0*a*c,zero(T))
	if a != 0
		z = max( (-b - sqrt(d)) / (T(2.0)*a), (-b + sqrt(d)) / (T(2.0)*a) )
	else
		z = max( (T(2.0)*c)/(-b - sqrt(d)), (T(2.0)*c)/(-b + sqrt(d)) )
	end
end

function roots(Coeff::AbstractArray{T}) where T
	z = @MArray zeros(T,2,1)
	roots!(z,Coeff)
	return z
end


text:
"""
    covdiff(xm, xd; include_diagonal=false)

Error of the covariance matrix.
"""
function covdiff(xm::AbstractArray, xd::AbstractArray; include_diagonal::Bool=false)
    @assert ndims(xm) == ndims(xd)
    @assert size(xm)[1:(end - 1)] == size(xd)[1:(end - 1)]
    Cm = cov(reshape(xm, :, size(xm)[end]); dims=2)
    Cd = cov(reshape(xd, :, size(xd)[end]); dims=2)
    if include_diagonal
        return norm(Cm - Cd) / length(Cm)
    else
        return norm(Cm - Diagonal(diag(Cm)) - Cd + Diagonal(diag(Cd))) / (length(Cm) - size(Cm, 1))
    end
end

"""
    adversarial_accuracy(xm, xd)

Adversarial accuracy indicator.
As defined in https://arxiv.org/abs/2105.13889.
"""
function adversarial_accuracy(xm::AbstractArray, xd::AbstractArray)

end


text:
@testset "ParallelJuliaSolver" begin
	include("testParallelJuliaSolver.jl")
end


text:
export computeTensorizedSP

function computeTensorizedSP(m::Int64,
                             α::Vector{Vector{Float64}},β::Vector{Vector{Float64}},
                             nodes::Vector{Vector{Float64}},weights::Vector{Vector{Float64}},
                             ind::Matrix{Int64};
                             issymmetric::BitArray=falses(length(α)))
    @assert m>=1 "`dimension` has to be positive"
    l,p = size(ind)
    l-=1
    T = spzeros(reduce(+,[l^i for i in 0:m]))
    # @show length(collect(with_replacement_combinations(0:l,m)))
    for tensor_ind in collect(with_replacement_combinations(0:l,m))
        index = 1 + reduce(+,[idx*l^(m-i) for (i,idx) in enumerate(tensor_ind)])
        T[index] = computeSP(tensor_ind,α,β,nodes,weights,ind;issymmetric=issymmetric)
    end
    return T
end
function computeTensorizedSP(m::Int64,op::Vector{OrthoPoly},q::Vector{Quad},ind::Matrix{Int64})
    @assert length(op)==length(q) "Inconsistent lengh of polynomials and quadrature rules"
    @assert issymmetric.(op)==issymmetric.(q)
    α,β = coeffs(op)
    n,w = nw(q)
    computeTensorizedSP(m,α,β,n,w,ind;issymmetric=issymmetric.(op))
end
function computeTensorizedSP(m::Int64,mOP::MultiOrthoPoly)
    @assert typeof(mOP.uni)==Vector{OrthoPolyQ} "no quadrature rules provided"
    ops = [ mOP.uni[i].op for i=1:length(mOP.uni) ]
    quads = [ mOP.uni[i].quad for i=1:length(mOP.uni) ]
    computeTensorizedSP(m,ops,quads,mOP.ind)
end

function computeTensorizedSP(m::Int64,opq::OrthoPolyQ)
    computeTensorizedSP(m,[opq.op],[opq.quad],calculateMultiIndices(1,deg(opq)) )
end


text:
function f_lm(x)
  [x[1], 2.0 - x[2]]
end
function g_lm(x)
  [1.0 0.0; 0.0 -1.0]
end

initial_x = [100.0, 100.0]

results = Optim.levenberg_marquardt(f_lm, g_lm, initial_x)
@assert norm(results.minimum - [0.0, 2.0]) < 0.01


text:
##  Copyright 2021 Myrl Marmarelis
##
##  Licensed under the Apache License, Version 2.0 (the "License");
##  you may not use this file except in compliance with the License.
##  You may obtain a copy of the License at
##
##    http://www.apache.org/licenses/LICENSE-2.0
##
##  Unless required by applicable law or agreed to in writing, software
##  distributed under the License is distributed on an "AS IS" BASIS,
##  WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
##  See the License for the specific language governing permissions and
##  limitations under the License.

#using Automa -- for a fast state-machine compiler

# remove duplicate entries, favoring earlier ones
function sanitize(choice::RankedChoice{R},
    n_candidates::Int)::RankedChoice{R} where R
  counts = @MVector zeros(Int, n_candidates)
  updated = MVector(choice)
  for r in 1:R
    c = choice[r]
    if c == 0
      continue
    end
    if counts[c] >= 1
      updated[r] = 0
    end
    counts[c] += 1
  end
  convert(SVector, updated)
end

function Base.vcat(issues::IssueVote{R}...)::IssueVote{R} where R
  choices = reduce(vcat, (issue.choices for issue in issues))
  n_candidates = issues[1].n_candidates
  @assert all( n_candidates == issue.n_candidates for issue in issues )
  IssueVote(choices, n_candidates)
end

function Base.filter(f::Function, issue::IssueVote{R})::IssueVote{R} where R
  IssueVote(filter(f, issue.choices), issue.n_candidates)
end

# generated functions to automatically wrap and unwrap?
function Base.map(f::Function, issue::IssueVote{R})::IssueVote where R
  IssueVote(map(f, issue.choices), issue.n_candidates)
end

function read_prm_file(filename, N, R)
  rankings = RankedChoice{R}[]
  for line in eachline(filename)
    header_content = split(line, "1) ")
    if length(header_content) != 2
      println("Line malformed: $line.")
      continue
    end
    choices = split(header_content[2], ',')
    equal_choices = filter(choices) do choice
      split(choice, '=') |> length > 1
    end
    # no support for tie votes...
    if length(equal_choices) > 0
      println("Can't set candidates equal...")
      continue
    end
    filtered_choices = filter(choices) do choice
      length(choice) > 0
    end
    candidates = map(filtered_choices) do choice
      candidate_string =
        split(split(choice, 'C')[2], '[')[1]
      candidate = parse(Int, candidate_string)
      if candidate < 1 || candidate > N
        println("Candidate number out of range: $candidate.")
      end
      candidate
    end
    n_choices = length(candidates)
    padded_candidates = vcat(candidates, zeros(Int, R-n_choices))
    push!(rankings, RankedChoice{R}(padded_candidates))
  end
  IssueVote(rankings, N)
end

using CSV, DataFrames

function read_csv(filename::String, candidates::AbstractVector{String})
  matrix = CSV.read(filename, DataFrame, header=false) |> Matrix
  parse_matrix(convert.(String, matrix), candidates)
end

using XLSX

function read_xlsx(filename::String, candidates::AbstractVector{String};
    sheet_index=1, columns::AbstractVector{Int}, first_row=2)
  sheet = XLSX.readxlsx(filename)
  matrix = sheet[sheet_index][:][first_row:end, columns]
  converted_matrix = map(matrix) do v
    v isa Missing ? "N/A" : convert(String, v)
  end # convert.(String, matrix) is such convenient notation
  parse_matrix(converted_matrix, candidates)
end

function parse_matrix(matrix::AbstractMatrix{S},
    candidates::AbstractVector{S}) where S <: AbstractString
  N = length(candidates)
  V, R = size(matrix)
  waste = Set{String}()
  sliced_votes = map(1:R) do rank
    slice = @view matrix[:, rank]
    slice_votes, slice_waste = quantize(slice, candidates)
    union!(waste, slice_waste)
    slice_votes
  end
  votes = zip(sliced_votes...) .|> RankedChoice{R}
  filtered_votes = filter(votes) do vote
    vote[1] != 0 # first entry can't be zero!! todo: ensure all other entries are zero too? if not, we could shift forward...
  end
  IssueVote(filtered_votes, N), waste
end

function quantize(items::AbstractVector{S}, keys::AbstractVector{S}) where S
  waste = Set{S}()
  quanta = map(items) do item
    for (index, key) in enumerate(keys)
      if item == key
        return index
      end
    end
    push!(waste, item)
    return 0
  end
  quanta, waste
end

# empty strings count as nothing
function parse_matrix(matrix::AbstractMatrix{S}) where S <: AbstractString
  candidates = filter(!isempty, unique(matrix)) # `!` operator creates a closure on the function automagically!
  sort!(candidates)
  votes, waste = parse_matrix(matrix, candidates)
  # (candidates, parse_matrix(matrix, candidates)...)
  votes, candidates
end


text:
# @inline tuplefirst(x) = x
# @inline tuplefirst(x::Tuple) = first(x)
# @inline tuplesecond(x) = x
# @inline tuplesecond(x::Tuple) = x[2]
# @inline tuplethird(x) = x
# @inline tuplethird(x::Tuple) = x[3]
# @inline tuplefourth(x) = x
# @inline tuplefourth(x::Tuple) = x[4]
# @inline tupletail(x) = x
# @inline tupletail(x::Tuple) = Base.tail(x)

# struct TupleLength{N} end
# @inline tuple_len(::Any) = nothing
# @inline tuple_len(::Tuple{Vararg{Any,N}}) where {N} = TupleLength{N}()
# @inline tuple_len(::TupleLength{N}, args...) where {N} = TupleLength{N}()
# @inline tuple_len(::Nothing, a, args...) = tuple_len(tuple_len(a), args...)
# @inline tuple_len(a, args...) = tuple_len(tuple_len(a), args...)

# unary # 2^2 - 2 = 2 definitions
@inline fmap(f::F, x::Tuple{X}) where {F,X} = (f(first(x)),)
@inline fmap(f::F, x::NTuple) where {F} = (f(first(x)), fmap(f, Base.tail(x))...)

# binary # 2^3 - 2 = 6 definitions
@inline fmap(f::F, x::Tuple{X}, y::Tuple{Y}) where {F,X,Y} = (f(first(x), first(y)),)
@inline fmap(f::F, x::Tuple{X}, y) where {F,X} = (f(first(x), y),)
@inline fmap(f::F, x, y::Tuple{Y}) where {F,Y} = (f(x, first(y)),)
@inline fmap(f::F, x::Tuple{Vararg{Any,N}}, y::Tuple{Vararg{Any,N}}) where {F,N} = (f(first(x), first(y)), fmap(f, Base.tail(x), Base.tail(y))...)
@inline fmap(f::F, x::Tuple, y) where {F} = (f(first(x), y), fmap(f, Base.tail(x), y)...)
@inline fmap(f::F, x, y::Tuple) where {F} = (f(x, first(y)), fmap(f, x, Base.tail(y))...)

fmap(f::F, x::Tuple{X}, y::Tuple) where {F,X} = throw("Dimension mismatch.")
fmap(f::F, x::Tuple, y::Tuple{Y}) where {F,Y} = throw("Dimension mismatch.")
fmap(f::F, x::Tuple, y::Tuple) where {F} = throw("Dimension mismatch.")

# ternary # 2^4 - 2 = 14 definitions, or 3
# @inline fmap(f::F, x, y, z) where {F} = fmap(f, tuple_len(x, y, z), x, y, z)
# @inline fmap(f::F, ::TupleLength{1}, x, y, z) where {F} = (f(tuplefirst(x), tuplefirst(y), tuplefirst(z)),)
# @inline fmap(f::F, ::TupleLength{2}, x, y, z) where {F} = (f(tuplefirst(x), tuplefirst(y), tuplefirst(z)), f(tuplesecond(x), tuplesecond(y), tuplesecond(z)))
# @inline function fmap(f::F, ::TupleLength{3}, x, y, z) where {F}
#     (
#         f(tuplefirst(x), tuplefirst(y), tuplefirst(z)),
#         f(tuplesecond(x), tuplesecond(y), tuplesecond(z)),
#         f(tuplethird(x), tuplethird(y), tuplethird(z))
#     )
# end
# @inline function fmap(f::F, ::TupleLength{4}, x, y, z) where {F}
#     (
#         f(tuplefirst(x), tuplefirst(y), tuplefirst(z)),
#         f(tuplesecond(x), tuplesecond(y), tuplesecond(z)),
#         f(tuplethird(x), tuplethird(y), tuplethird(z)),
#         f(tuplefourth(x), tuplefourth(y), tuplefourth(z))
#     )
# end
@generated function fmap(f::F, x::Vararg{Any,N}) where {F,N}
    q = Expr(:block, Expr(:meta, :inline))
    t = Expr(:tuple)
    U = 1
    call = Expr(:call, :f)
    syms = Vector{Symbol}(undef, N)
    istup = Vector{Bool}(undef, N)
    for n ∈ 1:N
        syms[n] = xₙ = Symbol(:x_, n)
        push!(q.args, Expr(:(=), xₙ, Expr(:ref, :x, n)))
        istup[n] = ist = (x[n] <: Tuple)
        if ist
            U = length(x[n].parameters)
            push!(call.args, Expr(:ref, xₙ, 1))
        else
            push!(call.args, xₙ)
        end
    end
    push!(t.args, call)
    for u ∈ 2:U
        call = Expr(:call, :f)
        for n ∈ 1:N
            xₙ = syms[n]
            if istup[n]
                push!(call.args, Expr(:ref, xₙ, u))
            else
                push!(call.args, xₙ)
            end
        end
        push!(t.args, call)
    end
    push!(q.args, t); q
end
# @inline function fmap(f::F, ::TupleLength, x, y, z) where {F}
#     (f(tuplefirst(x), tuplefirst(y), tuplefirst(z)), fmap(f, tupletail(x), tupletail(y), tupletail(z))...)
# end
# @inline fmap(f::F, ::TupleLength, x, y, z) where {F} = (f(tuplefirst(x), tuplefirst(y), tuplefirst(z)), fmap(f, tupletail(x), tupletail(y), tupletail(z))...)

# quaternary # 2^5 - 2 = 30 definitions, or 3
# @inline fmap(f::F, w, x, y, z) where {F} = fmap(f, tuple_len(w, x, y, z), w, x, y, z)
# @inline fmap(f::F, ::TupleLength{1}, w, x, y, z) where {F} = (f(tuplefirst(w), tuplefirst(x), tuplefirst(y), tuplefirst(z)),)
# @inline fmap(f::F, ::TupleLength, w, x, y, z) where {F} = (f(tuplefirst(w), tuplefirst(x), tuplefirst(y), tuplefirst(z)), fmap(f, tupletail(w), tupletail(x), tupletail(y), tupletail(z))...)

# @inline fmap(f::F, w) where {F} = VecUnroll(fmapt(f, unrolleddata(w)))
# @inline fmap(f::F, w, x) where {F} = VecUnroll(fmapt(f, unrolleddata(w), unrolleddata(x)))
# @inline fmap(f::F, w, x, y) where {F} = VecUnroll(fmapt(f, unrolleddata(w), unrolleddata(x), unrolleddata(y)))
# @inline fmap(f::F, w, x, y, z) where {F} = VecUnroll(fmapt(f, unrolleddata(w), unrolleddata(x), unrolleddata(y), unrolleddata(z)))

for op ∈ [:(-), :abs, :floor, :ceil, :trunc, :round, :sqrt, :!, :(~), :leading_zeros, :trailing_zeros]
    @eval @inline Base.$op(v1::VecUnroll{N,W,T}) where {N,W,T} = VecUnroll(fmap($op, v1.data))
end
@inline Base.inv(v::VecUnroll{N,W,T}) where {N,W,T<:Union{Float32,Float64}} = VecUnroll(fmap(vdiv, vbroadcast(Val{W}(), one(T)), v.data))
@inline Base.reinterpret(::Type{T}, v::VecUnroll) where {T<:Number} = VecUnroll(fmap(reinterpret, T, v.data))
for op ∈ [:(Base.:(+)),:(Base.:(-)),:(Base.:(*)),:(Base.:(/)),:(Base.:(%)),:(Base.:(<<)),:(Base.:(>>)),:(Base.:(>>>)),:(Base.:(&)),:(Base.:(|)),:(Base.:(⊻)),
          :(Base.:(÷)),:(Base.max),:(Base.min),:(Base.copysign),:(Base.:(<)),:(Base.:(≤)),:(Base.:(>)),:(Base.:(≥)),:(Base.:(==)),:(Base.:(≠)),:vadd,:vsub,:vmul]
    @eval begin
        # @inline $op(v1::VecUnroll{N,W,T}, v2::Real) where {N,W,T} = VecUnroll(fmap($op, v1.data, vbroadcast(Val{W}(), T, v2)))
        # @inline $op(v1::Real, v2::VecUnroll{N,W,T}) where {N,W,T} = VecUnroll(fmap($op, vbroadcast(Val{W}(), T, v1), v2.data))
        @inline $op(v1::VecUnroll{N,W,T}, v2::S) where {N,W,T,S<:Real} = VecUnroll(fmap($op, v1.data, vbroadcast(Val{W}(), v2)))
        @inline $op(v1::S, v2::VecUnroll{N,W,T}) where {N,W,T,S<:Real} = VecUnroll(fmap($op, vbroadcast(Val{W}(), v1), v2.data))
        # @inline $op(v1::VecUnroll{N,W,T}, v2::AbstractSIMD{W,S}) where {N,W,T,S<:Real} = VecUnroll(fmap($op, v1.data, promote_type(Vec{W,T}, Vec{W,S})(v2)))
        # @inline $op(v1::AbstractSIMD{W,S}, v2::VecUnroll{N,W,T}) where {N,W,T,S<:Real} = VecUnroll(fmap($op, promote_type(Vec{W,T}, Vec{W,S})(v1), v2.data))
        # @inline $op(v1::Real, v2::VecUnroll{N,W,T}) where {N,W,T} = VecUnroll(fmap($op, Vec{W,T}(v1), v2.data))
        @inline $op(v1::VecUnroll, v2::VecUnroll) = VecUnroll(fmap($op, v1.data, v2.data))
        @inline $op(v1::VecUnroll{N,W,T}, ::StaticInt{M}) where {N,W,T,M} = VecUnroll(fmap($op, v1.data, vbroadcast(Val{W}(), T(M))))
        @inline $op(::StaticInt{M}, v1::VecUnroll{N,W,T}) where {N,W,T,M} = VecUnroll(fmap($op, vbroadcast(Val{W}(), T(M)), v1.data))
    end
end
for op ∈ [:%, :&, :|, :⊻, :>>, :>>>, :<<,:(<),:(≤),:(>),:(≥),:(==),:(≠)]
    @eval begin
        @inline Base.$op(vu::VecUnroll, i::MM) = $op(vu, Vec(i))
        @inline Base.$op(i::MM, vu::VecUnroll) = $op(Vec(i), vu)
    end
end
for op ∈ [:>>, :>>>, :<<]
    @eval @inline Base.$op(m::Mask, vu::VecUnroll) = $op(Vec(m), vu)
end
# @inline Base.:(^)(v1::VecUnroll{N,W,T}, v2::Rational) where {N,W,T} = VecUnroll(fmap(^, v1.data, Vec{W,T}(v2)))
@inline Base.copysign(v1::Rational, v2::VecUnroll{N,W,T}) where {N,W,T} = VecUnroll(fmap(copysign, Vec{W,T}(v1), v2.data))
@inline Base.copysign(v1::Signed, v2::VecUnroll{N,W,T}) where {N,W,T} = VecUnroll(fmap(copysign, Vec{W,T}(v1), v2.data))
for op ∈ [:rotate_left,:rotate_right,:funnel_shift_left,:funnel_shift_right]
    @eval begin
        @inline $op(v1::VecUnroll{N,W,T}, v2::Real) where {N,W,T} = VecUnroll(fmap($op, v1.data, Vec{W,T}(v2)))
        @inline $op(v1::Real, v2::VecUnroll{N,W,T}) where {N,W,T} = VecUnroll(fmap($op, Vec{W,T}(v1), v2.data))
        @inline $op(v1::VecUnroll, v2::VecUnroll) = VecUnroll(fmap($op, v1.data, v2.data))
    end
end

for op ∈ [:(Base.muladd), :(Base.fma), :vfmadd, :vfnmadd, :vfmsub, :vfnmsub, :vfmadd231, :vfnmadd231, :vfmsub231, :vfnmsub231]
    @eval begin
        @inline $op(v1::VecUnroll{N,W,T}, v2::Real, v3::Real) where {N,W,T} = VecUnroll(fmap($op, v1.data, Vec{W,T}(v2), Vec{W,T}(v3)))
        @inline $op(v1::Real, v2::VecUnroll{N,W,T}, v3::Real) where {N,W,T} = VecUnroll(fmap($op, Vec{W,T}(v1), v2.data, Vec{W,T}(v3)))
        @inline $op(v1::Real, v2::Real, v3::VecUnroll{N,W,T}) where {N,W,T} = VecUnroll(fmap($op, Vec{W,T}(v1), Vec{W,T}(v2), v3.data))
        @inline $op(v1::VecUnroll{N,W,T}, v2::VecUnroll{N,W,T}, v3::Real) where {N,W,T} = VecUnroll(fmap($op, v1.data, v2.data, Vec{W,T}(v3)))
        @inline $op(v1::VecUnroll{N,W,T}, v2::Real, v3::VecUnroll{N,W,T}) where {N,W,T} = VecUnroll(fmap($op, v1.data, Vec{W,T}(v2), v3.data))
        @inline $op(v1::Real, v2::VecUnroll{N,W,T}, v3::VecUnroll{N,W,T}) where {N,W,T} = VecUnroll(fmap($op, Vec{W,T}(v1), v2.data, v3.data))
        $op(v1::VecUnroll, v2::VecUnroll, v3::Real) = throw("Size mismatch")
        $op(v1::VecUnroll, v2::Real, v3::VecUnroll) = throw("Size mismatch")
        $op(v1::Real, v2::VecUnroll, v3::VecUnroll) = throw("Size mismatch")
        # @inline $op(v1::VecUnroll, v2::VecUnroll, v3::VecUnroll) = VecUnroll(fmap($op, v1.data, v2.data, v3.data))
        $op(v1::VecUnroll, v2::VecUnroll, v3::VecUnroll) = throw("Size mismatch")
        @inline $op(v1::VecUnroll{N,W,T}, v2::VecUnroll{N,W,T}, v3::VecUnroll{N,W,T}) where {N,W,T} = VecUnroll(fmap($op, v1.data, v2.data, v3.data))
    end
end
@inline ifelse(v1::VecUnroll{N,W,<:Boolean}, v2::T, v3::T) where {N,W,T<:NativeTypes} = VecUnroll(fmap(ifelse, v1.data, Vec{W,T}(v2), Vec{W,T}(v3)))
@inline ifelse(v1::VecUnroll{N,W,<:Boolean}, v2::T, v3::T) where {N,W,T<:Real} = VecUnroll(fmap(ifelse, v1.data, v2, v3))
@inline ifelse(v1::Vec{W,Bool}, v2::VecUnroll{N,W,T}, v3::Union{NativeTypes,AbstractSIMDVector,StaticInt}) where {N,W,T} = VecUnroll(fmap(ifelse, Vec{W,T}(v1), v2.data, Vec{W,T}(v3)))
@inline ifelse(v1::Vec{W,Bool}, v2::Union{NativeTypes,AbstractSIMDVector,StaticInt}, v3::VecUnroll{N,W,T}) where {N,W,T} = VecUnroll(fmap(ifelse, Vec{W,T}(v1), Vec{W,T}(v2), v3.data))
@inline ifelse(v1::VecUnroll{N,W,<:Boolean}, v2::VecUnroll{N,W,T}, v3::Union{NativeTypes,AbstractSIMDVector,StaticInt}) where {N,W,T} = VecUnroll(fmap(ifelse, v1.data, v2.data, Vec{W,T}(v3)))
@inline ifelse(v1::VecUnroll{N,W,<:Boolean}, v2::Union{NativeTypes,AbstractSIMDVector,StaticInt}, v3::VecUnroll{N,W,T}) where {N,W,T} = VecUnroll(fmap(ifelse, v1.data, Vec{W,T}(v2), v3.data))
@inline ifelse(v1::Vec{W,Bool}, v2::VecUnroll{N,W,T}, v3::VecUnroll{N,W,T}) where {N,W,T} = VecUnroll(fmap(ifelse, Vec{W,T}(v1), v2.data, v3.data))
# ifelse(v1::VecUnroll, v2::VecUnroll, v3::Real) = throw("Size mismatch")
# ifelse(v1::VecUnroll, v2::Real, v3::VecUnroll) = throw("Size mismatch")
# ifelse(v1::Real, v2::VecUnroll, v3::VecUnroll) = throw("Size mismatch")
# ifelse(v1::VecUnroll, v2::VecUnroll, v3::VecUnroll) = throw("Size mismatch")
@inline ifelse(v1::VecUnroll{N,W,<:Boolean}, v2::VecUnroll{N,W,T}, v3::VecUnroll{N,W,T}) where {N,W,T} = VecUnroll(fmap(ifelse, v1.data, v2.data, v3.data))
@inline function ifelse(v1::VecUnroll{N,W,<:Boolean}, v2::VecUnroll{N,W}, v3::VecUnroll{N,W}) where {N,W}
    v4, v5 = promote(v2, v3)
    VecUnroll(fmap(ifelse, v1.data, v4.data, v5.data))
end

@inline Base.:(^)(v::VecUnroll, i::Integer) = VecUnroll(fmap(^, v.data, i))

@inline Base.:(==)(v::VecUnroll{N,W,T}, x::AbstractIrrational) where {N,W,T} = v == vbroadcast(Val{W}(), T(x))
@inline Base.:(==)(x::AbstractIrrational, v::VecUnroll{N,W,T}) where {N,W,T} = vbroadcast(Val{W}(), T(x)) == v

@inline Base.convert(::Type{T}, v::VecUnroll) where {T<:Real} = VecUnroll(fmap(convert, T, v.data))
@inline Base.convert(::Type{U}, v::VecUnroll) where {N,W,T,U<:VecUnroll{N,W,T}} = VecUnroll(fmap(convert, T, v.data))
@inline Base.convert(::Type{VU}, v::VU) where {N,W,T,VU <: VecUnroll{N,W,T}} = v
@inline Base.unsafe_trunc(::Type{T}, v::VecUnroll) where {T<:Real} = VecUnroll(fmap(unsafe_trunc, T, v.data))
@inline Base.:(%)(v::VecUnroll, ::Type{T}) where {T<:Real} = VecUnroll(fmap(%, v.data, T))

@inline (::Type{VecUnroll{N,W,T,V}})(vu::VecUnroll{N,W,T,V}) where {N,W,T,V<:AbstractSIMDVector{W,T}} = vu
@inline (::Type{VecUnroll{N,W,T,VT}})(vu::VecUnroll{N,W,S,VS})  where {N,W,T,VT<:AbstractSIMDVector{W,T},S,VS<:AbstractSIMDVector{W,S}} = VecUnroll(fmap(convert, Vec{W,T}, vu.data))

function collapse_expr(N, op)
    N += 1
    t = Expr(:tuple); s = Vector{Symbol}(undef, N)
    for n ∈ 1:N
        s_n = s[n] = Symbol(:v_, n)
        push!(t.args, s_n)
    end
    q = quote
        $(Expr(:meta,:inline))
        $t = data(vu)
    end
    while N > 1
        for n ∈ 1:N >>> 1
            push!(q.args, Expr(:(=), s[n], Expr(:call, op, s[n], s[n + (N >>> 1)])))
        end
        isodd(N) && push!(q.args, Expr(:(=), s[1], Expr(:call, op, s[1], s[N])))
        N >>>= 1
    end
    q
end
@generated collapse_add(vu::VecUnroll{N}) where {N} = collapse_expr(N, :vadd)
@generated collapse_mul(vu::VecUnroll{N}) where {N} = collapse_expr(N, :vmul)
@generated collapse_max(vu::VecUnroll{N}) where {N} = collapse_expr(N, :max)
@generated collapse_min(vu::VecUnroll{N}) where {N} = collapse_expr(N, :min)
@generated collapse_and(vu::VecUnroll{N}) where {N} = collapse_expr(N, :&)
@generated collapse_or(vu::VecUnroll{N}) where {N} = collapse_expr(N, :|)
@inline vsum(vu::VecUnroll) = vsum(collapse_add(vu))
@inline vprod(vu::VecUnroll) = vprod(collapse_mul(vu))
@inline vmaximum(vu::VecUnroll) = vmaximum(collapse_max(vu))
@inline vminimum(vu::VecUnroll) = vminimum(collapse_min(vu))
@inline vall(vu::VecUnroll) = vall(collapse_and(vu))
@inline vany(vu::VecUnroll) = vany(collapse_or(vu))


text:
"""
```
function LineNormals2D(vertices)
```
This function calculates the normals, of the line points
using the neighbouring points of each contour point, and
forward an backward differences on the end points
"""
function LineNormals2D(vertices)
    n = size(vertices,1)
    # Calculate tangent vectors
    DT = float.(vertices[1:n-1,:] .- vertices[2:n,:])

    # Make influence of tangent vector 1/Distance
    # (Weighted Central Differences. Points which are closer give a
    # more accurate estimate of the normal)
    LL=sqrt.(DT[ :,1].^2 .+ DT[:,2].^2);
    DT ./= max.(LL.^2, eps())

    D = vcat([0 0], DT) + vcat(DT, [0 0])
    #D=D1+D2;
    # Normalize the normal
    LL=sqrt.(D[:,1].^2+D[:,2].^2);
    N = zeros(size(vertices))
    N[:,1] = -D[:,2] ./ LL;
    N[:,2] = D[:, 1]./ LL;

    return N
end

struct BoundaryGradient
end

function boundary_map(b::BoundaryGradient, img::Matrix)
    (nr, nc) = size(img)
    gradX = zeros(size(img))
    for c = 2:nc-1
        gradX[:,c] = 0.5*(img[:,c+1] - img[:,c-1]);
    end
    gradX[:,1] = img[:,2] - img[:,1];
    gradX[:,nc] = img[:,nc] - img[:,nc-1];
    gradY = zeros(size(img))
    for r = 2:nr-1
        gradY[r,:] = 0.5*(img[r+1,:] - img[r-1,:]);
    end
    gradY[1,:] = img[2,:] - img[1,:];
    gradY[nr,:] = img[nr,:] - img[nr-1,:];
    return (abs.(gradX) .+ abs.(gradY)) .> 0;
end

struct BoundaryShift
end

""""
```
function boundary_map(seg::Matrix{T}) where T<:Integer
```
From a segmentation, compute a binary boundary map with 1 pixel wide  boundaries.
The boundary pixels are offset by 1/2 pixel towards the origin from the actual segment boundary.
"""
function boundary_map(cfg::BoundaryShift, seg::Matrix{T}) where T<:Integer
    (h,w) = size(seg)

    local ee = zeros(size(seg));
    local s = zeros(size(seg));
    local se = zeros(size(seg));

    ee[:,1:end-1] = seg[ :,2:end];
    s[1:end-1,:] = seg[2:end,:];
    se[1:end-1,1:end-1] = seg[2:end,2:end];

    local b = (seg.!=ee) .| (seg.!=s) .| (seg.!=se);
    b[end,:] = seg[end,:] .!= ee[end,:];
    b[:,end] = seg[:,end] .!= s[:,end];
    b[end,end] = 0;
    return b
end

struct BoundaryMooreTracing
end
function boundary_map(cfg::BoundaryMooreTracing, binary::Matrix)

    # Pad the input image with a 1-px border.
    ( rows, columns ) = size( binary );
    padded = falses(rows +2, columns +2)
    padded[ 2 : rows + 1, 2 : columns + 1 ] = binary;

    # Remove interior pixels with all 4-connected neighbors.
    N = circshift( padded, [  0  1 ] );
    S = circshift( padded, [  0 -1 ] );
    E = circshift( padded, [ -1  0 ] );
    W = circshift( padded, [  1  0 ] );
    boundary_image = padded .- (( padded .+ N .+ S .+ E .+ W) .== 5 );

    # To prevent reallocating boundary, we need to initialize it.
    boundary_size = sum( boundary_image);
    boundary = sizehint!(Vector(), boundary_size)

    # Scan for the first pixel, Left-to-Right & Top-to-Bottom.
    initial_entry = (0,0)
    for  j = 1:columns, i = 1:rows
        if binary[ i, j ] == 1
            initial_entry = (j,i) .+ 1
           break;
       end
    end
    # Set this pixel ( w/ padded offset ) as the initial entry point.

    # Designate a directional offset array for search positions.
    # [ 2 ][ 3 ][ 4 ]
    # [ 1 ][ X ][ 5 ]
    # [ 8 ][ 7 ][ 6 ]
    # Column 1: x-axis offset // Column 2: y-axis offset
    neighborhood = [ -1 0; -1 -1; 0 -1; 1 -1; 1 0; 1 1; 0 1; -1 1 ];
    exit_direction = [ 7 7 1 1 3 3 5 5 ];

    # Find the first point in the boundary, Moore-Neighbor of entry point.
    n = 0
    initial_position = -1
    for n = 1:8 # 8-connected neighborhood
        c = initial_entry .+ neighborhood[n,:];
        if padded[ c[ 2 ], c[ 1 ] ]
            initial_position = c;
            break;
        end
    end

    # Set next direction based on found pixel ( i.e. 3 -> 1).
    initial_direction = exit_direction[n];

    # Start the boundary set with this pixel.
    push!(boundary,  initial_position);

    # Initialize variables for boundary search.
    position = initial_position;
    direction = initial_direction;
    # return a list of the ordered boundary pixels.
    while true

        # Find the next neighbor with a clockwise search.
        for n in circshift( collect(1:8), [1-direction ] )
            c = position .+ neighborhood[ n, : ];
            if padded[ c[ 2 ], c[1 ] ]
                position = c;
                break;
            end
        end

        # Neighbor found, save its information.
        direction = exit_direction[n];
        push!(boundary, position)

        # Entered the initial pixel the same way twice, the end.
        if  (position == initial_position)  && ( direction == initial_direction )
           break;
        end
    end
    return [[j[2]-1;j[1]-1] for j in boundary]
end


text:
using WANPerformanceAnalyzer
using Plots

@show minimum(diff([1.0, 2.0, 3.0, 4.0]))
tr = Trace([1.0, 2.0, 3.0, 4.0], [3, 0, 0, 0])
c = 1.0
τ = 1.0
v = 0.5
cgf = CGF(tr, τ, c)
θ = collect(-2.0:0.1:2.0)
plot(θ, map(t -> rate(cgf, t, v), θ))
plot!(θ, map(t -> drate(cgf, t, v), θ))

max_stars_repo_name: BenLauwens/WANPerformanceAnalyzer

max_issues_repo_name: BenLauwens/WANPerformanceAnalyzer

max_forks_repo_name: BenLauwens/WANPerformanceAnalyzer


text:
using ITensors,
      LinearAlgebra, # For tr()
      Random,        # To set a seed
      Test

Random.seed!(12345)

digits(::Type{T},i,j,k) where {T} = T(i*10^2+j*10+k)

@testset "ITensor constructors" begin
  i = Index(2,"i")
  j = Index(2,"j")
  k = Index(2,"k")

  @testset "Default" begin
    A = ITensor()
    @test store(A) isa Dense{Nothing}
    @test isNull(A)
  end

  @testset "Undef with index" begin
    A = ITensor(undef, i)
    @test store(A) isa Dense{Float64}
    @test !isNull(A)
  end

  @testset "Default with indices" begin
    A = ITensor(i,j)
    @test store(A) isa Dense{Float64}
    @test !isNull(A)
  end

  @testset "Random" begin
    A = randomITensor(i,j)
    @test store(A) isa Dense{Float64}

    @test ndims(A) == order(A) == 2 == length(inds(A))
    @test size(A) == dims(A) == (2,2)

    @test !isNull(A)

    B = randomITensor(IndexSet(i,j))
    @test store(B) isa Dense{Float64}
    @test ndims(B) == order(B) == 2 == length(inds(B))
    @test size(B) == dims(B) == (2,2)
    @test !isNull(B)
  end

  @testset "From matrix" begin
    M = [1 2; 3 4]
    A = ITensor(M,i,j)
    @test store(A) isa Dense{Float64}

    @test M ≈ Matrix(A,i,j)
    @test M' ≈ Matrix(A,j,i)
    @test_throws DimensionMismatch Vector(A)

    @test size(A,1) == size(M,1) == 2
    @test size(A,3) == size(M,3) == 1
    @test_throws ErrorException size(A,0)
    @test_throws ErrorException size(M,0)

    M = [1 2 3; 4 5 6]
    @test_throws DimensionMismatch ITensor(M,i,j)
  end

  @testset "To Matrix" begin
    TM = randomITensor(i,j)

    M1 = Matrix(TM)
    for ni in i, nj in j
      @test M1[ni,nj] ≈ TM[i(ni),j(nj)]
    end

    M2 = Matrix(TM,j,i)
    for ni in i, nj in j
      @test M2[nj,ni] ≈ TM[i(ni),j(nj)]
    end

    T3 = randomITensor(i,j,k)
    @test_throws DimensionMismatch Matrix(T3,i,j)
  end

  @testset "To Vector" begin
    TV = randomITensor(i)

    V = Vector(TV)
    for ni in i
      @test V[ni] ≈ TV[i(ni)]
    end

    T2 = randomITensor(i,j)
    @test_throws DimensionMismatch Vector(T2)
  end

  @testset "Complex" begin
    A = ITensor(Complex,i,j)
    @test store(A) isa Dense{ComplexF64}
  end

  @testset "Random complex" begin
    A = randomITensor(Complex,i,j)
    @test store(A) isa Dense{ComplexF64}
  end

  @testset "From complex matrix" begin
    M = [1+2im 2; 3 4]
    A = ITensor(M,i,j)
    @test store(A) isa Dense{ComplexF64}
  end

end

@testset "Convert to complex" begin
  i = Index(2,"i")
  j = Index(2,"j")
  A = randomITensor(i,j)
  B = complex(A)
  for ii ∈ dim(i), jj ∈ dim(j)
    @test complex(A[i(ii),j(jj)]) == B[i(ii),j(jj)]
  end
end

@testset "similar" begin
  i = Index(2,"i")
  j = Index(2,"j")
  A = randomITensor(i,j)
  B = similar(A)
  @test inds(B) == inds(A)
  @test_throws ErrorException similar(A, ComplexF32)
end

@testset "fill!" begin
  i = Index(2,"i")
  j = Index(2,"j")
  A = randomITensor(i,j)
  fill!(A, 1.0)
  @test all(data(store(A)) .== 1.0)
end

@testset "copyto!" begin
  i = Index(2,"i")
  j = Index(2,"j")
  M = [1 2; 3 4]
  A = ITensor(M,i,j)
  N = 2*M
  B = ITensor(N,i,j)
  copyto!(A, B)
  @test A == B
  @test data(store(A)) == vec(N)
  A = ITensor(M,i,j)
  B = ITensor(N,j,i)
  copyto!(A, B)
  @test A == B
  @test data(store(A)) == vec(transpose(N))
end

@testset "Unary -" begin
  i = Index(2,"i")
  j = Index(2,"j")
  M = [1 2; 3 4]
  A = ITensor(M,i,j)
  @test -A == ITensor(-M, i, j)
end

@testset "dot" begin
  i = Index(2,"i")
  a = [1.0; 2.0]
  b = [3.0; 4.0]
  A = ITensor(a,i)
  B = ITensor(b,i)
  @test dot(A, B) == 11.0
end

@testset "add and axpy" begin
  i = Index(2,"i")
  a = [1.0; 2.0]
  b = [3.0; 4.0]
  A = ITensor(a,i)
  B = ITensor(b,i)
  c = [5.0; 8.0]
  @test A + B == ITensor([4.0; 6.0], i)
  @test axpy!(2.0, A, B) == ITensor(c, i)
  a = [1.0; 2.0]
  b = [3.0; 4.0]
  A = ITensor(a,i)
  B = ITensor(b,i)
  c = [8.0; 12.0]
  @test add!(A, 2.0, 2.0, B) == ITensor(c, i)

end

@testset "mul! and rmul!" begin
  i = Index(2,"i")
  a = [1.0; 2.0]
  b = [2.0; 4.0]
  A = ITensor(a,i)
  A2, A3 = copy(A), copy(A)
  B = ITensor(b,i)
  @test mul!(A2, A, 2.0) == B == ITensors.add!(A2, 0, 2, A)
  @test rmul!(A, 2.0) == B == ITensors.scale!(A3, 2)

  i = Index(2,"i")
  j = Index(2,"j")
  M = [1 2; 3 4]
  A = ITensor(M,i,j)
  N = 2*M
  B = ITensor(N,j,i)
  @test data(store(mul!(B, A, 2.0))) == 2.0*vec(transpose(M))
end

@testset "show" begin
  i = Index(2,"i")
  a = [1.0; 2.0]
  A = ITensor(a,i)
  s = split(sprint(show, A), '\n')
  @test s[1] == "ITensor ord=1 " * sprint(show, i) * " "
  @test s[2] == "Dense{Float64}"
  @test s[3] == " 1.0"
  @test s[4] == " 2.0"
end

@testset "Test isapprox for ITensors" begin
        m,n = rand(0:20,2)
        i = Index(m)
        j = Index(n)
        realData = rand(m,n)
        complexData = realData+ zeros(m,n)*1im
        A = ITensor(realData, i,j)
        B = ITensor(complexData, i,j)
        @test A≈B
        @test B≈A
        realDataT = Array(transpose(realData))
        A = ITensor(realDataT, j,i)
        @test A≈B
        @test B≈A
    end

@testset "ITensor tagging and priming" begin
  s1 = Index(2,"Site,s=1")
  s2 = Index(2,"Site,s=2")
  l = Index(3,"Link")
  A1 = randomITensor(s1,l,l')
  A2 = randomITensor(s2,l',l'')
  @testset "findindex(::ITensor,::String)" begin
    @test s1==findindex(A1,"Site")
    @test s1==findindex(A1,"s=1")
    @test s1==findindex(A1,"s=1,Site")
    @test l==findindex(A1,"Link,0")
    @test l'==findindex(A1,"1")
    @test l'==findindex(A1,"Link,1")
    @test s2==findindex(A2,"Site")
    @test s2==findindex(A2,"s=2")
    @test s2==findindex(A2,"Site")
    @test s2==findindex(A2,"0")
    @test s2==findindex(A2,"s=2,0")
    @test s2==findindex(A2,"Site,0")
    @test s2==findindex(A2,"s=2,Site,0")
    @test l'==findindex(A2,"1")
    @test l'==findindex(A2,"Link,1")
    @test l''==findindex(A2,"2")
    @test l''==findindex(A2,"Link,2")
  end
  @testset "addtags(::ITensor,::String,::String)" begin
    s1u = addtags(s1,"u")
    lu = addtags(l,"u")

    A1u = addtags(A1,"u")
    @test hasinds(A1u,s1u,lu,lu')

    A1u = addtags(A1,"u","Link")
    @test hasinds(A1u,s1,lu,lu')

    A1u = addtags(A1,"u","0")
    @test hasinds(A1u,s1u,lu,l')

    A1u = addtags(A1,"u","Link,0")
    @test hasinds(A1u,s1,lu,l')

    A1u = addtags(A1,"u","Link,1")
    @test hasinds(A1u,s1,l,lu')
  end
  @testset "removetags(::ITensor,::String,::String)" begin
    A2r = removetags(A2,"Site")
    @test hasinds(A2r,removetags(s2,"Site"),l',l'')

    A2r = removetags(A2,"Link","1")
    @test hasinds(A2r,s2,removetags(l,"Link")',l'')
  end
  @testset "replacetags(::ITensor,::String,::String)" begin
    s2tmp = replacetags(s2,"Site","Temp")
    ltmp = replacetags(l,"Link","Temp")

    A2r = replacetags(A2,"Site","Temp")
    @test hasinds(A2r,s2tmp,l',l'')

    A2r = replacetags(A2,"Link","Temp")
    @test hasinds(A2r,s2,ltmp',ltmp'')

    A2r = replacetags(A2,"Link","Temp","1")
    @test hasinds(A2r,s2,ltmp',l'')

    A2r = replacetags(A2,"Link,2","Temp,3")
    @test hasinds(A2r,s2,l',ltmp''')

    A2r = replacetags(A2,"1","5")
    @test hasinds(A2r,s2,prime(l,5),l'')

    #In-place version
    cA2 = copy(A2)
    replacetags!(cA2,"1","5")
    @test hasinds(cA2,s2,prime(l,5),l'')
  end
  @testset "prime(::ITensor,::String)" begin
    A2p = prime(A2)
    @test A2p==A2'
    @test hasinds(A2p,s2',l'',l''')

    A2p = prime(A2,2)
    A2p = A2''
    @test hasinds(A2p,s2'',l''',l'''')

    A2p = prime(A2,"s=2")
    @test hasinds(A2p,s2',l',l'')
  end

  @testset "mapprime" begin
    @test hasinds(mapprime(A2,1,7),s2,l^7,l'')
    @test hasinds(mapprime(A2,0,1),s2',l',l'')
  end
  @testset "setprime" begin
    @test hasinds(setprime(A2,2,s2),s2'',l',l'')
    @test hasinds(setprime(A2,0,l''),s2,l',l)
  end
  @testset "swapprime" begin
    @test hasinds(swapprime(A2,1,3),l''',s2,l'')
  end
end

@testset "Converting Real and Complex Storage" begin

  @testset "Add Real and Complex" begin
    i = Index(2,"i")
    j = Index(2,"j")
    TC = randomITensor(ComplexF64,i,j)
    TR = randomITensor(Float64,i,j)

    S1 = TC+TR
    S2 = TR+TC
    @test typeof(S1.store) == Dense{ComplexF64}
    @test typeof(S2.store) == Dense{ComplexF64}
    for ii=1:dim(i),jj=1:dim(j)
      @test S1[i(ii),j(jj)] ≈ TC[i(ii),j(jj)]+TR[i(ii),j(jj)]
      @test S2[i(ii),j(jj)] ≈ TC[i(ii),j(jj)]+TR[i(ii),j(jj)]
    end
  end

end

@testset "ITensor, Dense{$SType} storage" for SType ∈ (Float64,ComplexF64)
  mi,mj,mk,ml,mα = 2,3,4,5,6,7
  i = Index(mi,"i")
  j = Index(mj,"j")
  k = Index(mk,"k")
  l = Index(ml,"l")
  α = Index(mα,"alpha")
  @testset "Set and get values with IndexVals" begin
    A = ITensor(SType,i,j,k)
    for ii ∈ 1:dim(i), jj ∈ 1:dim(j), kk ∈ 1:dim(k)
      A[k(kk),j(jj),i(ii)] = digits(SType,ii,jj,kk)
    end
    for ii ∈ 1:dim(i), jj ∈ 1:dim(j), kk ∈ 1:dim(k)
      @test A[j(jj),k(kk),i(ii)]==digits(SType,ii,jj,kk)
    end
    @test_throws ErrorException A[1]
  end
  @testset "Test permute(ITensor,Index...)" begin
    A = randomITensor(SType,i,k,j)
    permA = permute(A,k,j,i)
    @test k==inds(permA)[1]
    @test j==inds(permA)[2]
    @test i==inds(permA)[3]
    for ii ∈ 1:dim(i), jj ∈ 1:dim(j), kk ∈ 1:dim(k)
      @test A[k(kk),i(ii),j(jj)]==permA[i(ii),j(jj),k(kk)]
    end
    for ii ∈ 1:dim(i), jj ∈ 1:dim(j), kk ∈ 1:dim(k)
      @test A[k(kk),i(ii),j(jj)]==permA[i(ii),j(jj),k(kk)]
    end
    @testset "getindex and setindex with vector of IndexVals" begin
        k_inds = [k(kk) for kk ∈ 1:dim(k)]
        for ii ∈ 1:dim(i), jj ∈ 1:dim(j)
          @test A[k_inds,i(ii),j(jj)]==permA[i(ii),j(jj),k_inds]
        end
        for ii ∈ 1:dim(i), jj ∈ 1:dim(j)
            A[k_inds,i(ii),j(jj)]=collect(1:length(k_inds))
        end
        permA = permute(A,k,j,i)
        for ii ∈ 1:dim(i), jj ∈ 1:dim(j)
          @test A[k_inds,i(ii),j(jj)]==permA[i(ii),j(jj),k_inds]
        end
    end
  end
  @testset "Set and get values with Ints" begin
    A = ITensor(SType,i,j,k)
    A = permute(A,k,i,j)
    for ii ∈ 1:dim(i), jj ∈ 1:dim(j), kk ∈ 1:dim(k)
      A[kk,ii,jj] = digits(SType,ii,jj,kk)
    end
    A = permute(A,i,j,k)
    for ii ∈ 1:dim(i), jj ∈ 1:dim(j), kk ∈ 1:dim(k)
      @test A[ii,jj,kk]==digits(SType,ii,jj,kk)
    end
  end
  @testset "Test scalar(ITensor)" begin
    x = SType(34)
    A = ITensor(x)
    @test x==scalar(A)
    A = ITensor(SType,i,j,k)
    @test_throws ArgumentError scalar(A)
    # test the storage_scalar error throw
    ds = Dense{Float64}(rand(10))
    @test_throws ErrorException ITensor.storage_scalar(ds)
  end
  @testset "Test norm(ITensor)" begin
    A = randomITensor(SType,i,j,k)
    @test norm(A)≈sqrt(scalar(dag(A)*A))
  end
  @testset "Test add ITensors" begin
    A = randomITensor(SType,i,j,k)
    B = randomITensor(SType,k,i,j)
    C = A+B
    for ii ∈ 1:dim(i), jj ∈ 1:dim(j), kk ∈ 1:dim(k)
      @test C[i(ii),j(jj),k(kk)]==A[j(jj),i(ii),k(kk)]+B[i(ii),k(kk),j(jj)]
    end
    @test Array(permute(C,i,j,k))==Array(permute(A,i,j,k))+Array(permute(B,i,j,k))
  end

  @testset "Test factorizations of an ITensor" begin

    A = randomITensor(SType,i,j,k,l)

    @testset "Test SVD of an ITensor" begin
      U,S,V,u,v = svd(A,(j,l))
      @test store(S) isa Diag{Vector{Float64}}
      @test A≈U*S*V
      @test U*dag(prime(U,u))≈δ(SType,u,u') atol=1e-14
      @test V*dag(prime(V,v))≈δ(SType,v,v') atol=1e-14
    end

    @testset "Test SVD truncation" begin
        M = randn(4,4) + randn(4,4)*1.0im
        (U,s,Vh) = svd(M)
        ii = Index(4)
        jj = Index(4)
        S = Diagonal(s)
        T = ITensor(IndexSet(ii,jj),Dense{ComplexF64}(vec(U*S*Vh)))
        (U,S,Vh) = svd(T,ii;maxdim=2)
        @test norm((U*S)*Vh-T)≈sqrt(s[3]^2+s[4]^2)
    end

    @testset "Test QR decomposition of an ITensor" begin
      Q,R,q = qr(A,(i,l))
      q = commonindex(Q,R)
      @test A≈Q*R
      @test Q*dag(prime(Q,q))≈δ(SType,q,q') atol=1e-14
    end

    @testset "Test polar decomposition of an ITensor" begin
      U,P,u = polar(A,(k,l))
      @test A≈U*P
      #Note: this is only satisfied when left dimensions
      #are greater than right dimensions
      UUᵀ =  U*dag(prime(U,u))
      for ii ∈ dim(u[1]), jj ∈ dim(u[2])
        @test UUᵀ[u[1](ii),u[2](jj),u[1]'(ii),u[2]'(jj)]≈one(SType) atol=1e-14
      end
    end

    # TODO: need to implement swapInds
    #@testset "Test eigen decomposition of an ITensor" begin
    #  A = A + swapInds(dag(A),(i,k),(j,l))
    #  U,D,u = eigen(A,(i,k),(j,l))
    #  @test A≈U*D*prime(dag(U))
    #  UUᵀ =  U*prime(dag(U),u)
    #  @test UUᴴ ≈ δ(u,u') atol=1e-14
    #end
  end # End ITensor factorization testset
end # End Dense storage test


text:
using NDPrototype
using Test
using LightGraphs

using NDPrototype: VertexBatch, parameter_range

@testset "NDPrototype.jl" begin
    @testset "constructor" begin
        g = complete_graph(10)
        vertexf = ODEVertex(; f=x->x^2, dim=1, pdim=2)
        edgef = StaticEdge(; f=x->x^2, dim=2, pdim=3, coupling=AntiSymmetric())

        nd = Network(g, vertexf, edgef; verbose=false);

        g = SimpleGraph(4)
        add_edge!(g, 1, 2)
        add_edge!(g, 1, 2)
    end

    @testset "Vertex batch" begin
        vb = VertexBatch([1,2,3,4], # vertices
                         sum, # function
                         3, # dimension
                         1, # first index in state vector
                         2, # p dim
                         4 # first index of p
                         )
        @test parameter_range(vb, 1) == 4:5
        @test parameter_range(vb, 2) == 6:7
        @test parameter_range(vb, 3) == 8:9
        @test parameter_range(vb, 4) == 10:11
    end
end

@testset "apply accumulation" begin
    using NDPrototype: apply_accumulation!
    acc = zeros(4)
    src_r, dst_r = 1:2, 3:4
    edge = [-1.0, 1.0]
    apply_accumulation!(Symmetric(), +, acc, src_r, dst_r, edge)
    @test acc == [-1, 1, -1, 1]

    acc = zeros(4)
    apply_accumulation!(AntiSymmetric(), +, acc, src_r, dst_r, edge)
    @test acc == [1, -1, -1, 1]

    acc = zeros(4)
    apply_accumulation!(Directed(), +, acc, src_r, dst_r, edge)
    @test acc == [0, 0, -1, 1]

    acc = zeros(4)
    edge = [1,2,3,4,5,6]
    apply_accumulation!(Fiducial(), +, acc, src_r, dst_r, edge)
    @test acc == [4, 5, 1, 2]
end

@testset "batch_identical" begin
    using NDPrototype: batch_identical
    v = :foo
    idx = [1,7,2,5]
    @test batch_identical(v, idx) == ([:foo], [idx])
    v = [:foo, :foo, :bar, :baz, :foo]
    idx = [5,4,3,2,1]
    @test batch_identical(v, idx) == ([:foo, :bar, :baz], [[5,4,1], [3], [2]])
end

@testset "test Index manager" begin
    using NDPrototype: _lastinrange, IndexManager, _lastindex_data, _lastindex_para, isdense
    im = IndexManager()
    @test isdense(im)

    @test _lastindex_data(im) == 0
    @test _lastindex_para(im) == 0
    im.v_data[1] = 1:6
    im.e_para[1] = 1:3
    @test isdense(im)
    @test _lastindex_data(im) == 6
    @test _lastindex_para(im) == 3
    im.e_data[2] = 8:8
    im.v_para[2] = 4:5
    @test !isdense(im)
    @test _lastindex_data(im) == 8
    @test _lastindex_para(im) == 5

    using NDPrototype: register_edges!, register_vertices!
    im = IndexManager()
    @test (1, 1) == register_edges!(im, [1,2,3,4], 2, 3)
    @test (9, 13) == register_edges!(im, [5,6], 1, 1)
    @test im.e_data.keys == [1,2,3,4,5,6]
    @test im.e_para.keys == [1,2,3,4,5,6]
    @test all(values(im.e_data) .== [1:2, 3:4, 5:6, 7:8, 9:9, 10:10])
    @test all(values(im.e_para) .== [1:3, 4:6, 7:9, 10:12, 13:13, 14:14])
    @test (11, 15) == register_vertices!(im, [1,2], 2, 3)
    @test (15, 21) == register_vertices!(im, [5,6], 1, 1)
    @test im.v_data.keys == [1,2,5,6]
    @test im.v_para.keys == [1,2,5,6]
    @test all(values(im.v_data) .== [11:12,13:14,15:15,16:16])
    @test all(values(im.v_para) .== [15:17,18:20,21:21,22:22])
    @test_throws ErrorException register_vertices!(im, [1,2], 1, 1)
    @test isdense(im)
end

@testset "greedy edge coloring" begin
    using NDPrototype: color_edges_greedy, isvalid
    for i in 1:20
        g = complete_graph(i)
        colors = color_edges_greedy(g)
        @test isvalid(g, colors)
    end
end


text:
Base.@deprecate method(x) summary(x)

const has_deprecated_fminbox = Ref(false)
function optimize(
        df::OnceDifferentiable,
        initial_x::Array{T},
        l::Array{T},
        u::Array{T},
        ::Type{Fminbox};
        x_tol::T = eps(T),
        f_tol::T = sqrt(eps(T)),
        g_tol::T = sqrt(eps(T)),
        allow_f_increases::Bool = true,
        iterations::Integer = 1_000,
        store_trace::Bool = false,
        show_trace::Bool = false,
        extended_trace::Bool = false,
        callback = nothing,
        show_every::Integer = 1,
        linesearch = LineSearches.HagerZhang{T}(),
        eta::Real = convert(T,0.4),
        mu0::T = convert(T, NaN),
        mufactor::T = convert(T, 0.001),
        precondprep = (P, x, l, u, mu) -> precondprepbox!(P, x, l, u, mu),
        optimizer = ConjugateGradient,
        optimizer_o = Options(store_trace = store_trace,
                                          show_trace = show_trace,
                                          extended_trace = extended_trace),
        nargs...) where T<:AbstractFloat
        if !has_deprecated_fminbox[]
            @warn("Fminbox with the optimizer keyword is deprecated, construct Fminbox{optimizer}() and pass it to optimize(...) instead.")
            has_deprecated_fminbox[] = true
        end
        optimize(df, initial_x, l, u, Fminbox{optimizer}();
                 allow_f_increases=allow_f_increases,
                 iterations=iterations,
                 store_trace=store_trace,
                 show_trace=show_trace,
                 extended_trace=extended_trace,
                 show_every=show_every,
                 callback=callback,
                 linesearch=linesearch,
                 eta=eta,
                 mu0=mu0,
                 mufactor=mufactor,
                 precondprep=precondprep,
                 optimizer_o=optimizer_o)
end

function optimize(::AbstractObjective)
    throw(ErrorException("Optimizing an objective `obj` without providing an initial `x` has been deprecated without backwards compatability. Please explicitly provide an `x`: `optimize(obj, x)``"))
end
function optimize(::AbstractObjective, ::Method)
    throw(ErrorException("Optimizing an objective `obj` without providing an initial `x` has been deprecated without backwards compatability. Please explicitly provide an `x`: `optimize(obj, x, method)``"))
end
function optimize(::AbstractObjective, ::Method, ::Options)
    throw(ErrorException("Optimizing an objective `obj` without providing an initial `x` has been deprecated without backwards compatability. Please explicitly provide an `x`: `optimize(obj, x, method, options)``"))
end
function optimize(::AbstractObjective, ::Options)
    throw(ErrorException("Optimizing an objective `obj` without providing an initial `x` has been deprecated without backwards compatability. Please explicitly provide an `x`: `optimize(obj, x, options)``"))
end

function optimize(df::OnceDifferentiable,
    l::Array{T},
    u::Array{T},
    F::Fminbox{O}; kwargs...) where {T<:AbstractFloat,O<:AbstractOptimizer}
    throw(ErrorException("Optimizing an objective `obj` without providing an initial `x` has been deprecated without backwards compatability. Please explicitly provide an `x`: `optimize(obj, x, l, u, method, options)``"))
end


text:
# %% Plot distribution of payouts when the global population plays a round of
# the game
using DelimitedFiles
using Gadfly
import Cairo, Fontconfig
using DataFrames
using Printf
using Compose
# %%
data = vec(readdlm("../../data/world.csv"))
# %%
df_wide = DataFrame(x=1:length(data), d=data, l=log2.(data))
df = stack(df_wide, [:d, :l])

# %%
theme = Theme(major_label_font="PragmataPro Mono",
              minor_label_font="PragmataPro Mono",
              key_title_font="PragmataPro Mono Bold",
              highlight_width=0.0025w,
              bar_highlight=color("black"))
Gadfly.push_theme(theme)

ylabs = Dict("d"=>"Linear", "l"=>"log₂")
p = plot(df,
         x=:x,
         y=:value,
         ygroup=:variable,
         Geom.subplot_grid(Geom.bar,
                           free_y_axis=true),
         Guide.xlabel("PAYOUT"),
         Guide.ylabel("NUMBER OF PEOPLE"),
         Guide.title("World Population Simulation Payouts"),
         Scale.ygroup(labels=i->ylabs[i]),
         Scale.x_continuous(labels = x -> @sprintf("2%d",x)),
         )
draw(PDF("../../figures/world.pdf"), p)
# %%

max_stars_repo_name: akoen/petersburg

max_issues_repo_name: akoen/petersburg

max_forks_repo_name: akoen/petersburg