text:
module Benchmarks

export run_performance_exps, show_performance_exps, visualize_benchmark
export run_variation_exps, show_variation_exps
export default_setting, variation_settings

using OneVision
using OneVision.Examples

using ProgressMeter
using DataFrames
using ThreadsX
using DrWatson
using CSV
using Plots
using Random
using Statistics
using Measurements
using Printf

function with_args(f; kws...)
    (args...; kws2...) -> f(args...; kws2..., kws...)
end

scenarios = [
    "1D Leader Linear" => with_args(Car1DExample.run_example,
        use_bang_bang = false, has_obstacle = false)
    "1D Leader With Obstacle" => with_args(Car1DExample.run_example,
        use_bang_bang = true, has_obstacle = true)
    "2D Formation Driving" => with_args(Car2DExamples.formation_example,
        switch_formation = false)
    "2D Formation Switching" => with_args(Car2DExamples.formation_example,
        switch_formation = true)
]

CFs = [
    :Naive => naive_cf,
    :Local => local_cf,
    :ConstU => const_u_cf,
    :OneVision => onevision_cf,
]

default_setting() = defaults

defaults = ExampleSetting(
    time_end = 20.0,
    freq = 100.0,
    noise = 0.05 / sqrt(100),
    sensor_noise = 0.05 / sqrt(100),
    delays = DelayModel(obs = 3, act = 4, com = 5, ΔT = 5),
    H = 20,
    model_error = 0.0,
)

variation_settings() = [
    "disturbance" =>
        [x => @set defaults.noise *= x
        for x in [0.0, 1.0, 3.0, 6.0, 10.0]]
    "noise" =>
        [x => @set defaults.sensor_noise *= x
        for x in [0.0, 1.0, 3.0, 6.0, 10.0]]
    "com delay" =>
        [x => @set defaults.delays.com = x
        for x in [1, 5, 15, 30, 50]]
    "horizon" =>
        [x => @set defaults.H = x
        for x in [1, 2, 4, 7, 10, 15, 20, 30]]
    "model error" =>
        [x => @set defaults.model_error = x
        for x in [0.0, 0.1, 0.3, 0.6, 1.0]]
]

basic_settings() = [
    "Default" => defaults
    "Larger Delays" =>  @set defaults.delays.com *= 4
    "No Delays" =>  @set defaults.delays = DelayModel(;obs = 0, act = 0, com = 1, ΔT = 5)
]

function run_setting(setting, seed)
    rows = []
    for (name, ex) in scenarios
        foreach(CFs) do (cf_str, CF)
            loss, metrics = ex(;plot_result = false, CF, setting, seed)
            metrics["loss"] = sum(loss) / setting.freq
            metrics["loss-state"] = sum(loss[:, :, :x]) / setting.freq
            for (metric, v) in metrics
                row = (task = name, cf = cf_str, metric = metric, value = v)
                push!(rows, row)
            end
        end
    end
    rows
end

"""
Rearrange the result DataFrame for easy comparison of the performance of
different CFs
"""
function rearrange_result_table(table::DataFrame)
    rows = []
    for df in groupby(table, :Task)
        task = df[1, :Task]
        cols = Symbol.(df[:, :CF]) .=> df[:, :value]
        best = sortperm(cols, by = x -> x[2])[1]

        push!(rows, (Task = task, cols..., Best = df[best, :CF]))
    end
    DataFrame(rows)
end

function display_table(rows, name)
    sep = " "^displaysize(stdout)[2]
    println("\r$sep")
    println("==== $name ====")
    display(rows)
end

"""
Generates `num` random seeds that are unlikely to be correlated.
"""
function independent_seeds(num; seed = 1)
    rng = MersenneTwister(seed)
    (1:num) .+ rand(rng, 1:num*100, num)
end

function run_performance_exps(; num_of_runs)
    OneVision.TrajPlanning.check_optimizer_converge[] = false
    seeds = independent_seeds(num_of_runs)

    @info "Running performance benchmarks ($num_of_runs tasks)..."
    prog = Progress(num_of_runs, 0.1, "Running benchmarks...")

    iolock = ReentrantLock()
    all_rows = []
    @time Threads.@threads for seed in seeds
        results = run_setting(defaults, seed)
        lock(iolock) do
            for result in results
                @unpack task, cf, metric, value = result
                push!(all_rows, @ntuple(task, cf, metric, seed, value))
            end
        end
        next!(prog)
    end

    return DataFrame(all_rows)
end

"""
Convert "loss" into "loss sqrt" and drop "loss-state".
Rename "avg deviation" to "avg distance".
"""
function prepare_metrics(df::AbstractDataFrame)
    df = filter(r -> r.metric != "loss-state", df)
    rows = df.metric .== "loss"
    df[rows, :value] = log10.(df[rows, :value])
    df[rows, :metric] .= Ref("log loss")
    df[df.metric .== "avg deviation", :metric] .= "avg distance"
    # df[df.metric .== "avg distance", :metric] .= "avg deviation"
    df
end

function show_performance_exps(df, out_path = nothing)
    comp_stats(vs) = mean(vs) # ± std(vs)
    display_num(v) = if v isa AbstractFloat; @sprintf("%.3f", v) else v end

    (out_path === nothing) || mkpath(out_path)
    df = prepare_metrics(df)
    df = groupby(df, Not([:seed, :value]))
    df = combine(df, :value => comp_stats => :value)
    cf_symbols = first.(CFs)
    for subdf in groupby(df, :metric)
        metric = subdf.metric[1]
        if metric == "loss-state" continue end
        println("=== Metric: $(metric) ===")
        table = @_ groupby(subdf, :task) |>
            combine(__, [:cf, :value] => ((cfs, vs) -> (; zip(cfs, vs)...)) => AsTable)
        display(table)
        (out_path === nothing) ||
            CSV.write(joinpath(out_path, "$(metric).csv"), table;
                transform = (_, val) -> display_num(val))
    end
end

function visualize_benchmark(scenario::String, setting::ExampleSetting, cf::String; seed = 1)
    ex = Dict(scenarios)[scenario]
    CF = Dict(CFs)[Symbol(cf)]
    ex(;plot_result = true, CF, setting, seed)
    nothing
end

function run_variation_exps(; num_of_runs)
    OneVision.TrajPlanning.check_optimizer_converge[] = false
    seeds = independent_seeds(num_of_runs)
    all_tasks = [
        @ntuple variation variable seed setting
        for seed in seeds  # this prevents different seeds from running in parallel and helps precompilation
        for (variation, tasks) in variation_settings()
        for (variable, setting) in tasks
    ]
    shuffle!(all_tasks)

    @info "Running benchmarks ($(length(all_tasks)) tasks)..."
    prog = Progress(length(all_tasks), 0.1, "Running benchmarks...")

    time_used = []
    iolock = ReentrantLock()
    all_rows = []
    @time Threads.@threads for task in all_tasks
        @unpack variation, variable, seed, setting = task
        stats = @timed run_setting(setting, seed)
        lock(iolock) do
            time = stats.time
            push!(time_used, @ntuple variation variable seed time)
            for result in stats.value
                @unpack task, cf, metric, value = result
                r = @ntuple variation variable task cf metric seed value
                push!(all_rows, r)
            end
        end
        next!(prog)
    end

    @_ DataFrame(time_used) |>
        groupby(__, Not(:seed)) |>
        combine(__, :time => mean => :avg_time) |>
        sort!(__, [:avg_time], rev = true) |>
        display

    sort!(all_rows)
    return DataFrame(all_rows)
end

function show_variation_exps(df::DataFrame; out_path = nothing, only_task = nothing)
    function comp_stats(vs)
        # quart1, quart2, quart3 = quantile(vs, [0.25, 0.5, 0.75])
        (mid = mean(vs), lb = minimum(vs), ub = maximum(vs))
    end

    data = prepare_metrics(df)
    data = groupby(data, Not([:seed, :value]))
    data = combine(data, :value => comp_stats => [:mid, :lb, :ub])
    (out_path === nothing) || mkpath(out_path)
    @showprogress map(d -> draw_variation_metric(d; out_path, only_task), groupby(data, [:variation, :metric], sort = true))
    nothing
end

function draw_variation_metric(subdata; out_path = nothing, only_task = nothing)
    var_name = subdata[1, :variation]
    metric_name = subdata[1, :metric]
    # yscale = startswith(metric_name, "loss") ? :log10 : :identity
    cf_names = hcat(first.(CFs)...)
    plots = []
    for taskdata in groupby(subdata, :task)
        task_name = taskdata.task[1]
        only_task === nothing || only_task == task_name || continue
        p = plot(;xlabel = var_name, ylabel = metric_name)
        for cfname in cf_names
            cfdata = filter(x -> x.cf == cfname, taskdata)
            plot!(cfdata.variable, cfdata.mid;
                title = task_name, #yscale = yscale,
                marker = true, markersize = 3,
                ribbon = (cfdata.mid .- cfdata.lb, cfdata.ub .- cfdata.mid))
        end
        push!(plots, p)
    end
    N = length(plots)
    plt = plot(plots...;
        labels=string.(cf_names), legend = :outertopright, fillalpha = 0.3,
        fontsize = 10, legendfontsize = 6, titlefontsize = 10,
        dpi=300,
        layout = (N, 1), size = (400, 250 * N))
    display(plt)
    if out_path !== nothing
        path = joinpath(out_path, "$(var_name)-$(metric_name).png")
        savefig(plt, path)
    end
end

end # module Benchmarks


text:
export Signal
export data
export fps
export seconds
export chunks
export chunk
export dft
export fft

import WAV

"""
analyse(signal)

Converts a signal from the time domain to the frequency domain
with the given signal.
"""
function analyse(signal)
    N = size(signal, 1)
    f(ω) = begin
        scale = 1 / sqrt(N)
        estimate = 0
        for t = 0:N-1
            estimate += signal[t+1] * exp(-im * 2π * ω / N * t)
        end
        scale * estimate
    end
    Ω = Array{Complex}(undef, N)
    for ω = 0:N-1
        Ω[ω+1] = f(ω)
    end
    Ω
end

"""
synthesize(signal)

Converts a signal from the frequency domain to the time domain
with the given signal.
"""
function synthesize(signal)
    N = size(signal, 1)
    f(t) = begin
        scale = 1 / sqrt(N)
        estimate = 0
        for ω = 0:N-1
            estimate += signal[ω+1] * exp(im * 2π * ω / N * t)
        end
        scale * estimate
    end
    T = Array{Complex}(undef, N)
    for t = 0:N-1
        D[t+1] = f(t)
    end
    T
end

struct Signal
    data::Array{Float64}
    fps::Float64
end

Signal(s::String; c=1) = begin
    y, fps = WAV.wavread(s)
    y = (y .+ 1) ./ 2
    Signal(y[:, c], fps)
end

data(s::Signal) = s.data

fps(s::Signal) = s.fps

seconds(s::Signal) = Int(size(data(s), 1) ÷ fps(s))

chunks(s::Signal, cps::Int) = seconds(s) * cps

"""
chunk(s::Signal, i::Int, cps::Int)

Gets a chunk of a signal such that there are a certain number of chunks per
second with the given signal, the chunk ordinal number and the number of chunks
per second.
"""
function chunk(s::Signal, i::Int, cps::Int)
    samples_per_chunk = Int(fps(s) ÷ cps)
    window = 2^11
    start = (i - 1) * samples_per_chunk + 1
    finish = start + window
    total_samples = size(data(s), 1)
    if finish < total_samples
        return Signal(data(s)[start:finish-1], fps(s))
    else
        return Signal(data(s)[total_samples-window+1:end], fps(s))
    end
end

function fft(s::Signal)
    x = data(s)
    N = size(x, 1)
    if N % 2 > 0
        println("Must be a power of 2.")
        return analyse(x)
    elseif N ≤ 2
        return analyse(x)
    else
        e = fft(Signal(view(x, 2:2:N), fps(s)))
        o = fft(Signal(view(x, 1:2:N-1), fps(s)))
        r = convert(Array{Int64}, floor.(range(0, stop=N-1, length=N)))
        v = exp.(-im .* 2pi .* r ./ N)
        i = Integer(N ÷ 2)
        return [e .+ v[1:i] .* o; e .+ v[i+1:end] .* o]
    end
end

# Discrete Fourier transform
dft(s::Signal) = analyse(data(s))

ext: jl


text:
module mbObserve
export Observer, observe, unobserve, notify

type Observer
	D::Dict{Any, Vector{Function}}
	Observer() = new( Dict{Any, Vector{Function}}() )
end

function observe(cb::Function, O::Observer, s)
	!haskey(O.D, s) && (O.D[s] = Function[])
	push!(O.D[s], cb)
	return Void
end

function unobserve(cb::Function, O::Observer, s)
	try
		delete!(O.D[s], cb)
		isempty(O.D[s]) && delete!(O.D, s)
	end
	return Void
end

#NOTE not exposed... conceptually dangerous... responsibility for callbacks from different
#	clients can only lie at observer. however clients should know their state, which is
#	not given if Observer can remove their callbacks.
function unobserveall(O::Observer, s)
	O.D[s] = Function[]
	return Void
end

#TODO DRYS
#TODO try... used this way for debugging

function notify(O::Observer, s, v...)
	for cb in get(O.D, s, [])
		# cb(v...)
		try cb(v...) catch e; println(e, "\nSymbol: $s.") end
	end
	return Void
end

function notify(O::Observer, s)
	for cb in get(O.D, s, [])
		# cb()
		try cb() catch e; println(e, "\nSymbol: $s.") end
	end
	return Void
end

end


