text:
using LoopVectorization, Test
# Tests for Tullio

@testset "Tullio Tests" begin
    A = (1:10) .^ 2; K = rand(10);

    function act!(ℛ::AbstractArray, A, 𝒶𝓍i = axes(A,1), 𝒶𝓍j = axes(ℛ,2))
        @avx for j in 𝒶𝓍j
            for i in 𝒶𝓍i
                ℛ[i, j] = A[i] / j
            end
        end
        ℛ
    end
    function act_noavx!(ℛ::AbstractArray, A, 𝒶𝓍i = axes(A,1), 𝒶𝓍j = axes(ℛ,2))
        for j in 𝒶𝓍j
            for i in 𝒶𝓍i
                ℛ[i, j] = A[i] / j
            end
        end
        ℛ
    end
    @test act!(rand(10,10), A) ≈ act_noavx!(rand(10,10), A)

    D = similar(A, 10, 10) .= 999;

    inds = [2,3,5,2];

    function two!(ℛ::AbstractArray, inds, A, 𝒶𝓍j = axes(ℛ,2), 𝒶𝓍i = axes(inds,1))
        @avx for i = 𝒶𝓍i
            for j = 𝒶𝓍j
                ℛ[inds[i], j] = A[j]
            end
        end
        ℛ
    end
    function two_noavx!(ℛ::AbstractArray, inds, A, 𝒶𝓍j = axes(ℛ,2), 𝒶𝓍i = axes(inds,1))
        for i = 𝒶𝓍i
            for j = 𝒶𝓍j
                ℛ[inds[i], j] = A[j]
            end
        end
        ℛ
    end
    @test two!(copy(D), inds, A) == two!(copy(D), inds, A)

    function three!(ℛ::AbstractArray, A, 𝒶𝓍i = axes(ℛ,1))
        @avx for i = 𝒶𝓍i
            ℛ[i] = A[2i + 1] + A[i]
        end
        ℛ
    end
    function three_noavx!(ℛ::AbstractArray, A, 𝒶𝓍i = axes(ℛ,1))
        for i = 𝒶𝓍i
            ℛ[i] = A[2i + 1] + A[i]
        end
        ℛ
    end
    @test three!(rand(4), A) == three_noavx!(rand(4), A)

    function and(A, 𝒶𝓍i = axes(A,1))
        𝒜𝒸𝒸 = true
        @avx for i = 𝒶𝓍i
            𝒜𝒸𝒸 = 𝒜𝒸𝒸 & (A[i] > 0)
        end
        𝒜𝒸𝒸
    end
    @test and(A)
    A[3] = -1
    @test !and(A)

    function dadada!(EV, Fu, V, ♻ = nothing)
        (ndims)(EV) == 5 || (throw)("expected a 5-array EV")
        (ndims)(Fu) == 2 || (throw)("expected a 2-array Fu")
        (ndims)(V) == 4 || (throw)("expected a 4-array V")

        local 𝒶𝓍a = (axes)(EV, 1)
        local 𝒶𝓍k = (axes)(EV, 2)
        (axes)(V, 1) == (axes)(EV, 2) || (throw)("range of index k must agree")
        local 𝒶𝓍iu = (axes)(Fu, 2)
        (axes)(V, 4) == (axes)(Fu, 2) || (throw)("range of index iu must agree")
        local 𝒶𝓍z = (axes)(EV, 4)
        (axes)(V, 3) == (axes)(EV, 4) || (throw)("range of index z must agree")
        local 𝒶𝓍u = (axes)(EV, 5)
        (axes)(Fu, 1) == (axes)(EV, 5) || (throw)("range of index u must agree")
        local 𝒶𝓍x = (axes)(EV, 3)
        (axes)(V, 2) == (axes)(EV, 3) || (throw)("range of index x must agree")

        ℛ = EV

        @avx for u = 𝒶𝓍u
            for z = 𝒶𝓍z
                for x = 𝒶𝓍x
                    for k = 𝒶𝓍k
                        for a = 𝒶𝓍a
                            𝒜𝒸𝒸 = zero(eltype(EV))  # simpler, same error
                            # 𝒜𝒸𝒸 = if ♻ === nothing
                            #         zero(𝒯)
                            #     else
                            #         ℛ[a, k, x, z, u]
                            #     end
                            for iu = 𝒶𝓍iu
                                𝒜𝒸𝒸 = 𝒜𝒸𝒸 + Fu[u, iu] * V[k, x, z, iu]
                            end
                            ℛ[a, k, x, z, u] = 𝒜𝒸𝒸
                        end
                    end
                end
            end
        end
    end
    function dadada_noavx!(EV, Fu, V, ♻ = nothing)
        (ndims)(EV) == 5 || (throw)("expected a 5-array EV")
        (ndims)(Fu) == 2 || (throw)("expected a 2-array Fu")
        (ndims)(V) == 4 || (throw)("expected a 4-array V")

        local 𝒶𝓍a = (axes)(EV, 1)
        local 𝒶𝓍k = (axes)(EV, 2)
        (axes)(V, 1) == (axes)(EV, 2) || (throw)("range of index k must agree")
        local 𝒶𝓍iu = (axes)(Fu, 2)
        (axes)(V, 4) == (axes)(Fu, 2) || (throw)("range of index iu must agree")
        local 𝒶𝓍z = (axes)(EV, 4)
        (axes)(V, 3) == (axes)(EV, 4) || (throw)("range of index z must agree")
        local 𝒶𝓍u = (axes)(EV, 5)
        (axes)(Fu, 1) == (axes)(EV, 5) || (throw)("range of index u must agree")
        local 𝒶𝓍x = (axes)(EV, 3)
        (axes)(V, 2) == (axes)(EV, 3) || (throw)("range of index x must agree")

        ℛ = EV

        @inbounds @fastmath for u = 𝒶𝓍u
            for z = 𝒶𝓍z
                for x = 𝒶𝓍x
                    for k = 𝒶𝓍k
                        for a = 𝒶𝓍a
                            𝒜𝒸𝒸 = zero(eltype(EV))  # simpler, same error
                            # 𝒜𝒸𝒸 = if ♻ === nothing
                            #         zero(𝒯)
                            #     else
                            #         ℛ[a, k, x, z, u]
                            #     end
                            for iu = 𝒶𝓍iu
                                𝒜𝒸𝒸 = 𝒜𝒸𝒸 + Fu[u, iu] * V[k, x, z, iu]
                            end
                            ℛ[a, k, x, z, u] = 𝒜𝒸𝒸
                        end
                    end
                end
            end
        end
    end
    EV, Fu, V = rand(3,3,3,3,3), rand(3,3), rand(3,3,3,3);
    EV2 = similar(EV);
    dadada!(EV, Fu, V)
    dadada_noavx!(EV2, Fu, V)
    @test EV ≈ EV2
end


text:
using Test, Agents, Random

include("api_tests.jl")
include("forest_fire_defs.jl")
include("space_test.jl")
include("data_collector_test.jl")
include("boltzmann_test.jl")
include("schelling_test.jl")
include("CA_test.jl")

ext: jl


text:
x = [-1, 1, 3, 3, -1]
y = [2, 0, -5, 2, -5]

@test_throws MethodError lineplot()
@test_throws MethodError lineplot(sin, x, y)
@test_throws ArgumentError lineplot(Function[], 0, 3)
@test_throws ArgumentError lineplot(Function[], x)
@test_throws ArgumentError lineplot(Function[])
@test_throws Union{BoundsError,DimensionMismatch} lineplot([1,2], [1,2,3])
@test_throws Union{BoundsError,DimensionMismatch} lineplot([1,2,3], [1,2])
@test_throws Union{BoundsError,DimensionMismatch} lineplot([1,2,3], 1:2)
@test_throws Union{BoundsError,DimensionMismatch} lineplot(1:3, [1,2])
@test_throws Union{BoundsError,DimensionMismatch} lineplot(1:3, 1:2)

@testset "Nan / Inf" begin
    lineplot([0, NaN], [0, 0])
    lineplot([0, 0], [+Inf, 0])
    lineplot([0, 0], [0, -Inf])
end

@testset "numeric array types" begin
    for p in (
            @inferred(lineplot(x, y)),
            @inferred(lineplot(float.(x), y)),
            @inferred(lineplot(x, float.(y))),
        )
        @test p isa Plot
        test_ref("references/lineplot/default.txt", @show_col(p))
    end

    for p in (
            @inferred(lineplot(y)),
            @inferred(lineplot(float.(y))),
        )
        @test p isa Plot
        test_ref("references/lineplot/y_only.txt", @show_col(p))
    end

    p = @inferred lineplot(6:10)
    @test p isa Plot
    test_ref("references/lineplot/range1.txt", @show_col(p))

    p = @inferred lineplot(11:15, 6:10)
    @test p isa Plot
    test_ref("references/lineplot/range2.txt", @show_col(p))
end

@testset "axis scaling and offsets" begin
    p = @inferred lineplot(x .* 1e3  .+ 15, y .* 1e-3 .- 15)
    test_ref("references/lineplot/scale1.txt", @show_col(p))

    p = @inferred lineplot(x .* 1e-3 .+ 15, y .* 1e3  .- 15)
    test_ref("references/lineplot/scale2.txt", @show_col(p))
    tx = [-1.,2, 3, 700000]
    ty = [1.,2, 9, 4000000]
    p = @inferred lineplot(tx, ty)
    test_ref("references/lineplot/scale3.txt", @show_col(p))
    p = @inferred lineplot(tx, ty, width=5, height=5)
    test_ref("references/lineplot/scale3_small.txt", @show_col(p))
end

@testset "dates" begin
    d = collect(Date(1999,12,31):Day(1):Date(2000,1,30))
    v = range(0, stop=3pi, length=31)
    p = @inferred lineplot(d, sin.(v), name="sin", height=5, xlabel="date")
    test_ref("references/lineplot/dates1.txt", @show_col(p))
    @test @inferred(lineplot!(p, d, cos.(v), name = "cos")) === p
    test_ref("references/lineplot/dates2.txt", @show_col(p))
end

@testset "line with intercept and slope" begin
    p = @inferred lineplot(y)
    @test @inferred(lineplot!(p, -3, 1)) === p
    test_ref("references/lineplot/slope1.txt", @show_col(p))
    @test @inferred(lineplot!(p, -4, .5, color = :cyan, name = "foo")) === p
    test_ref("references/lineplot/slope2.txt", @show_col(p))
end

@testset "functions" begin
    @test_throws ArgumentError lineplot(sin, 1:.5:12, color=:blue, ylim=(-1.,1., 2.))
    @test_throws ArgumentError lineplot(sin, 1:.5:12, color=(0,0,255), ylim=[-1.,1., 2.])
    p = @inferred lineplot(sin)
    test_ref("references/lineplot/sin.txt", @show_col(p))
    @test @inferred(lineplot!(p, cos)) === p
    test_ref("references/lineplot/sincos.txt", @show_col(p))
    p = @inferred lineplot([sin, cos])
    test_ref("references/lineplot/sincos.txt", @show_col(p))

    p = @inferred lineplot(sin, -.5, 6)
    test_ref("references/lineplot/sin2.txt", @show_col(p))
    @test @inferred(lineplot!(p, cos)) === p
    test_ref("references/lineplot/sincos2.txt", @show_col(p))
    @test @inferred(lineplot!(p, tan, 2.5, 3.5)) === p
    test_ref("references/lineplot/sincostan2.txt", @show_col(p))

    p = @inferred lineplot([sin, cos], -.5, 3)
    test_ref("references/lineplot/sincos3.txt", @show_col(p))

    tmp = [-.5, .6, 1.4, 2.5]
    p = @inferred lineplot(sin, tmp)
    test_ref("references/lineplot/sin4.txt", @show_col(p))
    @test @inferred(lineplot!(p, cos, tmp)) === p
    test_ref("references/lineplot/sincos4.txt", @show_col(p))
    p = @inferred lineplot([sin, cos], tmp)
    test_ref("references/lineplot/sincos4.txt", @show_col(p))

    @test_throws DimensionMismatch lineplot([sin, cos], -.5, 3, name = ["s", "c", "d"])
    @test_throws DimensionMismatch lineplot([sin, cos], -.5, 3, color = [:red])
    p = @inferred lineplot(
        [sin, cos], -.5, 3, name = ["s", "c"], color = [:red, :yellow],
        title = "Funs", ylabel = "f", xlabel = "num", xlim = (-.5, 2.5), ylim = (-.9, 1.2)
    )
    test_ref("references/lineplot/sincos_parameters.txt", @show_col(p))
    p = @inferred lineplot(
        [sin, cos], -.5, 3, name = ["s", "c"], color = [:red, :yellow],
        title = "Funs", ylabel = "f", xlabel = "num", xlim = [-.5, 2.5], ylim = [-.9, 1.2]
    )
    test_ref("references/lineplot/sincos_parameters.txt", @show_col(p))
end

