text:
module CanEvolve

using StatsBase
using Distributions
using LinearAlgebra
using DataFrames
using Random
using Plots
using StatsPlots
using PyCall

const np = PyNULL()
function __init__()
    copy!(np, pyimport_conda("numpy", "numpy"))
end

export

    # simulation.jl
    cell,
    mutate,
    primeSim,
    evolveTumour,
    BetaBin,
    sequence,
    detectableDrivers,
    processTumour,
    adjustProportion,
    mutationMetrics,
    simulation,
    full_synthetic_neutral,

    # automation.jl
    sampleParameters,
    autoSimulation,
    autoSimulation_fixed,
    autoSimulation_fixed_depth,
    plotVAF,
    frequencyThresholds,
    mapfrequency,
    uniformdensity,
    rescale,
    engineer,
    engineer_un

include("simulation.jl")
include("automation.jl")

end


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule glslang_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("glslang")
JLLWrappers.@generate_main_file("glslang", UUID("6fd5517d-459c-5c5a-9b1a-c968b4e37a81"))
end  # module glslang_jll

ext: jl


text:
"""
An evaluator specialised to dealing with regression problems.
"""
type RegressionEvaluator <: Evaluator
  stage::AbstractString
  num_samples::Int
  samples::Vector{Tuple{Float, Float}}

  RegressionEvaluator(stage::AbstractString, samples::Vector{Tuple{Float, Float}}) =
    new(stage, length(samples), samples)
end

"""
Provides a definition for a regression evaluator.
"""
type RegressionEvaluatorDefinition <: EvaluatorDefinition
  stage::AbstractString
  samples::Vector{Tuple{Float, Float}}
end

"""
Composes a regression evaluator from a provided definition.
"""
function compose!(def::RegressionEvaluatorDefinition, p::Population)
  if def.stage == ""
    def.stage = p.demes[1].species.genotype.label
  end
  RegressionEvaluator(def.stage, def.samples)
end

"""
The regression evaluator measures the fitness of a candidate solution against
the agreement with a series of provided sample observations from an oracle.

TODO: FOR NOW THESE OBSERVATIONS ARE FIXED.
"""
regression() =
  regression("", [
    (-10.0, 9090.0),
    (-9.0, 5904.0),
    (-8.0, 3640.0),
    (-7.0, 2100.0),
    (-6.0, 1110.0),
    (-5.0, 520.0),
    (-4.0, 204.0),
    (-3.0, 60.0),
    (-2.0, 10.0),
    (-1.0, 0.0),
    (0.0, 0.0),
    (1.0, 4.0),
    (2.0, 30.0),
    (3.0, 120.0),
    (4.0, 340.0),
    (5.0, 780.0),
    (6.0, 1554.0),
    (7.0, 2800.0),
    (8.0, 4680.0),
    (9.0, 7380.0)
  ])
regression(samples::Vector{Tuple{Float, Float}}) = regression("", samples)
regression(stage::AbstractString, samples::Vector{Tuple{Float, Float}}) =
  RegressionEvaluatorDefinition(stage, samples)

"""
Evaluates all unevaluated individuals within a provided deme according to
a given regression evaluator.
"""
function evaluate!(e::RegressionEvaluator, s::State, d::Deme)
  for (id, phenome) in enumerate(d.offspring.stages[e.stage])
    d.offspring.fitnesses[id] = evaluate!(e, s, d.species.fitness, get(phenome))
  end
  s.evaluations += length(d.offspring.stages[e.stage])
end

"""
Evaluates all unevaluated individuals contained within a given state according
to a given regression evaluator.
"""
evaluate!(e::RegressionEvaluator, s::State) =
  for deme in s.population.demes; evaluate!(e, s, deme); end

function evaluate!{T}(e::RegressionEvaluator, s::State, sc::FitnessScheme, candidate::T)
  sse = zero(Float)
  diff = zero(Float)
  for (x, y) in e.samples
    diff = y - execute(candidate, x)
    sse += diff * diff
  end
  assign(sc, sse)
end


text:
using Base.Dates
using Base.Test

using MarketData

using TimeSeries

@testset "modify" begin

@testset "update method works" begin
    new_cls  = update(cl, today(), 111.11)
    new_clv  = update(cl, today(), [111.11])
    new_ohlc = update(ohlc, today(), [111.11 222.22 333.33 444.44])
    empty1   = TimeArray(Vector{Date}(), Array{Int}(0,1))
    empty2   = TimeArray(Vector{Date}(), Array{Int}(0,2))

    @testset "empty time arrays can be constructed" begin
        @test length(empty1) == 0
        @test length(empty2) == 0
    end

    @testset "update an empty time array fails" begin
        @test_throws ArgumentError update(empty1, Date(2000,1,1), 10)
        @test_throws ArgumentError update(empty2, Date(2000,1,1), [10 11])
        # @test length(update1)      == 1
        # @test length(update2)      == 1
        # @test update1.timestamp[1] == Date(2000,1,1)
        # @test update2.timestamp[1] == Date(2000,1,1)
        # @test update1.values       == [12]
        # @test update2.values       == [10 11]
        # @test update1.colnames     == ["foo"]
        # @test update2.colnames     == ["foo", "bar"]
        # @test update1.meta         == "bar"
        # @test update2.meta         == "baz"
    end

    @testset "update a single column time array with single value vector" begin
        @test last(new_clv.values) == 111.11
        @test_throws DimensionMismatch update(cl, today(), [111.11, 222.22])
    end

    @testset "update a multi column time array" begin
        # @test last(new_ohlc).values == [111.11 222.22 333.33 444.44]
        @test tail(new_ohlc).values == [111.11 222.22 333.33 444.44]
        @test_throws MethodError update(ohlc, today(), [111.11, 222.22, 333.33])
    end

    @testset "cannot update more than one observation at a time" begin
        @test_throws(
            MethodError,
            update(cl, [Date(2002,1,1), Date(2002,1,2)], [111.11, 222,22]))
    end

    @testset "cannot update oldest observations" begin
        @test_throws ArgumentError update(cl, Date(1999,1,1), [111.11])
        @test_throws ArgumentError update(cl, Date(1999,1,1), 111.11)
    end

    @testset "cannot update in-between observations" begin
        @test_throws ArgumentError update(cl, Date(2000,1,8), [111.11])
        @test_throws ArgumentError update(cl, Date(2000,1,8), 111.11)
    end
end

@testset "rename method works" begin
    re_ohlc = rename(ohlc, ["a","b","c","d"])
    re_cl   = rename(cl, ["vector"])
    re_cls  = rename(cl, "string")

    @testset "change colnames with multi-member vector" begin
        @test colnames(re_ohlc) == ["a","b","c","d"]
        @test_throws DimensionMismatch rename(ohlc, ["a"])
    end

    @testset "change colnames with single-member vector" begin
        @test colnames(re_cl) == ["vector"]
        @test_throws DimensionMismatch rename(cl, ["a", "b"])
    end

    @testset "change colnames with string" begin
        @test colnames(re_cls) == ["string"]
        @test_throws MethodError rename(cl, "string_a", "string_b")
    end
end

end  # @testset "modify"


text:
using PowerSystems
using NLsolve
const PSY = PowerSystems

raw_file = joinpath(TEST_FILES_DIR, "benchmarks/psse/EXST1/TVC_System_32.raw")
dyr_file = joinpath(TEST_FILES_DIR, "benchmarks/psse/EXST1/TVC_System.dyr")

sys = System(raw_file, dyr_file)
source = first(get_components(Source, sys))
PSY.set_X_th!(source, 0.01)

for l in get_components(PSY.PowerLoad, sys)
    PSY.set_model!(l, PSY.LoadModels.ConstantImpedance)
end

for gen in get_components(ThermalStandard, sys)
    #Find the generator at bus 2
    if get_number(get_bus(gen)) == 2
        dynamic_injector = get_dynamic_injector(gen)
        machine = deepcopy(get_machine(dynamic_injector))
        shaft = deepcopy(get_shaft(dynamic_injector))
        tg = deepcopy(get_prime_mover(dynamic_injector))
        pss = deepcopy(get_pss(dynamic_injector))
        avr = deepcopy(get_avr(dynamic_injector))
        new_dynamic_injector = DynamicGenerator(
            name = get_name(gen),
            machine = machine,
            shaft = shaft,
            avr = avr_exst1(),
            prime_mover = tg,
            ω_ref = 1.0,
            pss = pss,
        )
        set_dynamic_injector!(gen, nothing)
        add_component!(sys, new_dynamic_injector, gen)
    end
end


text:
#=
Starting in the top left corner of a 2×2 grid, and only being able to move to the right and down, there are exactly 6 routes to the bottom right corner.
How many such routes are there through a 20×20 grid?
=#

# 1x1 = 2
# 2x2 = 6

# binomial coefficient is a way to calculate pascals triangle
# binomial(row, column)
# = row! / (column!*(row-column)!)
# reference: http://www.mathwords.com/b/binomial_coefficients_pascal.htm

function number_of_paths(x)
  x *= 2
  binomial(x,div(x,2))
end
@time println(number_of_paths(20))

max_stars_repo_name: daniel-beard/JuliaProjectEuler

max_issues_repo_name: daniel-beard/JuliaProjectEuler

max_forks_repo_name: daniel-beard/JuliaProjectEuler


text:
using Query
using Base.Test

# 3VL

@test Nullable(true) & Nullable(true) == Nullable(true)
@test Nullable(true) & Nullable(false) == Nullable(false)
@test Nullable(true) & Nullable{Bool}() == Nullable{Bool}()
@test Nullable(false) & Nullable(true) == Nullable(false)
@test Nullable(false) & Nullable(false) == Nullable(false)
@test Nullable(false) & Nullable{Bool}() == Nullable(false)
@test Nullable{Bool}() & Nullable(true) == Nullable{Bool}()
@test Nullable{Bool}() & Nullable(false) == Nullable(false)
@test Nullable{Bool}() & Nullable{Bool}() == Nullable{Bool}()

@test true & Nullable(true) == Nullable(true)
@test true & Nullable(false) == Nullable(false)
@test true & Nullable{Bool}() == Nullable{Bool}()
@test false & Nullable(true) == Nullable(false)
@test false & Nullable(false) == Nullable(false)
@test false & Nullable{Bool}() == Nullable(false)

@test Nullable(true) & true == Nullable(true)
@test Nullable(true) & false == Nullable(false)
@test Nullable(false) & true == Nullable(false)
@test Nullable(false) & false == Nullable(false)
@test Nullable{Bool}() & true == Nullable{Bool}()
@test Nullable{Bool}() & false == Nullable(false)

@test Nullable(true) | Nullable(true) == Nullable(true)
@test Nullable(true) | Nullable(false) == Nullable(true)
@test Nullable(true) | Nullable{Bool}() == Nullable(true)
@test Nullable(false) | Nullable(true) == Nullable(true)
@test Nullable(false) | Nullable(false) == Nullable(false)
@test Nullable(false) | Nullable{Bool}() == Nullable{Bool}()
@test Nullable{Bool}() | Nullable(true) == Nullable(true)
@test Nullable{Bool}() | Nullable(false) == Nullable{Bool}()
@test Nullable{Bool}() | Nullable{Bool}() == Nullable{Bool}()

@test true | Nullable(true) == Nullable(true)
@test true | Nullable(false) == Nullable(true)
@test true | Nullable{Bool}() == Nullable(true)
@test false | Nullable(true) == Nullable(true)
@test false | Nullable(false) == Nullable(false)
@test false | Nullable{Bool}() == Nullable{Bool}()

@test Nullable(true) | true == Nullable(true)
@test Nullable(true) | false == Nullable(true)
@test Nullable(false) | true == Nullable(true)
@test Nullable(false) | false == Nullable(false)
@test Nullable{Bool}() | true == Nullable(true)
@test Nullable{Bool}() | false == Nullable{Bool}()

@test !Nullable(true) == Nullable(false)
@test !Nullable(false) == Nullable(true)
@test !Nullable{Bool}() == Nullable{Bool}()

# null comparisons
@test (Nullable(5)==null) == false
@test (Nullable{Int}()==null) == true
@test (null==Nullable(5)) == false
@test (null==Nullable{Int}()) == true

@test (Nullable(5)!=null) == true
@test (Nullable{Int}()!=null) == false
@test (null!=Nullable(5)) == true
@test (null!=Nullable{Int}()) == false

:+, :-, :!, :~
@test +Nullable(1) == Nullable(+1)
@test +Nullable{Int}() == Nullable{Int}()
@test -Nullable(1) == Nullable(-1)
@test -Nullable{Int}() == Nullable{Int}()
@test ~Nullable(1) == Nullable(~1)
@test ~Nullable{Int}() == Nullable{Int}()

# TODO add ^, / back
for op in (:+, :-, :*, :%, :&, :|, :<<, :>>)
    @eval begin
        @test $op(Nullable(3), Nullable(5)) == Nullable($op(3, 5))
        @test $op(Nullable{Int}(), Nullable(5)) == Nullable{Int}()
        @test $op(Nullable(3), Nullable{Int}()) == Nullable{Int}()
        @test $op(Nullable{Int}(), Nullable{Int}()) == Nullable{Int}()

        @test $op(Nullable{Int}(3), 5) == Nullable($op(3, 5))
        @test $op(3, Nullable{Int}(5)) == Nullable($op(3, 5))
        @test $op(Nullable{Int}(), 5) == Nullable{Int}()
        @test $op(3, Nullable{Int}()) == Nullable{Int}()
    end
end

@test Nullable(3) == Nullable(3)
@test !(Nullable(3) == Nullable(4))
@test !(Nullable{Int}() == Nullable(3))
@test !(Nullable{Float64}() == Nullable(3))
@test !(Nullable(3) == Nullable{Int}())
@test !(Nullable(3) == Nullable{Float64}())
@test Nullable{Int}() == Nullable{Int}()
@test Nullable{Int}() == Nullable{Float64}()