text:
const templ = joinpath(splitdir(dirname(@__FILE__))[1], "scripts", "_template.jl")

const list_file = joinpath(dirname(@__FILE__), "data", "list.txt")
const file_one  = joinpath(dirname(@__FILE__), "data", "one.txt")
const file_two  = joinpath(dirname(@__FILE__), "data", "two.txt")

print("""

Test scripts
============
""")

print("""

STDIN -> STDOUT
---------------
""")

let out = readall( pipeline(`julia -e 'print(trues(2,2))'`, `julia $templ`) )

    @test ismatch(r"RUN : 0", out)
    @test length(matchall(r"true", out)) == 4
end

@unix_only begin

    print("""
    --arg
    """)

    let out = readall( pipeline(`cat $list_file`, `julia $templ --arg 42`) )

        @test ismatch(r"RUN : 42", out)
        @test length(matchall(r"\.txt", out)) == 2
    end

    print("""
    --list & -a
    """)

    let out = readall( pipeline(`cat $list_file`, `julia $templ --a 42 --list`) )

        @test ismatch(r"ONE", out) # Printed into STDOUT
        @test ismatch(r"TWO", out) # Printed into STDOUT
        @test !ismatch(r"RUN : 42", out) # Printed into the file
        @test filesize("one.mitos.txt") != 0 # Created file
        @test filesize("two.mitos.txt") != 0 # Created file
        @test ismatch(r"RUN : 42", readall("one.mitos.txt")) # Printed into the file
        @test ismatch(r"RUN : 42", readall("two.mitos.txt")) # Printed into the file

        # Clean up

        if isfile("one.mitos.txt")
            rm("one.mitos.txt")
        end

        if isfile("two.mitos.txt")
            rm("two.mitos.txt")
        end

    end

end

print("""

File -> File
------------
""")

print("""
    --list & -a
    """)

let out = readall( `julia $templ $list_file --a 42 --list` )

    @test ismatch(r"ONE", out) # Printed into STDOUT
    @test ismatch(r"TWO", out) # Printed into STDOUT
    @test !ismatch(r"RUN : 42", out) # Printed into the file
    @test filesize("one.mitos.txt") != 0 # Created file
    @test filesize("two.mitos.txt") != 0 # Created file
    @test ismatch(r"RUN : 42", readall("one.mitos.txt")) # Printed into the file
    @test ismatch(r"RUN : 42", readall("two.mitos.txt")) # Printed into the file

    # Clean up

    if isfile("one.mitos.txt")
        rm("one.mitos.txt")
    end

    if isfile("two.mitos.txt")
        rm("two.mitos.txt")
    end

end

print("""

-p (parallel)
-------------
""")

print("""
File -> File (--list)
""")

let out = readall( `julia $templ $list_file -p 2 --a 42 --list` )

    @test ismatch(r"ONE", out) # Printed into STDOUT
    @test ismatch(r"TWO", out) # Printed into STDOUT
    @test !ismatch(r"RUN : 42", out) # Printed into the file
    @test filesize("one.mitos.txt") != 0 # Created file
    @test filesize("two.mitos.txt") != 0 # Created file
    @test ismatch(r"RUN : 42", readall("one.mitos.txt")) # Printed into the file
    @test ismatch(r"RUN : 42", readall("two.mitos.txt")) # Printed into the file

    # Clean up

    if isfile("one.mitos.txt")
        rm("one.mitos.txt")
    end

    if isfile("two.mitos.txt")
        rm("two.mitos.txt")
    end

end

@unix_only begin

    print("""
    STDIN -> File (--list)
    """)

    let out = readall( pipeline(`cat $list_file`, `julia $templ -p 2 --a 42 --list` ) )

        @test ismatch(r"ONE", out) # Printed into STDOUT
        @test ismatch(r"TWO", out) # Printed into STDOUT
        @test !ismatch(r"RUN : 42", out) # Printed into the file
        @test filesize("one.mitos.txt") != 0 # Created file
        @test filesize("two.mitos.txt") != 0 # Created file
        @test ismatch(r"RUN : 42", readall("one.mitos.txt")) # Printed into the file
        @test ismatch(r"RUN : 42", readall("two.mitos.txt")) # Printed into the file

        # Clean up

        if isfile("one.mitos.txt")
            rm("one.mitos.txt")
        end

        if isfile("two.mitos.txt")
            rm("two.mitos.txt")
        end

    end

end

print("""

Distances.jl
============
""")

let path_script = joinpath(splitdir(dirname(@__FILE__))[1], "scripts", "Distances.jl"),
    path_file   = joinpath(dirname(@__FILE__), "data", "small.pdb"),
    out_intra  = readall(`julia $path_script $path_file -o STDOUT`),
    out_inter  = readall(`julia $path_script $path_file --inter -o STDOUT`)

    @test sum(diff(Float64[ parse(Float64,num) for num in  matchall(r"(\d\.\d+)\n",out_intra) ])) == 0.0
    @test length(matchall(r"\n1,A,", out_intra)) == 1
    @test length(matchall(r"\n1,B,", out_intra)) == 1

    @test sum(Float64[ parse(Float64,num) for num in  matchall(r"(\d\.\d+)\n",out_inter) ]) == 4.0*3.7836265671971385
    @test length(matchall(r"\n1,A,", out_inter)) == 5
    @test length(matchall(r"\n1,B,", out_inter)) == 1

end


text:
"""
    convolve(d1::Distribution, d2::Distribution)

Convolve two distributions and return the distribution corresponding to the sum of
independent random variables drawn from the underlying distributions.

Currently, the function is only defined in cases where the convolution has a closed form.
More precisely, the function is defined if the distributions of `d1` and `d2` are the same
and one of
* [`Bernoulli`](@ref)
* [`Binomial`](@ref)
* [`NegativeBinomial`](@ref)
* [`Geometric`](@ref)
* [`Poisson`](@ref)
* [`Normal`](@ref)
* [`Cauchy`](@ref)
* [`Chisq`](@ref)
* [`Exponential`](@ref)
* [`Gamma`](@ref)
* [`MvNormal`](@ref)

External links: [List of convolutions of probability distributions on Wikipedia](https://en.wikipedia.org/wiki/List_of_convolutions_of_probability_distributions)
"""
convolve(::Distribution, ::Distribution)

# discrete univariate
function convolve(d1::Bernoulli, d2::Bernoulli)
    _check_convolution_args(d1.p, d2.p)
    return Binomial(2, d1.p)
end

function convolve(d1::Binomial, d2::Binomial)
    _check_convolution_args(d1.p, d2.p)
    return Binomial(d1.n + d2.n, d1.p)
end

function convolve(d1::NegativeBinomial, d2::NegativeBinomial)
    _check_convolution_args(d1.p, d2.p)
    return NegativeBinomial(d1.r + d2.r, d1.p)
end

function convolve(d1::Geometric, d2::Geometric)
    _check_convolution_args(d1.p, d2.p)
    return NegativeBinomial(2, d1.p)
end

convolve(d1::Poisson, d2::Poisson) =  Poisson(d1.λ + d2.λ)

# continuous univariate
convolve(d1::Normal, d2::Normal) = Normal(d1.μ + d2.μ, hypot(d1.σ, d2.σ))
convolve(d1::Cauchy, d2::Cauchy) = Cauchy(d1.μ + d2.μ, d1.σ + d2.σ)
convolve(d1::Chisq, d2::Chisq) = Chisq(d1.ν + d2.ν)

function convolve(d1::Exponential, d2::Exponential)
    _check_convolution_args(d1.θ, d2.θ)
    return Gamma(2, d1.θ)
end

function convolve(d1::Gamma, d2::Gamma)
    _check_convolution_args(d1.θ, d2.θ)
    return Gamma(d1.α + d2.α, d1.θ)
end

# continuous multivariate
function convolve(d1::MvNormal, d2::MvNormal)
    _check_convolution_shape(d1, d2)
    return MvNormal(d1.μ + d2.μ, d1.Σ + d2.Σ)
end

function _check_convolution_args(p1, p2)
    p1 ≈ p2 || throw(ArgumentError(
    "$(p1) !≈ $(p2): distribution parameters must be approximately equal",
    ))
end

function _check_convolution_shape(d1, d2)
    length(d1) == length(d2) || throw(ArgumentError("$d1 and $d2 are not the same size"))
end


text:
similar(A::DEDataArray) = deepcopy(A)

done(A::DEDataArray, i::Int) = done(A.x,i)
eachindex(A::DEDataArray)      = eachindex(A.x)
next(A::DEDataArray, i::Int) = next(A.x,i)
start(A::DEDataArray)          = start(A.x)

length(A::DEDataArray) = length(A.x)
ndims(A::DEDataArray)  = ndims(A.x)
size(A::DEDataArray)   = size(A.x)

# Maybe should use fieldtype(typeof(dest), $i) ?
@generated function recursivecopy!{T<:DEDataArray}(dest::T, src::T)
   assignments = [:(typeof(getfield(dest,$i)) <: AbstractArray ? recursivecopy!(getfield(dest, $i),getfield(src, $i)) : setfield!(dest, $i, getfield(src, $i))) for i=1:nfields(T)]
   :($(assignments...); dest)
end

@generated function copy_non_array_fields{T<:DEDataArray}(previous::T, arr::AbstractArray)
  assignments = [:(getfield(previous,$i)) for i=2:nfields(T)]
  :(T(arr,$(assignments...)))
end

getindex( A::DEDataArray,    i::Int) = (A.x[i])
setindex!(A::DEDataArray, x, i::Int) = (A.x[i] = x)
getindex( A::DEDataArray,    i::Int...) = (A.x[i...])
setindex!(A::DEDataArray, x, i::Int...) = (A.x[i...] = x)
@compat Base.IndexStyle(::Type{DEDataArray}) = Base.IndexLinear()


text:
using ImageMagick, ColorTypes, FixedPointNumbers, IndirectArrays, FileIO, OffsetArrays
using ImageShow       # for show(io, ::MIME, img) & ImageMeta
using Test
using ImageCore

ontravis = haskey(ENV, "TRAVIS")

mutable struct TestType end

@testset "IO" begin
    workdir = joinpath(tempdir(), "Images")
    isdir(workdir) && rm(workdir, recursive=true)
    mkdir(workdir)

    a = rand(Bool,5,5,5,5)
    fn = joinpath(workdir, "5by5.png")
    @test_throws ErrorException ImageMagick.save(fn, a)

    a = [TestType() TestType()]
    fn = joinpath(workdir, "5by5.png")
    errfile, io = mktemp()  # suppress warning message
    redirect_stderr(io) do
        @test_throws MethodError ImageMagick.save(fn, a)
    end
    close(io)
    rm(errfile)

    @testset "Binary png" begin
        a = rand(Bool,5,5)
        fn = joinpath(workdir, "5by5.png")
        ImageMagick.save(fn, a)
        b = ImageMagick.load(fn)
        a8 = convert(Array{Gray{N0f8}}, a) # IM won't read back as Bool
        @test b == a8
        aim = colorview(Gray, a)
        ImageMagick.save(fn, aim)
        b = ImageMagick.load(fn)
        @test b == a8
        a = bitrand(5,5)
        fn = joinpath(workdir, "5by5.png")
        ImageMagick.save(fn, a)
        b = ImageMagick.load(fn)
        a8 = convert(Array{Gray{N0f8}}, a)
        @test b == a8
        aim = colorview(Gray, a)
        ImageMagick.save(fn, aim)
        b = ImageMagick.load(fn)
        @test b == a8

        @test ImageMagick.metadata(fn) == ((5,5), Gray{N0f8})
    end

    @testset "Gray png" begin
        a = [0 1/2^16 1/2^8; 1-1/2^8 1-1/2^16 1]
        aa = convert(Array{N0f8}, a)
        fn = joinpath(workdir, "2by3.png")
        ImageMagick.save(fn, a)
        b = ImageMagick.load(fn)
        @test b == aa
        ImageMagick.save(fn, aa)
        b = ImageMagick.load(fn)
        @test b == aa
        open(fn, "w") do io
            show(io, MIME("image/png"), b; minpixels=0)
        end
        bb = ImageMagick.load(fn)
        @test bb == b
        aaimg = Gray.(aa)
        ImageMagick.save(fn, aaimg)
        b = ImageMagick.load(fn)
        @test b == aaimg
        aa = convert(Array{N0f16}, a)
        ImageMagick.save(fn, aa)
        b = ImageMagick.load(fn)
        @test eltype(eltype(b)) == N0f16
        # Gets loaded as RGB{N0f16} on windows/osx
        @test N0f16.(Gray.(b)) == aa
        m = ImageMagick.metadata(fn)
        @test m[1]==(3,2)
    end

    @testset "Color" begin
        fn = joinpath(workdir, "2by2.png")
        A = rand(3,2,2)
        A[1] = 1
        img = colorview(RGB, A)
        img24 = convert(Array{RGB24}, img)
        ImageMagick.save(fn, img24)
        b = ImageMagick.load(fn)
        imgrgb8 = convert(Array{RGB{N0f8}}, img)
        @test imgrgb8 == b

        open(fn, "w") do io
            show(io, MIME("image/png"), imgrgb8; minpixels=0)
        end
        bb = ImageMagick.load(fn)
        @test bb == imgrgb8
    end

    @testset "Colormap usage" begin
        datafloat = reshape(range(0.5, stop=1.5, length=6), 2, 3)
        dataint = round.([UInt8], 254*(datafloat .- 0.5) .+ 1)  # ranges from 1 to 255
        # build our colormap
        b = RGB(0,0,1)
        w = RGB(1,1,1)
        r = RGB(1,0,0)
        cmaprgb = Array{RGB{Float64}}(undef, 255)
        f = range(0, stop=1, length=128)
        cmaprgb[1:128] = [(1-x)*b + x*w for x in f]
        cmaprgb[129:end] = [(1-x)*w + x*r for x in f[2:end]]
        img = IndirectArray(dataint, cmaprgb)
        ImageMagick.save(joinpath(workdir,"cmap.jpg"), img)
        cmaprgb = Array{RGB}(undef, 255) # poorly-typed cmap, Images issue #336
        cmaprgb[1:128] = [(1-x)*b + x*w for x in f]
        cmaprgb[129:end] = [(1-x)*w + x*r for x in f[2:end]]
        img = IndirectArray(dataint, cmaprgb)
        ImageMagick.save(joinpath(workdir, "cmap.png"), img)
    end

    @testset "Alpha" begin
        c = reinterpret(BGRA{N0f8}, [0xf0884422]'')
        fn = joinpath(workdir, "alpha.png")
    	ImageMagick.save(fn, c)
        C = ImageMagick.load(fn)
        if !ontravis || !Sys.islinux()
            # disabled on Linux Travis because it has a weird copy of
            # ImageMagick for which this fails (see Images#261)
            @test C[1] == c[1]
        end
        ImageMagick.save(fn, reinterpret(ARGB32, [0xf0884422]''))
        D = ImageMagick.load(fn)
        if !ontravis || !Sys.islinux()
            @test D[1] == c[1]
        end

        # Images#396
        c = colorview(RGBA, normedview(permuteddimsview(reshape(0x00:0x11:0xff, 2, 2, 4), (3,1,2))))
        ImageMagick.save(fn, c)
        D = ImageMagick.load(fn)
        if !ontravis || !Sys.islinux()
            @test D == c
        end
    end

    @testset "3D TIFF (issue #307)" begin
        Ar = rand(0x00:0xff, 2, 2, 4)
        Ar[1] = 0xff
        A = map(x->Gray(N0f8(x,0)), Ar)
        fn = joinpath(workdir, "3d.tif")
        ImageMagick.save(fn, A)
        B = ImageMagick.load(fn)

        @test A == B

        @test ImageMagick.metadata(fn) == ((2,2,4), Gray{N0f8})
    end

    @testset "16-bit TIFF (issue #49)" begin
        Ar = rand(0x0000:0xffff, 2, 2, 4)
        Ar[1] = 0xffff
        A = map(x->Gray(reinterpret(N0f16, x)), Ar)
        fn = joinpath(workdir, "3d16.tif")
        ImageMagick.save(fn, Ar)
        ImageMagick.save(fn, A)
        B = ImageMagick.load(fn)

        @test A == B
    end