@testset "keyword arguments" begin
    p = @inferred lineplot(x, y, xlim = (-1.5, 3.5), ylim = (-5.5, 2.5))
    test_ref("references/lineplot/limits.txt", @show_col(p))
    p = @inferred lineplot(x, y, xlim = [-1.5, 3.5], ylim = [-5.5, 2.5])
    test_ref("references/lineplot/limits.txt", @show_col(p))

    p = @inferred lineplot(x, y, grid = false)
    test_ref("references/lineplot/nogrid.txt", @show_col(p))

    p = @inferred lineplot(x, y, color = 20, name = "points1")
    test_ref("references/lineplot/blue.txt", @show_col(p))

    p = @inferred lineplot(x, y, name = "points1", title = "Scatter", xlabel = "x", ylabel = "y")
    @test p isa Plot
    test_ref("references/lineplot/parameters1.txt", @show_col(p))

    @test @inferred(lineplot!(p, [0.5, 1, 1.5], name = "points2")) === p
    test_ref("references/lineplot/parameters2.txt", @show_col(p))

    @test @inferred(lineplot!(p, [-0.5, 0.5, 1.5], [0.5, 1, 1.5], name = "points3")) === p
    test_ref("references/lineplot/parameters3.txt", @show_col(p))
    test_ref("references/lineplot/nocolor.txt", @show_nocol(p))

    p = lineplot(x, y, title = "Scatter", canvas = DotCanvas, width = 10, height = 5)
    @test p isa Plot
    test_ref("references/lineplot/canvassize.txt", @show_col(p))
end

@testset "stairs" begin
    sx = [1, 2, 4, 7, 8]
    sy = [1, 3, 4, 2, 7]

    p = @inferred stairs(sx, sy, style = :pre)
    @test p isa Plot
    test_ref("references/lineplot/stairs_pre.txt", @show_col(p))

    p = @inferred stairs(sx, sy)
    test_ref("references/lineplot/stairs_post.txt", @show_col(p))
    p = @inferred stairs(sx, sy, style = :post)
    test_ref("references/lineplot/stairs_post.txt", @show_col(p))

    p = @inferred stairs(sx, sy, title = "Foo", color = :red, xlabel = "x", name = "1")
    @test @inferred(stairs!(p, sx .- .2, sy .+ 1.5, name = "2")) === p
    test_ref("references/lineplot/stairs_parameters.txt", @show_col(p))
    @test @inferred(stairs!(p, sx, sy, name = "3", style = :pre)) === p
    test_ref("references/lineplot/stairs_parameters2.txt", @show_col(p))
    test_ref("references/lineplot/stairs_parameters2_nocolor.txt", @show_nocol(p))

    # special weird case
    p = stairs(sx, [1, 3, 4, 2, 7000])
    test_ref("references/lineplot/stairs_edgecase.txt", @show_col(p))

    p = stairs(sx, sy, width = 10, padding = 3)
    label!(p, :tl, "Hello")
    label!(p, :t, "how are")
    label!(p, :tr, "you?")
    label!(p, :bl, "Hello")
    label!(p, :b, "how are")
    label!(p, :br, "you?")
    lineplot!(p, 1, .5)
    test_ref("references/lineplot/squeeze_annotations.txt", @show_col(p))

    dx = [Date(2020,8,i) for i = 1:10]
    p = stairs(dx,1:10, xlim=(Date(2020,8,1), Date(2020,8,15)))
    test_ref("references/lineplot/stairs_date.txt", @show_col(p))
end

@testset "scales" begin
    x = y = collect(10:10:1000)
    tmp = tempname()
    for s ∈ (:ln, :log2, :log10)
        fscale = UnicodePlots.FSCALES[s]

        xs = fscale.(x)
        ys = fscale.(y)

        # xscale
        savefig(lineplot(xs, y, xlim=extrema(xs), labels=false), tmp; color=true)
        test_ref(tmp, @show_col(lineplot(x, y, xscale=s, labels=false)))

        # yscale
        savefig(lineplot(x, ys, ylim=extrema(ys), labels=false), tmp; color=true)
        test_ref(tmp, @show_col(lineplot(x, y, yscale=s, labels=false)))

        # xscale and yscale
        savefig(lineplot(xs, ys, xlim=extrema(xs), ylim=extrema(xs), labels=false), tmp; color=true)
        test_ref(tmp, @show_col(lineplot(x, y, xscale=s, yscale=s, labels=false)))

        # scale labels
        test_ref("references/lineplot/$(s)_scale.txt", @show_col(lineplot(x, y, xscale=s, yscale=s)))
    end

    lineplot(1:10, 1:10, xscale=x->log10(x))  # arbitrary scale function
    scl(x) = log(x)
    lineplot(1:10, 1:10, yscale=scl)
end


text:
module UsesPreferences
using Preferences

function set_backend(new_backend::String)
    if !(new_backend in ("OpenCL", "CUDA", "jlFPGA"))
        throw(ArgumentError("Invalid backend: \"$(new_backend)\""))
    end

    # Set it in our runtime values, as well as saving it to disk
    @set_preferences!("backend" => new_backend)
    @info("New backend set; restart your Julia session for this change to take effect!")
end

function clear_backend()
    @set_preferences!("backend" => nothing)
    @info("Backend cleared; restart your Julia session for this change to take effect!")
end

const backend = @load_preference("backend", "OpenCL")

# An example that helps us to prove that things are happening at compile-time
function do_computation()
    @static if backend == "OpenCL"
        return "OpenCL is the best!"
    elseif backend == "CUDA"
        return "CUDA; so fast, so fresh!"
    elseif backend == "jlFPGA"
        return "The Future is Now, jlFPGA online!"
    else
        return nothing
    end
end

# A non-compiletime preference
function set_username(username::String)
    @set_preferences!("username" => username)
end
function get_username()
    return @load_preference("username")
end
function has_username()
    return @has_preference("username")
end
function delete_username()
    @delete_preferences!("username")
end

end # module UsesPreferences


text:
"""
    Fermi.MøllerPlesset

Module for running Møller--Plesset perturbation theory computations.
"""
module MollerPlesset
# Import Fermi basics
using Fermi
using Fermi.Options
using Fermi: AbstractWavefunction
using Fermi.Integrals
using Fermi.Orbitals

function mp_header()
    output(repeat("=",80))
    output("|{:22}{:^34}{:22}|", "", "Møller-Plesset Perturbation Theory", "")
    output("|{:34}{:^9}{:34}|", "", "Module  by","")
    output("|{:25}{:^28}{:25}|", "", "G.J.R Aroeira and M.M. Davis", "")
    output(repeat("=",80))
end

"""
    Fermi.MollerPlesset.AbstractMPWavefunction

Abstract type common to all Møller-Plesset wave functions.

# Struct tree

**AbstractMPWavefunction** <: AbstractWavefunction
"""
abstract type AbstractMPWavefunction <: Fermi.AbstractWavefunction end

# Restricted MP2
include("RMP2/RMP2.jl")

end #module


text:
using SPICE
using Base.Test

this_path = dirname(Base.source_path())
data_path = joinpath(this_path,"data")
!isdir(data_path) && mkdir(data_path)
ngspice(file) = run(joinpath(this_path,file)|>((ignorestatus(setenv(`ngspice`;dir=data_path))|>DevNull).>DevNull))
ngspice("files/01.cir")
data1 = read_ascii_table(joinpath(data_path,"01-a.data"),["vo1","vo2","Iin"])
@test size(data1)[1] == 501


text:
export Kernel, cudims
const DEVICE_H = String(open(read, joinpath(@__DIR__,"device.h")))

mutable struct Kernel
    ptx::String
    funs::Vector{CuFunction}
end

function Kernel(kernel::String)
    kernel = "$DEVICE_H\n$kernel"
    ptx = NVRTC.compile(kernel)
    funs = Array{CuFunction}(undef, nthreads()*ndevices())
    Kernel(ptx, funs)
end

function (k::Kernel)(griddims, blockdims, args...; sharedmem=0, stream=C_NULL)
    id = getdevice() * nthreads() + threadid()
    if !isassigned(k.funs, id)
        k.funs[id] = CuFunction(k.ptx)
    end
    f = k.funs[id]

    argptrs = [cubox(args[i]) for i=1:length(args)]
    @apicall(:cuLaunchKernel, (
        Ptr{Cvoid},           # function
        Cuint,Cuint,Cuint,      # grid dimensions (x, y, z)
        Cuint,Cuint,Cuint,      # block dimensions (x, y, z)
        Cuint,                  # shared memory bytes,
        Ptr{Cvoid},             # stream
        Ptr{Ptr{Cvoid}},         # kernel parameters
        Ptr{Ptr{Cvoid}}),         # extra parameters
        f,
        griddims[1], griddims[2], griddims[3],
        blockdims[1], blockdims[2], blockdims[3],
        sharedmem, stream, pointer(argptrs), C_NULL)
end

function cudims(n::Int)
    bx = 512
    gx = n <= bx ? 1 : ceil(Int, n/bx)
    (gx,1,1), (bx,1,1)
end

cubox(x::Int) = cubox(Cint(x))
cubox(x::Ptr) = x
cubox(x::CuPtr) = x.ptr
cubox(x::Ref) = x
cubox(x) = Ref(x)
cubox(x::CuArray) = cubox(CuDeviceArray(x))
cubox(x::CuSubArray) = cubox(CuDeviceArray(x))


text:
# ------------------------------------------------------------
# What is type stability

using BenchmarkTools

# this function is not type stable
random_data(n) = isodd(n) ? rand(Int, n) : rand(Float64, n)

@code_warntype random_data(3)
#=
julia> @code_warntype random_data(3)
Variables
  #self#::Core.Compiler.Const(random_data, false)
  n::Int64

Body::Union{Array{Float64,1}, Array{Int64,1}}
1 ─ %1 = Main.isodd(n)::Bool
└──      goto #3 if not %1
2 ─ %3 = Main.rand(Main.Int, n)::Array{Int64,1}
└──      return %3
3 ─ %5 = Main.rand(Main.Float64, n)::Array{Float64,1}
└──      return %5
=#

function double_sum_of_random_data(n)
    data = random_data(n)
    total = 0
    for v in data
        total += 2 * v
    end
    return total
end

@btime double_sum_of_random_data(100000);
@btime double_sum_of_random_data(100001);
#=
julia> @btime double_sum_of_random_data(100000);
  347.050 μs (2 allocations: 781.33 KiB)

julia> @btime double_sum_of_random_data(100001);
  179.623 μs (2 allocations: 781.39 KiB)
=#

# See how may red marks with Union types?

@code_warntype double_sum_of_random_data(100000);
#=
julia> @code_warntype double_sum_of_random_data(100000);
Variables
  #self#::Core.Compiler.Const(double_sum_of_random_data, false)
  n::Int64
  data::Union{Array{Float64,1}, Array{Int64,1}}
  total::Union{Float64, Int64}
  @_5::Union{Nothing, Tuple{Union{Float64, Int64},Int64}}
  v::Union{Float64, Int64}

Body::Union{Float64, Int64}
1 ─       (data = Main.random_data(n))
│         (total = 0)
│   %3  = data::Union{Array{Float64,1}, Array{Int64,1}}
│         (@_5 = Base.iterate(%3))
│   %5  = (@_5 === nothing)::Bool
│   %6  = Base.not_int(%5)::Bool
└──       goto #4 if not %6
2 ┄ %8  = @_5::Tuple{Union{Float64, Int64},Int64}::Tuple{Union{Float64, Int64},Int64}
│         (v = Core.getfield(%8, 1))
│   %10 = Core.getfield(%8, 2)::Int64
│   %11 = total::Union{Float64, Int64}
│   %12 = (2 * v)::Union{Float64, Int64}
│         (total = %11 + %12)
│         (@_5 = Base.iterate(%3, %10))
│   %15 = (@_5 === nothing)::Bool
│   %16 = Base.not_int(%15)::Bool
└──       goto #4 if not %16
3 ─       goto #2
4 ┄       return total

=#


text:
# MIT License
#
# Copyright (c) 2018 Martin Biel
#
# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
#
# The above copyright notice and this permission notice shall be included in all
# copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
# SOFTWARE.

# VectorOfDecisions #
# ========================== #
struct VectorDecisionConstraintBridge{T, S <: MOI.AbstractVectorSet} <: MOIB.Constraint.AbstractBridge
    constraint::CI{MOI.VectorOfVariables, S}
    decisions::VectorOfDecisions
end

function MOIB.Constraint.bridge_constraint(::Type{VectorDecisionConstraintBridge{T,S}},
                                           model,
                                           f::VectorOfDecisions,
                                           set::S) where {T, S <: MOI.AbstractVectorSet}
    # Perform the bridge mapping manually
    mapped_variables = map(f.decisions) do decision
        bridged = MOIB.bridged_function(
            model,
            AffineDecisionFunction{T}(SingleDecision(decision)))
        bridged.decision_part.terms[1].variable_index
    end
    # Add the bridged constraint
    constraint = MOI.add_constraint(model,
                                    MOI.VectorOfVariables(mapped_variables),
                                    set)
    # Save the constraint indices and the decision to allow modifications
    return VectorDecisionConstraintBridge{T,S}(constraint, f)