@test Nullable(3) == 3
@test 3 == Nullable(3)
@test !(Nullable(3) == 4)
@test !(4 == Nullable(3))
@test !(Nullable{Int}() == 3)
@test !(3 == Nullable{Int}())

@test !(Nullable(3) != Nullable(3))
@test Nullable(3) != Nullable(4)
@test Nullable{Int}() != Nullable(3)
@test Nullable{Float64}() != Nullable(3)
@test Nullable(3) != Nullable{Int}()
@test Nullable(3) != Nullable{Float64}()
@test !(Nullable{Int}() != Nullable{Int}())
@test !(Nullable{Int}() != Nullable{Float64}())

@test !(Nullable(3) != 3)
@test !(3 != Nullable(3))
@test Nullable(3) != 4
@test 4 != Nullable(3)
@test Nullable{Int}() != 3
@test 3 != Nullable{Int}()

@test Nullable(4) > Nullable(3)
@test !(Nullable(3) > Nullable(4))
@test !(Nullable(4) > Nullable{Int}())
@test !(Nullable{Int}() > Nullable(3))
@test !(Nullable{Int}() > Nullable{Int}())

@test Nullable(4) > 3
@test !(Nullable(3) > 4)
@test !(Nullable{Int}() > 3)

@test 4 > Nullable(3)
@test !(3 > Nullable(4))
@test !(4 > Nullable{Int}())

@test lowercase(Nullable("TEST"))==Nullable("test")
@test lowercase(Nullable{String}())==Nullable{String}()

@test Nullable("TEST")[2:end]==Nullable("EST")
@test Nullable{String}()[2:end]==Nullable{String}()

@test length(Nullable("TEST"))==Nullable(4)
@test length(Nullable{String}())==Nullable{Int}()


text:
function _precompile_()
    #xlclip
    src_array = "1\t2\t3\t4\n5\t6\t7\t8\n9\t10\t11\t12"
    d = xlclip_reader(src_array)
    xlclip_writer(d)
end

max_issues_repo_name: alecloudenback/ExcelUtilities

max_forks_repo_name: alecloudenback/ExcelUtilities


text:
# Code for writing GraphViz Dot files from graphs.

md"""
# Generic Dot Code
"""

md"""
The `graph` argumrnt to these functions should support the
following methods:

  - nodes
  - edges
  - dotID
  - dotnode
  - dotedge

See the documentation for those methods.

Default *no-op* methods are provided here where reasonable.
"""

export dotID, dotgraph, dotnode, dotedge
export dotescape, diarc, dot_attributes_string
export rundot, DotError

"""
   nodes(graph)
Return a collection of all of the nodes of `graph`.
"""
function nodes end

"""
    edges(graph)
Return a collection of all of the edges of `graph`.
Each element of the collections is a Pair associating one
nodeof the graph with another.
"""
function edges end

"""
    dotescape(::AbstractString)::AbstractString
Escape an `ID` in the GraphViz Dot language.  'ID' is
the fundamental token in Dot.
"""
function dotescape(s::AbstractString)::AbstractString
    "\"" *
        replace(s, "\"" => "\\\"") *
        "\""
end

dotescape(s::Symbol) = dotescape(string(s))

md"""
Write a GraphViz dot file describing the panel cutting progression
described in a SearchState.
"""
function dotgraph end

struct DotError <: Exception
    command
    error
    stdout
    stderr
end

function Base.showerror(io::IO, e::DotError)
    write(io, "Error running dot:\n($e.command)\n$(e.stderr)\n")
end

md"""
    rundot(path)
Run the GraphViz dot command on the specified dot file
to produce an SVG file with the same basename.
"""
function rundot(path)
    out = IOBuffer()
    err = IOBuffer()
    cmd = `dot -Tsvg -O $path`
    try
        Base.run(pipeline(cmd; stdout=out, stderr=err))
    catch e
        throw(DotError(cmd, e, String(take!(out)), String(take!(err))))
    end
end

md"""
    dotgraph(path::String, graph, dotstyle)
Write the graph to the specified file.
If the file extension is `dot` then a GraphViz dot file is written.
Otherwise a dot description of `graph` is piped through the dot command,
the output of which will be written to `psth`.
"""
function dotgraph(path::String, graph, dotstyle)
    ext = last(splitext(path))[2:end]
    if ext == "dot"
        open(path, "w") do io
            dotgraph(io, graph, dotstyle)
        end
    else
        ### This code needs work.
        # If the file type is anything other than "dot" then run the
        # dot command to generate the file
        dot = IOBuffer()
        err = IOBuffer()
        cmd = `dot -T$ext -o$path`
        try
            proc = run(pipeline(cmd, stdin=dot,
                                stderr=err
                                # stderr=path*".stderr"
                                ))
            @info("dot exit status", status=proc.exitcode, cmd=cmd)
        catch e
            @warn("Error running dot: $e")
        end
        dotgraph(dot, graph)
        err = read(err, String)
        if length(err) > 0
            @warn("Error running dot", err=err)
        end
    end
    return path
end

md"""
    dotgraph(io::IO, graph, dotstyle)
Write a dot description of `graph` to `io`.
"""
function dotgraph(io::IO, graph, dotstyle)
    write(io, "digraph panels {\n")
    for (word, fun) in (("graph", graph_attributes),
                        ("node", node_attributes),
                        ("edge", edge_attributes))
        attrs = fun(dotstyle)
        if !isempty(attrs)
            write(io, "$word [$(dot_attributes_string(attrs))]\n")
        end
    end
    for node in nodes(graph)
        dotnode(io, graph, dotstyle, node)
    end
    for arc in edges(graph)
        dotedge(io, graph, dotstyle, arc.first, arc.second)
    end
    write(io, "}\n")
end

"""
    dotID(node)
Return a string to be used as the id of node in a GraphViz dot file.
"""
function dotID(x)
    string(x)
end

"""
    dotnode(io::IO, graph, dotstyle, node)
Write a Dot node statement to `io` describing `node`.
"""
function dotnode(io::IO, graph, dotstyle, node)
    id = dotescape(dotID(node))
    attrs = dot_attributes_string(node_attributes(dotstyle, node))
    if !isempty(attrs)
        attrs = " [" * attrs * "]"
    end
    write(io, """  $id$attrs\n""")
end

"""
    dotedge(io::IO, graph, dotstyle, from, to)
Write a Dot edge statement to `io` describing an edge from
`from` to `to`.
"""
function dotedge(io::IO, graph, dotstyle, from, to)
    diarc(io, from, to; edge_attributes(dotstyle, from, to)...)
end

function dot_attributes_string(attrs::Dict)::String
    a = []
    for (key, val) in attrs
        push!(a, "$(dotescape(key))=$(dotescape(val))")
    end
    join(a, "; ")
end

function dot_attributes_string(; kwargs...)::String
    attrs = []
    for (key, val) in kwargs
        push!(attrs, "$(dotescape(key))=$(dotescape(val))")
    end
    join(attrs, "; ")
end

function diarc(io::IO, arc::Pair, kwargs...)
    diarc(io, arc.first, arc.second; kwargs...)
end

"""
diark is a convenience function for `dotedge` methods to call.
"""
function diarc(io::IO, from, to; kwargs...)
    attrs = dot_attributes_string(; kwargs...)
    if !isempty(attrs)
        attrs = " [" * attrs * "]"
    end
    write(io, "  $(dotescape(dotID(from))) -> $(dotescape(dotID(to)))$attrs\n")
end


text:
function ok_year(year_str::String, lo::Int, hi::Int)
    if length(year_str) != 4
        return false
    end
    return lo <= parse(Int, year_str) <= hi
end

function ok_pid(pid_str::String)
    return !isnothing(match(r"^[0-9]{9}$", pid_str))
end

function ok_hcl(hcl_str::String)
    return !isnothing(match(r"^#[0-9a-f]{6}$", hcl_str))
end

function ok_hgt(hgt_str::String)
    if endswith(hgt_str, "in")
        hgt = tryparse(Int, hgt_str[1:end-2])
        lo, hi = 59, 76
    elseif endswith(hgt_str, "cm")
        hgt = tryparse(Int, hgt_str[1:end-2])
        lo, hi = 150, 193
    else
        hgt = Nothing()
    end
    if isnothing(hgt)
        return false
    else
        return lo <= hgt <= hi
    end
end

function ok_ecl(ecl_str::String)
    return ecl_str in Set(["amb", "blu", "brn", "gry", "grn", "hzl", "oth"])
end

expected_keys = Dict([
    ("byr", x -> ok_year(x, 1920, 2002)),
    ("iyr", x -> ok_year(x, 2010, 2020)),
    ("eyr", x -> ok_year(x, 2020, 2030)),
    ("hgt", ok_hgt),
    ("hcl", ok_hcl),
    ("ecl", ok_ecl),
    ("pid", ok_pid),
    ("cid", x -> true),
])

function solve_one()
    data = Dict([("cid", "123")])
    tot = 0

    for line in readlines("data/day_four.txt")
        if length(line) == 0
            if all(haskey(data, key) for key in keys(expected_keys))
                tot += 1
            end
            data = Dict([("cid", "123")])
        else
            for key in split(line, ' ')
                k, v = map(String, split(key, ':'))
                data[k] = v
            end
        end
    end
    return tot
end

function solve_two()
    data = Dict([("cid", "123")])
    tot = 0

    for line in readlines("data/day_four.txt")
        if length(line) == 0
            if all(haskey(data, key) for key in keys(expected_keys))
                tot += 1
            end
            data = Dict([("cid", "123")])
        else
            for key in split(line, ' ')
                k, v = map(String, split(key, ':'))
                if !haskey(expected_keys, k) | expected_keys[k](v)
                    data[k] = v
                end
            end
        end
    end
    return tot
end

println(solve_one())
println(solve_two())

max_stars_repo_name: ColCarroll/adventofcode

max_issues_repo_name: ColCarroll/adventofcode

max_forks_repo_name: ColCarroll/adventofcode


text:
function hans_solver(mdl::Union{ModelInit,Model}, lower::Vector{Float64}=[.975,.975], upper::Vector{Float64}=[1.075,1.075]; method_cube::String="CPU") # , lower = zeros(length(mdl.Prices)), upper = ones(length(mdl.Prices)))

    # initialise model

    mdl = init_model(mdl)
    results = Optim.optimize(prices -> objective_function(prices, mdl, method_cube), mdl.Prices, Optim.NelderMead())
    # results = optimize( prices -> objective_function(prices, mdl, method_cube), lower, upper, mdl.Prices , Fminbox( NelderMead() ) )
    # results = optimize( prices -> objective_function(prices, mdl, method_cube), mdl.Prices , ParticleSwarm(lower, upper, 12 ) )

    if (results.minimum < 10^(-8))

        printstyled("\n\n\nModel Solution found with precision 10^-8 with respect to value function.\nRunning last iteration to generate contents of Model\n\n)", color=28)
        mdl, new_prices = iterate_mdl(results.minimizer, mdl, method_cube)

    end

    if (results.minimum >= 10^(-8))

        printstyled("\n\n\nModel Solution NOT found with precision 10^-8 with respect to value function.\nRunning last iteration to generate contents of Model\n\n)", color=28)
        mdl, new_prices = iterate_mdl(results.minimizer, mdl, method_cube)

    end
    # mdl.Prices = new_prices
    # HA iteration loop
    # make use of dagger + arrayfire
    # also, make it as function.
    return(results, mdl)

end

function hans_solver_dynamic_griding(mdl::Union{ModelInit,Model}, grid_power::Int=2048, lower::Vector{Float64}=[.775,.775], upper::Vector{Float64}=[1.075,1.075]; method_cube::String="CPU") # , lower = zeros(length(mdl.Prices)), upper = ones(length(mdl.Prices)))

    # initialise model

    mdl = init_model(mdl)

    results = Optim.optimize(prices -> objective_function(prices, mdl, method_cube), mdl.Prices, Optim.NelderMead(), Optim.Options(g_tol=10^-4))
    # results = optimize( prices -> objective_function(prices, mdl, method_cube), lower, upper, mdl.Prices , Fminbox( NelderMead() ) )
    # results = optimize( prices -> objective_function(prices, mdl, method_cube), mdl.Prices , ParticleSwarm(lower, upper, 12 ) )

    printstyled("\n\n\nModel Solution found for initial Grid.\nRunning last iteration to generate contents of Model, than updating Grid\n\n", color=28)
    mdl, new_prices = iterate_mdl(results.minimizer, mdl, method_cube)
    # @load "mdl_tmp.jlo" results mdl
    print("\n\n")

    mdl.GridDef = update_grid(mdl, grid_power)
    # mdl = init_model(mdl)

    i = Int(0)

    while i < 12

        mdl = init_model(mdl)
        results = Optim.optimize(prices -> objective_function(prices, mdl, method_cube), mdl.Prices , Optim.NelderMead(),
         Optim.Options(g_tol=10.0^(-4 - i)) )
        mdl, new_prices = iterate_mdl(results.minimizer, mdl, method_cube)

        if (results.minimum < 10^(-12))

            printstyled("\n\n\nModel Solution found with precision 10^-3 with respect to value function.\nRunning last iteration to generate contents of Model\n\n\n", color=28)
            # mdl, new_prices = iterate_mdl( results.minimizer, mdl, method_cube )
            break

        else

            printstyled("\n\n\nModel Solution NOT found with precision 10^-3 with respect to value function.\nUpdating Grid and running next iteration\n\n\n\n\n\n", color=28)
            # @save "mdl_tmp.jlo" results mdl
            mdl.GridDef = update_grid(mdl, grid_power)
            i += 1

        end

    end
    # mdl.Prices = new_prices
    # HA iteration loop
    # make use of dagger + arrayfire
    # also, make it as function.
    return(results, mdl)

end


text:
### These benchmarks can be run from the top-level directory of this repo with
###
###     julia -e 'using PkgBenchmark; results=benchmarkpkg("Quaternionic"); export_markdown("benchmark/results.md", results)'
###
### This runs the benchmarks (possibly tuning them automatically first), and writes the results to a nice markdown file.