#= FAILS ON SCIENTIFIC LINUX 7.2 WITH IMAGEMAGICK 6.9.5, works for other combos
    @testset "32-bit TIFF (issue #49)" begin
        Ar = rand(0x00000000:0xffffffff, 2, 2, 4)
        Ar[1] = 0xffffffff
        A = map(x->Gray(reinterpret(N0f32, x)), Ar)
        fn = joinpath(workdir, "3d32.tif")
        ImageMagick.save(fn, A)
        B = ImageMagick.load(fn)

        @test A == B
    end
=#

    @testset "Clamping (issue #256)" begin
        Ar = rand(2,2)
        Ar[1] = 1
        A = colorview(Gray, Ar)
        A[1,1] = -0.4
        fn = joinpath(workdir, "2by2.png")
        errfile, io = mktemp()
        with_logger(SimpleLogger(io)) do
            @test_throws ArgumentError ImageMagick.save(fn, A)
        end
        close(io)
        @test occursin("out-of-range", readlines(errfile)[1])
        rm(errfile)
        ImageMagick.save(fn, A, mapi=clamp01nan)
        B = ImageMagick.load(fn)
        A[1,1] = 0
        @test B == map(Gray{N0f8}, A)
    end

    Sys.isunix() && @testset "Reading from a stream (issue #312)" begin
        fn = joinpath(workdir, "2by2.png")
        img = open(fn) do io
            ImageMagick.load(io)
        end
        @test size(img) == (2,2)
    end

    Sys.isunix() && @testset "Writing to a stream (PR #22)" begin
        orig_img = ImageMagick.load(joinpath(workdir, "2by2.png"))
        fn = joinpath(workdir, "2by2_fromstream.png")
        open(fn, "w") do f
            ImageMagick.save(Stream(format"PNG", f), orig_img)
        end
        img = ImageMagick.load(fn)
        @test img == orig_img
    end

    Sys.isunix() && @testset "Reading from a byte array (issue #279)" begin
        fn = joinpath(workdir, "2by2.png")
        io = open(fn)
        arr = read(io)
        close(io)
        img = readblob(arr)
        @test size(img) == (2,2)
    end

    @testset "show(MIME)" begin
        Ar = rand(UInt8, 3, 2, 2)
        Ar[1] = typemax(eltype(Ar))
        a = colorview(RGB, normedview(Ar))
        fn = joinpath(workdir, "2by2.png")
        open(fn, "w") do file
            show(file, MIME("image/png"), a, minpixels=0)
        end
        b = ImageMagick.load(fn)
        @test b == a

        Ar = rand(UInt8, 3, 1021, 1026)
        Ar[1] = typemax(eltype(Ar))
        abig = colorview(RGB, normedview(Ar))
        fn = joinpath(workdir, "big.png")
        open(fn, "w") do file
            show(file, MIME("image/png"), abig, maxpixels=10^6)
        end
        b = ImageMagick.load(fn)
        @test b == convert(Array{RGB{N0f8},2}, restrict(abig, (1,2)))

        # Issue #269
        Ar = rand(UInt16, 3, 1024, 1023)
        Ar[1] = typemax(eltype(Ar))
        abig = colorview(RGB, normedview(N0f16, Ar))
        open(fn, "w") do file
            show(file, MIME("image/png"), abig, maxpixels=10^6)
        end
        b = ImageMagick.load(fn)
        @test b == convert(Array{RGB{N0f8},2}, restrict(abig, (1,2)))
    end

    @testset "fps" begin
        A = rand(RGB{N0f8}, 100, 100, 5)
        fn = joinpath(workdir, "animated.gif")
        ImageMagick.save(fn, A, fps=2)
        wand = MagickWand()
        readimage(wand, fn)
        resetiterator(wand)
        @test ImageMagick.getimagedelay(wand) == 50
    end

    @testset "ImageMeta" begin
        # https://github.com/sisl/PGFPlots.jl/issues/5
        img = ImageMeta(rand(RGB{N0f8},3,5))
        fn = joinpath(workdir, "imagemeta.png")
        ImageMagick.save(fn, img)
        imgr = ImageMagick.load(fn)
        @test imgr == img
    end

    @testset "OffsetArrays" begin
        img = OffsetArray([true false; false true], 0:1, 3:4)
        fn = joinpath(workdir, "indices.png")
        ImageMagick.save(fn, img)
        imgr = ImageMagick.load(fn)
        @test imgr == parent(img)
    end
end

nothing


text:
@testset "Options" begin

    # Test Options
    opts = Options()
    @test typeof(opts) <: Options

end


text:
struct SimpleCellUpdater
end

function extensive2primitive(c::Conserved,
                             volume::Float64,
                             eos)
    density = c.mass/volume
    velocity = c.momentum/c.mass
    total_specific_energy = c.energy/c.mass
    kinetic_specific_energy = 0.5*velocity^2
    thermal_specific_energy = total_specific_energy-kinetic_specific_energy
    pressure = calcPressure(eos, density, thermal_specific_energy)
    return PrimitiveVars(density, pressure, velocity)
end

function updateCells!(cu::SimpleCellUpdater,
                      state::HydroSnapshot,
                      extensives::Array{Conserved, 1},
                      pg,
                      eos)
    cumulative_volumes = [calcVolume(pg, r) for r in state.grid]
    volumes = diff(cumulative_volumes)
    state.cells = [extensive2primitive(e,v,eos)
                    for (e,v) in zip(extensives, volumes)]
end

max_stars_repo_name: bolverk/Sirion

max_issues_repo_name: bolverk/Sirion

max_forks_repo_name: bolverk/Sirion


text:
using Colors

