text:
type _datalist_add{T}
    this::T
end

function (_this::_datalist_add)(xs...)
    push!(_this.this.value, xs...)
    _this.this
end

type _datalist_append{T}
    this::T
end

function (_this::_datalist_append)(x)
    append!(_this.this.value, x)
    _this.this
end

type _datalist_set{T}
    this::T
end

function (_this::_datalist_set)(x)
    _this.this.value = x
    _this.this
end

type _datalist_get{T}
    this::T
end

function (_this::_datalist_get)(x)
    _this.this.value
end

type DataList
    value::Vector

    add::_datalist_add{DataList}
    push::_datalist_add{DataList}
    push!::_datalist_add{DataList}

    append::_datalist_append{DataList}
    append!::_datalist_append{DataList}

    set::_datalist_set{DataList}
    set!::_datalist_set{DataList}

    get::_datalist_get{DataList}

    DataList(x=[]) = begin
        this = new(x)

        this.add = _datalist_add{DataList}(this)
        this.push = this.add
        this.push! = this.add

        this.append = _datalist_append{DataList}(this)
        this.append! = this.append

        this.set = _datalist_set{DataList}(this)
        this.set! = this.set

        this.get = _datalist_get{DataList}(this)

        this
    end
end

type _dataitem_set{T}
    this::T
end

function (_this::_dataitem_set)(x)
    _this.this.value = x
    _this.this
end

type _dataitem_get{T}
    this::T
end

function (_this::_dataitem_get)(x)
    _this.this.value
end

type DataItem
    value::Any

    DataItem() = begin
        this = new()

        this.set = _dataitem_set{DataItem}(this)
        this.set! = this.set

        this.get = _dataitem_get{DataItem}(this)

        this
    end
end

type DataSet
    x::DataList
    y::DataList
    z::DataList
    g::DataList
    s::DataList

    DataSet() = DataSet(DataList(), DataList(), DataList(), DataList(), DataList())
end


text:
struct EmptyBridge{T} <: MOIB.AbstractBridge
end

function EmptyBridge{T}(model::MOI.ModelLike, f::MOI.AbstractVectorFunction,
                        s::EmptyCone) where {T}
    @assert MOI.output_dimension(f) == MOI.dimension(s)
    return EmptyBridge{T}()
end

function MOI.supports_constraint(::Type{<:EmptyBridge},
                                 ::Type{<:MOI.AbstractVectorFunction},
                                 ::Type{<:EmptyCone})
    return true
end
function MOIB.added_constraint_types(::Type{<:EmptyBridge{T}}) where {T}
    # TODO remove vov-in-Nonneg when MOI v0.8.3 is released
    return Tuple{DataType, DataType}[(MOI.VectorOfVariables, MOI.Nonnegatives)]
end
function MOIB.concrete_bridge_type(::Type{<:EmptyBridge{T}},
                                   ::Type{<:MOI.AbstractVectorFunction},
                                   ::Type{EmptyCone}) where T
    return EmptyBridge{T}
end

# Indices
function MOI.delete(model::MOI.ModelLike, bridge::EmptyBridge) end

# TODO ConstraintPrimal
function MOI.get(::MOI.ModelLike,
                 ::Union{MOI.ConstraintDual, MOI.ConstraintPrimal},
                 bridge::EmptyBridge{T}) where T
    return T[]
end

ext: jl


text:
# # [Tonti Diagram tooling Examples](@id tonti_examples)
#
#md # [![](https://img.shields.io/badge/show-nbviewer-579ACA.svg)](@__NBVIEWER_ROOT_URL__/examples/Tonti_examples.ipynb)

using Catlab.Graphs
using Catlab.Graphics
using Catlab.CategoricalAlgebra
using CombinatorialSpaces
using CombinatorialSpaces: ∧
using Distributions
using CairoMakie
using DifferentialEquations
using Decapods.TontiDiagrams

using Logging: global_logger
using TerminalLoggers: TerminalLogger
global_logger(TerminalLogger())

function gen_form(s::EmbeddedDeltaSet2D, f::Function)
  map(f, point(s))
end

s = EmbeddedDeltaSet2D("meshes/naca0012_8.stl");
sp = Space(s)
sd = sp.sd

# Get boundary masks for BCs
boundary_e = findall(x -> x != 0, boundary(Val{2},s) * fill(1,ntriangles(s)))
boundary_v = unique(vcat(s[boundary_e,:src],s[boundary_e,:tgt]))

obj_v = filter(x->all((-15.0,-10.0,0.0) .<= s[x,:point] .<= (15.0,10.0,0.0)), boundary_v)
obj_e = unique(vcat(incident(s, obj_v,:src)..., incident(s, obj_v,:tgt)...));

left_v = filter(x->all((-50,-15.0,0.0) .<= s[x,:point] .<= (-50.0,15.0,0.0)), boundary_v)
left_e = unique(vcat(incident(s, left_v,:src)..., incident(s, left_v,:tgt)...));

right_v = filter(x->all((50,-15.0,0.0) .<= s[x,:point] .<= (50.0,15.0,0.0)), boundary_v)
right_e = unique(vcat(incident(s, right_v,:src)..., incident(s, right_v,:tgt)...));

top_v = filter(x->all((-50,15.0,0.0) .<= s[x,:point] .<= (50.0,15.0,0.0)), boundary_v)
top_e = unique(vcat(incident(s, top_v,:src)..., incident(s, top_v,:tgt)...));

bot_v = filter(x->all((-50,-15.0,0.0) .<= s[x,:point] .<= (50.0,-15.0,0.0)), boundary_v)
bot_e = unique(vcat(incident(s, bot_v,:src)..., incident(s, bot_v,:tgt)...));

center_v = filter(x->all((-15.0,-10.0,0.0) .<= s[x,:point] .<= (15.0,10.0,0.0)), 1:nv(s))

# Define DEC operators (will be moved to TontiDiagram tooling)

td = TontiDiagrams.TontiDiagram()

# Diffusion
add_variables!(td, (:C,0,true),(:ϕ,1,false),(:∑ϕ,2,false), (:∂C,0,true),
                   (:v,1,true), (:Cv,1,true))
add_derivatives!(td, sp, :ϕ=>:∑ϕ)
add_time_dep!(td, :∂C, :C)
add_laplacian!(td, sp, :C, :∂C; coef=0.1)

add_transition!(td, [:∑ϕ], (x,y)->(x.=sp.hodge[2,3]*y) ,[:∂C])

add_transition!(td, [:C, :v], (Cv,C,v)->(Cv .= ∧(Tuple{0,1},sd,C,v)) ,[:Cv])
add_transition!(td, [:Cv], (x,y)->(x.=sp.hodge[1,2]*y) ,[:ϕ])

# Flow
add_variables!(td, (:u,1,false),(:∂u,1,false), (:∂v,1,true))
add_transition!(td, [:v], (u,v)->(u .= sp.hodge[1,2]*v) ,[:u])
add_transition!(td, [:u,:v],
  (∂u,u,v)->( (sp.hodge[1,2]*∧(Tuple{1,0},sd,v,sp.hodge[2,3]*sp.boundary[2,2]*u) .+ sp.boundary[2,1]*sp.hodge[1,3]*∧(Tuple{1,1},sd,v,sp.hodge[2,2]*u)))
  ,[:∂u])
add_transition!(td, [:∂u], (∂v, ∂u)->(∂v .= sp.hodge[2,2]*∂u) ,[:∂v])
lap_arr = sp.boundary[1,1]*sp.hodge[2,3]*sp.boundary[2,2]*sp.hodge[1,2] + sp.hodge[2,2]*sp.boundary[2,1]*sp.hodge[1,3]*sp.boundary[1,2]
add_laplacian!(td, sp, :v, :∂v; coef = 1)
add_time_dep!(td, :∂v, :v)

add_bc!(td, :∂v, v->(v[vcat(left_e,right_e)].=0))
add_bc!(td, :v, v->(v[vcat(top_e,bot_e, obj_e)].=0))

# Pressure

add_variables!(td, (:p,0,true),(:∂p,0,true))
add_time_dep!(td, :∂p, :p)
pressure_op = 30 * sp.hodge[2,3] * sp.boundary[2,2]
add_transition!(td, [:u],  (∂p,u)->(∂p .= pressure_op*u), [:∂p])
depress_op = -1 .* sp.hodge[1,2]*sp.boundary[1,1]
add_transition!(td, [:p],  (∂u,p)->(∂u .= depress_op*p), [:∂u])

data, sim = vectorfield(td, sp);
@show data

c_range = range(data[:C]...,step=1)
v_range = range(data[:v]...,step=1)
p_range = range(data[:p]...,step=1)
u = zeros(Float64,maximum(last.(values(data))))

c = gen_form(s, x->pdf(MultivariateNormal([-25,0],[4.0,2.0]),[x[1],0]))
p = gen_form(s, x->0.0)

velocity(x) = begin
  amp = 2.0
  amp * Point{3,Float64}(-1,0,0)
end
v = ♭(sd, DualVectorField(velocity.(sd[triangle_center(sd),:dual_point])))

u[c_range] .= c
u[v_range] .= v
u[p_range] .= p

tspan=(0.0,15.0)
prob = ODEProblem(sim, u, tspan)
sol = solve(prob, Tsit5(), progress=true, progress_steps=1);

fig, ax, ob = mesh(s, color=sol(1)[1:nv(s)])
save("res.svg",fig)

t = tspan[2]
times = range(0,tspan[2], length=150)
colors = [sol(t)[c_range] for t in times]
figure, axis, scatter_thing = mesh(s, color=colors[1],
                                   colorrange=(minimum(vcat(colors...)),
                                               maximum(vcat(colors...))))
axis.aspect = AxisAspect(100.0/30.0)
framerate = 30

record(figure, "flow_conc.gif", collect(1:length(collect(times))); framerate = framerate) do i
  scatter_thing.color = colors[i]
end

times = range(0,tspan[2], length=150)
colors = [sol(t)[p_range] for t in times]

color_range = [sol(t)[p_range][center_v] for t in times]
figure, axis, scatter_thing = mesh(s, color=colors[1],
                                   colorrange=(minimum(vcat(color_range...)),
                                               maximum(vcat(color_range...))))
axis.aspect = AxisAspect(100.0/30.0)
framerate = 30

record(figure, "flow_press.gif", collect(1:length(collect(times))); framerate = framerate) do i
  scatter_thing.color = colors[i]
end


text:
# Generate a toy problem without any obvious structure in the mean, precision, or noise std.
# Important to ensure that the unit tests don't just pass for a special case by accident.
# Everything should be reasonably well conditioned.
function generate_toy_problem(rng, N, D, ::Type{<:Matrix})
    X = randn(rng, D, N)
    B = randn(rng, D, D)
    C = 0.1 * randn(rng, N, N)
    mw, Λw, Σy = randn(rng, D), B * B' + I, C * C' + I
    return X, BayesianLinearRegressor(mw, Λw), Σy
end

function generate_toy_problem(rng, N, D, ::Type{<:ColVecs})
    X, f, Σy = generate_toy_problem(rng, N, D, Matrix)
    return ColVecs(X), f, Σy
end

