text:
import ..Highlights.Compiler: NULL_RANGE, Context, nullmatch
using REPL
# Utilities.

function julia_is_identifier(ctx::Context, prefix = '\0')
    s = ctx.source
    i = ctx.pos[]
    (i > ncodeunits(s)) && return NULL_RANGE
    (c, i) = iterate(s, i)
    if prefix !== '\0'
        (c === prefix && !(i > ncodeunits(s))) || return NULL_RANGE
        (c, i) = iterate(s, i)
    end
    Base.is_id_start_char(c) || return NULL_RANGE
    prev_i = i

    while !(i > ncodeunits(s))
        prev_i = i
        (c, i) = iterate(s, i)
        Base.is_id_char(c) || break
    end
    return ctx.pos[]:prevind(s, Base.is_id_char(c) ? i : prev_i)
end
julia_is_macro_identifier(ctx::Context) = julia_is_identifier(ctx, '@')
julia_is_iterp_identifier(ctx::Context) = julia_is_identifier(ctx, '$')

function julia_is_operator(ctx::Context)
    ind = ctx.pos[]
    while ind <= ncodeunits(ctx.source)
        str = SubString(ctx.source, ctx.pos[], ind)
        Base.isoperator(Symbol(str)) || break
        ind = nextind(ctx.source, ind)
    end
    return ind == ctx.pos[] ? NULL_RANGE : (ctx.pos[]:(ind - 1))
end

function julia_is_symbol(ctx::Context)
    s = ctx.source
    i = ctx.pos[]
    valid = i > 1 ? iterate(s, prevind(s, i))[1] in ('(', '{', '[', ' ', '\n') : true
    return valid ? julia_is_identifier(ctx, ':') : NULL_RANGE
end

function julia_is_method_call(ctx::Context)
    range = julia_is_identifier(ctx)
    range === NULL_RANGE && return range
    i = nextind(ctx.source, last(range))
    (i > ncodeunits(ctx.source)) && return NULL_RANGE
    (c, i) = iterate(ctx.source, i)
    return (c === '(' || c === '{') ? range : NULL_RANGE
end

function julia_is_string_macro(ctx::Context, count::Integer = 1)
    range = julia_is_identifier(ctx)
    range == NULL_RANGE && return NULL_RANGE
    s = ctx.source
    epos = min(ctx.pos[] + length(range), length(s))
    i = prevind(s, epos + 1)
    num = 0
    while num < count && !(i > ncodeunits(s))
        (c, i) = iterate(s, i)
        c === '"' ? (num += 1) : break
    end
    num == count ? (ctx.pos[]:prevind(s, i)) : NULL_RANGE
end
julia_is_triple_string_macro(ctx::Context) = julia_is_string_macro(ctx, 3)

# Julia Script Lexer.

@lexer JuliaLexer let
    keywords = [kw.keyword for kw in REPL.REPLCompletions.complete_keyword("")]
    char_regex = [
        raw"'(\\.|\\[0-7]{1,3}|\\x[a-fA-F0-9]{1,3}|",
        raw"\\u[a-fA-F0-9]{1,4}|\\U[a-fA-F0-9]{1,6}|",
        raw"[^\\\'\n])'",
    ]
    Dict(
        :name => "Julia",
        :description => "A lexer for Julia source code.",
        :aliases => ["julia", "jl"],
        :filenames => ["*.jl"],
        :mimetypes => ["text/x-julia", "application/x-julia"],
        :tokens => Dict(
            :root => [
                (r"\n"m, TEXT),
                (r"[^\S\n]+"m, TEXT),
                (r"#=", COMMENT_MULTILINE, :block_comments),
                (r"#.*$"m, COMMENT_SINGLE),
                (r"[\[\]{}(),;]", PUNCTUATION),

                (julia_is_symbol, STRING_CHAR),
                (r"\b(?<![_.])in\b", KEYWORD_PSEUDO),
                (r"\b(?<![_.])end\b", KEYWORD),
                (r"\b(?<![_.])(true|false)\b", KEYWORD_CONSTANT),
                (r"\b(?<![_.])(local|global|const)\b", KEYWORD_DECLARATION),
                (words(keywords, prefix = "\\b(?<![_.])", suffix = "\\b"), KEYWORD),

                (Regex(join(char_regex)), STRING_CHAR),

                (r"\"\"\"", STRING, :triple_strings),
                (r"\"", STRING, :strings),

                (julia_is_triple_string_macro, STRING_OTHER, :triple_string_macros),
                (julia_is_string_macro, STRING_OTHER, :string_macros),

                (r"`", STRING_BACKTICK, :commands),

                (julia_is_method_call, NAME_FUNCTION),
                (julia_is_identifier, NAME),
                (julia_is_macro_identifier, NAME_DECORATOR),

                (r"(\d+(_\d+)+\.\d*|\d*\.\d+(_\d+)+)([eEf][+-]?[0-9]+)?", NUMBER_FLOAT),
                (r"(\d+\.\d*|\d*\.\d+)([eEf][+-]?[0-9]+)?", NUMBER_FLOAT),
                (r"\d+(_\d+)+[eEf][+-]?[0-9]+", NUMBER_FLOAT),
                (r"\d+[eEf][+-]?[0-9]+", NUMBER_FLOAT),
                (r"0b[01]+(_[01]+)+", NUMBER_BIN),
                (r"0b[01]+", NUMBER_BIN),
                (r"0o[0-7]+(_[0-7]+)+", NUMBER_OCT),
                (r"0o[0-7]+", NUMBER_OCT),
                (r"0x[a-fA-F0-9]+(_[a-fA-F0-9]+)+", NUMBER_HEX),
                (r"0x[a-fA-F0-9]+", NUMBER_HEX),
                (r"\d+(_\d+)+", NUMBER_INTEGER),
                (r"\d+", NUMBER_INTEGER),

                (julia_is_operator, OPERATOR),

                (r"."ms, TEXT),
            ],
            :block_comments => [
                (r"[^#=]", COMMENT_MULTILINE),
                (r"#=", COMMENT_MULTILINE, :__push__),
                (r"=#", COMMENT_MULTILINE, :__pop__),
                (r"[=#]", COMMENT_MULTILINE),
            ],
            :commands => [
                (r"`", STRING_BACKTICK, :__pop__),
                (julia_is_iterp_identifier, STRING_INTERPOL),
                (r"(\$)(\()", (STRING_INTERPOL, PUNCTUATION), :in_interpol),
                (r".|\s"ms, STRING_BACKTICK),
            ],
            :strings => [
                (r"\"", STRING, :__pop__),
                (r"\\([\"\\'\$nrbtfav]|(x|u|U)[a-fA-F0-9]+|\d+)", STRING_ESCAPE),
                (julia_is_iterp_identifier, STRING_INTERPOL),
                (r"(\$)(\()", (STRING_INTERPOL, PUNCTUATION), :in_interpol),
                (r".|\s"ms, STRING),
            ],
            :triple_strings => [
                (r"\"\"\"", STRING, :__pop__),
                (r"\\([\"\\'\$nrbtfav]|(x|u|U)[a-fA-F0-9]+|\d+)", STRING_ESCAPE),
                (julia_is_iterp_identifier, STRING_INTERPOL),
                (r"(\$)(\()", (STRING_INTERPOL, PUNCTUATION), :in_interpol),
                (r".|\s"ms, STRING),
            ],
            :string_macros => [
                (r"\"", STRING_OTHER, :__pop__),
                (r"\\\"", STRING_OTHER),
                (r".|\s"ms, STRING_OTHER),
            ],
            :triple_string_macros => [
                (r"\"\"\"", STRING_OTHER, :__pop__),
                (r".|\s"ms, STRING_OTHER),
            ],
            :in_interpol => [
                (r"\(", PUNCTUATION, :__push__),
                (r"\)", PUNCTUATION, :__pop__),
                :root,
            ],
        ),
    )
end

# Julia Console Lexer.

function julia_repl_splitter(ctx::Context)
    s = ctx.source
    i = ctx.pos[]
    while !(i > ncodeunits(s))
        (c, i) = iterate(s, i)
        if c == '\n'
            @label NEWLINE
            count = 0
            while !(i > ncodeunits(s))
                (c, i) = iterate(s, i)
                if c == '\n'
                    # Multiple blank lines in a row.
                    @goto NEWLINE
                elseif c == ' '
                    # Leading indent.
                    count += 1
                else
                    # When 7-space indented (width of `julia> `) we are still in the input,
                    # but when the indent is less than that we are now in the output.
                    count > 6 ? break : @goto FINISHED
                end
            end
        end
    end
    @label FINISHED
    return ctx.pos[]:prevind(s, i, 2)
end

@lexer JuliaConsoleLexer Dict(
    :name => "Julia Console",
    :description => "A lexer for Julia REPL sessions.",
    :aliases => ["jlcon"],
    :tokens => Dict(
        :root => [
            (r"^(julia|shell|help\?)> "m, NUMBER, :source),
            (r"."ms, TEXT),
        ],
        :source => [
            (julia_repl_splitter, JuliaLexer, :__pop__)
        ],
    ),
)


text:
using FFTW

function recenter_block(x::Block)
    return Block(x.data .- 128)
end

function decenter_block(x::Block)
    return Block(x.data .+ 128)
end

function dct(g::Block)::Matrix{Float32}
    rg = recenter_block(g).data
    return FFTW.dct(rg)
end

function invdct(g::Matrix{Float32})::Block
    rg = round.(FFTW.idct(g))
    b = decenter_block(Block(rg))
    return b
end

ext: jl


text:
#
# Copyright (c) 2021 Tobias Thummerer, Lars Mikelsons, Josef Kircher
# Licensed under the MIT license. See LICENSE file in the project root for details.
#

using FMI

pathToFMU = joinpath(dirname(@__FILE__), "../model/SpringFrictionPendulum1D.fmu")

myFMU = fmiLoad(pathToFMU)

fmiInstantiate!(myFMU; loggingOn=true)

fmiSetupExperiment(myFMU, 0.0)

fmiEnterInitializationMode(myFMU)
fmiExitInitializationMode(myFMU)

dt = 0.01
t_start = 0.0
t_stop = 8.0

data = fmiSimulateCS(myFMU, dt, t_start, t_stop, ["mass.s", "mass.v"])
fmiPlot(data)

fmiUnload(myFMU)


text:
using Base.Test

using Tick: Dag, node!, onfire!
using Tick: EvalDfs, EvalTSort, fire, tsort
using Tick: latest
using Tick: BVal, combine!, tr!
using Tick: add!, first!, last!, nth!, win!

# TODO FactCheck.jl

const testdir = dirname(@__FILE__)

@testset "Tick" begin
    include(joinpath(testdir, "evaldfs.jl"))
    include(joinpath(testdir, "evaltsort.jl"))
    include(joinpath(testdir, "combine.jl"))
    include(joinpath(testdir, "ops.jl"))
end


text:
module TensArrays

export TensArray

"""
    TensArray{T,N} <: AbstractArray{T,N}
    TensArray(x, cs, rs)

An Array with co- and contravariant dimensions.

TensArray wraps an AbstractArray with metadata that allows it to
be interpreted as a tensor.  Refer to the manual of Tensars.jl for
details.
"""
struct TensArray{T,N} <: AbstractArray{T,N}
    parent::AbstractArray{T,N}
    cs::Dims
    rs::Dims

    function TensArray(x::AbstractArray{T,N}, cs, rs) where T where N
        if match_shape(size(x), cs, rs)
            new{T,N}(Array(x), cs, rs)
        else
            "Array dimensions do not match tensor dimensions" |>
                DimensionMismatch |> throw
        end
    end
end

# 1D TensArray can be a matrix
function match_shape(as::Dims, cs::Dims, rs::Dims)
    cpas = cumprod([as...])
    prod(as) == prod(cs)*prod(rs) &&
        issubseq(cpas, cumprod([cs..., rs...])) &&
        (as ≠ () || cs == rs == ()) &&
        (cs == () || prod(cs) ∈ cpas)
end

function issubseq(xs, ys)
    if isempty(xs)
        true
    elseif isempty(ys)
        false
    elseif first(xs) == first(ys)
        issubseq(xs[2:end], ys[2:end])
    else
        issubseq(xs, ys[2:end])
    end
end

parent(x::TensArray) = x.parent

Base.size(x::TensArray) = size(parent(x))
Base.getindex(x::TensArray, j...) = x.parent[j...]
Base.setindex!(x::TensArray, v, j) = (x.parent[j] = v)

Base.similar(x::TensArray, T, dims) =
    TensArray(similar(parent(x), T, dims), x.cs, x.rs)

# TODO show with rowsize and colsize

# matmul

function Base.:*(x::TensArray, y::TensArray)
    a = parent(x)*parent(y)
    x.rs == y.cs ? TensArray(a, x.cs, y.rs) : a
end

# TODO reshaping

# Broadcasting, copied from the ArrayAndChar magic pixie dust.  For
# now, the tensor shape is discarded except for scalars and exactly
# matching matrices.

Base.BroadcastStyle(::Type{<:TensArray}) = Broadcast.ArrayStyle{TensArray}()

function Base.similar(bc::Broadcast.Broadcasted{Broadcast.ArrayStyle{TensArray}}, ::Type{ElType}) where ElType
    # Find the TensArray
    A = find_aac(bc)
    # Create the output
    pA = similar(Array{ElType}, axes(bc))
    (A == false) ? pA : TensArray(pA, A.cs, A.rs)
end