cmap = [
    RGB(0.267004, 0.004874, 0.329415),
    RGB(0.268510, 0.009605, 0.335427),
    RGB(0.269944, 0.014625, 0.341379),
    RGB(0.271305, 0.019942, 0.347269),
    RGB(0.272594, 0.025563, 0.353093),
    RGB(0.273809, 0.031497, 0.358853),
    RGB(0.274952, 0.037752, 0.364543),
    RGB(0.276022, 0.044167, 0.370164),
    RGB(0.277018, 0.050344, 0.375715),
    RGB(0.277941, 0.056324, 0.381191),
    RGB(0.278791, 0.062145, 0.386592),
    RGB(0.279566, 0.067836, 0.391917),
    RGB(0.280267, 0.073417, 0.397163),
    RGB(0.280894, 0.078907, 0.402329),
    RGB(0.281446, 0.084320, 0.407414),
    RGB(0.281924, 0.089666, 0.412415),
    RGB(0.282327, 0.094955, 0.417331),
    RGB(0.282656, 0.100196, 0.422160),
    RGB(0.282910, 0.105393, 0.426902),
    RGB(0.283091, 0.110553, 0.431554),
    RGB(0.283197, 0.115680, 0.436115),
    RGB(0.283229, 0.120777, 0.440584),
    RGB(0.283187, 0.125848, 0.444960),
    RGB(0.283072, 0.130895, 0.449241),
    RGB(0.282884, 0.135920, 0.453427),
    RGB(0.282623, 0.140926, 0.457517),
    RGB(0.282290, 0.145912, 0.461510),
    RGB(0.281887, 0.150881, 0.465405),
    RGB(0.281412, 0.155834, 0.469201),
    RGB(0.280868, 0.160771, 0.472899),
    RGB(0.280255, 0.165693, 0.476498),
    RGB(0.279574, 0.170599, 0.479997),
    RGB(0.278826, 0.175490, 0.483397),
    RGB(0.278012, 0.180367, 0.486697),
    RGB(0.277134, 0.185228, 0.489898),
    RGB(0.276194, 0.190074, 0.493001),
    RGB(0.275191, 0.194905, 0.496005),
    RGB(0.274128, 0.199721, 0.498911),
    RGB(0.273006, 0.204520, 0.501721),
    RGB(0.271828, 0.209303, 0.504434),
    RGB(0.270595, 0.214069, 0.507052),
    RGB(0.269308, 0.218818, 0.509577),
    RGB(0.267968, 0.223549, 0.512008),
    RGB(0.266580, 0.228262, 0.514349),
    RGB(0.265145, 0.232956, 0.516599),
    RGB(0.263663, 0.237631, 0.518762),
    RGB(0.262138, 0.242286, 0.520837),
    RGB(0.260571, 0.246922, 0.522828),
    RGB(0.258965, 0.251537, 0.524736),
    RGB(0.257322, 0.256130, 0.526563),
    RGB(0.255645, 0.260703, 0.528312),
    RGB(0.253935, 0.265254, 0.529983),
    RGB(0.252194, 0.269783, 0.531579),
    RGB(0.250425, 0.274290, 0.533103),
    RGB(0.248629, 0.278775, 0.534556),
    RGB(0.246811, 0.283237, 0.535941),
    RGB(0.244972, 0.287675, 0.537260),
    RGB(0.243113, 0.292092, 0.538516),
    RGB(0.241237, 0.296485, 0.539709),
    RGB(0.239346, 0.300855, 0.540844),
    RGB(0.237441, 0.305202, 0.541921),
    RGB(0.235526, 0.309527, 0.542944),
    RGB(0.233603, 0.313828, 0.543914),
    RGB(0.231674, 0.318106, 0.544834),
    RGB(0.229739, 0.322361, 0.545706),
    RGB(0.227802, 0.326594, 0.546532),
    RGB(0.225863, 0.330805, 0.547314),
    RGB(0.223925, 0.334994, 0.548053),
    RGB(0.221989, 0.339161, 0.548752),
    RGB(0.220057, 0.343307, 0.549413),
    RGB(0.218130, 0.347432, 0.550038),
    RGB(0.216210, 0.351535, 0.550627),
    RGB(0.214298, 0.355619, 0.551184),
    RGB(0.212395, 0.359683, 0.551710),
    RGB(0.210503, 0.363727, 0.552206),
    RGB(0.208623, 0.367752, 0.552675),
    RGB(0.206756, 0.371758, 0.553117),
    RGB(0.204903, 0.375746, 0.553533),
    RGB(0.203063, 0.379716, 0.553925),
    RGB(0.201239, 0.383670, 0.554294),
    RGB(0.199430, 0.387607, 0.554642),
    RGB(0.197636, 0.391528, 0.554969),
    RGB(0.195860, 0.395433, 0.555276),
    RGB(0.194100, 0.399323, 0.555565),
    RGB(0.192357, 0.403199, 0.555836),
    RGB(0.190631, 0.407061, 0.556089),
    RGB(0.188923, 0.410910, 0.556326),
    RGB(0.187231, 0.414746, 0.556547),
    RGB(0.185556, 0.418570, 0.556753),
    RGB(0.183898, 0.422383, 0.556944),
    RGB(0.182256, 0.426184, 0.557120),
    RGB(0.180629, 0.429975, 0.557282),
    RGB(0.179019, 0.433756, 0.557430),
    RGB(0.177423, 0.437527, 0.557565),
    RGB(0.175841, 0.441290, 0.557685),
    RGB(0.174274, 0.445044, 0.557792),
    RGB(0.172719, 0.448791, 0.557885),
    RGB(0.171176, 0.452530, 0.557965),
    RGB(0.169646, 0.456262, 0.558030),
    RGB(0.168126, 0.459988, 0.558082),
    RGB(0.166617, 0.463708, 0.558119),
    RGB(0.165117, 0.467423, 0.558141),
    RGB(0.163625, 0.471133, 0.558148),
    RGB(0.162142, 0.474838, 0.558140),
    RGB(0.160665, 0.478540, 0.558115),
    RGB(0.159194, 0.482237, 0.558073),
    RGB(0.157729, 0.485932, 0.558013),
    RGB(0.156270, 0.489624, 0.557936),
    RGB(0.154815, 0.493313, 0.557840),
    RGB(0.153364, 0.497000, 0.557724),
    RGB(0.151918, 0.500685, 0.557587),
    RGB(0.150476, 0.504369, 0.557430),
    RGB(0.149039, 0.508051, 0.557250),
    RGB(0.147607, 0.511733, 0.557049),
    RGB(0.146180, 0.515413, 0.556823),
    RGB(0.144759, 0.519093, 0.556572),
    RGB(0.143343, 0.522773, 0.556295),
    RGB(0.141935, 0.526453, 0.555991),
    RGB(0.140536, 0.530132, 0.555659),
    RGB(0.139147, 0.533812, 0.555298),
    RGB(0.137770, 0.537492, 0.554906),
    RGB(0.136408, 0.541173, 0.554483),
    RGB(0.135066, 0.544853, 0.554029),
    RGB(0.133743, 0.548535, 0.553541),
    RGB(0.132444, 0.552216, 0.553018),
    RGB(0.131172, 0.555899, 0.552459),
    RGB(0.129933, 0.559582, 0.551864),
    RGB(0.128729, 0.563265, 0.551229),
    RGB(0.127568, 0.566949, 0.550556),
    RGB(0.126453, 0.570633, 0.549841),
    RGB(0.125394, 0.574318, 0.549086),
    RGB(0.124395, 0.578002, 0.548287),
    RGB(0.123463, 0.581687, 0.547445),
    RGB(0.122606, 0.585371, 0.546557),
    RGB(0.121831, 0.589055, 0.545623),
    RGB(0.121148, 0.592739, 0.544641),
    RGB(0.120565, 0.596422, 0.543611),
    RGB(0.120092, 0.600104, 0.542530),
    RGB(0.119738, 0.603785, 0.541400),
    RGB(0.119512, 0.607464, 0.540218),
    RGB(0.119423, 0.611141, 0.538982),
    RGB(0.119483, 0.614817, 0.537692),
    RGB(0.119699, 0.618490, 0.536347),
    RGB(0.120081, 0.622161, 0.534946),
    RGB(0.120638, 0.625828, 0.533488),
    RGB(0.121380, 0.629492, 0.531973),
    RGB(0.122312, 0.633153, 0.530398),
    RGB(0.123444, 0.636809, 0.528763),
    RGB(0.124780, 0.640461, 0.527068),
    RGB(0.126326, 0.644107, 0.525311),
    RGB(0.128087, 0.647749, 0.523491),
    RGB(0.130067, 0.651384, 0.521608),
    RGB(0.132268, 0.655014, 0.519661),
    RGB(0.134692, 0.658636, 0.517649),
    RGB(0.137339, 0.662252, 0.515571),
    RGB(0.140210, 0.665859, 0.513427),
    RGB(0.143303, 0.669459, 0.511215),
    RGB(0.146616, 0.673050, 0.508936),
    RGB(0.150148, 0.676631, 0.506589),
    RGB(0.153894, 0.680203, 0.504172),
    RGB(0.157851, 0.683765, 0.501686),
    RGB(0.162016, 0.687316, 0.499129),
    RGB(0.166383, 0.690856, 0.496502),
    RGB(0.170948, 0.694384, 0.493803),
    RGB(0.175707, 0.697900, 0.491033),
    RGB(0.180653, 0.701402, 0.488189),
    RGB(0.185783, 0.704891, 0.485273),
    RGB(0.191090, 0.708366, 0.482284),
    RGB(0.196571, 0.711827, 0.479221),
    RGB(0.202219, 0.715272, 0.476084),
    RGB(0.208030, 0.718701, 0.472873),
    RGB(0.214000, 0.722114, 0.469588),
    RGB(0.220124, 0.725509, 0.466226),
    RGB(0.226397, 0.728888, 0.462789),
    RGB(0.232815, 0.732247, 0.459277),
    RGB(0.239374, 0.735588, 0.455688),
    RGB(0.246070, 0.738910, 0.452024),
    RGB(0.252899, 0.742211, 0.448284),
    RGB(0.259857, 0.745492, 0.444467),
    RGB(0.266941, 0.748751, 0.440573),
    RGB(0.274149, 0.751988, 0.436601),
    RGB(0.281477, 0.755203, 0.432552),
    RGB(0.288921, 0.758394, 0.428426),
    RGB(0.296479, 0.761561, 0.424223),
    RGB(0.304148, 0.764704, 0.419943),
    RGB(0.311925, 0.767822, 0.415586),
    RGB(0.319809, 0.770914, 0.411152),
    RGB(0.327796, 0.773980, 0.406640),
    RGB(0.335885, 0.777018, 0.402049),
    RGB(0.344074, 0.780029, 0.397381),
    RGB(0.352360, 0.783011, 0.392636),
    RGB(0.360741, 0.785964, 0.387814),
    RGB(0.369214, 0.788888, 0.382914),
    RGB(0.377779, 0.791781, 0.377939),
    RGB(0.386433, 0.794644, 0.372886),
    RGB(0.395174, 0.797475, 0.367757),
    RGB(0.404001, 0.800275, 0.362552),
    RGB(0.412913, 0.803041, 0.357269),
    RGB(0.421908, 0.805774, 0.351910),
    RGB(0.430983, 0.808473, 0.346476),
    RGB(0.440137, 0.811138, 0.340967),
    RGB(0.449368, 0.813768, 0.335384),
    RGB(0.458674, 0.816363, 0.329727),
    RGB(0.468053, 0.818921, 0.323998),
    RGB(0.477504, 0.821444, 0.318195),
    RGB(0.487026, 0.823929, 0.312321),
    RGB(0.496615, 0.826376, 0.306377),
    RGB(0.506271, 0.828786, 0.300362),
    RGB(0.515992, 0.831158, 0.294279),
    RGB(0.525776, 0.833491, 0.288127),
    RGB(0.535621, 0.835785, 0.281908),
    RGB(0.545524, 0.838039, 0.275626),
    RGB(0.555484, 0.840254, 0.269281),
    RGB(0.565498, 0.842430, 0.262877),
    RGB(0.575563, 0.844566, 0.256415),
    RGB(0.585678, 0.846661, 0.249897),
    RGB(0.595839, 0.848717, 0.243329),
    RGB(0.606045, 0.850733, 0.236712),
    RGB(0.616293, 0.852709, 0.230052),
    RGB(0.626579, 0.854645, 0.223353),
    RGB(0.636902, 0.856542, 0.216620),
    RGB(0.647257, 0.858400, 0.209861),
    RGB(0.657642, 0.860219, 0.203082),
    RGB(0.668054, 0.861999, 0.196293),
    RGB(0.678489, 0.863742, 0.189503),
    RGB(0.688944, 0.865448, 0.182725),
    RGB(0.699415, 0.867117, 0.175971),
    RGB(0.709898, 0.868751, 0.169257),
    RGB(0.720391, 0.870350, 0.162603),
    RGB(0.730889, 0.871916, 0.156029),
    RGB(0.741388, 0.873449, 0.149561),
    RGB(0.751884, 0.874951, 0.143228),
    RGB(0.762373, 0.876424, 0.137064),
    RGB(0.772852, 0.877868, 0.131109),
    RGB(0.783315, 0.879285, 0.125405),
    RGB(0.793760, 0.880678, 0.120005),
    RGB(0.804182, 0.882046, 0.114965),
    RGB(0.814576, 0.883393, 0.110347),
    RGB(0.824940, 0.884720, 0.106217),
    RGB(0.835270, 0.886029, 0.102646),
    RGB(0.845561, 0.887322, 0.099702),
    RGB(0.855810, 0.888601, 0.097452),
    RGB(0.866013, 0.889868, 0.095953),
    RGB(0.876168, 0.891125, 0.095250),
    RGB(0.886271, 0.892374, 0.095374),
    RGB(0.896320, 0.893616, 0.096335),
    RGB(0.906311, 0.894855, 0.098125),
    RGB(0.916242, 0.896091, 0.100717),
    RGB(0.926106, 0.897330, 0.104071),
    RGB(0.935904, 0.898570, 0.108131),
    RGB(0.945636, 0.899815, 0.112838),
    RGB(0.955300, 0.901065, 0.118128),
    RGB(0.964894, 0.902323, 0.123941),
    RGB(0.974417, 0.903590, 0.130215),
    RGB(0.983868, 0.904867, 0.136897),
    RGB(0.993248, 0.906157, 0.143936)
]

function viridis()
    return cmap
end


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

@testset "skipchars for IOStream" begin
    mktemp() do path, file
        function append_to_file(str)
            mark(file)
            print(file, str)
            flush(file)
            reset(file)
        end
        # test it doesn't error on eof
        @test eof(skipchars(isspace, file))

        # test it correctly skips
        append_to_file("    ")
        @test eof(skipchars(isspace, file))

        # test it correctly detects comment lines
        append_to_file("#    \n   ")
        @test eof(skipchars(isspace, file, linecomment='#'))

        # test it stops at the appropriate time
        append_to_file("   not a space")
        @test !eof(skipchars(isspace, file))
        @test read(file, Char) == 'n'

        # test it correctly ignores the contents of comment lines
        append_to_file("  #not a space \n   not a space")
        @test !eof(skipchars(isspace, file, linecomment='#'))
        @test read(file, Char) == 'n'

        # test it correctly handles unicode
        for (byte,char) in zip(1:4, ('@','߷','࿊','𐋺'))
            append_to_file("abcdef$char")
            @test Base.codelen(char) == byte
            @test !eof(skipchars(isalpha, file))
            @test read(file, Char) == char
        end
    end
end

@testset "readbytes_some! via readbytes!" begin
    mktemp() do path, file
        function append_to_file(str)
            mark(file)
            print(file, str)
            flush(file)
            reset(file)
        end
        append_to_file("aaaaaaaaaaaaaaaaa")
        b = UInt8[0]
        readbytes!(file, b, all=false)
        @test String(b) == "a"
        # with resizing of b
        b = UInt8[]
        readbytes!(file, b, 2, all=false)
        @test String(b) == "aa"
        # readbytes_all with resizing
        b = UInt8[]
        readbytes!(file, b, 15)
        @test String(b) == "aaaaaaaaaaaaaa"
    end
end
@testset "issue #18755" begin
    mktemp() do path, io
        write(io, zeros(UInt8, 131073))
        @test position(io) == 131073
        write(io, zeros(UInt8, 131073))
        @test position(io) == 262146
    end
end

@test Base.open_flags(read=false, write=true, append=false) == (read=false, write=true, create=true, truncate=true, append=false)

max_stars_repo_name: djsegal/julia-fork

max_issues_repo_name: djsegal/julia-fork

max_forks_repo_name: djsegal/julia-fork


text:
using DataDeps

register(DataDep(
  "chess-king-rook-vs-king-pawn",
  "http://archive.ics.uci.edu/ml/datasets/Chess+%28King-Rook+vs.+King-Pawn%29",
  "http://archive.ics.uci.edu/ml/machine-learning-databases/chess/king-rook-vs-king-pawn/kr-vs-kp.data",
  "70a5540902b4063d3977d3e820ac5b9fe5624eb5d93709eee5d50d7a94fca205",
  post_fetch_method=(path -> begin
    UCIData.process_dataset(path,
        target_index=37,
        feature_indices=1:36,
        categoric_indices=1:36,
        footerskip=2,
    )
  end),
))

max_stars_repo_name: JackDunnNZ/uci-data

max_issues_repo_name: JackDunnNZ/uci-data

max_forks_repo_name: JackDunnNZ/uci-data


text:
type Shader
    ptr::Ptr{Void}

    function Shader(ptr::Ptr{Void})
        s = new(ptr)
        # finalizer(s, destroy)
        # s
    end
end

function Shader(vertex_shader::AbstractString, frag_shader::AbstractString)
    vert = isempty(vertex_shader) ? C_NULL : vertex_shader
    frag = isempty(frag_shader) ? C_NULL : frag_shader
    Shader(ccall((:sfShader_createFromFile, libcsfml_graphics), Ptr{Void}, (Ptr{Cchar}, Ptr{Cchar},), vert, frag))
end
Shader() = Shader("", "")

function ShaderFromMemory(vertex_shader::AbstractString, frag_shader::AbstractString)
    vert = isempty(vertex_shader) ? C_NULL : vertex_shader
    frag = isempty(frag_shader) ? C_NULL : frag_shader
    Shader(ccall((:sfShader_createFromMemory, libcsfml_graphics), Ptr{Void}, (Ptr{Cchar}, Ptr{Cchar},), vert, frag))
end

function Shader(shader::AbstractString)
    shadername = basename(shader)
    if shadername[search(shadername, ".")+1] == "v"
        # This shader is a vertex shader
        Shader(shader, "")
    elseif shadername[search(shadername, ".") + 1] == "f"
        # This shader is a fragment shader
        Shader("", shader)
    else
        println("$shadername is not a valid shader file")
    end