end

function MOI.supports_constraint(::Type{<:VectorDecisionConstraintBridge{T}},
                                 ::Type{VectorOfDecisions},
                                 ::Type{<:MOI.AbstractVectorSet}) where T
    return true
end
function MOIB.added_constrained_variable_types(::Type{<:VectorDecisionConstraintBridge})
    return Tuple{DataType}[]
end
function MOIB.added_constraint_types(::Type{<:VectorDecisionConstraintBridge{T,S}}) where {T,S}
    return [(MOI.VectorOfVariables, S), (MOI.ScalarAffineFunction{T}, MOI.EqualTo{T})]
end
function MOIB.Constraint.concrete_bridge_type(::Type{<:VectorDecisionConstraintBridge{T}},
                              ::Type{VectorOfDecisions},
                              S::Type{<:MOI.AbstractVectorSet}) where T
    return VectorDecisionConstraintBridge{T,S}
end

MOI.get(b::VectorDecisionConstraintBridge{T,S}, ::MOI.NumberOfConstraints{MOI.VectorOfVariables, S}) where {T,S} = 1
MOI.get(b::VectorDecisionConstraintBridge{T}, ::MOI.ListOfConstraintIndices{MOI.VectorOfVariables, S}) where {T,S} = [b.constraint]

function MOI.get(model::MOI.ModelLike, attr::MOI.ConstraintFunction,
                 bridge::VectorDecisionConstraintBridge)
    return bridge.decisions
end

function MOI.get(model::MOI.ModelLike, attr::MOI.ConstraintSet,
                 bridge::VectorDecisionConstraintBridge)
    return MOI.get(model, MOI.ConstraintSet(), bridge.constraint)
end

function MOI.get(model::MOI.ModelLike, attr::MOI.ConstraintPrimal,
                 bridge::VectorDecisionConstraintBridge)
    return MOI.get(model, attr, bridge.constraint)
end

function MOI.get(model::MOI.ModelLike, attr::MOI.ConstraintDual,
                 bridge::VectorDecisionConstraintBridge)
    return MOI.get(model, attr, bridge.constraint)
end

function MOI.get(model::MOI.ModelLike, ::DecisionIndex,
                 bridge::VectorDecisionConstraintBridge{T}) where T
    return bridge.constraint
end

function MOI.delete(model::MOI.ModelLike, bridge::VectorDecisionConstraintBridge)
    MOI.delete(model, bridge.constraint)
    return nothing
end

function MOI.set(model::MOI.ModelLike, ::MOI.ConstraintSet,
                 bridge::VectorDecisionConstraintBridge{T,S}, change::S) where {T,S}
    MOI.set(model, MOI.ConstraintSet(), bridge.constraint, change)
    return nothing
end


text:
# key template data structure
@with_kw struct Template
    data::Dict{Symbol,Vector{_ATTR}} = Dict()
    layout::PlotlyAttribute = attr()
end
_json_lower(t::Template) = Dict(
    :data => _json_lower(t.data),
    :layout => _json_lower(t.layout)
)

==(t1::Template, t2::Template) =  t1.data == t2.data && t1.layout == t2.layout

# to allow Template(data = attr(...))
function Base.convert(T::Type{Dict{Symbol,Vector{_ATTR}}}, x::PlotlyAttribute)
    out = T()
    for (k, v) in pairs(x)
        if !(v isa Vector)
            @error "Could not convert to Template.data. Expected argument at $k to be Vector. Found $v"
        end
        want = _ATTR[]
        for item in v
            push!(want, item)
        end
        out[k] = want
    end
    out
end

function Template(data, layout::Layout)
    Template(data, attr(;layout.fields...))
end
_isempty(::Template) = false

Base.deepcopy(t::Template) = Template(deepcopy(t.data), deepcopy(t.layout))

Base.merge(t1::Template, _t2::Template) = merge!(deepcopy(t1), _t2)

function Base.merge!(out::Template, t2::Template)
    for k in keys(t2.data)
        t1_traces = Cycler(get(out.data, k, [attr()]))
        t2_traces = Cycler(t2.data[k])
        n_out = max(length(t1_traces), length(t2_traces))
        out.data[k] = [merge(t1_traces[ix], t2_traces[ix]) for ix in 1:n_out]
    end

    # now update layout
    merge!(out.layout, t2.layout)
    return out
end

# loading templates from disk
_template_dir() = template_dir = joinpath(artifact"plotly-base-artifacts", "templates")
_available_templates() = [replace(x, ".json" => "") for x in readdir(_template_dir())]
function _load_template(name::String)
    # first try to look it up
    name = Symbol(replace(name, ".json" => ""))
    if name in keys(templates) && !ismissing(templates.templates[name])
        return templates.templates[name]
    end

    # otherwise try to load
    template_dir = _template_dir()
    to_load = joinpath(template_dir, string(name, ".json"))

    # check if this template exists
    exists = isfile(to_load)
    if !exists
        @error "Unknown template $(name). Known values are $(_available_templates())"
    end

    # if it does, load and process it
    raw = _symbol_dict(open(JSON.parse, to_load))
    layout = attr(;_symbol_dict(get(raw, :layout, Dict()))...)
    data = Dict{Symbol,Vector{_ATTR}}()

    # data will be object of arrays. Each array has trace attributes
    for (trace_type, trace_val_array) in get(raw, :data, Dict())
        want = _ATTR[]
        for trace_data in trace_val_array
            push!(want, attr(;_symbol_dict(trace_data)...))
        end
        data[Symbol(trace_type)] = want
    end
    out = Template(data, layout)
    templates[name] = out
    return out
end

# Singleton object for handling templates
mutable struct _TemplatesConfig
    default::String
    templates::Dict{Symbol,_Maybe{Template}}
end

# custom set property to allow setting default with validation
Base.setproperty!(tc::_TemplatesConfig, k::Symbol, x) = setproperty!(tc, Val{k}(), x)
function Base.setproperty!(tc::_TemplatesConfig, ::Val{k}, x) where k
    if hasfield(_TemplatesConfig, k)
        return setfield!(tc, k, x)
    end
    return setindex!(tc, x, k)
end
function Base.setproperty!(tc::_TemplatesConfig, ::Val{:default}, x::String)
    # TODO: validation here
    setfield!(tc, :default, x)
end

# custom getproperty to add available
Base.getproperty(tc::_TemplatesConfig, k::Symbol) = getproperty(tc, Val{k}())
function Base.getproperty(tc::_TemplatesConfig, ::Val{k}) where k
    if hasfield(_TemplatesConfig, k)
        return getfield(tc, k)
    end
    return getindex(tc, k)
end
Base.getproperty(tc::_TemplatesConfig, ::Val{:available}) = collect(keys(tc.templates))
Base.propertynames(tc::_TemplatesConfig) = Symbol[fieldnames(_TemplatesConfig)...,  :available]

# setindex! and getindex, and keys work through the `.templates` dict
function Base.setindex!(tc::_TemplatesConfig, val::Template, k::Union{Symbol,String})
    setindex!(tc.templates, val, Symbol(k))
end
function Base.getindex(tc::_TemplatesConfig, k::Symbol)
    if haskey(tc.templates, k) && !ismissing(tc.templates[k])
        return tc.templates[k]
    end
    return _load_template(String(k))
end
function Base.getindex(tc::_TemplatesConfig, k::String)
    parts = strip.(split(k, "+"))
    reduce(merge, getindex.(Ref(tc), Symbol.(parts)))
end

Base.keys(tc::_TemplatesConfig) = keys(tc.templates)

function Base.show(io::IO, ::MIME"text/plain", tc::_TemplatesConfig)
    msg = """
    Templates configuration
    -----------------------
    Default template: $(tc.default)
    Available templates: $(tc.available)
    """
    println(io, msg)
end

const templates = _TemplatesConfig("plotly", Dict(Symbol(k) => missing for k in _available_templates()))


text:
"""
Return 6 surfaces of B-spline solid
"""
function _bsplinesurfaces(M::BSplineManifold{3})
    a = controlpoints(M)
    P1, P2, P3 = bsplinespaces(M)
    I1, I2, I3 = domain(P1), domain(P2), domain(P3)
    n1, n2, n3 = dim(P1), dim(P2), dim(P3)

    t_⚀ = minimum(I1)
    t_⚁ = minimum(I2)
    t_⚂ = minimum(I3)
    t_⚃ = maximum(I1)
    t_⚄ = maximum(I2)
    t_⚅ = maximum(I3)

    B_⚀ = bsplinebasis.(P1,1:n1,t_⚀)
    B_⚁ = bsplinebasis.(P2,1:n2,t_⚁)
    B_⚂ = bsplinebasis.(P3,1:n3,t_⚂)
    B_⚃ = bsplinebasis.(P1,1:n1,t_⚃)
    B_⚄ = bsplinebasis.(P2,1:n2,t_⚄)
    B_⚅ = bsplinebasis.(P3,1:n3,t_⚅)

    a_⚀ = sum(a[i1,:,:,]*B_⚀[i1] for i1 in 1:n1)
    a_⚁ = sum(a[:,i2,:,]*B_⚁[i2] for i2 in 1:n2)
    a_⚂ = sum(a[:,:,i3,]*B_⚂[i3] for i3 in 1:n3)
    a_⚃ = sum(a[i1,:,:,]*B_⚃[i1] for i1 in 1:n1)
    a_⚄ = sum(a[:,i2,:,]*B_⚄[i2] for i2 in 1:n2)
    a_⚅ = sum(a[:,:,i3,]*B_⚅[i3] for i3 in 1:n3)

    M_⚀ = BSplineManifold(a_⚀,(P2,P3))
    M_⚁ = BSplineManifold(a_⚁,(P1,P3))
    M_⚂ = BSplineManifold(a_⚂,(P1,P2))
    M_⚃ = BSplineManifold(a_⚃,(P2,P3))
    M_⚄ = BSplineManifold(a_⚄,(P1,P3))
    M_⚅ = BSplineManifold(a_⚅,(P1,P2))

    return (M_⚀, M_⚁, M_⚂, M_⚃, M_⚄, M_⚅)
end

function _arrayofvector2array(a::AbstractArray{<:AbstractVector{T},d})::Array{T,d+1} where {d, T<:Real}
    d̂ = length(a[1])
    s = size(a)
    N = prod(s)
    a_2dim = [a[i][j] for i in 1:N, j in 1:d̂]
    a′ = reshape(a_2dim, s..., d̂)
    return a′
end

function _array2arrayofvector(a::AbstractArray{T,d}) where {d, T<:Real}
    s = size(a)
    d̂ = s[end]
    N = s[1:end-1]
    a_flat = reshape(a,prod(N),d̂)
    a_vec = [SVector{d̂}(a_flat[i,:]) for i in 1:prod(N)]
    a′ = reshape(a_vec,N)
    return a′
end


text:
module TestVariableTypes

println("\nTestVariableTypes: Demonstrating the handling of various variable types")

@static if VERSION < v"0.7.0-DEV.2005"
    using Base.Test
    identity(m,n) = eye(m,n)
else
    using Test
    using LinearAlgebra
    identity(m, n) = Matrix{Float64}(I, m, n)
end

using Modia
using ModiaMath: plot

using StaticArrays
using Unitful
using Unitful.DefaultSymbols
import Unitful:
    nm, μm, mm, cm, m, km, inch, ft, mi,
    ac,
    mg, g, kg, A,
    °Ra, °F, °C, K,
    rad, °,
    ms, s, minute, hr,
    J, A, N, mol, cd, V,
    mW, W, Hz

#const Vec3 = Vec{3,Float64}
const Vec3 = SVector{3,Float64}