find_aac(bc::Base.Broadcast.Broadcasted) = find_aac(bc.args)
find_aac(::Tuple{}) = true
find_aac(args::Tuple) = find_aac(find_aac(args[1]), Base.tail(args))
find_aac(x) = isempty(size(x))
find_aac(a::TensArray) = a
find_aac(a, rest) = merge_aac(a, find_aac(rest))

merge_aac(a::Bool, b::TensArray) = a && b
merge_aac(b::TensArray, a::Bool) = a && b
merge_aac(a::Bool, b::Bool) = a && b
merge_aac(a::TensArray, b::TensArray) =
    (a.cs, a.rs) == (b.cs, b.rs) && a

end # module


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

using Markdown
using InteractiveUtils

# ╔═╡ c79638a6-c80d-11eb-3254-b33f7a382f78
begin
	using Pkg
	Pkg.activate(joinpath(@__DIR__, ".."))
	using Random
	using Agents
	using DataFrames
	using ImageFiltering
	using CairoMakie
	using InteractiveDynamics
end

# ╔═╡ 7852fa86-a513-49ac-92b1-7f0a47cedc17
function neighbors_at(n::Integer)
    k = 2 * n + 1
    kern = zeros(Int, k, k)
    center = k ÷ 2 + 1
    dist = map(CartesianIndices(kern)) do I
        i, j = Tuple(I)
        center - (abs(i - center) + abs(j - center))
    end
    dist[center, center] = 0
    map(findall(dist .> 0)) do I
        (I[1] - center, I[2] - center)
    end
end

# ╔═╡ 6ed9031f-8754-46bf-bf92-304d1ea462b8
const default_parameters = (#
    energy_miss=0.025,
    age_init=168,
    age_reproduce=504,
    age_old=600,
    age_die=720,
    pr_reproduce=0.176,
    pr_egg_death=0.0025,
    pr_old_death=0.04,
    offspring_max=12,
    offspring_min=5,
    energy_max=1.0,
    energy_transfer=0.1,
    energy_consume=0.025,
    energy_move=0.2,
    energy_reproduce=0.8,
    move_directions=neighbors_at(2)
)

# ╔═╡ 55d77475-590d-4c21-a47b-2012fba13d8f
Base.@kwdef mutable struct BPH <: AbstractAgent
    id::Int
    pos::Dims{2}
    energy::Float16
    age::Int16
end

# ╔═╡ cc080461-80c3-40cd-bc6e-d1c46f017b9a
function model_step!(model)
	T = eltype(model.food)
    @. model.food = min(model.food + model.food * T(0.008), one(T))
end

# ╔═╡ c7c02074-b970-425a-9586-48aa34e62c27
function agent_step!(agent, model)
    # position
    x, y = agent.pos

    # Older
    agent.age = agent.age + one(agent.age)

    # Step wise energy loss
    agent.energy = agent.energy - (agent.age ≥ model.age_init) * model.energy_consume

    # Move conditionally
    if agent.age ≥ model.age_init && agent.energy ≥ model.energy_move
        walk!(agent, rand(model.rng, model.move_directions), model)
    end

    # Eat conditionally
    if model.food[x, y] > 0 && agent.age ≥ model.age_init
        transfer = min(#
            model.energy_transfer,
            model.food[x, y],
            model.energy_max - agent.energy,
        )
        model.food[x, y] -= transfer
        agent.energy += transfer
    end

    # Reproduce conditionally
    if (
        agent.age ≥ model.age_reproduce && # Old enough
        agent.energy ≥ model.energy_reproduce && # Energy requirement
        rand(model.rng) ≤ model.pr_reproduce # Have RNG Jesus by your side
    )
        nb_offspring = rand(model.rng, (model.offspring_min):(model.offspring_max))
        for _ in 1:nb_offspring
            id = nextid(model)
            agent = BPH(;
                id=id,
                pos=agent.pos,
                energy=0.4,
                age=0,
            )
            add_agent_pos!(agent, model)
        end
        agent.energy -= one(agent.energy) / 10
    end

    # Die conditionally
    if (agent.energy ≤ 0) || # Exausted
       (agent.age ≥ model.age_die) || # Too old
       (rand(model.rng) ≤ model.pr_killed[x, y]) || # killed
       (
           model.age_die > agent.age ≥ model.age_old && # Old
           rand(model.rng) ≤ model.pr_old_death # And weak
       ) ||
       (
           agent.age < model.age_init && # Young
           rand(model.rng) ≤ model.pr_egg_death # And weak
       ) # then
        kill_agent!(agent, model)
        return nothing
    end
    # End of agent step
end

# ╔═╡ a7a6cbab-5b38-42d3-9ede-1e2b5550f2f1
function init_model(; envmap, pr_killed, seed, scheduler, kwargs...)
	nb_bph_init = 200
    rng = MersenneTwister(seed)
    food = collect(transpose(envmap))
    pr_killed = imfilter(isnan.(food) * pr_killed, Kernel.gaussian(2.5))

    # PROPERTIES

    params = merge(default_parameters, kwargs)
    properties = (#
        food=food,
        total_bph=nb_bph_init,
        death_natural=0,
        death_predator=0,
        pr_killed=pr_killed,
        pr_killed_positions=convert.(Tuple, findall(!iszero, pr_killed)),
        energy_full=1.0 - params.energy_transfer,
        params...,
    )

    # MODEL

    space = GridSpace(size(food); periodic=false)
    model = ABM(BPH, space; scheduler=scheduler, properties=properties, rng=rng)

    # AGENTS CREATION
    positions = Iterators.product(1:5, 1:5)
    positions = filter(pos -> !isnan(food[pos...]), collect(positions))
    for _ in 1:nb_bph_init
        bph = BPH(; #
            id=nextid(model),
            pos=rand(model.rng, positions),
            energy=rand(model.rng, 0.4:0.01:0.6),
            age=rand(model.rng, 0:300),
        )
        add_agent_pos!(bph, model)
    end

    # RETURN
    return model
end

# ╔═╡ 79076b76-917f-439c-b954-8046867614c0
heatkwargs=(
	nan_color=RGBAf0(1.0, 1.0, 0.0, 0.5),
	colormap=[RGBAf0(0, 1.0, 0, i) for i in 0:0.01:1],
	colorrange=(0, 1),
)

# ╔═╡ 11b14ef0-3554-498c-8e14-aafd5d026a49
adata, mdata = let food(model) = count(model.food .≥ 0.5),
	bph(agent) = agent.energy > 0
    [(bph, count)], [food]
end

# ╔═╡ 11362a40-572d-4019-9620-baecb02b01f9
begin
	function envmap(T::Type, d::Integer)
		env = ones(T, 100, 100)
		xstart, xend = (100 - d) ÷ 2, (100 + d) ÷ 2
		env[:, xstart+1:xend] .= convert(T, NaN)
		return env
	end
	envmap(d::Integer) = envmap(Float16, d)
end

# ╔═╡ 6cc30e3d-ca13-4ffb-b623-fcdb7c31936c
function run_and_collect_data(; seed, scheduler)
	model = init_model(;
		envmap=envmap(20),
		pr_killed=0.075,
		seed=seed,
		scheduler=scheduler)
	adf, mdf = run!(model, agent_step!, model_step!, 2880;
		adata=adata, mdata=mdata)
	leftjoin(adf, mdf, on=:step)
end

# ╔═╡ fe5717d3-dd4a-470a-82bb-453e9867dc3d
function run_and_collect_data2(; seed, scheduler)
	model = init_model(;
		envmap=envmap(20),
		pr_killed=0.075,
		seed=seed,
		scheduler=scheduler)
	adata = [:energy]
	adf_, mdf = run!(model, agent_step!, model_step!, 2880;
		adata=adata, mdata=mdata)
	adf = combine(groupby(adf_, :step),
		:energy => (energy -> count(x -> x > 0, energy)) => :count_bph)
	leftjoin(adf, mdf, on=:step)
end

# ╔═╡ f1bbd543-4f24-496a-804e-c594ebc22c5c
function run_and_make_video(videopath; seed, scheduler)
	model = init_model(;
		envmap=envmap(20),
		pr_killed=0.075,
		seed=seed,
		scheduler=scheduler)
	abm_video(videopath, model, agent_step!, model_step!;
		frames=2880,
		framerate=30,
		heatkwargs=heatkwargs,
		heatarray=model -> model.food,
		am=agent -> :circle
	)
	return videopath
end

# ╔═╡ dc557b3a-8f6d-4578-b723-60b8e02a02a7
let kwargs = (seed = rand(1:100), scheduler=Schedulers.randomly)
	df1=run_and_collect_data(;kwargs...)
	df2=run_and_collect_data2(;kwargs...)
	f = Figure()
	a1,a2 = Axis(f[1,1]), Axis(f[1,2])
	lines!(a1, df1.step, df1.count_bph)
	lines!(a1, df2.step, df2.count_bph)
	lines!(a2, df1.step, df1.food)
	lines!(a2, df2.step, df2.food)
	f
end

# ╔═╡ 8674a5e3-805a-4bfe-b8c1-1687a38df04e
let f = Figure()
	a = Axis(f[1,1])
	for seed in 1:20
		df = run_and_collect_data(;seed =seed, scheduler=Schedulers.randomly)
		kwargs = if df.food[end] ≤ df.food[1] ÷ 2
			lines!(a, df.step, df.count_bph; label = "Seed $seed")
		else
			lines!(a, df.step, df.count_bph;)
		end
	end
	f[1,2] = Legend(f, a, "Odd things", framevisible = false)
	f
end

# ╔═╡ 46abc4b7-6ed1-4adb-b31e-9973b06ef69d
run_and_make_video("Seed 4.mp4"; seed =4,scheduler = Schedulers.randomly)

# ╔═╡ f6fe4765-7a8f-4c48-907b-539610da783d
let f = Figure()
	a = Axis(f[1,1])
	haslegend = false
	for seed in 1:50
		df = run_and_collect_data(;seed =seed, scheduler=Schedulers.by_id)
		kwargs = if df.food[end] ≤ df.food[1] ÷ 2
			haslegend = true
			lines!(a, df.step, df.count_bph; label = "Seed $seed")
		else
			lines!(a, df.step, df.count_bph;)
		end
	end
	if haslegend
		f[1,2] = Legend(f, a, "Odd things", framevisible = false)
	end
	f
end

# ╔═╡ 18eaf72c-ec5f-4f18-9d87-e37415462643
run_and_make_video("Seed44-byid.mp4"; seed =44,scheduler = Schedulers.by_id)

# ╔═╡ 41f6b231-f8cf-46ca-a66f-ada0ce56e8b7
function run_and_make_video2(videopath; seed, scheduler)
    model = init_model(;
        envmap = envmap(20),
        pr_killed = 0.075,
        seed = seed,
        scheduler = scheduler,
    )
    abm_video(
        videopath,
        model,
        agent_step!,
        model_step!;
        frames = 2880,
        framerate = 30,
        heatkwargs = heatkwargs,
        heatarray = model -> model.food,
        am = agent -> :circle,
        scheduler = Schedulers.by_id
    )
    return videopath
end

# ╔═╡ 421492fc-72fc-446a-ba57-c2635065f7e5
run_and_make_video2("Seed4-new.mp4"; seed=4, scheduler=Schedulers.randomly)

# ╔═╡ Cell order:
# ╠═c79638a6-c80d-11eb-3254-b33f7a382f78
# ╠═7852fa86-a513-49ac-92b1-7f0a47cedc17
# ╠═6ed9031f-8754-46bf-bf92-304d1ea462b8
# ╠═55d77475-590d-4c21-a47b-2012fba13d8f
# ╠═cc080461-80c3-40cd-bc6e-d1c46f017b9a
# ╠═c7c02074-b970-425a-9586-48aa34e62c27
# ╠═a7a6cbab-5b38-42d3-9ede-1e2b5550f2f1
# ╠═79076b76-917f-439c-b954-8046867614c0
# ╠═11b14ef0-3554-498c-8e14-aafd5d026a49
# ╠═11362a40-572d-4019-9620-baecb02b01f9
# ╠═6cc30e3d-ca13-4ffb-b623-fcdb7c31936c
# ╠═fe5717d3-dd4a-470a-82bb-453e9867dc3d
# ╠═f1bbd543-4f24-496a-804e-c594ebc22c5c
# ╠═dc557b3a-8f6d-4578-b723-60b8e02a02a7
# ╠═8674a5e3-805a-4bfe-b8c1-1687a38df04e
# ╠═46abc4b7-6ed1-4adb-b31e-9973b06ef69d
# ╠═f6fe4765-7a8f-4c48-907b-539610da783d
# ╠═18eaf72c-ec5f-4f18-9d87-e37415462643
# ╠═41f6b231-f8cf-46ca-a66f-ada0ce56e8b7
# ╠═421492fc-72fc-446a-ba57-c2635065f7e5

max_stars_repo_name: ndgnuh/RiceBPH

max_issues_repo_name: ndgnuh/RiceBPH

max_forks_repo_name: ndgnuh/RiceBPH


text:
"""
This is a sample input file for the CloudScat code.  To use it run julia
and type

julia> include("wavelength.jl"); run()

The output file will be based on this file's name so to run different simulations
you can just copy this file, change the configuration and run it as above with
the new file name.  You can also run from the command line as e.g.

> julia wavelength.jl

- - -
2010 Alejandro Luque IAA-CSIC (aluque@iaa.es)

"""

using CloudScat

const co = CloudScat.constants