end

function destroy(shader::Shader)
    ccall((:sfShader_destroy, libcsfml_graphics), Void, (Ptr{Void},), shader.ptr)
end

function set_parameter(shader::Shader, name::AbstractString, x::Real)
    ccall((:sfShader_setFloatParameter, libcsfml_graphics), Void, (Ptr{Void}, Ptr{Cchar}, Cfloat,), shader.ptr, name, x)
end

function set_parameter(shader::Shader, name::AbstractString, x::Real, y::Real)
    ccall((:sfShader_setFloat2Parameter, libcsfml_graphics), Void, (Ptr{Void}, Ptr{Cchar}, Cfloat, Cfloat,), shader.ptr, name, x, y)
end

function set_parameter(shader::Shader, name::AbstractString, x::Real, y::Real, z::Real)
    ccall((:sfShader_setFloat3Parameter, libcsfml_graphics), Void, (Ptr{Void}, Ptr{Cchar}, Cfloat, Cfloat, Cfloat,), shader.ptr, name, x, y, z)
end

function set_parameter(shader::Shader, name::AbstractString, x::Real, y::Real, z::Real, w::Real)
    ccall((:sfShader_setFloat4Parameter, libcsfml_graphics), Void, (Ptr{Void}, Ptr{Cchar}, Cfloat, Cfloat, Cfloat, Cfloat), shader.ptr, name, x, y, z, w)
end

function set_parameter(shader::Shader, name::AbstractString, vector::Vector2f)
    ccall((:sfShader_setVector2Parameter, libcsfml_graphics), Void, (Ptr{Void}, Ptr{Cchar}, Vector2f,), shader.ptr, name, vector)
end

function set_parameter(shader::Shader, name::AbstractString, vector::Vector3f)
    ccall((:sfShader_setVector3Parameter, libcsfml_graphics), Void, (Ptr{Void}, Ptr{Cchar}, Vector3f), shader.ptr, name, vector)
end

function set_parameter(shader::Shader, name::AbstractString, color::Color)
    ccall((:sfShader_setColorParameter, libcsfml_graphics), Void, (Ptr{Void}, Ptr{Cchar}, Color), shader.ptr, name, color)
end

function set_parameter(shader::Shader, name::AbstractString, texture::Texture)
    ccall((:sfShader_setTextureParameter, libcsfml_graphics), Void, (Ptr{Void}, Ptr{Cchar}, Ptr{Void},), shader.ptr, name, texture.ptr)
end

function set_parameter(shader::Shader, name::AbstractString)
    ccall((:sfShader_setCurrentTextureParameter, libcsfml_graphics), Void, (Ptr{Void}, Ptr{Cchar},), shader.ptr, name)
end

function bind(shader::Shader)
    ccall((:sfShader_bind, libcsfml_graphics), Void, (Ptr{Void},), shader.ptr)
end

function unbind()
    ccall((:sfShader_bind, libcsfml_graphics), Void, (Ptr{Void},), C_NULL)
end

function shader_is_available()
    ccall((:sfShader_isAvailable, libcsfml_graphics), Bool, ())
end


text:
using LinearAlgebraicRepresentation
Lar = LinearAlgebraicRepresentation

function minimal_2cycles(V::Lar.Points, EV::Lar.ChainOp)

    function edge_angle(v::Int, e::Int)
        edge = EV[e, :]
        v2 = setdiff(edge.nzind, [v])[1]
        x, y = V[v2, :] - V[v, :]
        return atan(y, x)
    end

    for i in 1:EV.m
        j = min(EV[i,:].nzind...)
        EV[i, j] = -1
    end
    VE = convert(Lar.ChainOp, SparseArrays.transpose(EV))
    EF = Lar.Arrangement.minimal_cycles(edge_angle)(V, VE)

    return convert(Lar.ChainOp, SparseArrays.transpose(EF))
end

function minimal_3cycles(V::Lar.Points, EV::Lar.ChainOp, FE::Lar.ChainOp)

