text:
using IteratedProcessSimulations
using Documenter
using VegaLite, UUIDs

# Snippet source: https://github.com/queryverse/VegaLite.jl/blob/master/docs/make.jl
function Base.show(io::IO, m::MIME"text/html", v::VegaLite.VLSpec)
    divid = string("vl", replace(string(uuid4()), "-" => ""))
    print(io, "

")
    print(io, "")
end
# End snippet

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

makedocs(;
    modules=[IteratedProcessSimulations],
    authors="Jeremiah Lewis <4462211+jlewis91@users.noreply.github.com> and contributors",
    repo="https://github.com/jeremiahpslewis/IteratedProcessSimulations.jl/blob/{commit}{path}#{line}",
    sitename="IteratedProcessSimulations.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://jeremiahpslewis.github.io/IteratedProcessSimulations.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
        "Recommender Vignette" => "Vignette-Recommender.md",
    ],
)

deploydocs(;
    repo="github.com/jeremiahpslewis/IteratedProcessSimulations.jl",
    devbranch="main",
)


text:
# # Driving on a Straight Roadway

#md # [![](https://img.shields.io/badge/show-nbviewer-579ACA.svg)](@__NBVIEWER_ROOT_URL__/notebooks/straight_roadway.ipynb)

# This example demonstrates a simple, one-dimensional driving simulation in which
# cars drive along a straight roadway.
# The vehicles are represented using the `Entity` type with the followin attributes:
#
# - `state` - `VehicleState`, containing the vehicle position and speed
# - `def` - `VehicleDef`, containing length, width, and class
# - `id` - `Int`, a unique label for each vehicle`
#
# A driving situation with different vehicles at a given time is referred to as a scene. It is represented by the `Scene` object.
# A `Scene` can be thought of as a vector of vehicles. However, in addition to simple vectors it allows to query vehicles by ID using the `get_by_id` function.
# We use a straight roadway with 1 lane as the environment.

using AutomotiveSimulator
using AutomotiveVisualization # for rendering
AutomotiveVisualization.colortheme["background"] = colorant"white"; # hide

roadway = gen_straight_roadway(1, 2000.0)  # 200m long straight roadway with 1 lane
scene = Scene([
    Entity(VehicleState(VecSE2(10.0,0.0,0.0), roadway, 8.0), VehicleDef(), 1),
    Entity(VehicleState(VecSE2(50.0,0.0,0.0), roadway, 12.5), VehicleDef(), 2),
    Entity(VehicleState(VecSE2(150.0,0.0,0.0), roadway, 6.0), VehicleDef(), 3),
])

veh_1 = get_by_id(scene, 1) # note that the order of the vehicles in the scene does not necessarily match the id

camera = StaticCamera(position=VecE2(100.0,0.0), zoom=4.75, canvas_height=100)
snapshot = render([roadway, scene], camera=camera)
#md write("straight_roadway.svg", snapshot) # hide
#md # ![three cars on road](straight_roadway.svg)

# In this call to the `render` function, we use the default rendering behavior for
# entities. More advanced examples will show how the rendering of entities can be customized.

# We can add an overlay that displays the car id:

idoverlay = IDOverlay(scene=scene, color=colorant"black", font_size=20, y_off=1.)
snapshot = render([roadway, scene, idoverlay], camera=camera)
#md write("straight_roadway_with_id.svg", snapshot) # hide
#md # ![three cars with id](straight_roadway_with_id.svg)

# To run a simulation we need driving models that produce actions.
# For this we will use `LaneFollowingDriver`s that produce `LaneFollowingAccel`s.
# For this demo, we will give each car a different model.

models = Dict{Int, LaneFollowingDriver}(
    1 => StaticLaneFollowingDriver(0.0), # always produce zero acceleration
    2 => IntelligentDriverModel(v_des=12.0), # default IDM with a desired speed of 12 m/s
    3 => PrincetonDriver(v_des = 10.0), # default Princeton driver with a desired speed of 10m/s
)

nticks = 100
timestep = 0.1
scenes = simulate(scene, roadway, models, nticks, timestep)
#md nothing # hide

# We can visualize the simulation as a sequence of images, for example using the
# `Reel` package

using Reel

animation = roll(fps=1.0/timestep, duration=nticks*timestep) do t, dt
    i = Int(floor(t/dt)) + 1
    idoverlay.scene = scenes[i]
    renderables = [roadway, scenes[i], idoverlay]
    render(renderables, camera=camera)
end
#md write("straight_roadway_animated.gif", animation) # hide
#md nothing # hide
#md # ![three vehicles animated](straight_roadway_animated.gif)

# In order to inspect the simulation interactively, we can use the `Interact` package

#md using Interact
#md using Blink
#md using ElectronDisplay
#md
#md w = Window()
#md viz = @manipulate for step in 1 : length(scenes)
#md     renderables = [roadway, scenes[step], IDOverlay(scene=scenes[step], roadway=roadway)]
#md     render(renderables, camera=camera)
#md end
#md body!(w, viz)

ext: jl


text:
using ParticleFilters
using Parameters
using Distributions
using QMDP
using POMDPPolicies
using StaticArrays
using POMDPModelTools

import POMDPs: max_possible_weight, new_particle, reset_distribution

export
    SimpleLightDark,
    DSimpleLightDark,
    LDHSolver,
    LDSide,
    LDHeuristic,
    LDSidePolicy,

    max_possible_weight,
    new_particle,

    ContinuousLightDark,
    CLDPingPong,
    CLDPingPongPolicy,
    CLDInformationGain,
    CLDGoalPolicy,
    CLDHeuristic

include("simple_ld.jl")
include("continuous_ld.jl")


text:
@sumProductRule(:node_type     => GaussianWeightedMeanPrecision,
                :outbound_type => Message{GaussianWeightedMeanPrecision},
                :inbound_types => (Nothing, Message{PointMass}, Message{PointMass}),
                :name          => SPGaussianWeightedMeanPrecisionOutNPP)

@naiveVariationalRule(:node_type     => GaussianWeightedMeanPrecision,
                      :outbound_type => Message{GaussianWeightedMeanPrecision},
                      :inbound_types => (Nothing, ProbabilityDistribution, ProbabilityDistribution),
                      :name          => VBGaussianWeightedMeanPrecisionOut)


text:
using LinearAlgebra, SparseArrays, Plots
if 1==1
    if 1==1
        ρ = 0.05; δ = 0.10; A = 0.20;
        r(s,a) = A*s -a -0.5*(a)^2.0
        μ(s,a) = a - δ*s
        dr(s,a) = -1 -a
        FOC(s,Vs) = Vs - 1
        μ_inv(s,ṡ) = ṡ + δ*s
        s_ss = (A-ρ-δ)/(δ*(ρ+δ))
        a_ss = δ*s_ss
        v_ss = r(s_ss,a_ss)/ρ        # ∫exp(-ρt)*r(s_ss,a_ss) dt = r(s_ss,a_ss)/ρ
        #s_min = 0.00*s_ss
        s_min = -.5
        s_max = 2.000*s_ss
        # verify things are well defined @ corners
        s_max, μ_inv(s_max,0), dr(s_max,μ_inv(s_max,0))
        s_min, μ_inv(s_min,0), dr(s_min,μ_inv(s_min,0))
        H = 10_000;
        s = collect(LinRange(s_min, s_max, H))
        ds = (s_max-s_min)/(H-1)
        dVf, dVb            = [zeros(H,1) for i in 1:2]
        dV_Upwind, a_Upwind = [zeros(H,1) for i in 1:2]
        dVf_end       = dr(s_max,μ_inv(s_max,0))
        dVb_1         = dr(s_min,μ_inv(s_min,0))
        v0 = v_ss *ones(H)
        v0 = @. r(s, μ_inv(s,0))/ρ #initial guess for V
        v = v0
    end
    Δ = 1_000; maxit = 10_000;ε = 10e-8; dist=[];

    # Generic HJB Solver
    # Parimonious: dV_Upwind = dVf.*If + dVb.*Ib + dV0.*I0
    for n=1:maxit
        V=v
        dV = (V[2:end] - V[1:end-1])/ds

        # forward difference: if ṡ>0
        dVf = [dV; dVf_end]
        af  = FOC.(s,dVf)                 # choice w forward difference
        μ_f = μ.(s, af)                   # ṡ      w forward difference
        If  = μ_f .> 0

        # backward difference: if ṡ<0
        dVb = [dVb_1; dV]
        ab  = FOC.(s,dVb)                          # choice w backward difference
        μ_b = μ.(s, ab)
        Ib  = μ_b .< 0

        # neither difference: if ṡ=0
        a0  = μ_inv.(s,0)        # c if ṡ=0
        dV0 = dr.(s,a0)
        μ_0 = μ.(s, a0)          # μ_0 == zero(s)
        I0  = (1.0 .- If - Ib)   # choice betw forward & backward difference

        # I_concave = dVb .> dVf
        # scatter(I_concave) #1 everywhere EXCEPT @ last point H.

        dV_Upwind = dVf.*If + dVb.*Ib + dV0.*I0
        a_Upwind  = FOC.(s,dV_Upwind)
        μ_Upwind  = μ.(s, a_Upwind)
        u = r.(s,a_Upwind)

        # create the transition matrix AA
        X = -min.(μ_b,0)/ds
        Y = -max.(μ_f,0)/ds + min.(μ_b,0)/ds
        Z = max.(μ_f,0)/ds

        a1 = sparse(Diagonal((Y[:])))
        a2 = [zeros(1,H); sparse(Diagonal((X[2:H]))) zeros(H-1,1)]
        a3 = [zeros(H-1,1) sparse(Diagonal((Z[1:H-1]))); zeros(1,H)]
        AA = a1 + a2 + a3

        # Solve for new value
        B = (ρ + 1/Δ)*sparse(I,H,H) - AA
        b = u + V./Δ
        V = B \ b

        # 8: stopping criteria
        V_change = V-v
        v = V

        push!(dist,findmax(abs.(V_change))[1])
        println(n, " ", dist[n])
        if dist[n] .< ε
            println("Value Function Converged Iteration=")
            println(n)
            break
        end
    end
    s_dot = μ.(s, a_Upwind)
    v_err = r.(s, a_Upwind) + dV_Upwind.*s_dot - ρ.*v # approx @ borrowing constraint
    plot(s, v, xlabel="s", ylabel="V(s)",legend=false, title="")
    plot!([s_ss],  seriestype = :vline, lab="", color="grey", l=:dash)
    plot!([v_ss],  seriestype = :hline, lab="", color="grey", l=:dash)

    # Simulate.
    using Interpolations
    â = LinearInterpolation(s, a_Upwind[:], extrapolation_bc = Line())
    Δt = 0.01; T = 150; time = 0.0:Δt:T
    s_sim, a_sim, ṡ_sim = [zeros(length(time),1) for i in 1:3]
    s_0 =0.5*s_ss
    s_sim[1] = s_0
    for i in 2:length(time)
        a_sim[i-1] = â(s_sim[i-1])
        ṡ_sim[i-1] = μ(s_sim[i-1], a_sim[i-1])
        s_sim[i] = s_sim[i-1] + Δt * ṡ_sim[i-1]
        # s_sim[i] = s_sim[i-1] + Δt * μ(s_sim[i-1], a_sim[i-1])
    end
    ix = 1:(length(time)-1)
    plot()
    plot!(time[ix], s_sim[ix], lab = "k")
    plot!(time[ix], a_sim[ix], lab = "i")
    plot!(time[ix], ṡ_sim[ix], lab = "k̇")
    plot!([s_min],  seriestype = :hline, lab="", color="grey", l=:dash)
    plot!([s_ss a_ss],  seriestype = :hline, lab="", color="grey", l=:dash)
end


text:
#!/usr/bin/env julia

using GZip

"""
Remove variants from exonic (coding) or masked low resolution mapping regions of the genome
"""

wigf = ARGS[1]   # .wig file
bedf = ARGS[2]   # .bed file

function read(file)
    if ismatch(r".gz$", file)   # Gzipped .gz files
		f = GZip.open(file)
		lines = readlines(f)
		close(f)
    else
		f = open(file)
		lines = readlines(f)
		close(f)
	end
	return lines
end

"""
Expanding the dashed intervals into single nucleotide positions
example:
Y:2688071-2688073 =>
                    Y:2688071
                    Y:2688072
                    Y:2688073
"""
function expand(locci)
    list = []
    pos = split(locci, ':')
        chr = pos[1]
        region = split(pos[2], '-')
        n = parse(Int64, region[1])
        while n <= parse(Int64, region[end])
             out = "$chr" * ":" * "$n"
            push!(list, out)
            n +=  1
        end
    return list
end

"""
Read .wig file

track graphType=bar windowingFunction=none autoScale=off viewLimits=0:1
variableStep chrom=Y
2649597	1
2649610	1
2649664	1
"""
function readwig(fh)
    pos = []
    lines = read(fh)
    header = lines[1:2]
    lines  = lines[3:end]
    chr = split(chomp(header[2]), '=')
    for line in lines
        cols = split(line)
        out = chr[2] * ":" * cols[1]
        push!(pos, out)
    end
    chr = chr[2]
    unique(pos)
    return chr, header, pos
end

"""
Read .bed file

chr1	67000041	67000051	NM_032291_cds_0_0_chr1_67000042_f	0	+
chr1	67091529	67091593	NM_032291_cds_1_0_chr1_67091530_f	0	+
chr1	67098752	67098777	NM_032291_cds_2_0_chr1_67098753_f	0	+
"""
function readBED(bed_file, chr::AbstractString=chr)
    pos = []
    bed = read(bed_file)
    for line in bed[5265:5999]
        bed = split(line, '\t')
        chrm = replace(bed[1], "chr", "")
        if chrm == chr
            start, last = bed[2], bed[3]
            chroms = "$chrm" * ":" * "$start" * "-" * "$last"
            region = expand(chroms)
            append!(pos, region)
        end
    end
    pos = unique(pos)
    return pos
end