function run()
    for λ in [337 * co.nano, 777 * co.nano],
        R in [20 * co.micro]
        runone(λ, R)
    end
end

function runone(λ, R)
    params = Params(
        # Number of simulated photons
        # Revert to N = 10000000
        N = 10000000,

        # The wavelength of the scattered light
        λ = λ,

        # The source is defined as a straight line running from source_a
        # to source_b.  Within this line photons are distributed uniformly with
        # initially random, isotropic directions.  To use a point source set
        # source_b = source_a
        source_a = [0, 0, 10 * co.kilo],
        source_b = [0, 0, 12 * co.kilo])

    cth = 15 * co.kilo

    # World geometry
    # Describe the cloud geometry here.  It may consist in combinations of
    # geometrical figures (cylinders, spheres, cones, semi-planes).
    # See cloud_gemoetry.jl for a more detailed description
    #
    cloud = Cylinder(5 * co.kilo,  cth, 0, 0, 20 * co.kilo)

    # Geometry of the full domain.  Photons outside the domain are
    # immediately discarded.
    domain = shapediff(Plane(60 * co.kilo, false), Plane(7 * co.kilo, false))

    # Density of the scattering centers.
    nscat = 100 * co.centi^-3

    composition = Homogeneous(params.λ, nscat, R)

    # Define the full simulation world
    world = World(cloud, domain, composition)

    observers = [Observer(
        # Location of the observer
        position = [0, 0, 400 * co.kilo],

        # Sampling interval
        tsample = 1e-5,

        # Collected samples
        nsamples = 10000,

        # fov of the camera (half of the diagonal fov).
        fov = 40,

        # camera pixels
        pixels = 1024)]

    basename = splitext(splitdir(@__FILE__)[2])[1]

    # The saveto option specifies the output file.  Leave it as it is to base
    # the filename on the name of this file specifying the parameters.
    ofile = ("$(basename)_$(Int(λ/co.nano))nm_" *
             "$(Int(R/co.micro))um.h5")
    CloudScat.main(params, world, observers,
                   saveto=ofile)
end

# Check if the code has been 'included' or run from the shell.  In the latter
# case, run the simulation.
isinteractive() || run()


text:
using GraphOrdering, Test

G = Dict(
    1 => (3, 8, 9),
    2 => (3, 8, 7),
    3 => (1, 2),
    4 => (8, 9),
    5 => (7, 8),
    6 => (2, 7),
    7 => (5, 2, 6),
    8 => (1, 2, 4, 5),
    9 => (1, 4))

@test GraphOrdering.bandwidth(G) == 17


text:
f(x::Vector) = (x[1] - 5.0)^2
function g!(x::Vector, storage::Vector)
    storage[1] = 2.0 * (x[1] - 5.0)
end

initial_x = [0.0]

store_trace, show_trace = false, false
results = Optim.naive_gradient_descent(f, g!, initial_x,
	                                   0.1, 10e-8, 1_000,
	                                   store_trace, show_trace)
@assert length(results.trace.states) == 0
@assert results.converged
@assert norm(results.minimum - [5.0]) < 0.01

results = Optim.naive_gradient_descent(f, g!, initial_x)
@assert length(results.trace.states) == 0
@assert results.converged
@assert norm(results.minimum - [5.0]) < 0.01

eta = 0.9

f(x::Vector) = (1.0 / 2.0) * (x[1]^2 + eta * x[2]^2)

function g!(x::Vector, storage::Vector)
    storage[1] = x[1]
    storage[2] = eta * x[2]
end

store_trace, show_trace = true, false
results = Optim.naive_gradient_descent(f, g!, [1.0, 1.0],
		                               0.1, 10e-8, 1_000,
	                                   store_trace, show_trace)
@assert length(results.trace.states) > 0
@assert results.converged
@assert norm(results.minimum - [0.0, 0.0]) < 0.01


text:
@testset "DC Power Flow Models Monitored Line Flow Constraints and Static Unbounded" begin
    system = PSB.build_system(PSITestSystems, "c_sys5_ml")
    limits = PSY.get_flow_limits(PSY.get_component(MonitoredLine, system, "1"))
    for model in [DCPPowerModel, StandardPTDFModel]
        template = get_thermal_dispatch_template_network(model)
        op_problem_m = OperationsProblem(
            template,
            system;
            optimizer = OSQP_optimizer,
            PTDF = PSY.PTDF(system),
        )
        @test build!(op_problem_m; output_dir = mktempdir(cleanup = true)) ==
              PSI.BuildStatus.BUILT

        @test check_variable_bounded(op_problem_m, :Fp__MonitoredLine)
        @test check_variable_unbounded(op_problem_m, :Fp__Line)

        @test solve!(op_problem_m) == RunStatus.SUCCESSFUL
        @test check_flow_variable_values(
            op_problem_m,
            :Fp__MonitoredLine,
            "1",
            limits.from_to,
        )
    end
end

@testset "DC Power Flow Models Monitored Line Flow Constraints and Static with Bounds" begin
    system = PSB.build_system(PSITestSystems, "c_sys5_ml")
    set_rate!(PSY.get_component(Line, system, "2"), 1.5)
    for model in [DCPPowerModel, StandardPTDFModel]
        template = get_thermal_dispatch_template_network(model)
        set_device_model!(template, DeviceModel(Line, StaticBranch))
        set_device_model!(template, DeviceModel(MonitoredLine, StaticBranchUnbounded))
        op_problem_m = OperationsProblem(
            template,
            system;
            optimizer = OSQP_optimizer,
            PTDF = PSY.PTDF(system),
        )
        @test build!(op_problem_m; output_dir = mktempdir(cleanup = true)) ==
              PSI.BuildStatus.BUILT

        @test check_variable_unbounded(op_problem_m, :Fp__MonitoredLine)
        # Broken
        # @test check_variable_bounded(op_problem_m, :Fp__Line)

        @test solve!(op_problem_m) == RunStatus.SUCCESSFUL
        @test check_flow_variable_values(op_problem_m, :Fp__Line, "2", 1.5)
    end
end

@testset "AC Power Flow Monitored Line Flow Constraints" begin
    system = PSB.build_system(PSITestSystems, "c_sys5_ml")
    limits = PSY.get_flow_limits(PSY.get_component(MonitoredLine, system, "1"))
    template = get_thermal_dispatch_template_network(ACPPowerModel)
    op_problem_m = OperationsProblem(template, system; optimizer = ipopt_optimizer)
    @test build!(op_problem_m; output_dir = mktempdir(cleanup = true)) ==
          PSI.BuildStatus.BUILT

    @test check_variable_bounded(op_problem_m, :FpFT__MonitoredLine)
    @test check_variable_unbounded(op_problem_m, :FqFT__MonitoredLine)

    @test solve!(op_problem_m) == RunStatus.SUCCESSFUL
    @test check_flow_variable_values(
        op_problem_m,
        :FpFT__MonitoredLine,
        :FqFT__MonitoredLine,
        "1",
        0.0,
        limits.from_to,
    )
end

@testset "DC Power Flow Models for HVDCLine with with Line Flow Constraints, TapTransformer & Transformer2W Unbounded" begin
    ratelimit_constraint_names = [
        :RateLimit_ub__Transformer2W,
        :RateLimit_lb__Transformer2W,
        :RateLimit_ub__TapTransformer,
        :RateLimit_lb__TapTransformer,
    ]

    system = PSB.build_system(PSITestSystems, "c_sys14_dc")
    hvdc_line = PSY.get_component(HVDCLine, system, "DCLine3")
    limits_from = PSY.get_active_power_limits_from(hvdc_line)
    limits_to = PSY.get_active_power_limits_to(hvdc_line)
    limits_min = min(limits_from.min, limits_to.min)
    limits_max = min(limits_from.max, limits_to.max)

    tap_transformer = PSY.get_component(TapTransformer, system, "Trans3")
    rate_limit = PSY.get_rate(tap_transformer)

    transformer = PSY.get_component(Transformer2W, system, "Trans4")
    rate_limit2w = PSY.get_rate(tap_transformer)

    for model in [DCPPowerModel, StandardPTDFModel],
        hvdc_model in [HVDCDispatch, HVDCLossless]

        template = get_template_dispatch_with_network(model)
        set_device_model!(template, HVDCLine, hvdc_model)
        op_problem_m = OperationsProblem(
            template,
            system;
            optimizer = OSQP_optimizer,
            PTDF = PSY.PTDF(system),
        )
        @test build!(op_problem_m; output_dir = mktempdir(cleanup = true)) ==
              PSI.BuildStatus.BUILT

        @test check_variable_bounded(op_problem_m, :Fp__HVDCLine)
        @test check_variable_unbounded(op_problem_m, :Fp__TapTransformer)
        @test check_variable_unbounded(op_problem_m, :Fp__Transformer2W)

        psi_constraint_test(op_problem_m, ratelimit_constraint_names)

        @test solve!(op_problem_m) == RunStatus.SUCCESSFUL

        @test check_flow_variable_values(
            op_problem_m,
            :Fp__HVDCLine,
            "DCLine3",
            limits_min,
            limits_max,
        )
        @test check_flow_variable_values(
            op_problem_m,
            :Fp__TapTransformer,
            "Trans3",
            -rate_limit,
            rate_limit,
        )
        @test check_flow_variable_values(
            op_problem_m,
            :Fp__Transformer2W,
            "Trans4",
            -rate_limit2w,
            rate_limit2w,
        )
    end
end

@testset "DC Power Flow Models for Unbounded HVDCLine, and StaticBranchBounds for TapTransformer & Transformer2W" begin
    system = PSB.build_system(PSITestSystems, "c_sys14_dc")
    hvdc_line = PSY.get_component(HVDCLine, system, "DCLine3")
    limits_from = PSY.get_active_power_limits_from(hvdc_line)
    limits_to = PSY.get_active_power_limits_to(hvdc_line)
    limits_min = min(limits_from.min, limits_to.min)
    limits_max = min(limits_from.max, limits_to.max)

    tap_transformer = PSY.get_component(TapTransformer, system, "Trans3")
    rate_limit = PSY.get_rate(tap_transformer)

    transformer = PSY.get_component(Transformer2W, system, "Trans4")
    rate_limit2w = PSY.get_rate(tap_transformer)

    for model in [DCPPowerModel, StandardPTDFModel]
        template = get_template_dispatch_with_network(model)
        set_device_model!(template, DeviceModel(HVDCLine, HVDCUnbounded))
        set_device_model!(template, DeviceModel(TapTransformer, StaticBranchBounds))
        set_device_model!(template, DeviceModel(Transformer2W, StaticBranchBounds))
        op_problem_m = OperationsProblem(
            template,
            system;
            optimizer = OSQP_optimizer,
            PTDF = PSY.PTDF(system),
        )
        @test build!(op_problem_m; output_dir = mktempdir(cleanup = true)) ==
              PSI.BuildStatus.BUILT

        if model == DCPPowerModel
            # TODO: Currently Broken, remove variable bounds in HVDCUnbounded
            # @test check_variable_unbounded(op_problem_m, :Fp__HVDCLine)

            @test check_variable_bounded(op_problem_m, :Fp__TapTransformer)
            @test check_variable_bounded(op_problem_m, :Fp__TapTransformer)
        end

        @test solve!(op_problem_m) == RunStatus.SUCCESSFUL

        @test check_flow_variable_values(
            op_problem_m,
            :Fp__HVDCLine,
            "DCLine3",
            limits_min,
            limits_max,
        )
        @test check_flow_variable_values(
            op_problem_m,
            :Fp__TapTransformer,
            "Trans3",
            -rate_limit,
            rate_limit,
        )
        @test check_flow_variable_values(
            op_problem_m,
            :Fp__Transformer2W,
            "Trans4",
            -rate_limit2w,
            rate_limit2w,
        )
    end
end

@testset "AC Power Flow Models for HVDCLine Flow Constraints and TapTransformer & Transformer2W Unbounded" begin
    ratelimit_constraint_names = [
        :RateLimitFT__Transformer2W,
        :RateLimitTF__Transformer2W,
        :RateLimitFT__TapTransformer,
        :RateLimitTF__TapTransformer,
        :RateLimitFT__HVDCLine,
        :RateLimitTF__HVDCLine,
    ]

    system = PSB.build_system(PSITestSystems, "c_sys14_dc")

    hvdc_line = PSY.get_component(HVDCLine, system, "DCLine3")
    limits_from = PSY.get_active_power_limits_from(hvdc_line)
    limits_to = PSY.get_active_power_limits_to(hvdc_line)
    limits_min = min(limits_from.min, limits_to.min)
    limits_max = min(limits_from.max, limits_to.max)

    tap_transformer = PSY.get_component(TapTransformer, system, "Trans3")
    rate_limit = PSY.get_rate(tap_transformer)

    transformer = PSY.get_component(Transformer2W, system, "Trans4")
    rate_limit2w = PSY.get_rate(tap_transformer)

    template = get_template_dispatch_with_network(ACPPowerModel)
    set_device_model!(template, DeviceModel(HVDCLine, HVDCDispatch))
    op_problem_m = OperationsProblem(template, system; optimizer = ipopt_optimizer)
    @test build!(op_problem_m; output_dir = mktempdir(cleanup = true)) ==
          PSI.BuildStatus.BUILT

    check_variable_bounded(op_problem_m, :FqTF__HVDCLine)
    check_variable_bounded(op_problem_m, :Fp__HVDCLine)
    @test check_variable_bounded(op_problem_m, :FpFT__TapTransformer)
    @test check_variable_unbounded(op_problem_m, :FqFT__TapTransformer)
    @test check_variable_bounded(op_problem_m, :FpTF__Transformer2W)
    @test check_variable_unbounded(op_problem_m, :FqTF__Transformer2W)

    psi_constraint_test(op_problem_m, ratelimit_constraint_names)

    @test solve!(op_problem_m) == RunStatus.SUCCESSFUL

    @test check_flow_variable_values(
        op_problem_m,
        :Fp__HVDCLine,
        :FqTF__HVDCLine,
        "DCLine3",
        limits_max,
    )
    @test check_flow_variable_values(
        op_problem_m,
        :FpFT__TapTransformer,
        :FqFT__TapTransformer,
        "Trans3",
        rate_limit,
    )
    @test check_flow_variable_values(
        op_problem_m,
        :FpTF__Transformer2W,
        :FqTF__Transformer2W,
        "Trans4",
        rate_limit2w,
    )