using BenchmarkTools
using Random
using Quaternionic

Random.seed!(1234)

const SUITE = BenchmarkGroup()

const N = 1_000_000

SUITE["algebra"] = BenchmarkGroup()
for T in [Float64, Float32, Float16]
    a = collect(randn(Quaternion{T}, N))
    b = collect(randn(Quaternion{T}, N))
    for op in [:+, :-, :*, :/]
        SUITE["algebra"][op, T] = @benchmarkable $op.($a, $b)
    end
end

SUITE["math"] = BenchmarkGroup()
for T in [Float64, Float32, Float16]
    a = collect(randn(Quaternion{T}, N))
    for f in [abs, abs2, absvec, abs2vec, inv, log, exp, sqrt, angle]
        SUITE["algebra"][f, T] = @benchmarkable $f.($a)
    end
end


text:
if get(ENV, "CMBLENSING_TEST_CUDA", false) == "1"
    using Adapt, CUDA
    CUDA.allowscalar(false)
    maybegpu(x) = adapt(CuArray,x)
    storage = CuArray
    @info "Running tests on $(repr("text/plain", CUDA.device()))"
else
    maybegpu(x) = x
    storage = Array
    @info "Running tests on CPU"
end

##

using CMBLensing
using CMBLensing: @SMatrix, @SVector, AbstractCℓs, basis, Basis,
    LinearInterpolation, Measurement, RK4Solver, ±, typealias, BatchedReal

##

using FileIO
using FFTW
using FiniteDifferences
using LinearAlgebra
using Random
using Random: default_rng
using Requires
using Serialization
using SparseArrays
using Test
using Zygote

##

try
    push!(LOAD_PATH, "@v#.#") # assumes you have CirculantCov in your global environment
    using CirculantCov
catch
    @warn """CirculantCov.jl not found, not testing EquiRect fields.
    Run `pkg> add https://github.com/EthanAnderes/CirculantCov.jl` to add this package.
    """
end

##

macro test_real_gradient(f, x, tol=:(rtol=1e-3))
    esc(:(@test real(gradient($f,$x)[1]) ≈ central_fdm(5,1)($f,$x) $tol))
end

Nsides     = [(8,8), (4,8), (8,4)]
Nsides_big = [(128,128), (64,128), (128,64)]

Random.seed!(1)

has_batched_fft = (FFTW.fftw_provider != "mkl") || (storage != Array)

##

@testset "CMBLensing" begin

##

# basic printing sanity checks, which were super annoying to get right...
# see also: https://discourse.julialang.org/t/dispatching-on-the-result-of-unwrap-unionall-seems-weird/25677

@testset "Printing" begin

    # concrete types:
    for f in [maybegpu(FlatMap(rand(4,4))), maybegpu(FlatQUMap(rand(4,4),rand(4,4)))]
        @test occursin("pixel",sprint(show, MIME("text/plain"), f))
    end

end

##

@testset "Flat" begin

    Ny,Nx = Nside = first(Nsides)

    @testset "Constructors" begin

        @testset "batch=$D" for D in [(), (3,)]

            Ix = maybegpu(rand(Ny,Nx,D...))
            Il = maybegpu(rand(Ny÷2+1,Nx,D...))

            @testset "$(basis(F))" for (F,ks,args,kwargs) in [
                (FlatMap,        (:Ix,),        (Ix,),      ()),
                (FlatFourier,    (:Il,),        (Il,),      (;Ny,)),
                (FlatQUMap,      (:Qx,:Qx),     (Ix,Ix),    ()),
                (FlatQUFourier,  (:Ql,:Ql),     (Il,Il),    (;Ny,)),
                (FlatEBMap,      (:Ex,:Bx),     (Ix,Ix),    ()),
                (FlatEBFourier,  (:El,:Bl),     (Il,Il),    (;Ny,)),
                (FlatIQUMap,     (:Ix,:Qx,:Qx), (Ix,Ix,Ix), ()),
                (FlatIQUFourier, (:Il,:Ql,:Ql), (Il,Il,Il), (;Ny,)),
                (FlatIEBMap,     (:Ix,:Ex,:Bx), (Ix,Ix,Ix), ()),
                (FlatIEBFourier, (:Il,:El,:Bl), (Il,Il,Il), (;Ny,)),
            ]
                local f
                @test (f = F(args...; kwargs...)) isa F
                @test @inferred(F(getproperty.(Ref(f),ks)..., f.metadata)) == f
                @test (io=IOBuffer(); serialize(io,f); seekstart(io); deserialize(io) == f)
                @test (save(".test_field.jld2", "f", f); load(".test_field.jld2", "f") == f)
                @test_throws ErrorException F(args..., ProjLambert(Nx=Nx+1, Ny=Ny+1))

            end

        end

        rm(".test_field.jld2", force=true)

    end

    @testset "Basis conversions" begin
        @testset "$(typealias(Bin)) → $(typealias(Bout))" for (f,Bin,Bout) in [
            (f,Bin,Bout)
            for (f,Bs) in [
                (FlatMap(rand(Nside...)), (Map,Fourier)),
                (FlatQUMap(rand(Nside...),rand(Nside...)), (QUMap,QUFourier,EBMap,EBFourier))
            ]
            for Bin in Bs
            for Bout in Bs
        ]

            @test basis(@inferred(Bout(Bin(f)))) == Bout
            @test Bin(Bout(Bin(f))) ≈ f

        end
    end

end

##