@testset "VariableTypes" begin

    # Var(; args...) = Variable(; args...)
    # Float(; args...) = Variable(T=Float64; args...)

    FloatScalar(; args...) = Variable(T=Float64, size=(); args...)
    FloatArray(; args...) = Variable(T=Array{Float64,1}; args...)
    FloatMatrix(; args...) = Variable(T=Array{Float64,2}; args...)

    Float3(;args...) = Variable(T=Float64, size=(3,); args...)
    Float3x3(;args...) = Variable(T=Float64, size=(3, 3); args...)

    FixedFloat3(;args...) = Variable(T=Vec{3,Float64}; args...)
    FixedFloat3x3(;args...) = Variable(T=Mat{3,3,Float64}; args...)

    Voltage(; args...) = Variable(T=Unitful.V; args...)

    # Real = (;args...)->Variable(typ=Float64; args...)

    @model TestVariableTypes1 begin
        f = Float()
        b = Boolean()
        i = Integ()
        s = Str()
        c = Var(T=Complex)
        re = Float()
        im = Float()
        c1 = Var(start=Complex(2.0, 0.0))
        # c2 = Var(start=Complex(2.0, 3.0))  # InexactError
        @equations begin
            f = 1
            b = true
            i = 1
            s = "asdf"
            c = Complex(2.0, 3.0)
            re = real(c)
            im = imag(c)
            der(c1) = Complex(1.0, 0.0)
            # der(c2) = Complex(1.5, 2.5)  # InexactError # Allocation of Complex in state vector not supported yet.
        end
    end

    result = checkSimulation(TestVariableTypes1, 3, "i", 1)
    @show result["f"][end]
    @show result["b"][end]
    @show result["i"][end]
    @show result["s"][end]
    @show result["c"][end]
    @show result["re"][end]
    @show result["im"][end]
    # @show result["c2"][end]
    @test result["s"][end] == "asdf"

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

    @model TestArrays1 begin
        f = Float(start=[1,2,3])
        b = Boolean()
        i = Integ()
        s = Str()
        c1 = Var(start=[Complex(1.0, 0.0), Complex(2.0)])

        @equations begin
            der(f) = [2,4,6]
            b = [false,true]
            i = [1,2]
            s = ["asdf", "qwerty"]
            der(c1) = [Complex(2.0, 0.0), Complex(4.0, 0.0)]
            # der(c2) = Complex(1.0, 2.0)  # InexactError
        end
    end

    result = checkSimulation(TestArrays1, 1, "i", [1,2])
    # result = checkSimulation(TestArrays1, 1, "i", [1,2], expandArrayIncidence=true)
    if result != nothing
        @show keys(result)
        @show result["f"][end,:]
        @show result["der(f)"][end,:]
        @show result["b"][end]
        @show result["i"][end]
        @show result["s"][end]
        @show result["c1"][end,:]
        @show result["der(c1)"][end,:]
    end
    @test result["i"][end] == [1,2]
    println()

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

    @model TestVariableTypes2 begin
        scalar = FloatScalar()
        vector = FloatArray()
        matrix = FloatMatrix()
        anyVector = Variable(T=Array{Float64,1})
        v1 = Variable(T=SVector{3,Float64}, start=ones(3), min=3)
        v2 = Variable(T=Vec3, start=ones(3), min=3)
        v3 = Variable(T=Vec3)

        v4 = Float3(min=3)
        v5 = Var()
        v6 = Voltage(start=0.0)
        @equations begin
            scalar = 1
            vector = [1.0, 2.0, 3.0]
            matrix = identity(3, 3)
            anyVector = [1,2,3]
            v1 = ones(3)
            v2 = ones(3)
            v3 = 30.0

            v4 = [1,2,3]
            v5 = [1,2,3]
            v6 = 10.0u"V"
        end
    end

    result = simulate(TestVariableTypes2, 1, storeEliminated=false)

    println("Variable(T=Array{Float64,1}; args...) does not work with storeEliminated=true!")

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

    @model TestVariableUnits1 begin
        # Real3 := Variable{Float64,3}
        # typealias Voltage Variable{Volt}
        v1 = Variable()
        v2 = Var() # Voltage()  # Voltage should work!
        v3 = Var(T=typeof(1.0u"kg"), start=1.0u"kg")
        v4 = Var()
        m = 3.1u"kg" + 2.5kg + 3g
        F = 2.5u"kg*m/s^2" + 2u"N"
        a = Var() # (T=typeof(1.0u"m/s^2"))
        l = 2u"m" + 30cm + 55.7mm  # + 2m not possible due to variable m
        @equations begin
            # SIunits
            v1 = 10.0u"m"
            v2 = 5.5u"V"
            # Unitful
            v3 = 1u"kg" + 2u"g" + 10 * sin(2 * time) * u"hg"
            #  der(v3) = 1.0u"kg/s"  # Differential equations with units are not handled yet
            v4 = 3u"kg" * 0.5u"m/s^2" + time * u"N"
            a = F / m
        end
    end

    result = simulate(TestVariableUnits1, 1)
    #result = checkSimulation(TestVariableUnits1, 1, "v1", 10u"m")

    @show result["v2"][end]
    @show result["v3"][end]
    @show result["v4"][end]
    @show result["m"][end]
    @show result["F"][end]
    @show result["a"][end]
    @test result["v3"][end] == 1.9112974268256817u"kg"
    @test result["v4"][end] == 2.5u"kg*m/s^2" == 2.5N
    @test result["a"][end] == 0.8031411743708727u"m/s^2"  # Unit deduction
    @show result["l"][end]

    plot(result, ("v3", "v4", "a"), figure=1)

end

end


text:
"""
    simple_extr_val_sel(data::ClustData,
                             extreme_value_descr_ar::Array{SimpleExtremeValueDescr,1};
                             rep_mod_method::String="feasibility"
                             )
Selects simple extreme values and returns modified data, extreme values, and the corresponding indices.
"""
function simple_extr_val_sel(data::ClustData,
                             extr_value_descr_ar::Array{SimpleExtremeValueDescr,1};
                             rep_mod_method::String="feasibility"
                             )
  idcs = simple_extr_val_ident(data,extr_value_descr_ar)
  extr_vals = extreme_val_output(data,idcs;rep_mod_method=rep_mod_method)
  # for append method: modify data to be clustered to only contain the values that are not extreme values
  if rep_mod_method=="feasibility"
    data_mod = data
  elseif rep_mod_method=="append"
    data_mod = input_data_modification(data,idcs)
  else
    @error("rep_mod_method - "*rep_mod_method*" - does not exist")
  end
  return data_mod,extr_vals,idcs
end

"""
    simple_extr_val_sel(data::ClustData,
                             extreme_value_descr_ar::Array{SimpleExtremeValueDescr,1};
                             rep_mod_method::String="feasibility"
                             )

Wrapper function for only one simple extreme value.
Selects simple extreme values and returns modified data, extreme values, and the corresponding indices.
"""
function simple_extr_val_sel(data::ClustData,
                             extr_value_descr::SimpleExtremeValueDescr;
                             rep_mod_method::String="feasibility"
                             )
  return simple_extr_val_sel(data,[extr_value_descr];rep_mod_method=rep_mod_method)
end

"""
    simple_extr_val_ident(data::ClustData,extreme_value_descr::Array{SimpleExtremeValueDescr,1})

identifies multiple simple extreme values from the data and returns array of column indices of extreme value within data

data_type: any attribute from the attributes contained within *data*
extremum: "min" or "max"
peak_def: "absolute" or "integral"
"""
function simple_extr_val_ident(data::ClustData,
                               extreme_value_descr_ar::Array{SimpleExtremeValueDescr,1})
  idcs = Array{Int,1}()
  # for each desired extreme value description, finds index of that extreme value within data
  for i=1:length(extreme_value_descr_ar)
    append!(idcs,simple_extr_val_ident(data,extreme_value_descr_ar[i]))
  end
  return idcs
end

"""
    simple_extr_val_ident(data::ClustData,extreme_value_descr::SimpleExtremeValueDescr)

identifies a single simple extreme value from the data and returns column index of extreme value

data_type: any attribute from the attributes contained within *data*
extremum: "min" or "max"
peak_def: "absolute" or "integral"
periods: number of periods combined to analyze
"""
function simple_extr_val_ident(data::ClustData,
                               extreme_value_descr::SimpleExtremeValueDescr)
  return simple_extr_val_ident(data, extreme_value_descr.data_type; extremum=extreme_value_descr.extremum, peak_def=extreme_value_descr.peak_def, periods=extreme_value_descr.periods)
end

"""
    simple_extr_val_ident(data::ClustData,data_type::String;extremum="max",peak_def="absolute")

identifies a single simple extreme value from the data and returns column index of extreme value

data_type: any attribute from the attributes contained within *data*
extremum: "min" or "max"
peak_def: "absolute" or "integral"d
"""
function simple_extr_val_ident(data::ClustData,
                               data_type::String;
                               extremum::String="max",
                               peak_def::String="absolute",
                               periods::Int64=1)
  for name in keys(data.data)
    attr=split(name,"-")[1]
    if name==data_type
      return simple_extr_val_ident(data.data[data_type]; extremum=extremum, peak_def=peak_def, periods=periods)
    elseif attr==data_type
      return simple_extr_val_ident(data.data[name]; extremum=extremum, peak_def=peak_def, periods=periods)
    end
  end
  return throw(@error("the provided data type - "*data_type*" - is not contained in data"))
end

"""
    simple_extr_val_ident(data::Array{Float64};extremum="max",peak_def="absolute")
"""
function simple_extr_val_ident(data::Array{Float64};
                               extremum::String="max",
                               peak_def::String="absolute",
                               periods::Int64=1)

  # set data to be compared

  if peak_def=="absolute" && periods==1
      data_eval = data
  elseif peak_def=="integral"
    # The number of periods is substracted by one as k:k+period
    delta_period=periods-1
    data_eval=zeros(1,(size(data,2)-delta_period))
    for k in 1:(size(data,2)-delta_period)
      data_eval[1,k] = sum(data[:,k:(k+delta_period)])
    end
  else
    @error("peak_def - "*peak_def*" and periods $periods - not defined")
  end
  # find minimum or maximum index. Second argument returns cartesian indices, second argument of that is the column (period) index
  if extremum=="max"
    idx_k = findmax(data_eval)[2][2]
  elseif extremum=="min"
    idx_k = findmin(data_eval)[2][2]
  else
    @error("extremum - "*extremum*" - not defined")
  end
  idx=collect(idx_k:(idx_k+delta_period))
  return idx
end

"""
    input_data_modification(data::ClustData,extr_val_idcs::Array{Int,1})

returns ClustData structs with extreme vals and with remaining input data [data-extreme_vals].
Gives extreme vals the weight that they had in data.
This function is needed for the append method for representation modification
"""
function input_data_modification(data::ClustData,extr_val_idcs::Array{Int,1})
  unique_extr_val_idcs = unique(extr_val_idcs)
  K_dn = data.K- length(unique_extr_val_idcs)
  data_dn=Dict{String,Array}()
  index=setdiff(1:data.K,extr_val_idcs)
  for dt in keys(data.data)
    data_dn[dt] = data.data[dt][:,index] #take all columns but the ones that are extreme vals. If index occurs multiple times, setdiff only treats it as one.
  end
  weights_dn = data.weights[index]
  #QUESTION What with deltas and k_ids?
  deltas_dn= data.deltas[:,index]
  k_ids_dn=deepcopy(data.k_ids)
  k_ids_dn_data=k_ids_dn[findall(data.k_ids.!=0)]
  for k in extr_val_idcs
    k_ids_dn_data[k:end].-=1
    k_ids_dn_data[k]=0
  end
  k_ids_dn[findall(data.k_ids.!=0)]=k_ids_dn_data
  return ClustData(data.region,data.years,K_dn,data.T,data_dn,weights_dn,deltas_dn,k_ids_dn;mean=data.mean,sdv=data.sdv)
end

"""
    input_data_modification(data::ClustData,extr_val_idcs::Int)

wrapper function for a single extreme val.
returns ClustData structs with extreme vals and with remaining input data [data-extreme_vals].
Gives extreme vals the weight that they had in data.
"""
function input_data_modification(data::ClustData,extr_val_idcs::Int)
  return input_data_modification(data,[extr_val_idcs])
end

"""
   extreme_val_output(data::ClustData,
                            extr_val_idcs::Array{Int,1};
                            rep_mod_method="feasibility")

Takes indices as input and returns ClustData struct that contains the extreme vals from within data.
"""
function extreme_val_output(data::ClustData,
                            extr_val_idcs::Array{Int,1};
                            rep_mod_method="feasibility")
  unique_extr_val_idcs = unique(extr_val_idcs)
  K_ed = length(unique_extr_val_idcs)
  data_ed=Dict{String,Array}()
  for dt in keys(data.data)
    data_ed[dt] = data.data[dt][:,unique_extr_val_idcs]
  end
  weights_ed=[]
  if rep_mod_method == "feasibility"
    weights_ed = zeros(length(unique_extr_val_idcs))
  elseif rep_mod_method == "append"
    weights_ed = data.weights[unique_extr_val_idcs]
  else
    @error("rep_mod_method - "*rep_mod_method*" - does not exist")
  end
  deltas_ed=data.deltas[:,unique_extr_val_idcs]
  #QUESTION What with k_ids?
  # if original time series period isn't represented by any extreme period it has value 0
  k_ids_ed=zeros(Int64,size(data.k_ids))
  index_k_ids_data=findall(data.k_ids.!=0)
  k_ids_ed_data=k_ids_ed[index_k_ids_data]
  # each original time series period which is represented recieves the number of it's extreme period in this extreme value output
  k_ids_ed_data[unique_extr_val_idcs]=collect(1:K_ed)
  k_ids_ed[index_k_ids_data]=k_ids_ed_data
  extr_vals = ClustData(data.region,data.years,K_ed,data.T,data_ed,weights_ed,deltas_ed,k_ids_ed;mean=data.mean,sdv=data.sdv)
  return extr_vals
end

"""
   extreme_val_output(data::ClustData,
                            extr_val_idcs::Array{Int,1};
                            rep_mod_method="feasibility")

wrapper function for a single extreme val.
Takes indices as input and returns ClustData struct that contains the extreme vals from within data.
"""
function extreme_val_output(data::ClustData,
                            extr_val_idcs::Int;
                            rep_mod_method="feasibility")
  return extreme_val_output(data,[extr_val_idcs];rep_mod_method=rep_mod_method)