function printWIG(wigf, bedf)
    chr, header, ln = readwig(wigf)
    bed = readBED(bedf, chr)
    ncDNA = setdiff(Set(ln), Set(bed))

    for line in header
        println(chomp(line))
    end

   for i in ncDNA
       pattern = chr * ":"
       println(replace(i, pattern, ""), ' ', 1)
   end
end

printWIG(wigf, bedf)


text:
function KronGraph500NoPerm(SCALE,EdgePerVertex)

    # Set the number of vertices and edges to be generated
    N = 2^SCALE
    M = round(EdgePerVertex * N)

    # Set R_MAT (2x2 Kronecker) coefficients
    A = 0.57
    B = 0.19
    C = 0.19
    D = 1-(A+B+C)

    # Initialize index arrays
    ij = ones(Int,2,M)

    # Normalize coefficients
    ab = A+B
    c_norm = C/(1-(A+B))
    a_norm = A/(A+B)

    # Loop over each scale
    for ib = 1:SCALE
        ii_bit = rand(1,M) .> ab
        jj_bit = rand(1,M) .> (c_norm * ii_bit + a_norm * .!(ii_bit))
        ij = ij + 2^(ib-1) * [ii_bit; jj_bit]
    end

    StartVertex = ij[1,:]
    EndVertex = ij[2,:]

    return StartVertex, EndVertex
end


text:
import Base

function Base.haskey(m::RegexMatch, name::Symbol)
    idx = Base.PCRE.substring_number_from_name(m.regex.regex, name)
    return idx > 0
end
Base.haskey(m::RegexMatch, name::AbstractString) = haskey(m, Symbol(name))


text:
export @snoopl

using Serialization

"""
```
@snoopl "func_names.csv" "llvm_timings.yaml" begin
    # Commands to execute, in a new process
end
```
causes the julia compiler to log timing information for LLVM optimization during the
provided commands to the files "func_names.csv" and "llvm_timings.yaml". These files can
be used for the input to `SnoopCompile.read_snoopl("func_names.csv", "llvm_timings.yaml")`.

The logs contain the amount of time spent optimizing each "llvm module", and information
about each module, where a module is a collection of functions being optimized together.
"""
macro snoopl(flags, func_file, llvm_file, commands)
    return :(snoopl($(esc(flags)), $(esc(func_file)), $(esc(llvm_file)), $(QuoteNode(commands))))
end
macro snoopl(func_file, llvm_file, commands)
    return :(snoopl(String[], $(esc(func_file)), $(esc(llvm_file)), $(QuoteNode(commands))))
end

function snoopl(flags, func_file, llvm_file, commands)
    println("Launching new julia process to run commands...")
    # addprocs will run the unmodified version of julia, so we
    # launch it as a command.
    code_object = """
            using Serialization
            while !eof(stdin)
                Core.eval(Main, deserialize(stdin))
            end
            """
    process = open(`$(Base.julia_cmd()) $flags --eval $code_object`, stdout, write=true)
    serialize(process, quote
        let func_io = open($func_file, "w"), llvm_io = open($llvm_file, "w")
            ccall(:jl_dump_emitted_mi_name, Nothing, (Ptr{Nothing},), func_io.handle)
            ccall(:jl_dump_llvm_opt, Nothing, (Ptr{Nothing},), llvm_io.handle)
            try
                $commands
            finally
                ccall(:jl_dump_emitted_mi_name, Nothing, (Ptr{Nothing},), C_NULL)
                ccall(:jl_dump_llvm_opt, Nothing, (Ptr{Nothing},), C_NULL)
                close(func_io)
                close(llvm_io)
            end
        end
        exit()
    end)
    wait(process)
    println("done.")
    nothing
end


text:
using Test
using DFTK
using Logging
import DFTK: mpi_sum

include("testcases.jl")
using Random
Random.seed!(0)

function test_consistency_term(term; rtol=1e-4, atol=1e-8, ε=1e-6, kgrid=[1, 2, 3],
                               kshift=[0, 1, 0]/2, lattice=silicon.lattice,
                               Ecut=10, spin_polarization=:none)
    sspol = spin_polarization != :none ? " ($spin_polarization)" : ""
    xc    = term isa Xc ? "($(first(term.functionals)))" : ""
    @testset "$(typeof(term))$xc $sspol" begin
        n_dim = 3 - count(iszero, eachcol(lattice))
        Si = n_dim == 3 ? ElementPsp(14, psp=load_psp(silicon.psp)) : ElementCoulomb(:Si)
        atoms = [Si, Si]
        model = Model(lattice, atoms, silicon.positions; n_electrons=silicon.n_electrons,
                      terms=[term], spin_polarization, symmetries=true)
        basis = PlaneWaveBasis(model; Ecut, kgrid, kshift)

        n_electrons = silicon.n_electrons
        n_bands = div(n_electrons, 2, RoundUp)
        filled_occ = DFTK.filled_occupation(model)

        ψ = [Matrix(qr(randn(ComplexF64, length(G_vectors(basis, kpt)), n_bands)).Q)
             for kpt in basis.kpoints]
        occupation = [filled_occ * rand(n_bands) for _ in 1:length(basis.kpoints)]
        occ_scaling = n_electrons / sum(sum(occupation))
        occupation = [occ * occ_scaling for occ in occupation]
        ρ = with_logger(NullLogger()) do
            compute_density(basis, ψ, occupation)
        end
        E0, ham = energy_hamiltonian(basis, ψ, occupation; ρ=ρ)

        @assert length(basis.terms) == 1

        δψ = [randn(ComplexF64, size(ψ[ik])) for ik = 1:length(basis.kpoints)]
        function compute_E(ε)
            ψ_trial = ψ .+ ε .* δψ
            ρ_trial = with_logger(NullLogger()) do
                compute_density(basis, ψ_trial, occupation)
            end
            E, _ = energy_hamiltonian(basis, ψ_trial, occupation; ρ=ρ_trial)
            E.total
        end

        diff = (compute_E(ε) - compute_E(-ε)) / (2ε)

        diff_predicted = 0.0
        for ik in 1:length(basis.kpoints)
            Hψ = ham.blocks[ik]*ψ[ik]
            δψHψ = sum(occupation[ik][iband] * real(dot(δψ[ik][:, iband], Hψ[:, iband]))
                       for iband=1:n_bands)
            diff_predicted += 2 * basis.kweights[ik] * δψHψ
        end
        diff_predicted = mpi_sum(diff_predicted, basis.comm_kpts)

        err = abs(diff - diff_predicted)
        @test err < rtol * abs(E0.total) || err < atol
    end
end

@testset "Hamiltonian consistency" begin
    test_consistency_term(Kinetic())
    test_consistency_term(AtomicLocal())
    test_consistency_term(AtomicNonlocal())
    test_consistency_term(ExternalFromReal(X -> cos(X[1])))
    test_consistency_term(ExternalFromFourier(X -> randn()))
    test_consistency_term(LocalNonlinearity(ρ -> ρ^2))
    test_consistency_term(Hartree())
    test_consistency_term(Ewald())
    test_consistency_term(PspCorrection())
    test_consistency_term(Xc(:lda_xc_teter93))
    test_consistency_term(Xc(:lda_xc_teter93), spin_polarization=:collinear)
    test_consistency_term(Xc(:gga_x_pbe), spin_polarization=:collinear)
    test_consistency_term(Xc(:mgga_x_tpss))
    test_consistency_term(Xc(:mgga_x_scan))
    test_consistency_term(Xc(:mgga_c_scan), spin_polarization=:collinear)
    test_consistency_term(Xc(:mgga_x_b00))
    test_consistency_term(Xc(:mgga_c_b94), spin_polarization=:collinear)

    let
        a = 6
        pot(x, y, z) = (x - a/2)^2 + (y - a/2)^2
        Apot(x, y, z) = .2 * [y - a/2, -(x - a/2), 0]
        Apot(X) = Apot(X...)
        test_consistency_term(Magnetic(Apot); kgrid=[1, 1, 1], kshift=[0, 0, 0],
                              lattice=[a 0 0; 0 a 0; 0 0 0], Ecut=20)
        test_consistency_term(DFTK.Anyonic(2, 3.2); kgrid=[1, 1, 1], kshift=[0, 0, 0],
                              lattice=[a 0 0; 0 a 0; 0 0 0], Ecut=20)
    end
end


text:
module ForestBiometrics

const K = 0.005454154
const KMETRIC = 0.00007854

struct Metric end
const metric = Metric()

include("tree.jl")
include("stand.jl")
include("density.jl")
include("charts.jl")
include("height.jl")
include("expansion.jl")

export Tree
export Stand

export basal_area, basal_area!
export trees_per_area, trees_per_area!
export qmd, qmd!
export sdi

export gingrich_chart
export reineke_chart

export HeightModel
export estimate_height, estimate_height!

export curtis,
michailoff,
meyer,
micment,
micment2,
naslund,
naslund2,
naslund3,
naslund4,
power,
wyckoff

#3 parameter equations, mainly from LMFOR R package
#sorted alphabetically
export chapman,
gompertz,
hossfeldIV,
korf,
logistic,
monserud,
prodan,
ratkowsky,
sibbesen,
weibull

export /,
collapse, collapse!,
expand, expand!

end


text:
module OSMPBF
include("proto/OSMPBF_pb.jl")
include("model.jl")
include("scan.jl")
include("load.jl")

export scan_pbf, load_pbf, Node, Way, Relation, RelationMemberType
end


text:
module SauterSchwabQuadrature

using LinearAlgebra
using CompScienceMeshes
using StaticArrays

export sauterschwab_parameterized
export SauterSchwabStrategy, CommonFace, CommonEdge, CommonVertex, PositiveDistance
export generate_integrand_uv

"""
Compute interaction integrals using the quadrature introduced in [1].

    sauterschwab_parameterized(integrand, strategy)

Here, `integrand` is the pull-back of the integrand into the barycentric
parameter domain of the two triangles that define the integration domain. Note
that here we use (for a planar triangle) the representation:

    x = x[3] + u*(x[1]-x[3]) + v*(x[2]-x[3])

with `u` ranging from 0 to 1 and `v` ranging from 0 to 1-`u`. This parameter
domain and representation is different from the one used in [1].

The second argument 'strategy' is an object whose type is one of

- `CommonFace`
- `CommonEdge`
- `CommonVertex`
- `PositiveDistance`

according to the configuration of the two triangular patches defining the
domain of integration. The constructors of these classes take a single
argument `acc` that defines the number of quadrature points along each of the
four axes of the final rectangular (ξ,η) integration domain (see [1], Ch 5).

[1] Sauter. Schwwab, 'Boundary Element Methods', Springer Berlin Heidelberg, 2011
"""
function sauterschwab_parameterized end

include("pulled_back_integrands.jl")
include("sauterschwabintegral.jl")
include("reorder_vertices.jl")
include("parametric_kernel_generator.jl")

end


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

using Markdown
using InteractiveUtils

# ╔═╡ 63ba08cc-59a8-11eb-0a0f-27efac60d779
using Pkg, DrWatson

# ╔═╡ 6db218c6-59a8-11eb-2a8b-7107354cf590
begin
	#@quickactivate "StatisticalRethinkingStan"
	using StanSample, StanOptimize
	using StatisticalRethinking
end

# ╔═╡ 51fc19b8-59a8-11eb-2214-15aca59b807b
md" ## Figure 8.2s"

# ╔═╡ 8aaa4bcc-59a8-11eb-2003-f1213b116565
begin
	df = CSV.read(sr_datadir("rugged.csv"), DataFrame)
	df_africa = df[df.cont_africa .== 1, [:rgdppc_2000, :rugged]]
	dropmissing!(df_africa, :rgdppc_2000)
	dropmissing!(df_africa, :rugged)
	df_africa.log_gdp = log.(df_africa[:, :rgdppc_2000])
	scale!(df_africa, [:log_gdp, :rugged])
	PRECIS(df_africa)
end

# ╔═╡ 25ee6cd8-6a54-11eb-0a72-3fb09ee63b0e
begin
	df_non_africa = df[df.cont_africa .== 0, [:rgdppc_2000, :rugged]]
	dropmissing!(df_non_africa, :rgdppc_2000)
	dropmissing!(df_non_africa, :rugged)
	df_non_africa.log_gdp = log.(df_non_africa[:, :rgdppc_2000])
	scale!(df_non_africa, [:log_gdp, :rugged])
end;

# ╔═╡ d7d3e626-6a45-11eb-1820-a3ec98e556b1
stan8_0 = "
data {
	int N;
	vector[N] G;
	vector[N] R;
}

parameters {
	real a;
	real b;
	real<lower=0> sigma;
}

transformed parameters {
	vector[N] mu;
	mu = a + b * (R - 0.125);
}

model {
	a ~ normal(1, 1);
	b ~ normal(0, 1);
	sigma ~ exponential(1);
	G ~ normal(mu, sigma);
}
";

# ╔═╡ a3e7c070-6a46-11eb-072b-3943db854020
begin
	data1 = (N = size(df_africa, 1), G = df_africa.log_gdp_s,
		R = df_africa.rugged_s)
	m8_1s = SampleModel("m8.1s", stan8_0)
	rc8_1_1s = stan_sample(m8_1s; data=data1)
	if success(rc8_1_1s)
		post8_1_1s_df = read_samples(m8_1s, :dataframe)
		PRECIS(post8_1_1s_df[:, [:a, :b, :sigma]])
	end
end

# ╔═╡ cc11715a-6a54-11eb-1955-d38a021a3bb3
begin
	data2 = (N = size(df_non_africa, 1), G = df_non_africa.log_gdp_s,
		R = df_non_africa.rugged_s)
	rc8_1_2s = stan_sample(m8_1s; data=data2)
	if success(rc8_1_2s)
		post8_1_2s_df = read_samples(m8_1s, :dataframe)
		PRECIS(post8_1_2s_df[:, [:a, :b, :sigma]])
	end
end