end

#= This Test is Broken to do a missing implementation of add_variable
@testset "DC Power Flow Models for HVDCLine Dispatch and TapTransformer & Transformer2W Unbounded" begin
    ratelimit_constraint_names = [
        :RateLimitFT__Transformer2W,
        :RateLimitTF__Transformer2W,
        :RateLimitFT__TapTransformer,
        :RateLimitTF__TapTransformer,
        :RateLimitFT__HVDCLine,
        :RateLimitTF__HVDCLine,
    ]

    system = PSB.build_system(PSITestSystems, "c_sys14_dc")

    hvdc_line = PSY.get_component(HVDCLine, system, "DCLine3")
    limits_from = PSY.get_active_power_limits_from(hvdc_line)
    limits_to = PSY.get_active_power_limits_to(hvdc_line)
    limits_min = min(limits_from.min, limits_to.min)
    limits_max = min(limits_from.max, limits_to.max)

    tap_transformer = PSY.get_component(TapTransformer, system, "Trans3")
    rate_limit = PSY.get_rate(tap_transformer)

    transformer = PSY.get_component(Transformer2W, system, "Trans4")
    rate_limit2w = PSY.get_rate(tap_transformer)

    template = get_template_dispatch_with_network(StandardPTDFModel)
    set_device_model!(template, DeviceModel(HVDCLine, HVDCDispatch))
    op_problem_m = OperationsProblem(template, system; PTDF = PSY.PTDF(system), optimizer = ipopt_optimizer)
    @test build!(op_problem_m; output_dir = mktempdir(cleanup = true)) ==
          PSI.BuildStatus.BUILT

    check_variable_bounded(op_problem_m, :FqTF__HVDCLine)
    check_variable_bounded(op_problem_m, :FpTF__HVDCLine)
    @test check_variable_bounded(op_problem_m, :FpFT__TapTransformer)
    @test check_variable_unbounded(op_problem_m, :FqFT__TapTransformer)
    @test check_variable_bounded(op_problem_m, :FpTF__Transformer2W)
    @test check_variable_unbounded(op_problem_m, :FqTF__Transformer2W)

    psi_constraint_test(op_problem_m, ratelimit_constraint_names)

    @test solve!(op_problem_m) == RunStatus.SUCCESSFUL

    @test check_flow_variable_values(
        op_problem_m,
        :FpTF__HVDCLine,
        :FqTF__HVDCLine,
        "DCLine3",
        limits_max,
    )
    @test check_flow_variable_values(
        op_problem_m,
        :FpFT__TapTransformer,
        :FqFT__TapTransformer,
        "Trans3",
        rate_limit,
    )
    @test check_flow_variable_values(
        op_problem_m,
        :FpTF__Transformer2W,
        :FqTF__Transformer2W,
        "Trans4",
        rate_limit2w,
    )
end
=#


text:
#=------------------------------------------------------------------------------
   Core Operations associated with the COOTensor Class. This includes the
   getters/setters, addition, and contraction operations. More sophisticated
   numerical routines are saved for the NumericalRoutines.jl file.

  Operations Glossary
  -------------------
    * order - ()
    * add_edges! - (UNTESTED)
    * dense - (UNTESTED)
    * contract_edge_k_1 - ()
    * contract  - ()
    * dense_contract - ()
    * contract_k_1 - ()
    * contract_multi - ()
    * get_sub_tensor - ()
    * find_nnz - ()

------------------------------------------------------------------------------=#

export get_sub_tensor
export order
export save
export load
export flatten
export dense

#=------------------------------------------------------------------------------
						        Common Operators
------------------------------------------------------------------------------=#
function ==(A::Ten,B::Ten) where {Ten <: AbstractSSTen}
	if A.cubical_dimension != B.cubical_dimension
		return false
	end

	#TODO: add in sorting.
	for ((ind_A,v_A), (ind_B,v_B)) in zip(A,B)
	  if v_A != v_B || ind_A != ind_B
		  return false
	  end
	end
	true
end

function ==(A::SSSTensor,B::SSSTensor)
	if A.cubical_dimension != B.cubical_dimension
		return false
	end

	for (ind_A,v_A) in A
	  if get(B.edges,ind_A,nothing) != v_A
		  return false
	  end
	end
	true
end

#=------------------------------------------------------------------------------
						           File I/O
------------------------------------------------------------------------------=#
"""-----------------------------------------------------------------------------
    save(A, filepath)

  Saves the the file in the desired location as a ssten file. The ssten file
  format is a generalization of the smat format to the sparse symmetric tensor
  case.  The first line of the file is tab separated line indicated the order of
  the tensor, the cubical dimension, and the number of unique non-zeros in the
  symmetric tensor all stored in that order, as ints:

      (order::Int)\\t(cubical_dimension::Int)\\t(non-zero count)\\n  .

  The subsequent lines are tab separated lines which store the unique non-zeros
  which occur in the tensor. Each line contains the indices of the hyperedge in
  sorted order as ints, and the weight as a floating point number:

      '(v_{i_1}::Int)\\t...(v_{i_k}::Int)\\t(edgeweight::Float)\\n'  .

  If the filepath doesn't end in ".ssten", it will be appended to the filepath.

  Inputs:
  -------
  * A - (SSSTensor):
    Instance of the sparse symmetric tensor we wish to save.

  * filepath - (string):
    The location to save the file.
----------------------------------------------------------------------------"""
function save(A::Ten, filepath::String) where {Ten <: AbstractSSTen}
    file = open(filepath, "w")

    header = "$(order(A))\t$(A.cubical_dimension)\t$(length(A))\n"
    write(file,header);

    for (edge,weight) in A
		edge_string=""
        for v_i in edge
		   edge_string *= string(v_i,"\t")
        end
		edge_string *= string(weight,"\n")
		write(file,edge_string)
    end

	close(file)
end

#TODO: add in a loading bar
"""-----------------------------------------------------------------------------
    load(filepath, enforceFormatting, type, NoChecks)

  Loads in a sparse symmetric tensor from a .ssten file. See save for expected
  .ssten file formatting specifications. The routine will check to see if the
  vertex indices are 0 indexed or not, and will increment each by 1 if it is to
  adhere to Julia's use of indexing by 1.

  Inputs:
  -------
  * filepath - (String):
    The .ssten file to load from. Throws an error if not .ssten formatting.

  * enforceFormatting - (Optional Bool):
    Sorts each of the hyperedge indices before storing them. Can be used as a
    quick fix for a improperly formatted file. Not recommended as it will add
    to memory usage and run time.

  * type - (optional String):
    The subtype of AbstractSSTen to return

  * nochecks - (optional Bool):
    stops constructors from checking the validity of the input, only use if
    absolutely sure the file is properly formatted for the tensor (with the
    exception of 0 indexing).

  Outputs:
  --------
  * A - (SSSTensor):
    The Sparse symmetric tensor stored in the file.

-----------------------------------------------------------------------------"""
function load(filepath::String;enforceFormatting::Bool=false,
              type::String="DICTen",nochecks::Bool=false)
	#check path validity
	@assert filepath[end-5:end] == ".ssten"

    open(filepath) do file
		#preallocate from the header line
		order, n, m =
			[parse(Int,elem) for elem in split(chomp(readline(file)),'\t')]

		indices = Array{Int,2}(undef,m,order)
		values = Array{Float64,1}(undef,m)

		#hyperedges = Array{Tuple{Array{Int64,1},Float64},1}(undef, m)

		i = 1
		@inbounds for line in eachline(file)
			entries = split(chomp(line),'\t')
			indices[i,:] = [parse(Int,elem) for elem in entries[1:end-1]]
			if enforceFormatting
				sort!(indices[i,:])
			end
			values[i] = parse(Float64,entries[end])
			i += 1
		end

		#check for 0 indexing
		zero_indexed = false

		@inbounds for i in 1:m
		    if indices[i,1] == 0
    			zero_indexed = true
				break
			end
	    end

		if zero_indexed
			indices .+= 1
		end

		if type == "DICTen"
			if nochecks
				warn("nocheck unimplemented in DictTen currently")
			end
			return SSSTensor([(indices[i,:],values[i]) for i in 1:m],n)
		elseif type == "COOTen"
			return COOTen(indices,values,n;nocheck=nochecks)
		else
			error("type:$(type) not defined.\n Curently supporting:'DICTen','COOTen'.")
		end
	end
end

#=------------------------------------------------------------------------------
						           Getters
------------------------------------------------------------------------------=#

"""-----------------------------------------------------------------------------
 	  This function returns the order of the tensor passed in.

    Input:
    ------
    * A -(SSSTensor):

      The tensor in question.

    Output:
    -------
    * order - (Int)

      The order of the tensor.
-----------------------------------------------------------------------------"""
function order(A::SSSTensor)

    for (indices,_) in A.edges
        return length(indices)
    end
end

function order(A::COOTen)
    return A.order
end

"""-----------------------------------------------------------------------------
    find_nnz(A)

  Finds the non-zeros in a k-dimensional array and returns the list of the
indices associated along with a count of the non-zeros.
-----------------------------------------------------------------------------"""
@generated function find_nnz(A::Array{N,k}) where {N<:AbstractFloat,k}
    quote
        n = size(A)[1]
        y = Array{Tuple{Array{Int,1},N},1}(undef,n^k)
        nnz = 1

        @nloops $k i A begin
            val = @nref $k A i
            if val != 0.0
                y[nnz] = (collect(@ntuple $k i),val)
                nnz += 1
            end
        end
        nnz -= 1 #take care of overcount
        return y[1:nnz], nnz
    end
end

"""-----------------------------------------------------------------------------
    get_sub_tensor(A,indices,remap=false)

  This function produced a subtensor which only containes the vertices specified
in the array indices. Note that there must be hyper edges shared between the
vertices in indices, otherwise the SSSTensor constructor will throw an error
when given an empty list of hyperedges. Does not run the input checks on the
subtensor created.

Input:
------
* A - (SSSTensor)

    The tensor to produce a subtensor from.
* indices - (Array{Int,1} or Set{Int})

    The indices to build a subtensor from.
* remap - (optional bool)

    Indicates whether or not to remap the vertices from 1 to length(indices).
  Useful if the desired subtensor is to be used as a standalone tensor.
Output:
-------
* sub_tensor - (SSSTensor)

    The subtensor of A which only contains the hyperedges of A which all include
    each index in indices.
-----------------------------------------------------------------------------"""
function get_sub_tensor(A::Ten,indices::T,
                        remap::Bool=false) where {T <: Union{Array{Int,1},Set{Int}},
												  Ten <: AbstractSSTen}
	@assert 0 < length(indices) <= A.cubical_dimension
	@assert all(indices .> 0)

	if T == Array{Int,1}
		indices = Set(indices) #replace with a bitmap
	end

	incident_edges = find_edge_incidence(A)
	sub_tensor_edges = Array{Tuple{Array{Int,1},AbstractFloat}}(undef,0)

	for v_i in indices
		for (V,val) in get(incident_edges,v_i,[])
			edge = (V,val)

			#only include edge if all indices in hyperedge are desired
			if all(x -> x in indices,V)
				push!(sub_tensor_edges,edge)
			end

			#remove all other edges in incidence
			for v_j in V
				delete!(get(incident_edges,v_j,[]),edge)
			end
		end
	end

	if remap
		remap_indices!(sub_tensor_edges)
	end

	if Ten == COOTen

		nnz = length(sub_tensor_edges)
		k = length(sub_tensor_edges[1][1])

		sub_indices = Array{Int,2}(undef,nnz,k)
		sub_vals = Array{typeof(sub_tensor_edges[1][2]),1}(undef,nnz)

		for (i,(inds,val)) in zip(1:nnz,sub_tensor_edges)
			sub_indices[i,:] = inds
			sub_vals[i] = val
		end
		p = sort(1:nnz, by=i->sub_indices[i,:])

		COOTen(sub_indices[p,:],sub_vals[p];nocheck=true)#doesn't check input
	elseif Ten == SSSTensor
		SSSTensor(sub_tensor_edges,nocheck=true)#doesn't check input
	end
end

#=------------------------------------------------------------------------------
						             Setters
------------------------------------------------------------------------------=#