end

"""
    representation_modification(extr_vals::ClustData,
                                     clust_data::ClustData,
                                     )

Merges the clustered data and extreme vals into one ClustData struct. Weights are chosen according to the rep_mod_method
"""
function representation_modification(extr_vals::ClustData,
                                     clust_data::ClustData,
                                     )
  K_mod = clust_data.K + extr_vals.K
  data_mod=Dict{String,Array}()
  for dt in keys(clust_data.data)
    data_mod[dt] = [clust_data.data[dt] extr_vals.data[dt]]
  end
  weights_mod = [clust_data.weights; extr_vals.weights]
  deltas_mod = [clust_data.deltas extr_vals.deltas]
  # Question what with k_ids?
  # originial time series periods are regularly represented by periods in clust_data
  k_ids_mod=deepcopy(clust_data.k_ids)
  # if this particular original time series period is though represented in the extreme values, the new period number of the extreme value (clust_data.K+old number) is assigned to this original time series period QUESTION also true for feasibility? As Cost will be assumed zero, maybe there should be an if/else of representation_modification?
  k_ids_mod[findall(extr_vals.k_ids.!=0)]=extr_vals.k_ids[findall(extr_vals.k_ids.!=0)].+clust_data.K
  return ClustData(clust_data.region,clust_data.years,K_mod,clust_data.T,data_mod,weights_mod,deltas_mod,k_ids_mod;mean=clust_data.mean,sdv=clust_data.sdv)
end

"""
function representation_modification(extr_vals::ClustData,
                                     clust_data::ClustData,
                                     )

Merges the clustered data and extreme vals into one ClustData struct. Weights are chosen according to the rep_mod_method
"""
function representation_modification(extr_vals_array::Array{ClustData,1},
                                     clust_data::ClustData,
                                     )
    for extr_vals in extr_vals_array
      clust_data=representation_modification(extr_vals,clust_data)
    end
    return clust_data
end

"""
    representation_modification(full_data::ClustData,
                                     clust_data::ClustData,
                                     extr_val_idcs::Array{Int,1};
                                     rep_mod_method::String="feasibility")

Merges the clustered data and extreme vals into one ClustData struct. Weights are chosen according to the rep_mod_method
"""
function representation_modification(full_data::ClustData,
                                     clust_data::ClustData,
                                     extr_val_idcs::Array{Int,1};
                                     rep_mod_method::String="feasibility")
  extr_vals = extreme_val_output(full_data,extr_val_idcs;rep_mod_method=rep_mod_method)
  return representation_modification(extr_vals,clust_data;rep_mod_method=rep_mod_method)
end

"""
    representation_modification(full_data::ClustData,
                                     clust_data::ClustData,
                                     extr_val_idcs::Int;
                                     rep_mod_method::String="feasibility")

wrapper function for a single extreme val.
Merges the clustered data and extreme vals into one ClustData struct. Weights are chosen according to the rep_mod_method
"""
function representation_modification(full_data::ClustData,
                                     clust_data::ClustData,
                                     extr_val_idcs::Int;
                                     rep_mod_method::String="feasibility")
  return representation_modification(full_data,clust_data,[extr_val_idcs];rep_mod_method=rep_mod_method)
end


text:
"""
    partition(r, P)

Partition `r` into `P` pieces.

This will provide an equal partition when `N = length(r)` is divisible by `P`
and otherwise the ranges will have lengths of either `div(N, P)` or `cld(N, P)`.

# Examples

We can partition `10:20` into three pieces with
```jldoctest
julia> partition(10:20, 3)
3-element Vector{UnitRange{Int64}}:
 10:12
 13:16
 17:20
```

We can partition `10:2:20` into four pieces with
```jldoctest
julia> partition(10:2:20, 4)
4-element Vector{StepRange{Int64, Int64}}:
 10:2:10
 12:2:14
 16:2:16
 18:2:20
```
"""
partition(r, P) = [partition(r, p, P) for p in 1:P]

"""
    partition(r, p, P)

Partition `r` into `P` pieces and return the `p`th piece.

This will provide an equal partition when `N = length(r)` is divisible by `P`
and otherwise the ranges will have lengths of either `div(N, P)` or `cld(N, P)`.

# Examples

We can get the second piece of `10:20` partitioned into three pieces with
```jldoctest
julia> partition(10:20, 2, 3)
13:16
```
We can get the third piece of `10:2:20` partitioned into four pieces with
```jldoctest
julia> partition(10:2:20, 3, 4)
16:2:16
```
"""
partition(r, p, P) = r[(div((p - 1)*length(r), P) + 1):(div(p*length(r), P))]

"""
    hilbertcode(Y::AbstractArray{T}; bits=8sizeof(T)) where T

Given an array of axes coordinates `Y` stored as `bits`-bit integers
the function returns the associated Hilbert integer `H`.

The encoding of the Hilbert integer is best described by example.
If 5-bits are used from each of 3 coordinates then the function performs

     X[2]|                       H[0] = A B C D E
         | /X[1]       ------->  H[1] = F G H I J
    axes |/                      H[2] = K L M N O
         0------ X[0]                   high low

where the 15-bit Hilbert integer = `A B C D E F G H I J K L M N O` is stored
in `H`.

This function is based on public domain code from John Skilling which can be
found in <https://dx.doi.org/10.1063/1.1751381>.

# Examples
We can generate the two-dimensional Hilbert code for a 1-bit integers with
```jldoctest
julia> hilbertcode([0,0], bits=1)
2-element Vector{Int64}:
 0
 0

julia> hilbertcode([0,1], bits=1)
2-element Vector{Int64}:
 0
 1

julia> hilbertcode([1,1], bits=1)
2-element Vector{Int64}:
 1
 0

julia> hilbertcode([1,0], bits=1)
2-element Vector{Int64}:
 1
 1
```

# References
  John Skilling, "Programming the Hilbert curve",
  AIP Conference Proceedings 707, 381 (2004).
  <https://doi.org/10.1063/1.1751381>
"""
function hilbertcode(Y::AbstractArray{T}; bits=8sizeof(T)) where T
  # Below is Skilling's AxestoTranspose
  X = similar(Y)
  X .= Y
  N = length(X)
  M = one(T) << (bits-1)

  Q = M
  for j = 1:bits-1
    P = Q - one(T)
    for i = 1:N
      if X[i] & Q != zero(T)
        X[1] ⊻= P
      else
        t = (X[1] ⊻ X[i]) & P
        X[1] ⊻= t
        X[i] ⊻= t
      end
    end
    Q >>>= one(T)
  end

  for i = 2:N
    X[i] ⊻= X[i - 1]
  end

  t = zero(T)
  Q = M
  for j = 1:bits-1
    if X[N] & Q != zero(T)
      t ⊻= Q - one(T)
    end
    Q >>>= one(T)
  end

  for i = 1:N
    X[i] ⊻= t
  end

  # Below we transpose X and store it in H, i.e.:
  #
  #   X[0] = A D G J M               H[0] = A B C D E
  #   X[1] = B E H K N   <------->   H[1] = F G H I J
  #   X[2] = C F I L O               H[2] = K L M N O
  #
  # The 15-bit Hilbert integer is then = A B C D E F G H I J K L M N O
  H = zero(X)
  for i = 0:N-1, j = 0:bits-1
    k = i * bits + j
    bit = (X[N - mod(k, N)] >>> div(k, N)) & one(T)
    H[N - i] |= (bit << j)
  end

  return H
end

hilbertperm(indices) = hilbertperm(CartesianIndices(indices))
function hilbertperm(indices::CartesianIndices)
    return sortperm(vec(hilbertcode.(map(v->v.-1, SVector.(Tuple.(indices))))))
end

function hilbertindices(indices)
    CIs = CartesianIndices(indices)
    return reshape(invperm(hilbertperm(CIs)), size(CIs))
end

"""
    quantize([T=UInt64], x)

Quantize a number `x`, between `zero(x)` and `one(x)`, to an integer of type `T`
between `zero(T)` and `typemax(T)`.  If `x` is an array or tuple each element is
quantized.

# Examples
```jldoctest
julia> quantize(0.0)
0x0000000000000000

julia> quantize(UInt32, 0.5f0)
0x7fffffff

julia> quantize((1.0, 0.5))
(0xffffffffffffffff, 0x7fffffffffffffff)

```

"""
function quantize end

quantize(x) = quantize(UInt64, x)
quantize(::Type{T}, x::AbstractArray) where T = quantize.(T, x)
quantize(::Type{T}, x::Tuple) where T = quantize.(T, x)
quantize(::Type{T}, x::Number) where {T <: Integer} =
    convert(T, floor(typemax(T) * BigFloat(x; precision = 16sizeof(T))))


text:
using Utils

function convert_seat(s)
    replacements = [
        ('F', 'L'),
        ('B', 'R'),
    ]
    replaced = foldl((acc, (i, r)) -> replace(acc, r => i - 1), enumerate(replacements), init=s)
    parse(Int, replaced, base=2)
end

function load()
    re = []
    open(joinpath(@__DIR__, "05_plane-seats.txt"), "r") do f
        for line in eachline(f)
            push!(re, convert_seat(line))
        end
    end
    re
end

function find_missing_on_sorted(seats)
    for i in 2:length(seats)
        if seats[i-1] + 2 == seats[i]
            return seats[i-1] + 1
        end
    end
end
find_missing = find_missing_on_sorted ∘ sort

seats = load()
println("Highest seat number: $(maximum(seats))")
println("My seat is $(find_missing(seats))")

using Test

@testset "Binary Boarding" begin
    @testset "example" begin
        @test convert_seat("FBFBBFFRLR") == 357
        @test convert_seat("BFFFBBFRRR") == 567
        @test convert_seat("FFFBBBFRRR") == 119
        @test convert_seat("BBFFBBFRLL") == 820
    end

    @testset "results" begin
        @test maximum(seats) == 919
        @test find_missing(seats) == 642
    end
end

max_stars_repo_name: CmdQ/AoC


text:
using Revise

using BridgeLandmarks
const BL = BridgeLandmarks
using RCall
using Random
using Distributions
using DataFrames
using DelimitedFiles
using CSV
using StaticArrays
using LinearAlgebra
using JLD2
using FileIO
using Parameters
#using Debugger

Random.seed!(9)
workdir = @__DIR__
cd(workdir)
outdir = joinpath(workdir,"out")

## available update steps implemented:
println("available steps are [:mala_mom, :rmmala_mom, :rmrw_mom, :mala_pos, :rmmala_pos, :sgd_mom :parameter, :innov]")

####################### test landmarksmatching ######################
n = 6# 18
xobs0 = [PointF(2.0cos(t), sin(t))/4.0  for t in collect(0:(2pi/n):2pi)[2:end]]
θ, ψ =  π/3, 0.1
rot =  SMatrix{2,2}(cos(θ), sin(θ), -sin(θ), cos(θ))
stretch = SMatrix{2,2}(1.0 + ψ, 0.0, 0.0, 1.0 - ψ)
shift = [0.5, -1.0]
xobsT = [rot * stretch * xobs0[i] + shift  for i in 1:n]

@time landmarksmatching(xobs0,xobsT; ITER=300,pars=Pars_ms(δmom=0.01), outdir=outdir)
plotlandmarksmatching(outdir)

@time landmarksmatching(xobs0,xobsT; ITER=300,pars=Pars_ms(δmom=0.01), outdir=outdir,updatescheme=[:innov, :mala_mom, :parameter])

@time landmarksmatching(xobs0,xobsT; ITER=100,pars=Pars_ahs(δmom=0.01, db=[1.0,2.0],stdev=1.0), outdir=outdir)

@time landmarksmatching(xobs0,xobsT; ITER=300,pars=Pars_ms(δmom=0.001), outdir=outdir,updatescheme=[:innov, :sgd_mom, :parameter])

@time landmarksmatching(xobs0,xobsT; updatescheme=[:innov, :rmrw_mom], ITER=500,pars=Pars_ms(δmom=0.01), outdir=outdir)

landmarksmatching(xobs0,xobsT; outdir=outdir,ITER=500,
                            updatescheme = [:innov, :parameter,:mala_mom, :matching],
                            pars= Pars_ms(covθprop = Diagonal(fill(0.001,3)),η = n -> min(0.2, 10/(n))  ))

plotlandmarksmatching(outdir)
###################### test template estimation ######################
n = 6
nshapes = 8
q0 = [PointF(2.0cos(t), sin(t))  for t in collect(0:(2pi/n):2pi)[2:end]] # initial shape is an ellipse
x0 = State(q0, randn(PointF,n))
xobsT = Vector{PointF}[]
σobs = 0.05

T = 1.0; dt = 0.01; t = 0.0:dt:T
Ptrue = MarslandShardlow(2.0, 0.2, 1.7, 0.0, n)
for k in 1:nshapes
    Wf, Xf = BridgeLandmarks.landmarksforward(t, x0, Ptrue)
    push!(xobsT, Xf.yy[end].q + σobs * randn(PointF,n))
end

