text:
using Plots, Test, Colors

const PLOTS_DEFAULTS = Dict(:theme => :wong2, :fontfamily => :palantino)
Plots.__init__()

@testset "Loading theme" begin
    pl = plot(1:5)
    @test pl[1][1][:seriescolor] == RGBA(colorant"black")
    @test Plots.guidefont(pl[1][:xaxis]).family == "palantino"
end

empty!(PLOTS_DEFAULTS)
Plots.__init__()

@testset "Legend defaults" begin
    p = plot()
    @test p[1][:legend_font_family] == "sans-serif"
    @test p[1][:legend_font_pointsize] == 8
    @test p[1][:legend_font_halign] == :hcenter
    @test p[1][:legend_font_valign] == :vcenter
    @test p[1][:legend_font_rotation] == 0.0
    @test p[1][:legend_font_color] == RGB{Colors.N0f8}(0.0, 0.0, 0.0)
    @test p[1][:legend_position] == :best
    @test p[1][:legend_title] == nothing
    @test p[1][:legend_title_font_family] == "sans-serif"
    @test p[1][:legend_title_font_pointsize] == 11
    @test p[1][:legend_title_font_halign] == :hcenter
    @test p[1][:legend_title_font_valign] == :vcenter
    @test p[1][:legend_title_font_rotation] == 0.0
    @test p[1][:legend_title_font_color] == RGB{Colors.N0f8}(0.0, 0.0, 0.0)
    @test p[1][:legend_background_color] == RGBA{Float64}(1.0, 1.0, 1.0, 1.0)
    @test p[1][:legend_foreground_color] == RGB{Colors.N0f8}(0.0, 0.0, 0.0)
end # testset

@testset "Legend API" begin
    p = plot(;
        legendfontfamily = "serif",
        legendfontsize = 12,
        legendfonthalign = :left,
        legendfontvalign = :top,
        legendfontrotation = 1,
        legendfontcolor = :red,
        legend = :outertopleft,
        legendtitle = "The legend",
        legendtitlefontfamily = "helvetica",
        legendtitlefontsize = 3,
        legendtitlefonthalign = :right,
        legendtitlefontvalign = :bottom,
        legendtitlefontrotation = -5.2,
        legendtitlefontcolor = :blue,
        background_color_legend = :cyan,
        foreground_color_legend = :green,
    )
    @test p[1][:legend_font_family] == "serif"
    @test p[1][:legend_font_pointsize] == 12
    @test p[1][:legend_font_halign] == :left
    @test p[1][:legend_font_valign] == :top
    @test p[1][:legend_font_rotation] == 1.0
    @test p[1][:legend_font_color] == :red
    @test p[1][:legend_position] == :outertopleft
    @test p[1][:legend_title] == "The legend"
    @test p[1][:legend_title_font_family] == "helvetica"
    @test p[1][:legend_title_font_pointsize] == 3
    @test p[1][:legend_title_font_halign] == :right
    @test p[1][:legend_title_font_valign] == :bottom
    @test p[1][:legend_title_font_rotation] == -5.2
    @test p[1][:legend_title_font_color] == :blue
    @test p[1][:legend_background_color] == RGBA{Float64}(0.0, 1.0, 1.0, 1.0)
    @test p[1][:legend_foreground_color] == RGBA{Float64}(0.0, 0.5019607843137255, 0.0, 1.0)

    #setting whole font
    sp = plot(
        1:5,
        legendfont = font(12),
        legend_font_halign = :left,
        foreground_color_subplot = :red,
    )[1]
    @test Plots.legendfont(sp).pointsize == 12
    @test Plots.legendfont(sp).halign == :left
    # match mechanism
    @test sp[:legend_font_color] == sp[:foreground_color_subplot]
    @test Plots.legendfont(sp).color == sp[:foreground_color_subplot]
end # testset