@testset "Algebra" begin

    @testset "Nside = $Nside" for Nside in Nsides

        fs = ((B0,f0),(B2,f2),(Bt,ft)) = [
            (Fourier,    maybegpu(FlatMap(rand(Nside...)))),
            (EBFourier,  maybegpu(FlatQUMap(rand(Nside...),rand(Nside...)))),
            (Fourier,    maybegpu(FieldTuple(FlatMap(rand(Nside...)),FlatMap(rand(Nside...)))))
            # inference currently broken for this case:
            # (IEBFourier, maybegpu(FlatIQUMap(rand(Nside...),rand(Nside...),rand(Nside...)))),
        ]
        # MKL doesnt seem to support batched FFTs, not that theyre really useful on CPU
        has_batched_fft && append!(fs, [
            (Fourier,    maybegpu(FlatMap(rand(Nside...,2)))),
            (EBFourier,  maybegpu(FlatQUMap(rand(Nside...,2),rand(Nside...,2)))),
        ])

        @testset "f :: $(typeof(f)) " for (B,f) in fs

            local Ðf, Ðv, g, H

            @test similar(f) isa typeof(f)
            @test zero(f) isa typeof(f)
            @test similar(f,Float32) isa Field
            @test eltype(similar(f,Float32)) == Float32

            # used in lots of tests
            @test f ≈ f
            @test !(f ≈ 2f)

            # broadcasting
            @test (@inferred f .+ f) isa typeof(f)
            @test (@inferred f .+ Float32.(f)) isa typeof(f)

            # promotion
            @test (@inferred f + B(f)) isa typeof(f)
            @test (@inferred f + B(Float32.(f))) isa typeof(f)

            # Diagonal broadcasting
            @test (@inferred Diagonal(f) .* Diagonal(f) .* Diagonal(f)) isa typeof(Diagonal(f))

            # inverses
            @test (@inferred pinv(Diagonal(f))) isa Diagonal{<:Any,<:typeof(f)}
            @test_throws Exception inv(Diagonal(0*f))

            # trace
            @test all(tr(Diagonal(f)' * Diagonal(f)) ≈ f'f)
            @test all(tr(Diagonal(f) * Diagonal(f)') ≈ f'f)
            @test_broken all(tr(f*f') ≈ f'f) # broken by (intentionally) removing OuterProdOp

            # Field dot products
            D = Diagonal(f)
            @test (@inferred f' * f) isa Real
            @test (@inferred f' * B(f)) isa Real
            @test (@inferred f' * D * f) isa Real
            @test sum(f, dims=:) ≈ sum(f[:])

            # Explicit vs. lazy DiagOp algebra
            @test (Diagonal(Ð(f)) + Diagonal(Ð(f))) isa DiagOp{<:Field{basis(Ð(f))}}
            @test (Diagonal(Ł(f)) + Diagonal(Ð(f))) isa LazyBinaryOp
            @test (Diagonal(Ł(f)) + Diagonal(Ł(f))) isa DiagOp{<:Field{basis(Ł(f))}}

            if !(f isa FieldTuple)

                # gradients
                @test (Ðf = @inferred ∇[1]*f) isa Field
                @test all(∇[1]'*f ≈ -∇[1]*f)
                @test all(-∇[1]'*f ≈ ∇[1]*f)
                @test (@inferred mul!(Ðf,∇[1],Ð(f))) isa Field
                @test (Ðv = @inferred ∇*f) isa FieldVector
                @test (@inferred mul!(Ðv,∇,Ð(f))) isa FieldVector
                @test ((g,H) = map(Ł, (@inferred gradhess(f)))) isa NamedTuple{<:Any, <:Tuple{FieldVector, FieldMatrix}}

                # FieldVector dot product
                @test (@inferred Diagonal.(g)' * g) isa typeof(g[1])
                @test (@inferred mul!(similar(g[1]), Diagonal.(g)', g)) isa typeof(g[1])

                # FieldMatrix-FieldVector product
                @test (@inferred Diagonal.(H) * g) isa FieldVector
                @test (@inferred Diagonal.(H) * Diagonal.(g)) isa FieldOrOpVector
                @test (@inferred mul!(Diagonal.(similar.(g)), Diagonal.(H), Diagonal.(g))) isa FieldOrOpVector

            end

        end

        # # tuple adjoints
        # f0b = identity.(batch(f0, batch_length(f)))
        # v = similar.(@SVector[f0b, f0b])
        # @test (@inferred mul!(f0b, spin_adjoint(f), f)) isa Field{<:Any,S0}
        # @test (@inferred mul!(v, spin_adjoint(f), @SVector[f,f])) isa FieldVector{<:Field{<:Any,S0}}

        # mixed-spin
        @test (@inferred f0 .* f2) isa typeof(f2)

        # matrix type promotion
        @test (@inferred FlatMap(rand(Float64,2,2)) .+ FlatMap(view(rand(Float32,2,2),:,:))) isa FlatMap{<:Any,Float64,Matrix{Float64}}

    end

end

##

@testset "Log/Trace" begin

    @testset "logdet(Diagonal(::Map))" begin
        @test logdet(Diagonal(FlatMap([1 -2; 3 -4])))                                           ≈  log(24)
        @test logdet(Diagonal(FlatQUMap([1 -2; 3 -4], [1 -2; 3 -4])))                           ≈ 2log(24)
        @test logdet(Diagonal(FlatIQUMap([1 -2; 3 -4], [1 -2; 3 -4], [1 -2; 3 -4])))            ≈ 3log(24)
        @test logdet(Diagonal(FieldTuple(FlatMap([1 -2; 3 -4]), FlatMap([1 -2; 3 -4]))))        ≈ 2log(24)
        @test all(logdet(Diagonal(FlatMap(cat([1 -2; 3 -4],[1 -2; 3 -4],dims=3))))::BatchedReal ≈  log(24))
    end

    @testset "logdet(Diagonal(::Fourier)) Nside=$Nside" for Nside in Nsides_big
        x = rand(Nside...)
        @test logdet(Diagonal(Fourier(FlatMap(x))))                                                   ≈ real( logdet(Diagonal(fft(x)[:])))
        @test logdet(Diagonal(QUFourier(FlatQUMap(x,x))))                                             ≈ real(2logdet(Diagonal(fft(x)[:])))
        @test logdet(Diagonal(IQUFourier(FlatIQUMap(x,x,x))))                                         ≈ real(3logdet(Diagonal(fft(x)[:])))
        @test logdet(Diagonal(FieldTuple(Fourier(FlatMap(x)), Fourier(FlatMap(x)))))                  ≈ real(2logdet(Diagonal(fft(x)[:])))
        has_batched_fft && @test all(logdet(Diagonal(Fourier(FlatMap(cat(x,x,dims=3)))))::BatchedReal ≈ real( logdet(Diagonal(fft(x)[:]))))
    end

    @testset "tr(Diagonal(::Map))" begin
        @test tr(Diagonal(FlatMap([1 -2; 3 -4])))                                           ≈  -2
        @test tr(Diagonal(FlatQUMap([1 -2; 3 -4], [1 -2; 3 -4])))                           ≈  -4
        @test tr(Diagonal(FlatIQUMap([1 -2; 3 -4], [1 -2; 3 -4], [1 -2; 3 -4])))            ≈  -6
        @test tr(Diagonal(FieldTuple(FlatMap([1 -2; 3 -4]), FlatMap([1 -2; 3 -4]))))        ≈  -4
        @test all(tr(Diagonal(FlatMap(cat([1 -2; 3 -4],[1 -2; 3 -4],dims=3))))::BatchedReal ≈  -2)
    end

    @testset "tr(Diagonal(::Fourier)) Nside=$Nside" for Nside in Nsides_big
        x = rand(Nside...)
        @test tr(Diagonal(Fourier(FlatMap(x))))                                                   ≈  tr(Diagonal(fft(x)[:]))
        @test tr(Diagonal(QUFourier(FlatQUMap(x,x))))                                             ≈ 2tr(Diagonal(fft(x)[:]))
        @test tr(Diagonal(IQUFourier(FlatIQUMap(x,x,x))))                                         ≈ 3tr(Diagonal(fft(x)[:]))
        @test tr(Diagonal(FieldTuple(Fourier(FlatMap(x)), Fourier(FlatMap(x)))))                  ≈ 2tr(Diagonal(fft(x)[:]))
        has_batched_fft && @test all(tr(Diagonal(Fourier(FlatMap(cat(x,x,dims=3)))))::BatchedReal ≈ real(tr(Diagonal(fft(x)[:]))))
    end

end

##

@testset "FlatS2" begin
    @testset "Nside = $Nside" for Nside in Nsides
        C = maybegpu(Diagonal(EBFourier(FlatEBMap(rand(Nside...), rand(Nside...)))))
        f = maybegpu(FlatQUMap(rand(Nside...), rand(Nside...)))
        @test C*f ≈ FlatQUFourier(C[:QQ]*f[:Q]+C[:QU]*f[:U], C[:UU]*f[:U]+C[:UQ]*f[:Q])
    end
end

##

@testset "BatchedReal" begin

    @testset "Nside = $Nside" for Nside in Nsides

        r  = 1.
        rb = batch([1.,2])

        @testset "f :: $(typeof(f))" for (f,fb) in [
            (maybegpu(FlatMap(rand(Nside...))),                  maybegpu(FlatMap(rand(Nside...,2)))),
            (maybegpu(FlatQUMap(rand(Nside...),rand(Nside...))), maybegpu(FlatQUMap(rand(Nside...,2),rand(Nside...,2))))
        ]

            @test @inferred(r * f)  == f
            @test @inferred(r * fb) == fb
            @test unbatch(@inferred(rb * f)) == [f, 2f]
            @test unbatch(@inferred(rb * fb)) == [batch_index(fb,1), 2batch_index(fb,2)]

        end

    end

end

##

@testset "Misc" begin

    @testset "Nside = $Nside" for Nside in Nsides

        f = maybegpu(FlatMap(rand(Nside...)))

        @test                  (MidPass(100,200) .* Diagonal(Fourier(f))) isa Diagonal
        @test_throws Exception  MidPass(100,200) .* Diagonal(        f)

    end

end

##

@testset "Cℓs" begin

    @test InterpolatedCℓs(1:100, rand(100))      isa AbstractCℓs{Float64}
    @test InterpolatedCℓs(1:100, rand(100) .± 1) isa AbstractCℓs{Measurement{Float64}}
    @test (InterpolatedCℓs(1:100, 1:100) * ℓ²)[50] == 50^3

end

##

@testset "ParamDependentOp" begin

    D = Diagonal(maybegpu(FlatMap(rand(4,4))))

    @test ParamDependentOp((;x=1, y=1)->x*y*D)() ≈ D
    @test ParamDependentOp((;x=1, y=1)->x*y*D)(z=2) ≈ D
    @test ParamDependentOp((;x=1, y=1)->x*y*D)(x=2) ≈ 2D
    @test ParamDependentOp((;x=1, y=1)->x*y*D)((x=2,y=2)) ≈ 4D # tuple calling form
    @test_throws MethodError ParamDependentOp((;x=1, y=1)->x*y*D)(2) # only Tuple unnamed arg is OK

end

##

@testset "Chains" begin

    chains = CMBLensing.wrap_chains([
        [Dict(:i=>1, :b=>2), Dict(:i=>2       ), Dict(:i=>3, :b=>2)],
        [Dict(:i=>1, :b=>3), Dict(:i=>2, :b=>3), Dict(:i=>3, :b=>3)],
    ])

    # basic
    @test chains[1, 1, :i] == 1
    @test chains[:, 1, :i] == [1,1]
    @test chains[:, :, :i] == [[1, 2, 3], [1, 2, 3]]

    # slices
    @test chains[1, 1:2, :i] == [1, 2]
    @test chains[:, 1:2, :i] == [[1,2], [1,2]]

    # implied : in first dims
    @test chains[:i] == [[1,2,3],[1,2,3]]
    @test chains[1,:i] == [1,2,3]

    # missing
    @test all(chains[1,:b] .=== [2, missing, 2])

end;

##

@testset "Zygote" begin

    @testset "Nside = $Nside" for Nside in Nsides

        @testset "$FMap" for (FMap, FFourier, Npol) in [
            (FlatMap,   FlatFourier,   1),
            (FlatQUMap, FlatQUFourier, 2)
        ]

            Ny,Nx = Nside
            Ixs = collect(maybegpu(rand(Nside...)) for i=1:Npol)
            Ils = rfft.(Ixs)
            f,g,h = @repeated(maybegpu(FMap(Ixs...)),3)
            v = @SVector[f,f]
            D = Diagonal(f)
            A = 2

            @testset "Fields" begin

                @testset "sum" begin
                    @test ((δ = gradient(f -> sum(Map(f)),     Map(f))[1]); basis(δ)==basis(Map(f))     && δ ≈ one(Map(f)))
                    @test ((δ = gradient(f -> sum(Map(f)), Fourier(f))[1]); basis(δ)==basis(Fourier(f)) && δ ≈ Fourier(one(Map(f))))
                end

                @testset "B1=$B1, B2=$B2, B3=$B3" for B1=[Map,Fourier], B2=[Map,Fourier], B3=[Map,Fourier]

                    @test gradient(f -> dot(B1(f),B2(f)), B3(f))[1] ≈ 2f
                    @test gradient(f -> norm(B1(f)), B3(f))[1] ≈ f/norm(f)
                    @test gradient(f -> B1(f)' * B2(g), B3(f))[1] ≈ g
                    @test gradient(f -> sum(Diagonal(Map(f)) * B2(g)), B3(f))[1] ≈ g
                    @test gradient(f -> sum(Diagonal(Map(∇[1]*f)) * B2(g)), B3(f))[1] ≈ ∇[1]'*g
                    @test gradient(f -> B1(f)'*(D\B2(f)), B3(f))[1] ≈ D\f + D'\f
                    @test gradient(f -> (B1(f)'/D)*B2(f), B3(f))[1] ≈ D\f + D'\f
                    @test gradient(f -> B1(f)'*(D*B2(f)), B3(f))[1] ≈ D*f + D'*f
                    @test gradient(f -> (B1(f)'*D)*B2(f), B3(f))[1] ≈ D*f + D'*f
                    if B2==Map
                        @test gradient(f -> B1(f)'*Diagonal(B2(f))*f, B3(f))[1] ≈ @. 3*$B2(f)^2
                    else
                        @test_broken gradient(f -> B1(f)'*Diagonal(B2(f))*f, B3(f))[1] ≈ @. 3*$B2(f)^2
                    end

                end

                @testset "Broadcasting" begin
                    @test gradient(f -> sum(@. f*f + 2*f + 1), f)[1] ≈ 2*f+2
                    @test gradient(f -> sum(@. f^2 + 2*f + 1), f)[1] ≈ 2*f+2
                end

            end

            @testset "FieldVectors" begin

                @test gradient(f -> Map(∇[1]*f)' *     Map(v[1]) + Map(∇[2]*f)' *     Map(v[2]), f)[1] ≈ ∇' * v
                @test gradient(f -> Map(∇[1]*f)' * Fourier(v[1]) + Map(∇[2]*f)' * Fourier(v[2]), f)[1] ≈ ∇' * v
                @test gradient(f -> sum(Diagonal(Map(∇[1]*f)) * v[1] + Diagonal(Map(∇[2]*f)) * v[2]), f)[1] ≈ ∇' * v

                @test gradient(f -> @SVector[f,f]' * Map.(@SVector[g,g]), f)[1] ≈ 2g
                @test gradient(f -> @SVector[f,f]' * Fourier.(@SVector[g,g]), f)[1] ≈ 2g

                @test gradient(f -> sum(sum(@SVector[f,f])),                            f)[1] ≈ 2*one(f)
                @test gradient(f -> sum(sum(@SVector[f,f]      .+ @SVector[f,f])),      f)[1] ≈ 4*one(f)
                @test gradient(f -> sum(sum(@SMatrix[f f; f f] .+ @SMatrix[f f; f f])), f)[1] ≈ 8*one(f)

                # these were broken by (intentionally) removing OuterProdOp .they
                # seem like a fairly unusual, but keeping them here as broken for now...
                @test_broken gradient(f -> sum(Diagonal.(Map.(∇*f))' * Fourier.(v)), f)[1] ≈ ∇' * v
                @test_broken gradient(f -> sum(Diagonal.(Map.(∇*f))' * Map.(v)), f)[1] ≈ ∇' * v
                @test_broken gradient(f -> sum(sum(Diagonal.(@SMatrix[f f; f f]) * @SVector[f,f])), f)[1] ≈ 8*f

            end

            if f isa FlatS0

                @testset "logdet" begin
                    @test gradient(x->logdet(x*Diagonal(Map(f))),     1)[1] ≈ size(Map(f))[1]
                    @test gradient(x->logdet(x*Diagonal(Fourier(f))), 1)[1] ≈ size(Map(f))[1]
                    L = ParamDependentOp((;x=1)->x*Diagonal(Fourier(f)))
                    @test gradient(x->logdet(L(x=x)), 1)[1] ≈ size(Map(f))[1]
                end

                @test gradient(x -> norm(x*Fourier(f)), 1)[1] ≈ norm(f)
                @test gradient(x -> norm(x*Map(f)), 1)[1] ≈ norm(f)

                L₀ = Diagonal(Map(f))
                @test gradient(x -> norm((x*L₀)*f), 1)[1] ≈ norm(L₀*f)
                L₀ = Diagonal(Fourier(f))
                @test gradient(x -> norm((x*L₀)*f), 1)[1] ≈ norm(L₀*f)

                @test gradient(x -> norm((x*Diagonal(Map(f)))*f), 1)[1] ≈ norm(Diagonal(Map(f))*f)
                @test gradient(x -> norm((x*Diagonal(Fourier(f)))*f), 1)[1] ≈ norm(Diagonal(Fourier(f))*f)

            end

            @testset "Array Bailout" begin

                @testset "Fourier" begin

                    @test_real_gradient(A -> logdet(Diagonal(FFourier((A.*Ils)...; Ny))), A)
                    @test_real_gradient(A -> logdet(Diagonal(A*FFourier(Ils...; Ny))), A)
                    @test_real_gradient(A -> logdet(A*Diagonal(FFourier(Ils...; Ny))), A)

                    @test_real_gradient(A -> f' * (Diagonal(FFourier((A.*Ils)...; Ny))) * f, A)
                    @test_real_gradient(A -> f' * (Diagonal(A*FFourier(Ils...; Ny))) * f, A)
                    @test_real_gradient(A -> f' * (A*Diagonal(FFourier(Ils...; Ny))) * f, A)

                    @test_real_gradient(A -> norm(FFourier((A.*Ils)...; Ny)), A)
                    @test_real_gradient(A -> norm(A*FFourier(Ils...; Ny)), A)

                end

                @testset "Map" begin

                    @test_real_gradient(A -> logdet(Diagonal(FMap((A.*Ixs)...))), A)
                    @test_real_gradient(A -> logdet(Diagonal(A*FMap(Ixs...))), A)
                    @test_real_gradient(A -> logdet(A*Diagonal(FMap(Ixs...))), A)

                    @test_real_gradient(A -> f' * (Diagonal(FMap((A.*Ixs)...))) * f, A)
                    @test_real_gradient(A -> f' * (Diagonal(A*FMap(Ixs...))) * f, A)
                    @test_real_gradient(A -> f' * (A*Diagonal(FMap(Ixs...))) * f, A)

                    @test_real_gradient(A -> norm(FMap((A.*Ixs)...)), A)
                    @test_real_gradient(A -> norm(A*FMap(Ixs...)), A)

                end

            end

        end

    end

    @testset "LinearInterpolation" begin
        @test gradient(x->LinearInterpolation([1,2,3],[1,2,3])(x), 2)[1] == 1
        @test gradient(x->LinearInterpolation([1,2,3],[1,x,3])(2), 2)[1] == 1
        @test gradient(x->LinearInterpolation([1,x,3],[1,2,3])(2), 2)[1] == -1
    end

end

##

@testset "Lensing" begin

    local f, ϕ, Lϕ
    Cℓ = camb().unlensed_total

    @testset "$L" for (L,atol) in [(BilinearLens,300), (LenseFlow,0.1)]

        @testset "Nside = ($Ny,$Nx)" for (Ny,Nx) in Nsides_big

            @testset "T :: $T" for T in (Float32, Float64)

                proj = ProjLambert(;Ny,Nx,T,storage)

                Cϕ = maybegpu(Cℓ_to_Cov(:I, proj, Cℓ.ϕϕ))
                @test (ϕ = @inferred simulate(Cϕ)) isa FlatS0

                ## S0
                Cf = maybegpu(Cℓ_to_Cov(:I, proj, Cℓ.TT))
                @test (f = @inferred simulate(Cf)) isa FlatS0
                @test (Lϕ = cache(LenseFlow(ϕ),f)) isa CachedLenseFlow
                @test (@inferred Lϕ*f) isa FlatS0
                # adjoints
                f,g = simulate(Cf),simulate(Cf)
                @test f' * (Lϕ * g) ≈ (f' * Lϕ) * g
                # gradients
                δf, δϕ = simulate(Cf), simulate(Cϕ)
                @test_real_gradient(α -> norm(L(ϕ+α*δϕ)*(f+α*δf)), 0, atol=atol)

                # S2 lensing
                Cf = maybegpu(Cℓ_to_Cov(:P, proj, Cℓ.EE, Cℓ.BB))
                @test (f = @inferred simulate(Cf)) isa FlatS2
                @test (Lϕ = cache(LenseFlow(ϕ),f)) isa CachedLenseFlow
                @test (@inferred Lϕ*f) isa FlatS2
                # adjoints
                f,g = simulate(Cf),simulate(Cf)
                @test f' * (Lϕ * g) ≈ (f' * Lϕ) * g
                # gradients
                δf, δϕ = simulate(Cf), simulate(Cϕ)
                @test_real_gradient(α -> norm(L(ϕ+α*δϕ)*(f+α*δf)), 0, atol=atol)

            end

        end

    end

end

##

@testset "Posterior" begin

    Cℓ = camb()
    L = LenseFlow{RK4Solver{7}}
    T = Float64

    @testset "Nside = $Nside" for Nside in Nsides_big

        @testset "pol = $pol" for pol in (:I,:P)

            @unpack f,f̃,ϕ,ds,ds₀ = load_sim(
                Cℓ       = Cℓ,
                θpix     = 3,
                Nside    = Nside,
                T        = T,
                beamFWHM = 3,
                pol      = pol,
                storage  = storage,
                rng      = default_rng(),
                pixel_mask_kwargs = (edge_padding_deg=1,)
            )
            @unpack Cf,Cϕ = ds₀
            f°,ϕ° = mix(f,ϕ,ds)

            @test lnP(0,f,ϕ,ds) ≈ lnP(1,    f̃,  ϕ , ds) rtol=1e-4
            @test lnP(0,f,ϕ,ds) ≈ lnP(:mix, f°, ϕ°, ds) rtol=1e-4

            δf,δϕ = simulate(Cf, rng=default_rng()), simulate(Cϕ, rng=default_rng())

            @test_real_gradient(α->lnP(0,    f +α*δf, ϕ +α*δϕ, ds), 0, atol=0.5)
            @test_real_gradient(α->lnP(1,    f̃ +α*δf, ϕ +α*δϕ, ds), 0, atol=130)
            @test_real_gradient(α->lnP(:mix, f°+α*δf, ϕ°+α*δϕ, ds), 0, atol=0.5)

        end

    end

end

##

@require CirculantCov="edf8e0bb-e88b-4581-a03e-dda99a63c493" begin

@testset "EquiRect" begin

    θspan  = (π/2 .- deg2rad.((-40,-70)))
    φspan  = deg2rad.((-60, 60))
    φspan′ = deg2rad.((-50, 50))
    Cℓ     = camb()
    rtol   = 1e-5

    @testset "T = $T" for T in (Float32, Float64)

        @testset "Nside = $Nside" for Nside in [(32,64)]

            Ny, Nx = Nside

            # non-periodic
            proj′ = ProjEquiRect(;Ny, Nx, T, θspan, φspan=φspan′)
            @test (proj′.Ny == length(proj′.θ) == Ny) && (proj′.Nx == length(proj′.φ) == Nx)

            # Make a linear list `θedges[1], θ[1], θedges[2], θ[2], ..., θedges[n], θ[n], θedges[n+1]`
            # and test that it is strictly increasing.
            ∂θedges = vcat(vcat(proj′.θedges[1:end-1]', proj′.θ')[:],  proj′.θedges[end])
            @test all(diff(∂θedges) .> 0)

            proj′Δφpix   = rem2pi(proj′.φ[2]-proj′.φ[1],   RoundDown)
            proj′Δφspan  = rem2pi(proj′.φ[end]-proj′.φ[1], RoundDown) + proj′Δφpix
            inputΔφspan = φspan′ |> x->rem2pi(x[end]-x[1], RoundDown)
            @test proj′Δφspan ≈ inputΔφspan

            # with integer fraction of 2π (and Float64)
            proj = ProjEquiRect(;Ny, Nx, T, θspan, φspan)
            P = typeof(proj)
            f0   = EquiRectMap(randn(T, Nside...), proj)
            f2   = EquiRectQUMap(randn(T, Nside...), randn(T, Nside...), proj)

            @test f0 isa EquiRectMap
            @test f2 isa EquiRectQUMap

            # transform
            @test Map(AzFourier(f0)) ≈ f0
            @test QUMap(QUAzFourier(f2)) ≈ f2
            @test real(eltype(Map(AzFourier(f0)))) == T
            @test real(eltype(QUMap(QUAzFourier(f2)))) == T

            # transform (testing equality independent of dot)
            @test AzFourier(f0)[:Ix]   ≈ f0[:Ix]
            @test QUAzFourier(f2)[:Px] ≈ f2[:Px]
            @test Map(f0)[:Il]   ≈ f0[:Il]
            @test QUMap(f2)[:Pl] ≈ f2[:Pl]

            # dot product independent of basis
            @test dot(f0,f0) ≈ dot(AzFourier(f0), AzFourier(f0))
            @test dot(f2,f2) ≈ dot(QUAzFourier(f2), QUAzFourier(f2))

            # creating block-diagonal covariance operators
            Cf0 = Cℓ_to_Cov(:I, f0.proj, Cℓ.total.TT)
            Cf2 = Cℓ_to_Cov(:P, f2.proj, Cℓ.total.EE, Cℓ.total.BB)
            Bf0 = CMBLensing.Cℓ_to_Beam(:I, f0.proj, Cℓ.total.TT)
            Bf2 = CMBLensing.Cℓ_to_Beam(:P, f2.proj, Cℓ.total.TT)
            @test Cf0 isa BlockDiagEquiRect{AzFourier,T}
            @test Cf2 isa BlockDiagEquiRect{QUAzFourier,T}
            @test Bf0 isa BlockDiagEquiRect{AzFourier,T}
            @test Bf2 isa BlockDiagEquiRect{QUAzFourier,T}

            # sqrt
            @test (sqrt(Cf0) * sqrt(Cf0) * f0) ≈ (Cf0 * f0) rtol=rtol
            @test (sqrt(Cf2) * sqrt(Cf2) * f2) ≈ (Cf2 * f2) rtol=rtol

            # simulation
            @test real(eltype(simulate(Cf0) :: EquiRectS0)) == T
            @test real(eltype(simulate(Cf2) :: EquiRectS2)) == T

            # pinv
            @test (pinv(Cf0) * Cf0 * f0) ≈ f0 rtol=rtol
            @test (pinv(Cf2) * Cf2 * f2) ≈ f2 rtol=rtol

            @test (Cf0 \ Cf0 * f0) ≈ f0 rtol=rtol
            @test (Cf2 \ Cf2 * f2) ≈ f2 rtol=rtol

            @test (Cf0 / Cf0 * f0) ≈ f0 rtol=rtol
            @test (Cf2 / Cf2 * f2) ≈ f2 rtol=rtol

            # some operator algebra on ops
            @test (Cf0 + Cf0) * f0 ≈ Cf0 * (2 * f0) ≈ (2 * Cf0) * f0
            @test (Cf2 + Cf2) * f2 ≈ Cf2 * (2 * f2) ≈ (2 * Cf2) * f2

            # logdet
            @test logdet(Cf0) ≈ logabsdet(Cf0)[1]
            @test logdet(Cf2) ≈ logabsdet(Cf2)[1]

            # adjoint
            g0 = simulate(Cf0)
            g2 = simulate(Cf2)
            @test f0' * (Cf0 * g0) ≈ (f0' * Cf0) * g0 rtol=rtol
            @test f2' * (Cf2 * g2) ≈ (f2' * Cf2) * g2 rtol=rtol

            # Test for correct Fourier symmetry in monopole and nyquist f2
            let f2kk = f2[:Pl], f2xx = f2[:Px]
                v = f2kk[1:end÷2,1]
                w = f2kk[end÷2+1:end,1]
                @test v ≈ conj.(w)
                if iseven(size(f2xx,2))
                    v′ = f2kk[1:end÷2,end]
                    w′ = f2kk[end÷2+1:end,end]
                    @test v′ ≈ conj.(w′)
                end
            end

            # gradients w.r.t. fields
            @test gradient(f -> dot(f,f), f0)[1] ≈ 2*f0
            @test gradient(f -> dot(f,f), f2)[1] ≈ 2*f2
            @test gradient(f0 -> f0' * (pinv(Cf0) * f0), f0)[1] ≈ (2 * pinv(Cf0) * f0)
            @test gradient(f2 -> f2' * (pinv(Cf2) * f2), f2)[1] ≈ (2 * pinv(Cf2) * f2)

            # gradients through entries of BlockDiagEquiRect
            @test gradient(α -> logabsdet(α * Cf0)[1], 1)[1] ≈ size(Cf0,1)
            @test gradient(α -> logabsdet(α * Cf2)[1], 1)[1] ≈ size(Cf2,1)
            @test_broken gradient(α -> f0' * (pinv(α * Cf0) * f0), 1)[1] isa Real
            @test_broken gradient(α -> f2' * (pinv(α * Cf2) * f2), 1)[1] isa Real

        end

    end

end

end

##

end


text:
title = "Models"
weight = 10


text:
import Base: +, -, *, /, one, zero

"""
    DBit

A deterministic bit is a single bit representing ±1.

Fields:
- `bit::Bool`: a sample of a bitstream
"""
struct DBit <: AbstractBit
    bit::Bool
end
function DBit(bit::Real)
    if bit != 1 && bit != -1
        error("Cannot create DBit from value $bit (must be ±1).")
    end

    DBit((bit == 1) ? true : false)
end
DBit(bits::VecOrMat{Bool}) = DBit.(bits)

zero(::Type{DBit}) = DBit(false)
zero(::DBit) = zero(DBit)
one(::Type{DBit}) = DBit(true)
one(::DBit) = one(DBit)

"""
    float(b::DBit)

Map `b` to the underlying floating-point value using
``\\{0, 1\\} \\to \\{-1, 1\\}``
"""
float(b::DBit) = b.bit ? 1 : -1

"""
    DBitstream

A deterministic bitstream that looks like a PDM-encoded audio format.

Fields:
- `bits::Queue{DBit}`: the underlying bitstream
"""
struct DBitstream <: AbstractBitstream
    bits::Queue{DBit}
end
DBitstream() = DBitstream(Queue{DBit}())

for op in (:+, :-, :*, :/)
    @eval $(op)(x::DBit, y::DBit) = $(op)(float(x), float(y))
    @eval $(op)(x::DBit, y::Real) = $(op)(float(x), y)
    @eval $(op)(x::Real, y::DBit) = $(op)(x, float(y))
end


text:
"""
variant: test bloomfilter first
"""
function bloom_v4 end

function bloom_v4(t::SearchSequence)
    n = sizeof(t)
    skip = n
    tlast = _nthbyte(t,n)
    bloom_mask = UInt64(_search_bloom_mask(tlast))
    for j in 1:n-1
        bloom_mask |= _search_bloom_mask(_nthbyte(t,j))
        if _nthbyte(t,j) == tlast
            skip = n - j
        end
    end
    return t,bloom_mask,skip,tlast
end

function bloom_v4(s::SearchSequence, t::SearchSequence, i::Integer)
    bloom_v4(s,bloom_v4(t),i)
end

function bloom_v4(s::SearchSequence,p::Tuple,i::Integer,sv::MaybeVector=nothing)
    (t,bloom_mask,skip,tlast) = p # replaces original code, commented out
    DOSTATS = !(sv isa Nothing)
    n = sizeof(t)
    m = sizeof(s)
    if n == 0
        return 1 <= i <= m+1 ? max(1, i) : 0
    elseif m == 0
        return 0
    elseif n == 1
        return something(findnext(isequal(_nthbyte(t,1)), s, i), 0)
    end
    w = m - n
    if w < 0 || i - 1 > w
        return 0
    end

    if DOSTATS loops = 0 end
    if DOSTATS bloomtests = 0 end
    if DOSTATS bloomskips = 0 end
    i +=n-1
    while i <= m
        if DOSTATS loops += 1 end
        if DOSTATS bloomtests += 1 end
        silast = _nthbyte(s,i)
        if bloom_mask & _search_bloom_mask(silast) == 0
            if DOSTATS bloomskips += 1 end
            i += n
        else
            if silast == tlast
                # check candidate
                j = 1
                while _nthbyte(s,i-n+j) == _nthbyte(t,j)
                    j += 1
                    # match found?
                    if j == n
                        if DOSTATS recordcase(sv, loops, bloomtests, bloomskips, bitcount(bloom_mask), skip) end
                        return i-n+1
                    end
                end
                # no match: skip and test bloom
                i += skip
            else
                i +=1
            end
        end
    end
    if DOSTATS recordcase(sv, loops, bloomtests, bloomskips, bitcount(bloom_mask), skip) end
    0
end


text:
using POMDPModels
using POMDPToolbox
using GenerativeModels
using MCTS
using RobustMCTS

mdp = GridWorld()
rmdp = RobustAdapter(mdp)
solver = RobustMCTSSolver(rollout_nature=mdp, c=10.0)
policy = solve(solver, rmdp)
state = initial_state(mdp, MersenneTwister())
a = action(policy, state)

v = MCTS.TreeVisualizer(policy, state)
json = MCTS.create_json(v)


text:
__precompile__(true)

module Elly

using Compat
using ProtoBuf
using URIParser
using MbedTLS

# using from stdlib
using CRC32c
using Random
using UUIDs
using Base64
using Dates
using Sockets
using Serialization
import Sockets: connect
using Distributed
import Distributed: launch, manage

import Base: readdir, show, isfile, isdir, islink, stat, filesize, filemode, mtime, mkdir, mkpath,
        mv, rm, abspath, cd, pwd, touch, open, bytesavailable, cp, joinpath, dirname,
        eof, position, seek, seekend, seekstart, skip, read, write, read!, close, convert
import ProtoBuf: write_bytes, read_bytes, call_method
import URIParser: URI

export show, convert, URI

export UserGroupInformation, add_token, find_tokens, username

export HDFSClient, HDFSFile, HDFSFileInfo,
        hdfs_server_defaults, hdfs_default_block_size, hdfs_default_replication, hdfs_blocks, hdfs_set_replication,
        hdfs_status, hdfs_capacity, hdfs_capacity_used, hdfs_capacity_remaining, hdfs_renewlease,
        isfile, isdir, islink, stat, filesize, filemode, mtime, atime, du, exists, readdir,
        mkdir, mkpath, touch, mv, rm, abspath, cd, pwd, joinpath, dirname,
        eof, position, seek, seekend, seekstart, skip, bytesavailable,
        read!, read, write, readbytes, readall, open, close, cp

export YarnClient, YarnNode, YarnApp, YarnAppStatus, YarnAppAttemptStatus, nodecount, nodes, launchcontext, submit, kill, status, attempts, am_rm_token

export YarnAppMaster, register, unregister, kill, can_schedule_mem, can_schedule_cores, container_allocate, container_release, container_start, container_stop, callback

export YarnManager, launch, manage

const Lock = Channel
makelock() = Channel{Int}(1)

include("hadoop/hadoop.jl")
using Elly.hadoop
using Elly.hadoop.common
using Elly.hadoop.hdfs
using Elly.hadoop.yarn

const ELLY_CLIENTNAME = "elly"

include("ugi.jl")
include("rpc.jl")
include("sasl.jl")
include("envconsts.jl")

include("api_hdfs_base.jl")
include("api_hdfs_io.jl")

include("api_yarn_base.jl")
include("api_yarn_client.jl")
include("api_yarn_appmaster.jl")

include("cluster_manager.jl")

end # module


text:
export
ruleSPGaussianMeanVarianceOutNPP,
ruleSPGaussianMeanVarianceMPNP,
ruleSPGaussianMeanVarianceOutNGP,
ruleSPGaussianMeanVarianceMGNP,
ruleSPGaussianMeanVarianceVGGN,
ruleSPGaussianMeanVarianceVPGN,
ruleSPGaussianMeanVarianceOutNSP,
ruleSPGaussianMeanVarianceMSNP,
ruleSPGaussianMeanVarianceOutNGS,
ruleSPGaussianMeanVarianceMGNS,
ruleVBGaussianMeanVarianceM,
ruleVBGaussianMeanVarianceOut

ruleSPGaussianMeanVarianceOutNPP(   msg_out::Nothing,
                                    msg_mean::Message{PointMass, V},
                                    msg_var::Message{PointMass}) where V<:VariateType =
    Message(V, GaussianMeanVariance, m=deepcopy(msg_mean.dist.params[:m]), v=deepcopy(msg_var.dist.params[:m]))

ruleSPGaussianMeanVarianceMPNP(msg_out::Message{PointMass}, msg_mean::Nothing, msg_var::Message{PointMass}) =
    ruleSPGaussianMeanVarianceOutNPP(msg_mean, msg_out, msg_var)

function ruleSPGaussianMeanVarianceOutNGP(  msg_out::Nothing,
                                            msg_mean::Message{F, V},
                                            msg_var::Message{PointMass}) where {F<:Gaussian, V<:VariateType}

    d_mean = convert(ProbabilityDistribution{V, GaussianMeanVariance}, msg_mean.dist)

    Message(V, GaussianMeanVariance, m=d_mean.params[:m], v=d_mean.params[:v] + msg_var.dist.params[:m])
end

ruleSPGaussianMeanVarianceMGNP(msg_out::Message{F}, msg_mean::Nothing, msg_var::Message{PointMass}) where F<:Gaussian =
    ruleSPGaussianMeanVarianceOutNGP(msg_mean, msg_out, msg_var)

function ruleSPGaussianMeanVarianceVGGN(msg_out::Message{F1, Univariate},
                                        msg_mean::Message{F2, Univariate},
                                        msg_var::Nothing) where {F1<:Gaussian, F2<:Gaussian}

    d_out  = convert(ProbabilityDistribution{Univariate, GaussianMeanVariance}, msg_out.dist)
    d_mean = convert(ProbabilityDistribution{Univariate, GaussianMeanVariance}, msg_mean.dist)

    Message(Univariate, Function, log_pdf=(x)-> -0.5*log(d_out.params[:v] + d_mean.params[:v] + x) - 1/(2*x)*(d_out.params[:m] - d_mean.params[:m])^2)
end

function ruleSPGaussianMeanVarianceVPGN(msg_out::Message{PointMass, Univariate},
                                        msg_mean::Message{F, Univariate},
                                        msg_var::Nothing) where F<:Gaussian

    d_mean = convert(ProbabilityDistribution{Univariate, GaussianMeanVariance}, msg_mean.dist)

    Message(Univariate, Function, log_pdf=(x)-> -0.5*log(d_mean.params[:v] + x) - 1/(2*x)*(msg_out.dist.params[:m] - d_mean.params[:m])^2)
end

# Particle update
function ruleSPGaussianMeanVarianceOutNSP(msg_out::Nothing,
                                          msg_mean::Message{SampleList, V},
                                          msg_var::Message{PointMass}) where {V<:VariateType}

    samples = bootstrap(msg_mean.dist, msg_var.dist)
    weights = msg_mean.dist.params[:w]

    return Message(V, SampleList, s=samples, w=weights)
end

# Particle update
ruleSPGaussianMeanVarianceMSNP(msg_out::Message{SampleList},
                               msg_mean::Nothing,
                               msg_var::Message{PointMass}) = ruleSPGaussianMeanVarianceOutNSP(msg_mean, msg_out, msg_var)

# Particle update
function ruleSPGaussianMeanVarianceOutNGS(  msg_out::Nothing,
                                            msg_mean::Message{F, V},
                                            msg_var::Message{SampleList}) where {F<:Gaussian, V<:VariateType}

    samples = bootstrap(msg_mean.dist, msg_var.dist)
    weights = msg_var.dist.params[:w]

    Message(V, SampleList, s=samples, w=weights)
end

# Particle update
ruleSPGaussianMeanVarianceMGNS(msg_out::Message{F},
                               msg_mean::Nothing,
                               msg_var::Message{SampleList}) where F<:Gaussian = ruleSPGaussianMeanVarianceOutNGS(msg_mean, msg_out, msg_var)

ruleVBGaussianMeanVarianceM(dist_out::ProbabilityDistribution{V},
                            dist_mean::Any,
                            dist_var::ProbabilityDistribution) where V<:VariateType =
    Message(V, GaussianMeanVariance, m=unsafeMean(dist_out), v=unsafeMean(dist_var))

ruleVBGaussianMeanVarianceOut(  dist_out::Any,
                                dist_mean::ProbabilityDistribution{V},
                                dist_var::ProbabilityDistribution) where V<:VariateType =
    Message(V, GaussianMeanVariance, m=unsafeMean(dist_mean), v=unsafeMean(dist_var))


text:
@testset "default" begin
    p = @inferred boxplot([1,2,3,4,5])
    test_ref("references/boxplot/default.txt", @show_col(p))
    p = @inferred boxplot("series1", [1,2,3,4,5])
    test_ref("references/boxplot/default_name.txt", @show_col(p))
    p = @inferred boxplot("series1", [1,2,3,4,5], title = "Test", xlim = (-1,8), color = :blue, width = 50, border = :solid, xlabel="foo")
    test_ref("references/boxplot/default_parameters.txt", @show_col(p))
    p = @inferred boxplot("series1", [1,2,3,4,5], title = "Test", xlim = [-1,8], color = :blue, width = 50, border = :solid, xlabel="foo")
    test_ref("references/boxplot/default_parameters.txt", @show_col(p))
    test_ref("references/boxplot/default_parameters_nocolor.txt", @show_nocol(p))
end

@testset "scaling" begin
    for (i, max_x) in enumerate([5,6,10,20,40])
        p = @inferred boxplot([1,2,3,4,5], xlim=(0,max_x))
        test_ref("references/boxplot/scale$i.txt", @show_col(p))
        p = @inferred boxplot([1,2,3,4,5], xlim=[0,max_x])
        test_ref("references/boxplot/scale$i.txt", @show_col(p))
    end
end

@testset "multi-series" begin
    p = @inferred boxplot(["one", "two"], [[1,2,3,4,5], [2,3,4,5,6,7,8,9]], title="Multi-series", xlabel="foo", color=:yellow)
    test_ref("references/boxplot/multi1.txt", @show_col(p))
    @test @inferred(boxplot!(p, "one more", [-1,2,3,4,11])) === p
    test_ref("references/boxplot/multi2.txt", @show_col(p))
    @test @inferred(boxplot!(p, [4,2,2.5,4,14], name = "last one")) === p
    test_ref("references/boxplot/multi3.txt", @show_col(p))
end

@testset "colors" begin
    p = @inferred boxplot(["one", "two"], [[1,2,3], [4,5,6]], color=214)
    test_ref("references/boxplot/col1.txt", @show_col(p))

    p = @inferred boxplot(["one", "two"], [[1,2,3], [4,5,6]], color=(187,0,187))
    test_ref("references/boxplot/col2.txt", @show_col(p))
end


text:
@testset "ReadWrite" begin

    testdir = joinpath(@__DIR__, "data")
    files = map(x -> joinpath(testdir, x), readdir(testdir))

    for f in files
         res = featherread(f)
         columns, headers = res.columns, res.names

        ncols = length(columns)
        nrows = length(columns[1])

        temp = tempname()
        push!(temps, temp)

        featherwrite(temp, columns, headers, description=res.description, metadata=res.metadata)

         res2 = featherread(temp)
         columns2, headers2 = res2.columns, res2.names

        @test length(columns2) == ncols

        @test headers==headers2

        for (c1,c2) in zip(columns, columns2)
            @test length(c1)==nrows
            @test length(c2)==nrows
            for i = 1:nrows
                @test isequal(c1[i], c2[i])
            end
        end

    @test res.description == res2.description
    @test res.metadata == res2.metadata
    # for (col1,col2) in zip(source.ctable.columns,sink.ctable.columns)
    #     @test col1.name == col2.name
    #     @test col1.metadata_type == col2.metadata_type
    #     @test typeof(col1.metadata) == typeof(col2.metadata)
    #     @test col1.user_metadata == col2.user_metadata

    #     v1 = col1.values; v2 = col2.values
    #     @test v1.dtype == v2.dtype
    #     @test v1.encoding == v2.encoding
    #     # @test v1.offset == v2.offset # currently not python/R compatible due to wesm/feather#182
    #     @test v1.length == v2.length
    #     @test v1.null_count == v2.null_count
    #     # @test v1.total_bytes == v2.total_bytes
    # end
    end
end


text:
type MongoCursor
    _wrap_::Ptr{Void}

    MongoCursor(_wrap_::Ptr{Void}) = begin
        cursor = new(_wrap_)
        finalizer(cursor, destroy)
        return cursor
    end
end
export MongoCursor

# Iterator

start(cursor::MongoCursor) = nothing
export start

next(cursor::MongoCursor, _::Void) =
    (BSONObject(ccall(
        (:mongoc_cursor_current, libmongoc),
        Ptr{Void}, (Ptr{Void},),
        cursor._wrap_
        ), Union{}), _)
export next

done(cursor::MongoCursor, _::Void) = begin
    return !ccall(
        (:mongoc_cursor_next, libmongoc),
        Bool, (Ptr{Void}, Ptr{Ptr{Void}}),
        cursor._wrap_,
        Array(Ptr{Void}, 1)
        )
end
export done

destroy(collection::MongoCursor) =
    ccall(
        (:mongoc_cursor_destroy, libmongoc),
        Void, (Ptr{Void},),
        collection._wrap_
        )


text:
#
# Given electron density in real space, return Hartree potential in reciprocal
# space
#
function Poisson_solve( PW::PWGrid, rhoR )
    #
    G2 = PW.gvec.G2
    Ns = PW.Ns
    Npoints = prod(Ns)
    #
    ctmp = 4.0*pi*R_to_G( Ns, rhoR )
    #
    ctmp[1] = 0.0
    for ip = 2:Npoints
        ctmp[ip] = ctmp[ip]/G2[ip]
    end
    return ctmp
end


text:
"""
ДАНО: Робот - в произвольной клетке ограниченного прямоугольного поля (без внутренних перегородок)
РЕЗУЛЬТАТ: Робот - в исходном положении, в клетке с роботом стоит маркер, и все остальные клетки поля промаркированы в шахматном порядке
"""
function chesser(r::Robot)
    first=0 #Счётчики шагов
    second=0
    side=Ost
    while (isborder(r,Sud)&&isborder(r,West))==false
        first=moves!(r,West)
        second=moves!(r,Sud)
    end
    if isodd(first+second)==true #Проверка на нечётность (от неё зависит порядок расстановки маркеров!)
        while isborder(r,Nord)==false
            lets_mark_then(r,side)
            side=inverse(side)
            move!(r,Nord)
        end
        lets_mark_then(r,side)
    else
        while isborder(r,Nord)==false
            lets_mark_now(r,side)
            side=inverse(side)
            move!(r,Nord)
        end
        lets_mark_now(r,side)
    end
    go_back_SW(r)
    while first>0
        move!(r,Ost)
        first-=1
    end
    while second>0
        move!(r,Nord)
        second-=1
    end
end
function moves!(r::Robot,side::HorizonSide) #Движение и подсчёт шагов
    num_steps=0
    while isborder(r,side)==false
        move!(r,side)
        num_steps=num_steps+1
    end
    return num_steps
end
function lets_mark_now(r::Robot,side::HorizonSide) #Маркировка с данной клетки
    while isborder(r,side)==false
        putmarker!(r)
        move!(r,side)
        if isborder(r,side)==false
            move!(r,side)
        else break
        end
    end
end
function lets_mark_then(r::Robot,side::HorizonSide) #Маркировка со следующей клетки
    while isborder(r,side)==false
        move!(r,side)
        putmarker!(r)
        if isborder(r,side)==false
            move!(r,side)
        else break
        end
    end
end
function go_back_SW(r::Robot) #Возвращение в ЮЗ угол
    while isborder(r,Sud)==false
        move!(r,Sud)
    end
    while isborder(r,West)==false
        move!(r,West)
    end
end
function inverse(side::HorizonSide) #Возможно подключение библиотечной функции
    HorizonSide(mod(Int(side)+2, 4))
end
#=Возможна другая реализация функций
и программы в целом=#

max_stars_repo_name: alexgruzberg/mirea-progs

max_issues_repo_name: alexgruzberg/mirea-progs

max_forks_repo_name: alexgruzberg/mirea-progs


text:
allapprox(x,y,tol=1e-8) = all(isapprox.(x,y,atol=tol))
FDM = central_fdm(5,1)

function kappa_AD(::Val{:Zygote},k::Kernel,d::Real)
    first(Zygote.gradient(x->kappa(k,x),d))
end

function kappa_AD(::Val{:ForwardDiff},k::Kernel,d::Real)
    first(ForwardDiff.gradient(x->kappa(k,first(x)),[d]))
end

function kappa_fdm(k::Kernel,d::Real)
    first(FiniteDifferences.grad(FDM,x->kappa(k,x),d))
end

function transform_AD(::Val{:Zygote},t::Transform,A)
    ps = KernelFunctions.params(t)
    @test allapprox(first(Zygote.gradient(p->transform_with_duplicate(p,t,A),ps)),
        first(FiniteDifferences.grad(FDM,p->transform_with_duplicate(p,t,A),ps)))
    @test allapprox(first(Zygote.gradient(X->sum(transform(t,X,2)),A)),
            first(FiniteDifferences.grad(FDM,X->sum(transform(t,X,2)),A)))
end

function transform_AD(::Val{:ForwardDiff},t::Transform,A)
    ps = KernelFunctions.params(t)
    if t isa ScaleTransform
        @test allapprox(first(ForwardDiff.gradient(p->transform_with_duplicate(first(p),t,A),[ps])),
            first(FiniteDifferences.grad(FDM,p->transform_with_duplicate(p,t,A),ps)))
    else
        @test allapprox(ForwardDiff.gradient(p->transform_with_duplicate(p,t,A),ps),
            first(FiniteDifferences.grad(FDM,p->transform_with_duplicate(p,t,A),ps)))
    end
    @test allapprox(ForwardDiff.gradient(X->sum(transform(t,X,2)),A),
            first(FiniteDifferences.grad(FDM,X->sum(transform(t,X,2)),A)))
end

transform_with_duplicate(p,t,A) = sum(transform(KernelFunctions.duplicate(t,p),A,2))


text:
using Test

using MachineGa

@testset "Validation of byte arrays" begin
    @testset "NaN value" begin
        # 0x255, 0x255, ..., 0x255
        bytes = map(x -> UInt8(255), 1:8)
        @test !ByteWorks.validate(bytes)
    end

    @testset "NaN multiple values" begin
        # 0x255, 0x255, ..., 0x255
        bytes = map(x -> UInt8(255), 1:32)
        @test !ByteWorks.validate(bytes)
    end

    @testset "Valid value" begin
        bytes = ByteWorks.floattobytes(3.14159265)
        @test ByteWorks.validate(bytes)
    end

    @testset "Multiple valid values" begin
        bytes = ByteWorks.floatstobytes([3.14159265, 2.71828, 1234.5])
        @test ByteWorks.validate(bytes)
    end
end

@testset "Check bytes and restore if values are not valid" begin
    bytes = map(x -> UInt8(255), 1:16)
    newbytes = ByteWorks.checkandrestore(
        bytes,
        [0.0, 0.0],
        [10.0, 10.0]
    )
    floats = ByteWorks.bytestofloats(newbytes)
    @test length(newbytes) == 16
    @test ByteWorks.validate(newbytes)
    @test all(x -> x <= 10.0 && x >= 0, floats)
end

@testset "Floats to bytes" begin
    @testset "Single float to bytes" begin
        bytes = [0x1f, 0x85, 0xeb, 0x51, 0xb8, 0x1e, 0x01, 0x40]
        result = ByteWorks.floattobytes(2.14)
        @test bytes == result
    end

    @testset "Float array to bytes" begin
        bytes = [
            0x1f,
            0x85,
            0xeb,
            0x51,
            0xb8,
            0x1e,
            0x09,
            0x40,
            0xae,
            0x47,
            0xe1,
            0x7a,
            0x14,
            0xae,
            0x05,
            0x40,
        ]
        result = ByteWorks.floatstobytes([3.14, 2.71])
        @test bytes == result
    end

    @testset "Length of bytes" begin
        arr = [1.0, 2.0, 3.0, 4.0]
        result = ByteWorks.floatstobytes(arr)
        @test length(result) == length(arr) * sizeof(Float64)
    end
end

@testset "Bytes to floats" begin
    @testset "Bytes to floats" begin
        bytes = [
            0x1f,
            0x85,
            0xeb,
            0x51,
            0xb8,
            0x1e,
            0x09,
            0x40,
            0xae,
            0x47,
            0xe1,
            0x7a,
            0x14,
            0xae,
            0x05,
            0x40,
        ]
        result = ByteWorks.bytestofloats(bytes)
        @test result == [3.14, 2.71]
    end
end

@testset "Near byte mutation" begin
    bytes = ByteWorks.floatstobytes([3.14, 2.71])
    mutated = RawOperators.nearbytemutation(bytes, mutateprob = 1.0)
    difference = bytes .- mutated
    # Since the mutation probability is 1.0
    # Vector of differences of original and mutated bytes
    # can not include zeros
    @test all(x -> x != 0, difference)
end

@testset "Random byte mutation" begin
    bytes = ByteWorks.floatstobytes([3.14, 2.71])
    mutated = RawOperators.randombytemutation(bytes, mutateprob = 1.0)
    difference = bytes .- mutated
    # Since the mutation probability is 1.0
    # Vector of differences of original and mutated bytes
    # can not include zeros
    @test all(x -> x != 0, difference)
end

@testset "Crossover on bytes" begin
    @testset "One-point crossover on bytes -> x(C, C) = C" begin
        bytes1 = [0x1, 0x2, 0x3, 0x4, 0x5, 0x6, 0x7, 0x8, 0x9, 0x10]
        child = RawOperators.onepointcrossover(bytes1, bytes1)
        @test child == bytes1
    end

    @testset "One-point crossover on bytes -> x(C, D)" begin
        bytes1 = [0x1, 0x2, 0x3, 0x4, 0x5, 0x6, 0x7, 0x8, 0x9, 0x0a]
        bytes2 = [0x0a, 0x09, 0x08, 0x7, 0x6, 0x5, 0x4, 0x3, 0x2, 0x1]
        child = RawOperators.onepointcrossover(bytes1, bytes2)
        @test map((par1, par2, x) -> x in [par1, par2], bytes1, bytes2, child) |> all
    end

    @testset "Uniform crossover on bytes" begin
        bytes1 = [0x1, 0x2, 0x3, 0x4, 0x5, 0x6, 0x7, 0x8, 0x9, 0x0a]
        bytes2 = [0x0a, 0x09, 0x08, 0x7, 0x6, 0x5, 0x4, 0x3, 0x2, 0x1]
        child = RawOperators.uniformcrossover(bytes1, bytes2)
        @test map((par1, par2, x) -> x in [par1, par2], bytes1, bytes2, child) |> all
    end
end

@testset "Mutation on bytes of floats" begin

    @testset "Near byte mutation" begin

        floats_original = [3.14, 2.71, 5.0]
        bytes_original = ByteWorks.floatstobytes(floats_original)

        floats_mutated = FloatOperators.nearbytemutation(floats_original, mutateprob = 1.0)
        bytes_mutated = ByteWorks.floatstobytes(floats_mutated)

        differences = bytes_original .- bytes_mutated .|> abs
        @test all(x -> x != 0, differences)
    end

    @testset "Random byte mutation" begin

        floats_original = [3.14, 2.71, 5.0]
        bytes_original = ByteWorks.floatstobytes(floats_original)

        floats_mutated =
            FloatOperators.randombytemutation(floats_original, mutateprob = 1.0)
        bytes_mutated = ByteWorks.floatstobytes(floats_mutated)

        differences = bytes_original .- bytes_mutated .|> abs
        @test all(x -> x != 0, differences)
    end
end

@testset "Byte crossover on floats" begin
    @testset "One-point crossover" begin
        floats = [1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0]
        floatsnew = FloatOperators.onepointcrossover(floats, floats)
        difference = floats .- floatsnew .|> abs
        @test all(x -> x == 0, difference)
    end

    @testset "Uniform crossover" begin
        floats = [1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0]
        floatsnew = FloatOperators.uniformcrossover(floats, floats)
        difference = floats .- floatsnew .|> abs
        @test all(x -> x == 0, difference)
    end

    @testset "Uniform crossover enhanced" begin
        floats1 = [1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0]
        floats2 = floats1 * -1.0

        floatsnew = FloatOperators.uniformcrossover(floats1, floats2)
        eitherfirstorsecond =
            map((x, a, b) -> x == a || x == b, floatsnew, floats1, floats2)
        @test all(eitherfirstorsecond)
    end
end

@testset "Genetic algorithm operators - Crossover" begin
    @testset "One-point crossover" begin
        c1 = Chromosome([1.0, 2.0], 0.0)
        c2 = Chromosome([2.0, 3.0], -1.0)
        cnew = GaOperators.onepointcrossover(c1, c2)
        @test cnew isa Chromosome
        @test isinf(cnew.cost)
        @test length(cnew.genes) == length(c1.genes)
    end

    @testset "Uniform-point crossover" begin
        c1 = Chromosome([1.0, 2.0], 0.0)
        c2 = Chromosome([2.0, 3.0], -1.0)
        cnew = GaOperators.uniformcrossover(c1, c2)
        @test cnew isa Chromosome
        @test isinf(cnew.cost)
        @test length(cnew.genes) == length(c1.genes)
    end

    @testset "Uniform-point crossover enhanced" begin
        c1 = Chromosome([1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0], 0.0)
        c2 = Chromosome(c1.genes * -1, 100.0)
        cnew = GaOperators.uniformcrossover(c1, c2)
        @test cnew isa Chromosome
        @test isinf(cnew.cost)
        @test length(cnew.genes) == length(c1.genes)
        @test all(map((x, a, b) -> x == a || x == b, cnew.genes, c1.genes, c2.genes))
    end

end

@testset "Genetic algorithm operators - Mutation" begin
    @testset "Random mutation" begin
        c = Chromosome([1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0], 0.0)
        cnew = GaOperators.randombytemutation(c, mutateprob = 1.0)
        @test cnew isa Chromosome
        @test isinf(cnew.cost)
        @test length(c.genes) == length(cnew.genes)
    end

    @testset "Inc Dec mutation" begin
        c = Chromosome([1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0], 0.0)
        cnew = GaOperators.nearbytemutation(c, mutateprob = 1.0)
        @test cnew isa Chromosome
        @test isinf(cnew.cost)
        @test length(c.genes) == length(cnew.genes)
    end
end

@testset "Genetic algorithm operators - Tournament selection" begin
    @testset "The worst can be selected" begin
        cs = [
            Chromosome([], 10),
            Chromosome([], 20),
            Chromosome([], 30),
            Chromosome([], 40),
            Chromosome([], 50),
        ]
        c = GaOperators.tournamentselection(cs, tournaments = 2)
        @test c isa Chromosome
        @test c.cost in [10, 20, 30, 40, 50]
        @test length(c.genes) == 0
    end

    @testset "The worst should not be selected" begin
        cs = [
            Chromosome([], 10),
            Chromosome([], 20),
            Chromosome([], 30),
            Chromosome([], 40),
            Chromosome([], 50),
        ]
        c = GaOperators.tournamentselection(cs, tournaments = 1000)
        @test c isa Chromosome
        @test c.cost in [10, 20, 30, 40]
        @test length(c.genes) == 0
    end
end


text:
using IntervalLinearAlgebra
using Plots
using Documenter

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

makedocs(;
    modules=[IntervalLinearAlgebra],
    authors="Luca Ferranti",
    repo="https://github.com/JuliaIntervals/IntervalLinearAlgebra.jl/blob/{commit}{path}#{line}",
    sitename="IntervalLinearAlgebra.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://juliaintervals.github.io/IntervalLinearAlgebra.jl",
        assets=String[],
        collapselevel=1,
    ),
    pages=[
        "Home" => "index.md",
        "Tutorials" => [
            "Linear systems" => "tutorials/linear_systems.md",
            "Eigenvalue computations" => "tutorials/eigenvalues.md"
        ],
        "Applications" => "wip.md",
        "Explanations" => [
            "Interval system solution set" => "explanations/solution_set.md",
            "Preconditioning" => "explanations/preconditioning.md"
        ],
        "API" => [
            "Interval matrices classification" => "api/classify.md",
            "Solver interface" => "api/solve.md",
            "Interval linear systems" => "api/algorithms.md",
            "Preconditioners" => "api/precondition.md",
            "Verified real linear systems" => "api/epsilon_inflation.md",
            "Eigenvalues" => "api/eigenvalues.md",
            "Miscellaneous" => "api/misc.md"
        ],
        "References" => "references.md",
        "Contributing" => "CONTRIBUTING.md"
    ],
)

deploydocs(;
    repo="github.com/JuliaIntervals/IntervalLinearAlgebra.jl", devbranch = "main"
)


text:
"""
    FieldMapping

# Constructors

```julia
FieldMapping(name::String, units::Union{String, Nothing})
FieldMapping(; name::String, units::Union{String, Nothing} = nothing)
```

# Arguments

- `name`: The name of field in the SQL database table.
- `units`: The units of the field values in the SQL database table.

Provides the information needed to map a field understood by the BiTWeather.jl
module functions to a field in the SQL database table.
"""
struct FieldMapping
    name::String
    units::Union{String, Nothing}
end

function FieldMapping(; name::String, units::Union{String, Nothing} = nothing)
    return FieldMapping(name, units)
end

"""
    Configuration

# Constructors

```julia
Configuration(dsn::String, table::String, fieldMappings::Dict{Symbol, FieldMapping})
Configuration(; dsn::String, table::String, fieldMappings::Dict{Symbol, FieldMapping})
```

# Arguments

- `dsn`: A Data Source Name (DSN) that can connect to the SQL database
    and access `table`.
- `table`: The SQL database table containing the weather data.
- `fieldMappings`: A dictionary of [`BiTWeather.FieldMapping`](@ref) entries, each
    of which provides a mapping from a field understood by the BiTWeather.jl module
    methods to a field in the SQL database table. The BiTWeather.jl module fields and
    allowed units are

    - `:dateTime`
    - `:pressure`: *`"kPa"`*, `"bar"`, `"inHg"`, `"kgm/s2"`, `"mbar"`, `"mmHg"`, `"N/m2"`, `"Pa"`, `"Pascals"`
    - `:relativeHumidity`: *`"percent"`*
    - `:solarRadiation`: *`"W/m2"`*
    - `:temperature`: *`"C"`*, *`"Celsius"`*, `"F"`, `"Fahrenheit"`, `"K"`, `"Kelvin"`
    - `:windDirection`: *`"degrees"`*
    - `:windSpeed`:  *`"m/s"`*, `"Knots"`, `"k/hr"`, `"mph"`

The mapping for `:dateTime` is required. The SI units to which all other units
are converted appear at the beginning of the list in italic.

See [`configuration_example`](@ref) for an example `Configuration`.
"""
struct Configuration
    dsn::String
    table::String
    fieldMappings::Dict{Symbol, FieldMapping}

    function Configuration(dsn, table, fieldMappings)

        for (fieldName, fieldMapping) in fieldMappings
            if fieldName != :dateTime && !(fieldName in keys(unitMappings))
                error("""Unknown field $(fieldName).""")
            end
            if fieldName != :dateTime
                if isnothing(fieldMapping.units)
                    error(""""units" must be set for $(fieldName)".""")
                end
                units = keys(unitMappings[fieldName])
                if !(fieldMapping.units in units)
                    error("""The units of $(fieldName) must be one of $(units).""")
                end
            end
        end

        return new(dsn, table, fieldMappings)
    end
end

function Configuration(; dsn::String, table::String, fieldMappings::Dict{Symbol, FieldMapping})::Configuration

    return Configuration(dsn, table, fieldMappings)
end

max_stars_repo_name: pebender/BiTWeather

max_issues_repo_name: pebender/BiTWeather

max_forks_repo_name: pebender/BiTWeather


text:
module TestDay02

import AOC2021.Day02.part1, AOC2021.Day02.part2
using Test

test_input_raw = """
forward 5
down 5
forward 8
up 3
down 8
forward 2"""

test_input = [string(x) for x in split(test_input_raw, "\n")]

@testset "Day 02" begin
    @testset "part 1" begin
        @test part1(test_input) == 150
    end
    @testset "part 2" begin
        @test part2(test_input) == 900
    end
end

end # module


text:
export every, fps, fpswhen, throttle, debounce

"""
```
debounce(dt, input, f=(acc,x)->x, init=value(input), reinit=x->x;
            typ=typeof(init), name=auto_name!(string("debounce ",dt,"s"), input))
```

Creates a signal that will delay updating until `dt` seconds have passed since the last time `input` has updated. By default, the debounce signal holds the last update of the `input` signal since the debounce signal last updated.

This behavior can be changed by the `f`, `init` and `reinit` arguments. The `init` and `f` functions are similar to `init` and `f` in `foldp`. `reinit` is called after the debounce sends an update, to reinitialize the initial value for accumulation, it gets one argument, the previous accumulated value.

For example
    `y = debounce(0.2, x, push!, Int[], _->Int[])`
will accumulate a vector of updates to the integer signal `x` and push it after `x` is inactive (doesn't update) for 0.2 seconds.

"""
function debounce(dt, node::Signal{T}, f=(acc,x)->x, init=value(node),
        reinit=x->x; typ=typeof(init), name=auto_name!("debounce $(dt)s", node)) where T
    # we don't add `node` as a parent of throttle, as the action is added to the
    # `node` itself, which pushes to the output node at the appropriate time.
    output = Signal(typ, init, (); name=name)
    throttle_connect(dt, output, node, f, init, reinit, false, true)
    output
end

"""
```
throttle(dt, input, f=(acc,x)->x, init=value(input), reinit=x->x;
            typ=typeof(init), name=auto_name!(string("throttle ",dt,"s"), input), leading=false)
```

Throttle a signal to update at most once every dt seconds. By default, the throttled signal holds the last update of the `input` signal during each `dt` second time window.

This behavior can be changed by the `f`, `init` and `reinit` arguments. The `init` and `f` functions are similar to `init` and `f` in `foldp`. `reinit` is called when a new throttle time window opens to reinitialize the initial value for accumulation, it gets one argument, the previous accumulated value.

For example
    `y = throttle(0.2, x, push!, Int[], _->Int[])`
will create vectors of updates to the integer signal `x` which occur within 0.2 second time windows.

If `leading` is `true`, the first update from `input` will be sent immediately by the throttle signal. If it is false, the first update will happen `dt` seconds after `input`'s first update

New in v0.4.1: `throttle`'s behaviour from previous versions is now available with the `debounce` signal type.

"""
function throttle(dt, node::Signal{T}, f=(acc,x)->x, init=value(node),
        reinit=x->x; typ=typeof(init), name=auto_name!("throttle $(dt)s", node),
        leading=false) where T
    output = Signal(typ, init, (node,); name=name)
    throttle_connect(dt, output, node, f, init, reinit, leading, false)
    output
end

# Aggregate a signal producing an update at most once in dt seconds
function throttle_connect(dt, output, input, f, init, reinit, leading, debounce)
    collected = init
    timer = Timer(identity, interval=0) #dummy timer to initialise
    dopush(_) = begin
        push!(output, collected)
        collected = reinit(collected)
        prevpush = time()
    end

    # we add the do_throttle as a foreach on the input, so when input updates it
    # collects the input values and pushes to the output when the time is right.
    prevpush = 0 # immediate push of `input`'s first update (unless leading is false)
    function do_throttle(inpval)
        collected = f(collected, inpval)
        prevpush == 0 && !leading && (prevpush = time())
        elapsed = time() - prevpush
        debounce && (elapsed = 0) # for debounce, only the timer can trigger a push

        close(timer)
        if elapsed > dt
            # prevpush is reset in dopush, so that calls via the Timer also reset it
            dopush(elapsed)
        else
            timer = Timer(dopush, interval=dt-elapsed)
        end
        nothing
    end
    foreach(do_throttle, input; init=nothing, name="$(input.name) throttler")
end

"""
    every(dt)

A signal that updates every `dt` seconds to the current timestamp. Consider using `fpswhen` or `fps` if you want specify the timing signal by frequency, rather than delay.
"""
function every(dt; name=auto_name!("every $dt secs"))
    n = Signal(time(), (); name=name)
    every_connect(dt, n)
    n
end

function every_connect(dt, output)
    outputref = WeakRef(output)
    function onerror_close_timer(pushnode, val, error_node, ex)
        print_error(pushnode, val, error_node, ex)
        close(timer)
    end
    timer = Timer(x -> push!(outputref.value, time(), onerror_close_timer), dt; interval=dt)
    finalizer(_->close(timer), output)
    output
end

"""
    fpswhen(switch, rate)

returns a signal which when `switch` signal is true, updates `rate` times every second. If `rate` is not possible to attain because of slowness in computing dependent signal values, the signal will self adjust to provide the best possible rate.
"""
function fpswhen(switch, rate; name=auto_name!("$rate fpswhen", switch))
    # Creates a node and sets up a timer that pushes to the node every 1.0/rate
    # seconds.
    n = Signal(Float64, 0.0, (switch, ); name=name)
    fpswhen_connect(rate, switch, n, name)
    n
end

function setup_next_tick(outputref, switchref, dt, wait_dt)
    Timer(t -> begin
        if value(switchref.value)
            push!(outputref.value, dt)
        end
    end, interval=wait_dt)
end

function fpswhen_connect(rate, switch, output, name)
    prev_time = time()
    dt = 1.0/rate
    outputref = WeakRef(output)
    switchref = WeakRef(switch)
    timer = Timer(identity, interval=0) # dummy timer to initialise
    function fpswhen_runner()
        # this function will run if switch gets a new value (i.e. is "active")
        # and if output is pushed to (assumed to be by the timer)
        if switch.value
            start_time = time()
            timer = setup_next_tick(outputref, switchref, start_time-prev_time, dt)
            prev_time = start_time
        else
            close(timer)
        end
        switch.value
    end
    # the fpswhen_aux will start and stop the timer if switch's value updates, it'll
    # also setup the next tick once the first tick is pushed to output
    fpswhen_aux = Signal(switch.value, (switch, output); name="fpswhen runner switch: $(switch.name), output: $(output.name)")
    preserve(fpswhen_aux)
    add_action!(fpswhen_runner, fpswhen_aux)

    fpswhen_runner() # init
    # ensure timer stops when output node is garbage collected
    finalizer(_->close(timer), output)
end

"""
    fps(rate)

Same as `fpswhen(Input(true), rate)`
"""
fps(rate; name="$rate fps") = fpswhen(Signal(Bool, true, (); name="fps true"), rate; name=name)


text:
include("../src/TcrCore.jl")
using TcrCore, HypothesisTests

function differentMechanisms()
    srand(1234)
    for acc in 0:5:100
        applicationOnly = mean([simSimple(1000, setupAgentsWithFixedAccuracy(100, acc), [benchmarkRegistryMean]) for i in 1:20])[1]
        challenge = mean([simChallenge(1000, setupAgentsWithFixedAccuracy(100, acc), [benchmarkRegistryMean]) for i in 1:20])[1]
        token = mean([simToken(1000, setupAgentsWithFixedAccuracy(100, acc), noneRedistribution, [benchmarkRegistryMean]) for i in 1:20])[1]
        println("$acc, $applicationOnly, $challenge, $token")
    end
end

function tokensAndDiversifications()
    srand(1234)

    for acc in 0:5:100
        noTokenFixed = mean([simChallenge(1000, setupAgentsWithFixedAccuracy(100, acc), [benchmarkRegistryMean]) for i in 1:30])[1]
        noTokenDiv = mean([simChallenge(1000, setupRandomAgents(100, acc, 20), [benchmarkRegistryMean]) for i in 1:30])[1]
        tokenFixed = mean([simToken(1000, setupAgentsWithFixedAccuracy(100, acc), noneRedistribution, [benchmarkRegistryMean]) for i in 1:30])[1]
        tokenDiv = mean([simToken(1000, setupRandomAgents(100, acc, 20), noneRedistribution, [benchmarkRegistryMean]) for i in 1:30])[1]
        println("$acc, $noTokenFixed, $noTokenDiv, $tokenFixed, $tokenDiv")
    end
end

function tokensAndDiversificationsEffects()
    srand(1234)
    divEffects = []
    tokenEffect = []
    tokenDivEffects = []
    for acc in 0:5:100
        noTokenFixed = mean([simChallenge(1000, setupAgentsWithFixedAccuracy(100, acc), [benchmarkRegistryMean]) for i in 1:30])[1]
        noTokenDiv = mean([simChallenge(1000, setupRandomAgents(100, acc, 20), [benchmarkRegistryMean]) for i in 1:30])[1]
        tokenFixed = mean([simToken(1000, setupAgentsWithFixedAccuracy(100, acc), noneRedistribution, [benchmarkRegistryMean]) for i in 1:30])[1]
        tokenDiv = mean([simToken(1000, setupRandomAgents(100, acc, 20), noneRedistribution, [benchmarkRegistryMean]) for i in 1:30])[1]

        divEffect = mean([tokenDiv - tokenFixed, noTokenDiv - noTokenFixed])
        tokenEffect = mean([tokenDiv - noTokenDiv, tokenFixed - noTokenFixed])
        tokenDivEffect = tokenDiv - noTokenFixed
        println("$acc, $divEffect, $tokenEffect, $tokenDivEffect")
    end
end

function compareEfficiencyByAccuracy()
    for acc in 0:5:100
        noToken = [simChallenge(1000, setupRandomAgents(100, acc, 20), [benchmarkRegistryMean])[1] for i in 1:10]
        token = [simToken(1000, setupRandomAgents(100, acc, 20), noneRedistribution, [benchmarkRegistryMean])[1] for i in 1:10]
        boost = mean(token) - mean(noToken)
        test = pvalue(UnequalVarianceTTest(token, noToken))
        significant = test < 0.05
        println("$acc, $boost, $test, $significant")
    end
end

function efficiencyArea()
    for noAgents in 10:10:200
        for noSteps in 100:100:2000
            noToken = [simChallenge(noSteps, setupRandomAgents(noAgents, 50, 20), [benchmarkRegistryMean])[1] for i in 1:10]
            token = [simToken(noSteps, setupRandomAgents(noAgents, 50, 20), noneRedistribution, [benchmarkRegistryMean])[1] for i in 1:10]
            boost = mean(token) - mean(noToken)
            println("$noSteps, $noAgents, $boost")
        end
    end
end

function expertsInnerCirleSize()
    for expertCount in 1:1:50
        controlLow= mean([simChallenge(1000, setupAgentsWithFixedAccuracy(100, 20), [benchmarkRegistryMean])[1] for i in 1:20])
        expertTokens = mean([simToken(1000, setupMixedAgents(100-expertCount, 20, expertCount, 100), noneRedistribution, [benchmarkRegistryMean])[1] for i in 1:20])
        expertNoTokens = mean([simChallenge(1000, setupMixedAgents(100-expertCount, 20, expertCount, 100), [benchmarkRegistryMean])[1] for i in 1:20])
        controlHigh = mean([simChallenge(1000, setupAgentsWithFixedAccuracy(100, 100), [benchmarkRegistryMean])[1] for i in 1:20])
        println("$expertCount, $controlLow, $expertTokens, $expertNoTokens, $controlHigh")
    end
end

function compareRedistribution()
    for acc in 0:5:100
        challengerOnly = [simToken(1000, setupRandomAgents(100, acc, 20), onlyChallengerRewardRedistribution, [benchmarkRegistryMean])[1] for i in 1:20]
        punishment = [simToken(1000, setupRandomAgents(100, acc, 20), punishmentRedistribution, [benchmarkRegistryMean])[1] for i in 1:20]
        reward = [simToken(1000, setupRandomAgents(100, acc, 20), rewardRedistribution, [benchmarkRegistryMean])[1] for i in 1:20]
        punishmentAndReward = [simToken(1000, setupRandomAgents(100, acc, 20), punishmentAndRewardRedistribution, [benchmarkRegistryMean])[1] for i in 1:20]
        pBoost = mean(punishment) - mean(challengerOnly)
        rBoost = mean(reward) - mean(challengerOnly)
        prBoost = mean(punishmentAndReward) - mean(challengerOnly)
        # test = pvalue(UnequalVarianceTTest(full, challengerOnly))
        # significant = test < 0.05
        println("$acc, $pBoost, $rBoost, $prBoost")
    end
end

function runScenarios()
    compareRedistribution()
end

runScenarios()