# ╔═╡ 1968b688-6a49-11eb-2dec-9986830b4a7e
begin
	p1 = plotbounds(
		df_africa, :rugged, :log_gdp,
		post8_1_1s_df, [:a, :b, :sigma];
		bounds=[:none, :hpdi],
		colors=[:orange, :lightblue],
		title="African nations",
		xlab="ruggedness",
		ylab="log GDP"
	)

	df_afr = df[df.cont_africa .== 1, [:country, :rgdppc_2000, :rugged]]
	df_afr = df_afr[df_afr.rugged .> 4, :]
	for (ind, country) in enumerate(df_afr.country)
		annotate!([([df_afr.rugged[ind]+0.3], [log(df_afr.rgdppc_2000[ind])+0.15],
			Plots.text(df_afr.country[ind], 6, :red, :right))])
	end
end

# ╔═╡ 0f36fb46-6a5d-11eb-2d34-91266e770493
begin
	p2 = plotbounds(
		df_non_africa, :rugged, :log_gdp,
		post8_1_2s_df, [:a, :b, :sigma];
		bounds=[:none, :hpdi],
		colors=[:orange, :lightblue],
		title="Non-African nations",
		xlab="ruggedness",
		ylab="log GDP"
	)

	df_na = df[:, [:country, :rgdppc_2000, :rugged]]
	dropmissing!(df_na, :rgdppc_2000)
	dropmissing!(df_na, :rugged)
	df_na = df_na[df_na.rugged .> 4, :]
	for (ind, country) in enumerate(df_na.country)
		println(country)
		if !(country in df_afr.country)
			annotate!([([df_na.rugged[ind]+0.3],
						[log(df_na.rgdppc_2000[ind])+0.15],
				Plots.text(df_na.country[ind], 6, :red, :right))])
		end
	end
	plot(p1, p2, layout=(1,2))
end

# ╔═╡ 45767e2e-6a63-11eb-3e12-354a7e32a374
md" ## End of figure 8.2s"

# ╔═╡ Cell order:
# ╟─51fc19b8-59a8-11eb-2214-15aca59b807b
# ╠═63ba08cc-59a8-11eb-0a0f-27efac60d779
# ╠═6db218c6-59a8-11eb-2a8b-7107354cf590
# ╠═8aaa4bcc-59a8-11eb-2003-f1213b116565
# ╠═25ee6cd8-6a54-11eb-0a72-3fb09ee63b0e
# ╠═d7d3e626-6a45-11eb-1820-a3ec98e556b1
# ╠═a3e7c070-6a46-11eb-072b-3943db854020
# ╠═cc11715a-6a54-11eb-1955-d38a021a3bb3
# ╠═1968b688-6a49-11eb-2dec-9986830b4a7e
# ╠═0f36fb46-6a5d-11eb-2d34-91266e770493
# ╟─45767e2e-6a63-11eb-3e12-354a7e32a374


text:
# =================================================================
# Quadrotor model
# See https://easychair.org/publications/paper/gjfh
# =================================================================

using Reachability, MathematicalSystems, TaylorIntegration

# parameters of the model
const g = 9.81           # gravity constant in m/s^2
const R = 0.1            # radius of center mass in m
const l = 0.5            # distance of motors to center mass in m
const Mrotor = 0.1       # motor mass in kg
const M = 1.0            # center mass in kg
const m = M + 4*Mrotor   # total mass in kg
const mg = m*g

# moments of inertia
const Jx = (2/5)*M*R^2 + 2*l^2*Mrotor
const Jy = Jx
const Jz = (2/5)*M*R^2 + 4*l^2*Mrotor
const Cyzx = (Jy - Jz)/Jx
const Czxy = (Jz - Jx)/Jy
const Cxyz = 0.0 #(Jx - Jy)/Jz

# considering the control parameters as *parameters*
const u₁ = 1.0
const u₂ = 0.0
const u₃ = 0.0

@inline function quad_property(t, x)
    b1 = (x[3] < 1.4)
    b2 = t ≥ 1.0 ? (x[3] > 0.9) : true
    b3 = t ≥ 5.0 ? (0.98 ≤ x[3] ≤ 1.02) : true
    return b1 && b2 && b3
end

@taylorize function quadrotor!(dx, x, params, t)
    # unwrap the variables and the controllers; the last three are the controllers
    # x₁, x₂, x₃, x₄, x₅, x₆, x₇, x₈, x₉, x₁₀, x₁₁, x₁₂, u₁, u₂, u₃ = x
    x₁  = x[1]
    x₂  = x[2]
    x₃  = x[3]
    x₄  = x[4]
    x₅  = x[5]
    x₆  = x[6]
    x₇  = x[7]
    x₈  = x[8]
    x₉  = x[9]
    x₁₀ = x[10]
    x₁₁ = x[11]
    x₁₂ = x[12]

    # equations of the controllers
    F = (mg - 10*(x₃ - u₁)) + 3*x₆  # height control
    τϕ = -(x₇ - u₂) - x₁₀            # roll control
    τθ = -(x₈ - u₃) - x₁₁            # pitch control
    local τψ = 0.0                   # heading is uncontrolled
    #
    Tx = τϕ/Jx
    Ty = τθ/Jy
    Tz = τψ/Jz
    F_m = F/m

    # Some abbreviations
    sx7 = sin(x₇)
    cx7 = cos(x₇)
    sx8 = sin(x₈)
    cx8 = cos(x₈)
    sx9 = sin(x₉)
    cx9 = cos(x₉)
    #
    sx7sx9 = sx7*sx9
    sx7cx9 = sx7*cx9
    cx7sx9 = cx7*sx9
    cx7cx9 = cx7*cx9
    sx7cx8 = sx7*cx8
    cx7cx8 = cx7*cx8
    sx7_cx8 = sx7/cx8
    cx7_cx8 = cx7/cx8
    #
    x4cx8 = cx8*x₄
    #
    p11 = sx7_cx8*x₁₁
    p12 = cx7_cx8*x₁₂
    xdot9 = p11 + p12

    # differential equations for the quadrotor
    #
    dx[1] = (cx9*x4cx8 + (sx7cx9*sx8 - cx7sx9)*x₅) + (cx7cx9*sx8 + sx7sx9)*x₆
    dx[2] = (sx9*x4cx8 + (sx7sx9*sx8 + cx7cx9)*x₅) + (cx7sx9*sx8 - sx7cx9)*x₆
    dx[3] = (sx8*x₄ - sx7cx8*x₅) - cx7cx8*x₆
    dx[4] = (x₁₂*x₅ - x₁₁*x₆) - g*sx8
    dx[5] = (x₁₀*x₆ - x₁₂*x₄) + g*sx7cx8
    dx[6] = (x₁₁*x₄ - x₁₀*x₅) + (g*cx7cx8 - F_m)
    dx[7] = x₁₀ + sx8*xdot9
    dx[8] = cx7*x₁₁ - sx7*x₁₂
    dx[9] = xdot9
    dx[10] = Cyzx * (x₁₁ * x₁₂) + Tx
    dx[11] = Czxy * (x₁₀ * x₁₂) + Ty
    dx[12] = Cxyz * (x₁₀ * x₁₁) + Tz
     #
    return dx
end