@show V::Lar.Points;
@show SparseArrays.findnz(EV::Lar.ChainOp);
@show SparseArrays.findnz(FE::Lar.ChainOp);

	triangulated_faces = Array{Any, 1}(undef, FE.m)

    function face_angle(e::Int, f::Int)
        if !isassigned(triangulated_faces, f)
            vs_idxs = Array{Int64, 1}()
            edges_idxs = FE[f, :].nzind
            edge_num = length(edges_idxs)
            edges = zeros(Int64, edge_num, 2)

            for (i, ee) in enumerate(edges_idxs)
                edge = EV[ee, :].nzind
                edges[i, :] = edge
                vs_idxs = union(vs_idxs, edge)
            end

            #vs = V[vs_idxs, :]
			fv,edges = Lar.vcycle(EV, FE, f)

			vs = V[fv, :]

            v1 = LinearAlgebra.normalize(vs[2, :] - vs[1, :])
            v2 = [0 0 0]		# added for debug
            v3 = [0 0 0]
            err = 1e-8
            i = 3
            while -err < norm(v3) < err
                v2 = normalize(vs[i, :] - vs[1, :])
                v3 = cross(v1, v2)
                i = i + 1
            end
            M = reshape([v1; v2; v3], 3, 3)

            #vs = vs*M
			vs = (vs*M)[:, 1:2]

            # triangulated_faces[f] = Triangle.constrained_triangulation(
            #     Array{Float64,2}(vs), vs_idxs, edges, fill(true, edge_num))
			v = convert(Lar.Points, vs'[1:2,:])
			vmap = Dict(zip(fv,1:length(fv))) # vertex map
			mapv = Dict(zip(1:length(fv),fv)) # inverse vertex map
			trias = Lar.triangulate2d(v,edges)
			triangulated_faces[f] = [[mapv[v] for v in tria] for tria in trias]
        end
        edge_vs = EV[e, :].nzind

        t = findfirst(x->edge_vs[1] in x && edge_vs[2] in x, triangulated_faces[f])

        v1 = normalize(V[edge_vs[2], :] - V[edge_vs[1], :])

        if abs(v1[1]) > abs(v1[2])
            invlen = 1. / sqrt(v1[1]*v1[1] + v1[3]*v1[3])
            v2 = [-v1[3]*invlen, 0, v1[1]*invlen]
        else
            invlen = 1. / sqrt(v1[2]*v1[2] + v1[3]*v1[3])
            v2 = [0, -v1[3]*invlen, v1[2]*invlen]
        end

        v3 = cross(v1, v2)

        M = reshape([v1; v2; v3], 3, 3)

@show triangulated_faces;
@show f;
@show t;
        triangle = triangulated_faces[f][t]
        third_v = setdiff(triangle, edge_vs)[1]
        vs = V[[edge_vs..., third_v], :]*M

        v = vs[3, :] - vs[1, :]
        angle = atan(v[2], v[3])
        return angle
    end

    #EF = FE'
    EF = convert(Lar.ChainOp, LinearAlgebra.transpose(FE))
		println(">>>>>>>> going to loop?")
Verbose = true
@show Verbose
		FC = Lar.Arrangement.minimal_cycles(face_angle, true)(V, EF)  # , EV)

		#FC'
    return -convert(Lar.ChainOp, LinearAlgebra.transpose(FC))
end

function minimal_cycles(angles_fn::Function, verbose=true)
println(">>>>>>>>>>>>> eccomi !!")

    function _minimal_cycles(V::Lar.Points,
    		ld_bounds::Lar.ChainOp)  # , EV)

println(">>>>>>>>>>>>> eccomi !!!!")
Verbose = true
@show Verbose

        lld_cellsnum, ld_cellsnum = size(ld_bounds)
        count_marks = zeros(Int64, ld_cellsnum)
        dir_marks = zeros(Int64, ld_cellsnum)
        d_bounds = spzeros(Int64, ld_cellsnum, 0)

        angles = Array{Array{Int64, 1}, 1}(undef, lld_cellsnum)

        function get_seed_cell()
            s = -1
            for i in 1:ld_cellsnum
                if count_marks[i] == 0
                    return i
                elseif count_marks[i] == 1 && s < 0
                    s = i
                end
            end
            return s
        end

        for lld in 1:lld_cellsnum
            as = []
            for ld in ld_bounds[lld, :].nzind
                push!(as, (ld, angles_fn(lld, ld)))
            end
            sort!(as, lt=(a,b)->a[2]<b[2])
            as = map(a->a[1], as)
            angles[lld] = as
        end

if Verbose @show collect(enumerate(angles)) end
if Verbose @show collect(enumerate(Lar.cop2lar(ld_bounds))) end

        function nextprev(lld::Int64, ld::Int64, norp)
            as = angles[lld]
            #ne = findfirst(as, ld)  (findfirst(isequal(v), A), 0)[1]
            ne = (findfirst(isequal(ld), as), 0)[1]
            while true
                ne += norp # next or previous
                if ne > length(as)
                    ne = 1
                elseif ne < 1
                    ne = length(as)
                end

                if count_marks[as[ne]] < 2
                    break
                end
            end
            for k=1:length(count_marks)
            	if count_marks[k]>2  error("TGW is looping") end
            end

            as[ne]
        end

        while (sigma = get_seed_cell()) > 0
						if verbose
                print(Int(floor(50 * sum(count_marks) / ld_cellsnum)), "%\r") # <<<<<<<<<<<<<<<<<<<
            end

            c_ld = spzeros(Int8, ld_cellsnum)
            if count_marks[sigma] == 0
                c_ld[sigma] = 1
            else
                c_ld[sigma] = -dir_marks[sigma]
            end
            c_lld = ld_bounds*c_ld
            while c_lld.nzind != []
								corolla = spzeros(Int64, ld_cellsnum)
								#corolla = zeros(Int64, ld_cellsnum)

                for tau in c_lld.nzind # when looping, loops here !!
                    b_ld = ld_bounds[tau, :]
                    pivot = intersect(c_ld.nzind, b_ld.nzind)[1]
                    adj = nextprev(tau, pivot, sign(-c_lld[tau]))
                    corolla[adj] = c_ld[pivot]
                    if b_ld[adj] == b_ld[pivot]
                        corolla[adj] *= -1
                    end
                end

                c_ld += corolla
                c_lld = ld_bounds*c_ld
            end
            map(s->count_marks[s] += 1, c_ld.nzind)
            for k=1:length(count_marks)
            	if count_marks[k]>2  error("TGW is looping") end
            end
            map(s->dir_marks[s] = c_ld[s], c_ld.nzind)
            d_bounds = [d_bounds c_ld]

        end
        return d_bounds

    end

    return _minimal_cycles
end


text:
using Mimi

@defcomp radiativeforcing begin
    FORC = Variable(index=[time]) # Increase in radiative forcing (watts per m^2 from 1900)

    MAT = Parameter(index=[time]) # Carbon concentration increase in atmosphere (GtC from 1750)
    MATSUM = Parameter(index=[time]) #Sum of MAT[t] and MAT[t+1] to use in FORC[t] for radiativeforcing component
    mat1 = Parameter()
    forcoth = Parameter(index=[time]) # Exogenous forcing for other greenhouse gases
    fco22x = Parameter() # Forcings of equilibrium CO2 doubling (Wm-2)

    function run_timestep(p, v, d, t)
        #Define function for FORC
        if is_first(t)
            v.FORC[t] = p.fco22x * ((log10((((p.MAT[t] + p.mat1)/2)+0.000001)/596.4)/log10(2))) + p.forcoth[t]  #TEMP NOTE: Uses mat1 because it's given as a parameter...not calculated so couldn't use MATSUM
        else
            v.FORC[t] = p.fco22x * ((log10((((p.MATSUM[t])/2)+0.000001)/596.4)/log10(2))) + p.forcoth[t]
        end
    end
end


text:
using Revise
using AdFem
using PyCall
using LinearAlgebra
using PyPlot
using SparseArrays
using MAT
using ADCMEKit
np = pyimport("numpy")
using PyPlot
using SpecialFunctions

# Parameters setup
NT = 100
Δt = ones(NT)*0.0666
Δu = 1.0 # fixed Dirichlet boundary condition on the left hand side
n = 10
m = 5n
h = 0.4
left = bcnode("left", m, n, h)
right = bcnode("right", m, n, h)
dof = (m+1)*(n+1)
Force = zeros(NT, dof)

# Precompute the mass matrices, stiffness matrix, and damping matrix
H = diagm(0=>ones(2))
M = compute_fem_mass_matrix1(m, n, h)
K = compute_fem_stiffness_matrix1(H, m, n, h)
C = 0.0 * M + 0.0 * K

# Compute the force term
αt = αscheme_time(Δt)
ub = zeros(dof, NT)
for i = 1:NT
  ub[left, i] .= Δu
  ub[right, i] .= 0.0
  # Force at αt[i]
end
Force = (-K * ub)'|>Array

# Initial conditions
d0 = zeros(dof)
v0 = zeros(dof)
a0 = zeros(dof)

idof = ones(Bool, dof)
idof[[left;right]] .= false

a0[idof] = M[idof, idof]\(Force[1,idof])

# α scheme time stepping
function solver(A, rhs)
  rhs = scatter_update(rhs, [left;right], zeros(length(left)+length(right)))
  A, _ = fem_impose_Dirichlet_boundary_condition1(A, [left;right], m, n, h)
  A\rhs
end

d, v, a = αscheme(M, C, K, Force, d0, v0, a0, Δt; solve = solver)

# Simulation
sess = Session()
d_, v_, a_ = run(sess, [d, v, a])

d_[:, left] .+= Δu
# visualize_potential(permutedims(reshape(d_, NT, m+1, n+1), [1,3,2])[:,1:n,1:m], m, n, h)

close("all")
for (k,tid) in enumerate(LinRange{Int64}(1, NT+1, 5))
  t = (tid-1)*Δt[1]
  plot((d_[tid, :])[(1:m+1)],"C$k-", label="$t", markersize=3)
end
legend()

max_stars_repo_name: ADCMEMarket/ADCMEImages

max_issues_repo_name: ADCMEMarket/ADCMEImages

max_forks_repo_name: ADCMEMarket/ADCMEImages


text:
function makeTensorT{T}(factorW::Array{T,2}, measureOperator::Array{T,1} = ones(T, size(factorW,1)))
	leng = size(factorW,1)
	tensorT = zeros(T,leng,leng,leng,leng)
	for a = 1:leng, x = 1:leng, xp = 1:leng, y = 1:leng, yp = 1:leng
		tensorT[x,xp,y,yp] += measureOperator[a] * factorW[a,x] * factorW[a,xp] * factorW[a,y] * factorW[a,yp]
	end
	return tensorT
end

function makeTensorP{T}(factorW::Array{T,2})
	leng = size(factorW,1)
	tensorP = zeros(T,leng,leng,leng)
	for a = 1:leng, x = 1:leng, xp = 1:leng, s = 1:leng
		tensorP[x,xp,s] += factorW[a,x] * factorW[a,xp] * factorW[a,s]
	end
	return tensorP
end

function makeTensorQ{T}(factorW::Array{T,2})
	leng = size(factorW,1)
	tensorQ = zeros(T,leng,leng)
	for a = 1:leng, x = 1:leng, y = 1:leng
		tensorQ[x,y] += factorW[a,x] * factorW[a,y]
	end
	return tensorQ
end


text:
using LinearAlgebra
using Distributions
using Plots

include("./gen_data.jl")
include("./kalman_smoothers.jl")

"""Experimental parameters"""

# Time horizon
T = 100

# Experimental parameters
transition_coeffs = 1.0
emission_coeffs = 1.0

# Noises
process_noise = 0.5
measurement_noise = 0.1

# Prior state
state0 = ([0.], reshape([1e-12], 1,1))

"Generate data"

# Generate signal
observations, states = LGDS(transition_coeffs,
                            emission_coeffs,
                            process_noise,
                            measurement_noise,
                            state0;
                            time_horizon=T)

# Check signal visually
plot(1:T, states[2:end], linewidth=3, color="red", label="states")
scatter!(1:T, observations[1:end], color="blue", label="observations")

"""Basic Kalman smoother"""

# Call filter
msk, Psk = rts_smoother(observations,
                        transition_coeffs,
                        emission_coeffs,
                        process_noise,
                        measurement_noise,
                        state0)

# Visualize estimates
scatter(1:T, observations[1,:], color="black", label="observations")
plot!(1:T, states[1,2:end], color="red", label="latent states")
plot!(1:T, msk[:], color="purple", label="inferred")
plot!(1:T, msk[:],
      ribbon=[sqrt.(Psk[1,1,:]), sqrt.(Psk[1,1,:])],
      color="purple", alpha=0.1, label="")
xlabel!("time [k]")
ylabel!("signal")
savefig("./figures/LGDS_rts.png")

max_stars_repo_name: wmkouw/Bayesian-smoothers

max_issues_repo_name: wmkouw/Bayesian-smoothers

max_forks_repo_name: wmkouw/Bayesian-smoothers


text:
function hvdc_comparison(
      parent_folder::String,
      output_folder::String,
      objective::String;
      k_cond=[],
      k_ins=[],
      dc_subnet::Int64=-1,
      gen_areas::Array=[],
      area_transfer::Array=[],
      gen_zones::Array=[],
      zone_transfer::Array=[],
      enum_branches::Bool=false,
      plot_best_x::Int64=-1,
      print_results::Bool=false,
      series_labels::Array=[],
      results_folders::Array=[],
      scopf::Bool=false,
      no_converter_loss::Bool=false,
      output_to_files::Bool=true,
      regularize_f::Float64=0.0
   )
   if (length(k_cond)) > 0 && (length(k_ins) > 0)
      if scopf
         dc_params = ([["sn",contingency,"sn",dc_subnet,"k_ins"],["sn",contingency,"sn",dc_subnet,"k_cond"],["sn",contingency,"sn",dc_subnet,"f_max"]],[k_ins, k_cond, zeros(length(k_ins))])
      else
         dc_params = ([["sn",dc_subnet,"k_ins"],["sn",dc_subnet,"k_cond"],["sn",dc_subnet,"f_max"]],[k_ins, k_cond, zeros(length(k_ins))])
      end
   else
      dc_params = ()
   end
   num_series = 0
   if length(dc_params) > 0
      if length(dc_params[1]) > 0
         num_series = length(dc_params[2][1])
      end
   end
   println("num_series: $num_series")
   println("dc_params: $dc_params")
   # Run the opf for the base case then all dc configurations
   results_dict_allplots = []
   n_subnets = 0
   subnet_array = []
   idx_sorted = []
   if length(results_folders) == 0
      # Base case opf
      (results_dict,n_subnets,subnet_array,idx_sorted,series_output_folder,plot_best_x) = run_series(
         parent_folder,
         objective;
         gen_areas,
         area_transfer,
         gen_zones,
         zone_transfer,
         enum_branches,
         plot_best_x,
         print_results,
         suffix=series_labels[1],
         no_converter_loss=no_converter_loss,
         output_to_files=output_to_files,
         regularize_f=regularize_f
      )
      push!(results_dict_allplots, results_dict)

      results_dict_string = JSON.json(results_dict)
      println("Saving results at $(series_output_folder)/all_results_dict.json")
      open("$series_output_folder/all_results_dict.json", "w") do f
         write(f, results_dict_string)
      end
      push!(results_folders, series_output_folder)

      # Iterate through the parameter collections in order
      for param_i in 1:num_series
         override_param = Dict{Any, Any}()
         for (val_idx,value) in enumerate(dc_params[2])
            set_nested!(override_param, string.(dc_params[1][val_idx]), value[param_i])
         end
         (results_dict,n_subnets_dc,subnet_array_dc,idx_sorted_dc,series_output_folder,plot_best_x) = run_series(
            parent_folder,
            objective;
            gen_areas,
            area_transfer,
            gen_zones,
            zone_transfer,
            enum_branches,
            plot_best_x,
            print_results,
            override_param=override_param,
            suffix=series_labels[1+param_i],
            no_converter_loss=no_converter_loss,
            output_to_files=output_to_files,
            regularize_f=regularize_f
         )
         push!(results_dict_allplots, results_dict)

         results_dict_string = JSON.json(results_dict)
         println("Saving results at $(series_output_folder)/all_results_dict.json")
         open("$series_output_folder/all_results_dict.json", "w") do f
            write(f, results_dict_string)
         end
         push!(results_folders, series_output_folder)
      end
      println("results_folders: $results_folders")
   else
      # Use results data that has already been saved
      results_dict_allplots = []
      for (i,output_folder) in enumerate(results_folders)

         results_dict_parsed = JSON.parsefile("$(output_folder)/all_results_dict.json")
         if i == 1
            (n_subnets, subnet_array, idx_sorted) = get_sort_indices(results_dict_parsed, objective, gen_areas, area_transfer, gen_zones, zone_transfer)
         end

         push!(results_dict_allplots, results_dict_parsed)
      end
   end
   if output_to_files
      println("Saving results at $(output_folder)/results_dict_allplots.json")
      results_dict_string = JSON.json(results_dict_allplots)
      if !isdir(output_folder)
         mkpath(output_folder)
      end
      open("$(output_folder)/results_dict_allplots.json", "w") do f
         write(f, results_dict_string)
      end
   end
   plot_results_dicts_bar(results_dict_allplots,n_subnets,subnet_array,idx_sorted,output_folder,plot_best_x,series_labels)

end


text:
type Generator
  # specified
  t::Vector{Float64}
  n::Int64
  "number of dimensions."
  k::Int64

  # computed
  d::Tuple
  y::Matrix{Float64}
  dt::Vector{Float64}
  sdt::Vector{Float64}
  b::Matrix{Float64}
  db::Matrix{Float64}

  function Generator(t, s, k::Int)
    g      = new(t, length(t), k)
    g.d    = (s..., k)
    g.dt   = diff(t)
    g.sdt  = sqrt(g.dt)
    g.b    = zeros(g.d)
    g.db   = zeros(g.d)
    g.y    = zeros(g.d)
    return g
  end
end

@inline function next{P <: AProcess}(process::P, g::Generator, i::Int64)
  if i == 1
    return init!(process, g)
  end

  randn!(g.db)
  cmul!(g.db, g.sdt[i-1])
  result = next!(process, g, i)
  # Important to use previous b for process computation not to have bias.
  add!(g.b, g.db)
  result
end

function cumsim{P <: AProcess}(process::P, t, k::Int=1)
  s=size(vwrap(initial(process)))
  gen = Generator(vec(t), s, k)
  result = zeros(gen.n, s..., k)
  for i in 1:gen.n
    y = next(process, gen, i)
    result[i,:,:] .= y[:,:]
  end
  CumsimResult(collapse_but_first(result), t)
end

function sim{P <: AProcess}(process::P, t, k::Int=1)
  s=size(vwrap(initial(process)))
  gen = Generator(t, s, k)
  for i in 1:(gen.n-1)
    next(process, gen, i)
  end
  SimResult(collapse_but_first(next(process, gen, gen.n)), t[end], gen.b)
end


text:
export solve,
    Solver,
    solver_startsolutions,
    paths_to_track,
    parameter_homotopy,
    linear_subspace_homotopy

struct SolveStats
    regular::Threads.Atomic{Int}
    regular_real::Threads.Atomic{Int}
    singular::Threads.Atomic{Int}
    singular_real::Threads.Atomic{Int}
end
SolveStats() = SolveStats(Threads.Atomic{Int}.((0, 0, 0, 0))...)

function init!(SS::SolveStats)
    SS.regular[] = SS.regular_real[] = SS.singular[] = SS.singular_real[] = 0
    SS
end

function update!(stats::SolveStats, R::PathResult)
    is_success(R) || return stats

    if is_singular(R)
        Threads.atomic_add!(stats.singular_real, Int(is_real(R)))
        Threads.atomic_add!(stats.singular, 1)
    else
        Threads.atomic_add!(stats.regular_real, Int(is_real(R)))
        Threads.atomic_add!(stats.regular, 1)
    end
    stats
end

"""
    Solver(path_tracker; seed = nothing)

A struct containing multiple copies of `path_tracker`. This contains all pre-allocated
data structures to call [`solve`]. The most convenient way to construct a `Solver` is
via [`solver_startsolutions`](@ref).
"""
struct Solver{T<:AbstractPathTracker}
    trackers::Vector{T}
    seed::Union{Nothing,UInt32}
    stats::SolveStats
    start_system::Union{Nothing,Symbol}
end
Solver(
    tracker::AbstractPathTracker;
    seed::Union{Nothing,UInt32} = nothing,
    start_system = nothing,
) = Solver([tracker], seed, SolveStats(), start_system)

Base.show(io::IO, solver::Solver) = print(io, typeof(solver), "()")

"""
    solver_startsolutions(args...; kwargs...)

Takes the same input as [`solve`](@ref) but instead of directly solving the problem
returns a [`Solver`](@ref) struct and the start solutions.

## Example

Calling `solve(args..; kwargs...)` is equivalent to
```julia
solver, starts = solver_startsolutions(args...; kwargs...)
solve(solver, starts)
```
"""
function solver_startsolutions(
    F::AbstractVector{Expression},
    starts = nothing;
    parameters = Variable[],
    variables = setdiff(variables(F), parameters),
    variable_ordering = variables,
    variable_groups = nothing,
    kwargs...,
)
    sys = System(
        F,
        variables = variable_ordering,
        parameters = parameters,
        variable_groups = variable_groups,
    )
    solver_startsolutions(sys, starts; kwargs...)
end
function solver_startsolutions(
    F::AbstractVector{<:MP.AbstractPolynomial},
    starts = nothing;
    parameters = similar(MP.variables(F), 0),
    variables = setdiff(MP.variables(F), parameters),
    variable_ordering = variables,
    variable_groups = nothing,
    target_parameters = nothing,
    kwargs...,
)
    # handle special case that we have no parameters
    # to shift the coefficients of the polynomials to the parameters
    # this was the behaviour of HC.jl v1
    if isnothing(target_parameters) && isempty(parameters)
        sys, target_parameters = ModelKit.system_with_coefficents_as_params(
            F,
            variables = variable_ordering,
            variable_groups = variable_groups,
        )
    else
        sys = System(
            F,
            variables = variable_ordering,
            parameters = parameters,
            variable_groups = variable_groups,
        )
    end
    solver_startsolutions(sys, starts; target_parameters = target_parameters, kwargs...)
end
function solver_startsolutions(
    F::Union{System,AbstractSystem},
    starts = nothing;
    seed = rand(UInt32),
    start_system = isnothing(variable_groups(F)) ? :polyhedral : :total_degree,
    generic_parameters = nothing,
    p₁ = generic_parameters,
    start_parameters = p₁,
    p₀ = generic_parameters,
    target_parameters = p₀,
    compile::Union{Bool,Symbol} = COMPILE_DEFAULT[],
    start_subspace = nothing,
    target_subspace = nothing,
    intrinsic = nothing,
    kwargs...,
)
    !isnothing(seed) && Random.seed!(seed)

    used_start_system = nothing
    if start_subspace !== nothing
        if target_parameters !== nothing
            F = fix_parameters(F, target_parameters; compile = compile)
        end
        H = linear_subspace_homotopy(
            F,
            start_subspace,
            target_subspace;
            intrinsic = intrinsic,
            compile = compile,
        )
        tracker = EndgameTracker(H; kwargs...)
    else
        if target_subspace !== nothing
            F = slice(F, target_subspace, compile = compile)
        end
        if start_parameters !== nothing
            tracker = parameter_homotopy_tracker(
                F;
                start_parameters = start_parameters,
                target_parameters = target_parameters,
                compile = compile,
                kwargs...,
            )
        elseif start_system == :polyhedral
            used_start_system = :polyhedral
            tracker, starts = polyhedral(
                F;
                compile = compile,
                target_parameters = target_parameters,
                kwargs...,
            )
        elseif start_system == :total_degree
            used_start_system = :total_degree
            tracker, starts = total_degree(
                F;
                compile = compile,
                target_parameters = target_parameters,
                kwargs...,
            )
        else
            throw(
                KeywordArgumentException(
                    :start_system,
                    start_system,
                    "Possible values are: `:polyhedral` and `:total_degree`.",
                ),
            )
        end
    end

    Solver(tracker; seed = seed, start_system = used_start_system), starts
end

function solver_startsolutions(
    G::Union{System,AbstractSystem},
    F::Union{System,AbstractSystem},
    starts = nothing;
    seed = rand(UInt32),
    tracker_options = TrackerOptions(),
    endgame_options = EndgameOptions(),
    kwargs...,
)
    !isnothing(seed) && Random.seed!(seed)
    if F isa SlicedSystem && G isa SlicedSystem && system(F) == system(G)
        H = linear_subspace_homotopy(system(F), linear_subspace(F), linear_subspace(G))
    else
        H = start_target_homotopy(G, F; kwargs...)
    end
    tracker =
        EndgameTracker(H; tracker_options = tracker_options, options = endgame_options)

    Solver(tracker; seed = seed), starts
end

function parameter_homotopy_tracker(
    F::Union{System,AbstractSystem};
    tracker_options = TrackerOptions(),
    endgame_options = EndgameOptions(),
    kwargs...,
)
    H = parameter_homotopy(F; kwargs...)
    EndgameTracker(H; tracker_options = tracker_options, options = endgame_options)
end

"""
    parameter_homotopy(F; start_parameters, target_parameters)

Construct a [`ParameterHomotopy`](@ref). If `F` is homogeneous, then a random affine chart
is chosen (via [`AffineChartHomotopy`](@ref)).
"""
function parameter_homotopy(
    F::Union{System,AbstractSystem};
    generic_parameters = randn(ComplexF64, nparameters(F)),
    p₁ = generic_parameters,
    start_parameters = p₁,
    p₀ = generic_parameters,
    target_parameters = p₀,
    tracker_options = TrackerOptions(),
    endgame_options = EndgameOptions(),
    compile::Union{Bool,Symbol} = COMPILE_DEFAULT[],
    kwargs...,
)
    unsupported_kwargs(kwargs)
    m, n = size(F)
    H = ParameterHomotopy(fixed(F; compile = compile), start_parameters, target_parameters)
    f = System(F)
    if is_homogeneous(f)
        vargroups = variable_groups(f)
        if vargroups === nothing
            m ≥ (n - 1) || throw(FiniteException(n - 1 - m))
            H = on_affine_chart(H)
        else
            m ≥ (n - length(vargroups)) || throw(FiniteException(n - length(vargroups) - m))
            H = on_affine_chart(H, length.(vargroups,) .- 1)
        end
    else
        m ≥ n || throw(FiniteException(n - m))
    end

    H
end

"""
    linear_subspace_homotopy(F, V::LinearSubspace, W::LinearSubspace, intrinsic = nothing)

Constructs an [`IntrinsicSubspaceHomotopy`](@ref) (if `dim(V) < codim(V)` or
`intrinsic = true`) or [`ExtrinsicSubspaceHomotopy`](@ref).
Compared to the direct constructor, this also takes care of homogeneous systems.
"""
function linear_subspace_homotopy(
    F::Union{System,AbstractSystem},
    V::LinearSubspace,
    W::LinearSubspace;
    compile::Union{Bool,Symbol} = COMPILE_DEFAULT[],
    intrinsic = nothing,
)
    # Intrinsic: m+1,dim(V)+1
    # Extrinsic: m+codim(V), dim(V) + codim(V)
    if dim(V) < codim(V) || something(intrinsic, false)
        if is_linear(V) && is_linear(W) && is_homogeneous(System(F))
            IntrinsicSubspaceHomotopy(on_affine_chart(F; compile = compile), V, W)
        else
            IntrinsicSubspaceHomotopy(F, V, W; compile = compile)
        end
    else
        H = ExtrinsicSubspaceHomotopy(F, V, W; compile = compile)
        if is_linear(V) && is_linear(W) && is_homogeneous(System(F))
            on_affine_chart(H)
        else
            H
        end
    end
end

function start_target_homotopy(
    G::Union{System,AbstractSystem},
    F::Union{System,AbstractSystem};
    start_parameters = nothing,
    target_parameters = nothing,
    compile::Union{Bool,Symbol} = COMPILE_DEFAULT[],
    γ = 1.0,
    gamma = γ,
    kwargs...,
)
    unsupported_kwargs(kwargs)
    f, g = System(F), System(G)

    size(F) == size(G) || error("The provided systems don't have the same size.")
    is_homogeneous(f) == is_homogeneous(g) ||
        error("The provided systems are not both homogeneous.")
    variable_groups(f) == variable_groups(g) ||
        error("The provided systems don't decalare the same variable groups.")

    m, n = size(F)

    G = fixed(G; compile = compile)
    if !isnothing(start_parameters)
        G = FixedParameterSystem(G, start_parameters)
    end

    F = fixed(F; compile = compile)
    if !isnothing(target_parameters)
        F = FixedParameterSystem(F, target_parameters)
    end

    H = StraightLineHomotopy(G, F; gamma = gamma)
    if is_homogeneous(f)
        vargroups = variable_groups(f)
        if vargroups === nothing
            m ≥ (n - 1) || throw(FiniteException(n - 1 - m))
            H = on_affine_chart(H)
        else
            m ≥ (n - length(vargroups)) || throw(FiniteException(n - length(vargroups) - m))
            H = on_affine_chart(H, length.(vargroups,) .- 1)
        end
    else
        m ≥ n || throw(FiniteException(n - m))
    end

    H
end

function solver_startsolutions(
    H::Union{Homotopy,AbstractHomotopy},
    starts = nothing;
    compile::Union{Bool,Symbol} = COMPILE_DEFAULT[],
    seed = nothing,
    kwargs...,
)
    !isnothing(seed) && Random.seed!(seed)
    Solver(EndgameTracker(fixed(H; compile = compile)); seed = seed), starts
end

"""
    solve(f; options...)
    solve(f, start_solutions; start_parameters, target_parameters, options...)
    solve(f, start_solutions; start_subspace, target_subspace, options...)
    solve(g, f, start_solutions; options...)
    solve(homotopy, start_solutions; options...)

Solve the given problem. If only a single polynomial system `f` is given, then all
(complex) isolated solutions are computed.
If a system `f` depending on parameters together with start and target parameters is given
then a parameter homotopy is performed.
If two systems `g` and `f` with solutions of `g` are given then the solutions are tracked
during the deformation of `g` to `f`.
Similarly, for a given homotopy `homotopy` ``H(x,t)`` with solutions at ``t=1`` the solutions
at ``t=0`` are computed.
See the documentation for examples.
If the input is a *homogeneous* polynomial system, solutions on a random affine chart of
projective space are computed.

## General Options
The `solve` routines takes the following options:
* `catch_interrupt = true`: If this is `true`, the computation is gracefully stopped and a
  partial result is returned when the computation is interruped.
* `compile = $(COMPILE_DEFAULT[])`: If `true` then a `System` (resp. `Homotopy`) is compiled
  to a straight line program ([`CompiledSystem`](@ref) resp. [`CompiledHomotopy`](@ref))
  for evaluation. This induces a compilation overhead. If `false` then the generated program
  is only interpreted ([`InterpretedSystem`](@ref) resp. [`InterpretedHomotopy`](@ref)).
  This is slower than the compiled version, but does not introduce compilation overhead.
* `endgame_options`: The options and parameters for the endgame.
  See [`EndgameOptions`](@ref).
* `seed`: The random seed used during the computations. The seed is also reported in the
  result. For a given random seed the result is always identical.
* `show_progress= true`: Indicate whether a progress bar should be displayed.
* `stop_early_cb`: Here it is possible to provide a function (or any callable struct) which
  accepts a `PathResult` `r` as input and returns a `Bool`. If `stop_early_cb(r)` is `true`
  then no further paths are tracked and the computation is finished. This is only called
  for successfull paths. This is for example useful if you only want to compute one solution
  of a polynomial system. For this `stop_early_cb = _ -> true` would be sufficient.
* `threading = true`: Enable multi-threading for the computation. The number of
  available threads is controlled by the environment variable `JULIA_NUM_THREADS`.
* `tracker_options`: The options and parameters for the path tracker.
  See [`TrackerOptions`](@ref).

## Options depending on input

If only a polynomial system is given:
* `start_system`: Possible values are `:total_degree` and `:polyhedral`. Depending on the
  choice furhter options are possible. See also [`total_degree`](@ref) and
  [`polyhedral`](@ref).

If a system `f` depending on parameters together with start parameters (or start subspace),
start solutions and *multiple* target parameters (or target subspaces) then the following
options are also available:

* `flatten`: Flatten the output of `transform_result`. This is useful for example if
   `transform_result` returns a vector of solutions, and you only want a single vector of
   solutions as the result (instead of a vector of vector of solutions).
* `transform_parameters = identity`: Transform a parameters values `p` before passing it to
  `target_parameters = ...`.
* `transform_result`: A function taking two arguments, the `result` and the
  parameters `p`. By default this returns the tuple `(result, p)`.

## Basic example

```julia-repl
julia> @var x y;

julia> F = System([x^2+y^2+1, 2x+3y-1])
System of length 2
 2 variables: x, y

 1 + x^2 + y^2
 -1 + 2*x + 3*y

julia> solve(F)
Result with 2 solutions
=======================
• 2 non-singular solutions (0 real)
• 0 singular solutions (0 real)
• 2 paths tracked
• random seed: 0x75a6a462
• start_system: :polyhedral
```
"""
function solve(
    args...;
    show_progress::Bool = true,
    threading::Bool = Threads.nthreads() > 1,
    catch_interrupt::Bool = true,
    target_parameters = nothing,
    stop_early_cb = always_false,
    # many parameter options,
    transform_result = nothing,
    transform_parameters = identity,
    flatten = nothing,
    target_subspaces = nothing,
    kwargs...,
)
    many_parameters = false
    if target_subspaces !== nothing
        many_parameters = true
        solver, starts = solver_startsolutions(
            args...;
            target_subspace = first(target_subspaces),
            kwargs...,
        )
        target_parameters = target_subspaces
    elseif target_parameters !== nothing
        # check if we have many parameters solve
        if !isa(transform_parameters(first(target_parameters)), Number)
            many_parameters = true
            solver, starts = solver_startsolutions(
                args...;
                target_parameters = transform_parameters(first(target_parameters)),
                kwargs...,
            )
        else
            solver, starts = solver_startsolutions(
                args...;
                target_parameters = target_parameters,
                kwargs...,
            )
        end
    else
        solver, starts = solver_startsolutions(args...; kwargs...)
    end
    if many_parameters
        solve(
            solver,
            starts,
            target_parameters;
            show_progress = show_progress,
            threading = threading,
            catch_interrupt = catch_interrupt,
            transform_result = transform_result,
            transform_parameters = transform_parameters,
            flatten = flatten,
        )
    else
        solve(
            solver,
            starts;
            stop_early_cb = stop_early_cb,
            show_progress = show_progress,
            threading = threading,
            catch_interrupt = catch_interrupt,
        )
    end
end

solve(S::Solver, R::Result; kwargs...) =
    solve(S, solutions(R; only_nonsingular = true); kwargs...)
solve(S::Solver, s::AbstractVector{<:Number}; kwargs...) = solve(S, [s]; kwargs...)
solve(S::Solver, starts; kwargs...) = solve(S, collect(starts); kwargs...)
function solve(
    S::Solver,
    starts::AbstractArray;
    stop_early_cb = always_false,
    show_progress::Bool = true,
    threading::Bool = Threads.nthreads() > 1,
    catch_interrupt::Bool = true,
)
    n = length(starts)
    progress = show_progress ? make_progress(n; delay = 0.3) : nothing
    init!(S.stats)
    if threading
        threaded_solve(
            S,
            starts,
            progress,
            stop_early_cb;
            catch_interrupt = catch_interrupt,
        )
    else
        serial_solve(S, starts, progress, stop_early_cb; catch_interrupt = catch_interrupt)
    end
end
(solver::Solver)(starts; kwargs...) = solve(solver, starts; kwargs...)
track(solver::Solver, s; kwargs...) = track(solver.trackers[1], s; kwargs...)

function make_progress(n::Integer; delay::Float64 = 0.0)
    desc = "Tracking $n paths... "
    barlen = min(ProgressMeter.tty_width(desc, stdout, false), 40)
    progress =
        ProgressMeter.Progress(n; dt = 0.2, desc = desc, barlen = barlen, output = stdout)
    progress.tlast += delay
    progress
end
function update_progress!(progress, stats, ntracked)
    t = time()
    if ntracked == progress.n || t > progress.tlast + progress.dt
        showvalues = make_showvalues(stats, ntracked)
        ProgressMeter.update!(progress, ntracked; showvalues = showvalues)
    end
    nothing
end
@noinline function make_showvalues(stats, ntracked)
    showvalues = (("# paths tracked", ntracked),)
    nsols = stats.regular[] + stats.singular[]
    nreal = stats.regular_real[] + stats.singular_real[]
    (
        ("# paths tracked", ntracked),
        ("# non-singular solutions (real)", "$(stats.regular[]) ($(stats.regular_real[]))"),
        ("# singular endpoints (real)", "$(stats.singular[]) ($(stats.singular_real[]))"),
        ("# total solutions (real)", "$(nsols[]) ($(nreal[]))"),
    )
end
update_progress!(::Nothing, stats, ntracked) = nothing

function serial_solve(
    solver::Solver,
    starts,
    progress = nothing,
    stop_early_cb = always_false;
    catch_interrupt::Bool = true,
)
    path_results = Vector{PathResult}()
    tracker = solver.trackers[1]
    try
        for (k, s) in enumerate(starts)
            r = track(tracker, s; path_number = k)
            push!(path_results, r)
            update!(solver.stats, r)
            update_progress!(progress, solver.stats, k)
            if is_success(r) && stop_early_cb(r)
                break
            end
        end
    catch e
        (catch_interrupt && isa(e, InterruptException)) || rethrow(e)
    end

    Result(path_results; seed = solver.seed, start_system = solver.start_system)
end
function threaded_solve(
    solver::Solver,
    S::AbstractArray,
    progress = nothing,
    stop_early_cb = always_false;
    catch_interrupt::Bool = true,
)
    N = length(S)
    path_results = Vector{PathResult}(undef, N)
    interrupted = false
    started = Threads.Atomic{Int}(0)
    finished = Threads.Atomic{Int}(0)
    try
        Threads.resize_nthreads!(solver.trackers)
        tasks = map(enumerate(solver.trackers)) do (i, tracker)
            @tspawnat i begin
                while (k = Threads.atomic_add!(started, 1) + 1) ≤ N && !interrupted
                    r = track(tracker, S[k]; path_number = k)
                    path_results[k] = r
                    nfinished = Threads.atomic_add!(finished, 1) + 1
                    update!(solver.stats, r)
                    update_progress!(progress, solver.stats, nfinished[])
                    if is_success(r) && stop_early_cb(r)
                        interrupted = true
                    end
                end
            end
        end
        for task in tasks
            wait(task)
        end
    catch e
        if (
            isa(e, InterruptException) ||
            (isa(e, TaskFailedException) && isa(e.task.exception, InterruptException))
        )
            interrupted = true
        end
        if !interrupted || !catch_interrupt
            rethrow(e)
        end
    end
    # if we got interrupted we need to remove the unassigned filedds
    if interrupted
        assigned_results = Vector{PathResult}()
        for i = 1:started[]
            if isassigned(path_results, i)
                push!(assigned_results, path_results[i])
            end
        end
        Result(assigned_results; seed = solver.seed, start_system = solver.start_system)
    else
        Result(path_results; seed = solver.seed, start_system = solver.start_system)
    end
end

function start_parameters!(solver::Solver, p)
    for tracker in solver.trackers
        start_parameters!(tracker, p)
    end
    solver
end

function target_parameters!(solver::Solver, p)
    for tracker in solver.trackers
        target_parameters!(tracker, p)
    end
    solver
end
function parameters!(solver::Solver, p, q)
    for tracker in solver.trackers
        parameters!(tracker, p, q)
    end
    solver
end

"""
    paths_to_track(f; optopms..)

Returns the number of paths tracked when calling [`solve`](@ref) with the given arguments.
"""
function paths_to_track(
    f::Union{System,AbstractSystem};
    start_system::Symbol = :polyhedral,
    kwargs...,
)
    paths_to_track(f, Val(start_system); kwargs...)
end

#############################
### Many parameter solver ###
#############################

struct ManySolveStats
    solutions::Threads.Atomic{Int}
end

function solve(
    S::Solver,
    starts,
    target_parameters;
    show_progress::Bool = true,
    threading::Bool = Threads.nthreads() > 1,
    catch_interrupt::Bool = true,
    transform_result = nothing,
    transform_parameters = nothing,
    flatten = nothing,
)
    transform_result = something(transform_result, tuple) # (solutions ∘ first) ∘ tuple
    transform_parameters = something(transform_parameters, identity)
    flatten = something(flatten, false)

    n = length(target_parameters)

    progress = show_progress ? make_many_progress(n; delay = 0.3) : nothing
    many_solve(
        S,
        starts,
        target_parameters,
        progress,
        transform_result,
        transform_parameters,
        Val(flatten);
        catch_interrupt = catch_interrupt,
        threading = threading,
    )
end

function make_many_progress(n::Integer; delay::Float64 = 0.0)
    desc = "Solving for $n parameters... "
    barlen = min(ProgressMeter.tty_width(desc, stdout, false), 40)
    progress =
        ProgressMeter.Progress(n; dt = 0.3, desc = desc, barlen = barlen, output = stdout)
    progress.tlast += delay
    progress
end
function update_many_progress!(progress, results, k, paths_per_param; flatten::Bool)
    t = time()
    if k == progress.n || t > progress.tlast + progress.dt
        showvalues = make_many_showvalues(results, k, paths_per_param; flatten = flatten)
        ProgressMeter.update!(progress, k; showvalues = showvalues)
    end
    nothing
end
@noinline function make_many_showvalues(results, k, paths_per_param; flatten::Bool)
    if flatten
        [
            ("# parameters solved", k),
            ("# paths tracked", paths_per_param * k),
            ("# results", length(results)),
        ]
    else
        [("# parameters solved", k), ("# paths tracked", paths_per_param * k)]
    end
end
update_many_progress!(::Nothing, results, k, paths_per_param; kwargs...) = nothing

function many_solve(
    solver::Solver,
    starts,
    many_target_parameters,
    progress,
    transform_result,
    transform_parameters,
    ::Val{flatten};
    threading::Bool,
    catch_interrupt::Bool,
) where {flatten}
    q = first(many_target_parameters)
    target_parameters!(solver, transform_parameters(q))
    if threading
        res = threaded_solve(solver, starts; catch_interrupt = false)
    else
        res = serial_solve(solver, starts; catch_interrupt = false)
    end
    if flatten
        results = transform_result(res, q)
        if !(results isa AbstractArray)
            throw(ArgumentError("Cannot flatten arguments of type `$(typeof(results))`"))
        end
    else
        results = [transform_result(res, q)]
    end
    k = 1
    m = length(starts)
    update_many_progress!(progress, results, k, m; flatten = flatten)
    try
        for q in Iterators.drop(many_target_parameters, 1)
            target_parameters!(solver, transform_parameters(q))
            if threading
                res = threaded_solve(solver, starts; catch_interrupt = false)
            else
                res = serial_solve(solver, starts; catch_interrupt = false)
            end

            if flatten
                append!(results, transform_result(res, q))
            else
                push!(results, transform_result(res, q))
            end
            k += 1
            update_many_progress!(progress, results, k, m; flatten = flatten)
        end
    catch e
        if !(
            isa(e, InterruptException) ||
            (isa(e, TaskFailedException) && isa(e.task.exception, InterruptException))
        )
            rethrow(e)
        end
    end

    results
end


text: Dict{String, Any}("a" => Dict{String, Any}("b" => Dict{String, Any}("c" => Dict{String, Any}("answer" => Dict{String, Any}("value" => "42", "type" => "integer"))), "better" => Dict{String, Any}("value" => "43", "type" => "integer")))


text:
abstract type AbstractRating{T <: Number}
end

struct PredictRating{T <: Number} <: AbstractRating{T}
    value::Real
    target::T
    PredictRating(x::Real, preference::Preference{T}) where T <: Number = new{T}(correct(x, preference), round(x, preference))
end

struct Rating{T <: Number} <: AbstractRating{T}
    value::T
    Rating(x::T, preference::Preference{T}) where T <: Number = new{T}(correct(x, preference))
end

struct MissingRating{T <: Number} <: AbstractRating{T}
end

MissingRating() = MissingRating{Number}()

value(rating::Rating) = rating.value
value(rating::PredictRating) = rating.value
value(rating::MissingRating) = missing

Base.isnan(rating::AbstractRating{T}) where T <: Number = false
Base.isnan(rating::MissingRating{T}) where T <: Number = true

Base.zero(::Type{AbstractRating{T}}) where T <: Number = MissingRating{T}()

function convert(values::Array{T}, preference::Preference{T})::Array{AbstractRating{T}} where T <: Number
    [Rating(value, preference) for value in values]
end

function convert(values::Array{T})::Array{AbstractRating{T}} where T <: Number
    convert(values, Preference(values))
end

Base.:+(r1::AbstractRating{T}, r2::AbstractRating{T}) where {T <: Number} = value(r1) + value(r2)
Base.:+(r1::AbstractRating{T}, r2::MissingRating{T}) where {T <: Number} = value(r1)
Base.:+(r1::MissingRating{T}, r2::AbstractRating{T}) where {T <: Number} = value(r2)
Base.:+(r1::MissingRating{T}, r2::MissingRating{T}) where {T <: Number} = MissingRating{T}()

Base.:-(r1::AbstractRating{T}, r2::AbstractRating{T}) where {T <: Number} =value(r1) - value(r2)
Base.:-(r1::AbstractRating{T}, r2::MissingRating{T}) where {T <: Number} = value(r1)
Base.:-(r1::MissingRating{T}, r2::AbstractRating{T}) where {T <: Number} = value(r2)
Base.:-(r1::MissingRating{T}, r2::MissingRating{T}) where {T <: Number} = MissingRating{T}()

Base.:+(r1::AbstractRating{T}, x::Number) where {T <: Number} = value(r1) + x
Base.:+(x::Number, r1::AbstractRating{T}) where {T <: Number} = value(r1) + x

Base.:-(r1::AbstractRating{T}, x::Number) where {T <: Number} = value(r1) - x
Base.:-(x::Number, r1::AbstractRating{T}) where {T <: Number} = value(r1) - x

Base.:*(r1::AbstractRating{T}, x::Number) where {T <: Number} = value(r1) * x
Base.:*(x::Number, r1::AbstractRating{T}) where {T <: Number} = value(r1) * x

Base.:/(r1::AbstractRating{T}, x::Number) where {T <: Number} = value(r1) / x
Base.:/(x::Number, r1::AbstractRating{T}) where {T <: Number} = value(r1) / x


text:
module TaylorExponentialMatrix

include("expm2.jl")
include("expm3.jl")

end # module


text:
# Run JITL, SITL, and PITL first to generate their logs. Then run this script to
# compare them.

module Temp

using Plots
using HDF5

include(joinpath(Pkg.dir("Eratosthenes"), "src", "EratosthenesRotations.jl"))
using .EratosthenesRotations

# Nominal
t, q_BI = HDF5.h5open("out/basic-cube.h5", "r") do logs
    (read(logs, "/cube1/truth/time"),
     read(logs, "/cube1/truth/data/q_BI"))
end

# SITL
t_sitl, q_BI_sitl = HDF5.h5open("out/basic-cube-sitl.h5", "r") do logs
    (read(logs, "/cube1/truth/time"),
     read(logs, "/cube1/truth/data/q_BI"))
end

# PITL
t_pitl, q_BI_pitl = HDF5.h5open("out/basic-cube-pitl.h5", "r") do logs
    (read(logs, "/cube1/truth/time"),
     read(logs, "/cube1/truth/data/q_BI"))
end

# Scalars are logged as 1-by-n (don't know why). Convert to just n for
# the sake of plotting.
t      = dropdims(t, dims=1)
t_sitl = dropdims(t_sitl, dims=1)
t_pitl = dropdims(t_pitl, dims=1)

# Compare the literal values in the quaternion (not the underlying rotation);
# they should be identical.
Δq_sitl = q_BI - q_BI_sitl
display(maximum(Δq_sitl, 2))
Δq_pitl = q_BI - q_BI_pitl
display(maximum(Δq_pitl, 2))
Δq_pitl_wrt_sitl = q_BI_sitl - q_BI_pitl
display(maximum(Δq_pitl_wrt_sitl, 2))

# Ok, also compare the angular differences.
θ_sitl = map(k -> q2aa(qpos(qdiff(q_BI[:,k], q_BI_sitl[:,k])))[1], 1:size(q_BI, 2))
display(maximum(θ_sitl) * 180/pi)
θ_pitl = map(k -> q2aa(qpos(qdiff(q_BI[:,k], q_BI_pitl[:,k])))[1], 1:size(q_BI, 2))
display(maximum(θ_pitl) * 180/pi)
θ_pitl_wrt_sitl = map(k -> q2aa(qpos(qdiff(q_BI_sitl[:,k], q_BI_pitl[:,k])))[1], 1:size(q_BI_sitl, 2))
display(maximum(θ_pitl_wrt_sitl) * 180/pi)

# Show the divergence.
# pyplot()
plotlyjs()
# display(plot(t, transpose(Δq),
#              label  = ["Δq1" "Δq2" "Δq3" "Δq4"],
#              ylabel = "Δq"))
display(plot(t, θ_pitl * (180./pi * 3600.),
             xlabel = "Time (s)",
             ylabel = "θ (arcseconds)"))

end