template_estimation(xobsT; xinitq=2.0*xobsT[1],outdir=outdir, ITER=100, updatescheme = [:innov, :rmmala_pos, :parameter])  # deliberately initialise badly to show it works
plottemplate_estimation(outdir)

# cyclically shift landmarks for a few landmarks
xobsTshift = copy(xobsT)
xobsTshift[1] = Base.circshift(xobsT[1],2)
xobsTshift[2] = Base.circshift(xobsT[2],-1)
xobsTshift[3] = Base.circshift(xobsT[3],-3)
template_estimation(xobsTshift; xinitq=2.0*xobsT[3],outdir=outdir, updatescheme = [:innov, :rmmala_pos, :parameter, :matching],
    ITER=200, pars=Pars_ms(δpos=0.001,skip_saveITER=5))  # deliberately initialise badly to show it works

####################### below code for generating input data as matrices for landmarksmatching ######################
if false
    dat = load("../experiments/exp1/data_exp1.jld2")
    xobs0 = dat["xobs0"]
    xobsT = dat["xobsT"]
    if false # make example dataset
        writedlm("landmarks0.txt", hcat(extractcomp(xobs0,1), extractcomp(xobs0,2)))
        writedlm("landmarksT.txt", hcat(extractcomp(xobsT,1), extractcomp(xobsT,2)))
    end
    landmarks0 = readdlm("landmarks0.txt")
    landmarksT = readdlm("landmarksT.txt")

    landmarksmatching(xobs0,xobsT)
    landmarksmatching(landmarks0, landmarksT; outdir=outdir, ITER=10)
    landmarksmatching(landmarks0, landmarksT; outdir=outdir, ITER=10, pars=BL.Pars_ahs())
end

"""
een functie ``\\sin(x) + \\frac{x}{\\theta}``
"""
function frank()
    println("testje")
end


text:
using MetalheadUtils
using Test
using Images, Metalhead

image_urls = ["https://upload.wikimedia.org/wikipedia/commons/d/d7/Green_Sea_Turtle_grazing_seagrass.jpg", "https://upload.wikimedia.org/wikipedia/commons/0/0a/The_Great_Wave_off_Kanagawa.jpg",
"https://upload.wikimedia.org/wikipedia/commons/b/b4/Vassily_Kandinsky%2C_1913_-_Composition_7.jpg",
"https://upload.wikimedia.org/wikipedia/commons/0/00/Tuebingen_Neckarfront.jpg",
"https://upload.wikimedia.org/wikipedia/commons/6/68/Pillars_of_creation_2014_HST_WFC3-UVIS_full-res_denoised.jpg",
"https://upload.wikimedia.org/wikipedia/commons/thumb/e/ea/Van_Gogh_-_Starry_Night_-_Google_Art_Project.jpg/1024px-Van_Gogh_-_Starry_Night_-_Google_Art_Project.jpg"];

@testset "MetalheadUtils.jl" begin
    images = imgFromURL.(image_urls)
    @test length(images) == length(image_urls)
    @test images[1] isa Array{RGB{FixedPointNumbers.Normed{UInt8,8}},2}

    img1 = postprocess(Metalhead.preprocess(images[1]))
    @test img1 isa Array{RGB{FixedPointNumbers.Normed{UInt8,8}},2}
    r = img1 - impreprocess(images[1])
    @test maximum(abs.(channelview(r))) <= 1/500
end

using Flux
@testset "Features" begin
    r = Chain(x->x+1, x->x+2)
    @test Features(r, 20, [1]) == [21]
    @test Features(r, 20, [1,2]) == [21,23]

    vgg = VGG19()
    img1 = Metalhead.preprocess(imgFromURL(image_urls[1]))
    @test Features(vgg, img1, [1,2]) == [Features(vgg, img1, [1])[1], Features(vgg, img1, [2])[1]]
end


text:
module VPTrees

include("vantage_point_tree.jl")

export VPTree, find, find_nearest

end # module


text:
using Makie, VideoIO, Dates, Observables, DataStructures

include("testvideo.jl")

testvideo = joinpath(tempdir(), "testvideo.mp4")

# create the test video if it doesn't exist
if !isfile(testvideo)
    createtestvideo(testvideo)
end

avf = VideoIO.testvideo("annie_oakley")
testvideo = avf.io
# f = VideoIO.openvideo(avf)

# open the video
f = VideoIO.openvideo(testvideo)

duration = VideoIO.get_duration(f.avin.io)

img = read(f)
T = typeof(img)
bufferf = VideoIO.openvideo(testvideo)

# bufferradius = f.width*f.height*VideoIO.get_duration(f.avin.io)*f.framerate
bufferradius = 200
buffertime = bufferradius/f.framerate
bufferlength = 2bufferradius + 1
buffer = CircularBuffer{Pair{Float64, T}}(bufferlength)

function fillbuffer()
    for _ in 1:bufferlength
        push!(buffer, Pair(gettime(bufferf), read(bufferf)))
    end
end
fillbuffer()

current = Observable(first(buffer[bufferradius + 1]))

correctcurrent = lift(current) do c
    if c ∈ first.(buffer)
        c
    else
        seek(f, c)
        gettime(f)
    end
end

timestamp = lift(correctcurrent) do cc
    Time(0) + Millisecond(round(Int, 1000cc))
end

img = lift(correctcurrent) do cc
    i = findfirst(timg -> first(timg) > cc, buffer)
    if isnothing(i)
        read(f)
    else
        last(buffer[i])
    end
end
updatebuffer = on(correctcurrent) do cc
    m = max(cc - buffertime, 0.0)
    M = min(cc + buffertime, duration)
    mi = findfirst(timg -> first(timg) ≥ m, buffer)
    Mi = findfirst(timg -> first(timg) ≥ M, buffer)
    if mi == 1 && m > 0
        seek(bufferf, m)
        if Mi == 1
            fillbuffer()
        else
            # @assert Mi < bufferlength
            news = Vector{Pair{Float64, T}}()
            push!(news, Pair(gettime(bufferf), read(bufferf)))
            while first(news[end]) < first(first(buffer))
                push!(news, Pair(gettime(bufferf), read(bufferf)))
            end
            pop!(news)
            for timg in Iterators.reverse(news)
                pushfirst!(buffer, timg)
            end
        end
    elseif !isnothing(mi)
        # @assert isnothing(Mi) m, M, mi, Mi
        while first(buffer[end]) < M
            push!(buffer, Pair(gettime(bufferf), read(bufferf)))
        end
    else
        # @assert isnothing(Mi) m, M, mi, Mi
        seek(bufferf, m)
        fillbuffer()
    end
    0
end

# all the stuff needed for display
pixelaspectratio = VideoIO.aspect_ratio(f)
h = f.height
w = round(typeof(h), f.width * pixelaspectratio)
AbstractPlotting.set_theme!(backgroundcolor = :black)
scene = Scene(resolution = (w, h))
Makie.image!(scene, 1:h, 1:w, img, show_axis = false, scale_plot = false)
Makie.rotate!(scene, -0.5π)
update_limits!(scene)
update_cam!(scene)
Makie.update!(scene)

# the slider, updating currentframe and getting its string from timestamps
slider_h = slider(range(0, duration, step = 1/f.framerate), start = 0, valueprinter = i -> "no, this will be wrong", textcolor = :white)

lift(slider_h[end][:value]) do t
    current[] = t
end

# this is for a separate display of the time stamp, due to the fact that we can not calculate with the time stamp is from the frame number (which is what the slider has)
timestamp_h = text(lift(string, timestamp), color = :white, raw = true, camera = campixel!)

# I have to have a forward button becuase apparently the stuff in the lift gets evaluated, so to have the movie start in frame #1 I need to go forward once and the backwards...
fwdbutton = button(">", textcolor = :white)
lift(fwdbutton[end][:clicks]) do _
    i = findfirst(timg -> first(timg) == correctcurrent[], buffer)
    @assert !isnothing(i)
    correctcurrent[] = first(buffer[i])
    slider_h[end][:value][] += 1
    # so now the movie is in frame #2
    # end
end
#=bckbutton = button!(rsc(), "<", textcolor = :white)
lift(bckbutton[end][:clicks]) do _
if currentframe[] > 1 # check we're not in frame #1
currentframe[] -= 1
slider_h[end][:value][] -= 1
# and now we're back in frame #1
end
0
end=#

# done!
hbox(vbox(slider_h, fwdbutton, timestamp_h), scene)

# played = Observables.async_latest(play, play_bubtton[end][:clicks])

# hbox(vbox(play_button, next_button, slider_h), scene)


text:
using MAT
deltas = 10 .^LinRange(0,3.7,10)
seednum = 100
epsilon = 1.0
grownum = 10000

# M = matread("../data/processed/stackoverflow_answer_H.mat")
LabelMatrix = M["LabelMatrix"]
LabelNames = M["LabelNames"]
H = M["H"]
order = round.(Int64,vec(sum(H,dims=2)))
d = vec(sum(H,dims=1))
volA = sum(d)
m,n = size(H)

labels = [25849; 27596;28918;29386;43507]
lnum = length(labels)
v = [22943;28886;5713]

##
using Plots
plot()
s1 = 300
s2 = 250
ms = 4
lw = 2
using LaTeXStrings
x_label = L"\delta"
y_label = "F1 Scores"
for i = 1:3
    label = labels[i]
    outputmat = "Output_VaryDelta/Cluster_$(label)_$grownum"*"_$seednum"*"_$epsilon.mat"
    mat = matread(outputmat)
    hl_f1 = mat["hl_f1"]
    deltas = mat["deltas"]

    plot!(deltas, hl_f1[i,:],grid = false,label = LabelNames[label],
     markerstrokewidth = 0, markershape = :circle, linewidth = 2, xaxis = :log10,
     size = (s1,s2), markersize = ms,
     xlabel = x_label, ylabel = y_label, legend = :bottomright)
end

## These were run and stored and slightly differently
v = [22943;28886;5713]
for i = 1:3
    label = v[i]
    outputmat = "Output_VaryDelta/Cluster_$(label)_$grownum"*"_$seednum"*"_$epsilon.mat"
    mat = matread(outputmat)
    hl_f1 = mat["hl_f1"]

    plot!(deltas, vec(hl_f1[:,1]),grid = false,label = LabelNames[label],
     markerstrokewidth = 0, markershape = :circle, linewidth = 2, xaxis = :log10,
     size = (s1,s2), markersize = ms,xguidefont=font(18),legend = false,
     xlabel = x_label, ylabel = y_label, legendfont=font(7))

end

savefig("Plots/Stack_VaryDelta.pdf")

max_stars_repo_name: nveldt/HypergraphFlowClustering

max_issues_repo_name: nveldt/HypergraphFlowClustering

max_forks_repo_name: nveldt/HypergraphFlowClustering


text:
module TestContainers
	using Base.Test
	using DataArrays

	dv = @data ones(3)
	push!(dv, 3.0)
	push!(dv, NA)

	@assert isequal(dv, (@data [1.0, 1.0, 1.0, 3.0, NA]))

	a, b = pop!(dv), pop!(dv)
	@assert isna(a)
	@assert b == 3.0

	unshift!(dv, 3.0)
	unshift!(dv, NA)

	@assert isequal(dv, (@data [NA, 3.0, 1.0, 1.0, 1.0]))

	a, b = shift!(dv), shift!(dv)
	@assert isna(a)
	@assert b == 3.0
end


text:
# all other glasses should follow the format below, new glasses must be added to OTHER_GLASSES and OTHER_GLASS_NAMES where the index in the array matches the numeric part of the GlassID

module CARGILLE
using ..GlassCat: Glass, GlassID, OTHER
using StaticArrays: SVector

"""    CARGILLE.OG0608
```
ID:                      OTHER:1
RI @ 587nm:              1.457518
Abbe Number:             57.18978
ΔPgF:                    0.008
TCE (÷1e-6):             800.0
Density:                 0.878g/m³
Valid wavelengths:       0.32μm to 1.55μm
Reference Temp:          25°C
```
"""
const OG0608 = Glass(GlassID(OTHER, 1), -2, 1.4451400, 0.0043176, -1.80659e-5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.32, 1.55, -0.0009083144750540808, 0.0, 0.0, 0.0, 0.0, 0.0, 25.0, 0.008, -1.0, -1.0, 800.0, -1.0, 0, -1.0, [SVector(0.32, 0.03, 10.0), SVector(0.365, 0.16, 100.0), SVector(0.4047, 0.40, 100.0), SVector(0.480, 0.71, 100.0), SVector(0.4861, 0.72, 100.0), SVector(0.5461, 0.80, 100.0), SVector(0.5893, 0.90, 100.0), SVector(0.6328, 0.92, 100.0), SVector(0.6439, 0.95, 100.0), SVector(0.6563, 0.96, 100.0), SVector(0.6943, 0.99, 100.0), SVector(0.840, 0.99, 100.0), SVector(0.10648, 0.74, 100.0), SVector(0.1300, 0.39, 100.0), SVector(0.1550, 0.16, 100.0)], 1.457518, -1.0, -1.0, 0, 57.18978, 0, 0.878, -1)