function generate_toy_problem(rng, N, D, ::Type{<:RowVecs})
    X, f, Σy = generate_toy_problem(rng, N, D, Matrix)
    return RowVecs(collect(X')), f, Σy
end

# Some type-piracy.
Base.isapprox(dx::NamedTuple{(:X,)}, dy::ColVecs) = isapprox(dx.X, dy.X)

Base.isapprox(dx::NamedTuple{(:X,)}, dy::RowVecs) = isapprox(dx.X, dy.X)

# simple nonlinear mapping
ϕ(x::RowVecs) = RowVecs(hcat(ones(length(x)), prod.(x)))
ϕ(x::ColVecs) = ColVecs(hcat(ones(length(x)), prod.(x))')
ϕ(x::AbstractMatrix) = ϕ(ColVecs(x)).X


text: include("LinearFilter.jl")


text:
srand(1)

# TODO: Random links gets set up properly on instantiation.

context("landscape.jl") do
  context("NKLandscape") do
    l = NK.NKLandscape(4, 1)
    @fact length(l.links) --> l.n
    @fact length(l.contribs) --> l.n

    context("when near=true") do
      l = NK.NKLandscape(4, 2; near=true)
      @fact l.links[1] --> 0b1011
      @fact l.links[2] --> 0b0111
      @fact l.links[3] --> 0b1110
      @fact l.links[4] --> 0b1101
    end
  end

  context("NKqLandscape") do
    l = NK.NKqLandscape(4, 1, 2)
    @fact length(l.links) --> l.n
    @fact length(l.contribs) --> l.n

    context("when near=true") do
      l = NK.NKqLandscape(4, 2, 2; near=true)
      @fact l.links[1] --> 0b1011
      @fact l.links[2] --> 0b0111
      @fact l.links[3] --> 0b1110
      @fact l.links[4] --> 0b1101
    end
  end

  context("NKpLandscape") do
    l = NK.NKpLandscape(4, 1, 0.90)
    @fact length(l.links) --> l.n
    @fact length(l.contribs) --> l.n

    context("when near=true") do
      l = NK.NKpLandscape(4, 2, 0.90; near=true)
      @fact l.links[1] --> 0b1011
      @fact l.links[2] --> 0b0111
      @fact l.links[3] --> 0b1110
      @fact l.links[4] --> 0b1101
    end
  end
end


text:
using Knet
using CUDArt
using CUSPARSE
using Base.Test

a = rand(Float32,5,3)
b = convert(SparseMatrixCSC{Float32}, sprand(1000,3,.005))
b[10,:] = 1
c = a*b'
sc = sparse(c)
da = CudaArray(a)
db = CudaSparseMatrixCSC(b)
cc = similar(sc)
rand!(cc.nzval)
dc1 = CudaSparseMatrixCSR(cc)
Knet.A_mul_Bt_slow!(dc1, da, db)
sc1 = Knet.to_host2(dc1)
c1 = full(sc1)
@test @show c == c1
dc2 = CudaSparseMatrixCSR(cc)
Knet.A_mul_Bt!(dc2, da, db)
sc2 = Knet.to_host2(dc2)
c2 = full(sc2)
@test @show c == c2

T=1000
I=1000
J=20
K=10000

foo1(c,a,b)=(for t=1:T; Knet.A_mul_Bt_slow!(c,a,b); end)
foo2(c,a,b)=(for t=1:T; Knet.A_mul_Bt!(c,a,b); end)

a = rand(Float32, I, J)
b = convert(SparseMatrixCSC{Float32}, sprand(K,J,1/K))
c = a*b'
sc = sparse(c)
da = CudaArray(a)
db = CudaSparseMatrixCSC(b)
cc = similar(sc)
rand!(cc.nzval)
dc1 = CudaSparseMatrixCSR(cc)
dc2 = CudaSparseMatrixCSR(cc)
@date foo1(dc1, da, db)
@date foo1(dc1, da, db)
@date foo2(dc2, da, db)
@date foo2(dc2, da, db)

sc1 = Knet.to_host2(dc1)
c1 = full(sc1)
@test @show c == c1

sc2 = Knet.to_host2(dc2)
c2 = full(sc2)
@test @show c == c2


text:
export defaultmixin

makeresponse(r::Resource, req, id, res) = Response(res)

add_allow_header(r::Resource, req, id, res::Response) = if !haskey(res.headers, "Allow")
    res.headers["Allow"] = join(keys(r.methods), ", ")
end

response_option(r::Resource, req, id, res) = if res.status == 405 && req[:method] == :OPTIONS
    Response(200, help(r))
end

@mixin defaultmixin begin
    :onresponse => [makeresponse, response_option, add_allow_header]
end


text:
# Shared Memory (part of B.2)

export @cuStaticSharedMem, @cuDynamicSharedMem

shmem_id = 0

"""
    @cuStaticSharedMem(T::Type, dims) -> CuDeviceArray{T,AS.Shared}

Get an array of type `T` and dimensions `dims` (either an integer length or tuple shape)
pointing to a statically-allocated piece of shared memory. The type should be statically
inferable and the dimensions should be constant, or an error will be thrown and the
generator function will be called dynamically.
"""
macro cuStaticSharedMem(T, dims)
    # FIXME: generating a unique id in the macro is incorrect, as multiple parametrically typed
    #        functions will alias the id (and the size might be a parameter). but incrementing in
    #        the @generated function doesn't work, as it is supposed to be pure and identical
    #        invocations will erroneously share (and even cause multiple shmem globals).
    id = gensym("static_shmem")

    quote
        len = prod($(esc(dims)))
        ptr = _shmem(Val($(QuoteNode(id))), $(esc(T)), Val(len))
        CuDeviceArray($(esc(dims)), ptr)
    end
end

"""
    @cuDynamicSharedMem(T::Type, dims, offset::Integer=0) -> CuDeviceArray{T,AS.Shared}

Get an array of type `T` and dimensions `dims` (either an integer length or tuple shape)
pointing to a dynamically-allocated piece of shared memory. The type should be statically
inferable or an error will be thrown and the generator function will be called dynamically.

Note that the amount of dynamic shared memory needs to specified when launching the kernel.

Optionally, an offset parameter indicating how many bytes to add to the base shared memory
pointer can be specified. This is useful when dealing with a heterogeneous buffer of dynamic
shared memory; in the case of a homogeneous multi-part buffer it is preferred to use `view`.
"""
macro cuDynamicSharedMem(T, dims, offset=0)
    id = gensym("dynamic_shmem")

    # TODO: boundscheck against %dynamic_smem_size (currently unsupported by LLVM)

    quote
        len = prod($(esc(dims)))
        ptr = _shmem(Val($(QuoteNode(id))), $(esc(T))) + $(esc(offset))
        CuDeviceArray($(esc(dims)), ptr)
    end
end

# get a pointer to shared memory, with known (static) or zero length (dynamic shared memory)
@generated function _shmem(::Val{id}, ::Type{T}, ::Val{len}=Val(0)) where {id,T,len}
    eltyp = convert(LLVMType, T)

    # old versions of GPUArrays invoke _shmem with an integer id; make sure those are unique
    if !isa(id, String) || !isa(id, Symbol)
        id = "shmem$id"
    end

    T_ptr = convert(LLVMType, DevicePtr{T,AS.Shared})

    # create a function
    llvm_f, _ = create_function(T_ptr)

    # create the global variable
    mod = LLVM.parent(llvm_f)
    gv_typ = LLVM.ArrayType(eltyp, len)
    gv = GlobalVariable(mod, gv_typ, GPUCompiler.safe_name(string(id)), #=addrspace=# 3)
    if len > 0
        # static shared memory should be demoted to local variables, whenever possible.
        # this is done by the NVPTX ASM printer:
        # > Find out if a global variable can be demoted to local scope.
        # > Currently, this is valid for CUDA shared variables, which have local
        # > scope and global lifetime. So the conditions to check are :
        # > 1. Is the global variable in shared address space?
        # > 2. Does it have internal linkage?
        # > 3. Is the global variable referenced only in one function?
        linkage!(gv, LLVM.API.LLVMInternalLinkage)
        initializer!(gv, null(gv_typ))
    end
    # by requesting a larger-than-datatype alignment, we might be able to vectorize.
    # we pick 32 bytes here, since WMMA instructions require 32-byte alignment.
    # TODO: Make the alignment configurable
    alignment!(gv, Base.max(32, Base.datatype_alignment(T)))

    # generate IR
    Builder(JuliaContext()) do builder
        entry = BasicBlock(llvm_f, "entry", JuliaContext())
        position!(builder, entry)

        ptr = gep!(builder, gv, [ConstantInt(0, JuliaContext()),
                                 ConstantInt(0, JuliaContext())])

        val = ptrtoint!(builder, ptr, T_ptr)
        ret!(builder, val)
    end

    call_function(llvm_f, DevicePtr{T,AS.Shared})
end


text:
using ImageComponentAnalysis
using Test, Images

@testset "ImageComponentAnalysis.jl" begin
    include("test_images.jl")
    include("contour_tracing.jl")
end


text:
# Autogenerated wrapper script for RLIBM_32_jll for x86_64-w64-mingw32
export floatMathLib

JLLWrappers.@generate_wrapper_header("RLIBM_32")
JLLWrappers.@declare_library_product(floatMathLib, "floatMathLib.dll")
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_library_product(
        floatMathLib,
        "bin\\floatMathLib.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
using CSV
a = gf()
CSV.write("c:/data/a.csv", a)
;gzip c:/data/a.csv

using CodecZlib, BufferedStreams
io = open("c:/data/a.csv.gz") |> GzipDecompressorStream |> BufferedInputStream
@time a2 = CSV.read(io)
close(io)

using ZipFile, CSV, DataFrames, BufferedStreams

a = DataFrame(a = 1:3)
CSV.write("c:/data/a.csv", a)

# zip the file; Windows users who do not have zip available on the PATH can manual zip the CSV
;zip c:/data/a.zip c:/data/a.csv

io = open("c:/data/a.zip", "r")
z = ZipFile.Reader(io)

df = CSV.read(z.files[1])
close(io)


text:
using Test
using CommonSubexpressions
using CommonSubexpressions: binarize

@testset "basic usage" begin
    f_count = 0
    function f(x)
        f_count += 1
        1
    end

    g_count = 0
    function g(x)
        g_count += 1
        2
    end

    @test f_count == 0
    @test g_count == 0

    x = 1
    @cse [f(g(x)) for i in 1:10]
    @test f_count == 1
    @test g_count == 1

    f_count = 0
    g_count = 0
    [f(g(x)) for i in 1:10]
    @test f_count == 10
    @test g_count == 10

    h_count = 0
    function h(x, y)
        h_count += 1
        3
    end

    @cse function foo(x)
        [h(h(f(g(x)), g(x)), i) == i for i in 1:5]
    end
    f_count = 0
    g_count = 0
    h_count = 0
    @test foo(1) == [false, false, true, false, false]
    @test f_count == 1
    @test g_count == 1
    @test h_count == 6
end

@testset "function definition" begin
    f_count = 0
    function f(x)
        f_count += 1
        1
    end

    @cse function foo(x)
        [f(x) == i for i in 1:5]
    end

    @test foo(1) == [true, false, false, false, false]
    @test f_count == 1
end

@testset "algebra" begin
    H = [1 2; 3 4]
    W = [2 3; 4 5]
    G = [4 5; 6 7]

    @test (@cse inv(H) * G + W * (inv(H) + W)) == begin
        invH = inv(H)
        invH * G + W * (inv(H) + W)
    end
end

@testset "dict" begin
    @test (@cse begin
        x = Dict("foo" => sin(pi), "bar" => 2 + 2)
        x["bar"] = 100
        x["foo"], x["bar"]
    end) == begin
        x = Dict("foo" => sin(pi), "bar" => 2 + 2)
        x["bar"] = 100
        x["foo"], x["bar"]
    end
end

@testset "int" begin
    @test (@cse 1) == 1
end

@testset "nested usage" begin
    f_count = 0
    function f(x)
        f_count += 1
        x
    end

    g_count = 0
    function g(x)
        g_count += 1
        2x
    end

    @test f_count == 0
    @test g_count == 0

    result = @cse for i in 1:5
        x = f(2)
        @cse for j in 1:5
            y = g(i)
        end
    end
    @test f_count == 1
    @test g_count == 5

    @test result == for i in 1:5
        x = f(2)
        for j in 1:5
            y = g(i)
        end
    end
end

@testset "cse function" begin
    expr = cse(:((x + 1) * (x + 1)))
    @test expr.head == :block
    @test length(expr.args) == 3
    @test expr.args[1].head == :(=)
    @test expr.args[2].head == :(=)
    @test expr.args[1].args[2].head == :call
    @test expr.args[1].args[2].args == [:(+), :x, 1]
    @test expr.args[2].args[2].head == :call
    @test expr.args[2].args[2].args[1] == :(*)
end

@testset "warnings" begin
    @test_logs (:warn, "CommonSubexpressions can't yet handle expressions of this form: foo") cse(Expr(:foo, 1, 2, 3))

    @cse(1 + 2 + 3, false)
    @cse(1 + 2 + 3, warn=true)
end

@testset "inplace" begin
    x = 1
    @cse begin
        y = 2 + x
        x += y * x
    end
    @test y == 3
    @test x == 4
end

@testset "binarize" begin
    @test binarize(:a) == :a
    @test binarize(:(a)) == :(a)
    @test binarize(:(+a)) == :(+a)
    @test binarize(:(a + b)) == :(a + b)
    @test binarize(:(a + b + c)) == :((a + b) + c)
    @test binarize(:(a + b + c + d)) == :(((a + b) + c) + d)
    @test binarize(:(a + b + c + d + e)) == :((((a + b) + c) + d) + e)

    # Nested function calls
    @test binarize(:((a + b + c) + (d + e + f))) == :(((a + b) + c) + ((d + e) + f))

    # Arbirary binary functions should be left as-is
    @test binarize(:(f(a + b, c))) == :(f(a + b, c))
    @test binarize(:(f(a + b, g(c, d + e)))) == :(f(a + b, g(c, d + e)))

    # Make sure we can binarize functions whose arguments are expressions
    @test binarize(:(f(a + b, c + d, g(e, f)))) == :(f(f(a + b, c + d), g(e, f)))
end

# Test nesting `@cse` with a macro like `@fastmath`
module NestedMacroTest
    using CommonSubexpressions
    using MacroTools: postwalk, @capture
    using Test

    const special_plus_calls = Ref(0)

    function special_plus(args...)
        special_plus_calls[] += 1
        +(args...)
    end

    macro special_math(expr)
        esc(postwalk(expr) do ex
            if @capture(ex, +(x__))
                :(special_plus($(x...)))
            else
                ex
            end
        end)
    end

    @testset "Nested Macros" begin
        @test(@cse(@special_math(2 + 2)) == 4)
        @test special_plus_calls[] == 1

        special_plus_calls[] = 0
        @test(@cse(@binarize(@special_math(1 + 2 + 3 + 4))) == 10)
        # Three more calls to `special_plus`
        @test special_plus_calls[] == 3

        special_plus_calls[] = 0
        # Test that CSE is actually eliminated common subexpressions even
        # when nested with another macro.
        @test(@cse(@special_math((1 + 2) + (3 + 4))) == 10)
        # Three more calls: 1 + 2, 3 + 4, and 3 + 7
        @test special_plus_calls[] == 3

        special_plus_calls[] = 0
        @test(@cse(@special_math((2 + 2) + (2 + 2))) == 8)
        # Only two more calls because of CSE: 2 + 2, 4 + 4
        @test special_plus_calls[] == 2

        # These particular macros also commute:

        special_plus_calls[] = 0
        @test(@special_math(@cse((1 + 2) + (3 + 4))) == 10)
        @test special_plus_calls[] == 3

        special_plus_calls[] = 0
        @test(@special_math(@cse((2 + 2) + (2 + 2))) == 8)
        @test special_plus_calls[] == 2

        special_plus_calls[] = 0
        @test(@cse(@binarize(@special_math((1 + 2 + 3) + (1 + 2 + 4) + (1 + 2 + 5)))) == 21)
        # Test that the duplicate calls to `1 + 2` were eliminated
        @test special_plus_calls[] == 6
    end
end


text:
"""
A glue interaction that will apply to all atom pairs.
Finnis-Sinclair and similar interactions should sub-type this type.
This type should be a GeneralInteraction type. But due to the special
nature of glue interactions and the restriction to pair interactions
of the GeneralInteraction type, glue interactions are for now a sub-type of
SpecificInteraction.
"""
abstract type GlueInteraction <: SpecificInteraction end

struct FinnisSinclairPair{T<:Real}
    # container for two element interaction -> c c0 c1 c2
    c::T
    c₀::T
    c₁::T
    c₂::T
end

FinnisSinclairPair(c,c₀,c₁,c₂) = FinnisSinclairPair{typeof(c₀)}(c,c₀,c₁,c₂)

struct FinnisSinclairSingle{T<:Real}
    # container for single element stuff: glue density, glue energy -> A
    A::T
    β::T
    d::T
end

FinnisSinclairSingle(A,β,d) = FinnisSinclairPair{typeof(A)}(A,β,d)

"""
    FinnisSinclairInteraction(nl_only,pairs,singles,kb)

The Finnis-Sinclair interaction. This interaction expects units to be of
these https://lammps.sandia.gov/doc/units.html units (eV, Å, K, ps and so on).
"""
struct FinnisSinclair{T<:Real} <: GlueInteraction
    nl_only::Bool
    pairs::Dict{String,FinnisSinclairPair{T}}
    singles::Dict{String,FinnisSinclairSingle{T}}
    kb::T
end

FinnisSinclair(nl_only, pairs, singles, kb) = FinnisSinclair{typeof(kb)}(nl_only, pairs, singles, kb)

"""
    get_finnissinclair1984(nl_only)

Finnis and Sinclair 1984 parameterization: https://doi.org/10.1080/01418618408244210
"""
function get_finnissinclair1984(nl_only::Bool)

    elements = ["V", "Nb", "Ta", "Cr", "Mo", "W", "Fe"]

    d = [3.692767, 3.915354, 4.076980, 3.915720,
        4.114825, 4.400224, 3.699579] # (Å)
    A = [2.010637, 3.013789, 2.591061, 1.453418,
        1.887117, 1.896373, 1.889846] # (eV)
    β = [0, 0, 0, 1.8, 0, 0, 1.8] # (1)
    c = [3.8, 4.2, 4.2, 2.9, 3.25, 3.25, 3.4] # (Å)
    c₀ = [-0.8816318, -1.5640104, 1.2157373, 29.1429813,
        43.4475218, 47.1346499, 1.2110601] # (1)
    c₁ = [1.4907756, 2.0055779, 0.0271471, -23.3975027,
        -31.9332978, -33.7665655, -0.7510840] # (1)
    c₂ = [-0.3976370, -0.4663764, -0.1217350, 4.7578297,
        6.0804249, 6.2541999, 0.1380773] # (1)

    fs_singles = Dict()
    fs_pairs = Dict()
    for (i,el) in enumerate(elements)
        fs_singles[el] = FinnisSinclairSingle(A[i],β[i],d[i])
        fs_pairs[string(el,el)] = FinnisSinclairPair(c[i],c₀[i],c₁[i],c₂[i])
    end

    masses = Dict("V" => 50.9415, "Nb" => 92.9064, "Ta" => 180.9479,
                  "Cr" => 51.996, "Mo" => 95.94, "W" => 183.85,
                  "Fe" => 55.847) # g/mole

    # Å
    bcc_lattice_constants = Dict(
        "V" => 3.0399, "Nb" => 3.3008, "Ta" => 3.3058,
        "Cr" => 2.8845, "Mo" => 3.1472, "W" => 3.1652,
        "Fe" => 2.8665
    )

    kb = 8.617333262145e-5 # eV/K
    fs84 = FinnisSinclair(nl_only, fs_pairs, fs_singles, kb)

    u = [5.31, 7.57, 8.1, 4.1, 6.82, 8.9, 4.28] # eV/atom
    u_vac = [1.92, 2.64, 3.13, 1.97, 2.58, 3.71, 1.77] # eV/atom

    reference_energies = Dict(el=>Dict("u"=>u[i], "u_vac"=>u_vac[i])
                              for (i,el) in enumerate(elements))
    return fs84, elements, masses, bcc_lattice_constants, reference_energies
end

"""
    glue(r,β,d)

The core component of the Finnis-Sinclair type GlueInteraction.
Used to calculate contribution to the glue value of an atom, based on
its neighbouring atoms.
"""
function glue(r, β, d)
    return r > d ? 0 : (r-d)^2 + β*(r-d)^3/d
end

"""
    ∂glue_∂r(r, β, d)

Derivative of the glue density function.
"""
∂glue_∂r(r, β, d) = r > d ? 0 : 2*(r-d) + 3*β*(r-d)^2/d

"""
    Uglue(ρ, A)

Energy based on the glue density .
"""
function Uglue(ρ, A)
   return -A * √ρ
end

"""
    ∂Uglue_∂ρ(ρ, A)

Energy derivative given glue density.
"""
∂Uglue_∂ρ(ρ,A) = - A / (2 * √ρ)

"""
    Upair(r, c, c₀, c₁, c₂)

Energy contribution directly from atom pair distances.
"""
function Upair(r, c, c₀, c₁, c₂)
    return (r > c) ? 0 : (r - c)^2 * (c₀ + c₁*r + c₂*r^2)
end

"""
    ∂Upair_∂r(r, c, c₀, c₁, c₂)

Derivative of the pair potential.
"""
∂Upair_∂r(r, c, c₀,c₁, c₂) = (r > c) ? 0 : 2*(r - c) * (c₀ + c₁*r + c₂*r^2) + (r - c)^2 * (c₁ + 2*c₂*r)

"""
    update_glue_densities!(inter, coords, s, parallel)

Convenience function to update the densities before the forces are computed in serial/parallel.
"""
function update_glue_densities!(
        inter::FinnisSinclair,
        coords,
        s::Simulation;
        parallel::Bool=false
    )
    n_atoms = length(s.coords)

    # wiping old glue densities
    for i in 1:n_atoms
        s.glue_densities[i] = 0
    end

    # updating glue densities
    for (i,j) in s.neighbors

        # collecting parameters
        el_i = s.atoms[i].name
        el_j = s.atoms[j].name
        pi = inter.singles[el_i]
        pj = inter.singles[el_j]
        pij = inter.pairs[string(el_i,el_j)]

        # computing distance
        dr = vector(s.coords[i], s.coords[j], s.box_size)
        r = norm(dr)

        # updating glue densities
        s.glue_densities[i] += glue(r, pj.β, pj.d)
        s.glue_densities[j] += glue(r, pi.β, pi.d)
    end
    return nothing
end

@inline @inbounds function force!(forces, inter::FinnisSinclair, s::Simulation, i::Integer, j::Integer)
    fdr = force(inter, s, i, j)
    forces[i] -= fdr
    forces[j] += fdr
    return nothing
end

@inline @inbounds function force(inter::FinnisSinclair, s::Simulation, i::Integer, j::Integer)

    # parameters
    el_i = s.atoms[i].name
    el_j = s.atoms[j].name
    el_ij = string(el_i,el_j)
    pi = inter.singles[el_i]
    pj = inter.singles[el_j]
    pij = inter.pairs[el_ij]

    # distance i -> j
    dr = vector(s.coords[i], s.coords[j], s.box_size)
    r = norm(dr)
    dr = normalize(dr)

    # pair contribution
    f_pair = - ∂Upair_∂r(r, pij.c, pij.c₀, pij.c₁, pij.c₂)

    # glue contribution
    dudρ_i = ∂Uglue_∂ρ(s.glue_densities[i], pi.A)
    dudρ_j = ∂Uglue_∂ρ(s.glue_densities[j], pj.A)
    dΦdr_i = ∂glue_∂r(r, pi.β, pi.d)
    dΦdr_j = ∂glue_∂r(r, pj.β, pj.d)
    f_glue = - (dudρ_j * dΦdr_i + dudρ_i * dΦdr_j)

    # total force contribution
    f = f_pair + f_glue
    return f * dr
end

@inline @inbounds function potential_energy(inter::FinnisSinclair, s::Simulation, i)
    # logger - general inters - computes the glue energy part only for a single atom
    # note: assumes that the glue densities are up to date, only calcutes them when inter.glue_densities are all 0

    # check if densities are zero, if so calculate current, otherwise assume they are current
    no_glue = all(isapprox.(s.glue_densities, zeros(length(s.glue_densities)), atol=1e-4))
    if no_glue
        update_glue_densities!(inter, s.coords, s, parallel=false)
    end

    return Uglue(s.glue_densities[i], inter.singles[s.atoms[i].name].A)
end

@inline @inbounds function potential_energy(inter::FinnisSinclair, s::Simulation, i, j)
    # logger - general inters - computes the pair energy part only for a single atom pair
    dr = vector(s.coords[i], s.coords[j], s.box_size)
    r = norm(dr)
    pij = inter.pairs[string(s.atoms[i].name, s.atoms[j].name)]
    return Upair(r, pij.c, pij.c₀, pij.c₁, pij.c₂)
end


text:
#!/usr/bin/env julia

using Luxor, Colors

function test_pie(fname)
  pagewidth, pageheight = 1200, 1400
  Drawing(pagewidth, pageheight, fname)
  origin() # move 0/0 to center
  background("ivory")
  setopacity(0.8)
  setline(0.4)

  pagetiles = Tiler(pagewidth, pageheight, 4, 5, margin=50)
  rad = 60
  for (pos, n) in pagetiles
    randomhue()
    fromA, toA = rand() * 2pi, rand() * 2pi
    fromA_str, toA_str = string(convert(Int, round(rad2deg(fromA), 0))),
                         string(convert(Int, round(rad2deg(toA), 0)))
    pie(pos, rad, fromA, toA, :fill)

    #label

    text(fromA_str, pos.x + 1.1rad  * cos(fromA), pos.y + 1.1rad  * sin(fromA))
    text(toA_str, pos.x + 1.1rad  * cos(toA), pos.y + 1.1rad  * sin(toA))
    sethue("black")
    text(string(fromA_str, "°/", toA_str, "°"), pos.x, pos.y)
  end
  finish()
  println("pie-test saved in $(fname)")
end

fname = "/tmp/pie-test1.pdf"
test_pie(fname)


text:
export mouseover, mouse_selection, mouseposition, hovered_scene
export select_rectangle, select_line, select_point

"""
    mouseover(scene::SceneLike, plots::AbstractPlot...)

Returns true if the mouse currently hovers any of `plots`.
"""
function mouseover(scene::SceneLike, plots::AbstractPlot...)
    p, idx = mouse_selection(scene)
    return p in flatten_plots(plots)
end

"""
    onpick(f, scene::SceneLike, plots::AbstractPlot...)

Calls `f(idx)` whenever the mouse is over any of `plots`.
`idx` is an index, e.g. when over a scatter plot, it will be the index of the
hovered element
"""
function onpick(f, scene::SceneLike, plots::AbstractPlot...; range=1)
    fplots = flatten_plots(plots)
    args = range == 1 ? (scene,) : (scene, range)
    map_once(events(scene).mouseposition) do mp
        p, idx = mouse_selection(args...)
        (p in fplots) && f(idx)
        return
    end
end

"""
    mouse_selection(scene::Scene)
Returns the plot that is under the current mouse position
"""
function mouse_selection(scene::SceneLike)
    pick(scene, events(scene).mouseposition[])
end
function mouse_selection(scene::SceneLike, range)
    pick(scene, events(scene).mouseposition[], range)
end

function flatten_plots(x::Atomic, plots = AbstractPlot[])
    if isempty(x.plots)
        push!(plots, x)
    else
        flatten_plots(x.plots, plots)
    end
    plots
end

function flatten_plots(x::Combined, plots = AbstractPlot[])
    for elem in x.plots
        flatten_plots(elem, plots)
    end
    plots
end

function flatten_plots(array, plots = AbstractPlot[])
    for elem in array
        flatten_plots(elem, plots)
    end
    plots
end

"""
    mouse_in_scene(scene::Scene)
returns the mouseposition relative to `scene`
"""
function mouse_in_scene(scene::SceneLike)
    p = rootparent(scene)
    lift(pixelarea(p), pixelarea(scene), events(scene).mouseposition) do pa, sa, mp
        Vec(mp) .- minimum(sa)
    end
end

"""
Return the plot under pixel position x y
"""
function pick(scene::SceneLike, x::Number, y::Number)
    return pick(scene, Vec{2, Float64}(x, y))
end

"""
    pick(scene::Scene, xy::VecLike[, range])

Return the plot under pixel position xy
"""
function pick(scene::SceneLike, xy)
    screen = getscreen(scene)
    screen === nothing && return (nothing, 0)
    pick(scene, screen, Vec{2, Float64}(xy))
end
function pick(scene::SceneLike, xy, range)
    screen = getscreen(scene)
    screen === nothing && return (nothing, 0)
    pick(scene, screen, Vec{2, Float64}(xy), Float64(range))
end

"""
Normalizes mouse position relative to the screen rectangle
"""
function screen_relative(scene::Scene, mpos)
    return Point2f0(mpos) .- Point2f0(minimum(pixelarea(scene)[]))
end

"""
    mouseposition(scene = hovered_scene()) -> pos
Return the current position of the mouse `pos` in _data points_ of the
given `scene`.

By default uses the `scene` that the mouse is currently hovering over.
"""
function mouseposition(scene = hovered_scene())
    to_world(scene, mouseposition_px(scene))
end

function mouseposition_px(scene = hovered_scene())
    screen_relative(
        scene,
        events(scene).mouseposition[]
    )
end

"""
    hovered_scene()
Return the `scene` that the mouse is currently hovering over.

Properly identifies the scene for a plot with multiple sub-plots.
"""
hovered_scene() = error("hoevered_scene is not implemented yet.")

"""
    select_rectangle(scene; kwargs...) -> rect
Interactively select a rectangle on a 2D `scene` by clicking the left mouse button,
dragging and then un-clicking. The function returns an **observable** `rect` whose
value corresponds to the selected rectangle on the scene. In addition the function
_plots_ the selected rectangle on the scene as the user clicks and moves the mouse
around. When the button is not clicked any more, the plotted rectangle disappears.

The value of the returned observable is updated **only** when the user un-clicks
(i.e. when the final value of the rectangle has been decided) and only if the
rectangle has area > 0.

The `kwargs...` are propagated into `lines!` which plots the selected rectangle.
"""
function select_rectangle(scene; strokewidth = 3.0, kwargs...)
    key = Mouse.left
    waspressed = Node(false)
    rect = Node(FRect(0, 0, 1, 1)) # plotted rectangle
    rect_ret = Node(FRect(0, 0, 1, 1)) # returned rectangle

    # Create an initially hidden rectangle
    plotted_rect = poly!(
        scene, rect, raw = true, visible = false, color = RGBAf0(0, 0, 0, 0), strokecolor = RGBAf0(0.1, 0.1, 0.8, 0.5), strokewidth = strokewidth, kwargs...,
    )

    on(events(scene).mousedrag) do drag
        if ispressed(scene, key) && is_mouseinside(scene)
            mp = mouseposition(scene)
            if drag == Mouse.down
                waspressed[] = true
                plotted_rect[:visible] = true # start displaying
                rect[] = FRect(mp, 0.0, 0.0)
            elseif drag == Mouse.pressed
                mini = minimum(rect[])
                rect[] = FRect(mini, mp - mini)
            end
        else
            if drag == Mouse.up && waspressed[] # User has selected the rectangle
                waspressed[] = false
                r = absrect(rect[])
                w, h = widths(r)
                if w > 0.0 && h > 0.0 # Ensure that the rectangle has non0 size.
                    rect_ret[] = r
                end
            end
            # always hide if not the right key is pressed
            plotted_rect[:visible] = false # make the plotted rectangle invisible
        end
        return rect_ret
    end
    return rect_ret
end

"""
    select_line(scene; kwargs...) -> line
Interactively select a line (typically an arrow) on a 2D `scene` by clicking the left mouse button,
dragging and then un-clicking. Return an **observable** whose value corresponds
to the selected line on the scene. In addition the function
_plots_ the line on the scene as the user clicks and moves the mouse
around. When the button is not clicked any more, the plotted line disappears.

The value of the returned line is updated **only** when the user un-clicks
and only if the selected line has non-zero length.

The `kwargs...` are propagated into `lines!` which plots the selected line.
"""
function select_line(scene; kwargs...)
    key = Mouse.left
    waspressed = Node(false)
    line = Node([Point2f0(0,0), Point2f0(1,1)])
    line_ret = Node([Point2f0(0,0), Point2f0(1,1)])
    # Create an initially hidden  arrow
    plotted_line = lines!(
        scene, line; visible = false, color = RGBAf0(0.1, 0.1, 0.8, 0.5),
        linewidth = 4, kwargs...,
    )

    on(events(scene).mousedrag) do drag
        if ispressed(scene, key) && is_mouseinside(scene)
            mp = mouseposition(scene)
            if drag == Mouse.down
                waspressed[] = true
                plotted_line[:visible] = true  # start displaying
                line[][1] = mp
                line[][2] = mp
                line[] = line[]
            elseif drag == Mouse.pressed
                line[][2] = mp
                line[] = line[] # actually update observable
            end
        else
            if drag == Mouse.up && waspressed[] # User has selected the rectangle
                waspressed[] = false
                if line[][1] != line[][2]
                    line_ret[] = copy(line[])
                end
            end
            # always hide if not the right key is pressed
            plotted_line[:visible] = false
        end
        return line_ret
    end
    return line_ret
end

"""
    select_point(scene; kwargs...) -> point
Interactively select a point on a 2D `scene` by clicking the left mouse button,
dragging and then un-clicking. Return an **observable** whose value corresponds
to the selected point on the scene. In addition the function
_plots_ the point on the scene as the user clicks and moves the mouse
around. When the button is not clicked any more, the plotted point disappears.

The value of the returned point is updated **only** when the user un-clicks.

The `kwargs...` are propagated into `scatter!` which plots the selected point.
"""
function select_point(scene; kwargs...)
    key = Mouse.left
    pmarker = Circle(Point2f0(0, 0), Float32(1))
    waspressed = Node(false)
    point = Node([Point2f0(0,0)])
    point_ret = Node(Point2f0(0,0))
    # Create an initially hidden  arrow
    plotted_point = scatter!(
        scene, point; visible = false, marker = pmarker, markersize = 20px,
        color = RGBAf0(0.1, 0.1, 0.8, 0.5), kwargs...,
    )

    on(events(scene).mousedrag) do drag
        if ispressed(scene, key) && is_mouseinside(scene)
            mp = mouseposition(scene)
            if drag == Mouse.down
                waspressed[] = true
                plotted_point[:visible] = true  # start displaying
                point[][1] = mp
                point[] = point[]
            elseif drag == Mouse.pressed
                point[][1] = mp
                point[] = point[] # actually update observable
            end
        else
            if drag == Mouse.up && waspressed[] # User has selected the rectangle
                waspressed[] = false
                point_ret[] = copy(point[][1])
            end
            # always hide if not the right key is pressed
            plotted_point[:visible] = false
        end
        return point_ret
    end
    return point_ret
end


text:
# ## Data

import CSV
import TuringModels

using DataFrames
using StatsFuns

data_path = joinpath(TuringModels.project_root, "data", "chimpanzees.csv")
df = CSV.read(data_path, DataFrame; delim=';')
first(df, 10)

# ## Model

using Turing

@model m10_4(y, actors, x₁, x₂) = begin
    ## Number of unique actors in the data set
    N_actor = length(unique(actors))

    ## Set an TArray for the priors/param
    α ~ filldist(Normal(0, 10), N_actor)
    βp ~ Normal(0, 10)
    βpC ~ Normal(0, 10)

    logits = α[actors] .+ (βp .+ βpC * x₁) .* x₂
    y .~ BinomialLogit.(1, logits)
end

model = m10_4(df.pulled_left, df.actor, df.condition, df.prosoc_left);

# ## Output

chains = sample(model, NUTS(0.65), 1000)

# \defaultoutput{}

# ## Original output

m_10_04s_result = "
Iterations = 1:1000
Thinning interval = 1
Chains = 1,2,3,4
Samples per chain = 1000

Empirical Posterior Estimates:
        Mean        SD       Naive SE       MCSE      ESS
a.1 -0.74503184 0.26613979 0.0042080396 0.0060183398 1000
a.2 10.77955494 5.32538998 0.0842018089 0.1269148045 1000
a.3 -1.04982353 0.28535997 0.0045119373 0.0049074219 1000
a.4 -1.04898135 0.28129307 0.0044476339 0.0056325117 1000
a.5 -0.74390933 0.26949936 0.0042611590 0.0052178124 1000
a.6  0.21599365 0.26307574 0.0041595927 0.0045153523 1000
a.7  1.81090866 0.39318577 0.0062168129 0.0071483527 1000
bp  0.83979926 0.26284676 0.0041559722 0.0059795826 1000
bpC -0.12913322 0.29935741 0.0047332562 0.0049519863 1000
";


text:
# using DataStructures

function readinput(filename)
  lines = readlines(filename)
  M = zeros(Int, length(lines),length(lines[1]))
  for (i,l) in enumerate(lines)
    for (j,c) in enumerate(split(l,""))
      M[i,j] = parse(Int,String(c))
    end
  end
  M
end

function expanded_vals(n)
  S = [
    0 1 2 3 4
    1 2 3 4 5
    2 3 4 5 6
    3 4 5 6 7
    4 5 6 7 8
  ]
  N = zeros(Int, 5, 5)
  fill!(N,n)
  (x -> x>9 ? x-9 : x).(N + S)
end

function expandedmat(M)
  H = size(M,1)
  W = size(M,2)
  M2 = zeros(Int, H*5, W*5)
  for x in 1:H, y in 1:W
    M2[x:H:end,y:W:end] = expanded_vals(M[x,y])
  end
  M2
end

function cheapestpath(M)
  q = [(0,1,1)]
  costs = Dict()
  while true
    cost,x,y = q[1]
    if x==size(M,1) && y==size(M,2)
      return cost
      break
    end
    popfirst!(q)
    for (dx,dy) ∈ [(-1,0), (1,0), (0,1), (0,-1)]
      if x+dx in 1:size(M,1) && y+dy in 1:size(M,2)
        ncost = cost + M[x+dx,y+dy]
        if (x+dx,y+dy) ∈ keys(costs) && costs[(x+dx,y+dy)] <= ncost
          continue
        end
        costs[(x+dx,y+dy)] = ncost
        push!(q,(ncost,x+dx,y+dy))
      end
    end
    q = sort(q)
  end
end

function cheapestpath2(M,scalar)
  H = size(M,1)
  W = size(M,2)
  q = [(0,1,1)]
  costs = Dict()
  while true
    cost,x,y = q[1]
    if x==H*scalar && y==W*scalar
      return cost
    end
    popfirst!(q)
    # val = Msmall[mod(x,size(Msmall,1)),mod(y,size(Msmall,2))] + trunc(Int,x/size(Msmall,1)) + trunc(Int,y/size(Msmall,2))
    @show mod(x,)
    for (dx,dy) ∈ [(-1,0), (1,0), (0,1), (0,-1)]
      if x+dx in 1:H*scalar && y+dy in 1:W*scalar
        ncost = cost + M[x+dx,y+dy]
        if (x+dx,y+dy) ∈ keys(costs) && costs[(x+dx,y+dy)] <= ncost
          continue
        end
        costs[(x+dx,y+dy)] = ncost
        push!(q,(ncost,x+dx,y+dy))
      end
    end
    q = sort(q)
  end
end

M = readinput("inputs/11-15/15ex.txt")
# @show cheapestpath(M)

M2 = expandedmat(M)
# @show cheapestpath2(M2,M)
@show cheapestpath2(M,1)

using BenchmarkTools
# @btime cheapestpath(M)


text:
time_start = time()

# NOTE: requires JuMP v0.23.1, ComplexOptInterface v0.1.1
using PowerModels
using Ipopt
using JuMP
using ComplexOptInterface

pkg_load_time = time() - time_start

time_start = time()

file_name = "../data/pglib_opf_case5_pjm.m"
#file_name = "../data/pglib_opf_case118_ieee.m"

data = PowerModels.parse_file(file_name)
PowerModels.standardize_cost_terms!(data, order=2)
PowerModels.calc_thermal_limits!(data)
ref = PowerModels.build_ref(data)[:it][:pm][:nw][0]

data_load_time = time() - time_start

time_start = time()

model = Model(Ipopt.Optimizer)
#set_optimizer_attribute(model, "print_level", 0)

ComplexOptInterface.add_all_bridges(model)

@variable(model, -ref[:bus][i]["vmax"] <= vr[i in keys(ref[:bus])] <= ref[:bus][i]["vmax"], start=1.0)
@variable(model, -ref[:bus][i]["vmax"] <= vi[i in keys(ref[:bus])] <= ref[:bus][i]["vmax"], start=0.0)
V = Dict(i => vr[i] + vi[i]*im for i in keys(ref[:bus]))

@variable(model, ref[:gen][i]["pmin"] <= pg[i in keys(ref[:gen])] <= ref[:gen][i]["pmax"])
@variable(model, ref[:gen][i]["qmin"] <= qg[i in keys(ref[:gen])] <= ref[:gen][i]["qmax"])
G = Dict(i => pg[i] + qg[i]*im for i in keys(ref[:gen]))

@variable(model, -ref[:branch][l]["rate_a"] <= p[(l,i,j) in ref[:arcs]] <= ref[:branch][l]["rate_a"])
@variable(model, -ref[:branch][l]["rate_a"] <= q[(l,i,j) in ref[:arcs]] <= ref[:branch][l]["rate_a"])
S = Dict(k => p[k] + q[k]*im for k in ref[:arcs])

@objective(model, Min, sum(gen["cost"][1]*pg[i]^2 + gen["cost"][2]*pg[i] + gen["cost"][3] for (i,gen) in ref[:gen]))

for (i,bus) in ref[:ref_buses]
    @constraint(model, vi[i] == 0)
end

for (i,bus) in ref[:bus]
    bus_loads = [ref[:load][l] for l in ref[:bus_loads][i]]
    bus_shunts = [ref[:shunt][s] for s in ref[:bus_shunts][i]]

    @constraint(model, abs2(V[i]) >= bus["vmin"]^2)
    @constraint(model, abs2(V[i]) <= bus["vmax"]^2)

    @constraint(model,
        sum(S[a] for a in ref[:bus_arcs][i]) ==
        sum(G[g] for g in ref[:bus_gens][i]) -
        sum(load["pd"] + load["qd"]*im for load in bus_loads) -
        sum(shunt["gs"] - shunt["bs"]*im for shunt in bus_shunts)*abs2(V[i])
    )
end

# Branch power flow physics and limit constraints
for (i,branch) in ref[:branch]
    f_idx = (i, branch["f_bus"], branch["t_bus"])
    t_idx = (i, branch["t_bus"], branch["f_bus"])

    S_fr = S[f_idx]
    S_to = S[t_idx]

    V_fr = V[branch["f_bus"]]
    V_to = V[branch["t_bus"]]

    g, b = PowerModels.calc_branch_y(branch)
    Y = g + b*im

    tr, ti = PowerModels.calc_branch_t(branch)
    T = tr + ti*im

    Y_fr = branch["g_fr"] + branch["b_fr"]*im
    Y_to = branch["g_to"] + branch["b_to"]*im

    # Branch Flow
    @constraint(model, S_fr == conj(Y+Y_fr)/(T*conj(T))*V_fr*conj(V_fr) - conj(Y)/T*V_fr*conj(V_to))
    @constraint(model, S_to == conj(Y+Y_to)*V_to*conj(V_to) - conj(Y)/conj(T)*conj(V_fr)*V_to)

    # Voltage angle difference limit
    W = V_fr*conj(V_to)
    @constraint(model, imag(W) <= tan(branch["angmax"])*real(W))
    @constraint(model, imag(W) >= tan(branch["angmin"])*real(W))

    # Apparent power limit, from side and to side
    @constraint(model, abs2(S_fr) <= branch["rate_a"]^2)
    @constraint(model, abs2(S_to) <= branch["rate_a"]^2)
end

model_build_time = time() - time_start

#println("model building done")
#println(model)

time_start = time()

optimize!(model)
cost = objective_value(model)

solve_time = time() - time_start

println("")
println("\033[1mSummary\033[0m")
println("   case......: $(file_name)")
println("   cost......: $(round(Int, cost))")
println("   pkg time..: $(pkg_load_time)")
println("   data time.: $(data_load_time)")
println("   build time: $(model_build_time)")
println("   solve time: $(solve_time)")
println("")

max_stars_repo_name: lanl-ansi/rosetta-opf

max_issues_repo_name: lanl-ansi/rosetta-opf

max_forks_repo_name: lanl-ansi/rosetta-opf


text:
a = collect(3:3:30)
@test MultiScaleArrays.bisect_search(a,20) == 7
@test MultiScaleArrays.bisect_search(a,13) == 5
@test MultiScaleArrays.bisect_search(a,12) == 4
for i = 1:30
  @test MultiScaleArrays.bisect_search(a,i) == ((i-1)÷3)+1 #+1 for 1-based indexing
end


text:
# test_direct_discretization.jl

# direct discretization for elliptic and anisotropic problem in 2D/3D
function get_problem_2d_elliptic(n)
    A = NotSoSimpleMultigrid.MultigridMethod((n, m)->elliptic2d(f(n, m)), (n, n), V(3, 2))
    b = fill(1,(n-1)*(n-1))
    return (A,b)
end

@testset "Direct discretization of elliptic PDE, 2D" begin
    for n in [16 32 64 128 256 512]
        (A, b) = get_problem_2d_elliptic(n)
        A.grids[1].b .= b # copy rhs
        for item in Base.Iterators.take(A,15) end # iterate
        @test A.resnorm[end] < 1/n^2
        log(A,2)
    end
end

function get_problem_2d_anisotropic(n)
    A = NotSoSimpleMultigrid.MultigridMethod((n,m)->anisotropic2d(1e-8,n,m), (n, n), V(3, 2))
    b = fill(1,(n-1)*(n-1))
    return (A,b)
end

@testset "Direct discretization of anisotropic PDE, 2D" begin
    for n in [16 32 64 128 256 512]
        (A, b) = get_problem_2d_anisotropic(n)
        A.grids[1].b .= b # copy rhs
        for item in Base.Iterators.take(A,15) end # iterate
        @test A.resnorm[end] < 1/n^2
        log(A,2)
    end
end

function get_problem_2d_rotated_anisotropic(n)
    A = NotSoSimpleMultigrid.MultigridMethod((n,m)->anisotropic2d(1e-4,10*π/180,n,m), (n, n), V(3, 2))
    b = fill(1,(n-1)*(n-1))
    return (A,b)
end

@testset "Direct discretization of rotated anisotropic PDE, 2D" begin
    for n in [16 32 64 128 256 512]
        (A, b) = get_problem_2d_anisotropic(n)
        A.grids[1].b .= b # copy rhs
        for item in Base.Iterators.take(A,15) end # iterate
        @test A.resnorm[end] < 1/n^2
        log(A,2)
    end
end

function get_problem_3d_elliptic(n)
    A = NotSoSimpleMultigrid.MultigridMethod((n,m,k)->elliptic3d(f(n,m,k)), (n, n, n), V(3, 2))
    b = fill(1,(n-1)*(n-1)*(n-1))
    return (A,b)
end

@testset "Direct discretization of elliptic PDE, 3D" begin
    for n in [4 8 16 32 64]
        (A, b) = get_problem_3d_elliptic(n)
        A.grids[1].b .= b # copy rhs
        for item in Base.Iterators.take(A,15) end # iterate
        @test A.resnorm[end] < 1/n^3
        log(A,3)
    end
end

function get_problem_3d_anisotropic(n)
    A = NotSoSimpleMultigrid.MultigridMethod((n,m,k)->anisotropic3d(1e-8, 1, n, m, k), (n, n, n), V(3, 2))
    b = fill(1,(n-1)*(n-1)*(n-1))
    return (A,b)
end

@testset "Direct discretization of anisotropic PDE, 3D" begin
    for n in [4 8 16 32 64]
        (A, b) = get_problem_3d_anisotropic(n)
        A.grids[1].b .= b # copy rhs
        for item in Base.Iterators.take(A,15) end # iterate
        @test A.resnorm[end] < 1/n^3
        log(A,3)
    end
end


text:
using Revise
using LinearAlgebra
using Plots

includet("src/InternalArray.jl")
includet("src/FixedStepMethods.jl")
includet("src/MixedStepMethods.jl")
includet("src/Interpolation.jl")
includet("src/Arenstorf.jl")
IArray = InternalArray

const target_err = 1e-10
const closing_err = 1e-5
const stop_time = 40
const min_step = 1e-7
const max_step = 1.0
extra_info = Dict()

(t, x), time = @timed mixed_step_ode_solver(
    Arenstorf.vderivative,
    0,
    [Arenstorf.p0; Arenstorf.v0],
    (t, x) ->
        last(t) > 20 && norm(Arenstorf.p0 .- last(x)[1:2]) < closing_err || last(t) > stop_time,
    min_step = min_step,
    max_step = max_step,
    target_err = target_err,
    method = RK4_step,
    maxit = 5*10^6,
    debug_dict = extra_info
)

steps = [t[i]-t[i-1] for i in 2:length(t)]

p1 = Arenstorf.init_plot()
plot!(
    IArray(x, 1),
    IArray(x, 2),
    label = "satellite",
    title = "Arenstorf Orbit - RK4-Richardson",
)
p2 = plot(
    t[2:length(t)],
    steps,
    label = "Step Size",
    title = "Variation of step",
    legend = :bottomright,
)
plot(p1, p2)
savefig("media/arenstorf/richardson.png")

begin
    print("\033[32;1marenstorf>\033[0m \n")
    print("\tMethod: Richardson of RK4\n")
    print("\tMin. Step: $(min_step)\n")
    print("\tMax. Step: $(max_step)\n")
    print("\tTolerance: $(target_err)\n")
    print("\tTarget Error (closing condition): $(closing_err)\n")
    println()
    numit = extra_info[:numit]
    print("\tNumber of iterations: $(numit)\n")
    print("\tRunning time: $(time)s\n")
    print("\tTime (solution): $(last(t))\n")
    print("\tNumber of evaluations: $(11*numit)\n")
    print("\tOrbit closed with target tolerance: $(last(t) < stop_time)\n")
    print("\tDistance from last point to start point: ",
          "$(norm(last(x)[1:2] .- Arenstorf.p0))\n")
    n = length(x)
    last5x = [p[1] for p in x[n-4:n]]
    last5y = [p[2] for p in x[n-4:n]]
    closingx = Interpolation.hermite(last5y, last5x, 0)
    print("\tDistance from interpolated closing point to start point: ",
          "$(abs(closingx - Arenstorf.p0[1]))\n")
end

# Creation of Gif

using IJulia

len = length(t)÷2
p1 = plot(
    title = "Arenstorf orbit",
    aspect_ratio = :equal,
    framestyle = :box,
    xlims = (-1.4, 1.2),
    ylims = (-1.2, 1.2),
    legend = :bottomright,
)
plot!(p1, [x[1][1]], [x[1][2]], label = "satellite")
plot!(p1, [cos(0)], [sin(0)], color = :black, label = "moon")
scatter!((0, 0), label = "", color = :black, markersize = 1.5)
p2 = plot(
    title = "Variation of step",
    xlims = (0, t[len]),
    ylims = (0, maximum(steps[1:len-1])*1.2),
    legend = :bottomright,
)
plot!(p2, [0], [0], label = "step size")

anim = Animation()
i = 2
for j = 0:0.05:t[len]
    while i < len && t[i] < j
        i += 1
    end
    push!(p1, [(x[i][1], x[i][2]), (cos(2*π*t[i]/t[len]), sin(2*π*t[i]/t[len])), (NaN, NaN)])
    push!(p2, t[i], steps[i-1])
    plot(p1, p2)
    frame(anim)
end
gif(anim, "media/arenstorf/richardson.gif", fps = 15)

max_stars_repo_name: useredsa/numeric-differential-equations

max_issues_repo_name: useredsa/numeric-differential-equations

max_forks_repo_name: useredsa/numeric-differential-equations


text:
using LightningArt
using Luxor
using LinearAlgebra
using ForwardDiff
using DrWatson
@quickactivate "LightningArt"
import GeometryBasics
GB = GeometryBasics
using LightGraphs, MetaGraphs
include(scriptsdir("guitar_model.jl"))

@svg begin
    circle(Luxor.Point(0, 0), 200, :stroke)
end

@svg begin
    line(Luxor.Point(0, 0), Luxor.Point(100,100), :stroke)
end

graph = DiGraph()
add_vertex!(graph)
collect(vertices(graph))

ggplot(graph) = gplot(graph, nodelabel=1:nv(graph), edgelabel=1:ne(graph))
add_edge!(graph, 1, 2)

mg = MetaDiGraph(graph, 1.0)

abstract type AbstractCurve end
struct BasicCurve end
struct FunctionCurve{F} <: AbstractCurve
    curve::F
end
(curve::FunctionCurve)(t) = curve.curve(t)

add_vertex!(mg)
ggplot(mg)
set_props!(mg, Edge(1,2), Dict{Symbol, BasicCurve}(:curve=>BasicCurve()))
props(mg, Edge(1,2))
props(mg, Edge(1,3))
add_vertex!(mg)
add_edge!(mg, 1, 3)
ggplot(mg)

collect(edges(mg))
set_prop!(mg, 1, :position, Point2(0.0, -250.0))
set_prop!(mg, 2, :position, Point2(0.0, 10.0))
set_prop!(mg, 3, :position, Point2(10.0, 10.0))

function render_graph(metagraph; render_vertices=true, svg=true)
    if svg
        @svg begin
            position_array = get_prop(metagraph, :position)
            if render_vertices
                for vert in vertices(metagraph)
                    circle(Luxor.Point(position_array[vert]...), 0.1, :stroke)
                end
            end
            for edge in edges(metagraph)
                parent_pos = position_array[edge.src]
                child_pos = position_array[edge.dst]
                line(Luxor.Point(child_pos...), Luxor.Point(parent_pos...), :stroke)
            end
        end
    else
        @png begin
            position_array = get_prop(metagraph, :position)
            if render_vertices
                for vert in vertices(metagraph)
                    circle(Luxor.Point(position_array[vert]...), 0.1, :stroke)
                end
            end
            for edge in edges(metagraph)
                parent_pos = position_array[edge.src]
                child_pos = position_array[edge.dst]
                line(Luxor.Point(child_pos...), Luxor.Point(parent_pos...), :stroke)
            end
        end
    end
end

render_graph(mg)
set_prop!(mg, :position, [GeometryBasics.Point2(0.0, 0.0)])

normsq(v) = dot(v, v)

function create_diffeq_graph(n, initial_array, initial_gradients, potential_function, α, β, γ; expansion_rate=0.05, interior_func=x->true)
    ∇potential_function(p) = ForwardDiff.gradient(potential_function, p)
    metagraph = MetaDiGraph(n)
    position_array = initial_array
    potential_array = potential_function.(initial_array)
    ∇potential_array = initial_gradients
    ∇sign_array = [1.0, -1.0]
    set_prop!(metagraph, :position, position_array)
    topleft_corner = GeometryBasics.Point2(-300.0, -300.0)
    bottomright_corner = GeometryBasics.Point2(300.0, 300.0)
    difference_width = bottomright_corner .- topleft_corner
    for child_vertex in length(initial_array)+1:n
        # random parent, not very good
        #parent_vertex = rand(1:child_vertex-1)  # TODO: change this for more variety in sampling

        tries = 0
        while true
            target_position = GeometryBasics.Point2(rand(2)) .* difference_width .+ topleft_corner
            target_potential = potential_function(target_position)

            #@show target_potential, potential_array[2]
            #@show target_potential .- potential_array[2]
            distance_function(j) = α .* norm(target_position .- position_array[j]) .+ β .* norm(target_potential .- potential_array[j])

            parent_vertex = 1
            if child_vertex == 3
                parent_vertex = 1
            elseif child_vertex == 4
                parent_vertex = 2
            else
                parent_distance, parent_vertex = findmin(map(distance_function, 1:child_vertex-1))
            end

            gradient = ∇potential_array[parent_vertex]
            #@show parent_vertex
            #@show gradient

            parent_position = position_array[parent_vertex]
            child_position = expansion_rate .* normalize(normalize(target_position .- parent_position) .+ γ .* ∇sign_array[parent_vertex] .* normalize(gradient)) .+ parent_position
            if interior_func((child_position .- topleft_corner) ./ difference_width)
                push!(position_array, child_position)
                push!(potential_array, potential_function(child_position))
                push!(∇potential_array, ∇potential_function(child_position))
                push!(∇sign_array, ∇sign_array[parent_vertex])

                add_edge!(metagraph, parent_vertex, child_vertex, :curve, FunctionCurve(t->t * (child_position .- parent_position) .* t .+ parent_position))
                break
            end
            tries += 1
            if tries > 100
                throw(ErrorException("too many tries to find a connection :("))
            end
        end

    end
    metagraph

end

function top_bottom_proto_rrt(n; expansion_rate=5.0)
    top_middle_point = GeometryBasics.Point2(0.0, -220.0)
    bottom_middle_point = GeometryBasics.Point2(0.0, 220.0)
    initial_array = [top_middle_point, bottom_middle_point]
    initial_gradients = [GeometryBasics.Point2(0.0, 1.0), GeometryBasics.Point2(0.0, -1.0)]
    create_diffeq_graph(n, initial_array, initial_gradients, identity, 1.0, 0.0; expansion_rate=expansion_rate)
end

function top_bottom_proto_rrt_potential(n; expansion_rate=5.0)
    charge_constant = 100.0
    potential_function(p) = charge_constant .* (-1.0 ./ (norm(p .- top_middle_point) + 1e-1) .+ 1.0 ./ (norm(p .- bottom_middle_point) + 1e-1))
    top_middle_point = GeometryBasics.Point2(0.0, -200.0)
    bottom_middle_point = GeometryBasics.Point2(0.0, 200.0)
    initial_array = [top_middle_point, bottom_middle_point]
    initial_gradients = [GeometryBasics.Point2(0.0, 1.0), GeometryBasics.Point2(0.0, 1.0)]
    create_diffeq_graph(n, initial_array, initial_gradients, potential_function, 0.3, 1.0, 1.2; expansion_rate=expansion_rate)
end

middle_point = GeometryBasics.Point2(0.0,245.0)

mg = top_bottom_proto_rrt(3;expansion_rate=1.5)
mg = top_bottom_proto_rrt_potential(10000;expansion_rate=1.0)
render_graph(mg; render_vertices=true, svg=false)

function guitar_proto_rrt_potential(n, α, β, γ; expansion_rate=1.5)
    charge_constant = 100.0
    topleft_corner = GB.Point2(-300.0, -300.0)
    bottomright_corner = GB.Point2(300.0, 300.0)
    left_point = GB.Point2(0.5, 0.2) .* (bottomright_corner .- topleft_corner) .+ topleft_corner
    right_point = GB.Point2(0.25, 0.66) .* (bottomright_corner .- topleft_corner) .+ topleft_corner
    potential_function(p) = charge_constant .* (-1.0 ./ (norm(p .- left_point) + 1e-1) .+ 1.0 ./ (norm(p .- right_point) + 1e-1))
    initial_array = [left_point, right_point]
    initial_gradients = [GeometryBasics.Point2(0.0, 1.0), GeometryBasics.Point2(0.0, 1.0)]
    interior_func = inside_guitar_function
    create_diffeq_graph(n, initial_array, initial_gradients, potential_function, α, β, γ; expansion_rate=expansion_rate, interior_func)
end

mg = guitar_proto_rrt_potential(100000, 0.3, 1.0, 1.2;expansion_rate=1.0)
render_graph(mg; render_vertices=true, svg=true)

max_stars_repo_name: zoemcc/LightningArt

max_issues_repo_name: zoemcc/LightningArt

max_forks_repo_name: zoemcc/LightningArt


text:
"""
    UtilFuncs

This module collects many conventional utility functions
which are provided within the general form ``u(X|\\theta):\\mathbb{R}^{k}\\to\\mathbb{R}``.
where ``\\theta`` is `tuple` of parameters (in specific types of utility function, they are defined as keyword parameters), and ``k`` is the dimension of utilizable goods.

There is the contents of utility functions:
1. `u_log` - logarithm (risk netural, 1 good)
2. `u_crra` - CRRA
3. `u_cara` - CARA
4. `u_ces2` - CES (2 goods)
5. `u_cesk` - CES (k goods)
6. `u_cd1` - Cobb-Douglas (1 good)
7. `u_cd2` - Cobb-Douglas (2 goods)
8. `u_cdk` - Cobb-Douglas (k goods)
9. `u_logcd2` - Logarithm Cobb-Douglas (2 goods)
10. `u_logcdk` - Logarithm Cobb-Douglas (k goods)
11. `u_leontiefk` - Leontief (k goods)

Please refer to each specific function for detailed documentations.

This module depends on:
1. nothing

@author: Clpr @ GitHub
"""
module UtilFuncs
    export u_log, u_crra, u_cara
    export u_ces2, u_cesk
    export u_cd1, u_cd2, u_cdk
    export u_logcd2, u_logcdk
    export u_leontiefk
# =================================
"""
    u_log(x::Real)

Logarithm (risk-netural) utility function for one good
where `x` is naturally required to be positive.
The expression of logarithm utility function is:
``u(x)=\\log(x),x>0``
"""
function u_log(x::Real)
    return log(x)::Real
end # u_log

# ----------------
"""
    u_crra(x::Real ; η::Real)

CRRA (constant relative risk aversion) utility function of one good
whose expression is:
``u(x|\\eta>0) = \\frac{x^{1-\\eta}}{1-\\eta}``.
The relative risk aversion coefficient ``\\eta`` must be positive.
When ``\\eta`` is one, this function degenerates to logarithm utility function.
"""
function u_crra(x::Real ; η::Real = 1.1)
    if η<0; throw(ErrorException("requires positive eta (RRA)")); end
    η == 1 ? (return u_log(x)::Real)  :  (return (x^(1-η)/(1-η))::Real)
end # u_crra

# ----------------
"""
    u_cara(x::Real ; α::Real = 1.0)

CARA (constant absolute risk aversion) utility function of one good.
whose expression is:
``u(x|\\alpha>0) = -e^{-\\alpha x}``.
The absolute risk aversion coefficient ``\\alpha`` must be positive.
"""
function u_cara(x::Real ; α::Real = 1.0)
    if α<0; throw(ErrorException("requires positive alpha (ARA)")); end
    return (-exp(-α * x))::Real
end # u_cara

# ----------------
"""
    u_ces2(x1::Real, x2::Real ; ρ::Real = 1.0, β::Real = 1.0, α::Real = 0.5 )

CES (constant elasticity of substitution) utility function whose expression is:
``u(x_1,x_2)=(\\alpha x_1^{\\rho} + \\(1-alpha) x_2^{\\rho})^{\\beta/\\rho}``
The elasticity of substitution is ``\\frac{1}{1-\\rho}``.

Special cases:
1. when ``\\rho=1``, the function has linear/perfect substitution.
2. when ``\\rho=0``, the function becomes Cobb-Douglas function of two goods, i.e. ``u(x_1,x_2)=(x_1^\\alpha x_2^{1-\\alpha})^{\\beta}``
3. when ``\\rho=-\\infty``, the function becomes Leontief (perfect complimentary), i.e. ``u(x_1,x_2)=\\beta \\min\\{ x_1,x_2 \\}``
"""
function u_ces2(x1::Real, x2::Real ; ρ::Real = 1.0, β::Real = 1.0, α::Real = 0.5 )
    if ρ == 0
        return ( ( x1^α * x2^(1-α) )^β )::Real
    elseif ρ == -Inf
        return ( β * min(x1,x2) )::Real
    end # if
    return ( ( α*x1^ρ + (1-α)*x2^ρ )^(β/ρ) )::Real
end # u_ces2

# ----------------
"""
    u_cesk(X::Vector ; W::Vector = ones(length(X)) , β::Real = 1.0, ρ::Real = 1.0 )

CES utility function of k goods whose expression is:
``u(X)=(\\sum w_i \\log(x_i)^{\\rho})^{\\beta/\\rho}``.
requiring the length of `W` must be equal to the length of `X`.
This function does not require the summation of `W` must be 1.

Special cases:
1. when ``\\rho=1``, the function has linear/perfect substitution.
2. when ``\\rho=0``, the function becomes Cobb-Douglas function
3. when ``\\rho=-\\infty``, the function becomes Leontief (perfect complimentary)

"""
function u_cesk(X::Vector ; W::Vector = ones(length(X)) , β::Real = 1.0, ρ::Real = 1.0 )
    if ρ == 0
        return ( ( u_cdk(X, W=W) )^β )::Real
    elseif ρ == -Inf
        return ( β * u_leontiefk(X) )::Real
    end # if
    local U::Real = 0.0
    for (x,w) in zip(X,W)
        U += w * x ^ ρ
    end # for
    return ( U^(β/ρ) )::Real
end # u_cesk

# ----------------
"""
    u_cd1(x::Real ; θ::Real = 1.0)

Cobb-Douglas utility function of one good whose expression is:
``u(x)=x^\\theta``
"""
function u_cd1(x::Real ; θ::Real = 1.0)
    return (x^θ)::Real
end # u_cd1

# ----------------
"""
    u_cd2(x::Real ; θ::Real = 0.5)

Cobb-Douglas utility function of two good whose expression is:
``u(x_1,x_2)=x_1^\\theta\\cdot x_2^{1-\\theta}``
"""
function u_cd2(x1::Real, x2::Real ; θ::Real = 0.5)
    return (x1^θ * x2^(1-θ))::Real
end # u_cd2

# ----------------
"""
    u_cdk(X::Vector ; W::Vector = ones(length(X)))

Cobb-Douglas utility function of k goods whose expression is:
``u(X)=\\prod \\log(x_i)^{w_i}``.
requiring the length of `W` must be equal to the length of `X`.
This function does not require the summation of `W` must be 1.
"""
function u_cdk(X::Vector ; W::Vector = ones(length(X)) )
    return exp(u_logcdk(X, W = W))::Real
end # u_cdk

# ----------------
"""
    u_logcd2(x1::Real, x2::Real ; β::Real = 1.0 )

Logarithm Cobb-Douglas utility function of k goods whose expression is:
``u(x_1,x_2)=\\log(x_1)+\\beta\\log(x_2)``
"""
function u_logcd2(x1::Real, x2::Real ; β::Real = 1.0 )
    return ( log(x1) + β * log(x2) )::Real
end

# ----------------
"""
    u_logcdk(X::Vector ; W::Vector = ones(length(X)))

Logarithm Cobb-Douglas utility function whose expression is:
``u(X)=\\sum w_i\\log(x_i)``.
requiring the length of `W` must be equal to the length of `X`.
"""
function u_logcdk(X::Vector ; W::Vector = ones(length(X)) )
    local U::Real = 0.0
    for (x,w) in zip(X,W); U += w * x; end # for
    return U::Real
end # u_logcdk

# ----------------
"""
    u_leontiefk(X::Vector)

Leontief (perfect complimentary) utility function of k goods:
``u(\\mathbf{x})=\\min\\{x_1,x_2,\\dots\\}``
"""
function u_leontiefk(X::Vector)
    return findmin(X)
end # u_leontiefk

# =================================
end # UtilFuncs

max_stars_repo_name: Clpr/EconToolBox

max_issues_repo_name: Clpr/EconToolBox

max_forks_repo_name: Clpr/EconToolBox


text:
# Performs tests on a checker considering a bunch of different transforms of the
# needle and haystack
function test_checker(
    check::Function,
    haystack::AbstractString,
    right_needle,
    wrong_needle,
)
    @testset "$check" begin
        for (haystack_f, needle_f) in haystack_needle_transforms
            test_checker1(
                check,
                transform(haystack_f, haystack),
                transform(needle_f, right_needle),
                transform(needle_f, wrong_needle),
            )
        end
    end
end

const haystack_needle_transforms = Iterators.product(
    (identity, Symbol),  # haystack
    (identity, Symbol, first, Regex),  # needle
)

transform(f, x) = f(x)
transform(f, x::Tuple) = map(f, x)  # for testing _any varients

# Tests that the `check` function works correctly for these arguments
# this is a single case whereas `test_check` generates a bunch of cases for different argument types
function test_checker1(
    check::Function, haystack::H, right_needle::N, wrong_needle::N
) where {H,N}
    @testset "$H, $N" begin
        @testset "2 arg form" begin
            @test check(haystack, right_needle)
            @test !check(haystack, wrong_needle)
        end

        @testset "curried form" begin
            @test check(right_needle)(haystack)
            @test !check(wrong_needle)(haystack)

            @test (!check(wrong_needle))(haystack)  # check negated function
        end
    end
end

@testset "symbol_searching.jl" begin
    @testset "basic varients" begin
        test_checker(contains, "abc", "b", "x")
        test_checker(startswith, "abc", "a", "b")
        test_checker(endswith, "abc", "c", "b")
    end

    @testset "_any varients" begin
        @testset "singleton tuple" begin
            test_checker(contains_any, "abc", ("b",), ("x",))
            test_checker(startswith_any, "abc", ("a",), ("b",))
            test_checker(endswith_any, "abc", ("c",), ("b",))
        end

        @testset "right || wrong == right and wrong || wrong ==wrong" begin
            test_checker(contains_any, "abc", ("b", "q"), ("x", "q"))
            test_checker(startswith_any, "abc", ("a", "q"), ("b", "q"))
            test_checker(endswith_any, "abc", ("c", "q"), ("b", "q"))
        end
    end
end


text:
"""
	ternary(cmd0::String="", arg1=[]; kwargs...)

reads (x,y) pairs and generates PostScript code that will plot lines,
polygons, or symbols at those locations on a map.

Full option list at [`ternary`](http://gmt.soest.hawaii.edu/doc/latest/ternary.html)

Parameters
----------

- **A** : **straight_lines** : -- Str --

    By default, geographic line segments are drawn as great circle arcs. To draw them as straight
    lines, use the -A flag.
    [`-A`](http://gmt.soest.hawaii.edu/doc/latest/ternary.html#a)
- $(GMT.opt_J)
- $(GMT.opt_R)
- $(GMT.opt_B)
- $(GMT.opt_C)
- **G** : **fill** : **markerfacecolor** : **MarkerFaceColor** : -- Str --

    Select color or pattern for filling of symbols or polygons. BUT WARN: the alias 'fill' will set the
    color of polygons OR symbols but not the two together. If your plot has polygons and symbols, use
    'fill' for the polygons and 'markerfacecolor' for filling the symbols. Same applyies for W bellow
    [`-G`](http://gmt.soest.hawaii.edu/doc/latest/ternary.html#g)
- **L** : **labels** : -- Str --            Flags = a/b/c

    Set the labels for the three diagram vertices [none].
    [`-L`](http://gmt.soest.hawaii.edu/doc/latest/ternary.html#l)
- **M** : **no_plot** : -- Bool or [] --

    Do no plotting. Instead, convert the input (a,b,c[,*z*]) records to Cartesian (x,y,[,*z*]) records,
    where x, y are normalized coordinates on the triangle (i.e., 0-1 in xand 0-sqrt(3)/2 in y).
    [`-M`](http://gmt.soest.hawaii.edu/doc/latest/ternary.html#m)
- **N** : **no_clip** : -- Bool or [] --

    Do NOT clip symbols that fall outside map border
    [`-N`](http://gmt.soest.hawaii.edu/doc/latest/ternary.html#n)
- $(GMT.opt_P)
- **S** : **symbol** : **marker** : **Marker** : -- Str --

    Plot symbols (including vectors, pie slices, fronts, decorated or quoted lines).
    [`-S`](http://gmt.soest.hawaii.edu/doc/latest/ternary.html#s)
    Alternatively select a sub-set of symbols using the aliases: **marker** or **Marker** and values:

    + **-**, **x_dash**
    + **+**, **plus**
    + **a**, *, **star**
    + **c**, **circle**
    + **d**, **diamond**
    + **g**, **octagon**
    + **h**, **hexagon**
    + **i**, **v**, **inverted_tri**
    + **n**, **pentagon**
    + **p**, **.**, **point**
    + **r**, **rectangle**
    + **s**, **square**
    + **t**, **^**, **triangle**
    + **x**, **cross**
    + **y**, **y_dash**
- **W** : **line_attrib** : **markeredgecolor** : **MarkerEdgeColor** : -- Str --

    Set pen attributes for lines or the outline of symbols
    [`-W`](http://gmt.soest.hawaii.edu/doc/latest/ternary.html#w)
    WARNING: the pen attributes will set the pen of polygons OR symbols but not the two together.
    If your plot has polygons and symbols, use **W** or **line_attribs** for the polygons and
    **markeredgecolor** or **MarkerEdgeColor** for filling the symbols. Similar to S above.
- $(GMT.opt_U)
- $(GMT.opt_V)
- $(GMT.opt_X)
- $(GMT.opt_Y)
- **axis** : **aspect** : -- Str --
    When equal to "equal" makes a square plot.
- $(GMT.opt_a)
- $(GMT.opt_bi)
- $(GMT.opt_di)
- $(GMT.opt_e)
- $(GMT.opt_f)
- $(GMT.opt_g)
- $(GMT.opt_h)
- $(GMT.opt_i)
- $(GMT.opt_p)
- $(GMT.opt_t)
- $(GMT.opt_swap_xy)
"""
function ternary(cmd0::String="", arg1=[]; caller=[], K=false, O=false, first=true, kwargs...)

	arg2 = []		# May be needed if GMTcpt type is sent in via C
	N_args = isempty_(arg1) ? 0 : 1

	((isempty(cmd0) && isempty_(arg1) || occursin(" -", cmd0)) && return monolitic("ternary", cmd0, arg1)	# Monolitic mode

	d = KW(kwargs)
	output, opt_T, fname_ext = fname_out(d)		# OUTPUT may have been an extension only

	opt_J = " -JX12c/10.4c"                     # Equilateral triangle
	for sym in [:axis :aspect]
		if (haskey(d, sym))
			if (d[sym] == "equal")				# Need also a 'tight' option
				opt_J = " -JX12c"
			end
			break
		end
	end
	cmd, opt_B, opt_J, opt_R = parse_BJR(d, "", caller, O, opt_J)
	cmd = parse_UVXY(cmd, d)
	cmd, = parse_a(cmd, d)
	cmd, opt_bi = parse_bi(cmd, d)
	cmd, opt_di = parse_di(cmd, d)
	cmd, = parse_e(cmd, d)
	cmd, = parse_f(cmd, d)
	cmd, = parse_g(cmd, d)
	cmd, = parse_h(cmd, d)
	cmd, opt_i = parse_i(cmd, d)
	cmd, = parse_p(cmd, d)
	cmd, = parse_t(cmd, d)
	cmd, = parse_swap_xy(cmd, d)
	cmd = parse_params(cmd, d)

	cmd, K, O, opt_B = set_KO(cmd, opt_B, first, K, O)		# Set the K O dance

	# If file name sent in, read it and compute a tight -R if this was not provided
	cmd, arg1, = read_data(d, cmd0, cmd, arg1, opt_R, opt_i, opt_bi, opt_di, false)

	cmd, arg1, arg2, = add_opt_cpt(d, cmd, [:C :color :cmap], 'C', N_args, arg1, arg2)

	cmd = add_opt(cmd, 'L', d, [:L :straight_lines])
	cmd = add_opt(cmd, 'M', d, [:M :offset])
	cmd = add_opt(cmd, 'N', d, [:N :error_bars])

	cmd = add_opt(cmd, 'G', d, [:G :fill])
	opt_Gsymb = ""			# Filling color for symbols
	for sym in [:G :markerfacecolor :MarkerFaceColor]
		if (haskey(d, sym))
			opt_Gsymb = " -G" * arg2str(d[sym])
			break
		end
	end

	opt_Wmarker = ""
	for sym in [:markeredgecolor :MarkerEdgeColor]
		if (haskey(d, sym))
			opt_Wmarker = "0.5p," * arg2str(d[sym])		# 0.25p is so thin
			break
		end
	end

	cmd = add_opt(cmd, 'L', d, [:L :labels])
	cmd = add_opt(cmd, 'M', d, [:M :no_plot])
	cmd = add_opt(cmd, 'N', d, [:N :no_clip])

	opt_W = ""
	pen = build_pen(d)						# Either a full pen string or empty ("")
	if (!isempty(pen))
		opt_W = " -W" * pen
	else
		for sym in [:W :line_attrib]
			if (haskey(d, sym))
				if (isa(d[sym], Tuple))		# Like this it can hold the pen, not extended atts
					opt_W = " -W" * parse_pen(d[sym])
				else
					opt_W = " -W" * arg2str(d[sym])
				end
				break
			end
		end
	end

	opt_S = ""
	for sym in [:S :symbol]
		if (haskey(d, sym))
			opt_S = " -S" * arg2str(d[sym])
			break
		end
	end
	if (isempty(opt_S))			# OK, no symbol given via the -S option. So fish in aliases
		marca = ""
		for sym in [:marker :Marker]
			if (haskey(d, sym))
				t = d[sym]
				if (isa(t, Symbol))	t = string(t)	end
				if (t == "-"     || t == "x-dash")   marca = "-"
				elseif (t == "+" || t == "plus")     marca = "+"
				elseif (t == "a" || t == "*" || t == "star")     marca = "a"
				elseif (t == "c" || t == "circle")   marca = "c"
				elseif (t == "d" || t == "diamond")  marca = "d"
				elseif (t == "g" || t == "octagon")  marca = "g"
				elseif (t == "h" || t == "hexagon")  marca = "h"
				elseif (t == "i" || t == "v" || t == "inverted_tri")  marca = "i"
				elseif (t == "n" || t == "pentagon")  marca = "n"
				elseif (t == "p" || t == "." || t == "point")     marca = "p"
				elseif (t == "r" || t == "rectangle") marca = "r"
				elseif (t == "s" || t == "square")    marca = "s"
				elseif (t == "t" || t == "^" || t == "triangle")  marca = "t"
				elseif (t == "x" || t == "cross")     marca = "x"
				elseif (is3D && (t == "u" || t == "cube"))  marca = "u"
				elseif (t == "y" || t == "y-dash")    marca = "y"
				end
				break
			end
		end
		if (!isempty(marca))
			done = false
			for sym in [:markersize :MarkerSize :size]
				if (haskey(d, sym))
					marca = marca * arg2str(d[sym])
					done = true
					break
				end
			end
			if (!done)  marca = marca * "8p"  end			# Default to 8p
		end
		if (!isempty(marca))  opt_S = " -S" * marca  end
	end

	if (!isempty(opt_S))			#
		opt_ML = ""
		for sym in [:markerline :MarkerLine]
			if (haskey(d, sym))
				if (isa(d[sym], Tuple))	# Like this it can hold the pen, not extended atts
					opt_ML = " -W" * parse_pen(d[sym])
				else
					opt_ML = " -W" * arg2str(d[sym])
				end
				if (!isempty(opt_Wmarker))
					opt_Wmarker = ""
					@warn("markerline overrides markeredgecolor")
				end
				break
			end
		end
		if (!isempty(opt_W) && !isempty(opt_ML))
			@warn("You cannot use both markeredgecolor and W or line_attrib keys.")
		end
	end

	if (!isempty(opt_W) && isempty(opt_S)) 			# We have a line/polygon request
		cmd = [finish_PS(d, cmd * opt_W, output, K, O)]
	elseif (isempty(opt_W) && !isempty(opt_S))		# We have a symbol request
		if (!isempty(opt_Wmarker) && isempty(opt_W))
			opt_Gsymb = opt_Gsymb * " -W" * opt_Wmarker	# Piggy back in this option string
		end
		if (!isempty(opt_ML))						# If we have a symbol outline pen
			cmd = cmd * opt_ML
		end
		cmd = [finish_PS(d, cmd * opt_S * opt_Gsymb, output, K, O)]
	elseif (!isempty(opt_W) && !isempty(opt_S))		# We have both line/polygon and a symbol
		# that is not a vector (because Vector width is set by -W)
		if (opt_S[4] == 'v' || opt_S[4] == 'V' || opt_S[4] == '=')
			cmd = [finish_PS(d, cmd * opt_W * opt_S * opt_Gsymb, output, K, O)]
		else
			if (!isempty(opt_Wmarker))
				opt_Wmarker = " -W" * opt_Wmarker	# Set Symbol edge color
			end
			cmd1 = cmd * opt_W
			cmd2 = replace(cmd, opt_B => "") * opt_S * opt_Gsymb * opt_Wmarker	# Don't repeat option -B
			if (!isempty(opt_ML))					# If we have a symbol outline pen
				cmd1 = cmd1 * opt_ML
			end
			cmd = [finish_PS(d, cmd1, output, true, O)
			       finish_PS(d, cmd2, output, K, true)]
		end
	elseif (!isempty(opt_S) && !isempty(opt_ML))		# We have a symbol outline pen
		cmd = [finish_PS(d, cmd * opt_ML * opt_S * opt_Gsymb, output, K, O)]
	else
		cmd = [finish_PS(d, cmd, output, K, O)]
	end

	return finish_PS_module(d, cmd, "", output, fname_ext, opt_T, K, "ternary", arg1, arg2)
end

# ---------------------------------------------------------------------------------------------------
ternary!(cmd0::String="", arg1=[]; caller=[], K=true, O=true,  first=false, kw...) =
	ternary(cmd0, arg1; caller=caller, K=K, O=O,  first=first, kw...)
ternary!(arg1=[]; caller=[], K=true, O=true,  first=false, kw...) =
	ternary("", arg1; caller=caller, K=K, O=O,  first=first, kw...)

psternary  = ternary            # Aliases
psternary! = ternary!           # Aliases


text:
# Implementation of the "bwsj" method, from
# Sheather and Jones, "A Reliable Data-based Bandwidth Selection Method for Kernel Density Estimation"
# Notation below follows roughly the notation in the paper.

struct Block
	depth::Int
	k1::Int
	k2::Int
	lb::Float64
	ub::Float64
end
Base.isless(b1::Block,b2::Block) = isless(b1.ub-b1.lb, b2.ub-b2.lb)

function ϕsum(::Val{N}, ::Type{T}, k1, k2, α, X, leafSize) where {N,T}
	n = length(X)

	i1 = (k1-1)*leafSize+1
	i2 = min(k1*leafSize, n)

	s = zero(T)
	if k1==k2
		@inbounds for i in i1:i2
			for j in i+1:i2 # upper triangular part
				s += ϕ(Val{N}(), (X[j]-X[i])/α)
			end
		end
	else
		j1 = (k2-1)*leafSize+1
		j2 = min(k2*leafSize, n)

		@inbounds for i in i1:i2
			for j in j1:j2
				s += ϕ(Val{N}(), (X[j]-X[i])/α)
			end
		end
	end
	s # exact
end

function bounds(::Val{N}, ::Type{T}, d, k1, k2, α, X, tree::SumTree)::Tuple{T,T} where {N,T}
	n = length(X)
	intervalSize = tree.leafSize * 2^(length(tree.intervalSums)-d)

	i1 = (k1-1)*intervalSize+1
	i2 = min(k1*intervalSize, n)

	if k1==k2 # block on the diagonal
		npoints = div((i2-i1+1)*(i2-i1),2)
		npoints==0 && return zero(T),zero(T) # degenerate case
		x = tree.diagonalSums[d][k1]/(npoints*α)
		npoints.*ϕbounds(Val{N}(), 0.0, (X[i2]-X[i1])/α, x)
	else
		j1 = (k2-1)*intervalSize+1
		j2 = min(k2*intervalSize, n)
		ni = intervalSize
		nj = (j2-j1+1)
		npoints = ni*nj
		meanI = tree.intervalSums[d][k1]/(ni*α)
		meanJ = tree.intervalSums[d][k2]/(nj*α)
		npoints.*ϕbounds(Val{N}(), (X[j1]-X[i2])/α, (X[j2]-X[i1])/α, meanJ-meanI)
	end
end

function ssumstep!(heap::BinaryMaxHeap{Block}, ::Val{N}, ::Type{T},α,X,tree::SumTree)::Tuple{T,T} where {N,T}
	block = pop!(heap)

	# remove existing bounds
	lb = -block.lb
	ub = -block.ub

	if block.depth >= depth(tree) # Fallback to exact sum
		s = ϕsum(Val{N}(), T, block.k1, block.k2, α, X, tree.leafSize)
		lb += s
		ub += s
	else # create child blocks
		d = block.depth+1
		k1 = 2block.k1-1 # first child in first direction
		k2 = 2block.k2-1 # first child in second direction
		nbrBlocks = length(tree.intervalSums[d])

		lb11,ub11 = bounds(Val{N}(), T, d, k1, k2, α, X, tree)
		lb += lb11
		ub += ub11
		lb11!=ub11 && push!(heap, Block(d, k1, k2, lb11, ub11))

		if k1<k2<=nbrBlocks # only care about upper triangular part
			lb21,ub21 = bounds(Val{N}(), T, d, k1+1, k2, α, X, tree)
			lb += lb21
			ub += ub21
			lb21!=ub21 && push!(heap, Block(d, k1+1, k2, lb21, ub21))
		end

		if k2<nbrBlocks
			lb12,ub12 = bounds(Val{N}(), T, d, k1, k2+1, α, X, tree)
			lb += lb12
			ub += ub12
			lb12!=ub12 && push!(heap, Block(d, k1, k2+1, lb12, ub12))

			lb22,ub22 = bounds(Val{N}(), T, d, k1+1, k2+1, α, X, tree)
			lb += lb22
			ub += ub22
			lb22!=ub22 && push!(heap, Block(d, k1+1, k2+1, lb22, ub22))
		end
	end
	lb,ub
end

# Let s = ∑ᵢ∑ⱼϕ⁽ᴺ⁾(α⁻¹(Xᵢ-Xⱼ)), where the sum is taken over i<j.
# returns lower and upper bounds for s, terminating when the tolerance is good enough
function sumapprox(::Val{N},::Type{T},α,X,tree::SumTree;rtol)::Tuple{T,T} where {N,T}
	heap = BinaryMaxHeap{Block}()
	lb,ub = bounds(Val{N}(),T,1,1,1,α,X,tree)
	push!(heap,Block(1,1,1,lb,ub))
	while !isempty(heap) && !isapprox(lb,ub;rtol=rtol)
		lb,ub = (lb,ub) .+ ssumstep!(heap,Val{N}(),T,α,X,tree)
	end
	lb,ub
end

# Let s = ∑ᵢ∑ⱼϕⁱᵛ(α⁻¹(Xᵢ-Xⱼ)), where the sum is taken over i<j.
# returns sign(s-C)
function ssign(::Type{T},α,X,tree::SumTree,C)::Int where T
	heap = BinaryMaxHeap{Block}()
	lb,ub = bounds(Val{4}(),T,1,1,1,α,X,tree)
	push!(heap,Block(1,1,1,lb,ub))
	while !isempty(heap) && lb<=C<=ub
		lb,ub = (lb,ub) .+ ssumstep!(heap,Val{4}(),T,α,X,tree)
	end
	lb > C && return 1
	ub < C && return -1
	@warn "Numerical issues for upper/lower bounds, ub=$ub, lb=$lb, C=$C."
	return 0
end

# Computes the sign of the objective function by gradually approximating using lower/upper bounds
function objectivesign(h, X::AbstractArray{T}, tree::SumTree, α2Constant) where T
	n = length(X)
	α2 = α2Constant*h^(5/7) # 1.357[SD(a)/TD(b)]^(1/7) * h^(5/7)

	C = (n-1)*α2Constant^5*h^(-10/7)/√2
	C -= n*ϕ4(0) # get rid of the diagonal entries
	C/=2 # because of symmetry, we can sum over j>i (upper triangular part), effectively halving the sum

	ssign(promote_type(Float64,T), α2, X, tree, C)
end

function SD(α, X::AbstractArray{T}, tree::SumTree; rtol=0.05) where T
	# α⁻⁵∑ᵢ∑ⱼϕⁱᵛ(α⁻¹(Xᵢ-Xⱼ))
	lb,ub = sumapprox(Val{4}(), promote_type(Float64,T), α, X, tree; rtol=rtol)
	s = (lb+ub) + length(X)*ϕ4(0) # 2*(lb+ub)/2 + diagonal entries
	s/α^5
end

function TD(b, X::AbstractArray{T}, tree::SumTree; rtol=0.05) where T
	# b⁻⁷∑ᵢ∑ⱼϕᵛⁱ(b⁻¹(Xᵢ-Xⱼ))
	lb,ub = sumapprox(Val{6}(), promote_type(Float64,T), b, X, tree; rtol=rtol)
	s = (lb+ub) + length(X)*ϕ6(0) # 2*(lb+ub)/2 + diagonal entries
	-s/b^7
end

"""
	bwsj(X; rtol=0.1, leafsize=10, lower=λ*n^(-1/5)*1e-9, upper=λ*n^(-1/5)*1e9)

Bandwidth selection for kernel density estimation using Gaussian kernels.
Based on the paper Sheather, S. J., & Jones, M. C. (1991). A reliable data‐based bandwidth selection method for kernel density estimation. Journal of the Royal Statistical Society: Series B (Methodological), 53(3), 683-690.

Solves the equation [ R(K) / (nσ⁴ₖS_D(α₂(h))) ]^(1/5) - h = 0 by root finding, bounding the objective such that each iteration can be terminated when the sign of the objective is known.

* `rtol` - the relative tolerance for the bandwidth.
* `leafsize` - Below this size, exact computations are used instead of bounds. Does not impact the accuracy of the result, only the execution time.
* `lower` - Bandwidth lower bound. There is usally no need to change the bounds.
* `upper` - Bandwidth upper bound.
"""
function bwsj(X; rtol=0.1, leafsize=10, lower=nothing, upper=nothing)
	issorted(X) || (X=sort(X))
	(length(X)<=1 || X[1]==X[end]) && return NaN # degenerate cases
	n = length(X)
	tree = SumTree(X,leafsize)

	q25,q75 = quantile(X,(0.25,0.75))
	λ = min(q75-q25, std(X)*1.349) # Givivng the same scale estimate as e.g. the KernSmooth R package.
	a = 0.920λ*n^(-1/7)
	b = 0.912λ*n^(-1/9)
	α2Constant = 1.357*(SD(a,X,tree;rtol=0.5rtol)/TD(b,X,tree;rtol=0.5rtol))^(1/7)

	# Decide if we need to deal with bounds better.
	# But note that execution is very fast for bad guesses for h.

	# The Bandwidth scales roughly with n^(-1/5)
	isnothing(lower) && (lower = λ*n^(-1/5)*1e-9)
	isnothing(upper) && (upper = λ*n^(-1/5)*1e9)

	find_zero(h->objectivesign(h,X,tree,α2Constant), (lower,upper), Bisection(), xrtol=rtol)
end


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

using Markdown
using InteractiveUtils

# ╔═╡ db064048-3ec4-11eb-36a6-2de2e93c3b98
mutable struct State
	step
	spoken
	last
end

# ╔═╡ 21d25d36-3ec5-11eb-123c-6179223a23c9
function speak(state, num)
	state.step += 1
	state.spoken[state.last] = state.step - 1
	state.last = num
end

# ╔═╡ 49925646-3ec5-11eb-20a0-8dc9eb4f160d
function nextnum(state)
	if !(state.last in keys(state.spoken))
		return 0
	else
		return state.step - state.spoken[state.last]
	end
end

# ╔═╡ 308d4de0-3ec5-11eb-0556-813bef55256d
function runit(starting, steps; return_history=true)
	history = []
	state = State(0, Dict(), -1)
	for n in starting
		speak(state, n)
		return_history && push!(history, n)
	end
	for i=length(starting)+1:steps
		n = nextnum(state)
		speak(state, n)
		return_history && push!(history, n)
	end
	if return_history
		history
	else
		state.last
	end
end

# ╔═╡ 223129f8-3ee6-11eb-0711-f5b9c5b005ae
starting = [7,14,0,17,11,1,2]

# ╔═╡ 1b4b07d0-3ee6-11eb-0550-69e3dcfc6914
runit(starting, 2020)[end]

# ╔═╡ f6e372ec-3ee5-11eb-2406-31cda672cc86
runit(starting, 30000000, return_history=false)

# ╔═╡ Cell order:
# ╠═db064048-3ec4-11eb-36a6-2de2e93c3b98
# ╠═21d25d36-3ec5-11eb-123c-6179223a23c9
# ╠═49925646-3ec5-11eb-20a0-8dc9eb4f160d
# ╠═308d4de0-3ec5-11eb-0556-813bef55256d
# ╠═223129f8-3ee6-11eb-0711-f5b9c5b005ae
# ╠═1b4b07d0-3ee6-11eb-0550-69e3dcfc6914
# ╠═f6e372ec-3ee5-11eb-2406-31cda672cc86

max_stars_repo_name: st--/advent-of-code

max_issues_repo_name: st--/advent-of-code

max_forks_repo_name: st--/advent-of-code


text:
# NOTE none of this is exported, and follow upstreaming at JuliaRobotics/DistributedFactorGraphs.jl/issues/727
# copied from original in DFG/src/LightDFG/FactorGraphs/BiMaps.jl

# will be overloaded
import Base: getindex, setindex!, firstindex, lastindex, iterate, keys, isempty

struct _BiDictMap{T <: Integer}
    int_sym::Dict{T,Symbol}
    sym_int::Dict{Symbol,T}
end

_BiDictMap{T}(;sizehint=100) where T<:Integer = begin
    int_sym = Dict{T,Symbol}()
    sizehint!(int_sym, sizehint)
    sym_int = Dict{Symbol,T}()
    sizehint!(sym_int, sizehint)
    _BiDictMap{T}(int_sym, sym_int)
end

_BiDictMap(;sizehint=100) = _BiDictMap{Int}(;sizehint=sizehint)

Base.getindex(b::_BiDictMap, key::Int) = b.int_sym[key]
Base.getindex(b::_BiDictMap, key::Symbol) = b.sym_int[key]

# setindex!(b, value, key) = b[key] = value
function Base.setindex!(b::_BiDictMap, s::Symbol, i::Int)
    haskey(b.sym_int, s) && delete!(b.int_sym, b[s])
    haskey(b.int_sym, i) && delete!(b.sym_int, b[i])

    b.int_sym[i] = s
    b.sym_int[s] = i
end

function Base.setindex!(b::_BiDictMap, i::Int, s::Symbol)
    haskey(b.int_sym, i) && delete!(b.sym_int, b[i])
    haskey(b.sym_int, s) && delete!(b.int_sym, b[s])

    b.int_sym[i] = s
    b.sym_int[s] = i
end

function Base.delete!(b::_BiDictMap, i::Int)
    s = b[i]
    delete!(b.int_sym, i)
    delete!(b.sym_int, s)
    return b
end

Base.haskey(b::_BiDictMap, s::Symbol) = haskey(b.sym_int, s)
Base.haskey(b::_BiDictMap, i::Int) = haskey(b.int_sym, i)

Base.length(b::_BiDictMap) = length(b.int_sym)
#NOTE This will work only with LightGraphs that assumes indices 1:nv(g)
Base.firstindex(v::_BiDictMap) = 1
Base.lastindex(v::_BiDictMap) = length(v.int_sym)
Base.iterate(v::_BiDictMap, i=1) = (length(v.int_sym) < i ? nothing : (v.int_sym[i], i + 1))
Base.keys(v::_BiDictMap) = Base.OneTo(length(v.int_sym))
Base.isempty(v::_BiDictMap) = (length(v.int_sym) == 0)


text:
module Example

using Glimmer, Glimmer.FlexUI
using Images, ImageSegmentation

println("Start [$(splitext(basename(@__FILE__))[1])]")

#---------------------------------------------------------------
# define the application and some basic properties such as title and initial window size
#---------------------------------------------------------------
app = App()
prop!(app, :title, "Glimmer Example - Images")
prop!(app, :winInitWidth, 1200)
prop!(app, :winInitHeight, 800)

#---------------------------------------------------------------
# Define the Variables we will be using in this example
#---------------------------------------------------------------

# image file name
filename = addVariable!(app, Variable(name="filename", type="string",
        value=abspath(joinpath(dirname(pathof(Glimmer)), "Data", "horse.jpg"))))
# size of segments
size = addVariable!(app, Variable(name="size", type="int",value=10))
# contains the result of the segmentation
details = addVariable!(app, Variable(name="details", type="string",value=10))
# source and destination images
src_image = addVariable!(app, Variable(name="src_image", type="image", value=""))
dst_image = addVariable!(app, Variable(name="dst_image", type="image", value=""))

#---------------------------------------------------------------
# Define the UI Controls
#---------------------------------------------------------------
ui = VContainer(
    H1Label("Segmantation Example"),
    HContainer(
        Button(
            text="Load a file",
            fileTypes=".png,.jpg",
            variable="filename",
        ),
        Label(text="Current file [\$(filename)]"),
    ),
    Container(
        direction = "row warp",
        align = "space-between center",
        children = [
            Slider(
                text="Segments Size",
                trailingText="[\$()]",
                min=2,
                max=100,
                value=8,
                variable="size"
            ),
            Label(text="\$(details)"),
        ],
    ),
    H1Label("Image Viewer allowing Pan (left-drag) and Zoom (wheel)"),
    PanZoom(
        style="width: 100%; height: 60vh;  border: 1px solid black;",
        content = Container(
            direction = "row",
            children = [
                Card(
                    title="Source Image",
                    # style=card_style,
                    content=VContainer(
                        Image(
                            source="\$(src_image)",
                        ),
                    ),
                ),
                Card(
                    title="Segmented Image",
                    # style=card_style,
                    content=VContainer(
                        Image(
                            source="\$(dst_image)",
                        ),
                    ),
                ),
            ]
        ),
    ),

    Glimmer.exampleSourceAsCard(@__FILE__),     # add the source code of the example as the last control

)
# set the controls for the application
controls!(app, ui)

#---------------------------------------------------------------
# utility function to load the current image from disk
#---------------------------------------------------------------
function loadImage()
    @info "Loading Image [$(filename[])]"
    global src_img = load(filename[])
end

#---------------------------------------------------------------
# define an event that will trigger when the vriable "filename" is changed
#---------------------------------------------------------------
on(filename) do val
    loadImage()
end

#---------------------------------------------------------------
# the render function will be caulled when any of the UI controls changes value
#---------------------------------------------------------------
src_img = nothing
function processImages()
    global src_img

    segments = felzenszwalb(src_img, size[])
    dst_img = map(i->segment_mean(segments,i), labels_map(segments))
    src_image["png"] = src_img
    dst_image["png"] = dst_img
    details[] = "$segments"
end

function render()
    processImages();
end
renderFunction!(app, render)

#---------------------------------------------------------------
# load the initial image and Run the application
#---------------------------------------------------------------
loadImage()
run(app)

println("End [$(splitext(basename(@__FILE__))[1])]")

end # module