#UNTESTED
"""-----------------------------------------------------------------------------
    add_edges!(A,edges)

  This function takes in a list of edges and adds them into the SSSTensor. If an
edge is already present in the tensor, and the value is added in at that index.

Input:
------
* A - (SSSTensor)

    The tensor to add hyper edges to.
* edges - (Array{Tuple{Array{Int,1},Float},1})

    An array of pairs which contain the indices in the first element, and the
    value in the second element. Note each edge's indices must be in range.
-----------------------------------------------------------------------------"""
function add_edges!(A::SSSTensor,edges::Array{Tuple{Array{Int,1},N},1}) where N <: AbstractFloat
    #check edges' validity
    for (indices,_) in edges
        sort!(indices)
        @assert indices[1] > 0
        @assert indices[end] <= A.cubical_dimension
    end

    for (indices,v) in edges
        if haskey(A.edges, indices)
            A.edges[indices] += v
        else
            A.edges[indices] = v
        end
    end
end

#=------------------------------------------------------------------------------
						   Data Representation Converters
------------------------------------------------------------------------------=#

#UNTESTED
"""-----------------------------------------------------------------------------
    dense(A)

  This function returns a dense representation of the SSSTensor passed in.

Inputs:
-------
* A - (SSSTensor)

    The sparse tensor to be converted.

Output:
-------
* B - (Array{AbstractFloat,k})

    The corresponding dense tensor representation of A.
-----------------------------------------------------------------------------"""
function dense(A::SSSTensor)
    B = zeros(eltype(A.edges.vals),Tuple(repeat([A.cubical_dimension],order(A))))

    for (indices,_) in A.edges
        for p in permutations(indices)
            B[CartesianIndex(Tuple(p))] = A.edges[indices]
        end
    end

    return B
end

#UNTESTED
"""-----------------------------------------------------------------------------
    flatten(A,colIndex)

  This function returns a sparse matrix from flattening k-2 modes into the
  columns of the matrix.  Defaults to col major formatting.

  TODO: add colIndex function to parameters.

Inputs:
-------
* A - (SSSTensor)

    The sparse tensor to be converted.
*
Output:
-------
* B - (Array{AbstractFloat,k})

    The corresponding dense tensor representation of A.
-----------------------------------------------------------------------------"""
function flatten(A::Ten) where {Ten <: AbstractSSTen}

	ord = order(A)
	nnz = length(A.edges)*factorial(ord)
	I = Array{Int,1}(undef,nnz)
	J = Array{Int,1}(undef,nnz)
	V = Array{Float64,1}(undef,nnz)

	i = 0
	for (edge,val) in A

		for indices in unique(permutations(edge))
			i+=1
			I[i] = indices[1]
			V[i] = val

			#adjust back to 0 indexing
			J[i] = foldl((x,y)-> x+A.cubical_dimension*(y-1),indices[2:end])

		end
	end

	return SparseArrays.sparse(I[1:i],J[1:i],V[1:i],A.cubical_dimension,A.cubical_dimension^(ord-1))
end


text:
mutable struct Node{Data}
    data::Data
    children::Vector{Node{Data}}
    foldchildren::Bool
    lastbranch::Bool   # used internally for working storage during `iterate`
    parent::Node{Data}

    """
        Node(data, foldchildren::Bool=false)

    Create the root of a folded tree. `data` holds the data associated with the root node,
    and `foldchildren` specifies whether its children are initially folded.
    """
    function Node{Data}(data, foldchildren::Bool=false) where Data
        new(data, Node[], foldchildren, false)
    end

    """
        Node(data, parent::Node, foldchildren::Bool=false)

    Add a new child node to `parent`. `data` holds the data associated with the child node,
    and `foldchildren` specifies whether its children are initially folded.
    """
    function Node{Data}(data, parent::Node{Data}, foldchildren::Bool=false) where Data
        child = new(data, Node[], foldchildren, false, parent)
        push!(parent.children, child)
        return child
    end
end

Node(data, foldchildren::Bool=false) = Node{typeof(data)}(data, foldchildren)
Node(data, parent::Node{Data}, foldchildren::Bool=false) where Data = Node{Data}(data, parent, foldchildren)

Base.eltype(::Type{Node{Data}}) where Data = Data
Base.IteratorSize(::Type{<:Node}) = Base.SizeUnknown()

"""
    isroot(node)

Return `true` if `node` is the root node (meaning, it has no parent).
"""
isroot(node::Node) = !isdefined(node, :parent)

"""
    toggle!(node)

Change the folding state of `node`. See also [`fold!`](@ref) and [`unfold!`](@ref).
"""
toggle!(node::Node) = (node.foldchildren = !node.foldchildren)

"""
    fold!(node)

Fold the children of `node`. See also [`unfold!`](@ref) and [`toggle!`](@ref).
"""
fold!(node::Node) = node.foldchildren = true

"""
    unfold!(node)

Unfold the children of `node`. See also [`fold!`](@ref) and [`toggle!`](@ref).
"""
unfold!(node::Node) = node.foldchildren = false

"""
    count_open_leaves(node)

Return the number of unfolded descendants of `node`.
"""
function count_open_leaves(node::Node, count::Int=0)
    count += 1    # count self
    if !node.foldchildren
        # count children (and their children...)
        for child in node.children
            count = count_open_leaves(child, count)
        end
    end
    return count
end

"""
    newnode, newdepth = next(node, depth::Int=0)

Return the next node in a depth-first search.
`depth` counts the number of levels below the root.
The parent is visited before any children.
"""
function next(node, depth::Int=0)
    if !node.foldchildren && !isempty(node.children)
        return first(node.children), depth+1
    end
    return upnext(node, depth)
end
function upnext(node, depth)
    # node.parent must be defined if we're calling this
    p = node.parent
    myidx = findfirst((==)(node), p.children)
    if myidx < length(p.children)
        return p.children[myidx+1], depth
    end
    return upnext(p, depth-1)
end

"""
    newnode, newdepth = prev(node, depth::Int=0)

Return the previous node in a depth-first search.
`depth` counts the number of levels below the root.

This traverses in the opposite direction as [`next`](@ref),
so last, deepest children are visited before their parents.
"""
function prev(node, depth::Int=0)
    p = node.parent
    myidx = findfirst((==)(node), p.children)
    if myidx > 1
        return lastchild(p.children[myidx-1], depth)
    end
    return p, depth-1
end
function lastchild(node, depth)
    if node.foldchildren || isempty(node.children)
        return node, depth
    end
    return lastchild(node.children[end], depth+1)
end

# During iteration, we mark each node as to whether it's on the terminal branch of all ancestors.
function Base.iterate(root::Node)
    root.lastbranch = true   # root has no ancestors so it is last by definition
    return root.data, (root, 0)
end

function Base.iterate(root::Node, state)
    node, depth = state
    lb = node.lastbranch
    thislb = node.foldchildren | isempty(node.children)
    (lb & thislb) && return nothing
    # We can use `next` safely because parents are visited before children,
    # so `lastbranch` is guaranteed to be set properly.
    node, depth = next(node, depth)
    p = node.parent
    node.lastbranch = p.lastbranch && node == p.children[end]
    return node.data, (node, depth)
end

struct NodeWrapper{N}
    node::N
end

function Base.iterate(rootw::NodeWrapper)
    root = rootw.node
    root.lastbranch = true   # root has no ancestors so it is last by definition
    return root, (root, 0)
end

function Base.iterate(rootw::NodeWrapper, state)
    node, depth = state
    lb = node.lastbranch
    thislb = node.foldchildren | isempty(node.children)
    (lb & thislb) && return nothing
    # We can use `next` safely because parents are visited before children,
    # so `lastbranch` is guaranteed to be set properly.
    node, depth = next(node, depth)
    p = node.parent
    node.lastbranch = p.lastbranch && node == p.children[end]
    return node, (node, depth)
end

"""
    itr = nodes(node)

Create an iterator `itr` that will return the nodes, rather than the node data.

# Example

```julia
julia> foreach(unfold!, nodes(root))
```

would ensure that each node in the tree is unfolded.
"""
nodes(node::Node) = NodeWrapper(node)


text:
using DrWatson
@quickactivate "ConstrainedFWIExamples"
DrWatson.greet()

max_stars_repo_name: slimgroup/ConstrainedFWIExamples

max_issues_repo_name: slimgroup/ConstrainedFWIExamples

max_forks_repo_name: slimgroup/ConstrainedFWIExamples


text:
__precompile__()

module ImageMorphology

using ImageCore

include("dilation_and_erosion.jl")
include("thinning.jl")

export
    dilate,
    erode,
    opening,
    closing,
    tophat,
    bothat,
    morphogradient,
    morpholaplace,

    thinning,
    GuoAlgo

end # module


text:
"""
m: minibatch size
N: full data size
"""
function subsample_index(m::Integer, N::Integer)
  @assert m < N
  shuffled_idx = Random.shuffle(1:N)
  return shuffled_idx[1:m], shuffled_idx[(m+1):end]
end

function sample_minibatch(batchsizes::Vector{S},
                          y::Vector{Matrix{T}}) where {T <: AbstractFloat,
                                                       S <: Integer}
  I = length(y)
  @assert length(batchsizes) == I
  N = size.(y, 1)
  J = size(y[1], 2)

  idx_mini, idx_mega = zip([subsample_index(batchsizes[i], N[i]) for i in 1:I]...)
  ymini = [y[i][idx_mini[i], :] for i in 1:I]
  ymega = [y[i][idx_mega[i], :] for i in 1:I]

  return idx_mini, idx_mega, Data(ymini), Data(ymega)
end


text:
using Images
using ImageView
using LinearAlgebra
using Plots

n = 1000 # 1000 x 1000 square
X = zeros(n, n)
X[250:750, 250:750] .= 1
imshow(X)

Xrot = imrotate(X, -pi/18, axes(X))
Xrot[isnan.(Xrot)] .= 0
imshow(Xrot)

U, S, V = svd(X) # SVD well-aligned square
U2, S2, V2 = svd(Xrot) # SVD rotated square
p1 = plot(S, yaxis = :log)
p2 = plot(S2, yaxis = :log)


text:
using Test
using BenchmarkTools
using polyBindingModel
using ForwardDiff
using LinearAlgebra

"""
multinomial(k...)
Multinomial coefficient where `n = sum(k)`.
"""
function multinomial(k...)
    s = 0
    result = 1
    @inbounds for i in k
        s += i
        result *= binomial(s, i)
    end
    result
end

function vec_comb(length, rsum, resid)
    if length <= 1
        return fill(rsum, (1, 1))
    end
    enum = Array{Int64}(undef, length, 0)
    for i = max(0, rsum - sum(resid[2:end])):min(resid[1], rsum)
        other_ele = vec_comb(length - 1, rsum - i, resid[2:end])
        first_ele = ones(Int64, 1, size(other_ele, 2)) * i
        enum = hcat(enum, vcat(first_ele, other_ele))
    end
    return enum
end

function polyfc_via_polyc(Req, L0::Real, KxStar::Real, f::Number, Rtot::Vector, LigC::Vector, Kav::Matrix)
    LigC /= sum(LigC)
    Cplx = vec_comb(length(LigC), f, repeat([f], length(LigC)))'
    Ctheta = exp.(Cplx * log.(LigC)) .* [multinomial(x...) for x in eachrow(Cplx)]
    @assert sum(Ctheta) ≈ 1.0 "Ctheta is $(Ctheta) with sum $(sum(Ctheta)) != 1.0"

    return polyc_Req(Req, L0, KxStar, Rtot, Cplx, Ctheta, Kav)
end

L0 = 1.0e-8
KxStar = 1.0e-10
Rtot = [100.0, 1000.0, 10.0, 1000000.0]
Kav = ones((3, 4)) * 1.0e6
Kav[1, 2] = 1.0e9