function quad(; T=5.0, plot_vars=[0, 3],
                property=quad_property,
                project_reachset=true)

    # equations, x' = f(x(t))
    𝐹 = BlackBoxContinuousSystem(quadrotor!, 12)

    # initial conditions
    Wpos = 0.4
    Wvel = 0.4
    X0c = zeros(12)
    ΔX0 = [Wpos, Wpos, Wpos, Wvel, Wvel, Wvel, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
    X0 = Hyperrectangle(X0c, ΔX0)

    # instantiate the IVP
    𝑃 = InitialValueProblem(𝐹, X0)

    # general options
    𝑂 = Options(:T=>T, :plot_vars=>plot_vars, :property=>property,
                :project_reachset=>project_reachset, :mode=>"check")

    return (𝑃, 𝑂)
end

max_stars_repo_name: JuliaReach/ReachabilityBenchmarks

max_issues_repo_name: JuliaReach/ReachabilityBenchmarks

max_forks_repo_name: JuliaReach/ReachabilityBenchmarks


text:
module PTGen
    export PTVer
    include("Types.jl")
    include("Functions.jl")
end


text:
const indent = 4

mutable struct Logger
    level::Int
    verbose::Bool
    Logger(verbose) = new(0, verbose)
end

function print_message(l::Logger, msg, eq::AbstractArray; color=:green)
    if !l.verbose return end
    print(repeat(" ", l.level * indent))
    print(msg, ": ")

    printstyled("{"; color=:red)
    for i = 1:length(eq)-1
        printstyled(eq[i]; color)
        printstyled(", "; color=:red)
    end
    printstyled(eq[end]; color)
    printstyled("}\n"; color=:red)
end

function print_message(l::Logger, msg, eq; color=:red)
    if !l.verbose return end
    print(repeat(" ", l.level * indent))
    print(msg, ": ")
    printstyled(eq, '\n'; color)
end

function print_message(l::Logger, msg)
    if !l.verbose return end
    print(repeat(" ", l.level * indent))
    print(msg, '\n')
end

function print_message(l::Logger, msg, solved, unsolved; color1=:green, color2=:red)
    if !l.verbose return end
    print(repeat(" ", l.level * indent))
    print(msg, ": ")
    printstyled(solved; color=color1)
    if !isequal(unsolved, 0)
        print(" + ")
        printstyled("∫ ", unsolved; color=color2)
    end
    println()
end

##############################################################################

function attempt(l::Logger, msg, eq)
    print_message(l, msg, eq)
    l.level += 1
end

function attempt(l::Logger, msg)
    print_message(l, msg)
    l.level += 1
end

function inform(l::Logger, msg, eq)
    print_message(l, msg, eq)
end

function inform(l::Logger, msg)
    print_message(l, msg)
end

function result(l::Logger, msg, eq)
    if l.level > 0 l.level -= 1 end
    print_message(l, msg, eq)
end

function result(l::Logger, msg)
    if l.level > 0 l.level -= 1 end
    print_message(l, msg)
end

function result(l::Logger)
    if l.level > 0 l.level -= 1 end
end

function result(l::Logger, msg, solved, unsolved)
    if l.level > 0 l.level -= 1 end
    print_message(l, msg, solved, unsolved)
end


text:
py"""
class DiffDelage(torch.nn.Module):
    '''
    DiffDelage

Differentiable Delage model.

Atributes:
- `D::Int`: dimension (i.e. number of assets)
- `w::PyObject`: differentiable weights
- `c::Vector{PyObject}`: convex constraints on weights
- `μ::Array{Float64,1}`: predicted mean vector
- `sqrtΣ::Array{Float64,2}`: predicted upper triangular Cholesky factor of covariance matrix
- `solver_args::Dict=Dict("use_indirect" => true, "eps" => 1e-5, "max_iters" => 20000}`: solver arguments
    '''
    def __init__(
        self,
        D,
        w,
        c,
        μ,
        sqrtΣ,
        solver_args={'use_indirect': True, 'eps': 1e-5, 'max_iters': 20000}
    ):
        super(DiffDelage, self).__init__()
        self.D = D
        self.μ = μ
        Σ = sqrtΣ.T @ sqrtΣ
        self.Σ = Σ
        self.solver_args=solver_args

        P = cp.Variable((D+1, D+1), PSD=True)
        Q = cp.Variable((D+1, D+1), PSD=True)

        params = cp.Parameter(2)

        constraints = [
            P[D, 0:D] == -Q[D, 0:D] + w / 2 - Q[0:D, 0:D] @ μ,
            Q[0:D, 0:D] >> 0,
        ]
        objective = cp.Minimize(params[1] * cp.sum(cp.multiply(Σ, Q[0:D, 0:D])) - cp.quad_form(μ, Q[0:D, 0:D]) + Q[D, D] + cp.sum(cp.multiply(Σ, P[0:D, 0:D])) - 2 * μ.T @ P[D, 0:D] + params[0] * P[D, D])
        problem = cp.Problem(objective, constraints + c)

        self.cvxpylayer = CvxpyLayer(problem, parameters=[params], variables=[w, P, Q])

    def forward(self, params):
        return self.cvxpylayer(params, solver_args=self.solver_args)[0]

    def get_variables(self, params):
        params_ = torch.tensor(params, dtype=torch.float64, requires_grad=True)
        return self.cvxpylayer(params_, solver_args=self.solver_args)

    def get_w(self, params):
        params_ = torch.tensor(params, dtype=torch.float64, requires_grad=True)
        return self.get_variables(params_)[0]

    def get_lagrangian(self, params):
        params_ = torch.tensor(params, dtype=torch.float64, requires_grad=True)
        D = self.D
        μ_ = torch.tensor(self.μ, dtype=torch.float64)
        Σ_ = torch.tensor(self.Σ, dtype=torch.float64)
        w_, P_, Q_ = self.cvxpylayer(params_, solver_args=self.solver_args)
        return params_[1] * torch.sum(Σ_ * Q_[0:D, 0:D]) - torch.dot(μ_, Q_[0:D, 0:D] @ μ_) + Q_[D, D] + torch.sum(Σ_ * P_[0:D, 0:D]) - 2 * μ_.T @ P_[D, 0:D] + params_[0] * P_[D, D]
"""

DiffDelage = py"DiffDelage"

max_stars_repo_name: molet/DiffDRO

max_issues_repo_name: molet/DiffDRO

max_forks_repo_name: molet/DiffDRO


text:
#
# load and analyze Kaggle photo quality prediction data
#
using Ladder
using Losses
using Gadfly
using DataFrames
using Distributions
using DataFrames

# Required files:
# location of directory containing all submissions in .csv format
all_submissions_dir = "../photo/ALLSUBS"
# location of csv file containing lines of form: "team","submission name" for all submissions
team_file = "../photo/teams.csv"
# location of csv file containing solution labels
solution_file = "../photo/PhotoSolution.csv"

function sub2labels(sub)
    # Parse labels in file sub
    f = open(sub)
    lines = readlines(f)
    close(f)
    # dealing with missing line breaks
    if length(lines) == 1
        lines = split(lines[1],"\r")
    end
    # removing header line if present
    if lines[1][1:2] == "id" || length(lines) == 12001
        lines = lines[2:end]
    end
    if length(lines) != 12000
        # fixes some submissions that include training labels
        lines = lines[1:12000]
    end
    labels = Float64[]
    for line in lines
        # fixes two common formatting errors
        line = replace(line,"\"","")
        line = replace(line,";",",")
        tokens = split(line,",")
        #handles submissions with only one column
        if length(tokens)==1
            # extra string() conversion fixed julia bug in float
            push!(labels,float(string(tokens[1],"")))
        else
            push!(labels,float(string(tokens[2],"")))
        end
    end
    return labels
end

function privpart(solution)
    # parse which labels are private
    f = open(solution)
    lines = readlines(f)
    close(f)
    lines = lines[2:end]
    private = Int[]
    for line in lines
        if split(line,",")[3][1:2]=="Pr"
            push!(private,1)
        else
            push!(private,0)
        end
    end
    return private
end

function build_teams()
    # Assumes team.csv contains lines of form: team,submission name
    # for all submissions
    team2id = (String => Array{String,1})[]
    id2team = (String => String)[]
    teams = readcsv(team_file)
    t = size(teams)[1]
    for i in 1:t
        team2id[teams[i,1]] = String[]
    end
    for i in 1:t
        push!(team2id[teams[i,1]],idfromname(teams[i,2]))
        id2team[idfromname(teams[i,2])] = teams[i,1]
    end
    return team2id,id2team
end

idfromname(name) = split(name,' ')[1]

function build_allsubs()
    allsubs = Submission[]
    badsubs = Submission[]
    csvfiles = readdir(all_submissions_dir)
    for sub in csvfiles
        id = idfromname(sub)
        try
            predictions = sub2labels(string(all_submissions_dir,"/",sub))
            # test if submission is valid by scoring it
            empirical_loss(solution[private],predictions[private],cliplogloss)
            empirical_loss(solution[public],predictions[public],cliplogloss)
            push!(allsubs,Submission(id,id2team[id],predictions))
        catch
            push!(badsubs,Submission(id,id2team[id],Float64[]))
        end
    end
    return allsubs,badsubs
end

function lb(privs=private,pubs=public,p=kaggle_precision,lossfct=cliplogloss)
    # Kaggle leaderboard per team is equivalent to running fixed step size ladder
    # for each team with Kaggle accuracy
    prl = (String => LadderState)[]
    pbl = (String => LadderState)[]
    for team in keys(team2id)
        prl[team] = ladder(solution[privs],10.0^(-p))
        pbl[team] = ladder(solution[pubs],10.0^(-p))
    end
    for sub in allsubs
        score!(prl[sub.team],Submission(sub.id,sub.team,sub.labels[privs]),lossfct)
        score!(pbl[sub.team],Submission(sub.id,sub.team,sub.labels[pubs]),lossfct)
    end
    return prl,pbl
end

function lbpf(privs=private,pubs=public,lossfct=cliplogloss)
    #
    # Evaluate parameter free ladder, using one instance per team
    #
    prl = (String => LadderState)[]
    pbl = (String => LadderState)[]
    for team in keys(team2id)
        prl[team] = ladder(solution[privs])
        pbl[team] = ladder(solution[pubs])
    end
    for sub in allsubs
        score!(prl[sub.team],Submission(sub.id,sub.team,sub.labels[privs]),lossfct)
        score!(pbl[sub.team],Submission(sub.id,sub.team,sub.labels[pubs]),lossfct)
    end
    return prl,pbl
end

function lbdiff(prl::Dict{String,LadderState},pbl::Dict{String,LadderState})
    prm = (Float64,String)[]
    pbm = (Float64,String)[]
    for team in keys(team2id)
        push!(prm,(prl[team].current,team))
        push!(pbm,(pbl[team].current,team))
    end
    return prm,pbm
end

#
# Experiments
#

function standard_splits()
    println("Plot: Score distribution")
    # compute kaggle leaderboards for each team
    prl, pbl = lb()
    # compute best sub for each team
    prm, pbm = lbdiff(prl,pbl)
    # compute leaderboards according to parameter-free ladder
    prlpf,pblpf = lbpf()
    prmpf, pbmpf = lbdiff(prlpf,pblpf)

    function plotscores(df,title="Deviation of scores")
        return Gadfly.plot(df, x="x", y="y", Geom.point,
                           Coord.Cartesian(ymin=0.182,ymax=0.202,xmin=0.182,xmax=0.202),
                           Guide.yticks(ticks=[0.19,0.20]),
                           Guide.xticks(ticks=[0.18,0.19,0.20]),
                           Guide.xlabel("Public score"), Guide.ylabel("Private score"),
                           Guide.title(title),
                           Theme(minor_label_font_size=15pt,major_label_font_size=15pt,
                                 key_label_font_size=14pt,key_title_font_size=15pt),
    layer(x=[0.182,0.202],y=[0.182,0.202],Geom.line,Theme(default_color=color("red"))))
    end

    top50 = sort(prm)[1:50]
    top50priv = [x for (x,_) in top50]
    top50pub = [pbl[t].current for (_,t) in top50]
    top50pf = sort(prmpf)[1:50]
    top50pfpriv = [x for (x,_) in top50pf]
    top50pfpub = [pblpf[t].current for (_,t) in top50pf]
    df1 = DataFrame(x=top50pub,y=top50priv)
    df2 = DataFrame(x=top50pfpub,y=top50pfpriv)
    top50kaggle = plotscores(df1,"Top 50 Kaggle scores on private data")
    top50ladder = plotscores(df2,"Top 50 Ladder scores on private data")
    draw(PDF(string("top50kaggle-private.pdf"), 5inch, 4inch),top50kaggle)
    draw(PDF(string("top50ladder-private.pdf"), 5inch, 4inch),top50ladder)
    writetable(string("top50kaggle-private.csv"),df1)
    writetable(string("top50ladder-private.csv"),df2)

    println("Plot: Top 50 score comparison")
    dftop50test = DataFrame(x=[1:50],y=top50priv,label="Test")
    dftop50kaggle = DataFrame(x=[1:50],y=top50pub,label="Kaggle")
    dftop50ladder = DataFrame(x=[1:50],y=top50pfpub,label="Ladder")
    dftop50 = vcat(dftop50test,dftop50kaggle,dftop50ladder)
    top50plot = Gadfly.plot(dftop50, x="x", y="y", Geom.point,
                           Coord.Cartesian(ymin=0.182,ymax=0.202),
                            Geom.line,color="label",Scale.color_discrete_manual("green","red","blue"),
                           Guide.yticks(ticks=[0.19,0.20]),
                           Guide.xticks(ticks=[1,10,20,30,40,50]),
                           Guide.xlabel("Reported score"), Guide.ylabel("Rank by test score"),
                           Guide.title("Top 50 submissions by test score"),
                           Theme(minor_label_font_size=15pt,major_label_font_size=15pt,
                                 key_label_font_size=14pt,key_title_font_size=15pt))
    draw(PDF(string("publicscore-top50.pdf"), 5inch, 4inch),top50plot)
    writetable(string("publicscore-top50.csv"),dftop50)

    println("Plot: Significance analysis")
    #function ttest(m1,s1,m2,s2)
    function ttest(n,m,s)
        # one-sided Welch's t-test
        t = m/(s/sqrt(n))
        T = TDist(n-1)
        return 2.0* min(cdf(T,t),1-cdf(T,t))
    end

    function pvalues()
        ranking = sort(pbm)
        n = 8400
        k = 50
        losses = Array{Float64,1}[]
        pvals = Float64[]
        for i in 1:k
            labels = subs[pbl[ranking[i][2]].leader.id].labels
            push!(losses,map(cliplogloss,solution[private],labels[private]))
        end
        for i in 1:k
            push!(pvals,ttest(n,mean(losses[1]-losses[i]),std(losses[1]-losses[i])))
        end
        return pvals
    end

    function plotvals(pvals)
        df = DataFrame(x=[1:length(pvals)],y=pvals)
        pvalplot = Gadfly.plot(df, x="x", y="y", Geom.point,
                           Guide.yticks(ticks=[0.05,0.2,0.4,0.6,0.8]),
                           Guide.xticks(ticks=[2:length(pvals)]),
                           Coord.Cartesian(ymin=0.0,ymax=1.0),
                           Guide.xlabel("rank"), Guide.ylabel("p value"),
                           Theme(minor_label_font_size=15pt,major_label_font_size=15pt,
                                 key_label_font_size=14pt,key_title_font_size=15pt),
                                layer(x=[2,length(pvals)],y=[0.05,0.05],Geom.line,Theme(default_color=color("red"))))
        return df,pvalplot
    end

    pvals = pvalues()
    pvaldf,pvalplot = plotvals(pvals[1:10])
    draw(PDF(string("top10pvals.pdf"), 5inch, 4inch),pvalplot)
    writetable(string("top10pvals.csv"),pvaldf)
    # bonferroni correction for 9 comparisons
    pvalsbon = pvals[1:10]*9
    pvaldfbon,pvalplotbon = plotvals(pvalsbon)
    draw(PDF(string("top10pvalsbon.pdf"), 5inch, 4inch),pvalplotbon)
    writetable(string("top10pvalsbon.csv"),pvaldfbon)
end

function fresh_splits()
    #
    # Plot 3: Score deviations on new splits
    #
    println("Plot: Score distribution on fresh splits")

    function newsplit()
        rand(1:n,n)
        randperm = shuffle([1:n][private])
        split1 = bool(zeros(n))
        split2 = bool(zeros(n))
        split1[randperm[1:4200]]=true
        split2[randperm[4201:end]]=true

        # compute kaggle leaderboards for each team
        prl, pbl = lb(split1,split2)
        # compute best sub for each team
        prm, pbm = lbdiff(prl,pbl)
        # compute leaderboards according to parameter-free ladder
        prlpf,pblpf = lbpf(split1,split2)
        prmpf, pbmpf = lbdiff(prlpf,pblpf)

        top50 = sort(prm)[1:50]
        top50priv = [x for (x,_) in top50]
        top50pub = [pbl[t].current for (_,t) in top50]

        top50pf = sort(prmpf)[1:50]
        top50pfpriv = [x for (x,_) in top50pf]
        top50pfpub = [pblpf[t].current for (_,t) in top50pf]

        df1 = DataFrame(x=top50pub,y=top50priv)
        df2 = DataFrame(x=top50pfpub,y=top50pfpriv)

        return df1,df2
    end

    function plotscoresbars(means,stds,title="Top 50 scores on fresh split")
        return plot(x=means, y=means, ymin=means-stds,ymax=means+stds,Geom.point,Geom.errorbar,
                           Coord.Cartesian(ymin=0.182,ymax=0.202,xmin=0.182,xmax=0.202),
                           Guide.yticks(ticks=[0.19,0.20]),
                           Guide.xticks(ticks=[0.18,0.19,0.20]),
                           Guide.xlabel("Public score"), Guide.ylabel("Private score"),
                           Guide.title(title),
                           Theme(minor_label_font_size=15pt,major_label_font_size=15pt,
                                 key_label_font_size=13pt,key_title_font_size=15pt),
    layer(x=[0.182,0.202],y=[0.182,0.202],Geom.line,Theme(default_color=color("red"))))
    end

    # compute statistics across 20 independent iterations
    dfs = [newsplit() for _ in 1:20]
    dfs1 = [ x[1] for x in dfs]
    dfs2 = [ x[2] for x in dfs]
    scores1 = [ float([dfs1[i][j,1] for i in 1:20]) for j in 1:50]
    means1 = float([mean(scores1[i]) for i in 1:50])
    stds1 = float([std(scores1[i]) for i in 1:50])
    scores2 = [ float([dfs2[i][j,1] for i in 1:20]) for j in 1:50]
    means2 = float([mean(scores2[i]) for i in 1:50])
    stds2 = float([std(scores2[i]) for i in 1:50])
    plotbars1 = plotscoresbars(means1,stds1,"Kaggle top 50 scores on fresh splits")
    plotbars2 = plotscoresbars(means2,stds2,"Ladder top 50 scores on fresh splits")
    draw(PDF(string("top50kaggle-fresh.pdf"), 6inch, 4.5inch),plotbars1)
    draw(PDF(string("top50ladder-fresh.pdf"), 6inch, 4.5inch),plotbars2)
end

# process all teams and submission
solution = sub2labels(solution_file)
n = length(solution)
private = bool(privpart(solution_file))
# Kaggle used 5 bits of precision
kaggle_precision = 5
public = ~private

println("Building teams.")
team2id, id2team = build_teams()
teams = [ team for team in keys(team2id)]
println("Building submissions.")
allsubs, badsubs = build_allsubs()
subs = (String => Submission)[]
for sub in allsubs
    subs[sub.id] = sub
end
println("Done.")


text:
function setup(backend::GPUBackend, layer::GaussianKLLossLayer, inputs::Vector{Blob}, diffs::Vector{Blob})
  data_type = eltype(inputs[1])
  aux_ones = make_blob(backend, data_type, size(inputs[2])...) # for summing
  fill!(aux_ones, 1.0)
  tmp = make_blob(backend, data_type, size(inputs[2])...)
  state = GaussianKLLossLayerState(layer, zero(data_type), zero(data_type), 0,
                                   @compat(Dict(:aux_ones => aux_ones, :tmp => tmp)))
  return state
end

function forward(backend::GPUBackend, state::GaussianKLLossLayerState, inputs::Vector{Blob})
  mu  = inputs[1]
  sigma = inputs[2]
  num = get_num(mu)
  data_type = eltype(mu)
  n = length(mu) # length or num?

  aux_ones  = state.tmp_blobs[:aux_ones]
  log_sigma_tmp = state.tmp_blobs[:tmp]
  copy!(log_sigma_tmp, sigma)
  CuVec.log!(backend, log_sigma_tmp)

  Σμ²  = CuBLAS.dot(backend.cublas_ctx, data_type, n, mu.ptr, 1, mu.ptr, 1)
  Σσ²  = CuBLAS.dot(backend.cublas_ctx, data_type, n, sigma.ptr, 1, sigma.ptr, 1)
  logΣ = CuBLAS.dot(backend.cublas_ctx, data_type, n, log_sigma_tmp.ptr, 1, aux_ones.ptr, 1)
  state.loss = 0.5(Σμ² + Σσ² - 2logΣ - n) * state.layer.weight / num

  # accumulate statistics
  state.loss_accum *= state.n_accum
  state.loss_accum += state.loss * n
  state.loss_accum /= state.n_accum + n

  state.n_accum += n
end

function backward(backend::GPUBackend, state::GaussianKLLossLayerState, inputs::Vector{Blob}, diffs::Vector{Blob})

  mu  = inputs[1]
  sigma = inputs[2]
  data_type = eltype(mu)
  num = get_num(mu)

  # diff = df/dmu[i] = mu
  diff = diffs[1]
  if isa(diff, CuTensorBlob)
    copy!(diff, mu)
    CuVec.mul_scal!(backend, diff, convert(data_type, state.layer.weight/num))
  end

  # diff = df/dsigma[i] = sigma - 1/sigma
  diff = diffs[2]

  tmp = state.tmp_blobs[:tmp]
  copy!(tmp, sigma)
  if isa(diff, CuTensorBlob)
    copy!(diff, sigma)
    copy!(tmp, sigma)
    CuVec.pow!(backend, tmp, convert(data_type, -1.0))
    CuVec.sub!(backend, diff, tmp)
    CuVec.mul_scal!(backend, diff, convert(data_type, state.layer.weight/num))
  end
end


text:
using Parsers
export print_smart, parse_operator_string

"""
Wrapper for julia's base display() function
to make it more user friendly.
"""
function print_smart(M::AbstractVecOrMat, header::String="Data")
    println(header * ": ")
    display(M)
    println()
end

"""
Parses strings of operators seperated by spaces into a
list of tuples. creation operator denoted with '^' is
represented by '1' and annhilation operators are '0'
Ex: "3^ 1^ 2 4" => [(4, 0), (2, 0), (1, 1), (3, 1)]
"""
function parse_operator_string(s::String)
    sites = split(s, ' ', keepempty=false)

    operators = Tuple[]

    # create array of tuples
    # Second value denotes annhiliate/create: 0/1
    for num in sites
        if occursin("^", num)
            pushfirst!(operators, (parse(Int, chop(num)), 1))
        else
            pushfirst!(operators, (parse(Int, num), 0))
        end
    end

    return operators
end


text:
# This file was generated, do not modify it. # hide
res = evaluate!(multi; opts...)
round(res.measurement[1], sigdigits=3)


text:
"""
More general version of hypergeometric2F1 function than is provided by GSL library, based on 'Computation of Hypergeometric Functions' (Pearson 2009)
"""
function hypergeom_2F1_expanded(a::G, b::G, c::G, z::G) where G <: AbstractFloat
    if z < -1.0
        return (1.0 - z) ^ (-1.0 * a) * hypergeom_2F1_expanded(a, c - b, c, z / (z - 1.0))
    elseif z < 0.5
        return hypergeom([a, b], c, z)
    elseif z < 1.0
        coef1 = gamma(c) * gamma(c - a - b) / (gamma(c - a) * gamma(c - b))
        hypergeom1 = hypergeom_2F1_expanded(a, b, a + b - c + 1.0, 1.0 - z)
        coef2 = (1.0 - z) ^ (c - a - b) * gamma(c) * gamma(a + b - c) / (gamma(a) * gamma(b))
        hypergeom2 = hypergeom_2F1_expanded(c - a, c - b, c - a - b + 1.0, 1.0 - z)
        return coef1 * hypergeom1 + coef2 * hypergeom2
    else #if z >= 1.0
        coef1 = z ^ (-1.0 * a) * gamma(c) * gamma(c - a - b) / (gamma(c - a) * gamma(c - b))
        hypergeom1 = hypergeom_2F1_expanded(a, a - c + 1.0, a + b - c + 1.0, 1.0 - 1.0 / z)
        coef2 = z ^ (a - c) * (1.0 - z) ^ (c - a - b) * gamma(c) * gamma(a + b - c) / (gamma(a) * gamma(b))
        hypergeom2 = hypergeom_2F1_expanded(c - a, 1.0 - a, c - a - b + 1.0, 1.0 - 1.0 / z)
        return coef1 * hypergeom1 + coef2 * hypergeom2
    end
end

"""
Dark matter Density profile given in units of ρs
"""
function density_profile(x::G) where G <: AbstractFloat
    return (x * (1.0 + x) ^ 2) ^ -1
end

function density_profile(x::G, p::P) where {G <: AbstractFloat, P <: SFWParameters}
    return (x ^ -p.γ) * (1.0 + x ^ p.α) ^ ((p.γ - p.β) / p.α)
end

"""
Define gravitationla potential function, results are in units of Φs (i.e. are often multiplied by it)
"""
function gravitational_potential(x::G) where G <: AbstractFloat
    return 1.0 - (log(1.0 + x) / x)
end

function gravitational_potential(x::G, p::NFWParameters{G}) where G <: AbstractFloat
    return p.Φs * (1.0 - (log(1.0 + x) / x))
end

"""
Slower version of function gravitational_potential0() but works in cases where a analytic form for solution does not exist
"""
function gravitational_potential_integral(x::G, p::P) where {G <: AbstractFloat, P <: SFWParameters}
    function integrand1(xi::G)
        return  density_profile(xi, p) * xi * xi
    end
    function integrand2(xi::G)
        return density_profile(xi, p) * xi
    end
    Φ1 = quadgk(integrand1, 0.0, x)[1]
    Φ2 = quadgk(integrand2, x, Inf)[1]

    return p.Φs * (1.0 - (Φ1 / x + Φ2))
end

"""
Gravitational potential for SFW profile where center of galaxy is not forced to zero, use to compute normalization factor Φ0
"""
function gravitational_potential0(x::G, p::P) where {G <: AbstractFloat, P <: SFWParameters}
    if ((3.0 - p.β) / p.α % 1.0) == 0.0 || ((p.γ - p.β) / p.α % 1.0) == 0.0
        return gravitational_potential_integral(x, p)
    end
    p1b = hypergeom_2F1_expanded((3.0 - p.γ) / p.α, (p.β - p.γ) / p.α, (3.0 + p.α - p.γ) / p.α, -x ^ p.α)
    I1 = -(x ^ (3.0 - p.γ) * p1b) / (x * (p.γ - 3.0))
    p2 = hypergeom_2F1_expanded((p.β - 2.0) / p.α, (p.β - p.γ) / p.α, (p.α + p.β - 2.0) / p.α, -x ^ (-p.α))
    I2 = x ^ (2.0 - p.β) * p2 / (p.β - 2.0)
    return p.Φs * (1.0 - (I1 + I2))
end

#complete if re-encouter the error with hypergeom...
function gravitational_potential0_robust(x::G, p::P) where {G <: AbstractFloat, P <: SFWParameters}
    return try
        gravitational_potential0(x, p) #determine if should be 0 or not
    catch
        gravitational_potential_integral(x, p)
    end
end

"""
Gravitational potential for SFW profile where center of galaxy is forced to have a potential of 0.
"""
function gravitational_potential(x::G, p::P) where {G <: AbstractFloat, P <: SFWParameters}
    return gravitational_potential0(x, p) -  p.Φ0
end


text:
let
    gw = SimpleGridWorld(size=(2,2))

    pvec = fill(:left, length(states(gw)))

    solver = VectorSolver(pvec)

    p = solve(solver, gw)

    # test default
    @test sprint(showpolicy, gw, p) == " [1, 1] -> :left\n [2, 1] -> :left\n [1, 2] -> :left\n [2, 2] -> :left\n [-1, -1] -> :left"

    # test with small display
    iob = IOBuffer()
    io = IOContext(iob, :limit=>true, :displaysize=>(7, 7))
    showpolicy(io, gw, p, pre="@ ")
    @test String(take!(iob)) == "@ [1, 1] -> :left\n@ [2, 1] -> :left\n@ [1, 2] -> :left\n@ …"

    # test very long policy with small display
    struct M <: MDP{Int, Int}
        n::Int
    end
    iob = IOBuffer()
    io = IOContext(iob, :limit=>true, :displaysize=>(7, 7))
    POMDPs.states(m::MDP) = 1:m.n
    POMDPs.actions(m::MDP) = 1:3
    m = M(1_000_000_000)
    showpolicy(io, m, RandomPolicy(m))
    # Below, the actual values could be different because of the RandomPolicy, but length should be the same
    @test length(String(take!(iob))) == length(" 1 -> 2\n 2 -> 1\n 3 -> 3\n …")
end

max_stars_repo_name: JuliaPOMDP/POMDPPolicies

max_issues_repo_name: JuliaPOMDP/POMDPPolicies

max_forks_repo_name: JuliaPOMDP/POMDPPolicies


text:
#-----------------------------------------------------------
# Wrapper for ensuring that in case function call fails,
# a default value (e.g. Inf) is returned
#-----------------------------------------------------------

function safewrapper(f; E = PosDefException, inform::Bool=false, retvalue = Inf)

    function safe_f(x)

        try

            f(x)

        catch err

            if isa(err, E)

                inform ? @warn("Caught exception!") : nothing

                return retvalue

            else

                throw(err)

            end

        end

    end

end


text:
"""
    FixedLoad

Fixed load, a.k.a. Uncontrollable load
"""
struct FixedLoad <: Resource
    index::Int           # Index of the resource
    num_timesteps::Int   # Number of time-steps in the battery's operation

    load::Vector{Float64}     # Fixed load at each time step

    function FixedLoad(;index::Int=0, T::Int=0, load::Vector{Float64}=Float64[])
        T == length(load) || throw(DimensionMismatch(
            "T=$T but load has size $(length(load))"
        ))
        return new(index, T, load)
    end
end

function add_resource_to_model!(
    m::MOI.ModelLike,
    h::House,
    l::FixedLoad,
    var2idx, con2idx
)
    T = l.num_timesteps

    # ==========================================
    #    I. Add local variables
    # ==========================================
    # Create variables in model
    pload = MOI.add_variables(m, T)
    # Set variable bounds
    # Indices of variable bounds are not tracked
    for t in 1:T
        MOI.add_constraint(m, MOI.SingleVariable(pload[t]), MOI.EqualTo(l.load[t]))
    end

    # Update list of variable indices
    for t in 1:T
        var2idx[(:fixed, l.index, :pnet, t)] = pload[t]  # net load
    end

    # ==========================================
    #    II. Update linking constraints
    # ==========================================
    for t in 1:T
        cidx = con2idx[(:house, h.index, :link, t)]
        MOI.modify(m, cidx, MOI.ScalarCoefficientChange(pload[t], -1.0))
    end

    # ==========================================
    #    III. Add local constraints
    # ==========================================
    # (None here)

    # ==========================================
    #    IV. Add sub-resources to current model
    # ==========================================
    # (None here)

    return nothing
end


text:
"""
    LookingGlass.jl

This package contains a collection of code reflection and investigation utilities that can
be (potentially) useful for understanding, analyzing, and investigating julia code.

Most of the reflection functions are named hierarchically by what they reflect over. For
example:
- `func_specializations(f)` lists all specializations for a given function.
- `module_functions(m::Module)` lists all the functions in a given module.
- `module_globals_names(m::Module)` lists the names of all global variables in a given module.

There are also other non-reflection-focused utilities, such as `@quot`.
"""
module LookingGlass
import Base.Iterators

"""
    @quot 2 + 3

Macro for debugging during metaprogramming; it simply returns the exact expression that a
macro gets as input. This can be helpful if you're wondering what inputs a macro can expect.
"""
macro quot(e) QuoteNode(e) end

# ---------------------------------------------------------------------------
# -- Utilities for Functions
# ---------------------------------------------------------------------------

"""
    func_specializations(f) -> Dict((method,signature) => MethodInstance)

Return all the specializations (MethodInstances) of each method of the given function (or
callable object), if there are any.

Specializations (also called MethodInstances) are the compiler-generated code for a method,
generated by specializing the method for specific argument types.

Most often, it's just the keys that are interesting, if you want to know _what_ specializations
a function has, rather than accessing the MethodInstances themselves (they're
compiler-internal and sort of gross):

```julia-repl
julia> foo(x) = 2x
foo (generic function with 1 method)

julia> foo(2) + foo(2.0)
8.0

julia> keys(LookingGlass.func_specializations(foo))
Base.KeySet for a Dict{Tuple{Method,DataType},Core.TypeMapEntry} with 2 entries. Keys:
  (foo(x) in Main at none:1, Tuple{typeof(foo),Int64})
  (foo(x) in Main at none:1, Tuple{typeof(foo),Float64})
```
"""
func_specializations(f) =
    Dict((m, s.sig) => s
        for m in methods(f).ms
        # m.specializations returns a custom iterate, which we collect here via `s.next`
        for s in Channel() do ch
            s = m.specializations
            while s isa Core.TypeMapEntry  # TODO: Can there be other things in here besides TypeMapEntry?
                put!(ch, s)
                s=s.next
            end
        end)

"""
    func_backedges(f) -> Dict((method, signature) => [backedges...])

Return all the backedges on all specializations of a function, if there are any.
Backedges are a compiler-internal mechanism used for recursively invalidating the cached
compilations of a method -- its MethodInstances -- if the inlined code for any functions it
calls changes. This is the mechanism that fixes the infamous Julia #265 Issue.

Since these are _back_ edges, not _forward_ edges, if a function `foo` has a backedge listed
to function `bar`, this means `bar` has _inlined `foo`_ when it was compiled, so if the code
for `foo` is redefined by the user, `bar` needs to be recompiled the next time it's called.

```julia-repl
julia> foo(x) = 2x
foo (generic function with 1 method)

julia> bar(x) = foo(x) + 1
bar (generic function with 1 method)

julia> bar(1)
3

julia> LookingGlass.func_backedges(foo)
Dict{Any,Array{Any,1}} with 2 entries:
  (foo(x) in Main at none:1, Tuple{typeof(foo),Int64}) => Any[MethodInstance for bar(::Int64)]
  :MethodTable                                         => Any[]
```

Backedges on a function's MethodTable means changes to _any methods_ (including adding
a method) will trigger recompilation of the target function.
"""
function func_backedges(f)
    out = Dict{Tuple{Any,Type}, Vector{Any}}(
        # Method backedges
        k => try s.func.backedges catch ; [] end
        for (k,s) in func_specializations(f)
    )
    # MethodTable backedges
    for (typ, method) in
            # MethodTable edges are pairs of (type, method_instance) (but stored as a flat array)
            try Tuple.(Iterators.partition(methods(f).mt.backedges, 2)) catch ; [] end
        methods = get!(out, (:MethodTable, typ), [])
        push!(methods, method)
    end
    out
end

# ---------------------------------------------------------------------------
# -- Utilities for Modules
# ---------------------------------------------------------------------------

"""
    module_submodules(m::Module; recursive=true, base=false) -> Vector{Module}
Return a list of all Modules that are submodules of `m`. If `recursive=true`, this returns
all recursive submodules. If `recursive=false`, it returns only direct children of `m`.

By default, it skips Julia's `Base` and `Core` modules, but you can enable those with
`base=true`.
NOTE: currently, `module_submodules(Main, recursive=true, base=true)` will trigger a
StackOverflowError, since recursively listing submodules of `Core` or `Base` infinite loops.
"""
module_submodules(m::Module; recursive=true, base=false) =
    if recursive
        _module_recursive_submodules(m, base=base)
    else
        _module_direct_submodules(m, base=base)
    end

_module_recursive_submodules(m; base) = collect(Iterators.flatten(
    [x, _module_recursive_submodules(x, base=base)...]
        for x in _module_direct_submodules(m, base=base)
    ))
_module_direct_submodules(m; base) =
    Module[submodule
        for x in filter(x->name_is_submodule(m,x), names(m, all=true))
        for submodule in (Core.eval(m, x),)  # assign to temporary variable (comprehensions are weird)
        if submodule ∉ (Base, Core)]

name_is_submodule(m::Module, s::Symbol) =
    isdefined(m, s) && isa(Core.eval(m,s), Module) && nameof(m) != s
#module_is_submodule(m::Module, s::Module) = issubmodule(m, nameof(s))  # BROKEN: might just be a name collision. This is not the right way to check this.

"""
    module_functions_names(m::Module) -> Vector{Symbol}
Return a list of the names of all functions defined in Module `m`.
"""
module_functions_names(m::Module) =
    [name
     for name in names(m, all=true)
     if Core.eval(m, name) isa Function && name ∉ (:include, :eval)]
"""
    module_functions(m::Module) -> Dict{Symbol, Function}
Return a list of all the Function objects defined in Module `m`.
"""
module_functions(m::Module) = [Core.eval(m, n) for n in module_functions_names(m)]

"""
    module_globals_names(m; constness=:all, mutability=:all) -> Vector{Symbol}
Return a list of the names of all global variables in Module `m`.

To return only const globals or only non-const globals, set the `constness=` keyword
argument to one of `constness=:const` or `constness=:nonconst`.
To return only mutable globals or only immutable globals, set the `mutability=` keyword
argument to one of `mutability=:mutable` or `mutability=:immutable`.
"""
module_globals_names(m::Module; constness=:all, mutability=:all) =
    [n for n in names(m, all=true)
        if module_name_isglobal(m, n; constness=constness, mutability=mutability)]
function module_name_isglobal(m::Module, n::Symbol; constness, mutability)
    @assert constness ∈ (:all, :const, :nonconst)
    @assert mutability ∈ (:all, :mutable, :immutable)
    try
        if String(n)[1] == '#' return false end
        v = Core.eval(m, n)
        return !isa(v, Union{DataType, UnionAll, Function, Module}) &&
            (constness == :all || (constness == :const && _isconst_global(m, n) ||
                                   constness == :nonconst && !_isconst_global(m, n))) &&
            (mutability == :all || (mutability == :mutable && !isimmutable(v) ||
                                    mutability == :immutable && isimmutable(v)))
    catch
        false
    end
end
"""
    module_globals(m::Module) -> Dict{Symbol, Any}
Return a Dict mapping the name to the value of all global variables in Module `m`.
"""
module_globals(m::Module) = Dict(n => Core.eval(m, n) for n in module_globals_names(m))

function _isconst_global(m::Module, s::Symbol)
    b = _getbinding(m,s)
    b != nothing && b.constp
end
#--- Compiler-Internal Julia Binding struct ------------------------------------------------
# Stores information about a variable binding, which we can interrogate for information
# such as whether a variable is a global `const`.
# NOTE: This is very specific to a given version of Julia.
struct _Binding_t
    name::Ptr{Nothing}
    value::Ptr{Nothing}
    globalref::Ptr{Nothing}
    owner::Ptr{Nothing}
    constp::Bool
end
function _getbinding(m, s)
    p = ccall(:jl_get_binding, Ptr{_Binding_t}, (Any, Any), m, s)
    p == Ptr{Nothing}(0) ? nothing : unsafe_load(p)
end
#-------------------------------------------------------------------------------------------

"""
    module_recursive_globals_names(m; constness=:all, mutability=:all) -> Dict(submodule => name)
Return a list of the names of all global variables for each submodule in Module `m`.

Defaults to all globals, can toggle only const-globals or nonconst-globals via keyword args.
See [`module_recursive_globals`](@ref).
"""
module_recursive_globals_names(m::Module; constness=:all, mutability=:all) =
    merge!(
        Dict(m => module_globals_names(m, constness=constness, mutability=mutability)),
        Dict(
            sm => names
            for sm in module_submodules(m, recursive=true)
            for names in (module_globals_names(sm; constness=constness, mutability=mutability),)
            if !isempty(names)
        ))

end # module


text:
"""
	rates(param::parameters,iterations::Int64,divergence::Array,sfs::Array)

Function to solve randomly *N* scenarios. The function will create *N* models, defined by ```Analytical.parameters()```, to estimate analytically fixation and polymorphic rates for each model. The rates will be used to compute summary statistics required at ABC. The function output a HDF5 file containing the solved models, the selected DAC and the analytical rates.

If ρ and/or theta are set to ```nothing```, the function will input random values given the range 0.0005:0.0005:0.01. Otherwise you can fix the values.

If gL is set to ```nothing```, the function will not account the role of the weakly selected alleles in the estimation.

# Arguments
 - `param::parameters`: mutable structure containing the model
 - `binom::binomialDict` : structure containing the binomial convolution
 - `gH::Array{Int64,1}` : Range of strong selection coefficients
 - `gL::Union{Array{Int64,1},Nothing}`: Range of weak selection coefficients
 - `gam_neg::Array{Int64,1}` : Range of deleterious selection coefficients
 - `theta::Union{Float64,Nothing}` : Population-scaled mutation rate on coding region
 - `ρ::Union{Float64,Nothing}` : Population-scaled recombination rate
 - `shape::Float64=0.184` : DFE shape parameter
 - `iterations::Int64` : Number of solutions
 - `output::String` : File to output HDF5 file
# Returns
 - `Array`: summary statistics
 - `Output`: HDF5 file containing models solved and rates.
"""
function rates(;param::parameters,
				gH::S,
				gL::Union{S,Nothing},
				gam_neg::S,
				alpha::Vector{Float64}=[0.1,0.9],
				theta::Union{Float64,Nothing}=0.001,
				rho::Union{Float64,Nothing}=0.001,
				iterations::Int64,
				output::String,
				threads::Bool=false) where S <: Union{Vector{Int64},UnitRange{Int64}}

	# @assert (alpha[1] >= 0) & (alpha[end] <= 1) "α values must be at the interval [0,1]"

	assertion_params(param)

	# temp = mktempdir(homedir());
	# temp = temp .* "/tmp_" .* string.(1:iterations) .* ".txt"
	# Iterations = models to solve

	# Factor to modify input Γ(shape) parameter. Flexible Γ distribution over negative alleles
	fac     = rand(-2:0.05:2,iterations)
	afac    = @. param.al*(2^fac)

	# Deleting shape > 1. Negative alpha_x values
	idx = findall(afac .> 1)
	if !isempty(idx)
		afac[idx] = rand(afac[afac .< 1],size(idx,1))
	end

	# Random α values
	nTot    = rand(alpha[1]:0.01:alpha[end],iterations)

	# Defining αW. It is possible to solve non-accounting for weak fixations
	if isnothing(gL)
		# Setting αW to 0 for all estimations
		nLow    = fill(0.0,iterations)
		# Random strong selection coefficients
		ngl     = rand(repeat([1],iterations),iterations);
	else
		# Setting αW as proportion of α
		lfac    = rand(0.0:0.05:0.9,iterations)
		nLow    = @. nTot * lfac
		# Random weak selection coefficients
		ngl     = rand(repeat(gL,iterations),iterations);
	end

	# Random strong selection coefficients
	ngh     = rand(repeat(gH,iterations),iterations);
	# Random negative selection coefficients
	ngam_neg = rand(repeat(gam_neg,iterations),iterations);

	if !isnothing(theta)
		θ = fill(theta,iterations)
	else
		θ = rand(0.0005:0.0005:0.01,iterations)
	end

	# Random rho on coding regions
	if !isnothing(rho)
		ρ = fill(rho,iterations)
	else
		ρ = rand(0.0005:0.0005:0.05,iterations)
	end

	# Creating N models to iter in threads. Set N models (paramerters) and sampling probabilites (binomialDict)
	binom  = binom_op!(param)
	# Estimations to distributed workers
	if threads
		@time out = ThreadsX.mapi( (α_strong,α_weak,γ_strong,γ_weak,γ_neg,shape,θ_n,ρ_n) -> iter_rates(param, binom, α_strong,α_weak,γ_strong,γ_weak,γ_neg,shape,θ_n,ρ_n),nTot, nLow, ngh, ngl, ngam_neg, afac, θ, ρ);
	else
		@time m,r_ps,r_pn,r_f = unzip(ParallelUtilities.pmapbatch( (α_strong,α_weak,γ_strong,γ_weak,γ_neg,shape,θ_n,ρ_n) -> iter_rates(param, binom, α_strong,α_weak,γ_strong,γ_weak,γ_neg,shape,θ_n,ρ_n),nTot, nLow, ngh, ngl, ngam_neg, afac, θ, ρ));
	end

	# Reducing models array
	df  = vcat(m...);
	# Filtering cases where deleterious DFE cannot be solved.
	idx = sum.(eachrow(df)).!=0

	# Saving models and rates
	models = @view DataFrame(df,[:B,:al_low,:al_tot,:gam_neg,:gL,:gH,:al,:ρ])[idx,:];
	neut   = @view vcat(r_ps...)[idx,:];
	sel    = @view vcat(r_pn...)[idx,:];
	dsdn   = @view vcat(r_f...)[idx,:];

	# Saving multiple summary statistics
	n = OrderedDict{Int,Array}()
	s = OrderedDict{Int,Array}()
	for i ∈ eachindex(param.dac)
		n[param.dac[i]] = neut[:,i]
		s[param.dac[i]] = sel[:,i]
	end

	# Writting HDF5 file
	string_cutoff = "cutoff=[" * string(param.cutoff[1]) * "," * string(param.cutoff[end]) * "]"
	JLD2.jldopen(output, "a+") do file
		file[string(param.N)* "/" * string(param.n) * "/" * string_cutoff *  "/models"] = models;
		file[string(param.N)* "/" * string(param.n) * "/" * string_cutoff *  "/neut"]   = n;
		file[string(param.N)* "/" * string(param.n) * "/" * string_cutoff *  "/sel"]    = s;
		file[string(param.N)* "/" * string(param.n) * "/" * string_cutoff *  "/dsdn"]   = dsdn;
		file[string(param.N)* "/" * string(param.n) * "/" * string_cutoff *  "/dac"]    = param.dac;
	end;
end

"""
	iter_rates(param::parameters,afac::Float64,bfac::Float64,al_tot::Float64,al_low::Float64,divergence::Array,sfs::Array)

Estimating rates given a model for all B range.

# Arguments
 - `param::parameters`
 - `binom::Dict`
 - `al_tot::Float64`
 - `al_low::Float64`
 - `gH::Int64`
 - `gL::Int64`
 - `gam_neg::Int64`
 - `afac::Float64`
 - `ρ::Float64`
 - `θ::Float64`
# Output
 - `Array{Float64,2}`
"""
function iter_rates(param::parameters,binom::Dict{Float64, SparseMatrixCSC{Float64, Int64}},al_tot::Float64,al_low::Float64,gH::Int64,gL::Int64,gam_neg::Int64,afac::Float64,θ::Float64,ρ::Float64)

	# Creating model to solve
	# Γ distribution
	param.al    = afac; param.be = abs(afac/gam_neg); param.gam_neg = gam_neg
	# α, αW
	param.al_low = al_low; param.al_tot = al_tot;
	# Positive selection coefficients
	param.gH    = gH;param.gL = gL
	# Mutation rate and recomb
	param.θ_coding = θ; param.ρ = ρ
	#=param.θ_coding = θ; param.θᵣ .= θᵣ; param.ρ = ρ=#

	# Allocate array to solve the model for all B values
	m    = zeros(size(param.B_bins,1),8)
	r_ps = zeros(size(param.B_bins,1),size(param.dac,1))
	r_pn = zeros(size(param.B_bins,1),size(param.dac,1))
	r_f  = zeros(size(param.B_bins,1),4)
	# r  = zeros(size(param.B_bins,1),(size(param.dac,1) * 2) + 12)

	# Solving θ on non-coding region and probabilites to get α value without BGS
	param.B = 0.999
	set_θ!(param)
	try
		set_ppos!(param)
		for j ∈ eachindex(param.B_bins)
			# Set B value
			param.B = param.B_bins[j]
			# Solve θ non-coding for the B value.
			set_θ!(param)
			# Solve model for the B value
			x,y,z,w = try
				getting_rates(param,binom[param.B])
			catch e
				zeros(size(m,2)),zeros(size(r_ps,2)),zeros(size(r_pn,2)),zeros(size(r_pf,2))
			end
			@inbounds m[j,:] = x
			@inbounds r_ps[j,:] = y
			@inbounds r_pn[j,:] = z
			@inbounds r_f[j,:] = w
		end
		return (m,r_ps,r_pn,r_f)
	catch
		return (m,r_ps,r_pn,r_f)
	end

end

"""
	getting_rates(gammaL,gammaH,ppos_l,ppos_h,observedData,nopos)

Estimating analytical rates of fixation and polymorphism to approach α value accouting for background selection, weakly and strong positive selection. Output values will be used to sample from a Poisson distribution the total counts of polymorphism and divergence using observed data.

# Arguments
 - `param::parameters`
 - `binom::SparseMatrixCSC{Float64,Int64}`
# Returns
 - `Array{Float64,2}` containing solved model, fixation and polymorphic rates
"""
function getting_rates(param::parameters,binom::SparseMatrixCSC{Float64,Int64})

	################################################
	# Subset rates accounting for positive alleles #
	################################################

	# Fixation
	fN       = param.B*fix_neut(param)
	fNeg     = param.B*fix_neg(param,0.5*param.ppos_h+0.5*param.ppos_l)
	fPosL    = fix_pos_sim(param,param.gL,0.5*param.ppos_l)
	fPosH    = fix_pos_sim(param,param.gH,0.5*param.ppos_h)

	ds       = fN
	dn       = fNeg + fPosL + fPosH

	# Polymorphism
	neut::Array{Float64,1} = sfs_neut(param,binom)
	selH::Array{Float64,1} = if isinf(exp(param.gH * 2))
		sfs_pos_float(param,param.gH,param.ppos_h,binom)
	else
		sfs_pos(param,param.gH,param.ppos_h,binom)
	end
	selL::Array{Float64,1} = sfs_pos(param,param.gL,param.ppos_l,binom)
	selN::Array{Float64,1} = sfs_neg(param,param.ppos_h+param.ppos_l,binom)
	# Cumulative rates
	if(sum(param.cutoff)!=1)
		freq = round.(collect(1:(param.nn-1))./param.nn,digits=4)

		tmp = hcat(freq,neut,selH,selL,selN)
		tmp = @view tmp[(tmp[:,1] .>= param.cutoff[1]) .& (tmp[:,1] .<= param.cutoff[end]),2:end]

		tmp = cumulative_sfs(tmp,false)
	else
		tmp = cumulative_sfs(hcat(neut,selH,selL,selN),false)
	end

	split_columns(matrix::Array{Float64,2}) = (view(matrix, :, i) for i ∈ 1:size(matrix, 2));
	neut, selH, selL, selN = split_columns(tmp)
	sel = (selH+selL)+selN

	##########
	# Output #
	##########
	analytical_m::Matrix{Float64}  = vcat(param.B,param.al_low,param.al_tot,param.gam_neg,param.gL,param.gH,param.al,param.θ_coding)'
	analytical_ps::Matrix{Float64} = neut[param.dac]'
	analytical_pn::Matrix{Float64} = sel[param.dac]'
	analytical_f::Matrix{Float64}  = hcat(ds,dn,fPosL,fPosH)
	# analytical_values::Array{Float64,2} = vcat(param.B,param.al_low,param.al_tot,param.gam_neg,param.gL,param.gH,param.al,param.θ_coding,neut[param.dac],sel[param.dac],ds,dn,fPosL,fPosH)'

	return (analytical_m,analytical_ps,analytical_pn,analytical_f)
end


text:
println(".... run examples and notebooks ....")

using NBInclude

const ex_dir = "../examples/src/"
const nb_dir = "../examples/nb/"

dir = pwd()
cd(@__DIR__)

ex = [x for x in readdir(ex_dir) if occursin(r".jl", x)]
nb = [x for x in readdir(nb_dir) if occursin(r".ipynb", x)]

for x in ex
    println("... include $ex_dir$x:")
    include(ex_dir*x)
end

for x in nb
    println("... @nbinclude $nb_dir$x:")
    @nbinclude(nb_dir*x)
end

cd(dir)


text:
module Observations

export SyntheticObservations, observation_times

using ..Utils: prettyvector

using Oceananigans
using Oceananigans: fields
using Oceananigans.Grids: AbstractGrid
using Oceananigans.Grids: cpu_face_constructor_x, cpu_face_constructor_y, cpu_face_constructor_z
using Oceananigans.Grids: pop_flat_elements, topology, halo_size, on_architecture
using Oceananigans.TimeSteppers: update_state!, reset!
using Oceananigans.Fields
using Oceananigans.Utils: SpecifiedTimes, prettytime
using Oceananigans.Architectures
using Oceananigans.Architectures: arch_array, architecture
using JLD2

import Oceananigans.Fields: set!

using ParameterEstimocean.Utils: field_name_pairs
using ParameterEstimocean.Transformations: Transformation, compute_transformation

abstract type AbstractObservation end

struct SyntheticObservations{F, R, G, T, P, M, Þ} <: AbstractObservation
    field_time_serieses :: F
    forward_map_names :: R
    grid :: G
    times :: T
    path :: P
    metadata :: M
    transformation :: Þ
end

"""
    SyntheticObservations(path=nothing;
                          field_names,
                          forward_map_names = field_names,
                          transformation = Transformation(),
                          times = nothing,
                          field_time_serieses = nothing,
                          architecture = CPU(),
                          regrid = nothing)

Return a time series of synthetic observations generated by Oceananigans.jl's simulations
gridded as Oceananigans.jl fields.
"""
function SyntheticObservations(path=nothing;
                               field_names,
                               forward_map_names = field_names,
                               transformation = Transformation(),
                               times = nothing,
                               field_time_serieses = nothing,
                               architecture = CPU(),
                               regrid = nothing)

    field_names = tupleit(field_names)
    forward_map_names = tupleit(forward_map_names)

    all(n ∈ field_names for n in forward_map_names) ||
        throw(ArgumentError("All of the forward map names $forward_map_names must be in field names $field_names"))

    if isnothing(field_time_serieses)
        raw_time_serieses = NamedTuple(name => FieldTimeSeries(path, string(name); times, architecture) for name in field_names)
    else
        raw_time_serieses = field_time_serieses
    end

    raw_grid = first(raw_time_serieses).grid
    times = first(raw_time_serieses).times
    boundary_conditions = first(raw_time_serieses).boundary_conditions

    if isnothing(regrid)
        field_time_serieses = raw_time_serieses
        grid = raw_grid

    else # Well, we're gonna regrid stuff
        if regrid isa Tuple
            grid = with_size(regrid, raw_grid)
        elseif regrid isa AbstractGrid
            grid = regrid
        else
            error("$regrid must be a Tuple specifying a new grid size, or a new grid!")
        end

        @info string("Regridding synthetic observations...", '\n',
                     "    original grid: ", summary(raw_grid), '\n',
                     "         new grid: ", summary(grid))

        field_time_serieses = Dict()

        # Re-grid the data in `field_time_serieses`
        for (field_name, ts) in zip(keys(raw_time_serieses), raw_time_serieses)

            #LX, LY, LZ = location(ts[1])
            LX, LY, LZ = infer_location(field_name)

            new_ts = FieldTimeSeries{LX, LY, LZ}(grid, times; boundary_conditions)

            # Loop over time steps to re-grid each constituent field in `field_time_series`
            for n = 1:length(times)
                regrid!(new_ts[n], ts[n])
            end

            field_time_serieses[field_name] = new_ts
        end

        field_time_serieses = NamedTuple(field_time_serieses)
    end

    # validate_data(fields, grid, times) # might be a good idea to validate the data...
    if !isnothing(path)
        file = jldopen(path)
        metadata = NamedTuple(Symbol(group) => read_group(file[group])
                              for group in filter(n -> n ∉ not_metadata_names, keys(file)))
        close(file)
    else
        metadata = nothing
    end

    transformation = field_name_pairs(transformation, forward_map_names, "transformation")
    transformation = Dict(name => compute_transformation(transformation[name], field_time_serieses[name])
                          for name in forward_map_names)

    return SyntheticObservations(field_time_serieses, forward_map_names,
                                 grid, times, path, metadata, transformation)
end

observation_names(observations::SyntheticObservations) = keys(observations.field_time_serieses)
forward_map_names(observations::SyntheticObservations) = observations.forward_map_names

function get_names(observations, name_getter)
    names = Set()
    for obs in observations
        push!(names, name_getter(obs)...)
    end

    return names
end

"""
    observation_names(observations::Vector{<:SyntheticObservations})

Return a Set representing the union of all names in `observations`.
"""
observation_names(observations::Vector{<:SyntheticObservations}) = get_names(observations, observation_names)

"""
    forward_map_names(observations::Vector{<:SyntheticObservations})

Return a Set representing the union of all names in `observations`.
"""
forward_map_names(observations::Vector{<:SyntheticObservations}) = get_names(observations, forward_map_names)

Base.summary(observations::SyntheticObservations) =
    "SyntheticObservations of $(keys(observations.field_time_serieses)) on $(summary(observations.grid))"

Base.summary(observations::Vector{<:SyntheticObservations}) =
    "Vector{<:SyntheticObservations} of $(keys(first(observations).field_time_serieses)) on $(summary(first(observations).grid))"

tupleit(t) = try
    Tuple(t)
catch
    tuple(t)
end

const not_metadata_names = ("serialized", "timeseries")

read_group(group::JLD2.Group) = NamedTuple(Symbol(subgroup) => read_group(group[subgroup]) for subgroup in keys(group))
read_group(group) = group

using Oceananigans.Grids: ZRegRectilinearGrid

function with_size(new_size, old_grid)

    old_grid isa ZRegRectilinearGrid ||
        error("Cannot remake stretched grid \n $grid \n with a new size!")

    topo = topology(old_grid)

    x = cpu_face_constructor_x(old_grid)
    y = cpu_face_constructor_y(old_grid)
    z = cpu_face_constructor_z(old_grid)

    # Remove elements of size and new_halo in Flat directions as expected by grid
    # constructor
    new_size = pop_flat_elements(new_size, topo)
    halo = pop_flat_elements(halo_size(old_grid), topo)

    new_grid = RectilinearGrid(architecture(old_grid), eltype(old_grid);
        size = new_size,
        x = x, y = y, z = z,
        topology = topo,
        halo = halo)

    return new_grid
end

location_guide = Dict(:u => (Face, Center, Center),
                      :v => (Center, Face, Center),
                      :w => (Center, Center, Face))

function infer_location(field_name)
    if field_name in keys(location_guide)
        return location_guide[field_name]
    else
        return (Center, Center, Center)
    end
end

function observation_times(data_path::String)
    file = jldopen(data_path)
    iterations = parse.(Int, keys(file["timeseries/t"]))
    times = [file["timeseries/t/$i"] for i in iterations]
    close(file)
    return times
end

observation_times(observation::SyntheticObservations) = observation.times

function observation_times(obs::Vector)
    @assert all([o.times ≈ obs[1].times for o in obs]) "Observations must have the same times."
    return observation_times(first(obs))
end

#####
##### set! for simulation models and observations
#####

"""
    column_ensemble_interior(observations::Vector{<:SyntheticObservations},
                             field_name, time_index, (Nensemble, Nbatch, Nz))

Return an `Nensemble × Nbatch × Nz` Array of `(1, 1, Nz)` `field_name` data,
given `Nbatch` `SyntheticObservations` objects. The `Nbatch × Nz` data for `field_name`
is copied `Nensemble` times to form a 3D Array.
"""
function column_ensemble_interior(observations::Vector{<:SyntheticObservations},
                                  field_name, time_index, (Nensemble, Nbatch, Nz))

    zeros_column = zeros(1, 1, Nz)
    Nt = length(first(observations).times)

    batched_data = []
    for observation in observations
        fts = observation.field_time_serieses
        if field_name in keys(fts) && time_index <= Nt
            field_column = interior(fts[field_name][time_index])
            push!(batched_data, interior(fts[field_name][time_index]))
        else
            push!(batched_data, zeros_column)
        end
    end

    # Make a Vector of 1D Array into a 3D Array
    flattened_data = cat(batched_data..., dims = 2) # (Nbatch, Nz)
    ensemble_interior = cat((flattened_data for i = 1:Nensemble)..., dims = 1) # (Nensemble, Nbatch, Nz)

    return ensemble_interior
end

function set!(model, obs::SyntheticObservations, time_index=1)
    for field_name in keys(fields(model))
        model_field = fields(model)[field_name]

        if field_name ∈ keys(obs.field_time_serieses)
            obs_field = obs.field_time_serieses[field_name][time_index]
            set!(model_field, obs_field)
        else
            fill!(parent(model_field), 0)
        end
    end

    update_state!(model)

    return nothing
end

function set!(model, observations::Vector{<:SyntheticObservations}, time_index=1)
    for field_name in keys(fields(model))
        model_field = fields(model)[field_name]
        model_field_size = size(model_field)
        Nensemble = model.grid.Nx

        observations_data = column_ensemble_interior(observations, field_name, time_index, model_field_size)

        # Reshape `observations_data` to the size of `model_field`'s interior
        reshaped_data = arch_array(architecture(model_field), reshape(observations_data, size(model_field)))

        # Sets the interior of field `model_field` to values of `reshaped_data`
        model_field .= reshaped_data
    end

    update_state!(model)

    return nothing
end

#####
##### FieldTimeSeriesCollector for collecting data while a simulation runs
#####

struct FieldTimeSeriesCollector{G, D, F}
    grid :: G
    times :: Vector{Float64}
    field_time_serieses :: D
    collected_fields :: F
    collection_times :: Vector{Float64}
end

"""
    FieldTimeSeriesCollector(collected_fields, times;
                             architecture = Architectures.architecture(first(collected_fields)))

Return a `FieldTimeSeriesCollector` for `fields` of `simulation`.
`fields` is a `NamedTuple` of `AbstractField`s that are to be collected.
"""
function FieldTimeSeriesCollector(collected_fields, times;
                                  architecture = Architectures.architecture(first(collected_fields)))

    grid = on_architecture(architecture, first(collected_fields).grid)
    field_time_serieses = Dict{Symbol, Any}()

    for field_name in keys(collected_fields)
        field = collected_fields[field_name]
        LX, LY, LZ = location(field)
        field_time_series = FieldTimeSeries{LX, LY, LZ}(grid, times)
        field_time_serieses[field_name] = field_time_series
    end

    # Convert to NamedTuple
    field_time_serieses = NamedTuple(name => field_time_serieses[name] for name in keys(collected_fields))

    collection_times = similar(times)

    return FieldTimeSeriesCollector(grid, times, field_time_serieses, collected_fields, collection_times)
end

function (collector::FieldTimeSeriesCollector)(simulation)
    for field in collector.collected_fields
        compute!(field)
    end

    current_time = simulation.model.clock.time
    time_index = findfirst(t -> t ≈ current_time, collector.times)
    if isnothing(time_index)
        @warn string("Current time ", prettytime(current_time), " not found in
                     time collector times ", prettytime.(collector.times))
        return nothing
    end

    for field_name in keys(collector.collected_fields)
        field_time_series = collector.field_time_serieses[field_name]
        if architecture(collector.grid) != architecture(simulation.model.grid)
            arch = architecture(collector.grid)
            device_collected_field_data = arch_array(arch, parent(collector.collected_fields[field_name]))
            parent(field_time_series[time_index]) .= device_collected_field_data
        else
            set!(field_time_series[time_index], collector.collected_fields[field_name])
        end
    end

    # We _have_ collected data
    collector.collection_times[time_index] = current_time

    return nothing
end

#####
##### Initializing simulations
#####

nothingfunction(simulation) = nothing

function initialize_forward_run!(simulation, observations,
                                 time_series_collector, initialize_simulation! = nothingfunction)

    reset!(simulation)

    times = observation_times(observations)
    initial_time = times[1]
    simulation.model.clock.time = initial_time

    collected_fields = time_series_collector.collected_fields
    arch = architecture(time_series_collector.grid)

    # Clear potential NaNs from timestepper data.
    # Particularly important for Adams-Bashforth timestepping scheme.
    # Oceananigans ≤ v0.71 initializes the Adams-Bashforth scheme with an Euler step by *multiplying* the tendency
    # at time-step n-1 by 0. Because 0 * NaN = NaN, this fails when the tendency at n-1 contains NaNs.
    timestepper = simulation.model.timestepper
    for field in tuple(timestepper.Gⁿ..., timestepper.G⁻...)
        if !isnothing(field)
            parent(field) .= 0
        end
    end

    # Initialize FieldTimeSeriesCollector
    time_series_collector.collection_times .= 0
    for time_series in time_series_collector.field_time_serieses
        parent(time_series) .= 0
    end

    :nan_checker ∈ keys(simulation.callbacks) && pop!(simulation.callbacks, :nan_checker)
    :data_collector ∈ keys(simulation.callbacks) && pop!(simulation.callbacks, :data_collector)
    simulation.callbacks[:data_collector] = Callback(time_series_collector, SpecifiedTimes(times...))
    simulation.stop_time = times[end]

    set!(simulation.model, observations, 1)

    initialize_simulation!(simulation)

    return nothing
end

summarize_metadata(::Nothing) = ""
summarize_metadata(metadata) = keys(metadata)

function Base.show(io::IO, obs::SyntheticObservations)
    times_str = prettyvector(prettytime.(obs.times, false))

    print(io, "SyntheticObservations with fields $(propertynames(obs.field_time_serieses))", '\n',
              "├── times: $times_str", '\n',
              "├── grid: $(summary(obs.grid))", '\n',
              "├── path: \"$(obs.path)\"", '\n',
              "├── metadata: ", summarize_metadata(obs.metadata), '\n',
              "└── transformation: $(summary(obs.transformation))")
end

end # module


text:
using Base: min, max

export
    Rectangle,
    set!, intersect!, intersects,
    bounds!, contains_point

# A two-dimensional axis-aligned rectangle.
#
# X-axis directed towards the right
#
# Y-Axis directed downward.
#
# left(X)/Top(Y)
#       *------------.  --> X
#       |            |  |
#       |            |  v Y
#       |            |
#       |            |
#       .------------*
#              Right(X)/Bottom(Y)
#
mutable struct Rectangle{T <: AbstractFloat}
    # Top-left corner
    min::Point{T}
    # Bottom-right corner
    max::Point{T}

    width::T
    height::T

    function Rectangle{T}() where {T <: AbstractFloat}
        new(Point{T}(), Point{T}(1.0, 1.0), 1.0, 1.0)
    end

    function Rectangle{T}(minx::T, miny::T, maxx::T, maxy::T) where {T <: AbstractFloat}
        new(Point{T}(minx, miny), Point{T}(maxx, maxy), maxx - minx, maxy - miny)
    end

    function Rectangle{T}(min::Point{T}, max::Point{T}) where {T <: AbstractFloat}
        new(min, max, max.x - min.x, max.y - min.y)
    end

end

# setters/getters
function set!(rect::Rectangle{T}, minx::T, miny::T, maxx::T, maxy::T) where {T <: AbstractFloat}
    set!(rect.min, minx, miny);
    set!(rect.max, maxx, maxy);
    rect.width = maxx - minx;
    rect.height = maxy - miny;
end

# algorithms
function intersect!(intersect::Rectangle{T}, rectA::Rectangle{T}, rectB::Rectangle{T}) where {T <: AbstractFloat}
    x0 = max(rectA.min.x, rectB.min.x);
    x1 = min(rectA.max.x, rectB.max.x);

    if x0 <= x1
        y0 = max(rectA.min.y, rectB.min.y);
        y1 = min(rectA.max.y, rectB.max.y);

        if y0 <= y1
            set!(intersect, x0, y0, x1, y1);
        end
    end
end

function intersects(rectA::Rectangle{T}, rectB::Rectangle{T}) where {T <: AbstractFloat}
    rectA.min.x <= rectB.min.x + rectB.width &&
    rectB.min.x <= rectA.min.x + rectA.width &&
    rectA.max.y <= rectB.max.y + rectB.height &&
    rectB.max.y <= rectA.max.y + rectA.height
end

# Returns a new rectangle which completely contains `rectA` and `rectB`.
function bounds!(bounds::Rectangle{T}, rectA::Rectangle{T}, rectB::Rectangle{T}) where {T <: AbstractFloat}
    right = max(rectA.max.x, rectB.max.x);
    bottom = max(rectA.max.y, rectB.max.y);
    left = min(rectA.min.x, rectB.min.x);
    top = min(rectA.min.y, rectB.min.y);

    set!(bounds.min, left, top);
    set!(bounds.max, right, bottom);
    bounds.width = right - left;
    bounds.height = bottom - top;
end

function contains_point(rect::Rectangle{T}, p::Point{T}) where {T <: AbstractFloat}
    p.x >= rect.min.x &&
  p.x <= rect.max.x &&
  p.y >= rect.min.y &&
  p.y <= rect.max.y
end

max_stars_repo_name: wdevore/Ranger-Julia-SDL

max_issues_repo_name: wdevore/Ranger-Julia-SDL

max_forks_repo_name: wdevore/Ranger-Julia-SDL


text:
struct ModelConfiguration{T<:AbstractFloat}
    κh::T  # Horizontal Laplacian heat diffusion [m²/s]. diffKhT in MITgcm.
    κv::T  # Vertical Laplacian heat diffusion [m²/s]. diffKzT in MITgcm.
    𝜈h::T  # Horizontal eddy viscosity [m²/s]. viscAh in MITgcm.
    𝜈v::T  # Vertical eddy viscosity [m²/s]. viscAz in MITgcm.

    function ModelConfiguration(κh, κv, 𝜈h, 𝜈v)
        @assert κh >= 0
        @assert κv >= 0
        @assert 𝜈h >= 0
        @assert 𝜈v >= 0
        new{Float64}(κh, κv, 𝜈h, 𝜈v)
    end
end


text:
# include
include("Include.jl")

# load data and model files -
path_model_file = "Coagulation.net"
path_data_file = joinpath(_PATH_TO_DATA, "Thrombin-TF.csv")

# load the experimental data -
data_table = load(path_data_file)
model_buffer = read_model_file(path_model_file)

# get size -
(ℛ, 𝒞) = size(data_table)
SF = 1e9 # all concentrations are in nmol/L

# build the default dictionary -
model = build_default_model_dictionary(model_buffer)

# setup -
FIIa_archive = Array{Array{Float64,1},1}(undef, ℛ)
T_archive = Array{Array{Float64,1},1}(undef, ℛ)

# main loop -
for i ∈ 1:ℛ

    # build new model -
    dd = deepcopy(model)

    # setup static -
    sfa = dd["static_factors_array"]
    sfa[1] = (2.5e-9) * SF        # 1 TFPI
    sfa[2] = (3.4e-6) * SF        # 2 AT
    sfa[3] = (5e-12) * SF         # 3 TF
    sfa[6] = 0.005                # 6 TRAUMA

    # grab the multiplier from the data -
    ℳ = dd["number_of_dynamic_states"]
    xₒ = zeros(ℳ)
    xₒ[1] = (1.4e-6) * (data_table[i, :II] / 100.0)     # 1 FII
    xₒ[2] = (1e-8) * (data_table[i, :VII] / 100.0)      # 2 FVII
    xₒ[3] = (2e-8) * (data_table[i, :V] / 100.0)        # 3 FV
    xₒ[4] = (1.6e-7) * (data_table[i, :X] / 100.0)      # 4 FX
    xₒ[5] = (7e-10) * (data_table[i, :VIII] / 100.0)    # 5 FVIII
    xₒ[6] = (9e-8) * (data_table[i, :IX] / 100.0)       # 6 FIX
    xₒ[7] = (1e-8) * (data_table[i, :XI] / 100.0)       # 7 FXI
    xₒ[8] = (1e-8) * (data_table[i, :XII] / 100.0)      # 8 FXII
    xₒ[9] = (1e-13)  # 9 FIIa
    xₒ = SF * xₒ # convert to nmol

    # update α -
    α = dd["α"]
    α[1] = 0.061
    α[9] = 0.70

    # setup -
    G = dd["G"]
    idx = indexin(dd, "FVIIa")
    G[idx, 4] = 0.1

    # what is the index of TRAUMA?
    idx = indexin(dd, "AT")
    G[idx, 9] = 0.045

    # what is the index of TFPI?
    idx = indexin(dd, "TFPI")
    G[idx, 1] = -0.65

    # Phase 3: solve the model -
    # setup the solver -
    tspan = (0.0, 20.0)
    prob = ODEProblem(balances, xₒ, tspan, dd; saveat = 0.01)
    soln = solve(prob)

    T = soln.t
    X = hcat(soln.u...)

    FIIa_archive[i] = X[9, :]
    T_archive[i] = T
end

# put simulation archives -
X = hcat(FIIa_archive...)
T = hcat(T_archive...)

# compute the stylized fact -
lagtime_array = τ_lag(x -> x > 2.0, T, X)
peak_FIIa = max_FIIa(X)
pta = peak_time(T, X)

max_stars_repo_name: varnerlab/UVM-TopDown-LegacyModel

max_issues_repo_name: varnerlab/UVM-TopDown-LegacyModel

max_forks_repo_name: varnerlab/UVM-TopDown-LegacyModel


text:
"""
    bar(x, y)

Creates an `EChart` where values plotted vertically as rectangular columns.

## Methods
```julia
bar(x::AbstractVector, y::AbstractVector{<:Union{Missing, Real}})
bar(x::AbstractVector, y::AbstractArray{<:Union{Missing, Real},2})
bar(df::AbstractDataFrame, x::Symbol, y::Symbol)
bar(df::AbstractDataFrame, x::Symbol, y::Symbol, group::Symbol)
```

## Arguments
* `mark::Union{String, AbstractVector} = "bar"` : how to display plotted points
* `stack::Union{Bool, AbstractVector, Void} = nothing` : stack (add together) when multple series present?
* `legend::Bool` : display legend?
* `scale::Bool = false` : show full Y-axis or truncated
* `kwargs` : varargs to set any field of resulting `EChart` struct

## Notes

Reasonable defaults set for different methods of `bar`, such as displaying a legend when two or more series present.

## Examples
```julia
x = ["Monday","Tuesday","Wednesday","Thursday","Friday","Saturday","Sunday"]
y = [11, 11, 15, 13, 12, 13, 10]
b = bar(x, y)
```
"""
function bar(x::AbstractVector, y::AbstractVector{<:Union{Missing, Real}};
			mark::Union{String, AbstractVector} = "bar",
			stack::Union{Bool, AbstractVector, Nothing} = nothing,
			legend::Bool = false,
			scale::Bool = false,
			kwargs...)

	 return xy_plot(x, y; mark = mark, stack = stack, legend = legend, scale = scale, kwargs...)

end

#multiple y
function bar(x::AbstractVector, y::AbstractArray{<:Union{Missing, Real},2};
			mark::Union{String, AbstractVector} = "bar",
			stack::Union{Bool, AbstractVector, Nothing} = nothing,
			legend::Bool = true,
			scale::Bool = false,
			kwargs...)

	 return xy_plot(x, y; mark = mark, stack = stack, legend = legend, scale = scale, kwargs...)

end

# dataframe, single series
function bar(df::AbstractDataFrame, x::Symbol, y::Symbol;
			mark::Union{String, AbstractVector} = "bar",
			stack::Union{Bool, AbstractVector, Nothing} = nothing,
			legend::Bool = false,
			scale::Bool = false,
			kwargs...)

	 return xy_plot(df, x, y; mark = mark, stack = stack, legend = legend, scale = scale, kwargs...)

end

# dataframe, group argument
function bar(df::AbstractDataFrame, x::Symbol, y::Symbol, group::Symbol;
			mark::Union{String, AbstractVector} = "bar",
			stack::Union{Bool, AbstractVector, Nothing} = nothing,
			legend::Bool = true,
			scale::Bool = false,
			kwargs...)

	 return xy_plot(df, x, y, group; mark = mark, stack = stack, legend = legend, scale = scale, kwargs...)

end