"""    CARGILLE.OG0607
```
ID:                      OTHER:2
RI @ 587nm:              1.457587
Abbe Number:             57.19833
ΔPgF:                    0.008
TCE (÷1e-6):             700.0
Density:                 0.878g/m³
Valid wavelengths:       0.32μm to 1.55μm
Reference Temp:          25°C
```
"""
const OG0607 = Glass(GlassID(OTHER, 2), -2, 1.44503, 0.0044096, -2.85878e-5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.32, 1.55, -0.0009083144750540808, 0.0, 0.0, 0.0, 0.0, 0.0, 25.0, 0.008, -1.0, -1.0, 700.0, -1.0, 0, -1.0, [SVector(0.32, 0.15, 10.0), SVector(0.365, 0.12, 100.0), SVector(0.4047, 0.42, 100.0), SVector(0.480, 0.78, 100.0), SVector(0.4861, 0.79, 100.0), SVector(0.5461, 0.86, 100.0), SVector(0.5893, 0.90, 100.0), SVector(0.6328, 0.92, 100.0), SVector(0.6439, 0.90, 100.0), SVector(0.6563, 0.92, 100.0), SVector(0.6943, 0.98, 100.0), SVector(0.840, 0.99, 100.0), SVector(0.10648, 0.61, 100.0), SVector(0.1300, 0.39, 100.0), SVector(0.1550, 0.11, 100.0)], 1.457587, -1.0, -1.0, 0, 57.19833, 0, 0.878, -1)

"""    CARGILLE.OG081160
```
ID:                      OTHER:3
RI @ 587nm:              1.515549
Abbe Number:             36.82493
ΔPgF:                    0.014
TCE (÷1e-6):             700.0
Density:                 1.11g/m³
Valid wavelengths:       0.32μm to 1.55μm
Reference Temp:          25°C
```
"""
const OG081160 = Glass(GlassID(OTHER, 3), -2, 1.49614, 0.00692199, -8.07052e-5, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.32, 1.55, -0.000885983052189022, 0.0, 0.0, 0.0, 0.0, 0.0, 25.0, 0.014, -1.0, -1.0, 700.0, -1.0, 0, -1.0, [SVector(0.32, 0.04, 100.0), SVector(0.365, 0.13, 100.0), SVector(0.4047, 0.26, 100.0), SVector(0.480, 0.48, 100.0), SVector(0.4861, 0.49, 100.0), SVector(0.5461, 0.60, 100.0), SVector(0.5893, 0.68, 100.0), SVector(0.6328, 0.71, 100.0), SVector(0.6439, 0.73, 100.0), SVector(0.6563, 0.74, 100.0), SVector(0.6943, 0.76, 100.0), SVector(0.840, 0.83, 100.0), SVector(0.10648, 0.86, 100.0), SVector(0.1300, 0.89, 100.0), SVector(0.1550, 0.90, 100.0)], 1.515549, -1.0, -1.0, 0, 36.82493, 0, 1.11, -1)

end

const OTHER_GLASSES = [CARGILLE.OG0607, CARGILLE.OG0608, CARGILLE.OG081160]
const OTHER_GLASS_NAMES = ["CARGILLE.OG0607", "CARGILLE.OG0608", "CARGILLE.OG081160"]


text:
using Base.Test
using MIPVerify: MaskedReLU, mip

@testset "masked_relu.jl" begin

@testset "MaskedReLU" begin
    mask = rand(MersenneTwister(0), [-1, 0, 1], 10)
    @testset "Initialize without tightening algorithm" begin
        p = MaskedReLU(mask)
        @test isnull(p.tightening_algorithm)
    end

    @testset "Initialize with tightening algorithm" begin
        p = MaskedReLU(mask, mip)
        @test !isnull(p.tightening_algorithm)
        @test get(p.tightening_algorithm) == mip
    end

    @testset "Base.show" begin
        p = MaskedReLU(mask)
        io = IOBuffer()
        Base.show(io, p)
        @test String(take!(io)) == "MaskedReLU with expected input size (10,). (4 zeroed, 4 as-is, 2 rectified)."
    end
end

end


text:
using LinearAlgebra
function Vectorout(x1,y1,z1,x2,y2,z2,x3,y3,z3,x4,y4,z4)
    v = zeros(Float64,4,3)
    AB = [x2-x1;y2-y1;z2-z1]
    AC = [x3-x1;y3-y1;z3-z1]
    na = LinearAlgebra.cross(AB,AC)
    DA = [x1-x4;y1-y4;z1-z4]
    j1 = dot(na,DA)
    if j1 > 0
        n1 = na
    else
        n1 = -na
    end
    union1 = n1/sqrt(n1[1]^2+n1[2]^2+n1[3]^2)
    v[1,1] = union1[1]
    v[1,2] = union1[2]
    v[1,3] = union1[3]
    AB = [x2-x1;y2-y1;z2-z1]
    AD = [x4-x1;y4-y1;z4-z1]
    nb = LinearAlgebra.cross(AB,AD)
    CA = [x1-x3;y1-y3;z1-z3]
    j2 = dot(nb,CA)
    if j2 > 0
        n2 = nb
    else
        n2 = -nb
    end
    union2 = n2/sqrt(n2[1]^2+n2[2]^2+n2[3]^2)
    v[2,1] = union2[1]
    v[2,2] = union2[2]
    v[2,3] = union2[3]
    BC = [x3-x2;y3-y2;z3-z2]
    BD = [x4-x2;y4-y2;z4-z2]
    nc = LinearAlgebra.cross(BC,BD)
    AB = [x2-x1;y2-y1;z2-z1]
    j3 = dot(nc,AB)
    if j3 > 0
        n3 = nc
    else
        n3 = -nc
    end
    union3 = n3/sqrt(n3[1]^2+n3[2]^2+n3[3]^2)
    v[3,1] = union3[1]
    v[3,2] = union3[2]
    v[3,3] = union3[3]
    DA = [x1-x4;y1-y4;z1-z4]
    DC = [x3-x4;y3-y4;z3-z4]
    nd = LinearAlgebra.cross(DA,DC)
    BD = [x4-x2;y4-y2;z4-z2]
    j4 = dot(nd,BD)
    if j4 > 0
        n4 = nd
    else
        n4 = -nd
    end
    union4 = n4/sqrt(n4[1]^2+n4[2]^2+n4[3]^2)
    v[4,1] = union4[1]
    v[4,2] = union4[2]
    v[4,3] = union4[3]

    return v
end

max_stars_repo_path: src/single-core juSFEM/vectorout.jl

max_stars_repo_name: Elliothuo/juSFEM

max_issues_repo_path: src/single-core juSFEM/vectorout.jl

max_issues_repo_name: Elliothuo/juSFEM

max_forks_repo_path: src/single-core juSFEM/vectorout.jl

max_forks_repo_name: Elliothuo/juSFEM


text:
println("hello world")
println("this doesn't work")
println("hopefully this works!")
println("La terre nous aimait un peu je me souviens")


text:
#####################################################################################################
# Define the basic type system
#####################################################################################################
"Most abstract supertype for all SEMs"
abstract type AbstractSem end

"Supertype for all single SEMs, e.g. SEMs that have the fields `observed`, `imply`, loss, diff"
abstract type AbstractSemSingle <: AbstractSem end

"Supertype for all collections of multiple SEMs"
abstract type AbstractSemCollection <: AbstractSem end

"Supertype for all loss functions of SEMs. If you want to implement a custom loss function, it should be of this type."
abstract type SemLossFunction end

"""
    SemLoss(args...; ...)

Constructs the loss field of a SEM. Can contain multiple `SemLossFunction`s, the model is optimized over their sum.
See also [`SemLossFunction`](@ref).

# Arguments
- `args...`: Multiple `SemLossFunction`s.

# Examples
```julia
my_ml_loss = SemML(...)
my_ridge_loss = SemRidge(...)
my_loss = SemLoss(SemML, SemRidge)
```
"""
mutable struct SemLoss{F <: Tuple, T, FT, GT, HT}
    functions::F
    weights::T

    objective::FT
    gradient::GT
    hessian::HT
end

function SemLoss(functions...; loss_weights = nothing, parameter_type = Float64, kwargs...)

    n_par = length(gradient(functions[1]))
    !isnothing(loss_weights) || (loss_weights = Tuple(nothing for _ in 1:length(functions)))

    return SemLoss(
        functions,
        loss_weights,

        zeros(parameter_type, 1),
        zeros(parameter_type, n_par),
        zeros(parameter_type, n_par, n_par))
end

objective(loss::SemLoss) = loss.objective
gradient(loss::SemLoss) = loss.gradient
hessian(loss::SemLoss) = loss.hessian

"""
Supertype of all objects that can serve as the diff field of a SEM.
Connects the SEM to its optimization backend and controls options like the optimization algorithm.
If you want to connect the SEM package to a new optimization backend, you should implement a subtype of SemDiff.
"""
abstract type SemDiff end

"""
Supertype of all objects that can serve as the observed field of a SEM.
Pre-processes data and computes sufficient statistics for example.
If you have a special kind of data, e.g. ordinal data, you should implement a subtype of SemObs.
"""
abstract type SemObs end

"""
Supertype of all objects that can serve as the imply field of a SEM.
Computed model-implied values that should be compared with the observed data to find parameter estimates,
e. g. the model implied covariance or mean.
If you would like to implement a different notation, e.g. LISREL, you should implement a subtype of SemImply.
"""
abstract type SemImply end

"Subtype of SemImply for all objects that can serve as the imply field of a SEM and use some form of symbolic precomputation."
abstract type SemImplySymbolic <: SemImply end

"""
    Sem(;observed = SemObsCommon, imply = RAM, loss = (SemML,), diff = SemDiffOptim, kwargs...)

Constructor for the basic `Sem` type.
All additional kwargs are passed down to the constructors for the fields.

# Arguments
- `observed`: isa `SemObs` or a constructor.
- `imply`: isa `SemImply` or a constructor.
- `loss`: Tuple of objects that are `SemLossFunction`s or constructors.
- `observed`: isa `SemObs` or a constructor.

Returns a struct of type Sem with fields
- `observed::SemObs`: Stores observed data, sample statistics, etc. See also [`SemObs`](@ref).
- `imply::SemImply`: Computes Σ, μ, etc. See also [`SemImply`](@ref).
- `loss::SemLoss`: Computes the objective and gradient of a sum of loss functions. See also [`SemLoss`](@ref).
- `diff::SemDiff`: Connects the model to the optimizer. See also [`SemDiff`](@ref).
"""
mutable struct Sem{O <: SemObs, I <: SemImply, L <: SemLoss, D <: SemDiff} <: AbstractSemSingle
    observed::O
    imply::I
    loss::L
    diff::D
end

function (model::Sem)(par, F, G, H)
    model.imply(par, F, G, H, model)
    model.loss(par, F, G, H, model)
end

function (loss::SemLoss)(par, F, G, H, model)
    for lossfun in loss.functions lossfun(par, F, G, H, model) end
    if H
        loss.hessian .= 0.0
        for (lossfun, c) in zip(loss.functions, loss.weights)
            if isnothing(c)
                loss.hessian .+= hessian(lossfun)
            else
                loss.hessian .+= c*hessian(lossfun)
            end
        end
    end
    if G
        loss.gradient .= 0.0
        for (lossfun, c) in zip(loss.functions, loss.weights)
            if isnothing(c)
                loss.gradient .+= gradient(lossfun)
            else
                loss.gradient .+= c*gradient(lossfun)
            end
        end
    end
    if F
        loss.objective[1] = 0.0
        for (lossfun, c) in zip(loss.functions, loss.weights)
            if isnothing(c)
                loss.objective[1] += objective(lossfun)[1]
            else
                loss.objective[1] += c*objective(lossfun)[1]
            end
        end
    end
end

#####################################################################################################
# automatic differentiation
#####################################################################################################

struct SemFiniteDiff{O <: SemObs, I <: SemImply, L <: SemLoss, D <: SemDiff, G} <: AbstractSemSingle
    observed::O
    imply::I
    loss::L
    diff::D
    has_gradient::G
end

function (model::SemFiniteDiff)(par, F, G, H)

    if H
        model.loss.hessian .= FiniteDiff.finite_difference_hessian(x -> objective!(model, x), par)
    end

    if model.has_gradient
        model.imply(par, F, G, false, model)
        model.loss(par, F, G, false, model)
    else

        if G
            model.loss.gradient .= FiniteDiff.finite_difference_gradient(x -> objective!(model, x), par)
        end

        if F
            model.imply(par, F, false, false, model)
            model.loss(par, F, false, false, model)
        end

    end

end

struct SemForwardDiff{O <: SemObs, I <: SemImply, L <: SemLoss, D <: SemDiff, G} <: AbstractSemSingle
    observed::O
    imply::I
    loss::L
    diff::D
    has_gradient::G
end