text:
function tvlqr(A, B, Q, R)
    T = length(Q)

    P = [zero(A[1]) for t = 1:T]
    K = [zero(B[1]') for t = 1:T-1]
    P[T] = Q[T]

    for t = T-1:-1:1
        K[t] = (R[t] + B[t]' * P[t+1] *  B[t]) \ (B[t]' * P[t+1] * A[t])
        P[t] = (Q[t] + K[t]' * R[t] * K[t]
                + (A[t] - B[t] * K[t])' * P[t+1] * (A[t] - B[t] * K[t]))
    end

    return K, P
end

function reference_gains(s::Simulation, traj::ContactTraj, obj::TrackingVelocityObjective;
		N::Int=10, κ=2e-4, U_scaling=100, V_scaling=100)

	model = s.model
	H = traj.H
	z = [[deepcopy(traj.z) for i = 1:N]...;]
	θ = [[deepcopy(traj.θ) for i = 1:N]...;]

	rz = [LinearizedStep(s, z[t], θ[t], κ).rz for t = 1:N*H]
	rθ = [LinearizedStep(s, z[t], θ[t], κ).rθ for t = 1:N*H]

	idx_z = indices_z(s)
	idx_θ = indices_θ(model, nf=1)

	∂z∂θ = [-rz[t] \ rθ[t] for t = 1:N*H]

	∂q2∂q1 = [zeros(model.nq, model.nq) for t = 1:N*H]
	∂q2∂q2 = [I(model.nq) for t = 1:N*H]
	∂q3∂q1 = [∂z∂θ[t][idx_z.q, idx_θ.q1] for t = 1:N*H]
	∂q3∂q2 = [∂z∂θ[t][idx_z.q, idx_θ.q2] for t = 1:N*H]

	∂q2∂u = [zeros(model.nq, model.nu) for t = 1:N*H]
	∂q3∂u = [∂z∂θ[t][idx_z.q, idx_θ.u] for t = 1:N*H]

	A = [[∂q2∂q1[t] ∂q2∂q2[t]; ∂q3∂q1[t] ∂q3∂q2[t]] for t = 1:N*H]
	B = [[∂q2∂u[t]; ∂q3∂u[t]] for t = 1:N*H]
	Q = [
		[obj.q[1]+V_scaling*obj.v[1]    -V_scaling*obj.v[1];
		 -V_scaling*obj.v[1]    obj.q[1]+V_scaling*obj.v[1]
		 ] for i = 1:N*H+1]
	R = [U_scaling * obj.u[1] for i = 1:N*H]

	K, P = tvlqr(A, B, Q, R)
	return K[1:H]
end

ext: jl


text:
using PLext
using Test
using TestSetExtensions

@testset ExtendedTestSet "PLext.jl" begin
    @includetests ARGS
end

max_stars_repo_name: mlxd/PLext

max_issues_repo_name: mlxd/PLext

max_forks_repo_name: mlxd/PLext


text:
struct LoadConfig
    k_range::UnitRange
    o_range::UnitRange
    blacklist::Vector{Chain_ID}
    whitelist::Vector{Chain_ID}
    function LoadConfig(k_range, o_range; blacklist=Vector{Chain_ID}(), whitelist=Vector{Chain_ID}())
        assert_loadconfig(k_range, o_range) && new(k_range, o_range, blacklist, whitelist)
    end
end

function assert_loadconfig(k_range,o_range)
    k_range[1] < 1 && throw(ArgumentError("Minimum value of LoadConfig k_range is 1!"))
    o_range[1] < 0 && throw(ArgumentError("Minimum value of LoadConfig o_range is 0!"))
    return true
end

#subfunc to handle balancing memory load on dissimilar machines in cluster
function load_balancer(no_models::Integer, hmm_jobs::RemoteChannel, config::LoadConfig, timeout::Integer=300)
    jobid::Chain_ID, start_iterate::Integer, hmm::BHMM, job_norm::AbstractFloat, observations::Matrix = take!(hmm_jobs)

    starttime=time()
    while (jobid.K < config.k_range[1] || jobid.K > config.k_range[end] || jobid.order < config.o_range[1] || jobid.order > config.o_range[end] || jobid in config.blacklist || (length(config.whitelist)>0 && !(jobid in config.whitelist))) && (time()-starttime <= timeout) #while a job prohibited by load table, keep putting the job back and drawing a new one
        put!(hmm_jobs, (jobid, start_iterate, hmm, job_norm, observations))
        jobid, start_iterate, hmm, job_norm, observations = take!(hmm_jobs)
    end
    (time()-starttime > timeout) ? (return 0,0,0,0,0) : (return jobid, start_iterate, hmm, job_norm, observations)
end


text:
# push!(LOAD_PATH, ".")
import JSON
using NNS
using ArgParse
import NNS.save, NNS.load

include("src/io.jl")
include("src/textmodel.jl")
include("src/tfidf.jl")
include("src/diskdata.jl")

function create_model(args)
    if length(args["qlist"]) > 0
        qlist = parse.(split(args["qlist"], ','))
    else
        qlist = []
    end

    if length(args["nlist"]) > 0
        nlist = parse.(split(args["nlist"], ','))
    else
        nlist = []
    end

    key = args["key"]
    input = args["input"]

    corpus = Vector{String}()
    itertweets(input == "-" ? STDIN : input) do tweet
        push!(corpus, tweet[key])
    end

    model = TfidfModel(corpus, TextConfig(qlist, nlist))
    path = joinpath(args["env"], "model")
    savepath(path, model)
end

function convert_database(args)
    model = loadpath(joinpath(args["env"], "model"), TfidfModel)
    input = args["input"]
    output = args["output"]
    out = output == "-" ? STDOUT : open(output, "w")
    key = args["key"]
    itertweets(input == "-" ? STDIN : input) do tweet
        save(out, vectorize(tweet[key], model))
    end

    output != "-" && close(out)
end

### index related functions
function set_locality_oracle_handler(args, index, model)
    env = args["env"]
    fastlinks = KnnResult[]
    filename = joinpath(env, "fastlinks")
    if isfile(filename)
        open(filename) do f
            while !eof(f)
                links = load(f, KnnResult)
                push!(fastlinks, links)
            end
        end
    else
        num = args["fastlinks"]
        for i in 1:length(model.voc)
            push!(fastlinks, KnnResult(num))
        end
    end

    function locality_oracle(q::DocumentType)::Vector{Int32}
        L = Int32[]
        for term in q.terms
          for p in fastlinks[term.id]
            push!(L, p.objID)
          end
        end
        L
    end

    index.locality_oracle = locality_oracle
    fastlinks
end

function append_index(args)
    env = args["env"]
    key = args["key"]
    model = loadpath(joinpath(env, "model"), TfidfModel)
    index = loadindex(args)
    fastlinks = set_locality_oracle_handler(args, index, model)
    input = args["database"]
    wcorpus = DiskWriter(joinpath(env, "corpus"))
    fdb = open(joinpath(env, "db"), "a")
    istream = input == "-" ? STDIN : open(input)
    count = 0

    while !eof(istream)
        line = readline(istream)
        tweet = parsetweet(line)
        v = vectorize(tweet[key], model)
        length(v.terms) == 0 && continue
        links = push!(index, v)
        count += 1
	save(fdb, v)
        save(fdb, links)
        push!(wcorpus, line)
        if log(2, count) % 1 == 0
            info("saving a snapshot of the index")
            flush(fdb)
            flush(wcorpus)
            saveindex(args, index, fastlinks)
        end

        objID = length(index.db) |> Int32

        for i in 1:length(v.terms)
            termID = v.terms[i].id
            push!(fastlinks[termID], objID, -v.terms[i].weight)
        end
    end

    input != "-" && close(istream)
    close(wcorpus)
    close(fdb)
    saveindex(args, index, fastlinks)
end

function search_index(args)
    env = args["env"]
    key = args["key"]

    index = loadindex(args)
    model = loadpath(joinpath(env, "model"), TfidfModel)
    fastlinks = set_locality_oracle_handler(args, index, model)

    k = args["knn"]
    if args["opt-recall"] > 0
        _k = args["opt-knn"] > 0 ? args["opt-knn"] : k
        _db = args["opt-queries"] == "" ? index.db : loaddatabase(args["opt-queries"])
        perf = Performance(_db, index.dist)
        optimize!(index, recall=args["opt-recall"], k=_k, perf=perf)
    end

    # t = time()
    raw_output = args["raw-output"]
    rcorpus = DiskReader(joinpath(env, "corpus"))
    input = args["queries"]
    istream = input == "-" ? STDIN : open(input)

    while !eof(istream)
        tweet = readline(istream) |> parsetweet
        q = vectorize(tweet[key], model)
        if length(q.terms) == 0
            warn("no terms were found:", tweet[args["key"]])
            continue
        end

        res = search(index, q, KnnResult(k))
        if raw_output
            tweet["knn"] = res
        else
            L = []
            for p in res
                s = rcorpus[p.objID] |> JSON.parse
                push!(L, (p, s))
            end
            tweet["knn"] = L
        end
        write(STDOUT, JSON.json(tweet, 2))

        push!(L, tweet)
    end
    # t = time() - t
    # info("Solved $count queries in ", t, "seconds, avg search time: ", t/count)

    input != "-" && close(istream)
end

function main()
    s = ArgParseSettings(description = "Ranix command line")

    @add_arg_table s begin
        "model"
            action = :command
            help = "creates a new model"
        "convert"
            action = :command
            help = "converts an json-like input to a sparse vector output"
        "append"
            action = :command
            help = "append items to the index"
        "search"
            action = :command
            help = "opens a repl for searching over a given index"
    end

    @add_arg_table s["model"] begin
      "env"
        required = true
        arg_type = String
        help = "The output environment"
      "--input"
        arg_type = String
        default = "-"
        help = "The input dataset"
      "--qlist"
        required = true
        help = "a comma separated list of q sizes (character q-grams)"
        arg_type = String
      "--nlist"
        required = true
        help = "a comma separated list of n sizes (word n-grams)"
        arg_type = String
      "--key"
        help = "The keyword that contains the text data to be modeled"
        arg_type = String
        default = "text"
    end

    @add_arg_table s["convert"] begin
        "env"
            required = true
            arg_type = String
            help = "index's environment"
        "--input"
            arg_type = String
            default = "-"
            help = "The input dataset"
        "--output"
            arg_type = String
            default = "-"
            help = "filename to store the output"
        "--key"
            arg_type = String
            help = "The keyword that contains the text data"
            default = "text"
    end

    @add_arg_table s["append"] begin
        "env"
            required = true
            arg_type = String
            help = "index's environment"
        "--database"
            arg_type = String
            default = "-"
            help = "input database"
        "--key"
            arg_type = String
            help = "The keyword that contains the text data"
            default = "text"
        "--knn"
            arg_type = Int
            default = 10
            help = "the index should be optimized to solve k queries"
        "--recall"
            arg_type = Float64
            default = 0.90
            help = "the index will be optimized to solve queries with this expected recall"
        "--fastlinks"
            arg_type = Int
            default = 3
            help = "the number of fast links"
    end

    @add_arg_table s["search"] begin
        "env"
            required = true
            arg_type = String
            help = "index's environment"
        "--queries"
            arg_type = String
            default = "-"
            help = "the input queries, in database format"
        "--knn"
            arg_type = Int
            default = 10
            help = "the number of neighbors to be retrieved"
        "--key"
            arg_type = String
            help = "The keyword that contains the text data"
            default = "text"
        "--raw-output"
            action = :store_true
            help = "The output is shown without accessing the corpus"
        "--opt-recall"
            arg_type = Float64
            default = 0.0
            help = "if given, the search parameters will be adjusted to try to reach this recall"
        "--opt-knn"
            arg_type = Int
            default = -1
            help = "if given, the search parameters will be adjusted to try to reach --opt-recall for this number of neighbors"
        "--opt-queries"
            arg_type = String
            default = ""
            help = "a small sample of items following the query set distribution (useful to optimize parameters --opt-*)"
    end
    args = parse_args(s)
    cmd = args["%COMMAND%"]

    if cmd == "model"
        create_model(args["model"])
    elseif cmd == "convert"
        convert_database(args["convert"])
    elseif cmd == "append"
        append_index(args["append"])
    elseif cmd == "search"
        search_index(args["search"])
    end
end

if !isinteractive()
    textindex = main()
end

max_stars_repo_name: INGEOTEC/ranix

max_issues_repo_name: INGEOTEC/ranix

max_forks_repo_name: INGEOTEC/ranix


text:
# Graph → Syntax

toexpr(f, xs...) = :($f($(xs...)))

binding(bindings::AbstractDict, v) =
  haskey(bindings, v) ? bindings[v] : (bindings[v] = gensym("edge"))

function syntax(head::DVertex; bindconst = !isfinal(head))
  vs = topo(head)
  ex, bs = :(;), IdDict{Any ,Any}()
  for v in vs
    x = toexpr(value(v), [binding(bs, n) for n in inputs(v)]...)
    if !bindconst && isconstant(v) && nout(v) > 1
      bs[v] = v.value.value
    elseif nout(v) > 1 || (!isfinal(head) && v ≡ head)
      edge = binding(bs, v)
      push!(ex.args, :($edge = $x))
    elseif haskey(bs, v)
      if MacroTools.inexpr(ex, bs[v])
        ex = MacroTools.replace(ex, bs[v], x)
      else
        push!(ex.args, :($(bs[v]) = $x))
      end
    else
      isfinal(v) ? push!(ex.args, x) : (bs[v] = x)
    end
  end
  head ≢ vs[end] && push!(ex.args, binding(bs, head))
  return ex
end

# TODO: this is butt ugly

function constructor(g)
  vertex = isa(g, DVertex) ? :(DataFlow.dvertex) : :(DataFlow.vertex)
  g = mapv(g) do v
    value(v) isa Lambda && (v.value = :(DataFlow.Lambda($(value(v).args), $(constructor(value(v).body)))))
    value(v) isa Constant && (v.value = :(DataFlow.Constant($(value(v).value))))
    prethread!(v, typeof(v)(Constant(value(v))))
    prethread!(v, typeof(v)(Constant(vertex)))
    v.value = Call()
    v
  end
  ex = syntax(g)
  decls, exs = [], []
  for x in block(ex).args
    if @capture(x, v_ = $vertex(f_, a__))
      push!(decls, :($v = $vertex($f)))
      push!(exs, :(DataFlow.thread!($v, $(a...))))
    else
      push!(exs, x)
    end
  end
  return :($(decls...);$(exs...))
end


text:
using StippleUI
using Test

@testset "StippleUI.jl" begin

    @test 1 == 1

end


text:
#using Distributed
#@everywhere using MambaModels
using MambaModels
gr(size=(400,400))

# Data
line = Dict{Symbol, Any}()

howell1 = CSV.read(rel_path("..", "data", "Howell1.csv"), delim=';')
df = convert(DataFrame, howell1);

df2 = filter(row -> row[:age] >= 18, df);
mean_weight = mean(df2[:weight])
df2[:weight_c] = convert(Vector{Float64}, df2[:weight]) .- mean_weight ;
line[:x] = convert(Array{Float64,1}, df2[:weight_c]);
line[:y] = convert(Array{Float64,1}, df2[:height]);
line[:xmat] = convert(Array{Float64,2}, [ones(length(line[:x])) line[:x]])

model = Model(
  y = Stochastic(1,
    (xmat, beta, s2) -> MvNormal(xmat * beta, sqrt(s2)),
    false
  ),
  beta = Stochastic(1, () -> MvNormal([178, 0], [sqrt(10000), sqrt(100)])),
  s2 = Stochastic(() -> Uniform(0, 50))
)

inits = [
  Dict{Symbol, Any}(
    :y => line[:y],
    :beta => [rand(Normal(178, 100)), rand(Normal(0, 10))],
    :s2 => rand(Uniform(0, 50))
  )
  for i in 1:3
]

scale1 = [0.5, 0.25]
summary1 = identity
eps1 = 0.5

scale2 = 0.5
summary2 = x -> [mean(x); sqrt(var(x))]
eps2 = 0.1

scheme = [
  Mamba.NUTS([:beta]),
  Mamba.Slice([:s2], 10)
]

setsamplers!(model, scheme)

chn = mcmc(model, line, inits, 10000, burnin=1000, chains=3)

describe(chn)

chn2 = MCMCChains.Chains(chn.value, String.(chn.names))

MCMCChains.describe(chn2)

MCMCChains.plot(chn2)

# This file was generated using Literate.jl, https://github.com/fredrikekre/Literate.jl


text:
using Revise
using MetaGen
using Gen
using Profile
using StatProfilerHTML
using GenRFS

#Profile.init(; n = 10^4, delay = 1e-5)

#GenRFS.modify_partition_ctx!(1000)

#call it
#@profilehtml gt_trace,_ = Gen.generate(metacog, (possible_objects,))
#@profilehtml gt_trace,_ = Gen.generate(metacog, (possible_objects,))
num_frames = 1
num_videos = 2
params = Video_Params(n_possible_objects = 2)

@time gt_trace,_ = Gen.generate(main, (false, num_videos, num_frames, params));

@profilehtml Gen.generate(main, (num_videos, num_frames, params));

#Profile.init(; n = 10^7, delay = 1e-6)
@profilehtml gt_trace,_ = Gen.generate(main, (num_videos, num_frames, params));
#println(gt_trace)
#gt_choices = get_choices(gt_trace)

max_stars_repo_name: zhangir-azerbayev/MetaGen

max_issues_repo_name: zhangir-azerbayev/MetaGen

max_forks_repo_name: zhangir-azerbayev/MetaGen


text:
# Run Julia Script in Project Environment

```jl
julia --startup-file=no --project=~/amzn-repo/axam/axam_wbx src/SIM/AXAM-37/M2-M3-report.jl

```

max_stars_repo_name: mahiki/julia-is-great

max_issues_repo_name: mahiki/julia-is-great

max_forks_repo_name: mahiki/julia-is-great


text:
@testset "sizing" begin
    for dims in (
        (0, 0), (1, 1), (2, 1), (1, 2), (2, 2), (3, 6), (6, 3), (9, 4), (4, 9),
        (5, 0), (0, 5), (5, 1), (1, 5), (10, 0), (0, 10), (6, 8), (8, 6),
        (10, 10), (10, 15), (15, 10), (2000, 200), (200, 2000), (2000, 2000)
    )
        seed!(RNG, 1337)
        p = @inferred heatmap(randn(RNG, dims))
        @test p isa Plot
        test_ref("references/heatmap/default_$(dims[1])x$(dims[2]).txt", @show_col(p, :displaysize=>T_SZ))
    end
end

@testset "axis scaling" begin
    x = repeat(collect(0:10), outer=(1, 11))
    p = @inferred heatmap(x, xfact=0.1)
    test_ref("references/heatmap/scaling_$(size(x, 1))x$(size(x, 2))_xfact_0.1.txt", @show_col(p, :displaysize=>T_SZ))
    p = @inferred heatmap(x, xfact=0.1, xoffset=-0.5)
    test_ref("references/heatmap/scaling_$(size(x, 1))x$(size(x, 2))_xfact_0.1_xoffset-0.5.txt", @show_col(p, :displaysize=>T_SZ))
    p = @inferred heatmap(x, yfact=0.1)
    test_ref("references/heatmap/scaling_$(size(x, 1))x$(size(x, 2))_yfact_0.1.txt", @show_col(p, :displaysize=>T_SZ))
    p = @inferred heatmap(x, yfact=0.1, yoffset=-0.5)
    test_ref("references/heatmap/scaling_$(size(x, 1))x$(size(x, 2))_yfact_0.1_yoffset-0.5.txt", @show_col(p, :displaysize=>T_SZ))
    p = @inferred heatmap(x, xfact=0.1, xoffset=-0.5, yfact=10, yoffset=-50)
    test_ref(
        "references/heatmap/scaling_$(size(x, 1))x$(size(x, 2))_xfact_0.1_xoffset-0.5_yfact_10_yoffset_-50.txt",
        @show_col(p, :displaysize=>T_SZ)
    )
end

@testset "integer values" begin
    x = repeat(collect(1:20), outer=(1, 20))
    p = @inferred heatmap(x)
    test_ref("references/heatmap/integers_$(size(x, 1))x$(size(x, 2)).txt", @show_col(p, :displaysize=>T_SZ))
end

@testset "color maps" begin
    x = collect(0:30) * collect(0:30)'
    for cmap in keys(UnicodePlots.COLOR_MAP_DATA)
        p = @inferred heatmap(x, colormap=cmap)
        test_ref("references/heatmap/colormap_$(size(x, 1))x$(size(x, 2))_$cmap.txt", @show_col(p, :displaysize=>T_SZ))
    end
    p = @inferred heatmap(x, colormap=reverse(UnicodePlots.COLOR_MAP_DATA[:viridis]))
    test_ref("references/heatmap/colormap_$(size(x, 1))x$(size(x, 2))_reverse_viridis.txt", @show_col(p, :displaysize=>T_SZ))

    seed!(RNG, 1337)
    rgb = rand(RNG, 20, 20, 3)
    img = RGB.(rgb[:, :, 1], rgb[:, :, 2], rgb[:, :, 3])
    p = @inferred heatmap(img)  # RGB Matrix
    test_ref("references/heatmap/rgb.txt", @show_col(p, :displaysize=>T_SZ))
end

@testset "axis limits" begin
    x = collect(0:30) * collect(0:30)'
    p = @inferred heatmap(x, ylim=[10, 20])
    test_ref("references/heatmap/limits_$(size(x, 1))x$(size(x, 2))_ylim_10_20.txt", @show_col(p, :displaysize=>T_SZ))
    p = @inferred heatmap(x, ylim=[50, 50])
    test_ref("references/heatmap/limits_$(size(x, 1))x$(size(x, 2))_ylim_50_50.txt", @show_col(p, :displaysize=>T_SZ))
    p = @inferred heatmap(x, xlim=[10, 20])
    test_ref("references/heatmap/limits_$(size(x, 1))x$(size(x, 2))_xlim_10_20.txt", @show_col(p, :displaysize=>T_SZ))
    p = @inferred heatmap(x, xlim=[50, 50])
    test_ref("references/heatmap/limits_$(size(x, 1))x$(size(x, 2))_xlim_50_50.txt", @show_col(p, :displaysize=>T_SZ))
    p = @inferred heatmap(x, xlim=[10, 20], ylim=[10, 20])
    test_ref("references/heatmap/limits_$(size(x, 1))x$(size(x, 2))_xlim_10_20_ylim_10_20.txt", @show_col(p, :displaysize=>T_SZ))
    p = @inferred heatmap(x, xlim=[0, 0.1], xfact=0.01)
    test_ref("references/heatmap/limits_$(size(x, 1))x$(size(x, 2))_xfact_0.01_xlim_0_0.1.txt", @show_col(p, :displaysize=>T_SZ))
    p = @inferred heatmap(x, ylim=[0, 0.1], yfact=0.01)
    test_ref("references/heatmap/limits_$(size(x, 1))x$(size(x, 2))_yfact_0.01_ylim_0_0.1.txt", @show_col(p, :displaysize=>T_SZ))
    p = @inferred heatmap(x, xlim=[0, 0.1], xfact=0.01, ylim=[0.1, 0.2], yfact=0.01)
    test_ref(
        "references/heatmap/limits_$(size(x, 1))x$(size(x, 2))_xfact_0.01_xlim_0_0.1_yfact_0.01_ylim_0.1_0.2.txt",
        @show_col(p, :displaysize=>T_SZ)
    )
    p = @inferred heatmap(x, ylim=[1, 50])
    test_ref("references/heatmap/limits_$(size(x, 1))x$(size(x, 2))_ylim_1_50.txt", @show_col(p, :displaysize=>T_SZ))
    p = @inferred heatmap(x, xlim=[1, 50])
    test_ref("references/heatmap/limits_$(size(x, 1))x$(size(x, 2))_xlim_1_50.txt", @show_col(p, :displaysize=>T_SZ))
    p = @inferred heatmap(x, xlim=[1, 50], ylim=[1, 50])
    test_ref("references/heatmap/limits_$(size(x, 1))x$(size(x, 2))_xlim_1_50_ylim_1_50.txt", @show_col(p, :displaysize=>T_SZ))
    zmax = 2 * maximum(x)
    p = @inferred heatmap(x, zlim=[0, zmax])
    test_ref("references/heatmap/limits_$(size(x, 1))x$(size(x, 2))_zlim_0_$zmax.txt", @show_col(p, :displaysize=>T_SZ))
end

@testset "all zero values" begin
    x = zeros(20, 20)
    p = @inferred heatmap(x)
    test_ref("references/heatmap/zeros_$(size(x, 1))x$(size(x, 2)).txt", @show_col(p, :displaysize=>T_SZ))
end

@testset "parameters" begin
    seed!(RNG, 1337)
    p = @inferred heatmap(randn(RNG, 200, 200), colorbar=false)
    test_ref("references/heatmap/parameters_200x200_no_colorbar.txt", @show_col(p, :displaysize=>T_SZ))
    seed!(RNG, 1337)
    p = @inferred heatmap(randn(RNG, 200, 200), labels=false)
    test_ref("references/heatmap/parameters_200x200_no_labels.txt", @show_col(p, :displaysize=>T_SZ))
    seed!(RNG, 1337)
    p = heatmap(randn(RNG, 200, 200), title="Custom Title", zlabel="Custom Label", colorbar_border=:ascii, colormap=:inferno)
    test_ref("references/heatmap/parameters_200x200_zlabel_ascii_border.txt", @show_col(p, :displaysize=>T_SZ))
    seed!(RNG, 1337)
    p = @inferred heatmap(randn(RNG, 10, 10), xfact=0.1, colormap=:inferno)
    test_ref("references/heatmap/parameters_10x10_xfact_inferno.txt", @show_col(p, :displaysize=>T_SZ))
    seed!(RNG, 1337)
    p = @inferred heatmap(randn(RNG, 10, 11), xfact=0.1, colormap=:inferno)
    test_ref("references/heatmap/parameters_10x11_xfact_inferno.txt", @show_col(p, :displaysize=>T_SZ))
    seed!(RNG, 1337)
    p = @inferred heatmap(randn(RNG, 10, 10), yfact=1, colormap=:inferno)
    test_ref("references/heatmap/parameters_10x10_yfact_inferno.txt", @show_col(p, :displaysize=>T_SZ))
end


text:
using TropicalTensors
using Test

@testset "dump load" begin
    jt = Randn()
    ht = Randn()
    rsg = rand_spinglass(Float64, SquareLattice(3,3); jt=jt, ht=ht, dumpto=".")
    rsg2 = rand_spinglass(Float64, SquareLattice(3,3))
    load_params(".", rsg2, jt, ht, 2)
    @test rsg.Js ≈ rsg2.Js
    @test rsg.hs ≈ rsg2.hs
end


text:
# Michiel Stock
# Example of a source code file implementing a module.

# all your code is part of the module you are implementing
module Example

# you have to import everything you need for your module to work
# if you use a new package, don't forget to add it in the package manager
using LinearAlgebra

# export all functions that are relevant for the user
export solve_quadratic_system, quadratic_function

"""
    solve_quadratic_system(P, q, r=0; testPD=false)

Find the minimizer of a canonical quadratic system:

``f(x) = 0.5xᵀ  P  x + x ⋅ q + r``

with `P` a positive scalar or a positive-definite n x n matrix and `q` a scalar or
a n x n vector. The intercept `r` can optionally be given but does not influence
the minimizer. One can use the keyword argument `testPD` to test if `P` is positive definite,
though this comes at a computational cost.

## Examples

Scalar case:

```julia-repl
julia> solve_quadratic_system(8, -4, 3)
0.5
```

Vector case:

```julia-repl
julia> P = [3 1; 1 2];

julia> q = [0.5, -2];

julia> solve_quadratic_system(P, q)
2-element Array{Float64,1}:
 -0.6
  1.2999999999999998
```
"""
function solve_quadratic_system(P, q, r=0; testPD=false)
    testPD && @assert isposdef(P) "Provided P is not positive definite"
    return - P \ q
end

"""
    quadratic_function(P, q, r)

Give the paramters of a canonical quatratic function and returns a function.
"""
quadratic_function(P, q, r) = x -> 0.5x' * P * x + q ⋅ x + r

end

max_stars_repo_name: kirstvh/STMOZOO

max_issues_repo_name: kirstvh/STMOZOO

max_forks_repo_name: kirstvh/STMOZOO


text:
"""
save resonance groupings of a compound.
Choices for `f`` are: `real()`, `imag()`, or `abs()`. These corresponds to real part, imaginary part, and magnitude spectrum, respectively.
"""
function plotgroups(title_string::String,
    P,
    U,
    q,
    qs,
    q_singlets,
    f::Function,
    return_val_type::T;
    canvas_size::Tuple{Int,Int} = (1000, 400)) where T

    U_rad = U .* (2*π)

    q_U = q.(U_rad)
    plot_obj = Plots.plot( P,
        f.(q_U),
        title = title_string,
        label = "sum of sub-models",
        seriestype = :line,
        ticks = :native,
        xlims = (P[1],P[end]),
        hover = P,
        linewidth = 4,
        xflip = true,
        size = canvas_size)

    qs_U = Vector{Vector{Vector{Complex{T}}}}(undef, length(qs))
    for i = 1:length(qs)

        qs_U[i] = Vector{Vector{Complex{T}}}(undef, length(qs[i]))
        for k = 1:length(qs[i])

            qs_U[i][k] = qs[i][k].(U_rad)

            Plots.plot!(plot_obj, P, f.(qs_U[i][k]), label = "sub-group ($(i),$(k))",
                seriestype = :line,
                linestyle = :dot,
                xflip = true,
                linewidth = 4)
        end
    end

    q_singlets_U = q_singlets.(U_rad)
    Plots.plot!(plot_obj, P, f.(q_singlets_U), label = "group of all singlets",
        seriestype = :line,
        linestyle = :dot,
        xflip = true,
        linewidth = 4)
    #
    #Plots.plot!(plot_obj, P, f.(q_U),
    #markershape = :circle,
    #seriestype = :scatter,
    #xflip = true)

    return plot_obj, q_U, qs_U, q_singlets_U
end

"""
only plots the first compound in molecule_names.
"""
function plotgroupsfullscript(plot_title, molecule_names,
    H_params_path, project_path, tol_coherence, α_relative_threshold,
    Δc_partition_radius, Δcs_max, κ_λ_lb, κ_λ_ub;
    display_flag = false,
    prune_low_signal_for_display_flag::Bool = false,
    display_reduction_factor = 100,
    display_threshold_factor =  α_relative_threshold/10,
    canvas_size = (1600, 900),
    reset_plot_dir_flag = true,
    plot_imag_and_mag_flag = false,
    save_plot_flag = true)

    load_path = joinpath(project_path, "experiment.bson")
    save_folder = joinpath(project_path, "plots")
    isdir(save_folder) || mkpath(save_folder)

    ### load block.
    dic = BSON.load(load_path)
    fs = dic[:fs]
    SW = dic[:SW]
    ν_0ppm = dic[:ν_0ppm]
    λ_0ppm = dic[:λ_0ppm]

    hz2ppmfunc = uu->(uu - ν_0ppm)*SW/fs
    ppm2hzfunc = pp->(ν_0ppm + pp*fs/SW)

    Δcs_max_mixture = collect( Δcs_max for i = 1:length(molecule_names))

    # get a surrogate where K_{n,i} is encouraged to be no larger than `early_exit_part_size`.
    # println("Timing: setupmixtureproxies()")
    # @time mixture_params = NMRSpectraSimulator.setupmixtureproxies(molecule_names,
    mixture_params = NMRSpectraSimulator.setupmixtureproxies(molecule_names,
        H_params_path, ppm2hzfunc, fs, SW,
        λ_0ppm, ν_0ppm, dummy_SSFID;
        tol_coherence = tol_coherence,
        α_relative_threshold = α_relative_threshold,
        Δc_partition_radius = Δc_partition_radius)
    As = mixture_params
    A = As[1]

    ## frequency locations. For plotting.
    ΩS_ppm = NMRSpectraSimulator.getPsnospininfo(mixture_params, hz2ppmfunc)
    ΩS_ppm_sorted = sort(NMRSpectraSimulator.combinevectors(ΩS_ppm))

    u_offset = 0.5 # in units ppm.
    u_min = ppm2hzfunc(ΩS_ppm_sorted[1] - u_offset)
    u_max = ppm2hzfunc(ΩS_ppm_sorted[end] + u_offset)

    P = LinRange(hz2ppmfunc(u_min), hz2ppmfunc(u_max), 50000)
    U = ppm2hzfunc.(P)
    U_rad = U .* (2*π)

    ## end frequency locations.

    # println("Timing: fitproxies()")
    # @time NMRSpectraSimulator.fitproxies!(As;
    NMRSpectraSimulator.fitproxies!(As;
        κ_λ_lb = κ_λ_lb,
        κ_λ_ub = κ_λ_ub,
        u_min = u_min,
        u_max = u_max,
        Δr = 1.0,
        Δκ_λ = 0.05)

    # store fitted surrogate for future use.
    save_simulation_path = joinpath(project_path, "simulation.bson")

    #NMRSpectraSimulator.removeauxinfo!(mixture_params)
    #BSON.bson(save_simulation_path, mixture_params = mixture_params)

    # create the functions for each resonance group.
    qs = collect( collect( ωω->A.qs[i][k](ωω-A.ss_params.d[i], A.ss_params.κs_λ[i]) for k = 1:length(A.qs[i]) ) for i = 1:length(A.qs) )
    q_singlets = ωω->NMRSpectraSimulator.evalsinglets(ωω, A.d_singlets, A.αs_singlets, A.Ωs_singlets, A.β_singlets, A.λ0, A.κs_λ_singlets)

    # create the function for the entire compound.
    q = uu->NMRSpectraSimulator.evalitpproxymixture(uu, As[1:1])

    # evaluate at the plotting positions.
    q_U = q.(U_rad)

    qs_U = collect( collect( qs[i][k].(U_rad) for k = 1:length(qs[i]) ) for i = 1:length(qs) )
    q_singlets_U = q_singlets.(U_rad)

    # sanity check.
    q_check_U = q_singlets_U
    if !isempty(qs) # some compounds only have singlets.
        q_check_U += sum( sum( qs[i][k].(U_rad) for k = 1:length(qs[i]) ) for i = 1:length(qs) )
    end
    discrepancy = norm(q_check_U- q_U)

    #println("sanity check. This should be numerically zero: ", discrepancy)
    @assert discrepancy < 1e-14

    # plot.
    P_display = P
    U_display = U
    if prune_low_signal_for_display_flag

        inds, _ = NMRSpectraSimulator.prunelowsignalentries(q_U, display_threshold_factor, display_reduction_factor)
        P_display = P[inds]
        U_display = U[inds]
    end

    #println("length(P) = ", length(P))
    #println("length(P_display) = ", length(P_display))

    # should be clear the plotting directory before saving plots?
    if reset_plot_dir_flag && save_plot_flag
        paths_to_delete = readdir(save_folder, join = true)
        rm.(paths_to_delete)
    end

    plots_save_path = joinpath(save_folder, "groups_real.html")
    title_string = "$(plot_title), $(molecule_names), real spectrum"
    plot_obj, q_U, qs_U, q_singlets_U = plotgroups(title_string, P_display, U_display, q, qs, q_singlets, real, P[1]; canvas_size = canvas_size)
    if save_plot_flag
        Plots.savefig(plot_obj, plots_save_path)
    end
    if display_flag
        display(plot_obj)
    end
    #println("length(q_U) = ", length(q_U))
    #@assert 1==44

    if plot_imag_and_mag_flag

        plots_save_path = joinpath(save_folder, "groups_imag.html")
        title_string = "$(plot_title), $(molecule_names), imaginary spectrum"
        plot_obj, q_U, qs_U, q_singlets_U = plotgroups(title_string, P_display, U_display, q, qs, q_singlets, imag, P[1]; canvas_size = canvas_size)
        if save_plot_flag
            Plots.savefig(plot_obj, plots_save_path)
        end
        if display_flag
            display(plot_obj)
        end

        plots_save_path = joinpath(save_folder, "groups_magnitude.html")
        title_string = "$(plot_title), $(molecule_names), magnitude spectrum"
        plot_obj, q_U, qs_U, q_singlets_U = plotgroups(title_string, P_display, U_display, q, qs, q_singlets, abs, P[1]; canvas_size = canvas_size)
        if save_plot_flag
            Plots.savefig(plot_obj, plots_save_path)
        end
        if display_flag
            display(plot_obj)
        end
    end
end

"""
root_folder should contain only folders of 1D 1H NMR compound standard experiments, with the compound name being the subfolder name.
I.e., a folder containing an experiment of L-Serine and L-Histidine should have two folders in it with those names that contain the corresponding experimental data.
"""
function batchplotgroups(plot_title, root_path,
    H_params_path, tol_coherence,
    α_relative_threshold, Δc_partition_radius, Δcs_max, κ_λ_lb, κ_λ_ub;
    display_flag = false,
    prune_low_signal_for_display_flag::Bool = false,
    display_reduction_factor = 100,
    display_threshold_factor =  α_relative_threshold/10,
    canvas_size = (1600, 900),
    reset_plot_dir_flag = true,
    plot_imag_and_mag_flag = false,
    save_plot_flag = true)

    # get a list of experiment paths and compounds names. Assume the folder name of the experiment contains the compound.
    tmp = readdir(root_path, join = true)
    inds = findall(xx->isdir(xx), tmp)

    project_paths = tmp[inds]
    experiment_names = readdir(root_path)[inds]
    project_names = collect( "$(experiment_names[i])" for i = 1:length(experiment_names))

    for i = 1:length(project_names)

        project_path = project_paths[i]
        compound_name = experiment_names[i]

        println("Now on ", project_path)

        # make sure this compound is in our library.
        records = GISSMOReader.getGISSMOentriesall()
        record_names = collect( records[i].molecule_name for i = 1:length(records) )

        k = findfirst(xx->xx==compound_name, record_names)

        # decide if we should simulate and plot.
        if typeof(k) != Nothing
            plotgroupsfullscript(plot_title, [compound_name;],
                H_params_path, project_path, tol_coherence, α_relative_threshold,
                Δc_partition_radius, Δcs_max, κ_λ_lb, κ_λ_ub;
                display_flag = display_flag,
                prune_low_signal_for_display_flag = prune_low_signal_for_display_flag,
                display_reduction_factor = display_reduction_factor,
                display_threshold_factor = display_threshold_factor,
                canvas_size = canvas_size,
                reset_plot_dir_flag = reset_plot_dir_flag,
                plot_imag_and_mag_flag = plot_imag_and_mag_flag,
                save_plot_flag = save_plot_flag)
        else
            println("Compound not in the local GISSMO library. Skip.")
        end
    end

end

####


text:
# Single first-order reaction
function foc(τ, (q, α, τ₋,  τ₊, k))
    exp(-k*τ)
end

# Competing first-order production and consumption reactions
function fopac(τ, kprod, kcons, substrate0, product0)
    (substrate0*kprod/(kcons - kprod))*(exp(-kprod*τ) - exp(-kcons*τ)) + product0*exp(-kcons*τ)
end


text:
# # Testing the stability of a CKSVAR model

# Here, we test the stability of the CKSVAR model estimated for Japan in:
# - Ikeda, Li, Mavroeidis and Zanetti (2021), "[Testing the effectiveness of monetary policy in Japan and the United States](https://arxiv.org/abs/2012.15158)." Working paper.
# This paper is hereafter referred to as ILMZ21.
##
# In ILMZ21, the effectiveness of unconventional monetary policy in Japan is investigated via estimation of a CKSVAR(2) model in inflation, output gap and the Call Rate for Japan, over the period 1985q3-2019q1, with shadow call rate as the latent variable.
#
# Converting this model into TAR form, we determine upper bounds on the state-constrained joint spectral radius (SCJSR)
##
using DataFrames, DataFramesMeta, CSV, HTTP
using ThresholdStability

extract_est(parname, df) = @subset(df, in([parname]).(:Parameter)).:Estimate[1]

JP = CSV.read(HTTP.get("https://raw.githubusercontent.com/samwycherley/ThresholdStability.jl/master/examples/src/estimates/estimates_JP.csv").body, DataFrame, copycols=true)  # Japan data from CSV
##
##
# This loads the Japan estimates as a dataframe and defines a function to extract estimates. Next, we build the matrices:
##
##
JPβtilde = zeros(2)
for i in 1:2
    JPβtilde[i] = extract_est("tildebeta_$i", JP)
end

JPCstar = zeros(3, 2)
for i in 1:3
    for j in 1:2
        JPCstar[i, j] = extract_est("Eq.$i lSR_$j", JP)
    end
end

Cbars = []
for j in 1:2
    Cbar_j = zeros(3, 3)
    for i in 1:3
        Cbar_j[i, 1] = extract_est("Eq.$i INFL_$j", JP)
        Cbar_j[i, 2] = extract_est("Eq.$i YGAP_BOJ_$j", JP)
        Cbar_j[i, 3] = extract_est("Eq.$i SR_$j", JP)
    end
    push!(Cbars, Cbar_j)
end
JPCbar = hcat(Cbars[1], Cbars[2])

JPC = copy(JPCbar); k = size(JPC, 1)  # convert C̄ to C
for i in 1:2
    JPC[:, i*k] -= JPCstar[:, i]
end

Σ_JP, X_JP = CKSVAR_to_TAR(JPC, JPCstar, JPβtilde, 2)
##
# Write this as a state-dependent switched linear system:
##
##
G_JP = automaton_constructor(Σ_JP)
s_JP = discreteswitchedsystem(Σ_JP, G_JP, X_JP)
##
@show sosbound_γ(s_JP, 2)
##
# We see that the model estimated for Japan is stable


text:
"""
```
relaxation!(ral, xₙ₋₁, Ψₙ₋₁; tol = 1e-10, max_iters = 1000, damping = .5, pnorm = Inf,
            schur_fnct = schur!, autodiff = :central, use_anderson = false, m = 5,
            verbose = :none, kwargs...)
```

solves for the coefficients ``(z, y, \\Psi)`` of a risk-adjusted linearization by the following relaxation algorithm:

1. Initialize guesses for ``(z, y, \\Psi)``
2. Do until convergence

    a) Solve for ``(z, y)`` using the expectational and state transition equations and fixing ``\\Psi``.

    b) Use a QZ decomposition to solve for ``\\Psi`` while fixing ``(z, y)``.

### Types:
- `S1 <: Number`
- `S2 <: Real`
- `S3 <: Real`

### Inputs
- `m::RiskAdjustedLinearization`: object holding functions needed to calculate
    the risk-adjusted linearization
- `xₙ₋₁::AbstractVector{S1}`: initial guess for ``(z, y)``
- `Ψₙ₋₁::AbstractVector{S1}`: initial guess for ``\\Psi``

### Keywords
- `tol::S2`: convergence tolerance of residual norm for relaxation algorithm
- `max_iters::Int`: maximumm number of iterations
- `damping::S2`: guesses are updated as the weighted average
    `xₙ = damping * proposal + (1 - damping) * xₙ₋₁`.
- `pnorm::S3`: norm for residual tolerance
- `schur_fnct::Function`: function for calculating the Schur factorization during QZ decomposition
- `autodiff::Symbol`: specifies whether to use autoamtic differentiation in `nlsolve`
    (and is the same keyword as the `autodiff` keyword for `nlsolve`)
- `use_anderson::Bool`: set to true to apply Anderson acceleration to the
    fixed point iteration of the relaxation algorithm
- `m::Int`: `m` coefficient if using Anderson acceleration
- `sparse_jacobian::Bool = false`: if true, exploit sparsity in the Jacobian in calls to `nlsolve` using SparseDiffTools.jl.
    If `jac_cache` and `sparsity` are `nothing`, then `relaxation!` will attempt to determine the sparsity pattern.
- `sparsity::Union{AbstractArray, Nothing} = nothing`: sparsity pattern for the Jacobian in calls to `nlsolve`
- `colorvec = nothing`: matrix coloring vector for sparse Jacobian in calls to `nlsolve`
- `jac_cache = nothing`: pre-allocated Jacobian cache for calls to `nlsolve` during the numerical algorithms
- `sparsity_detection::Bool = false`: If true, use SparsityDetection.jl to detect sparsity pattern (only relevant if
    both `jac_cache` and `sparsity` are `nothing`). If false,  then the sparsity pattern is
    determined by using finite differences to calculate a Jacobian and assuming any zeros will always be zero.
    Currently, SparsityDetection.jl fails to work.
- `verbose::Symbol`: verbosity of information printed out during solution.
    a) `:low` -> statement when homotopy continuation succeeds
    b) `:high` -> statement when homotopy continuation succeeds and for each successful iteration
"""
function relaxation!(ral::RiskAdjustedLinearization, xₙ₋₁::AbstractVector{S1}, Ψₙ₋₁::AbstractMatrix{S1};
                     tol::S2 = 1e-10, max_iters::Int = 1000, damping::S2 = .5, pnorm::S3 = Inf,
                     schur_fnct::Function = schur!, autodiff::Symbol = :central,
                     use_anderson::Bool = false, m::Int = 5,
                     sparse_jacobian::Bool = false,  sparsity::Union{AbstractArray, Nothing} = nothing,
                     colorvec = nothing, jac_cache = nothing,
                     sparsity_detection::Bool = true, verbose::Symbol = :none,
                     kwargs...) where {S1 <: Number, S2 <: Real, S3 <: Real}
    # Set up
    err = 1.
    nl  = nonlinear_system(ral)
    li  = linearized_system(ral)
    Nzy = ral.Nz + ral.Ny
    AA  = Matrix{Complex{S1}}(undef, Nzy, Nzy) # pre-allocate these matrices to calculate QZ decomp for Ψ
    BB  = similar(AA)

    # Initialize system of equations
    _my_eqn = (F, x, Ψ, 𝒱) -> _relaxation_equations(F, x, ral, Ψ, 𝒱)

    if use_anderson
        # Some aliases/views will be useful
        zₙ    = ral.z
        yₙ    = ral.y
        Ψₙ    = ral.Ψ
        𝒱ₙ₋₁  = nl[:𝒱_sss]
        J𝒱ₙ₋₁ = li[:JV]

        _anderson_f = function _my_anderson(F::AbstractArray{T}, xₙ₋₁::AbstractVector{T}) where {T <: Number}
            zₙ₋₁  = @view xₙ₋₁[1:ral.Nz]
            yₙ₋₁  = @view xₙ₋₁[(ral.Nz + 1):Nzy]
            Ψₙ₋₁  = @view xₙ₋₁[(Nzy + 1):end]
            Ψₙ₋₁  = reshape(Ψₙ₋₁, ral.Ny, ral.Nz)

            # Calculate entropy terms 𝒱ₙ₋₁, J𝒱ₙ₋₁
            update!(nl, zₙ₋₁, yₙ₋₁, Ψₙ₋₁; select = Symbol[:𝒱]) # updates nl.𝒱_sss
            update!(li, zₙ₋₁, yₙ₋₁, Ψₙ₋₁; select = Symbol[:JV]) # updates li.JV

            # Solve state transition and expectational equations for (zₙ, yₙ), taking 𝒱ₙ₋₁ and Ψₙ₋₁ as given
            solve_steadystate!(ral, vcat(zₙ₋₁, yₙ₋₁), _my_eqn, Ψₙ₋₁, 𝒱ₙ₋₁; autodiff = autodiff, # updates ral.z and ral.y
                               sparse_jacobian = sparse_jacobian,
                               sparsity = sparsity, colorvec = colorvec,
                               jac_cache = jac_cache,
                               sparsity_detection = sparsity_detection,
                               verbose = verbose, kwargs...)

            # Update Γ₁, Γ₂, Γ₃, Γ₄, given (zₙ, yₙ)
            update!(li, zₙ, yₙ, Ψₙ₋₁; select = Symbol[:Γ₁, :Γ₂, :Γ₃, :Γ₄]) # updates li.Γᵢ

            # QZ decomposition to get Ψₙ, taking Γ₁, Γ₂, Γ₃, Γ₄, and J𝒱ₙ₋₁ as given
            Ψₙ .= compute_Ψ!(AA, BB, li; schur_fnct = schur_fnct)

            # Update zₙ, yₙ, and Ψₙ; then calculate error for convergence check
            zₙ .= (1 - damping) .* zₙ + damping .* zₙ₋₁
            yₙ .= (1 - damping) .* yₙ + damping .* yₙ₋₁
            Ψₙ .= (1 - damping) .* Ψₙ + damping .* Ψₙ₋₁
            err = norm(vcat(zₙ - zₙ₋₁, yₙ - yₙ₋₁, vec(Ψₙ - Ψₙ₋₁)), pnorm)

            # Calculate residual
            F[1:ral.Nz] = zₙ - zₙ₋₁
            F[(ral.Nz + 1):Nzy] = yₙ - yₙ₋₁
            F[(Nzy + 1):end] = vec(Ψₙ - Ψₙ₋₁)

            return F
        end

        out   = nlsolve(_anderson_f, vcat(xₙ₋₁, vec(Ψₙ₋₁)); m = m, ftol = tol, iterations = max_iters)
        count = out.iterations
        if out.f_converged
            update!(ral, out.zero[1:ral.Nz], out.zero[(ral.Nz + 1):Nzy],
                    reshape(out.zero[(Nzy + 1):end], ral.Ny, ral.Nz); update_cache = false)
        end
    else
        count = 1

        # Some aliases/views will be useful
        zₙ₋₁  = @view xₙ₋₁[1:ral.Nz]
        yₙ₋₁  = @view xₙ₋₁[(ral.Nz + 1):end]
        zₙ    = ral.z
        yₙ    = ral.y
        Ψₙ    = ral.Ψ
        𝒱ₙ₋₁  = nl[:𝒱_sss]
        J𝒱ₙ₋₁ = li[:JV]

        while (err > tol) && (count < max_iters)

            # Calculate entropy terms 𝒱ₙ₋₁, J𝒱ₙ₋₁
            update!(nl, zₙ₋₁, yₙ₋₁, Ψₙ₋₁; select = Symbol[:𝒱]) # updates nl.𝒱_sss
            update!(li, zₙ₋₁, yₙ₋₁, Ψₙ₋₁; select = Symbol[:JV]) # updates li.JV

            # Solve state transition and expectational equations for (zₙ, yₙ), taking 𝒱ₙ₋₁ and Ψₙ₋₁ as given
            solve_steadystate!(ral, xₙ₋₁, _my_eqn, Ψₙ₋₁, 𝒱ₙ₋₁; autodiff = autodiff, # updates ral.z and ral.y
                               sparse_jacobian = sparse_jacobian,
                               sparsity = sparsity, colorvec = colorvec,
                               jac_cache = jac_cache,
                               sparsity_detection = sparsity_detection,
                               verbose = verbose, kwargs...)

            # Update Γ₁, Γ₂, Γ₃, Γ₄, given (zₙ, yₙ)
            update!(li, zₙ, yₙ, Ψₙ₋₁; select = Symbol[:Γ₁, :Γ₂, :Γ₃, :Γ₄]) # updates li.Γᵢ

            # QZ decomposition to get Ψₙ, taking Γ₁, Γ₂, Γ₃, Γ₄, and J𝒱ₙ₋₁ as given
            Ψₙ .= compute_Ψ!(AA, BB, li; schur_fnct = schur_fnct)

            # Update zₙ, yₙ, and Ψₙ; then calculate error for convergence check
            zₙ .= (1 - damping) .* zₙ + damping .* zₙ₋₁
            yₙ .= (1 - damping) .* yₙ + damping .* yₙ₋₁
            Ψₙ .= (1 - damping) .* Ψₙ + damping .* Ψₙ₋₁
            err = norm(vcat(zₙ - zₙ₋₁, yₙ - yₙ₋₁, vec(Ψₙ - Ψₙ₋₁)), pnorm)

            # Update zₙ₋₁, yₙ₋₁, and Ψₙ₋₁ (without reallocating them)
            zₙ₋₁ .= zₙ
            yₙ₋₁ .= yₙ
            Ψₙ₋₁ .= Ψₙ

            if verbose == :high
                println("Iteration $(count): error under norm=$(pnorm) is $(err)")
            end

            count += 1
        end
    end

    if count == max_iters
        throw(RALRelaxationError("Relaxation method to find the risk-adjusted linearization did not converge."))
    else
        update!(ral)

        if verbose == :low
            errvec = steady_state_errors(ral)
            println("Convergence achieved after $(count) iterations! Error under norm = $(pnorm) is " *
                    "$(norm(errvec, pnorm)).")
        elseif verbose == :high
            errvec = steady_state_errors(ral)
            println("")
            println("Convergence achieved after $(count) iterations! Error under norm = $(pnorm) is " *
                    "$(norm(errvec, pnorm)).")
        end

        return ral
    end
end

function solve_steadystate!(m::RiskAdjustedLinearization, x0::AbstractVector{S1},
                            f::Function, Ψ::AbstractMatrix{<: Number}, 𝒱::AbstractVector{<: Number};
                            sparse_jacobian::Bool = false, sparsity::Union{AbstractArray, Nothing} = nothing,
                            colorvec = nothing, jac_cache = nothing,
                            sparsity_detection::Bool = true, autodiff::Symbol = :central,
                            verbose::Symbol = :none, kwargs...) where {S1 <: Real, S2 <: Real}

    # Exploit sparsity?
    if sparse_jacobian
        nlsolve_jacobian!, jac =
            construct_sparse_jacobian_function(m, (F, x) -> f(F, x, Ψ, 𝒱), :relaxation, autodiff;
                                               sparsity = sparsity, colorvec = colorvec,
                                               jac_cache = jac_cache, sparsity_detection = sparsity_detection)
        out = nlsolve(OnceDifferentiable((F, x) -> f(F, x, Ψ, 𝒱), nlsolve_jacobian!, x0, copy(x0), jac), x0; kwargs...)
    else
        out = nlsolve(OnceDifferentiable((F, x) -> f(F, x, Ψ, 𝒱), x0, copy(x0), autodiff,
                                         ForwardDiff.Chunk(ForwardDiff.pickchunksize(min(m.Nz, m.Ny)))), x0; kwargs...)
    end

    if out.f_converged
        m.z .= out.zero[1:m.Nz]
        m.y .= out.zero[(m.Nz + 1):end]
    else
        if verbose == :high
            println(out)
        end
        throw(RALRelaxationError())
    end
end

function _relaxation_equations(F::AbstractArray, x::AbstractArray, m::RiskAdjustedLinearization,
                               Ψ::AbstractMatrix{<: Number}, 𝒱::AbstractVector{<: Number})
    # Unpack
    z = @view x[1:m.Nz]
    y = @view x[(m.Nz + 1):end]

    # Update μ(z, y) and ξ(z, y)
    update!(m.nonlinear, z, y, Ψ; select = Symbol[:μ, :ξ])

    # Calculate residuals
    μ_sss             = get_tmp(m.nonlinear.μ.cache, z, y, (1, 1)) # select the first DiffCache b/c that one
    ξ_sss             = get_tmp(m.nonlinear.ξ.cache, z, y, (1, 1)) # corresponds to autodiffing both z and y
    F[1:m.Nz]         = μ_sss - z
    F[(m.Nz + 1):end] = ξ_sss + m.linearization[:Γ₅] * z + m.linearization[:Γ₆] * y + 𝒱
end

mutable struct RALRelaxationError <: Exception
    msg::String
end
RALRelaxationError() =
    RALRelaxationError("A solution for (z, y), given Ψ and 𝒱, to the state transition and expectational equations could not be found.")
Base.showerror(io::IO, ex::RALRelaxationError) = print(io, ex.msg)


text:
s = Set([11, 14, 13, 7, 14, 11])
#> Set([7, 14, 13, 11])

# empty set:
Set() #> Set(Any[])

s1 = Set([11, 25])
s2 = Set([25, 3.14])
union(s1, s2) #> Set([3.14, 25.0, 11.0])
intersect(s1, s2) #> Set([25])
setdiff(s1, s2) #> Set([11])
setdiff(s2, s1) #> Set([3.14])
issubset(s1, s2) #> false
issubset(s1, Set([11, 25, 36])) #> true

push!(s1, 32) #> Set([25,32,11])
in(32, s1) #> true
in(100, s1) #> false

s1 = Set([1,2,3]) #> Set([2,3,1])
typeof(s1) #> Set{Int64}
Set([[1,2,3]]) #> Set(Array{Int64,1}[[1, 2, 3]])

x = Set(collect(1:100))
@time 2 in x
#> 0.003186 seconds (33 allocations: 2.078 KiB)
x2 = Set(collect(1:1000000))
@time 2 in x2
# 0.000003 seconds (4 allocations: 160 bytes)


text:
@inline getindex{T <: FixedVector}(x::T, i::Union{Range, Integer}) = x.(1)[i]
@inline getindex{T <: FixedVectorNoTuple}(x::T, i::Integer) = x.(i)
@inline getindex{N, M, T}(a::Mat{N, M, T}, i::Range, j::Int) = ntuple(IndexFunc(a, j), Val{length(i)})::NTuple{length(i), T}
@inline getindex{N, M, T}(a::Mat{N, M, T}, i::Int, j::Union{Range, Int}) = a.(1)[j][i]
@inline getindex{N, M, T}(a::Mat{N, M, T}, i::Int) = a[ind2sub((N,M), i)...]
@inline getindex(A::FixedArray, I::Tuple) = map(IndexFunctor(A), I)

@inline setindex(a::FixedArray, value, index::Int...) = map(SetindexFunctor(a, value, index), typeof(a))

@inline column{N, T}(v::FixedVector{N, T}) = v
@inline column{R, C, T}(a::Mat{R, C, T}, i::Union{Range, Int}) = a.(1)[i]

@inline row{N, T}(v::FixedVector{N, T}) = Mat{1, N, T}(v...)
@inline row{N, T}(v::FixedVector{N, T}, i::Int) = (v[i],)
@inline row{R, C, T}(a::Mat{R, C, T}, j::Int) = ntuple(IndexFunc(a, j), Val{C})::NTuple{C, T}
@inline row{R, T}(a::Mat{R, 1, T}, j::Int) = (a.(1)[1][j],)
@inline row{R, T}(a::Mat{R, 2, T}, j::Int) = (a.(1)[1][j], a.(1)[2][j])
@inline row{R, T}(a::Mat{R, 3, T}, j::Int) = (a.(1)[1][j], a.(1)[2][j], a.(1)[3][j])
@inline row{R, T}(a::Mat{R, 4, T}, j::Int) = (a.(1)[1][j], a.(1)[2][j], a.(1)[3][j], a.(1)[4][j],)

# the columns of the ctranspose are the complex conjugate rows
@inline crow{R, C, T}(a::Mat{R, C, T}, j::Int) = ntuple(IndexFunc(a, j)', Val{C})::NTuple{C, T}
@inline crow{R, T}(a::Mat{R, 1, T}, j::Int) = (a.(1)[1][j]',)
@inline crow{R, T}(a::Mat{R, 2, T}, j::Int) = (a.(1)[1][j]', a.(1)[2][j]')
@inline crow{R, T}(a::Mat{R, 3, T}, j::Int) = (a.(1)[1][j]', a.(1)[2][j]', a.(1)[3][j]')
@inline crow{R, T}(a::Mat{R, 4, T}, j::Int) = (a.(1)[1][j]', a.(1)[2][j]', a.(1)[3][j]', a.(1)[4][j]',)

# Unified slicing along combinations of fixed and variable dimensions

"Get index of a field with `name` in type `T`.  Should this be in Base?"
@generated function fieldindex{T}(::Type{T}, name::Symbol)
    # Expand all field names inline to allow this to become a constant
    # expression.  Simple alternative is `findfirst(fieldnames(T), name)`
    exprs = [:($(Expr(:quote, n)) == name && return $i) for (i,n) in enumerate(fieldnames(T))]
    quote
        $(Expr(:meta, :inline))
        $(exprs...)
        error("No field \"$name\" in type $T")
    end
end

# Turn `A[1,2,...]` into `destructure(A)[1,2,3,4]`
#
# Turn `A[:x,:y,1,2,...]` into
# `destructure(A)[fieldindex(eltype(A),:x), fieldindex(eltype(A),:y), 1,2,...]`
#
# Also works on the left side of an assignment
function fixed_slice_expr(expr)
    assignrhs = nothing
    if expr.head == :(=)
        @assert length(expr.args) == 2
        assignrhs = expr.args[2]
        expr = expr.args[1]
    end
    expr.head == :ref || error("Array reference not found in expression $expr")
    inds = Any[]
    name = expr.args[1]
    for i = 2:length(expr.args)
        ind = expr.args[i]
        if isa(ind,Expr) && ind.head == :quote
            push!(inds, :(fieldindex(eltype(tmp), $(esc(ind)))))
        else
            push!(inds, esc(ind))
        end
    end
    if assignrhs === nothing
        quote
            tmp = $(esc(name))
            destructure(tmp)[$(inds...)]
        end
    else
        quote
            tmp = $(esc(name))
            destructure(tmp)[$(inds...)] = $(esc(assignrhs))
        end
    end
end

"""
Slice across both fixed and variable size dimensions of
`Array{F<:FixedArray,M}`.  Before slicing, the array is reshaped to the natural
memory ordering as in `destructure()`: the `N=ndims(F)` fixed dimensions come
first, followed by the `M` dimensions of variable length.

Examples:

    # Array of fixed size vectors
    a = [Vec(i,j) for i=1:5, j=1:5]
    @fslice a[2,:,:] = 10
    xcomps = @fslice a[1,:,:]

    # Vector of fixed size matrices
    m = [@fsa([i 0; 0 i^2]) for i=1:4]
    @fslice m[1,1,:]

    # Slice immutables by field name
    immutable MyVec <: FixedVectorNoTuple{2,Float64}
        x::Float64
        y::Float64
    end
    v = MyVec[MyVec(i,-i) for i=1:5]
    @fslice v[:x, :]
"""
macro fslice(expr)
    fixed_slice_expr(expr)
end


text:
using Behavior
using Behavior.Gherkin
using Niancat.Languages
using Niancat.Games.NiancatGames
using Niancat.Games.NiancatGames: Incorrect, Correct
using Niancat.Games
import Niancat.Gameface: score!

struct ByColumnTable
    table::Gherkin.DataTable
end

column(t::ByColumnTable, n::Int = 1) = [row[n] for row in t.table]

@given("a Swedish dictionary") do context
    table = ByColumnTable(context.datatable)
    words = column(table)
    context[:dictionary] = SwedishDictionary(words)
end

@given("a new game") do context
    dictionary = context[:dictionary]
    game = NiancatGame(dictionary, TestingGameService())
    context[:game] = game
end

@given("a puzzle {String}") do context, puzzle
    team = Team(1, "defaultteam", "")
    user = User(1, "name", "displayname", team)
    response = gamecommand(context[:game], user, SetPuzzle(puzzle))
    context[:response] = response
end

@when("a guess {String} is made") do context, word
    team = Team(1, "defaultteam", "")
    user = User(1, "name", "displayname", team)
    response = gamecommand(context[:game], user, Guess(word))
    context[:response] = response
end

@then("the response is that the guess is wrong") do context
    @expect context[:response] isa Incorrect
end

@then("the response is that the word is correct") do context
    @expect context[:response] isa Correct
end


text:
import ModelSanitizer
import Test

a = [1, 2, 3]
b = [4 5 6; 7 8 9; 10 11 12]
c = Array{Int, 3}(undef, 3, 3, 3)
c[:] .= 100

d = ["foo", "bar", "baz"]
e = ["foo" "bar" "baz"; "foo" "bar" "baz"]
f = Array{String, 3}(undef, 3, 3, 3)
f[:] .= "foo"

g = Any["foo", "bar", "baz"]
h = Any["foo" "bar" "baz"; "foo" "bar" "baz"]
i = Array{Any, 3}(undef, 3, 3, 3)
i[:] .= "foo"

j = Union{String, Missing}["foo", missing, "baz"]
k = Union{Missing, Nothing}[missing, nothing, missing]
l = Nothing[nothing, nothing, nothing]

Test.@test !any(a .== 0)
Test.@test !any(b .== 0)
Test.@test !any(c .== 0)
Test.@test !any(d .== "")
Test.@test !any(e .== "")
Test.@test !any(f .== "")
Test.@test !any(g .== 0)
Test.@test !any(h .== 0)
Test.@test !any(i .== 0)
Test.@test !any(skipmissing(j .== ""))
Test.@test !(skipmissing(k) == skipmissing([missing, nothing, missing]))
Test.@test all(l .== nothing)

ModelSanitizer.zero!(a)
ModelSanitizer.zero!(b)
ModelSanitizer.zero!(c)
ModelSanitizer.zero!(d)
ModelSanitizer.zero!(e)
ModelSanitizer.zero!(f)
ModelSanitizer.zero!(g)
ModelSanitizer.zero!(h)
ModelSanitizer.zero!(i)
ModelSanitizer.zero!(j)
ModelSanitizer.zero!(k)
ModelSanitizer.zero!(l)

Test.@test all(a .== 0)
Test.@test all(b .== 0)
Test.@test all(c .== 0)
Test.@test all(d .== "")
Test.@test all(e .== "")
Test.@test all(f .== "")
Test.@test all(g .== 0)
Test.@test all(h .== 0)
Test.@test all(i .== 0)
Test.@test all(j .== "")
Test.@test all(k .== nothing)
Test.@test all(l .== nothing)

struct Z end

Test.@test zero(Any) == 0
Test.@test zero(Z) == 0
Test.@test zero(String) == ""
Test.@test zero(Nothing) == nothing
Test.@test zero(Missing) === missing
Test.@test ismissing(zero(Missing))


text:
module AnalyticsService

using HTTP, Sockets, Dates

include("postgres.jl")
include("utils.jl")

function getPoints(req::HTTP.Request)::HTTP.Response
    result = getPgData()

    @info "Marshall time"
    resp = @time df2json(result)

    return HTTP.Response(200, resp)
    end;

const ROUTER = HTTP.Router()
HTTP.@register(ROUTER, "GET", "/api", getPoints)

@info "Run server"
HTTP.serve(ROUTER, Sockets.localhost, 8081)

end # module

max_stars_repo_name: kuznetsovin/julia-first-app

max_issues_repo_name: kuznetsovin/julia-first-app

max_forks_repo_name: kuznetsovin/julia-first-app


text:
using Test
using BlockSystems
using ModelingToolkit
using ModelingToolkit: get_iv, get_eqs, get_states
using Graphs

# ordering and simplification should not matter for equity of equations!
import Base.==
function ==(A::Vector{Equation}, B::Vector{Equation})
    Set(simplify.(A)) == Set(simplify.(B))
end

@info "Testes of transformations.jl"

@testset "transformations.jl" begin
    @testset "pairwise cycle free" begin
        using BlockSystems: _pairwise_cycle_free
        g = SimpleDiGraph(5)
        add_edge!(g, 1=>2)
        add_edge!(g, 2=>4)
        add_edge!(g, 4=>5)
        add_edge!(g, 3=>5)
        add_edge!(g, 3=>1)
        add_edge!(g, 2=>3)
        add_edge!(g, 5=>4)
        cycles = simplecycles(g)
        @test _pairwise_cycle_free(g) == [3,5]
    end

    @testset "remove superfluous states" begin
        using BlockSystems: remove_superfluous_states
        @parameters t a b i(t)
        @variables x(t) y(t) o(t) o1(t) o2(t)
        D = Differential(t)
        eqs = [D(x) ~ x,
               D(y) ~ y,
               o ~ a^b]
        reqs = remove_superfluous_states(IOBlock(eqs, [], [x]; iv=t)) |> equations
        @test reqs == eqs[[1]]
        reqs = remove_superfluous_states(IOBlock(eqs, [], [y]; iv=t)) |> equations
        @test reqs == eqs[[2]]
        reqs = remove_superfluous_states(IOBlock(eqs, [], [o]; iv=t)) |> equations
        @test reqs == eqs[[3]]
        eqs = [D(x) ~ y,
               D(y) ~ x+o,
               o ~ a^b]
        reqs = remove_superfluous_states(IOBlock(eqs, [], [y]; iv=t)) |> equations
        @test reqs == eqs
        reqs = remove_superfluous_states(IOBlock(eqs, [], [x]; iv=t)) |> equations
        @test reqs == eqs
        reqs = remove_superfluous_states(IOBlock(eqs, [], [o]; iv=t)) |> equations
        @test reqs == eqs[[3]]
        eqs = [D(x) ~ x+o,
               D(y) ~ y,
               o ~ a^b]
        reqs = remove_superfluous_states(IOBlock(eqs, [], [x]; iv=t)) |> equations
        @test reqs == eqs[[1,3]]
        reqs = remove_superfluous_states(IOBlock(eqs, [], [y]; iv=t)) |> equations
        @test reqs == eqs[[2]]
        reqs = remove_superfluous_states(IOBlock(eqs, [], [o]; iv=t)) |> equations
        @test reqs == eqs[[3]]
    end

    @testset "remove algebraic states" begin
        @parameters t a b i(t)
        @variables x(t) y(t) o(t) o1(t) o2(t)
        D = Differential(t)
        # variable o can be removed
        eqs = [D(x) ~ x + o,
               D(y) ~ y + o,
               o ~ a^b]
        reqs = substitute_algebraic_states(IOBlock(eqs,[],[x,y]))

        @test equations(reqs) == [D(x) ~ x + a^b,
                          D(y) ~ y + a^b]
        @test reqs.removed_eqs == [o ~ a^b]

        # variable o1 and o2 can be removed
        eqs = [D(x) ~ x + o1+o2,
               D(y) ~ y + o2,
               o1 ~ a^b + o2,
               o2 ~ a-b]
        reqs = substitute_algebraic_states(IOBlock(eqs,[],[x,y]))
        @test equations(reqs) == [D(x) ~ x + a^b + a-b + a-b,
                          D(y) ~ y + a-b]
        @test reqs.removed_eqs == [o1 ~ a^b + a - b,
                          o2 ~ a - b]

        eqs = [D(x) ~ i + o,
               o ~ x + i]
        reqs = substitute_algebraic_states(IOBlock(eqs,[],[x]))
        @test equations(reqs) == [D(x) ~ i + (x + i)]
        @test reqs.removed_eqs == [o ~ x + i]

        eqs = [D(x) ~ i,
               o1 ~ x + o2,
               D(y) ~ i,
               o2 ~ y + 2*o1]
        reqs = substitute_algebraic_states(IOBlock(eqs,[],[x,y]))
        @test equations(reqs) == [D(x) ~ i,
                          D(y) ~ i,
                          0 ~ x + (y+2*o1) - o1]
        @test reqs.removed_eqs == [o2 ~ y + 2*o1]

        rreqs = substitute_algebraic_states(IOBlock(equations(reqs), [],[x,y]))
        @test equations(rreqs) == equations(reqs)
        @test reqs.removed_eqs == reqs.removed_eqs

        # test skip condition
        eqs = [D(x) ~ i + o,
               o ~ x + i]
        reqs = substitute_algebraic_states(IOBlock(eqs,[],[x,o]))
        @test equations(reqs) == eqs
        @test reqs.removed_eqs == []

        eqs = [D(x) ~ i,
               o1 ~ x + o2,
               D(y) ~ i,
               o2 ~ y + 2*o1]
        reqs = substitute_algebraic_states(IOBlock(eqs,[],[x,y,o1]))
        @test equations(reqs) == [D(x) ~ i,
                          D(y) ~ i,
                          0 ~ x + (y + 2*o1) - o1]
        @test reqs.removed_eqs == [o2 ~ y + 2*o1]

        reqs = substitute_algebraic_states(IOBlock(eqs,[],[x,y,o2]))
        @test equations(reqs) == [D(x) ~ i,
                          D(y) ~ i,
                          0 ~ y + 2*(x + o2) - o2]
        @test reqs.removed_eqs == [o1 ~ x + o2]
    end

    @testset "test of connect_system" begin
        #=
                    +------------+
        in1 --> i1 -|    iob1    |
        in2 --> i2 -|(x1, x2)(a) |-o--+     +-----+
                    +------------+    +-ina-| add |- add ---> out
                    +------------+    +-inb-|     |
        in3 --> i1 -|    iob2    |-o--+     +-----+
        in4 --> i2 -|(x1, x2)(b) |
                    +------------+
        =#
        # same system as in testset before
        @parameters t i1(t) i2(t) a b ina(t) inb(t)
        @variables x1(t) x2(t) o(t) add(t)
        D = Differential(t)
        eqs1  = [D(x1) ~ a*i1, D(x2)~i2, o~x1+x2]
        iob1 = IOBlock(eqs1, [i1, i2], [o], name=:A)
        eqs2  = [D(x1) ~ b*i1, D(x2)~i2, o~x1+x2]
        iob2 = IOBlock(eqs2, [i1, i2], [o], name=:B)
        ioadd = IOBlock([add ~ ina + inb], [ina, inb], [add], name=:add)
        @parameters in1(t) in2(t) in3(t) in4(t)
        @variables out(t)
        sys = IOSystem([iob1.o => ioadd.ina, iob2.o => ioadd.inb],
                       [iob1, iob2, ioadd],
                       namespace_map = Dict(iob1.i1 => in1,
                                            iob1.i2 => in2,
                                            iob2.i1 => in3,
                                            iob2.i2 => in4,
                                            ioadd.add => out),
                       outputs = [out],
                       name=:sys)
        iob = connect_system(sys)

        @test Set(sys.inputs) == Set(iob.inputs)
        @test Set(sys.iparams) == Set(iob.iparams)
        # the connection will delete the states iob1.o and iob2.o
        @test Set(sys.istates) == Set(iob.istates) ∪ Set([iob1.o, iob2.o])
        @test Set(sys.outputs) == Set(iob.outputs)
        @test iob.name == sys.name

        @variables A₊x1(t) A₊x2(t) B₊x1(t) B₊x2(t) A₊o(t) B₊o(t)
        eqs = [D(A₊x1) ~ a * in1,
               D(A₊x2) ~ in2,
               D(B₊x1) ~ b * in3,
               D(B₊x2) ~ in4,
               out ~ (A₊x1 + A₊x2) + (B₊x1 + B₊x2)]
        @test eqs == get_eqs(iob.system)
        @test iob.removed_eqs == [A₊o ~ A₊x1 + A₊x2, B₊o ~ B₊x1 + B₊x2]

        @testset "test rename_vars" begin
            @test_throws ArgumentError new = rename_vars(iob, in2=:in1)
            new = rename_vars(iob, in2=:in2N, out=:outN, A₊x1=:y, a=:c)

            @parameters in2N(t) c
            @variables outN(t) y(t)
            eqs = [D(y) ~ c * in1,
                   D(A₊x2) ~ in2N,
                   D(B₊x1) ~ b * in3,
                   D(B₊x2) ~ in4,
                   outN ~ (y + A₊x2) + (B₊x1 + B₊x2)]
            @test eqs == get_eqs(new.system)
            @test new.removed_eqs == [A₊o ~ y + A₊x2, B₊o ~ B₊x1 + B₊x2]
        end
    end

    @testset "system with removed equations in subsystem" begin
        #=
           +------+         +---------+
        +->| x'=i |--+------| o=a1+a2 |--o
        |  +------+  |   +--|         |
        |  +------+  |   |  +---------+
        +--| o = i|<-+   |
           +------+      |
           +------+      |
        +->| x'=i |--+---+
        |  +------+  |
        |  +------+  |
        +--| o = i|<-+
           +------+
        =#
        @parameters t i(t)
        @variables x(t) o(t)
        D = Differential(t)

        b1 = IOBlock([D(x) ~ i], [i], [x], name=:a)
        b2 = IOBlock([o ~ i], [i], [o], name=:b)
        subsys = IOSystem([b2.o => b1.i, b1.x => b2.i], [b1, b2], name=:A,
                          namespace_map = [b1.x=>x],
                          outputs = [x])
        @test Set(subsys.removed_states) == Set()

        subsysA = connect_system(subsys)
        @test get_eqs(subsysA.system) == [D(x) ~ x]
        @test Set(subsysA.removed_states) == Set([o])
        @test subsysA.removed_eqs == [o ~ x]

        subsysB = IOBlock(subsysA, name=:B)
        @test subsysB.name == :B
        @test subsysB.system.name == :B
        @test get_eqs(subsysB.system) == [D(x) ~ x]
        @test Set(subsysB.removed_states) == Set([o])
        @test subsysB.removed_eqs == [o ~ x]

        @parameters a1(t) a2(t)
        adder = IOBlock([o ~ a1 + a2], [a1, a2], [o], name=:add)

        system = IOSystem([subsysA.x => adder.a1, subsysB.x => adder.a2],
                          [adder, subsysA, subsysB])

        @test Set(system.removed_states) == Set([subsysA.o, subsysB.o])

        systemblock = connect_system(system)
        @test systemblock.system.name == systemblock.name == system.name
        @test Set(systemblock.removed_states) == Set(system.removed_states)

        using BlockSystems: namespace_rem_eqs
        @test systemblock.removed_eqs == vcat(namespace_rem_eqs(subsysA), namespace_rem_eqs(subsysB))

        # psst, i'm putting a test for function_generation here, don't tell ma
        gen = generate_io_function(systemblock,
                                   f_states=[systemblock.A₊x, systemblock.B₊x],
                                   f_rem_states=[systemblock.A₊o, systemblock.B₊o],
                                   warn=false)
        @test Set(gen.states) == Set(Sym{Real}.([:B₊x, :A₊x, :add₊o]))
        @test Set(gen.inputs) == Set()
        @test Set(gen.params) == Set()
        @test Set(gen.rem_states) == Set(Sym{Real}.([:B₊o, :A₊o]))

        out = zeros(2)
        st = rand(3)
        gen.g_ip(out, st, (), (), (), 0.0)
        @test out == st[1:2]
        @test gen.g_oop(st, (), (), (), 0.0) == st[1:2]
    end

    @testset "substitution of differential states" begin
        # first block
        @parameters t i(t)
        @variables o(t)
        D = Differential(t)
        A = IOBlock([D(o) ~ 1 + D(i)], [i], [o])
        @test length(BlockSystems.rhs_differentials(A)) == 1

        # second block, does not contain derivative of output
        @variables x(t)
        @parameters a(t)
        B = IOBlock([x ~ a + 1], [], [x])
        @test length(BlockSystems.rhs_differentials(B)) == 0

        # third block, contains derivative of output
        @variables y(t)
        @parameters b(t)
        C = IOBlock([D(y) ~ b + 1], [], [y])
        @test length(BlockSystems.rhs_differentials(B)) == 0

        sysAB = IOSystem([B.x => A.i], [A, B]) |> connect_system
        @test length(BlockSystems.rhs_differentials(sysAB)) == 1

        sysAC = IOSystem([C.y => A.i], [A, C]) |> connect_system
        @test length(BlockSystems.rhs_differentials(sysAC)) == 0

        @test equations(sysAC.system) == [D(o) ~ 2 + b, D(y) ~ 1 + b]

        # fouth block with more complex stuff
        @variables x(t) y(t)
        B = IOBlock([D(x) ~ a, y ~ x + 3], [], [y])

        sysAB = IOSystem([B.y => A.i], [A, B]) |> connect_system
        @test length(BlockSystems.rhs_differentials(sysAB)) == 0
    end

    @testset "subs of differential states" begin
        @parameters t a(t) b
        @variables x(t) y(t) z(t)
        D = Differential(t)

        blk = IOBlock([D(x) ~ 1 + D(y),
                       D(y) ~ a + b], [], [])
        blk2 = substitute_derivatives(blk, verbose=true)
        @test equations(blk2) == [D(x) ~ 1 + a + b,
                                  D(y) ~ a + b]

        blk = IOBlock([D(x) ~ 1 + D(y),
                       y ~ 2 + b], [], [])
        blk2 = substitute_derivatives(blk, verbose=true)
        @test equations(blk2) == [D(x) ~ 1 ,
                                  y ~ 2 + b]

        blk = IOBlock([D(x) ~ 1 + D(y),
                       y ~ 2 + a], [], [])
        blk2 = substitute_derivatives(blk, verbose=true)
        @test equations(blk2) == [D(x) ~ 1 + D(a),
                                  y ~ 2 + a]

        blk = IOBlock([D(x) ~ 1 + D(y),
                       y ~ 2 + z,
                       D(z) ~ 5*x + b], [], [])
        blk2 = substitute_derivatives(blk, verbose=true)
        @test equations(blk2) == [D(x) ~ 1 + 5*x + b,
                                  y ~ 2 + z,
                                  D(z) ~ 5*x + b]
    end

    @testset "set p" begin
        @parameters t a(t) b
        @variables x(t) y(t) z(t)
        D = Differential(t)

        blk = IOBlock([D(x) ~ 1 + D(y) + a*b], [], [])

        blk2 = set_p(blk, :a=>0)
        @test equations(blk2) == [D(x) ~ 1 + D(y)]
        blk2 = set_p(blk, blk.a=>0)
        @test equations(blk2) == [D(x) ~ 1 + D(y)]

        blk2 = set_p(blk, :b=>1)
        @test equations(blk2) == [D(x) ~ 1 + D(y) + a]
        blk2 = set_p(blk, blk.b=>1)
        @test equations(blk2) == [D(x) ~ 1 + D(y) + a]

        blk2 = set_p(blk, Dict(blk.a=>2, b=>4))
        @test equations(blk2) == [D(x) ~ 9 + D(y)]

        blk2 = set_p(blk, blk.a=>2, b=>4; warn=false)
        @test equations(blk2) == [D(x) ~ 9 + D(y)]

        @test_throws ArgumentError blk2 = set_p(blk, :a=>:bla)
        @test_throws ArgumentError blk2 = set_p(blk, x=>2.0)
    end
end


text:
using BlockHankelMatrices
using Documenter

DocMeta.setdocmeta!(BlockHankelMatrices, :DocTestSetup, :(using BlockHankelMatrices); recursive=true)

makedocs(;
    modules=[BlockHankelMatrices],
    authors="Zdeněk Hurák <hurak@fel.cvut.cz>",
    repo="https://github.com/hurak/BlockHankelMatrices.jl/blob/{commit}{path}#{line}",
    sitename="BlockHankelMatrices.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://hurak.github.io/BlockHankelMatrices.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
)

deploydocs(;
    repo="github.com/hurak/BlockHankelMatrices.jl",
    devbranch="main",
)


text:
# function green_fn(sim::Simulation{1},ω::Number,x::Real)
#     inds = IrosBase.Tessellations.get_surrounding_sites(sim.tessellation,x,y)
#     weights = IrosBase.Simulations.generate_weights(x,y,sim.x[inds],sim.y[inds])
#     j = zeros(sim.n)
#     x = sim.x[inds]; y = sim.y[inds]
#     v1 = [x[2]-x[1],y[2]-y[1]]
#     v2 = [x[3]-x[2],y[3]-y[2]]
#     cosθ = (norm(v1+v2)^2 - norm(v1)^2 - norm(v2)^2)/2/norm(v1)/norm(v2)
#     sinθ = sqrt(1-cosθ^2)
#     area = norm(v1)*norm(v2)*sinθ/2
#     j[inds] = weights./area
#     return solve_scatter(sim,k,j)
# end

function scattering(
            sim::Simulation{1},
            ω::Number,
            j::AbstractArray;
            ky::Number=0,
            kz::Number=0)

    A,B = maxwell_lep(sim;ky=ky,kz=kz)(ω)
    return scattering_core(A-B*ω^2,j)
end
#
#
# function solve_scatter(sim::Simulation,kx::Number,ky::Number,shape::AbstractShape,f::Function)
#     k = hypot(kx,ky)
#     x = sim.x[sim.interior]
#     y = sim.y[sim.interior]
#     L = Laplacian(sim,k,0,0)
#     φ = (f.(kx,ky,x,y).*shape.(x,y))
#     j = (L.D+sparse(I*k^2,sim.n,sim.n))*φ
#     A = L.D + L.M*k^2
#     scattered_field = solve_scatter_core(A,-j)
#     total_field = scattered_field + φ
#     return total_field, scattered_field, φ
# end
# solve_scatter(sim::Simulation,k::Vector,shape::AbstractShape,f::Function) = solve_scatter(sim,k[1],k[2],shape,f)


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct SshPublicKey <: SwaggerModel
    path::Any # spec type: Union{ Nothing, String } # spec name: path
    keyData::Any # spec type: Union{ Nothing, String } # spec name: keyData

    function SshPublicKey(;path=nothing, keyData=nothing)
        o = new()
        validate_property(SshPublicKey, Symbol("path"), path)
        setfield!(o, Symbol("path"), path)
        validate_property(SshPublicKey, Symbol("keyData"), keyData)
        setfield!(o, Symbol("keyData"), keyData)
        o
    end
end # type SshPublicKey

const _property_map_SshPublicKey = Dict{Symbol,Symbol}(Symbol("path")=>Symbol("path"), Symbol("keyData")=>Symbol("keyData"))
const _property_types_SshPublicKey = Dict{Symbol,String}(Symbol("path")=>"String", Symbol("keyData")=>"String")
Base.propertynames(::Type{ SshPublicKey }) = collect(keys(_property_map_SshPublicKey))
Swagger.property_type(::Type{ SshPublicKey }, name::Symbol) = Union{Nothing,eval(Meta.parse(_property_types_SshPublicKey[name]))}
Swagger.field_name(::Type{ SshPublicKey }, property_name::Symbol) =  _property_map_SshPublicKey[property_name]

function check_required(o::SshPublicKey)
    true
end

function validate_property(::Type{ SshPublicKey }, name::Symbol, val)
end


text:
abstract QuODESolvers <: QuPropagatorMethod

const type_to_method_ode = @compat Dict{Any, Any}(:QuODE45 => ODE.ode45, :QuODE78 => ODE.ode78, :QuODE23s => ODE.ode23s)

for qu_ode_type in keys(type_to_method_ode)
    method_name = type_to_method_ode[qu_ode_type]
    @eval  begin
        @doc """
        ODE solver, using $($(qu_ode_type))

        Step Propagation using the ode solver $($(method_name)).

        ### Fields :

        * options :: Dict

        Dictionary to set the relative tolerance and absolute tolerance by using
        keys as `:reltol` and `:abstol`.
        """ ->
        immutable $qu_ode_type <: QuODESolvers
            options::Dict{Symbol, Any}
        end
        $qu_ode_type() = $qu_ode_type(Dict())
    end
end

for (qu_ode_type,ode_solver) in type_to_method_ode
    @eval begin
        function propagate(prob::$qu_ode_type, eq::QuEquation, t, current_t, current_qustate)
            op = operator(eq, t)
            dims = size(current_qustate)
            # Convert the current_qustate to complex as it might result in a Inexact Error. After complex is in QuBase.jl (PR #38)
            # we could just do a complex(vec(current_qustate)) avoiding the coeffs(coeffs(vec(current_qustate))).
            next_state = $ode_solver((t,y)-> -im*coeffs(op)*y, complex(coeffs(vec(current_qustate))), [current_t, t], points=:specified,
                                  reltol = get(prob.options, :reltol, 1.0e-5), abstol = get(prob.options, :abstol, 1.0e-8))[2][end]
            CQST = QuBase.similar_type(current_qustate)
            return CQST(reshape(next_state, dims), bases(current_qustate))
        end
    end
end

export QuODE45,
      QuODE78,
      QuODE23s


text:
"""
    @type_commutative 2 3 f(op, x::Type1, y::Type2) = ...

Wrap a method definition, resulting in the definition of another method
which is identical except that the types of x and y are swapped.
Note that, any expression valid for the type of x must be valid for the type of y,
since the body of the function stays exactly the same.
"""
macro type_commutative(ix, iy, expr)
    additional_decl = deepcopy(expr)

    fargs = extract_fargs(ix, iy, additional_decl)
    ix, iy = handle_ix_iy(ix, iy, fargs)

    if length(fargs[ix].args) == length(fargs[iy].args) == 1
        fargs[ix].args[1], fargs[iy].args[1] = fargs[iy].args[1], fargs[ix].args[1]
    else
        fargs[ix].args[2], fargs[iy].args[2] = fargs[iy].args[2], fargs[ix].args[2]
    end

    quote
        $(esc(expr))
        $(esc(additional_decl))
    end
end

macro type_commutative(expr) :(@type_commutative($nothing, $nothing, $expr)) end

"""
    @commutative 2 3 f(op, x::Type1, y::Type2) = ...

Wrap a method definition, resulting in the definition of another method
which is identical except that x and y are swapped.
"""
macro commutative(ix, iy, expr)
    additional_decl = deepcopy(expr)

    fargs = extract_fargs(ix, iy, additional_decl)
    ix, iy = handle_ix_iy(ix, iy, fargs)

    fargs[ix], fargs[iy] = fargs[iy], fargs[ix]

    quote
        $(esc(expr))
        $(esc(additional_decl))
    end
end

macro commutative(expr) :(@commutative($nothing, $nothing, $expr)) end

function extract_fargs(ix, iy, decl)
    fdecl, _ = decl.args
    if fdecl.head == :where
        fdecl = first(fdecl.args)
    end

    fargs = @view fdecl.args[2:end]
    @assert length(fargs) >= 2 "The method must have at least two arguments"

    fargs
end

function handle_ix_iy(ix, iy, fargs)
    if isnothing(ix) || isnothing(iy)
        @assert length(fargs) == 2 "`ix` and `iy` must be provided for a method with more than 2 arguments"
        ix = 1
        iy = 2
    end

    ix, iy
end

"""
Return `val` as a subscript, used for printing `UnitBlade` and `Blade` instances.
"""
function subscript(val)
    r = div(val, 10)
    subscript_char(x) = Char(8320 + x)
    r > 0 ? string(subscript_char(r), subscript_char(mod(val, 10))) : string(subscript_char(val))
end


text:
using OpenTelemetryProto: OpentelemetryClients
using Test

@testset "OpenTelemetryProto.jl" begin
    req =
        OpentelemetryClients.opentelemetry.proto.collector.trace.v1.ExportTraceServiceRequest(
            resource_spans = [
                OpentelemetryClients.opentelemetry.proto.trace.v1.ResourceSpans(
                    resource = OpentelemetryClients.opentelemetry.proto.resource.v1.Resource(
                        dropped_attributes_count = 0,
                    ),
                    instrumentation_library_spans = [
                        OpentelemetryClients.opentelemetry.proto.trace.v1.InstrumentationLibrarySpans(
                            instrumentation_library = OpentelemetryClients.opentelemetry.proto.common.v1.InstrumentationLibrary(;
                                name = "julia",
                                version = "1.0.0",
                            ),
                            spans = [
                                OpentelemetryClients.opentelemetry.proto.trace.v1.Span(
                                    trace_id = reinterpret(UInt8, rand(UInt128, 1)),
                                    span_id = reinterpret(UInt8, rand(UInt64, 1)),
                                    trace_state = "",
                                    parent_span_id = [],
                                    name = "JJJJJJJJJJJJJJJJJJ",
                                    kind = 0,
                                    start_time_unix_nano = 0xe3c6aac9ebe5b6f8,
                                    end_time_unix_nano = 0xe3c6aac9ebe5b6f9,
                                ),
                            ],
                            schema_url = "http://localhost:7777",
                        ),
                    ],
                    schema_url = "http://localhost:8888",
                ),
            ],
        )

    client = OpentelemetryClients.TraceServiceBlockingClient("http://localhost:4317")

    # OpentelemetryClients.Export(client, req)

end


text:
m = [1 2 3; 3 4 6]
dm = DataArray(m)
df = DataFrame(m)
@assert isequal(DataFrame([2.0 3.0 4.5;], colnames(df)), colmeans(df))
@assert isequal(DataFrame([4 6 9;], colnames(df)), colsums(df))
@assert isequal(DataArray([2.0, 3.0, 4.5]), colmeans(dm))


text:
# # Helmholtz decomposition

#md # ```@meta
#md # CurrentModule = ImmersedLayers
#md # ```

#=
For vector fields on the grid, we often prefer to work with potentials
of that field and compose the vector field from these potentials. This
is the so-called Helmholtz decomposition

$$\mathbf{v} = \nabla\phi + \nabla\times\psi$$

where $\psi$ is a vector potential (or streamfunction in two dimensions)
and $\phi$ is a scalar potential. Note that, in an unbounded domain, $\phi$
satisfies a Poisson equation:

$$\nabla^2\phi = \nabla\cdot\mathbf{v}$$

and, assuming that $\psi$ has zero 'gauge' (i.e., its divergence is zero), then
$\psi$ also satisfies a Poisson equation:

$$\nabla^2\psi = -\nabla\times\mathbf{v}$$

So we can recover the velocity field from the divergence and curl of the
vector field (in fluid dynamics, these are the rate of dilatation and the vorticity,
respectively) by solving these equations for the potentials and then
evaluating their gradient and curl.

With immersed layers, these are extendable to domains with surfaces
and their associated jump in the vector field $[\mathbf{v}] = \mathbf{v}^+ - \mathbf{v}^-$.
The crucial relationships are

$$\nabla\cdot{\overline{\mathbf{v}}} = \overline{\nabla\cdot\mathbf{v}} + \delta(\chi)\mathbf{n}\cdot[\mathbf{v}]$$

and

$$\nabla\times{\overline{\mathbf{v}}} = \overline{\nabla\times\mathbf{v}} + \delta(\chi)\mathbf{n}\times[\mathbf{v}]$$

These are, respectively, the divergence and curl of the masked vector field, $\overline{\mathbf{v}}$.
In `ImmersedLayers.jl`, we are always working with the masked form of a field,
whether we acknowledge it or not. In contrast, the first terms on the right-hand sides of these are the masked divergence and curl
of the vector field in each side of the surface. The difference between these is
a type of single-layer surface term---a 'source sheet' or a 'vortex sheet', respectively---
due to the jump in vector field. If we omit the overbar notation, then the Helmholtz
decomposition and the Poisson equations all still hold.

This package has a number of routines that can be used to perform the
operations described here. They rely on two caches, [`ScalarPotentialCache`](@ref)
and [`VectorPotentialCache`](@ref), each of which supports one of the
respective potentials. The function [`vecfield_helmholtz!`](@ref) carries
out all of these operations in order to assemble the vector field $\mathbf{v}$
from the masked divergence and curl fields and the jump in the vector field.
(We can also add an additional irrotational, divergence-free vector field during
this assembly.
=#

#md # ## Helmholtz decomposition functions

#md # ```@docs
#md # ScalarPotentialCache
#md # VectorPotentialCache
#md # vectorpotential_from_masked_curlv!
#md # scalarpotential_from_masked_divv!
#md # vectorpotential_from_curlv!
#md # vecfield_from_vectorpotential!
#md # masked_curlv_from_curlv_masked!
#md # curlv_masked_from_masked_curlv!
#md # scalarpotential_from_divv!
#md # masked_divv_from_divv_masked!
#md # divv_masked_from_masked_divv!
#md # vecfield_from_scalarpotential!
#md # vecfield_helmholtz!
#md # vectorpotential_uniformvecfield!
#md # scalarpotential_uniformvecfield!
#md # vecfield_uniformvecfield!
#md # ```


text:
using JuliaCon2020RANSAC
using RANSAC
using StaticArrays
using FileIO

using RANSACVisualizer
using Makie

fname = joinpath(@__DIR__, "models", "furcsa32k.obj")

m = load(fname)
pc = RANSACCloud(m.position, m.normals, 4)

## Without translational
defpars = ransacparameters([FittedSphere, FittedPlane, FittedCylinder])
planep = (ϵ = 0.2, α = deg2rad(4),)
spherep = (ϵ = 0.3, α = deg2rad(1),)
cylinderp = (ϵ = 0.4, α = deg2rad(7),)
conep = (ϵ = 0.2, α = deg2rad(1),)
iterp = (minsubsetN=15, itermax=10000, τ=1200, prob_det=0.99,)
params = ransacparameters(defpars, iteration=iterp, plane=planep, sphere=spherep, cylinder=cylinderp, cone=conep,)

extr2, t = ransac(pc, params, true; reset_rand=true)
showshapes(pc, extr2, show_axis = false)
sc = showbytype(pc, extr2, show_axis = false, resolution=(1920,1080))
sc.center = false
save("no_transl.png", sc)

## With translational
defpars = ransacparameters([FittedSphere, FittedPlane, FittedCylinder, FittedTranslational])
translp = (ϵ = 0.3, α = deg2rad(3), α_perpend = cosd(89.5), thinning_par=0.3/4, force_transl=false, thin_method=:slow, min_normal_num=0.99,)

params = ransacparameters(defpars, iteration=iterp, plane=planep, sphere=spherep, cylinder=cylinderp, cone=conep, translational=translp,)

extr, t = ransac(pc, params, true; reset_rand=true)

showshapes(pc, extr)
sc2 = showbytype(pc, extr, show_axis = false, resolution=(1920,1080))
sc2.center = false
save("transl.png", sc2)


text:
function destroy!(itr::AttributesIterator)
    if itr.first_ber_element != C_NULL
        ber_free(itr.first_ber_element, 0)
        itr.first_ber_element = C_NULL
    end
    nothing
end

function AttributesIterator(entry::Entry)
    first_ber_element_ref = Ref{Ptr{Cvoid}}()
    cstr = ldap_first_attribute(entry.chain.ldap.handle, entry.handle, first_ber_element_ref)

    local str::Union{Nothing, String}
    if cstr == C_NULL
        str = nothing
    else
        str = unsafe_string(cstr)
        ldap_memfree(cstr)
    end

    return AttributesIterator(entry, first_ber_element_ref[], str)
end

each_attribute(entry::Entry) = AttributesIterator(entry)

function Base.collect(itr::AttributesIterator)
    result = Vector{Attribute}()
    for a in itr
        push!(result, a)
    end
    return result
end

struct Attribute
    name::String
    vals::Vector{String}
end

function Base.getindex(entry::Entry, attribute_name::AbstractString)
    vals_ptr_vec = ldap_get_values(entry.chain.ldap.handle, entry.handle, attribute_name)
    vals_string_vec = parse_null_terminated_string_vector(vals_ptr_vec)
    ldap_value_free(vals_ptr_vec)
    return vals_string_vec
end

Attribute(entry::Entry, attribute_name::AbstractString) = Attribute(attribute_name, entry[attribute_name])

function Base.iterate(itr::AttributesIterator)
    if itr.first_ber_element == C_NULL
        @assert itr.first_attribute_name == nothing
        return nothing
    else
        return Attribute(itr.entry, itr.first_attribute_name), itr.first_ber_element
    end
end

function Base.iterate(itr::AttributesIterator, state::Ptr{Cvoid})
    next_attribute_cstr = ldap_next_attribute(itr.entry.chain.ldap.handle, itr.entry.handle, state)
    if next_attribute_cstr == C_NULL
        return nothing
    else
        str = unsafe_string(next_attribute_cstr)
        ldap_memfree(next_attribute_cstr)
        return Attribute(itr.entry, str), itr.first_ber_element
    end
end