@testset "randomBinding tests" begin
    IgGC = [1.0, 0.1, 3.0]

    @testset "Can successfully assemble the parameters and get a sane result." begin
        @btime polyfc($L0, $KxStar, 4, $Rtot, $IgGC, $Kav)

        for ff in 1:40
            out = polyfc(L0, KxStar, ff, Rtot, IgGC, Kav)
            # Mass balance
            @test all(out.Rbound_n .<= Rtot)
            @test all(out.Rbound_n .>= 0.0)
            @test all(out.Req .<= Rtot)
            @test all(isapprox.(out.Rbound_n .+ out.Req, Rtot))

            nres = polyfc_via_polyc(out.Req, L0, KxStar, ff, Rtot, IgGC, Kav)
            @test out.Lbound ≈ sum(nres[1])
            @test isapprox(out.Rbound, sum(nres[2]), rtol = 1e-6)
        end
    end

    @testset "Test monovalent case." begin
        out = polyfc(L0, KxStar, 1, Rtot, [1], Kav[1, :]')

        # Note f is not used
        comp = vec(Kav[1, :]' .* L0 .* Rtot' ./ (1 .+ (Kav[1, :]' .* L0)))

        @test all(out.Lbound .≈ sum(comp))
        @test all(out.Rbound_n .≈ comp)
        @test all(out.Rbound_n .+ out.Req .≈ Rtot)
        @test out.Rmulti == 0.0
        @test all(out.Rmulti_n .== 0.0)
    end
end

@testset "Test Lfbnd in complexBinding when f = 1, 2, 3" begin
    # f = 1
    Cplx = [1 0 0; 0 1 0; 0 0 1]
    θ = rand(3)
    Lbounds, Rbounds, Lfbnds, Lmbnds = polyc(L0, KxStar, Rtot, Cplx, θ / sum(θ), Kav)
    @test all(Lbounds .≈ Lfbnds)
    @test all(Lbounds .≈ sum(Rbounds, dims = 2))
    @test all(Lbounds .≈ Lmbnds)

    # f = 3
    Cplx = [3 0 0; 2 1 0; 2 0 1; 1 2 0; 1 1 1; 1 0 2; 0 3 0; 0 2 1; 0 1 2; 0 0 3]
    θ = rand(10)
    Lbounds, Rbounds, Lfbnds, Lmbnds = polyc(L0, KxStar, Rtot, Cplx, θ / sum(θ), Kav)
    Lbbnds = Lbounds - Lfbnds - Lmbnds
    @test all(sum(Rbounds, dims = 2) .≈ Lfbnds * 3 + Lbbnds * 2 + Lmbnds)
end

include("testForwardDiff.jl")


text:
const SGA = TransferOperatorGenerator{SingleGrid{RectangularBinning}}
@recipe function f(iv::InvariantMeasure{TransferOperatorApproximation{<:SGA}})
    heatmap(iv.to.M)
end


text:
# separable.jl : separable Gaussian random fields

## SeparableCovarianceFunction ##
struct SeparableCovarianceFunction{d,V} <: AbstractCovarianceFunction{d}
    cov::V
end

"""
    SeparableCovarianceFunction(cov...)

Create a separable covariance function in `length(cov)` dimensions for the covariance structures `cov`. Usefull for defining anisotropic covariance functions, or if the usual KL expansion is too expensive.

# Examples
```
julia> e = Exponential(0.1)
exponential (λ=0.1, σ=1.0, p=2.0)

julia> m = Matern(0.01,1.0)
Matérn (λ=0.01, ν=1.0, σ=1.0, p=2.0)

julia> c = SeparableCovarianceFunction([e,m])
2d separable covariance function [ exponential (λ=0.1, σ=1.0, p=2.0), Matérn (λ=0.01, ν=1.0, σ=1.0, p=2.0) ]

```
See also: [`CovarianceFunction`](@ref), [`KarhunenLoeve`](@ref)
"""
SeparableCovarianceFunction(cov::Vector{<:CovarianceStructure}) = SeparableCovarianceFunction{length(cov),typeof(cov)}(cov)
SeparableCovarianceFunction(cov::CovarianceStructure...) = SeparableCovarianceFunction([cov...])

Base.eltype(cov::SeparableCovarianceFunction) = mapreduce(eltype, promote_type, cov.cov)

function GaussianRandomField(mean::Array{<:Real}, cov::SeparableCovarianceFunction{d}, kl::KarhunenLoeve{n}, pts::Vararg{AbstractVector,d}; kwargs...) where {d,n}
    size(mean) == length.(pts) || throw(DimensionMismatch("size of the mean does not correspond to the dimension of the points"))

    # generate the 1d grf's
    data = ntuple(length(cov.cov)) do i
        cov_ = CovarianceFunction(1, cov.cov[i])
        if isa(cov.cov[i], Exponential) && isone(cov.cov[i].p)
            compute_analytic(cov_,n,pts[i])
        else
            _GaussianRandomField(mean,cov_,kl,pts[i];kwargs...).data
        end
    end

    # determine n-d eigenvalues
    m = map(prod, Iterators.product((data[i].eigenval for i in 1:d)...))
    idx = sortperm(vec(m), rev=true)
    pidx = Iterators.product((1:length(data[i].eigenval) for i in 1:d)...)

    alldata = (collect(pidx)[idx], data)

    GaussianRandomField{typeof(kl),typeof(cov),typeof(pts),typeof(mean),typeof(alldata)}(mean,cov,pts,alldata)
end

# zero-mean GRF
GaussianRandomField(cov::SeparableCovarianceFunction{d}, kl, pts::Vararg{AbstractVector,d}; kwargs...) where d =
GaussianRandomField(zeros(eltype(cov), length.(pts)), cov, kl, pts...; kwargs...)

# constant mean GRF
GaussianRandomField(mean::Real, cov::SeparableCovarianceFunction{d}, kl, pts::Vararg{AbstractVector,d}; kwargs...) where d =
GaussianRandomField(fill(convert(eltype(cov), mean), length.(pts)), cov, kl, pts...; kwargs...)

# sample function
function sample(grf::GaussianRandomField{KarhunenLoeve{n},<:SeparableCovarianceFunction}; xi::Vector{<:Real} = randn(n)) where n
    length(xi) == n || throw(DimensionMismatch("length of random points vector must be equal to $(n)"))
    order, data = grf.data
    x = zeros(eltype(xi), length(grf.mean))
    d = length(grf.cov.cov)
    for i in 1:n
        orderi = order[i]
        ev = prod(data[j].eigenval[orderi[j]] for j = 1:d)
        ef = kron((data[j].eigenfunc[:,orderi[j]] for j = 1:d)...)
        x .+= (xi[i] * ev) .* ef
    end
    grf.mean + prod(grf.cov.cov[i].σ for i = 1:d) * reshape(x, size(grf.mean))
end

function sample(grf::GaussianRandomField{KarhunenLoeve{n},<:SeparableCovarianceFunction{1}}; xi::Vector{<:Real} = randn(n)) where n
    length(xi) == n || throw(DimensionMismatch("length of random points vector must be equal to $(n)"))
    order, data = grf.data
    x = data[1].eigenfunc * (data[1].eigenval .* xi)
    grf.mean + grf.cov.cov[1].σ*x
end

function Base.show(io::IO, s::SeparableCovarianceFunction)
    print(io, ndims(s), "d separable covariance function [ ", s.cov[1])
    for i in 2:length(s.cov)
        print(io, ", ", s.cov[i])
    end
    print(io, " ]")
end


text:
#=
Plasticity between 2 "virtual" neurons
The neurons are input neurons, i.e. their spikes are
entirely driven and specified externally
=#
push!(LOAD_PATH, abspath(@__DIR__,".."))

using Test
using LinearAlgebra,Statistics,StatsBase,Distributions
using Plots,NamedColors ; theme(:dark)
using SparseArrays
using SpikingRNNs; const global S = SpikingRNNs
using BenchmarkTools
function oneDSparse(w::Real)
  mat=Matrix{Float64}(undef,1,1) ; mat[1,1]=w
  return sparse(mat)
end

##
dt = 1E-3
Ttot = 10.0
# start simple: two Poisson neurons
# 1 is connected to 2 and has some plasticity
ps1 = S.PSPoisson(12.,1E4,1)
ps2 = S.PSPoisson(20.,-1E4,1)

myplasticity = let τplus = 0.2,
  τminus = 0.2
  τx = 0.3
  τy = 0.3
  A2plus = 0.01
  A3plus = 0.01
  A2minus = 0.02
  A3minus = 0.02
  (n_post,n_pre) = (2,2)
  S.PlasticityTriplets(τplus,τminus,τx,τy,A2plus,A3plus,
    A2minus,A3minus,n_post,n_pre)
end

conn_2_1 = let w = oneDSparse(10.)
  S.ConnectionPlasticityTest(w,myplasticity)
end

# neuron 1 receives no connections
pop1 = S.UnconnectedPopulation(ps1)
# neuron 2 is connected to 1
pop2 = S.Population(ps2,(conn_2_1,ps1))

# that's it, let's make the network
myntw = S.RecurrentNetwork(dt,pop1,pop2)

## recorders : spikes
krec=1
rec_spikes1 = S.RecSpikes(ps1,50.0,Ttot)
rec_spikes2 =  S.RecSpikes(ps2,50.0,Ttot)
# all weights (but there is only one)
rec_weights = S.RecWeights(conn_2_1,krec,dt,Ttot)

## Run

times = (0:myntw.dt:Ttot)
nt = length(times)
# clean up
S.reset!.([rec_spikes1,rec_spikes2,rec_weights])
S.reset!.([ps1,ps2])
S.reset!.(conn_2_1.plasticities)
# initial conditions

for (k,t) in enumerate(times)
  rec_spikes1(t,k,myntw)
  rec_spikes2(t,k,myntw)
  rec_weights(t,k,myntw)
  S.dynamics_step!(t,myntw)
end

##

_ = let plt=plot(leg=false)
  scatter!(plt,_->1,rec_spikes1.spiketimes,marker=:vline,color=:white,ylims=(0.5,3))
  scatter!(plt,_->2,rec_spikes2.spiketimes,marker=:vline,color=:white,ylims=(0.5,3))
end

##

_ = let x=rec_weights.times
 y = rec_weights.weights_now[:]
 spkdict = S.get_spiketimes_dictionary(rec_spikes1)
 plt=plot(leg=false)
 plot!(plt,x,y;linewidth=2,leg=false)
 for (neu,spkt) in pairs(spkdict)
  scatter!(twinx(plt),_->neu,spkt ;
    marker=:vline,color=:green,ylims=(0.5,3))
 end
 plt
end


text:
using PkgSRA
using Plots, Random, OrdinaryDiffEq, Statistics
using Distributions, Turing
using StatsPlots, DSP
pyplot()
Random.seed!(11)

#####
##### Define the controlled dynamical system
#####
# Load example problem
include("paper_settings.jl")
include(EXAMPLE_FOLDERNAME*"example_lorenz.jl")

# Define the forcing functions
num_ctr = 5;
    U_starts = rand(num_ctr,1) .* tspan[2]
    U_width = 0.05;
    my_U_func_time(t, u) = U_func_time(t, u, U_width, U_starts, F_dim=3)

#####
##### Produce data
#####
sol = solve_lorenz_system(my_U_func_time)
dat_clean = Array(sol)

# NEW: ADD NOISE
rng = MersenneTwister(128);
noise = 0.1
dat = dat_clean .+ noise.*randn(rng, eltype(dat_clean), size(dat_clean))

plot3d(dat[1, :], dat[2, :], dat[3, :])
    title!("Noisy Data")

# Derivatives
numerical_grad = numerical_derivative(dat, ts)
true_grad = zeros(size(dat))
for i in 1:size(dat,2)
    true_grad[:,i] = lorenz_system(true_grad[:,i], dat[:,i], p, [0])
end
# plot(numerical_grad[3,:], label="numerical")
#     plot!(true_grad[3,:], label="true uncontrolled")

# True control signal
U_true = zeros(size(dat))
for (i, t) in enumerate(ts)
    U_true[:,i] = my_U_func_time(t, dat[:,i])
end

# plot(U_true')

#####
##### Define the Bayesian modeling framework
#####

# Define the model
function lorenz_system(u, p, t)
    u = convert.(eltype(p),u)

    ρ, σ, β = p
    x, y, z = u
    du = [σ*(y-x); x*(ρ-z) - y; x*y - β*z]
    return du
end

@model lorenz_grad_residual(y) = begin
    # Lorenz parameters
    ρ ~ Normal(p[1], 1.0)
    σ ~ Normal(p[2], 0.1)
    β ~ Normal(p[3], 0.1)
    params = [ρ, σ, β]

    t = [0] # Not time dependent
    for i in 1:size(y,2)
        u = dat[:,i]
        du_params = lorenz_system(u, params, t)
        # TODO: Hard-coded noise
        y[:, i] ~ MvNormal(du_params, [10.0, 10.0, 10.0])
    end
end;

# Settings of the Hamiltonian Monte Carlo (HMC) sampler.
iterations = 1000
    ϵ = 0.02
    τ = 10
# Try to predict the GRADIENT from data
num_training_pts = 500
    train_ind = 201:num_training_pts+200
    y = numerical_grad[:,train_ind]
chain = sample(lorenz_grad_residual(y),
                HMC(iterations, ϵ, τ));
# plot3d(dat[1, train_ind], dat[2, train_ind], dat[3, train_ind])
#     title!("Training data")
# plot(chain)

# Generate test trajectories from the posterior
num_samples = 100
    param_samples = sample(chain, num_samples)
    save_ind = 3
    t = [0]
    num_test_pts = 500
    vars = [:ρ, :σ, :β]
    all_vals = zeros(num_samples, num_test_pts)
    plot_ind = 201:num_test_pts+200
for i in 1:num_samples
    these_params = [param_samples[v].value[i] for v in vars]
    for (i_save, i_dat) in enumerate(plot_ind)
        all_vals[i, i_save] = lorenz_system(dat[:,i_dat],
                                these_params, t)[save_ind]
    end
end

# Calculate the residuals
this_std = std(all_vals, dims=1)
    this_mean = mean(all_vals, dims=1)

# Align the signals
dat_grad = numerical_grad[save_ind, plot_ind]

# Function for processing the residual into a true guess
function process_residual(mean, std)
    low = mean .- std
    high = mean .+ std
    real_ind = .!((low .< 0.0) .& (high .> 0.0))

    ctr = zeros(size(mean))
    ctr[real_ind] = mean[real_ind]

    return ctr
end

#####
##### Plot everything
#####

this_plot_ind = 1:300

p1 = plot(vec(dat_grad[this_plot_ind]),
            linecolor=COLOR_DICT["true"], label="Data", lw=3);
    plot!(vec(this_mean[this_plot_ind]),
            linecolor=COLOR_DICT["intrinsic"], ribbon=vec(this_std),
            fillalpha=0.5, label="Models", lw=1);
    title!("Reconstruction");

residual = vec(dat_grad[this_plot_ind]) .- vec(this_mean[this_plot_ind])
p2 = plot(residual, lw=2, ribbon=vec(this_std), fillalpha=0.5,
            linecolor=COLOR_DICT["intrinsic"], legend=false);
    title!("Residual");

ctr_guess = process_residual(residual, this_std[this_plot_ind])
p3 = plot(ctr_guess, lw=2, label="Processed Residual",
            linecolor=COLOR_DICT["control_time"]);
    plot!(vec(U_true[3,plot_ind][this_plot_ind]), label="True control signal",
                linecolor=COLOR_DICT["control_true"]);
    title!("Learned signal");

## Start: Test plot
# residual_ratio = abs.(residual) ./ vec(this_std[this_plot_ind])
# plot(residual_ratio, label="Residual")
#     plot!(vec(U_true[3,plot_ind][this_plot_ind]), label="True control signal");
#     title!("??")
## End: Test plot

# All together
l = @layout [p1; p2; p3]
    p_all = plot(p1, p2, p3, layout=l)

fname = FIGURE_FOLDERNAME * "fig_bayesian_sra2.png";
savefig(p_all, fname)


text:
mutable struct Int64Particle
  x::Int64
  Int64Particle() = new()
end

struct RPIO{Particle}
  N::Int64
  n::Int64
  threads::Int64
  potentials::Vector{Float64}
  particles::Vector{Particle}
  ks::Vector{Int64}
  scratch::Vector{Float64}
  scratchInt::Vector{Int64}
  scratchInt64Particle::Vector{Int64Particle}
end

function RPIO{Particle}(N::Int64, n::Int64, threads::Int64) where Particle
  potentials::Vector{Float64} = Vector{Float64}(undef, N)
  particles::Vector{Particle} = Vector{Particle}(undef, N)
  for i in 1:N
    particles[i] = Particle()
  end
  ks::Vector{Int64} = Vector{Int64}(undef, n-1)
  scratch::Vector{Float64} = Vector{Float64}(undef, threads)
  scratchInt::Vector{Int64} = Vector{Int64}(undef, threads)
  scratchInt64Particle::Vector{Int64Particle} =
    Vector{Int64Particle}(undef, threads)
  for i in 1:threads
    scratchInt64Particle[i] = Int64Particle()
  end
  return RPIO{Particle}(N, n, threads, potentials, particles, ks, scratch,
    scratchInt, scratchInt64Particle)
end


text:
using FourVectors
using FourVectors.StaticArrays
using Test

@test FourVector(1,2,3; t=4) == FourVector(MVector(1,2,3,4.0))
@test FourVector(1,2,3; t=4) == Particle(p=[1,2,3], E=4.0)
@test Particle(0,0,3; msq=4^2) == Particle(p=[0,0,3], E=5.0)
#
@test mass(Particle(1,2,3; msq=4)) ≈ 2

@test collect(FourVector(1,2,3; t=4)) == [1,2,3,4]

x = Particle(1.1,2.2,3.0; msq=4)
newp = [1.1,1.1,1.1]
x.p .= newp
@test invmasssq(x) == x.E^2 - sum(abs2,newp)

@test x.E == x.T == x[0]== x[4]
@test x.Px == x.X == x.x == x.px
@test x.Py == x.Y == x.y == x.py
@test x.Pz == x.Z == x.z == x.pz

@test length(x) == 4
@test mass(x * 1000) ≈ 1000*mass(x)

@test sphericalangles(Particle(0,4,3; msq=4^2)) == (3/5,π/2)
@test boostfactor(Particle(4,0,0; msq=3^2)) == 5/3

x1,x2 = x[[1,2]]
x[[1,2]] .= x[[2,1]]
@test x[[1,2]] == [x2,x1]


text:
function oneOrderHigher(diag::Diagram{W}, ::Type{Id}, subdiagram=[]) where {W,Id}
    if diag.id isa PropagatorId && (diag.id isa Id) == false
        #for bare propagator, a derivative of different id vanishes
        return nothing
    end
    id = deepcopy(diag.id)
    if Id == BareGreenId
        id.order[1] += 1
    elseif Id == BareInteractionId
        id.order[2] += 1
    else
        error("not implemented!")
    end

    d = Diagram{W}(id, diag.operator, subdiagram; name=diag.name, factor=diag.factor, weight=diag.weight)
    return d
end

"""
    function derivative(diags::Union{Diagram,Tuple,AbstractVector}, ::Type{ID}) where {ID<:PropagatorId}

    Automatic differentiation derivative on the diagrams

# Arguments
- diags     : diagrams to take derivative
- ID        : DiagramId to apply the differentiation
"""
function derivative(diags::Union{Diagram,Tuple,AbstractVector}, ::Type{ID}) where {ID<:PropagatorId}
    # use a dictionary to host the dual diagram of a diagram for a given hash number
    # a dual diagram is defined as the derivative of the original diagram

    single = false
    if diags isa Diagram
        diags = [diags,]
        single = true
    end
    dual = Dict{Int,Any}()
    for diag in diags
        for d in PostOrderDFS(diag)
            if haskey(dual, d.hash)
                continue
            end
            id = d.id
            if id isa PropagatorId
                # for propagators like bare Green's function and interaction, derivative simply means increase an order by one
                dual[d.hash] = oneOrderHigher(d, ID)
            else # composite diagram
                if d.operator isa Sum
                    # for a diagram which is a sum of subdiagrams, derivative means a sub of derivative subdiagrams
                    children = [dual[sub.hash] for sub in d.subdiagram if isnothing(dual[sub.hash]) == false]
                    if isempty(children)
                        dual[d.hash] = nothing
                    else
                        dual[d.hash] = oneOrderHigher(d, ID, children)
                    end
                elseif d.operator isa Prod
                    # d = s1xs2x... = s1'xs2x... + s1xs2'x... + ...
                    terms = []
                    for (si, sub) in enumerate(d.subdiagram)
                        if isnothing(dual[sub.hash])
                            continue
                        end
                        children = deepcopy(d.subdiagram)
                        children[si] = dual[sub.hash]
                        dd = oneOrderHigher(d, ID, children)
                        if isnothing(dd) == false
                            push!(terms, dd)
                        end
                    end
                    if isempty(terms)
                        dual[d.hash] = nothing
                    else
                        newid = deepcopy(id)
                        newid.order .= terms[1].id.order
                        dual[d.hash] = Diagram{id.para.weightType}(newid, Sum(), terms, name=Symbol("$(d.name)'"))
                    end
                else
                    error("not implemented!")
                end
            end
        end
    end
    if single
        return dual[diags[1].hash]
    else
        return [dual[diag.hash] for diag in diags]
    end
end

"""
    function derivative(diags::Union{Diagram,Tuple,AbstractVector}, ::Type{ID}, order::Int) where {ID<:PropagatorId}

    Automatic differentiation derivative on the diagrams

# Arguments
- diags      : diagrams to take derivative
- ID         : DiagramId to apply the differentiation
- order::Int : derivative order
"""
function derivative(diags::Union{Diagram,Tuple,AbstractVector}, ::Type{ID}, order::Int) where {ID<:PropagatorId}
    @assert order >= 0
    if order == 0
        return diags
    end
    result = diags
    for o in 1:order
        result = derivative(result, ID)
    end
    return result
end

"""
    function removeHatreeFock!(diags::Union{Diagram,Tuple,AbstractVector})

    Remove the Hatree-Fock insertions that without any derivatives on the propagator and the interaction.

# Arguments
- diags      : diagrams to remove the Fock insertion

# Remarks
- The operations removeHatreeFock! and taking derivatives doesn't commute with each other!
- If the input diagram is a Hatree-Fock diagram, then the overall weight will become zero!
"""
function removeHatreeFock!(diags::Union{Diagram,Tuple,AbstractVector})
    single = false
    if diags isa Diagram
        diags = [diags,]
        single = true
    end
    for diag in diags
        for d in PreOrderDFS(diag)
            # for subdiag in d.subdiagram
            if d.id isa SigmaId
                # println(d, " with ", d.id.para.innerLoopNum)
                if isempty(d.id.order) || all(x -> x == 0, d.id.order) #eithr order is empty or a vector of zeros
                    if d.id.para.innerLoopNum == 1
                        d.factor = 0.0
                    end
                end
            end
            # end
        end
    end
    if single
        return diags[1]
    else
        return diags
    end
end


text:
## distribution.jl : unit testing for distribution.jl
#
# This file is part of MultilevelEstimators.jl - A Julia toolbox for
# Multilevel Monte Carlo Methods (c) Pieterjan Robbe, 2019

@testset "Uniform                      " begin

distr1 = Uniform()
distr2 = Uniform(0.,1.)
distr3 = Uniform(0,1.)
@test typeof(distr1) <: Uniform
@test typeof(distr2) <: Uniform
@test typeof(distr3) <: Uniform
@test_throws ArgumentError Uniform(Inf, 1.)
@test_throws ArgumentError Uniform(0., NaN)
@test_throws ArgumentError Uniform(1., 0.)
@sprintf "%s" distr1
for xᵢ in 0.1:0.1:0.9
    @test transform(Uniform(), xᵢ) == xᵢ
end

end

@testset "Normal                       " begin

distr1 = Normal()
distr2 = Normal(0.,1.)
distr3 = Normal(0,1.)
@test typeof(distr1) <: Normal
@test typeof(distr2) <: Normal
@test typeof(distr3) <: Normal
@test_throws ArgumentError Normal(Inf, 1.)
@test_throws ArgumentError Normal(0., NaN)
@test_throws ArgumentError Normal(0., -1.)
@sprintf "%s" distr1
y = [-1.2815515655446004, -0.8416212335729143, -0.5244005127080409, -0.2533471031357997, 0.0,
      0.2533471031357997,  0.5244005127080407,  0.8416212335729143,  1.2815515655446004]
for (xᵢ, yᵢ) in Base.Iterators.zip(0.1:0.1:0.9, y)
    @test transform(Normal(), xᵢ) ≈ yᵢ
end

end

@testset "TruncatedNormal              " begin

distr1 = TruncatedNormal()
distr2 = TruncatedNormal(0.,1.)
distr3 = TruncatedNormal(0.,1.,-1.,1.)
distr4 = TruncatedNormal(0.,1.,-2//3,10)
@test typeof(distr1) <: TruncatedNormal
@test typeof(distr2) <: TruncatedNormal
@test typeof(distr3) <: TruncatedNormal
@test typeof(distr4) <: TruncatedNormal
@test_throws ArgumentError TruncatedNormal(Inf, 1.)
@test_throws ArgumentError TruncatedNormal(0., NaN)
@test_throws ArgumentError TruncatedNormal(0., 1., Inf, 2.)
@test_throws ArgumentError TruncatedNormal(0., 1., -2., Inf)
@test_throws ArgumentError TruncatedNormal(0., -1.)
@test_throws ArgumentError TruncatedNormal(1., 0.)
@sprintf "%s" distr1
y = [ -1.184032466693905, -0.7938201191289552, -0.4984028824219386, -0.2415871851410768, 0.0,
      0.2415871851410768,  0.4984028824219383,  0.7938201191289559,  1.1840324666939055]
for (xᵢ, yᵢ) in Base.Iterators.zip(0.1:0.1:0.9, y)
    @test transform(TruncatedNormal(), xᵢ) ≈ yᵢ
end

end

@testset "Weibull                      " begin

distr1 = Weibull()
distr2 = Weibull(2.,1.)
distr3 = Weibull(2,1.)
@test typeof(distr1) <: Weibull
@test typeof(distr2) <: Weibull
@test typeof(distr3) <: Weibull
@test_throws ArgumentError Weibull(Inf, 1.)
@test_throws ArgumentError Weibull(0., NaN)
@test_throws ArgumentError Weibull(1., -1)
@test_throws ArgumentError Weibull(-1., 1)
@sprintf "%s" distr1
y = [  0.4590436050264207, 0.6680472308365775,  0.8446004309005917,  1.0107676525947897, 1.1774100225154747,
       1.3537287260556712, 1.5517556536555206,  1.7941225779941017,  2.1459660262893476]
for (xᵢ, yᵢ) in Base.Iterators.zip(0.1:0.1:0.9, y)
    @test transform(Weibull(), xᵢ) ≈ yᵢ
end

end


text:
module TestRollCall
    using Base.Test
    import RollCallDataIO

    encoding = RollCallDataIO.RollCallEncoding()
    votes = [1 2 3; 4 5 6; 7 8 9; 0 0 0]
    roll_call = RollCallDataIO.RollCall(votes, encoding)

    @test isa(roll_call, RollCallDataIO.RollCall)
    @test isa(roll_call.votes, Matrix{Int})
    @test isa(roll_call.encoding, RollCallDataIO.RollCallEncoding)
end

max_stars_repo_name: johnmyleswhite/RollCallDataIO

max_issues_repo_name: johnmyleswhite/RollCallDataIO

max_forks_repo_name: johnmyleswhite/RollCallDataIO


text:
# ---
# title: 761. Special Binary String
# id: problem761
# author: Tian Jun
# date: 2020-10-31
# difficulty: Hard
# categories: String, Recursion
# link: <https://leetcode.com/problems/special-binary-string/description/>
# hidden: true
# ---
#
# _Special_ binary strings are binary strings with the following two properties:
#
# * The number of 0's is equal to the number of 1's.
# * Every prefix of the binary string has at least as many 1's as 0's.
#
# Given a special string `S`, a _move_ consists of choosing two consecutive,
# non-empty, special substrings of `S`, and swapping them. _(Two strings are
# consecutive if the last character of the first string is exactly one index
# before the first character of the second string.)_
#
# At the end of any number of moves, what is the lexicographically largest
# resulting string possible?
#
# **Example 1:**
#
#
#
#     Input: S = "11011000"
#     Output: "11100100"
#     Explanation:
#     The strings "10" [occuring at S[1]] and "1100" [at S[3]] are swapped.
#     This is the lexicographically largest string possible after some number of swaps.
#
#
# **Note:**
#
#   1. `S` has length at most `50`.
#   2. `S` is guaranteed to be a _special_ binary string as defined above.
#
#
## @lc code=start
using LeetCode

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


text:
## Based on PyTorch tutorial: https://pytorch.org/tutorials/beginner/blitz/neural_networks_tutorial.html#sphx-glr-beginner-blitz-neural-networks-tutorial-py
## Note that that tutorial and this file use dummy inputs just to demonstrate API,
## for fully working example see cnn_mnist2.jl

using Avalon

# include("../src/core.jl")
# __init__()

mutable struct Net
    conv1::Conv2d
    conv2::Conv2d
    fc1::Linear
    fc2::Linear
    fc3::Linear
end

Net() = Net(
    Conv2d(1, 6, 3),
    Conv2d(6, 16, 3),
    Linear(16 * 6 * 6, 120),
    Linear(120, 84),
    Linear(84, 10)
)

function (m::Net)(x::AbstractArray)
    x = maxpool2d(relu.(m.conv1(x)), (2, 2))
    x = maxpool2d(relu.(m.conv2(x)), (2, 2))
    # prod(size(x)[1:3])
    x = reshape(x, 576, :)
    x = relu.(m.fc1(x))
    x = relu.(m.fc2(x))
    x = m.fc3(x)
    return x
end

function main()
    m = Net()
    X = rand(32, 32, 1, 1000);
    Y = rand(10, 1000);
    loss_fn = MSELoss()
    fit!(m, X, Y, loss_fn; lr=1e-3)
end


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

@testset "termination" begin
  # Different infeasibility information scenarios.
  infeas_info_no_infeasibility1 = FirstOrderLp.InfeasibilityInformation()

  infeas_info_no_infeasibility2 = FirstOrderLp.InfeasibilityInformation()
  infeas_info_no_infeasibility2.primal_ray_linear_objective = -1.0
  infeas_info_no_infeasibility2.primal_ray_quadratic_norm = 1.0
  infeas_info_no_infeasibility2.max_dual_ray_infeasibility = 1.0

  infeas_info_dual_infeasible = deepcopy(infeas_info_no_infeasibility1)
  infeas_info_dual_infeasible.primal_ray_linear_objective = -1.0
  infeas_info_primal_infeasible = deepcopy(infeas_info_no_infeasibility1)
  infeas_info_primal_infeasible.dual_ray_objective = 1.0

  # Test infeasibility criteria.
  eps_primal_infeasible = 1e-6
  eps_dual_infeasible = 1e-6

  @testset "primal infeasibility" begin
    @test !FirstOrderLp.primal_infeasibility_criteria_met(
      eps_primal_infeasible,
      infeas_info_no_infeasibility1,
    )
    @test !FirstOrderLp.primal_infeasibility_criteria_met(
      eps_primal_infeasible,
      infeas_info_no_infeasibility2,
    )
    @test !FirstOrderLp.primal_infeasibility_criteria_met(
      eps_primal_infeasible,
      infeas_info_dual_infeasible,
    )
    @test FirstOrderLp.primal_infeasibility_criteria_met(
      eps_primal_infeasible,
      infeas_info_primal_infeasible,
    )
  end

  @testset "dual infeasibility" begin
    @test !FirstOrderLp.dual_infeasibility_criteria_met(
      eps_dual_infeasible,
      infeas_info_no_infeasibility1,
    )
    @test !FirstOrderLp.dual_infeasibility_criteria_met(
      eps_dual_infeasible,
      infeas_info_no_infeasibility2,
    )
    @test FirstOrderLp.dual_infeasibility_criteria_met(
      eps_dual_infeasible,
      infeas_info_dual_infeasible,
    )
    @test !FirstOrderLp.dual_infeasibility_criteria_met(
      eps_dual_infeasible,
      infeas_info_primal_infeasible,
    )
  end

  # Different convergence information  stats scenarios.
  convergence_information_optimal = FirstOrderLp.ConvergenceInformation()
  convergence_information_optimal.primal_objective = 1.0
  convergence_information_optimal.dual_objective = 1.0
  convergence_information_optimal.l_inf_primal_variable = 1.0
  convergence_information_optimal.l2_primal_variable = 1.0
  convergence_information_optimal.l_inf_dual_variable = 2.0
  convergence_information_optimal.l2_dual_variable = 2.0

  convergence_info_dont_terminate1 = deepcopy(convergence_information_optimal)
  convergence_info_dont_terminate1.primal_objective = 10.0
  convergence_info_dont_terminate2 = deepcopy(convergence_information_optimal)
  convergence_info_dont_terminate2.l_inf_primal_residual = 1.0
  convergence_info_dont_terminate2.l2_primal_residual = 1.0
  convergence_info_dont_terminate3 = deepcopy(convergence_information_optimal)
  convergence_info_dont_terminate3.l_inf_dual_residual = 1.0
  convergence_info_dont_terminate3.l2_dual_residual = 1.0

  iteration_stats_optimal = FirstOrderLp.IterationStats()
  iteration_stats_optimal.iteration_number = 5
  iteration_stats_optimal.cumulative_kkt_matrix_passes = 100.5
  iteration_stats_optimal.cumulative_time_sec = 5.0
  iteration_stats_optimal.convergence_information =
    [convergence_information_optimal]
  iteration_stats_optimal.infeasibility_information =
    [infeas_info_no_infeasibility1]

  iteration_stats_dont_terminate1 = FirstOrderLp.IterationStats()
  iteration_stats_dont_terminate1.iteration_number = 5
  iteration_stats_dont_terminate1.cumulative_kkt_matrix_passes = 100.5
  iteration_stats_dont_terminate1.cumulative_time_sec = 5.0
  iteration_stats_dont_terminate1.convergence_information =
    [convergence_info_dont_terminate1]
  iteration_stats_dont_terminate1.infeasibility_information =
    [infeas_info_no_infeasibility1]

  qp_cache = FirstOrderLp.cached_quadratic_program_info(example_qp())

  for optimality_norm in [FirstOrderLp.L_INF, FirstOrderLp.L2]
    eps_optimal_relative = 1e-4
    eps_optimal_absolute = 1e-4

    @testset "optimality" begin
      # Test optimality criteria.
      @test !FirstOrderLp.optimality_criteria_met(
        optimality_norm,
        eps_optimal_relative,
        eps_optimal_absolute,
        convergence_info_dont_terminate1,
        qp_cache,
      )
      @test !FirstOrderLp.optimality_criteria_met(
        optimality_norm,
        eps_optimal_relative,
        eps_optimal_absolute,
        convergence_info_dont_terminate2,
        qp_cache,
      )
      @test !FirstOrderLp.optimality_criteria_met(
        optimality_norm,
        eps_optimal_relative,
        eps_optimal_absolute,
        convergence_info_dont_terminate3,
        qp_cache,
      )
      @test FirstOrderLp.optimality_criteria_met(
        optimality_norm,
        eps_optimal_relative,
        eps_optimal_absolute,
        convergence_information_optimal,
        qp_cache,
      )
    end

    full_termination_criteria = FirstOrderLp.construct_termination_criteria(
      optimality_norm = optimality_norm,
      eps_optimal_absolute = eps_optimal_absolute,
      eps_optimal_relative = eps_optimal_relative,
      eps_primal_infeasible = eps_primal_infeasible,
      eps_dual_infeasible = eps_dual_infeasible,
      time_sec_limit = 100.0,
      iteration_limit = 10,
      kkt_matrix_pass_limit = 10000.0,
    )

    @testset "check_termination_criteria" begin
      @test FirstOrderLp.check_termination_criteria(
        full_termination_criteria,
        qp_cache,
        iteration_stats_optimal,
      ) == FirstOrderLp.TERMINATION_REASON_OPTIMAL

      @test !FirstOrderLp.check_termination_criteria(
        full_termination_criteria,
        qp_cache,
        iteration_stats_dont_terminate1,
      )

      full_termination_criteria.time_sec_limit = 1.0
      @test FirstOrderLp.check_termination_criteria(
        full_termination_criteria,
        qp_cache,
        iteration_stats_dont_terminate1,
      ) == FirstOrderLp.TERMINATION_REASON_TIME_LIMIT

      full_termination_criteria.time_sec_limit = 10.0
      full_termination_criteria.iteration_limit = 1

      @test FirstOrderLp.check_termination_criteria(
        full_termination_criteria,
        qp_cache,
        iteration_stats_dont_terminate1,
      ) == FirstOrderLp.TERMINATION_REASON_ITERATION_LIMIT

      full_termination_criteria.iteration_limit = 10
      full_termination_criteria.kkt_matrix_pass_limit = 40.0
      @test FirstOrderLp.check_termination_criteria(
        full_termination_criteria,
        qp_cache,
        iteration_stats_dont_terminate1,
      ) == FirstOrderLp.TERMINATION_REASON_KKT_MATRIX_PASS_LIMIT
    end
  end
end


text:
module Xsum

export xsum, XAccumulator

# Load xsum libraries from our deps.jl
const depsjl_path = joinpath(dirname(@__FILE__), "..", "deps", "deps.jl")
if !isfile(depsjl_path)
    error("Xsum not installed properly, run Pkg.build(\"Xsum\"), restart Julia and try again")
end
include(depsjl_path)

function __init__()
    check_deps()
end

include("xsum_h.jl")

"""
`acc = XAccumulator()` creates an object `acc` that can be used to accumulate an
exactly rounded double-precision sum.   The sum is accumulated by repeatedly
calling `accumulate!(acc, x)` on summands `x`, and then the final exactly-rounded
result is obtained by calling `float(acc)`.
"""
struct XAccumulator <: Real
    acc::Base.RefValue{xsum_small_accumulator}
    XAccumulator(acc::Base.RefValue{xsum_small_accumulator}) = new(acc)
    function XAccumulator()
        acc = Base.RefValue{xsum_small_accumulator}()
        xsum_small_init(acc)
        return new(acc)
    end
end
XAccumulator(s::XAccumulator) = XAccumulator(s.acc)
Base.float(s::XAccumulator) = xsum_small_round(s.acc)
XAccumulator(init::Real) = accumulate!(XAccumulator(), init)
Base.show(io::IO, s::XAccumulator) = print(io, "XAccumulator(", float(s), ")")

"""
    accumulate!(s::XAccumulator, x::Real)

Add `Float64(x)` to the accumulated sum `s`, with the sum performed in (effectively) infinite precision.
"""
Base.accumulate!(s::XAccumulator, x::Real) = begin; xsum_small_add1(s.acc, x); s; end

Base.accumulate!(s::Complex{XAccumulator}, x::Complex) = Complex(accumulate!(real(s), real(x)), accumulate!(imag(s), imag(x)))
Base.accumulate!(s::Complex{XAccumulator}, x::Real) = Complex(accumulate!(real(s), x), imag(s))
Base.accumulate!(s::XAccumulator, x::Complex) = Complex(accumulate!(s, x), XAccumulator())
XAccumulator(init::Complex) = Complex(XAccumulator(real(init)), XAccumulator(imag(init)))

"""
    xsum([f,] itr)

Compute the exactly rounded double-precision sum of calling `f` (which defaults to `identity`) on
the elements of the iterable collection `itr`. That is, return the sum as if (1) the elements of `f.(itr)`
were converted to double (`Float64`) precision, (2) summed in *infinite* precision, and
(3) rounded to the closest double-precision value in the final result.    Both real and complex
sums are supported.
"""
function xsum(f, itr)
    i = iterate(itr)
    if i === nothing
        zero = Base.mapreduce_empty_iter(f, +, itr, Base.IteratorEltype(itr))
        return zero isa Real ? Float64(zero) : zero isa Complex ? ComplexF64(zero) : zero
    end
    v, state = i
    s = XAccumulator(f(v))
    while true
        i = iterate(itr, state)
        i === nothing && return float(s)
        v, state = i
        accumulate!(s, f(v))
    end
end

xsum(itr) = xsum(identity, itr)

function xsum(a::StridedVector{Float64})
    stride(a,1) != 1 && return xsum(identity, x)
    n = length(a)
    if n < 256 # empirical threshold for small accumulator to be faster
        let acc = Base.RefValue{xsum_small_accumulator}()
            xsum_small_init(acc)
            xsum_small_addv(acc, a)
            return xsum_small_round(acc)
        end
    else
        let acc = Base.RefValue{xsum_large_accumulator}()
            xsum_large_init(acc)
            if n ≤ typemax(Cint)
                xsum_large_addv(acc, a)
            else # libxsum needs arrays in typemax(Cint) chunks
                i = firstindex(a)
                l = i + typemax(Cint) - 1
                e = lastindex(a)
                while l ≤ e
                    xsum_large_addv(acc, @view a[i:l])
                    i = l + 1
                    l = i + typemax(Cint) - 1
                end
                xsum_large_addv(acc, @view a[i:e])
            end
            return xsum_large_round(acc)
        end
    end
end

end # module


text:
# This file was generated, do not modify it. # hide
r = range(stack, :(regressor2.lambda), lower = 1, upper = 20, scale=:log)
tuned_stack = TunedModel(model=stack,
                         ranges=r,
                         tuning=Grid(),
                         measure=rms,
                         resampling=Holdout())

mach = fit!(machine(tuned_stack,  X, y), verbosity=0)
best_stack = fitted_params(mach).best_model
best_stack.regressor2.lambda


text:
## ----------------------------------------------------------------------------
function update_remote_backup(;force = false, loglevel = :verbose)

    bkdir = _remote_backup_dir()
    remote_url = _remote_url()

    _update_remote_backup(bkdir, remote_url; force, loglevel)

end