function (model::SemForwardDiff)(par, F, G, H)

    if H
        model.loss.hessian .= ForwardDiff.hessian(x -> objective!(model, x), par)
    end

    if model.has_gradient
        model.imply(par, F, G, false, model)
        model.loss(par, F, G, false, model)
    else

        if G
            model.loss.gradient .= ForwardDiff.gradient(x -> objective!(model, x), par)
        end

        if F
            model.imply(par, F, false, false, model)
            model.loss(par, F, false, false, model)
        end

    end

end

#####################################################################################################
# ensemble models
#####################################################################################################

struct SemEnsemble{N, T <: Tuple, V <: AbstractVector, D, I, FT, GT, HT} <: AbstractSemCollection
    n::N
    sems::T
    weights::V
    diff::D
    identifier::I

    objective::FT
    gradient::GT
    hessian::HT
end

function SemEnsemble(models...; diff = SemDiffOptim, weights = nothing, parameter_type = Float64, kwargs...)
    n = length(models)

    # default weights

    if isnothing(weights)
        nobs_total = sum(n_obs.(models))
        weights = [n_obs(model)/nobs_total for model in models]
    end

    # check identifier equality
    id = identifier(models[1])
    for model in models
        if id != identifier(model)
            throw(ErrorException("The identifier of your models do not match. \n
            Maybe you tried to specify models of an ensemble via ParameterTables. \n
            In that case, you may use RAMMatrices instead."))
        end
    end

    npar = n_par(models[1])

    # diff
    if !isa(diff, SemDiff)
        diff = diff(;kwargs...)
    end

    return SemEnsemble(
        n,
        models,
        weights,
        diff,
        id,

        zeros(parameter_type, 1),
        zeros(parameter_type, npar),
        zeros(parameter_type, npar, npar))
end

function (ensemble::SemEnsemble)(par, F, G, H)

    if H ensemble.hessian .= 0.0 end
    if G ensemble.gradient .= 0.0 end
    if F ensemble.objective .= 0.0 end

    for (model, weight) in zip(ensemble.sems, ensemble.weights)
        model(par, F, G, H)
        if H ensemble.hessian .+= weight*hessian(model) end
        if G ensemble.gradient .+= weight*gradient(model) end
        if F ensemble.objective .+= weight*objective(model) end
    end

end

n_models(ensemble::SemEnsemble) = ensemble.n
models(ensemble::SemEnsemble) = ensemble.sems
weights(ensemble::SemEnsemble) = ensemble.weights
diff(ensemble::SemEnsemble) = ensemble.diff

#####################################################################################################
# additional methods
#####################################################################################################

observed(model::Sem) = model.observed
imply(model::Sem) = model.imply
loss(model::Sem) = model.loss
diff(model::Sem) = model.diff

observed(model::SemForwardDiff) = model.observed
imply(model::SemForwardDiff) = model.imply
loss(model::SemForwardDiff) = model.loss
diff(model::SemForwardDiff) = model.diff
has_gradient(model::SemForwardDiff) = model.has_gradient

observed(model::SemFiniteDiff) = model.observed
imply(model::SemFiniteDiff) = model.imply
loss(model::SemFiniteDiff) = model.loss
diff(model::SemFiniteDiff) = model.diff
has_gradient(model::SemFiniteDiff) = model.has_gradient

#####################################################################################################
# gradient, objective, hessian helpers
#####################################################################################################

objective(lossfun::SemLossFunction) = lossfun.objective
gradient(lossfun::SemLossFunction) = lossfun.gradient
hessian(lossfun::SemLossFunction) = lossfun.hessian

objective(model::AbstractSem) = model.loss.objective
gradient(model::AbstractSem) = model.loss.gradient
hessian(model::AbstractSem) = model.loss.hessian

objective(model::SemEnsemble) = model.objective
gradient(model::SemEnsemble) = model.gradient
hessian(model::SemEnsemble) = model.hessian

function objective!(model::AbstractSem, parameters)
    model(parameters, true, false, false)
    return objective(model)[1]
end

function gradient!(model::AbstractSem, parameters)
    model(parameters, false, true, false)
    return gradient(model)
end

function gradient!(grad, model::AbstractSem, parameters)
    model(parameters, false, true, false)
    copyto!(grad, gradient(model))
    return gradient(model)
end

function hessian!(model::AbstractSem, parameters)
    model(parameters, false, false, true)
    return hessian(model)
end

function hessian!(hessian, model::AbstractSem, parameters)
    model(parameters, false, false, true)
    copyto!(hessian, hessian(model))
    return hessian(model)
end

function objective_gradient!(model::AbstractSem, parameters)
    model(parameters, true, true, false)
    return objective(model)[1], copy(gradient(model))
end

function objective_gradient!(grad, model::AbstractSem, parameters)
    model(parameters, true, true, false)
    copyto!(grad, gradient(model))
    return objective(model)[1]
end


text:
module PointSampler
using LinearAlgebra: norm
using Memoize: @memoize
using Parameters: @with_kw

export monte_carlo_th, MonteCarloThDesign

const RVec = AbstractVector{<:Real}
const RVecArr = AbstractVector{<:RVec};

abstract type PointIterator end

include("utils.jl")

include("monte_carlo_thresholded.jl");

end


text:
using RecipesBase

@recipe plot(ts::TSeries...) = begin
    for t in ts
        @series begin
            (Float64[mitrange(t)...], t.values)
        end
    end
end


text:
module Introspection

using InteractiveUtils

import MyterialColors: orange, grey_dark, light_green

import Term:
    highlight,
    escape_brackets,
    join_lines,
    unescape_brackets,
    split_lines,
    do_by_line,
    expr2string

import ..Panels: Panel
import ..Dendograms: Dendogram
import ..Trees: Tree

include("_inspect.jl")

export inspect, typestree, expressiontree

# ---------------------------------------------------------------------------- #
#                                TYPES HIERARCHY                               #
# ---------------------------------------------------------------------------- #

function typestree(io::IO, T::DataType)
    return print(
        io,
        Panel(
            Tree(T);
            title = "Types hierarchy",
            style = "blue dim",
            title_style = orange * " default",
            title_justify = :right,
            fit = true,
            padding = (1, 4, 1, 1),
        ),
    )
end
typestree(T::DataType) = typestree(stdout, T)

function expressiontree(io::IO, e::Expr)
    _expr = expr2string(e)
    tree = Tree(e)

    return print(
        io,
        Panel(
            tree;
            title = _expr,
            title_style = "$light_green default bold",
            title_justify = :center,
            style = grey_dark,
            fit = tree.measure.w > 88,
            width = max(tree.measure.w, 88),
            subtitle = "inspect",
            subtitle_justify = :right,
            justify = :center,
            padding = (1, 1, 1, 1),
        ),
    )
end
expressiontree(e::Expr) = expressiontree(stdout, e)

# ---------------------------------------------------------------------------- #
#                                EXPR. DENDOGRAM                               #
# ---------------------------------------------------------------------------- #

function inspect(io::IO, expr::Expr)
    _expr = expr2string(expr)
    dendo = Dendogram(expr)

    return print(
        io,
        Panel(
            dendo;
            title = _expr,
            title_style = "$light_green default bold",
            title_justify = :center,
            style = grey_dark,
            fit = true,
            # width=dendo.measure.w,
            subtitle = "inspect",
            subtitle_justify = :right,
            justify = :center,
            padding = (1, 1, 1, 1),
        ),
    )
end

end


text:
function parse_sequence(str::AbstractString)
    sequence = split.(split(str, ";"), ",")
    sequence = [uppercase(join(first.(s))) for s in sequence]
    return join(sequence, "#")
end

function parse_strokes(str::AbstractString)
    lines = filter(x -> x[1] != '#', split(str))
    sequence = parse_sequence(lines[1])
    data = map(x -> parse.(Float64, split(x, ",")), lines[2:end])
    return sequence, reduce(hcat, data)
end

function load_strokes(fn::AbstractString; transform=false, jitter=false)
    sequence, data = read(fn, String) |> parse_strokes
    transform && transform_strokes!(data)
    jitter && add_jitter!(data)
    return sequence, data
end

function transform_strokes!(data::AbstractMatrix;
                            invert_x::Bool=false, invert_y::Bool=true,
                            in_scale::Real=176.389, out_scale::Real=5.0)
    data .*= (out_scale / in_scale)
    if invert_x data[1,:] .*= -1 end
    if invert_y data[2,:] .*= -1 end
    min_bounds, max_bounds = minimum(data, dims=2), maximum(data, dims=2)
    mid_x, mid_y = 0.5*(min_bounds[1] + max_bounds[1]), 0.5*(min_bounds[2] + max_bounds[2])
    data[1,:] .+= (out_scale / 2 - mid_x)
    data[2,:] .+= (out_scale / 2 - mid_y)
    return data
end

function add_jitter!(data::AbstractMatrix; x_std::Real=0.05, y_std::Real=0.05)
    data[1,:] .+= (randn(size(data)[2]) .* x_std)
    data[2,:] .+= (randn(size(data)[2]) .* y_std)
    return data
end

function load_dataset(dir::AbstractString; transform=true, jitter=false)
    fns = readdir(dir)
    fns = filter(s -> match(r".*\.txt", s) !== nothing, fns)
    dataset = load_strokes.(joinpath.(dir, fns); transform=transform, jitter=jitter)
    return dataset
end


text:
"""
    PyIterable{T=PyObject}(o)

Wrap the Python object `o` into a Julia object which iterates values of type `T`.
"""
struct PyIterable{T}
    ref::PyRef
    PyIterable{T}(o) where {T} = new{T}(PyRef(o))
end
PyIterable(o) = PyIterable{PyObject}(o)
export PyIterable

ispyreftype(::Type{<:PyIterable}) = true
pyptr(x::PyIterable) = pyptr(x.ref)
Base.unsafe_convert(::Type{CPyPtr}, x::PyIterable) = checknull(pyptr(x))
C.PyObject_TryConvert__initial(o, ::Type{T}) where {T<:PyIterable} =
    C.putresult(T(pyborrowedref(o)))

Base.length(x::PyIterable) = Int(pylen(x))

Base.IteratorSize(::Type{<:PyIterable}) = Base.SizeUnknown()

Base.IteratorEltype(::Type{<:PyIterable}) = Base.HasEltype()

Base.eltype(::Type{PyIterable{T}}) where {T} = T

function Base.iterate(x::PyIterable{T}, it = pyiter(PyRef, x)) where {T}
    vo = C.PyIter_Next(it)
    if !isnull(vo)
        r = C.PyObject_Convert(vo, T)
        C.Py_DecRef(vo)
        checkm1(r)
        (takeresult(T), it)
    elseif C.PyErr_IsSet()
        pythrow()
    else
        nothing
    end
end


text:
import JuliaDB.ML

@testset "feature extraction" begin
    @testset "schema" begin
        @test ML.schema(1:10).stat == fit!(Variance(), 1:10)
        x = repeat([1,2], inner=5)
        @test ML.schema(x, ML.Categorical).stat == fit!(CountMap(Int), x) == ML.schema(PooledArray(x)).stat
        m = ML.schema(Vector{Union{Int,Missing}}(x))
        @test m isa ML.Maybe
        @test m.feature.stat == ML.schema(x).stat
        t = table(PooledArray([1,2,1,2]), [1,2,3,4], names=[:x, :y])
        sch = ML.schema(t)
        @test sch[:x].stat == ML.schema(PooledArray([1,2,1,2])).stat
        @test sch[:y].stat == ML.schema([1,2,3,4]).stat
    end

    @testset "featuremat" begin
        @test ML.featuremat([1,3,5]) == [-1.0, 0, 1.0]'
        Δ = 0.5 / std([1,2])
        @test ML.featuremat([1,2,missing]) ≈ [0 0 1; -Δ Δ 0]
        @test isempty(ML.featuremat(["x","y","x"]))
        @test ML.featuremat(PooledArray(["x","y","x"])) == [1  0  1; 0  1  0]

        t = table(PooledArray([1,1,2,2]), [1,2,3,4], [1,2,3,4], ["x", "y", "z", "a"], names=[:a,:b,:c,:d])
        @test ML.featuremat(t) == collect(ML.featuremat(distribute(t, 2)))

        x = randn(100)
        @test ML.featuremat(x) ≈ ((x .- mean(x)) ./ std(x))'
    end

    @testset "Issue 149" begin
        x = rand(10), rand(1:4, 10), rand(1:4, 10)
        t  = table(x..., names = [:x1, :x2, :x3])
        td = table(x..., names = [:x1, :x2, :x3], chunks = 2)
        sch = ML.schema(t, hints=Dict(
            :x1 => ML.Continuous,
            :x2 => ML.Categorical,
            :x3 => ML.Categorical,
            )
        )
        schd = ML.schema(td, hints=Dict(
            :x1 => ML.Continuous,
            :x2 => ML.Categorical,
            :x3 => ML.Categorical,
            )
        )

        @test schd[:x1].stat.σ2 ≈ sch[:x1].stat.σ2
        @test schd[:x2].stat == sch[:x2].stat
        @test schd[:x3].stat == sch[:x3].stat
    end
end


text:
#!/usr/bin/env julia

using EchogramPyPlot
using Test

a = rand(100,100)
eg(a)