text:
using CImGui
using CImGui: ImVec2, FLT_MAX, ImGuiSizeCallbackData
using CImGui.CSyntax
using CImGui.CSyntax.CStatic

"""
    ShowExampleAppConstrainedResize(p_open::Ref{Bool})
Create a window with custom resize constraints.
"""
global function ShowExampleAppConstrainedResize(p_open::Ref{Bool})
    function Square(data::Ptr{ImGuiSizeCallbackData})::Cvoid
        desired_size = CImGui.Get_DesiredSize(data)
        max_size = max(desired_size.x, desired_size.y)
        CImGui.Set_DesiredSize(data, ImVec2(max_size, max_size))
    end

    function Step(data::Ptr{ImGuiSizeCallbackData})::Cvoid
        desired_size = CImGui.Get_DesiredSize(data)
        step::Cfloat = CImGui.Get_UserData(data) |> Int
        size_x = trunc(desired_size.x / step + 0.5) * step
        size_y = trunc(desired_size.y / step + 0.5) * step
        CImGui.Set_DesiredSize(data, ImVec2(size_x, size_y))
    end

    square_fptr = @cfunction($Square, Cvoid, (Ptr{ImGuiSizeCallbackData},))
    step_fptr = @cfunction($Step, Cvoid, (Ptr{ImGuiSizeCallbackData},))

    @cstatic auto_resize=false type=Cint(0) display_lines=Cint(10) begin
        type == 0 && CImGui.SetNextWindowSizeConstraints(ImVec2(-1, 0),    ImVec2(-1, FLT_MAX))      # vertical only
        type == 1 && CImGui.SetNextWindowSizeConstraints(ImVec2(0, -1),    ImVec2(FLT_MAX, -1))      # horizontal only
        type == 2 && CImGui.SetNextWindowSizeConstraints(ImVec2(100, 100), ImVec2(FLT_MAX, FLT_MAX)) # width > 100, Height > 100
        type == 3 && CImGui.SetNextWindowSizeConstraints(ImVec2(400, -1),  ImVec2(500, -1))          # width 400-500
        type == 4 && CImGui.SetNextWindowSizeConstraints(ImVec2(-1, 400),  ImVec2(-1, 500))          # height 400-500
        type == 5 && CImGui.SetNextWindowSizeConstraints(ImVec2(0, 0),     ImVec2(FLT_MAX, FLT_MAX), square_fptr) # always Square
        type == 6 && CImGui.SetNextWindowSizeConstraints(ImVec2(0, 0),     ImVec2(FLT_MAX, FLT_MAX), step_fptr, Ptr{Cvoid}(100)) # fixed Step

        flags = auto_resize ? CImGui.ImGuiWindowFlags_AlwaysAutoResize : 0
        if CImGui.Begin("Example: Constrained Resize", p_open, flags)
            desc = [
                "Resize vertical only",
                "Resize horizontal only",
                "Width > 100, Height > 100",
                "Width 400-500",
                "Height 400-500",
                "Custom: Always Square",
                "Custom: Fixed Steps (100)",
            ]
            CImGui.Button("200x200") && CImGui.SetWindowSize(ImVec2(200, 200))
            CImGui.SameLine()
            CImGui.Button("500x500") && CImGui.SetWindowSize(ImVec2(500, 500))
            CImGui.SameLine()
            CImGui.Button("800x200") && CImGui.SetWindowSize(ImVec2(800, 200))
            CImGui.PushItemWidth(200)
            @c CImGui.Combo("Constraint", &type, desc, length(desc))
            @c CImGui.DragInt("Lines", &display_lines, 0.2, 1, 100)
            CImGui.PopItemWidth()
            @c CImGui.Checkbox("Auto-resize", &auto_resize)
            foreach(i->CImGui.Text("$(" "^4i)Hello, sailor! Making this line long enough for the example."), 1:display_lines)
        end
        CImGui.End()
    end # @cstatic
end


text:
# Effective perceived noise level
function f_epnl(t_o, pnlt)

    # Number of time steps
    n_t = size(t_o)[1]

    # Interpolate time, pnlt and C
    dt = 0.5
    n_ip = Int64(ceil((t_o[end]-t_o[1])/dt))

    T = eltype(pnlt)
    t_ip = zeros(T, n_ip)
    for i in range(1, n_ip, step=1)
        t_ip[i] = t_o[1] + (i-1)*dt
    end

    # Interpolate the data
    #f_ipnlterp = PCHIPInterpolation.Interpolator(t_o, reshape(pnlt, (n_t, )))
    f_ipnlterp = LinearInterpolation(t_o, reshape(pnlt, (n_t, )))
    pnlt_ip = f_ipnlterp.(t_ip)

    # Compute max. PNLT
    pnltm = maximum(pnlt_ip)

    # Compute max. PNLT point (k_m)
    I = findall(x->x>pnltm - 10.,pnlt_ip)

    # ICAO Annex 16 procedures (p132 - App. 2-18)
    f_int = 10 .^(pnlt_ip / 10.)

    # Compute duration correction
    if pnltm > 10
        # Compute integration bounds: lower bound
        i_1 = I[1]
        if i_1 != 1
            if abs(pnlt_ip[i_1] - (pnltm - 10)) > abs(pnlt_ip[i_1 - 1] - (pnltm - 10))
                i_1 = i_1 - 1
            end
        end

        # Compute integration bounds: upper bound
        i_2 = I[end]
        if i_2 != size(pnlt_ip)[1]
            if (abs(pnlt_ip[i_2] - (pnltm - 10))) > (abs(pnlt_ip[i_2 + 1] - (pnltm - 10)))
                i_2 = i_2 + 1
            end
        end

        # Duration correction
        D = 10 * log10(sum(f_int[i_1:i_2])) - pnltm - 10 * log10(10. / dt)
    else
        # Duration correction
        D = 10 * log10(sum(f_int)) - pnltm - 10 * log10(10. / dt)
    end

    # Compute EPNL
    epnl = pnltm + D

    return epnl
end

ext: jl

max_stars_repo_name: MIT-LAE/pyNA

max_issues_repo_name: MIT-LAE/pyNA

max_forks_repo_name: MIT-LAE/pyNA


text:
module FeldtLib

export setgooglemapsapikey!, swedish_postnr2position, address2position
export ggplot2_distribution

#include("ranks.jl")
#include("baumgartner_weis_schindler.jl")
#include("fitall.jl")
include("cacheddict.jl")
include("geocoding.jl")
include("geonames.jl")
include("ggplot2_graphs.jl")

end


text:
using Documenter, SBML

makedocs(
    modules = [SBML],
    clean = false,
    format = Documenter.HTML(
        prettyurls = !("local" in ARGS),
        canonical = "https://lcsb-biocore.github.io/DistributedData.jl/stable/",
    ),
    sitename = "SBML.jl",
    authors = "The developers of SBML.jl",
    linkcheck = !("skiplinks" in ARGS),
    pages = ["Documentation" => "index.md"],
)

deploydocs(
    repo = "github.com/LCSB-BioCore/SBML.jl.git",
    target = "build",
    branch = "gh-pages",
    push_preview = true,
)


text:
@doc raw"""
# Massless charged particle in 2D

The Lagrangian is given by
```math
L(x, \dot{x}) = A(x) \cdot \dot{x} - \phi (x) ,
```
with magnetic vector potential
```math
A(x) = \frac{A_0}{2} \big( 1 + x_1^2 + x_2^2 \big) \begin{pmatrix}
- x_2 \\
+ x_1 \\
\end{pmatrix} ,
```
electrostatic potential
```math
\phi(x) =  E_0 \, \big( \cos (x_1) + \sin(x_2) \big) ,
```
and magnetic and electric fields
```math
\begin{aligned}
B(x) &= \nabla \times A(x) = A_0 \, (1 + 2 x_1^2 + 2 x_2^2) , \\
E(x) &= - \nabla \phi(x) = E_0 \, \big( \sin x_1, \, - \cos x_2 \big)^T .
\end{aligned}
```

The Hamiltonian form of the equations of motion reads
```math
\dot{x} = \frac{1}{B(x)} \begin{pmatrix}
\hphantom{-} 0 & + 1 \\
- 1 & \hphantom{+} 0 \\
\end{pmatrix} \nabla \phi (x) .
```

"""
module MasslessChargedParticle

    using GeometricEquations

    using ..Diagnostics

    import ..Diagnostics: compute_invariant_error, compute_momentum_error

    export ϑ, A, B, ϕ, E, hamiltonian
    export massless_charged_particle_ode, massless_charged_particle_iode,
           massless_charged_particle_idae, massless_charged_particle_idae_spark
    export compute_energy_error, compute_momentum_error

    # default simulation parameters
    Δt = 0.2
    nt = 5000

    # default initial conditions and parameters
    q₀ = [1.0, 1.0]
    parameters = (A₀ = 1.0, E₀ = 1.0)

    # components of the vector potential
    A₁(q, params) = - params[:A₀] * q[2] * (1 + q[1]^2 + q[2]^2) / 2
    A₂(q, params) = + params[:A₀] * q[1] * (1 + q[1]^2 + q[2]^2) / 2

    A(q, params) = [A₁(q, params), A₂(q, params)]

    # z-componend of the magnetic field
    B(q, params) = params[:A₀] * (1 + 2 * q[1]^2 + 2 * q[2]^2)

    # electrostatic potential
    ϕ(q, params) = params[:E₀] * (cos(q[1]) + sin(q[2]))
    # ϕ(q, params) = E₀ * (q[1]^2 + q[2]^2)

    # components of the electric field
    E₁(q, params) = + params[:E₀] * sin(q[1])
    E₂(q, params) = - params[:E₀] * cos(q[2])
    # E₁(q, params) = + params[:E₀] * q[1]
    # E₂(q, params) = + params[:E₀] * q[2]

    E(q, params) = [E₁(q, params), E₂(q, params)]

    # components of the velocity
    v₁(t, q, params) = - E₂(q, params) / B(q, params)
    v₂(t, q, params) = + E₁(q, params) / B(q, params)

    # components of the one-form (symplectic potential)
    ϑ₁(t, q, params) = A₁(q, params)
    ϑ₂(t, q, params) = A₂(q, params)

    ϑ(t, q, params) = [ϑ₁(t, q, params), ϑ₂(t, q, params)]

    function ϑ(t::Number, q::AbstractVector, params::NamedTuple, k::Int)
        if k == 1
            ϑ₁(t, q, params)
        elseif k == 2
            ϑ₂(t, q, params)
        else
            throw(BoundsError(ϑ,k))
        end
    end

    # derivatives of the one-form components
    dϑ₁dx₁(t, q, params) = - params[:A₀] * q[1] * q[2]
    dϑ₁dx₂(t, q, params) = - params[:A₀] * (1 + q[1]^2 + 3 * q[2]^2) / 2
    dϑ₂dx₁(t, q, params) = + params[:A₀] * (1 + 3 * q[1]^2 + q[2]^2) / 2
    dϑ₂dx₂(t, q, params) = + params[:A₀] * q[1] * q[2]

    # components of the force
    f₁(t, q, v, params) = dϑ₁dx₁(t, q, params) * v[1] + dϑ₂dx₁(t, q, params) * v[2]
    f₂(t, q, v, params) = dϑ₁dx₂(t, q, params) * v[1] + dϑ₂dx₂(t, q, params) * v[2]

    g₁(t, q, v, params) = dϑ₁dx₁(t, q, params) * v[1] + dϑ₁dx₂(t, q, params) * v[2]
    g₂(t, q, v, params) = dϑ₂dx₁(t, q, params) * v[1] + dϑ₂dx₂(t, q, params) * v[2]

    # Hamiltonian (total energy)
    hamiltonian(t, q, params) = ϕ(q, params)

    # components of the gradient of the Hamiltonian
    dHd₁(t, q, params) = - E₁(q, params)
    dHd₂(t, q, params) = - E₂(q, params)

    function massless_charged_particle_dH(t, q, dH, params)
        dH[1] = dHd₁(t, q, params)
        dH[2] = dHd₂(t, q, params)
        nothing
    end

    function massless_charged_particle_v(t, q, v, params)
        v[1] = v₁(t, q, params)
        v[2] = v₂(t, q, params)
        nothing
    end

    function massless_charged_particle_ϑ(t, q, Θ, params)
        Θ[1] = ϑ₁(t, q, params)
        Θ[2] = ϑ₂(t, q, params)
        nothing
    end

    massless_charged_particle_ϑ(t, q, v, Θ, params) = massless_charged_particle_ϑ(t, q, Θ, params)

    function massless_charged_particle_f(t, q, v, f, params)
        f[1] = f₁(t, q, v, params) - dHd₁(t, q, params)
        f[2] = f₂(t, q, v, params) - dHd₂(t, q, params)
        nothing
    end

    function massless_charged_particle_f̄(t, q, v, f, params)
        f[1] = - dHd₁(t, q, params)
        f[2] = - dHd₂(t, q, params)
        nothing
    end

    function massless_charged_particle_g(t, q, v, g, params)
        g[1] = f₁(t, q, v, params)
        g[2] = f₂(t, q, v, params)
        nothing
    end

    massless_charged_particle_g(t, q, p, v, g, params) = massless_charged_particle_g(t, q, v, g, params)

    function massless_charged_particle_u(t, q, v, u, params)
        u .= v
        nothing
    end

    massless_charged_particle_u(t, q, p, v, u, params) = massless_charged_particle_u(t, q, v, u, params)

    function massless_charged_particle_ϕ(t, q, p, ϕ, params)
        ϕ[1] = p[1] - ϑ₁(t,q,params)
        ϕ[2] = p[2] - ϑ₂(t,q,params)
        nothing
    end

    function massless_charged_particle_ψ(t, q, p, v, f, ψ, params)
        ψ[1] = f[1] - g₁(t,q,v,params)
        ψ[2] = f[2] - g₂(t,q,v,params)
        nothing
    end

    "Creates an ODE object for the massless charged particle in 2D."
    function massless_charged_particle_ode(q₀=q₀; params=parameters)
        ODE(massless_charged_particle_v, q₀; invariants=(h=hamiltonian,), parameters=params)
    end

    "Creates an implicit ODE object for the massless charged particle in 2D."
    function massless_charged_particle_iode(q₀=q₀; params=parameters)
        IODE(massless_charged_particle_ϑ, massless_charged_particle_f,
                massless_charged_particle_g, q₀, ϑ(0., q₀, params);
                v̄=massless_charged_particle_v, f̄=massless_charged_particle_f,
                invariants=(h=hamiltonian,), parameters=params)
    end

    "Creates an implicit DAE object for the massless charged particle in 2D."
    function massless_charged_particle_idae(q₀=q₀; params=parameters)
        IDAE(massless_charged_particle_ϑ, massless_charged_particle_f,
                massless_charged_particle_u, massless_charged_particle_g,
                massless_charged_particle_ϕ, q₀, ϑ(0., q₀, params), zero(q₀);
                v̄=massless_charged_particle_v, f̄=massless_charged_particle_f,
                invariants=(h=hamiltonian,), parameters=params)
    end

    "Creates an implicit DAE object for the massless charged particle in 2D."
    function massless_charged_particle_idae_spark(q₀=q₀; params=parameters)
        IDAE(massless_charged_particle_ϑ, massless_charged_particle_f̄,
                massless_charged_particle_u, massless_charged_particle_g,
                massless_charged_particle_ϕ, q₀, ϑ(0., q₀, params), zero(q₀);
                v̄=massless_charged_particle_v, f̄=massless_charged_particle_f,
                invariants=(h=hamiltonian,), parameters=params)
    end

    compute_energy_error(t,q,params) = compute_invariant_error(t,q, (t,q) -> hamiltonian(t,q,params))
    compute_momentum_error(t,q,p,params::NamedTuple) = compute_momentum_error(t, q, p, (t,q,k) -> ϑ(t,q,params,k))

end


text:
# test memory allocation

using MicrochipRTDM
using FTD2XX

projectdirectory = "c:/Users/Chris/MPLABXProjects/lbcmcp/PWM_RTDM.x"
mapdict = microchip_parsemap(projectdirectory,"production")
uartconfig1 = UARTConfiguration(921600,8,1,"n",readtimeout = 1000,writetimeout = 1000)
io = open(FT_DeviceIndex(0),uartconfig1)
setlatencytimer(io,2); #Should speed things up a little
rtdminterface = RTDMInterface(mapdict,io)

const address = UInt32(mapdict[:offsetP3])
const buffer = [0x1234]
function f()
  for i in 1:1000
    #rtdm_write(rtdminterface, buffer, address)
    rtdm_write(rtdminterface, 0x1234, address)
  end
  return nothing
end
f()

Profile.clear_malloc_data()
f()

close(rtdminterface)
quit()


text:
function accumulate!(
    ocn :: Ocean,
)

    for f in fieldnames(AccumulativeVariables)
        v = getfield(ocn.acc_vars, f)
        v .+= getfield(ocn, f)
    end

end

function avg_accumulate!(
    ocn :: Ocean;
    count :: Real,
)
    for f in fieldnames(AccumulativeVariables)

        ocn_field = getfield(ocn, f)
        acc_field = getfield(ocn.acc_vars, f)

        for i in eachindex(ocn_field)
            ocn_field[i] = acc_field[i] / count
        end

        # clear variable by default
        acc_field .= 0.0

    end

end

max_stars_repo_name: meteorologytoday/IOM

max_issues_repo_name: meteorologytoday/IOM

max_forks_repo_name: meteorologytoday/IOM


text:
using FractionalDelayFilter
using Test

@testset "Max filter order" begin
    delaystep = 0.25
    delays = 0.0:delaystep:5.0
    numdelays = length(delays)
    fdfcoefs = Vector{Vector{Float64}}(undef, numdelays)
    nonzeroindices = Vector{Int}(undef, numdelays)
    filterorders = Vector{Int}(undef, numdelays)
    sinewave = sin.(0:π/40:4π)

    maxorder = 3
    @testset "Odd value: $maxorder" begin
        for delayindex = 1:numdelays
            filterorders[delayindex] = filtord(delays[delayindex], maxorder=maxorder)
            (fdfcoefs[delayindex], nonzeroindices[delayindex]) = getfdfcoef(filterorders[delayindex], delays[delayindex])
        end

        delayedwave = zeros(Float64, length(sinewave)*2, numdelays)
        for delayindex = 1:numdelays
            delayedwave[1+nonzeroindices[delayindex]:length(sinewave)+nonzeroindices[delayindex]+length(fdfcoefs[delayindex])-2,delayindex] = fdfilter(sinewave, fdfcoefs[delayindex])
        end
        @test sortperm(delayedwave[7,:], rev=true) == collect(1:numdelays)
    end

    maxorder = 4
    @testset "Even value: $maxorder" begin
        for delayindex = 1:numdelays
            filterorders[delayindex] = filtord(delays[delayindex], maxorder=maxorder)
            (fdfcoefs[delayindex], nonzeroindices[delayindex]) = getfdfcoef(filterorders[delayindex], delays[delayindex])
        end

        delayedwave = zeros(Float64, length(sinewave)*2, numdelays)
        for delayindex = 1:numdelays
            delayedwave[1+nonzeroindices[delayindex]:length(sinewave)+nonzeroindices[delayindex]+length(fdfcoefs[delayindex])-2,delayindex] = fdfilter(sinewave, fdfcoefs[delayindex])
        end
        @test sortperm(delayedwave[7,:], rev=true) == collect(1:numdelays)
    end

end


text:
using Db
using Test

@testset "Db.jl" begin
    # Write your tests here.
end


text:
include("../src/BlueMesh7.jl")
using .BlueMesh7

using Random
using Plots
using Statistics: mean
using ProgressMeter: @showprogress
using BSON: @save, @load

ngraphs = 4
graphs = [generate_positions(dims=(30, 30), n=96) for _ = 1:ngraphs]
# ■

emit_rates = collect(10:10:50) |> reverse
setup_nrelays = collect(0:16:96) |> reverse

# ■
pdrs = []

@showprogress for emit_rate in emit_rates
    for nrelays in setup_nrelays
        localpdrs = []

        Threads.@threads for graph in graphs
            node_roles = zeros(length(graph))

            for idx = 1:nrelays
                node_roles[idx] = 1
            end

            mesh = initialize_mesh(graph, Random.shuffle(node_roles))
            mesh.packet_emit_rate = emit_rate / 1000

            push!(localpdrs, start(mesh, minutes=1).pdr)
        end

        push!(pdrs, mean(localpdrs))
    end
end

@save "surface.bson" pdrs
# ■

p = plot(reverse(setup_nrelays), reverse(emit_rates), pdrs, st = :surface, camera=(60, 40),
         xticks=([reverse(setup_nrelays);], setup_nrelays), yticks=([reverse(emit_rates);], emit_rates),
         xlabel="relays", ylabel="emit", zlabel="pdr")

savefig("resurfaced-more.png")
println("wrote resurfaced-more.png")


text:
using CSV, Distributions, HypothesisTests

data1 = CSV.read("../data/machine1.csv", header=false)[:,1]
data2 = CSV.read("../data/machine2.csv", header=false)[:,1]
xBar1, n1 = mean(data1), length(data1)
xBar2, n2 = mean(data2), length(data2)
sig1, sig2 = 1.2, 1.6
delta0 = 0

testStatistic = ( xBar1-xBar2 - delta0 ) / ( sqrt( sig1^2 / n1 + sig2^2 / n2 ) )
pVal = 2*ccdf(Normal(), abs(testStatistic))

println("Sample mean machine 1: ",xBar1)
println("Sample mean machine 2: ",xBar2)
println("Manually calculated test statistc: ", testStatistic)
println("Manually calculated p-value: ", pVal)

max_stars_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook

max_issues_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook

max_forks_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook


text:
Base.show(io::IO, coll::Collection) = print(io, "Collection($(coll.database), \"", coll.name, "\")")

# docstring for this function is at database.jl
function command_simple(collection::Collection, command::BSON) :: BSON
    reply = BSON()
    err_ref = Ref{BSONError}()
    ok = mongoc_collection_command_simple(collection.handle, command.handle, C_NULL, reply.handle,
        err_ref)

    if !ok
        throw(err_ref[])
    end
    return reply
end

function read_command(collection::Collection, command::BSON;
        options::Union{Nothing, BSON}=nothing) :: BSON

    reply = BSON()
    err_ref = Ref{BSONError}()
    options_handle = options == nothing ? C_NULL : options.handle
    ok = mongoc_collection_read_command_with_opts(collection.handle,
            command.handle, C_NULL, options_handle, reply.handle, err_ref)

    if !ok
        throw(err_ref[])
    end
    return reply
end

# Aux function to add _id field to document if it does not exist.
# If `_id_` is not present, creates a copy of `document`.
# This function returns `document, inserted_oid`.
function _new_id(document::BSON)
    if haskey(document, "_id")
        return document, nothing
    else
        inserted_oid = BSONObjectId()

        # copies it so this function doesn't have side effects
        document = deepcopy(document)

        document["_id"] = inserted_oid
        return document, inserted_oid
    end
end

function insert_one(collection::Collection, document::BSON;
        options::Union{Nothing, BSON}=nothing) :: InsertOneResult

    document, inserted_oid = _new_id(document)
    reply = BSON()
    err_ref = Ref{BSONError}()
    options_handle = options == nothing ? C_NULL : options.handle
    ok = mongoc_collection_insert_one(collection.handle, document.handle, options_handle,
                                      reply.handle, err_ref)
    if !ok
        throw(err_ref[])
    end
    return InsertOneResult(reply, inserted_oid)
end

function delete_one(collection::Collection, selector::BSON; options::Union{Nothing, BSON}=nothing)
    reply = BSON()
    err_ref = Ref{BSONError}()
    options_handle = options == nothing ? C_NULL : options.handle
    ok = mongoc_collection_delete_one(collection.handle, selector.handle, options_handle,
                                      reply.handle, err_ref)
    if !ok
        throw(err_ref[])
    end
    return reply
end

function delete_many(collection::Collection, selector::BSON; options::Union{Nothing, BSON}=nothing)
    reply = BSON()
    err_ref = Ref{BSONError}()
    options_handle = options == nothing ? C_NULL : options.handle
    ok = mongoc_collection_delete_many(collection.handle, selector.handle, options_handle,
                                       reply.handle, err_ref)
    if !ok
        throw(err_ref[])
    end
    return reply
end

function update_one(collection::Collection, selector::BSON, update::BSON;
        options::Union{Nothing, BSON}=nothing)

    reply = BSON()
    err_ref = Ref{BSONError}()
    options_handle = options == nothing ? C_NULL : options.handle
    ok = mongoc_collection_update_one(collection.handle, selector.handle, update.handle,
                                      options_handle, reply.handle, err_ref)
    if !ok
        throw(err_ref[])
    end
    return reply
end

function update_many(collection::Collection, selector::BSON, update::BSON;
        options::Union{Nothing, BSON}=nothing)

    reply = BSON()
    err_ref = Ref{BSONError}()
    options_handle = options == nothing ? C_NULL : options.handle

    ok = mongoc_collection_update_many(collection.handle, selector.handle, update.handle,
        options_handle, reply.handle, err_ref)

    if !ok
        throw(err_ref[])
    end
    return reply
end

"""
    replace_one(
            collection::Collection,
            selector::BSON,
            replacement::BSON;
            options::Union{Nothing, BSON}=nothing
        )

Replace the result of querying by `selector` with `replacement`.
"""
function replace_one(
            collection::Collection,
            selector::BSON,
            replacement::BSON;
            options::Union{Nothing, BSON}=nothing
        )

    reply = BSON()
    err_ref = Ref{BSONError}()
    options_handle = options == nothing ? C_NULL : options.handle
    ok = mongoc_collection_replace_one(collection.handle, selector.handle, replacement.handle,
                                       options_handle, reply.handle, err_ref)
    if !ok
        throw(err_ref[])
    end

    return reply
end

function BulkOperationResult(reply::BSON, server_id::UInt32)
    BulkOperationResult(reply, server_id, Vector{Union{Nothing, BSONObjectId}}())
end

function execute!(bulk_operation::BulkOperation) :: BulkOperationResult
    if bulk_operation.executed
        error("Bulk operation was already executed.")
    end

    try
        reply = BSON()
        err_ref = Ref{BSONError}()

        bulk_operation_result = mongoc_bulk_operation_execute(bulk_operation.handle,
            reply.handle, err_ref)

        if bulk_operation_result == 0
            throw(err_ref[])
        end
        return BulkOperationResult(reply, bulk_operation_result)
    finally
        destroy!(bulk_operation)
    end
end

function bulk_insert!(bulk_operation::BulkOperation, document::BSON;
        options::Union{Nothing, BSON}=nothing)

    err_ref = Ref{BSONError}()
    options_handle = options == nothing ? C_NULL : options.handle
    ok = mongoc_bulk_operation_insert_with_opts(bulk_operation.handle, document.handle,
        options_handle, err_ref)

    if !ok
        throw(err_ref[])
    end
    nothing
end

function bulk_replace_one!(bulk_operation::BulkOperation, selector::BSON, replacement::BSON;
        options::Union{Nothing, BSON}=nothing)

    err_ref = Ref{BSONError}()
    options_handle = options === nothing ? C_NULL : options.handle
    ok = mongoc_bulk_operation_replace_one_with_opts(
        bulk_operation.handle,
        selector.handle,
        replacement.handle,
        options_handle,
        err_ref
    )
    if !ok
        throw(err_ref[])
    end
    nothing
end

function insert_many(collection::Collection, documents::Vector{BSON};
        bulk_options::Union{Nothing, BSON}=nothing, insert_options::Union{Nothing, BSON}=nothing)

    inserted_oids = Vector{Union{Nothing, BSONObjectId}}()

    bulk_operation = BulkOperation(collection, options=bulk_options)
    for doc in documents
        doc, inserted_oid = _new_id(doc)
        bulk_insert!(bulk_operation, doc, options=insert_options)
        push!(inserted_oids, inserted_oid)
    end
    result = execute!(bulk_operation)
    append!(result.inserted_oids, inserted_oids)
    return result
end

"""
    find(collection::Collection, bson_filter::BSON=BSON();
        options::Union{Nothing, BSON}=nothing) :: Cursor

Executes a query on `collection` and returns an iterable `Cursor`.

# Example

```julia
function find_contract_codes(collection, criteria::Dict=Dict()) :: Vector{String}
    result = Vector{String}()

    let
        bson_filter = Mongoc.BSON(criteria)
        bson_options = Mongoc.BSON(\"\"\"{ "projection" : { "_id" : true }, "sort" : { "_id" : 1 } }\"\"\")
        for bson_document in Mongoc.find(collection, bson_filter, options=bson_options)
            push!(result, bson_document["_id"])
        end
    end

    return result
end
```

Check the [libmongoc documentation](http://mongoc.org/libmongoc/current/mongoc_collection_find_with_opts.html)
for more information.
"""
function find(collection::Collection, bson_filter::BSON=BSON();
        options::Union{Nothing, BSON}=nothing) :: Cursor

    options_handle = options == nothing ? C_NULL : options.handle
    cursor_handle = mongoc_collection_find_with_opts(collection.handle, bson_filter.handle,
        options_handle, C_NULL)

    if cursor_handle == C_NULL
        error("Couldn't execute query.")
    end
    return Cursor(collection, cursor_handle)
end

"""
    count_documents(collection::Collection, bson_filter::BSON=BSON();
        options::Union{Nothing, BSON}=nothing) :: Int

Returns the number of documents on a `collection`,
with an optional filter given by `bson_filter`.

`length(collection)` and `Mongoc.count_documents(collection)`
produces the same output.

# Example

```julia
result = length(collection, Mongoc.BSON("_id" => oid))
```
"""
function count_documents(collection::Collection, bson_filter::BSON=BSON();
        options::Union{Nothing, BSON}=nothing) :: Int

    err_ref = Ref{BSONError}()
    options_handle = options == nothing ? C_NULL : options.handle
    len = mongoc_collection_count_documents(collection.handle, bson_filter.handle,
        options_handle, C_NULL, C_NULL, err_ref)

    if len == -1
        throw(err_ref[])
    end
    return Int(len)
end

function set_limit!(cursor::Cursor, limit::Int)
    ok = mongoc_cursor_set_limit(cursor.handle, limit)
    if !ok
        error("Couldn't set cursor limit to $limit.")
    end
    nothing
end

"""
    find_one(collection::Collection, bson_filter::BSON=BSON();
        options::Union{Nothing, BSON}=nothing) :: Union{Nothing, BSON}

Execute a query to a collection and returns the first element of the result set.

Returns `nothing` if the result set is empty.
"""
function find_one(collection::Collection, bson_filter::BSON=BSON();
                  options::Union{Nothing, BSON}=nothing) :: Union{Nothing, BSON}

    cursor = find(collection, bson_filter, options=options)
    set_limit!(cursor, 1)
    next = _iterate(cursor)
    if next == nothing
        return nothing
    else
        bson_document, _state = next
        return bson_document
    end
end

"""
    aggregate(collection::Collection, bson_pipeline::BSON;
        flags::QueryFlags=QUERY_FLAG_NONE,
        options::Union{Nothing, BSON}=nothing) :: Cursor

Use `Mongoc.aggregate` to execute an aggregation command.

# Example

The following reproduces the example from the
[MongoDB Tutorial](https://docs.mongodb.com/manual/aggregation/).

```julia
docs = [
    Mongoc.BSON(\"\"\"{ "cust_id" : "A123", "amount" : 500, "status" : "A" }\"\"\"),
    Mongoc.BSON(\"\"\"{ "cust_id" : "A123", "amount" : 250, "status" : "A" }\"\"\"),
    Mongoc.BSON(\"\"\"{ "cust_id" : "B212", "amount" : 200, "status" : "A" }\"\"\"),
    Mongoc.BSON(\"\"\"{ "cust_id" : "A123", "amount" : 300, "status" : "D" }\"\"\")
]

collection = client["my-database"]["aggregation-collection"]
append!(collection, docs)

# Sets the pipeline command
bson_pipeline = Mongoc.BSON(\"\"\"
    [
        { "\$match" : { "status" : "A" } },
        { "\$group" : { "_id" : "\$cust_id", "total" : { "\$sum" : "\$amount" } } }
    ]
\"\"\")

for doc in Mongoc.aggregate(collection, bson_pipeline)
  println(doc)
end
```

The result of the script above is:

```julia
BSON("{ "_id" : "B212", "total" : 200 }")
BSON("{ "_id" : "A123", "total" : 750 }")
```
"""
function aggregate(collection::Collection, bson_pipeline::BSON;
                   flags::QueryFlags=QUERY_FLAG_NONE,
                   options::Union{Nothing, BSON}=nothing) :: Cursor

    options_handle = options == nothing ? C_NULL : options.handle
    cursor_handle = mongoc_collection_aggregate(collection.handle, flags,
                        bson_pipeline.handle, options_handle, C_NULL)

    if cursor_handle == C_NULL
        error("Couldn't execute aggregate command.")
    end
    return Cursor(collection, cursor_handle)
end

"""
    drop(database::Database, opts::Union{Nothing, BSON}=nothing)
    drop(collection::Collection, opts::Union{Nothing, BSON}=nothing)

Drops `database` or `collection`.

For information about `opts` argument, check the libmongoc documentation for
[database drop](http://mongoc.org/libmongoc/current/mongoc_database_drop_with_opts.html)
or
[collection drop](http://mongoc.org/libmongoc/current/mongoc_collection_drop_with_opts.html).
"""
function drop(collection::Collection, opts::Union{Nothing, BSON}=nothing)
    opts_handle = opts == nothing ? C_NULL : opts.handle
    err_ref = Ref{BSONError}()
    ok = mongoc_collection_drop_with_opts(collection.handle, opts_handle, err_ref)

    if !ok
        throw(err_ref[])
    end
    nothing
end

# findAndModify
function Base.setproperty!(builder::FindAndModifyOptsBuilder, opt::Symbol, val::BSON)
    if opt == :update
        set_opt_update!(builder, val)

    elseif opt == :sort
        set_opt_sort!(builder, val)

    elseif opt == :fields
        set_opt_fields!(builder, val)

    else
        error("Unknown option for FindAndModifyOptsBuilder: $opt")
    end
end

function Base.setproperty!(builder::FindAndModifyOptsBuilder, opt::Symbol, val::FindAndModifyFlags)
    @assert opt == :flags "Can't set $val to field $opt."
    set_opt_flags!(builder, val)
end

function set_opt_update!(builder::FindAndModifyOptsBuilder, val::BSON)
    ok = mongoc_find_and_modify_opts_set_update(builder.handle, val.handle)
    if !ok
        error("Couldn't set option update $val for FindAndModifyOptsBuilder.")
    end
    nothing
end

function set_opt_sort!(builder::FindAndModifyOptsBuilder, val::BSON)
    ok = mongoc_find_and_modify_opts_set_sort(builder.handle, val.handle)
    if !ok
        error("Couldn't set option sort $val for FindAndModifyOptsBuilder.")
    end
    nothing
end

function set_opt_fields!(builder::FindAndModifyOptsBuilder, val::BSON)
    ok = mongoc_find_and_modify_opts_set_fields(builder.handle, val.handle)
    if !ok
        error("Couldn't set option fields $val for FindAndModifyOptsBuilder.")
    end
    nothing
end

function set_opt_flags!(builder::FindAndModifyOptsBuilder, val::FindAndModifyFlags)
    ok = mongoc_find_and_modify_opts_set_flags(builder.handle, val)
    if !ok
        error("Couldn't set option flags $val for FindAndModifyOptsBuilder.")
    end
    nothing
end

function Base.setproperty!(builder::FindAndModifyOptsBuilder, opt::Symbol, val::Bool)
    @assert opt == :bypass_document_validation "Can't set $val to field $opt."
    set_opt_bypass_document_validation!(builder, val)
end

function set_opt_bypass_document_validation!(builder::FindAndModifyOptsBuilder, bypass::Bool)
    ok = mongoc_find_and_modify_opts_set_bypass_document_validation(builder.handle, bypass)
    if !ok
        error("Couldn't set option `bypass document validation` for FindAndModifyOptsBuilder.")
    end
    nothing
end

"""
    find_and_modify(collection::Collection, query::BSON;
        update::Union{Nothing, BSON}=nothing,
        sort::Union{Nothing, BSON}=nothing,
        fields::Union{Nothing, BSON}=nothing,
        flags::Union{Nothing, FindAndModifyFlags}=nothing,
        bypass_document_validation::Bool=false,
    ) :: BSON

Find documents and updates them in one go.

See [`Mongoc.FindAndModifyFlags`](@ref) for a list of accepted values
for `flags` argument.

# C API

* [mongoc_collection_find_and_modify](http://mongoc.org/libmongoc/current/mongoc_collection_find_and_modify.html).
"""
function find_and_modify(collection::Collection, query::BSON;
            update::Union{Nothing, BSON}=nothing,
            sort::Union{Nothing, BSON}=nothing,
            fields::Union{Nothing, BSON}=nothing,
            flags::Union{Nothing, FindAndModifyFlags}=nothing,
            bypass_document_validation::Bool=false,
        ) :: BSON

    opts = FindAndModifyOptsBuilder(
            update=update,
            sort=sort,
            fields=fields,
            flags=flags,
            bypass_document_validation=bypass_document_validation,
    )

    reply = BSON()
    err_ref = Ref{BSONError}()

    ok = mongoc_collection_find_and_modify_with_opts(
            collection.handle,
            query.handle, opts.handle,
            reply.handle, err_ref)

    if !ok
        throw(err_ref[])
    end

    return reply
end

const FIND_ONE_AND_DELETE_FIELDS = Dict(
    "collation" => "collation",
    "maxTimeMS" => "maxTimeMS",
    "projection" => "fields",
    "sort" => "sort",
    "writeConcern" => "writeConcern",
)

"""
    find_one_and_delete(collection::Collection, bson_filter::BSON;
        options::Union{Nothing, BSON}=nothing
    ) :: Union{Nothing, BSON}

Delete a document and return it.

See [db.collection.findOneAndDelete](https://docs.mongodb.com/manual/reference/method/db.collection.findOneAndDelete/)
for a list of accepted options.
"""
function find_one_and_delete(collection::Collection, bson_filter::BSON;
        options::Union{Nothing, BSON}=nothing) :: Union{Nothing, BSON}

    command_bson = BSON(
        "findAndModify" => collection.name,
        "query" => bson_filter,
        "remove" => true,
    )

    options_bson = BSON()

    if !isnothing(options)
        for (opt_key, comm_key) in FIND_ONE_AND_DELETE_FIELDS
            if haskey(options, opt_key)
                command_bson[comm_key] = options[opt_key]
            end
        end

        if haskey(options, "sessionId")
            options_bson["sessionId"] = options["sessionId"]
        end
    end

    reply = write_command(
        collection.database,
        command_bson,
        options=options_bson,
    )

    return reply["value"]
end

const FIND_ONE_AND_MODIFY_FIELDS = Dict(
    "arrayFilters" => "arrayFilters",
    "bypassDocumentValidation" .=> "bypassDocumentValidation",
    "collation" => "collation",
    "maxTimeMS" => "maxTimeMS",
    "projection" => "fields",
    "returnNewDocument" => "new",
    "sort" => "sort",
    "upsert" => "upsert",
    "writeConcern" => "writeConcern",
)

function _find_one_and_modify(collection::Collection, bson_filter::BSON, bson_modify::BSON;
        options::Union{Nothing, BSON}=nothing) :: Union{Nothing, BSON}

    command_bson = BSON(
        "findAndModify" => collection.name,
        "query" => bson_filter,
        "update" => bson_modify,
    )

    options_bson = BSON()

    if !isnothing(options)
        for (opt_key, comm_key) in FIND_ONE_AND_MODIFY_FIELDS
            if haskey(options, opt_key)
                command_bson[comm_key] = options[opt_key]
            end
        end

        if haskey(options, "sessionId")
            options_bson["sessionId"] = options["sessionId"]
        end
    end

    reply = write_command(
        collection.database,
        command_bson,
        options=options_bson,
    )

    return reply["value"]
end

"""
    find_one_and_replace(collection::Collection, bson_filter::BSON, bson_replacement::BSON;
        options::Union{Nothing, BSON}=nothing
    ) :: Union{Nothing, BSON}

Replace a document and return the original.

See [db.collection.findOneAndReplace](https://docs.mongodb.com/manual/reference/method/db.collection.findOneAndReplace/)
for a list of accepted options.
"""
function find_one_and_replace(collection::Collection, bson_filter::BSON, bson_replacement::BSON;
        options::Union{Nothing, BSON}=nothing) :: Union{Nothing, BSON}

    # Ensure the document has no operators
    for key in keys(bson_replacement)
        if startswith(key, '$')
            throw(ArgumentError("replacement must not contain update operators"))
        end
    end

    return _find_one_and_modify(collection, bson_filter, bson_replacement, options=options)
end

"""
    find_one_and_update(collection::Collection, bson_filter::BSON, bson_update::BSON;
        options::Union{Nothing, BSON}=nothing
    ) :: Union{Nothing, BSON}

Update a document and return the original.

See [db.collection.findOneAndUpdate](https://docs.mongodb.com/manual/reference/method/db.collection.findOneAndUpdate/)
for a list of accepted options.
"""
function find_one_and_update(collection::Collection, bson_filter::BSON, bson_update::BSON;
        options::Union{Nothing, BSON}=nothing) :: Union{Nothing, BSON}

    # Ensure the document only contains operators
    for key in keys(bson_update)
        if !startswith(key, '$')
            throw(ArgumentError("update must only contain update operators"))
        end
    end

    return _find_one_and_modify(collection, bson_filter, bson_update, options=options)
end

#
# High-level API
#

function _iterate(cursor::Cursor, state::Nothing=nothing)
    bson_handle_ref = Ref{Ptr{Cvoid}}()
    has_next = mongoc_cursor_next(cursor.handle, bson_handle_ref)

    if has_next
        # The bson document is valid only until the next call to mongoc_cursor_next.
        # So we should return a deepcopy.
        return deepcopy(BSON(bson_handle_ref[], enable_finalizer=false)), nothing
    else
        err_ref = Ref{BSONError}()
        if mongoc_cursor_error(cursor.handle, err_ref)
            throw(err_ref[])
        end

        return nothing
    end
end

Base.iterate(cursor::Cursor, state::Nothing=nothing) = _iterate(cursor, state)

function Base.iterate(coll::Collection)
    cursor = find(coll)
    return iterate(coll, cursor)
end

function Base.iterate(coll::Collection, state::Cursor)
    next = _iterate(state)
    if next == nothing
        return nothing
    else
        doc, _ = next
        return doc, state
    end
end

function Base.push!(collection::AbstractCollection, document::BSON;
                    options::Union{Nothing, BSON}=nothing)
    insert_one(collection, document; options=options)
end

function Base.append!(collection::AbstractCollection, documents::Vector{BSON};
                      bulk_options::Union{Nothing, BSON}=nothing,
                      insert_options::Union{Nothing, BSON}=nothing)
    insert_many(collection, documents; bulk_options=bulk_options, insert_options=insert_options)
end

function Base.length(collection::AbstractCollection, bson_filter::BSON=BSON();
                     options::Union{Nothing, BSON}=nothing)
    count_documents(collection, bson_filter; options=options)
end

function Base.isempty(collection::AbstractCollection, bson_filter::BSON=BSON();
                      options::Union{Nothing, BSON}=nothing)
    count_documents(collection, bson_filter; options=options) == 0
end

Base.empty!(collection::C) where {C<:AbstractCollection} = delete_many(collection, BSON())

function Base.collect(cursor::Cursor) :: Vector{BSON}
    result = Vector{BSON}()
    for doc in cursor
        push!(result, doc)
    end
    return result
end

function Base.collect(collection::AbstractCollection, bson_filter::BSON=BSON())
    return collect(find(collection, bson_filter))
end


text:
"""
    struct SkeletonPair{L,R} <: GridapType
      left::L
      right::R
    end
"""
struct SkeletonPair{L,R} <: GridapType
  left::L
  right::R
end

function reindex(a::AbstractArray,b::SkeletonPair)
  left = reindex(a,b.left)
  right = reindex(a,b.right)
  SkeletonPair(left,right)
end


text:
# MadNLP.jl.
# Created by Sungho Shin (sungho.shin@wisc.edu)

# Build info
default_linear_solver() = MadNLPUmfpack
default_dense_solver() = MadNLPLapackCPU
default_iterator() = MadNLPRichardson

# Dummy module
module DummyModule end

# Logger
@kwdef mutable struct Logger
    print_level::LogLevels = INFO
    file_print_level::LogLevels = INFO
    file::Union{IOStream,Nothing} = nothing
end

get_level(logger::Logger) = logger.print_level
get_file_level(logger::Logger) = logger.file_print_level
get_file(logger::Logger) = logger.file
finalize(logger::Logger) = logger.file != nothing && close(logger.file)

for (name,level,color) in [(:trace,TRACE,7),(:debug,DEBUG,6),(:info,INFO,256),(:notice,NOTICE,256),(:warn,WARN,5),(:error,ERROR,9)]
    @eval begin
        macro $name(logger,str)
            gl = $get_level
            gfl= $get_file_level
            gf = $get_file
            l = $level
            c = $color
            code = quote
                if $gl($logger) <= $l
                    if $c == 256
                        println($str)
                    else
                        printstyled($str,"\n",color=$c)
                    end
                end
                if $gf($logger) != nothing && $gfl($logger) <= $l
                    println($gf($logger),$str)
                end
            end
            esc(code)
        end
    end
end

# BLAS
const blas_num_threads = Ref{Int}(1)
function set_blas_num_threads(n::Integer;permanent::Bool=false)
    permanent && (blas_num_threads[]=n)
    BLAS.set_num_threads(n)
end
macro blas_safe_threads(args...)
    code = quote
        set_blas_num_threads(1)
        Threads.@threads($(args...))
        set_blas_num_threads(blas_num_threads[])
    end
    return esc(code)
end

# Type definitions
SubVector{Tv}=SubArray{Tv, 1, Vector{Tv}, Tuple{Vector{Int}}, false}
StrideOneVector{Tv}=Union{
    Vector{Tv},
    SubArray{Tv,1,Vector{Tv},Tuple{UnitRange{Int}},true},
    SubArray{Tv, 1, Matrix{Tv}, Tuple{Base.Slice{Base.OneTo{Int}}, Int}, true}
}


text:
module Day9
export day9
import ..Day1

function readInput(inputfile="inputs/input09.txt")
    [parse(Int, line) for line in eachline(inputfile)]
end

function day9(inputfile="inputs/input09.txt")
    values = readInput(inputfile)
    part1 = day9Part1(values)
    part2 = day9Part2(values, part1)
    (part1, part2)
end

function day9Part1(values)
    preamble = 25
    for (i, val) in enumerate(values[27:end])
        front,back = Day1.doublePointerSearch(sort(values[i:i+25]), val)
        if front == back # failure
            return val
        end
    end
end

function day9Part2(values, part1)
    lower = 1
    upper = 2
    s = values[lower]+values[upper]
    while s !== part1 && upper <= length(values)
        if s > part1
            lower += 1
        else
            upper += 1
        end
        s = sum(@view values[lower:upper])
    end
    if upper > length(values)
        error("[Day09 Part2]No contiguous range in input summing to $part1 found!")
    end
    return sum(extrema(@view values[lower:upper]))
end

end #module


text:
include("_common.jl")

using StatsPlots, OnlineStats

calib_start = 1826  # 1826
outflows = Dict()
climate = FULL_CLIMATE
calib_obs = CALIB_OBS[calib_start+1:end]

baseline_path = joinpath(FIG_PATH, "calib_baseline")

for app in APPROACHES
    display_name = replace(app, "_"=>" ")
    sn, n_id = setup_network(joinpath("baselines", "cotter_baseline_IHACRES_$(app).yml"))
    node = sn[n_id]
    run_basin!(sn, climate)

    calib_outflow = node.outflow[calib_start+1:CALIB_LN]
    gw_store = node.gw_store[calib_start+1:CALIB_LN]

    cmd = node.storage[calib_start+1:CALIB_LN]
    log_cmd = log.(cmd)
    max_cmd = maximum(log_cmd)
    quantiles = [0.2, 0.80]

    normal_dist = fit(Normal, log_cmd)
    μ = mean(normal_dist)
    σ = std(normal_dist)
    d_min = minimum(log_cmd)
    d_max = maximum(log_cmd)
    trunc_dist = truncated(Normal(μ, σ), d_min, d_max)

    thresholds = quantile(log_cmd, quantiles)
    @info "Log CMD thresholds: [0.1, 0.9 / min / max]" thresholds d_min d_max
    histogram(log_cmd, alpha=0.4, normalize=true,
              title="Distribution of Log CMD values\n($(display_name) calibrated baseline)",
              legend=:topleft,
              label=false,
              xlabel="Log CMD",
              ylabel="Density")
    StatsPlots.plot!(trunc_dist, linewidth=2, label="Fitted normal distribution")
    vline!(thresholds; color="orange", linewidth=1.5, label="State thresholds")
    for (thres, annot) in zip(thresholds, quantiles)
        annotate!(thres, 0.0015, text("q$annot", 8))
    end

    savefig(joinpath(baseline_path, "baseline_cmd_dist_$(app).png"))

    # Save performance over calibration period
    rmse = round(Streamfall.RMSE(calib_obs, calib_outflow), digits=2)
    mkge = round(Streamfall.mKGE(calib_obs, calib_outflow), digits=2)
    nse = round(Streamfall.NSE(calib_obs, calib_outflow), digits=2)
    nnse = round(Streamfall.NNSE(calib_obs, calib_outflow), digits=2)
    nmkge = round(Streamfall.NmKGE(calib_obs, calib_outflow), digits=2)
    mean_nmkge = round(Streamfall.mean_NmKGE(calib_obs, calib_outflow), digits=2)
    @info "$app calibration: "
    report_metrics(calib_obs, calib_outflow)

    plot(CALIB_DATES[calib_start+1:end], calib_obs, label="Historic", title="$app Calibration\n(RMSE: $rmse; NSE: $nse; KGE': $mkge)", legend=:topright)
    plot!(CALIB_DATES[calib_start+1:end], calib_outflow, label="Simulated", alpha=0.7, ylabel="Stream flow [ML]")
    savefig(joinpath(baseline_path, "baseline_calib_perf_$(app).png"))

    qqplot(calib_obs, calib_outflow, title="Calibration Period\n$display_name")
    xlabel!("Historic")
    ylabel!("Simulated")
    savefig(joinpath(baseline_path, "baseline_calib_perf_qq_$(display_name).png"))

    @info "CMD" app minimum(node.storage[2:end]) maximum(node.storage[2:end])
    @info "GW Store" app minimum(node.gw_store[2:end]) maximum(node.gw_store[2:end])
    @info "Quick" app minimum(node.quick_store[2:end]) maximum(node.quick_store[2:end])
    @info "Slow" app minimum(node.slow_store[2:end]) maximum(node.slow_store[2:end])

    # Save performance over validation period
    # valid_cmd = node.storage[CALIB_LN+1:end]
    valid_outflow = node.outflow[CALIB_LN+1:end]
    outflows[app] = valid_outflow

    rmse = round(Streamfall.RMSE(VALID_OBS, valid_outflow), digits=2)
    mkge = round(Streamfall.mKGE(VALID_OBS, valid_outflow), digits=2)
    nse = round(Streamfall.NSE(VALID_OBS, valid_outflow), digits=2)
    nnse = round(Streamfall.NNSE(VALID_OBS, valid_outflow), digits=2)
    nmkge = round(Streamfall.NmKGE(VALID_OBS, valid_outflow), digits=2)
    mean_nmkge = round(Streamfall.mean_NmKGE(VALID_OBS, valid_outflow), digits=2)
    @info "$app validation: "
    report_metrics(VALID_OBS, valid_outflow)

    plot(VALID_DATES, VALID_OBS, label="Historic", title="$display_name Validation\n(RMSE: $rmse; NSE: $nse; KGE': $mkge)", legend=:topright)
    plot!(VALID_DATES, valid_outflow, label="Simulated", alpha=0.7, ylabel="Stream flow [ML]")
    savefig(joinpath(baseline_path, "baseline_valid_perf_$(app).png"))

    qqplot(VALID_OBS, valid_outflow, title="Validation Period\n$display_name")
    xlabel!("Historic")
    ylabel!("Simulated")
    savefig(joinpath(baseline_path, "baseline_valid_perf_qq_$(app).png"))

    full_dates = vcat(CALIB_DATES, VALID_DATES)
    plot(full_dates[calib_start:CALIB_LN], CALIB_OBS[calib_start:end], label="Calibration", title="$display_name C&V", legend=:topright)
    plot!(full_dates[CALIB_LN+1:end], VALID_OBS, label="Validation", color="green")
    plot!(full_dates[calib_start:end], node.outflow[calib_start:end], label="Simulated", linestyle=:dashdot, color="red", alpha=0.7, ylabel="Stream flow [ML]")
    savefig(joinpath(baseline_path, "baseline_full_perf_$(app).png"))
end

# Ensemble results
ensemble_outflow = mean.(outflows["NNSE"] .+ outflows["mean_NmKGE"])
rmse = round(Streamfall.RMSE(VALID_OBS, ensemble_outflow), digits=2)
mkge = round(Streamfall.mKGE(VALID_OBS, ensemble_outflow), digits=2)
nse = round(Streamfall.NSE(VALID_OBS, ensemble_outflow), digits=2)

plot(VALID_DATES, VALID_OBS, label="Historic", title="Baseline Ensemble\n(RMSE: $rmse; NSE: $nse; KGE': $mkge)", legend=:topright)
plot!(VALID_DATES, ensemble_outflow, label="Ensemble", alpha=0.8)
savefig(joinpath(FIG_PATH, "baseline_valid_perf_ensemble_NNSE_meanNmKGE.png"))

# keeping the below as it plots out simulation for each bin
# climate = Climate(training, "_rain", "_temp")

# # Set up stream network
# network = YAML.load_file(joinpath(DATA_PATH, "cotter_baseline_calibrated_NNSE.yml"))
# nnse_sn = create_network("Calibrated NNSE", network)

# run_basin!(nnse_sn, climate)

# node = nnse_sn[1]
# nnse_outflow = node.outflow
# nnse_cmd = node.storage[2:end]  # skip the first entry as it is the initial condition

# normal_dist = fit(Normal, nnse_cmd)
# trunc_dist = truncated(Normal(mean(normal_dist), std(normal_dist)), 0.0, maximum(nnse_cmd))

# # save to file
# dist_data = Dict(
#     "model" => ["Baseline_NNSE", "Baseline_NmKGE", "Baseline_RMSE"],
#     "mean" => [mean(normal_dist), ],
#     "std" => [std(normal_dist), ],
#     "max" => [maximum(nnse_cmd), ],
# )

# nnse_bin_bounds = quantile.(trunc_dist, [0.5-0.25, 0.5+0.25, 1.0])
# histogram(nnse_cmd, alpha=0.4, normalize=true,
#           title="Distribution of CMD values\n(NNSE calibrated baseline)",
#           legend=:left,
#           label="CMD",
#           xlabel="Catchment Moisture Deficit",
#           ylabel="Probability")
# StatsPlots.plot!(trunc_dist, linewidth=2, label="Fitted normal distribution")
# vline!(nnse_bin_bounds; color="orange", linewidth=1.5, label="Ensemble bounds")

# savefig(joinpath(FIG_PATH, "NNSE_calibrated_dist.png"))

# ######

# # Set up stream network
# nmkge_network = YAML.load_file(joinpath(DATA_PATH, "cotter_baseline_calibrated_NmKGE.yml"))
# nmkge_sn = create_network("Calibrated NmKGE", nmkge_network)

# run_basin!(nmkge_sn, climate)

# node = nmkge_sn[1]
# nmkge_outflow = node.outflow
# nmkge_cmd = node.storage[2:end]

# normal_dist = fit(Normal, nmkge_cmd)
# trunc_dist = truncated(Normal(mean(normal_dist), std(normal_dist)), 0.0, maximum(nmkge_cmd))

# # Save data
# append!(dist_data["mean"], mean(normal_dist))
# append!(dist_data["std"], std(normal_dist))
# append!(dist_data["max"], maximum(nmkge_cmd))

# nmkge_bin_bounds = quantile.(trunc_dist, [0.5-0.25, 0.5+0.25, 1.0])
# histogram(nmkge_cmd, alpha=0.4, normalize=true,
#           title="Distribution of CMD values\n(NmKGE calibrated baseline)",
#           legend=:left,
#           label="CMD",
#           xlabel="Catchment Moisture Deficit",
#           ylabel="Probability")
# StatsPlots.plot!(trunc_dist, linewidth=2, label="Fitted normal distribution")
# vline!(nmkge_bin_bounds; color="orange", linewidth=1.5, label="Ensemble bounds")

# savefig(joinpath(FIG_PATH, "NmKGE_calibrated_dist.png"))

# ######

# # Set up stream network
# rmse_network = YAML.load_file(joinpath(DATA_PATH, "cotter_baseline_calibrated_RMSE.yml"))
# rmse_sn = create_network("Calibrated RMSE", rmse_network)

# run_basin!(rmse_sn, climate)

# node = rmse_sn[1]
# rmse_outflow = node.outflow
# rmse_cmd = node.storage[2:end]

# normal_dist = fit(Normal, rmse_cmd)
# trunc_dist = truncated(Normal(mean(normal_dist), std(normal_dist)), 0.0, maximum(rmse_cmd))

# # Save data
# append!(dist_data["mean"], mean(normal_dist))
# append!(dist_data["std"], std(normal_dist))
# append!(dist_data["max"], maximum(rmse_cmd))
# CSV.write(joinpath(DATA_PATH, "dist_data.csv"), DataFrame(dist_data))

# nmkge_bin_bounds = quantile.(trunc_dist, [0.5-0.25, 0.5+0.25, 1.0])
# histogram(rmse_cmd, alpha=0.4, normalize=true,
#           title="Distribution of CMD values\n(RMSE calibrated baseline)",
#           legend=:left,
#           label="CMD",
#           xlabel="Catchment Moisture Deficit",
#           ylabel="Probability")
# StatsPlots.plot!(trunc_dist, linewidth=2, label="Fitted normal distribution")
# vline!(nmkge_bin_bounds; color="orange", linewidth=1.5, label="Ensemble bounds")

# savefig(joinpath(FIG_PATH, "RMSE_calibrated_dist.png"))

# #######

# # Performance of baseline and ensemble of baselines
# hist_flow = training_flow["410730"]
# mean_outflow = mean([nnse_outflow, nmkge_outflow], dims=1)[1]
# @info "Performance of baseline and ensemble of baselines"
# @info "NNSE Calibrated" Streamfall.RMSE(hist_flow, nnse_outflow)
# @info "NmKGE Calibrated" Streamfall.RMSE(hist_flow, nmkge_outflow)
# @info "RMSE Calibrated" Streamfall.RMSE(hist_flow, rmse_outflow)
# @info "Mean of Ensemble" Streamfall.RMSE(hist_flow, mean_outflow)

# # Performance within each bin
# nnse_bin_membership = bin_membership(nnse_cmd, nnse_bin_bounds)
# @info "Performance within bins (NNSE)"
# for idx in 1:length(nnse_bin_bounds)+1
#     bin_member = findall(x -> x == idx, nnse_bin_membership)
#     @info "Bin $idx ($(length(bin_member)))"
#     @info "NNSE" Streamfall.RMSE(hist_flow[bin_member], nnse_outflow[bin_member])
#     @info "Ensemble" Streamfall.RMSE(hist_flow[bin_member], mean_outflow[bin_member])
# end

# ######
# println("------------")

# nmkge_bin_membership = bin_membership(nmkge_cmd, nmkge_bin_bounds)
# @info "Performance within bins (NmKGE)"
# for idx in 1:length(nmkge_bin_bounds)+1
#     bin_member = findall(x -> x == idx, nmkge_bin_membership)
#     @info "Bin $idx ($(length(bin_member)))"
#     if length(bin_member) == 0
#         @info "No members"
#         continue
#     end
#     @info "NmKGE" Streamfall.RMSE(hist_flow[bin_member], nmkge_outflow[bin_member])
#     @info "Ensemble" Streamfall.RMSE(hist_flow[bin_member], mean_outflow[bin_member])
# end

# plot(hist_flow, title="Baseline results", label="Observed",
#      xlabel="Time step", ylabel="Streamflow [ML]")
# num_bounds = length(nnse_bin_bounds)+1
# for idx in 1:num_bounds
#     bin_member = findall(x -> x == idx, nnse_bin_membership)
#     if length(bin_member) == 0
#         continue
#     end
#     scatter!(bin_member, nnse_outflow[bin_member],
#              markersize=2, label="Bin $idx",
#              markerstrokewidth=0,
#              alpha=0.6;
#              palette=cgrad(:matter, num_bounds, categorical = true))
# end

# savefig(joinpath(FIG_PATH, "NNSE_baseline_perf.png"))

# plot(hist_flow, title="Baseline results", label="Observed",
#      xlabel="Time step", ylabel="Streamflow [ML]")
# for idx in 1:num_bounds
#     bin_member = findall(x -> x == idx, nmkge_bin_membership)
#     if length(bin_member) == 0
#         continue
#     end
#     scatter!(bin_member, nmkge_outflow[bin_member],
#              markersize=2, label="Bin $idx",
#              markerstrokewidth=0,
#              alpha=0.6;
#              palette=cgrad(:matter, num_bounds, categorical = true))
# end

# savefig(joinpath(FIG_PATH, "NmKGE_baseline_perf.png"))

# v_flow = validation_flow["410730"]
# reset!(nnse_sn)
# v_climate = Climate(validation, "_rain", "_temp")
# run_basin!(nnse_sn, v_climate)

# plot(v_flow, label="Observed", alpha=0.5, title="Validation Period")
# nnse_outflow = nnse_sn[1].outflow
# rmse = round(Streamfall.RMSE(v_flow, nnse_outflow), digits=4)
# plot!(nnse_outflow, label="NNSE calibrated: RMSE = $rmse", linestyle=:dash, alpha=0.9)

# ###
# reset!(nmkge_sn)
# run_basin!(nmkge_sn, v_climate)
# nmkge_outflow = nmkge_sn[1].outflow
# rmse = round(Streamfall.RMSE(v_flow, nmkge_outflow), digits=4)
# plot!(nmkge_outflow, label="NmKGE calibrated: RMSE = $rmse", linestyle=:dashdotdot, alpha=0.8)

# ####
# reset!(rmse_sn)
# run_basin!(rmse_sn, v_climate)
# rmse_outflow = rmse_sn[1].outflow
# rmse = round(Streamfall.RMSE(v_flow, rmse_outflow), digits=4)
# plot!(rmse_outflow, label="RMSE calibrated: RMSE = $rmse", linestyle=:dashdotdot, alpha=0.7)

# # mean_outflow = mean([nnse_outflow, nmkge_outflow], dims=1)[1]
# # rmse = round(Streamfall.RMSE(v_flow, mean_outflow), digits=4)
# # plot!(mean_outflow, label="Ensemble mean: RMSE = $rmse")

# savefig(joinpath(FIG_PATH, "calibrated_performance.png"))


max_stars_repo_name: MechaJulia/MechaJulia

max_issues_repo_name: MechaJulia/MechaJulia

max_forks_repo_name: MechaJulia/MechaJulia


text:
int main() {

  double[] vector = new double[4];
  double[][] matrix = new double[3][4];

  int i=0;
  while (i<matrix.length) {
    int j=0;
    while (j<matrix[0].length) {
       matrix[i][j] = 5.0;
       j++;
    }
    i++;
  }

  i=0;
  while (i<vector.length-1) {
     vector[i] = 3.0;
     i++;
  }

  matrix[0] = vector;

  int j = 0;
  while (j<vector.length) {
    matrix[1][j] = vector[j]+1.0;
    j++;
  }

  for (double[] x : matrix)
    for (double y : x)
      printDouble(y);

  return 0 ;
}

max_stars_repo_name: carlostome/Javalette-Compiler

max_issues_repo_name: carlostome/Javalette-Compiler

max_forks_repo_name: carlostome/Javalette-Compiler


text:
# SPDX-License-Identifier: MIT

using ArgParse
using HDF5

function parse_commandline()
    s = ArgParseSettings()
    @add_arg_table s begin
        "--compression-level", "-c"
            arg_type = Int
            default = 0
            help = "compression level of the HDF5 archive"
        "out"
            help = "output HDF5 archive"
            required = true
        "archives"
            nargs = '+'
            help = "archives to concatenate"
            required = true
    end
    parse_args(s)
end

function main(args)
    archives = args["archives"]
    h5open(args["out"], "w") do f

        # Use the first archive of the list to add the attributes to
        # the output archive.
        h5open(archives[1], "r") do farchive
            for k in keys(attributes(farchive))
                attributes(f)[k] = read(attributes(farchive)[k])
            end
        end

        # Copy the data of each archive in the new one.
        for arch in archives
            h5open(arch, "r") do farchive
                for k in keys(farchive["/"])
                    write(f, k, read(farchive[k]))
                end
            end
        end
    end
end

args = parse_commandline()
main(args)

max_stars_repo_name: MartinKocour/SpeechLab

max_issues_repo_name: MartinKocour/SpeechLab

max_forks_repo_name: MartinKocour/SpeechLab


text:
module JSONRenderer

import DataFrames
import JSON

import Dance.Configuration
import Dance.Logger
import Dance.PayloadUtils

const HTTP_STATUS_UNAUTHORIZED = "Unauthorized"
const HTTP_STATUS_PAYMENT_REQUIRED = "Payment Required"
const HTTP_STATUS_FORBIDDEN = "Forbidden"
const HTTP_STATUS_REQUEST_TIMEOUT = "Request Timeout"

const HTTP_STATUS_ERROR_CODES = Dict{String, Int16}(
    HTTP_STATUS_UNAUTHORIZED => 401,
    HTTP_STATUS_PAYMENT_REQUIRED => 402,
    HTTP_STATUS_FORBIDDEN => 403,
    HTTP_STATUS_REQUEST_TIMEOUT => 408
)

"""
    render(;headers::Dict{String, String}, status_code::Int16, data::Union{DataFrames.DataFrame, Dict}) :: Dict{Symbol, Union{Dict{String, String}, Int16, String}}

Output JSON renderer

Though `status code` is pre-supplied, can become 500

Cannot set `headers type to `Dict{String, String} here, as can be blank (status_code !=200)
"""
function render(;headers::Dict{String, String}, status_code::Int16, data::Union{DataFrames.DataFrame, Dict}) :: Dict{Symbol, Union{Dict{String, String}, Int16, String}}
    body::Union{Array{Any,1}, Dict} = Dict()

    # One of HTTP_STATUS_ERROR_CODES is supplied as data{Dict}
    if status_code==200 && isa(data, Dict)
        if length(data)==1 && :error in keys(data)
            try
                status_code = HTTP_STATUS_ERROR_CODES[data[:error]]
                body = Dict(
                    :error => string(status_code, ": ", data[:error])
                )
            catch e
                status_code = 500
                body = Dict(
                    :error => "500:  Internal Server Error"
                )

                data_string::String = string(data)
                Logger.log("500 Internal Server Error when rendering page with $data_string: $e")
            end
        else
            body = data
        end
    else
        if isa(data, DataFrames.DataFrame)
            body = PayloadUtils.convert_dataframe_to_array(data)
        else
            body = data
        end
    end

    headers["Access-Control-Allow-Origin"] = Configuration.Settings[:api_access_control_allow_origin]
    return Dict(
        :headers => headers,
        :status_code => status_code,
        :content_type => "application/json",
        :body => JSON.json(body)
    )
end

end


text:
"""
Proportional-integral-derivative (PID) controller.
# Notes
The derivative term is obtained via second-order filter.
For more details, see [1] and singular perturbation theory.
# References
[1] J. A. Farrell, M. Polycarpou, M. Sharma, and W. Dong, “Command Filtered Backstepping,” IEEE Transactions on Automatic Control, vol. 54, no. 6, pp. 1391–1395, Jun. 2009, doi: 10.1109/TAC.2009.2015562.
"""
struct PID <: AbstractEnv
    K_P
    K_I
    K_D
    ω_n
    ζ
    windup_limit
    function PID(K_P, K_I, K_D; ω_n=1e1, ζ=5e-1, windup_limit=Inf)
        @assert windup_limit > 0.0
        @assert ω_n > 0.0
        @assert ζ > 0.0
        new(K_P, K_I, K_D, ω_n, ζ, windup_limit)
    end
end

"""
∫e: integral of the error
ê: the estimate of the error
ė̂: the estimate of the error rate
"""
function State(controller::PID)
    function state(; ∫e=zeros(1), ê=zeros(1), ė̂=zeros(1))
        ComponentArray(∫e=∫e, ê=ê, ė̂=ė̂)
    end
end

function Dynamics!(controller::PID)
    @unpack ω_n, ζ, windup_limit = controller
    @Loggable function dynamics!(dx, x, p, t; e)
        @assert !(typeof(e) <: Number)  # make sure that `e` is a 1-d array
        @unpack ∫e, ê, ė̂ = x
        @onlylog e, ∫e, ê, ė̂
        if norm(∫e) < windup_limit
            dx.∫e .= e
        elseif norm(∫e) > windup_limit && all(sign.(e) .* sign.(∫e) .< zero.(e))
            dx.∫e .= e
        else
            dx.∫e .= zero.(e)
        end
        dx.ê .= ω_n*ė̂
        dx.ė̂ = -2*ζ*ω_n*ė̂ - ω_n*(ê - e)
    end
end

function Command(controller::PID)
    @unpack K_P, K_I, K_D = controller
    function command(e, ∫e, ė̂)
        -(K_P*e + K_I*∫e + K_D*ė̂)
    end
end


text:
# OptimalTransport.jl -- optimal transportation algorithms for Julia
# See prettyprinted documentation at https://juliaoptimaltransport.github.io/OptimalTransport.jl/dev/
#

module OptimalTransport

using Distances
using LinearAlgebra
using IterativeSolvers, SparseArrays
using LogExpFunctions: LogExpFunctions
using MathOptInterface
using Distributions
using PDMats
using QuadGK
using NNlib: NNlib
using StatsBase: StatsBase

export SinkhornGibbs, SinkhornStabilized, SinkhornEpsilonScaling
export SinkhornBarycenterGibbs
export QuadraticOTNewton

export sinkhorn, sinkhorn2
export emd, emd2
export sinkhorn_stabilized, sinkhorn_stabilized_epsscaling, sinkhorn_barycenter
export sinkhorn_unbalanced, sinkhorn_unbalanced2
export quadreg
export ot_cost, ot_plan, wasserstein, squared2wasserstein

const MOI = MathOptInterface

include("distances/bures.jl")
include("utils.jl")
include("exact.jl")
include("wasserstein.jl")

include("entropic/sinkhorn.jl")
include("entropic/sinkhorn_gibbs.jl")
include("entropic/sinkhorn_stabilized.jl")
include("entropic/sinkhorn_epsscaling.jl")
include("entropic/sinkhorn_unbalanced.jl")
include("entropic/sinkhorn_barycenter.jl")
include("entropic/sinkhorn_barycenter_gibbs.jl")
include("entropic/sinkhorn_solve.jl")

include("quadratic.jl")
include("quadratic_newton.jl")

end


text:
using Base: @deprecate

@deprecate by{T,N}(ta::TimeArray{T,N}, t::Int; period::Function=day) when(ta, period, t)
@deprecate by{T,N}(ta::TimeArray{T,N}, t::String; period::Function=day) when(ta, period, t)

@deprecate to(ta::TimeArray, y::Int, m::Int, d::Int) to(ta, Date(y, m, d))
@deprecate from(ta::TimeArray, y::Int, m::Int, d::Int) from(ta, Date(y, m, d))

@deprecate findall(ta::TimeArray) find(ta)

@deprecate collapse{T,N,D}(ta::TimeArray{T,N,D}, timestamp::Function; period::Function=week) collapse(ta, period, timestamp)

# since julia 0.6

# deprecate non-dot function due to 0.6 syntactic loop fusion
for f ∈ (:^, :/, :abs, :sign, :sqrt, :cbrt,
         :log, :log2, :log10, :log1p,
         :exp, :exp2, :exp10, :expm1,
         :cos, :sin, :tan, :cosd, :sind, :tand,
         :acos, :asin, :atan, :acosd, :asind, :atand,
         :isnan, :isinf)
    @eval import Base: $f
    @eval @deprecate $f(ta::TimeArray, args...) $f.(ta, args...)
end

for f ∈ (:+, :-, :*, :%,
         :|, :&, :<, :>, :(==), :(!=), :>=, :<=)
    @eval import Base: $f
    @eval @deprecate $f(ta::TimeArray, args...) $f.(ta, args...)
    @eval @deprecate $f(n::Number, ta::TimeArray) $f.(n, ta)
end

# non-dot operators
import Base: $, !, ~

@deprecate ($)(ta1::TimeArray, ta2::TimeArray) xor.(ta1, ta2)
@deprecate ($)(n::Integer, ta::TimeArray) xor.(n, ta)
@deprecate ($)(ta::TimeArray, n::Integer) xor.(ta, n)

@deprecate ~(ta::TimeArray) .~(ta)
@deprecate !(ta::TimeArray) .!(ta)


text:
using TeeStreams, Test, SimpleBufferStream, HTTP,
      CodecZlib, CodecZstd, SHA, Tar

function write_things(io)
    # String
    write(io, "hello, string\n")
    # Single byte and byte arrays
    for b in Vector{UInt8}("hello, ")
        write(io, b)
    end
    write(io, Vector{UInt8}("bytes\n"))
    # Write from file (IOStream)
    f, s = mktemp()
    write(s, "hello, file\n"); close(s)
    open(io1 -> write(io, io1), f)
    # Read from external process (Base.Process)
    open(`$(Base.julia_cmd()[1]) -e 'println("hello, process")'`) do proc
        write(io, proc)
    end
    # print(ln)
    print(io, "hello, "); println(io, "print")
    return io
end

@testset "TeeStreams" begin
    # write(::TeeStream, x) for different x and tee'd streams

    correct = String(take!(write_things(IOBuffer())))

    # tee = TeeStream(...)
    iob = IOBuffer()
    ioc = IOContext(IOBuffer())
    f, iof = mktemp()
    bs = BufferStream()
    tee = TeeStream(iob, ioc, iof, bs)
    write_things(tee)
    close(iof); close(bs)
    @test String(take!(iob)) == String(take!(ioc.io)) ==
          read(f, String) == read(bs, String) == correct
    try
        close(tee)
        write_things(tee)
    catch err
        @test err isa CompositeException
        @test length(err.exceptions) == 4
        @test all(x -> x isa TaskFailedException, err.exceptions)
    end

    # TeeStream(...) do tee
    mktempdir() do tmpd;
        f1 = joinpath(tmpd, "file")
        io1 = open(f1, "w")
        f2 = joinpath(tmpd, "file2")
        io2 = open(f2, "w")
        @test isopen(io1)
        @test isopen(io2)
        TeeStream(io1, io2) do tee
            write_things(tee)
        end
        @test !isopen(io1)
        @test !isopen(io2)
        @test read(f1, String) == read(f2, String) == correct
    end

    # Test some integration with other packages
    mktempdir() do tmpd
        url = "https://julialang-s3.julialang.org/bin/linux/x64/1.5/julia-1.5.3-linux-x86_64.tar.gz"
        expected_shasum = "f190c938dd6fed97021953240523c9db448ec0a6760b574afd4e9924ab5615f1"

        #          HTTP.Stream
        #               +
        #               |
        #      +--------+----------+
        #      |                   |
        #      v                   v
        # BufferStream      GzipDecompressor
        #      +                   +
        #      |                   |
        #      v                   v
        # SHA.sha256          BufferStream
        #                          +
        #                          |
        #                          v
        #                     Tar.rewrite
        #                          +
        #                          |
        #                 +--------+---------+
        #                 |                  |
        #                 v                  v
        #            GzipCompressor   ZstdCompressor

        buffer_shasum = BufferStream()
        buffer_tar = BufferStream()
        decompressor = GzipDecompressorStream(buffer_tar)
        tee = TeeStream(buffer_shasum, decompressor)
        compressors = TeeStream(
            GzipCompressorStream(open(joinpath(tmpd, "julia.tar.gz"), "w")),
            ZstdCompressorStream(open(joinpath(tmpd, "julia.tar.zst"), "w")),
        )
        @sync begin
            dl_task = @async (HTTP.get(url; response_stream=tee); close(tee))
            sha_task = @async bytes2hex(SHA.sha256(buffer_shasum))
            tar_task = @async Tar.rewrite(buffer_tar, compressors)
            @test fetch(sha_task) == expected_shasum
        end
    end

end #testset


text:
module Config
using YAML
using Logging
using PaddedViews
using Distributions
using ArgParse

include("functions.jl")
functions = Array{Function}(0)

function init(config::Dict)
    for k in keys(config)
        if k == "functions"
            append!(functions, load_functions(config["functions"]))
        else
            if config[k] != nothing
                eval(parse(string(k, "=", config[k])))
            end
        end
    end
end

function bloat()
    ((mutate_method in [:mixed_node_mutate, :adaptive_node_mutate,
                        :mixed_subtree_mutate, :adaptive_subtree_mutate]) ||
     (crossover_method in [:output_graph_crossover, :subgraph_crossover]))
end

function init(file::String)
    init(YAML.load_file(file))
end

function reset()
    empty!(functions)
end

function add_arg_settings!(s::ArgParseSettings)

    mutations = [":gene_mutate", ":mixed_node_mutate", ":mixed_subtree_mutate",
                 ":adaptive_node_mutate", ":adaptive_subtree_mutate"]

    crossovers = [":single_point_crossover", ":random_node_crossover",
                  ":aligned_node_crossover", ":proportional_crossover",
                  ":output_graph_crossover", ":subgraph_crossover"]

    distances = [":genetic_distance", ":positional_distance",
                 ":constant_functional_distance", ":random_functional_distance",
                 ":active_distance"]

    @add_arg_table s begin
        "--mutate_method"
            default = nothing
            range_tester = (x->x ∈ mutations)
            help = "mutation method; must be one of " * join(mutations, ", ", " or ")
        "--crossover_method"
            default = nothing
            range_tester = (x->x ∈ crossovers)
            help = "crossover method; must be one of " * join(crossovers, ", ", " or ")
        "--distance_method"
            default = nothing
            range_tester = (x->x ∈ distances)
            help = "distance method; must be one of " * join(distances, ", ", " or ")
    end

    params = ["input_start", "recurrency", "input_mutation_rate",
        "output_mutation_rate", "node_mutation_rate", "node_size_delta",
        "modify_mutation_rate", "ga_elitism_rate", "ga_crossover_rate",
        "ga_mutation_rate"]

    for p in params
        add_arg_table(s, ["--$p"], Dict(:help=>"Parameter: $p", :arg_type=>Float64))
    end

    for p in ["lambda", "ga_population", "starting_nodes", "static_node_size",
              "node_size_cap", "total_evals"]
        add_arg_table(s, ["--$p"], Dict(:help=>"Parameter: $p", :arg_type=>Int64))
    end

    for p in ["active_mutate", "weights", "save_best"]
        add_arg_table(s, ["--$p"], Dict(:help=>"Parameter: $p", :arg_type=>Bool))
    end
    s
end

function get_arg_settings()
    s = ArgParseSettings()
    add_arg_settings!(s)
    s
end

function to_string()
    @sprintf(
        "%s %s %s %s %d %d %d %0.3f %0.3f %0.3f %0.3f %0.3f %0.3f %0.3f %d %d %0.3f %0.3f %0.3f",
        string(mutate_method), string(active_mutate), string(crossover_method),
        string(weights), starting_nodes, static_node_size, node_size_cap,
        input_start, recurrency, input_mutation_rate,
        output_mutation_rate, node_mutation_rate, node_size_delta,
        modify_mutation_rate, lambda, ga_population, ga_elitism_rate,
        ga_crossover_rate, ga_mutation_rate)

end

# append!(functions, [f_input])
export init, reset
end


text:
import Base.hcat, Base.vcat, Base.hvcat

hcat(A::AbstractLinearOperator, B::AbstractMatrix) = hcat(A, LinearOperator(B))

hcat(A::AbstractMatrix, B::AbstractLinearOperator) = hcat(LinearOperator(A), B)

function hcat_prod!(res::AbstractVector, A::AbstractLinearOperator{T}, B::AbstractLinearOperator{T},
                    Ancol::I, nV::I, v::AbstractVector, α, β) where {T,I<:Integer}
  mul!(res, A, view(v, 1:Ancol), α, β)
  mul!(res, B, view(v, (Ancol+1): nV), α, one(T))
end

function hcat_tprod!(res::AbstractVector, A::AbstractLinearOperator{T}, B::AbstractLinearOperator{T},
                     Ancol::I, nV::I, u::AbstractVector, α, β) where {T,I<:Integer}
  mul!(view(res, 1:Ancol), transpose(A), u, α, β)
  mul!(view(res, (Ancol+1): nV), transpose(B), u, α, β)
end

function hcat_ctprod!(res::AbstractVector, A::AbstractLinearOperator{T}, B::AbstractLinearOperator{T},
                      Ancol::I, nV::I, w::AbstractVector, α, β) where {T,I<:Integer}
  mul!(view(res, 1:Ancol), adjoint(A), w, α, β)
  mul!(view(res, (Ancol+1): nV), adjoint(B), w, α, β)
end

function hcat(A::AbstractLinearOperator, B::AbstractLinearOperator)
  size(A, 1) == size(B, 1) || throw(LinearOperatorException("hcat: inconsistent row sizes"))

  nrow = size(A, 1)
  Ancol, Bncol = size(A, 2), size(B, 2)
  ncol = Ancol + Bncol
  S = promote_type(eltype(A), eltype(B))
  MvAB = res = Vector{S}(undef, nrow)
  MtuAB = res = Vector{S}(undef, nrow)
  MawAB = Vector{S}(undef, ncol)

  prod = @closure (res, v, α, β) -> hcat_prod!(res, A, B, Ancol, Ancol+Bncol, v, α, β)
  tprod = @closure (res, u, α, β) -> hcat_tprod!(res, A, B, Ancol, Ancol+Bncol, u, α, β)
  ctprod = @closure (res, w, α, β) -> hcat_ctprod!(res, A, B, Ancol, Ancol+Bncol, w, α, β)
  LinearOperator{S}(nrow, ncol, false, false, prod, tprod, ctprod, MvAB, MtuAB, MawAB)
end

function hcat(ops::AbstractLinearOperator...)
  op = ops[1]
  for i = 2:length(ops)
    op = [op ops[i]]
  end
  return op
end

vcat(A::AbstractLinearOperator, B::AbstractMatrix) = vcat(A, LinearOperator(B))

vcat(A::AbstractMatrix, B::AbstractLinearOperator) = vcat(LinearOperator(A), B)

function vcat_prod!(res::AbstractVector, A::AbstractLinearOperator{T}, B::AbstractLinearOperator{T},
                    Anrow::I, nV::I, u::AbstractVector, α, β) where {T,I<:Integer}
  mul!(view(res, 1:Anrow), A, u, α, β)
  mul!(view(res, (Anrow+1): nV), B, u, α, β)
end

function vcat_tprod!(res::AbstractVector, A::AbstractLinearOperator{T}, B::AbstractLinearOperator{T},
                     Anrow::I, nV::I, v::AbstractVector, α, β) where {T,I<:Integer}
  mul!(res, transpose(A), view(v, 1:Anrow), α, β)
  mul!(res, transpose(B), view(v, (Anrow+1): nV), α, one(T))
end

function vcat_ctprod!(res::AbstractVector, A::AbstractLinearOperator{T}, B::AbstractLinearOperator{T},
                     Anrow::I, nV::I, v::AbstractVector, α, β) where {T,I<:Integer}
  mul!(res, adjoint(A), view(v, 1:Anrow), α, β)
  mul!(res, adjoint(B), view(v, (Anrow+1): nV), α, one(T))
end

function vcat(A::AbstractLinearOperator, B::AbstractLinearOperator)
  size(A, 2) == size(B, 2) || throw(LinearOperatorException("vcat: inconsistent column sizes"))

  Anrow, Bnrow = size(A, 1), size(B, 1)
  nrow = Anrow + Bnrow
  ncol = size(A, 2)
  S = promote_type(eltype(A), eltype(B))
  MvAB = Vector{S}(undef, nrow)
  MtuAB = Vector{S}(undef, ncol)
  MawAB = Vector{S}(undef, ncol)

  prod! = @closure (res, v, α, β) -> vcat_prod!(res, A, B, Anrow, Anrow+Bnrow, v, α, β)
  tprod! = @closure (res, u, α, β) -> vcat_tprod!(res, A, B, Anrow, Anrow+Bnrow, u, α, β)
  ctprod! = @closure (res, w, α, β) -> vcat_ctprod!(res, A, B, Anrow, Anrow+Bnrow, w, α, β)
  return LinearOperator{S}(nrow, ncol, false, false, prod!, tprod!, ctprod!, MvAB, MtuAB, MawAB)
end

function vcat(ops::AbstractLinearOperator...)
  op = ops[1]
  for i = 2:length(ops)
    op = [op; ops[i]]
  end
  return op
end

# Removed by https://github.com/JuliaLang/julia/pull/24017
function hvcat(rows::Tuple{Vararg{Int}}, ops::AbstractLinearOperator...)
  nbr = length(rows)
  rs = Array{AbstractLinearOperator, 1}(undef, nbr)
  a = 1
  for i = 1:nbr
    rs[i] = hcat(ops[a:(a - 1 + rows[i])]...)
    a += rows[i]
  end
  vcat(rs...)
end


text:
using Distributed
using Plots
addprocs(48)
nprocs()
using SharedArrays
@everywhere using FRGRealTime, DelimitedFiles, Dierckx

kdata=sol.t[1:1:end]|>reverse
lamdata=sol[1,:][1:1:end]|>reverse
m2data=-sol[2,:][1:1:end]|>reverse

@eval @everywhere kdata=$kdata
@eval @everywhere lamdata=$lamdata
@eval @everywhere m2data=$m2data

@everywhere kdata=collect(1:1.0:1000)
@everywhere lamdata=-10*one.(kdata)
@everywhere m2data=10.0^2 *one.(kdata)

@everywhere lamfun=Spline1D(kdata,lamdata)
@everywhere m2fun=Spline1D(kdata,m2data)
v1=SharedArray(collect(0.1:4:50.0))

@everywhere kdata=readdlm("/home/wjfu1/tyy/CTP-fRG-Test/Real-Time-data/realtime_zero_k.dat")[:,1]|>reverse
@everywhere lamdata=readdlm("/home/wjfu1/tyy/CTP-fRG-Test/Real-Time-data/realtime_zero_lamdak.dat")[:,1]|>reverse
@everywhere m2data=readdlm("/home/wjfu1/tyy/CTP-fRG-Test/Real-Time-data/realtime_zero_m2k.dat")[:,1]|>reverse
@everywhere lamfun=Spline1D(kdata,0.5*lamdata)
@everywhere m2fun=Spline1D(kdata,m2data)

plot(x->m2fun(x),1.0,800.0)
plot(x->lamfun(x),1.0,800.0)

p0v1=SharedArray(collect(0.01:0.5:400.0))
psv1=SharedArray(collect(0.01:2.0:300.0))

Rev1=[(i,j) for j in psv1 for i in p0v1]

@time outvRe2=pmap(
    pv -> FRGRealTime.fastpropReintqs_All(
        pv[1],
        pv[2],
        temp_T,
        1.0,
        800.0,
        4.0,
        m2fun,
        lamfun,
        rtol = 1e-5,
        atol = 1e-10,
        maxevals=100000000,
        initdiv=20,
    ),Imv1
)
# interrupt(collect(1:200))
# p0v1=readdlm("$dir/p0.dat")[:,1]
# outvRe=readdlm("$dir/Re.dat")[:,1]

FRGRealTime.fastpropReintqs_All(
        1.01,
        0.000001,
        20.44,
        1.0,
        800.0,
        4.0,
        m2fun,
        lamfun,
        rtol = 1e-4,
        atol = 1e-5,
        maxevals=10000000,
        initdiv=40,
    )

plot(p0v1,outvRe2)

interrupt(collect(1:100))

plot!(p0v1,p0v1.^2)

interrupt(collect(1:400))

dir="/home/wjfu1/tyy/FRG-RealTime-data/Gamma2Im"
dir2="/home/wjfu1/tyy/FRG-RealTime-data/Re/ps=0"

# writedlm("$dir/Rep0.dat",v1)
# writedlm("$dir/ReGamm2_ps=0.dat",outvRe)

Imgamm2fun_All=Spline1D(
    readdlm("$dir/Imp0.dat")[:, 1],
    readdlm("$dir/propImI1p0+q0.dat")[:, 1] +
    readdlm("$dir/propImI1p0-q0.dat")[:, 1] +
    readdlm("$dir/propImI2p0+q0.dat")[:, 1] +
    readdlm("$dir/propImI2p0-q0.dat")[:, 1],
)

Imgamm2fun_1=Spline1D(
    readdlm("$dir/Imp0.dat")[:, 1],
    readdlm("$dir/propImI1p0+q0.dat")[:, 1]
)

Imgamm2fun_2=Spline1D(
    readdlm("$dir/Imp0.dat")[:, 1],
    readdlm("$dir/propImI1p0-q0.dat")[:, 1]
)
Imgamm2fun_3=Spline1D(
    readdlm("$dir/Imp0.dat")[:, 1],
    readdlm("$dir/propImI2p0+q0.dat")[:, 1]
)
Imgamm2fun_4=Spline1D(
    readdlm("$dir/Imp0.dat")[:, 1],
    readdlm("$dir/propImI2p0-q0.dat")[:, 1]
)

Regamm2fun=Spline1D(readdlm("$dir2/Rep0.dat")[:, 1],
    readdlm("$dir2/ReGamm2_ps=0.dat")[:, 1],
)

plot(x->Imgamm2fun_1(x),1.0,50.0)

plot(x->-(Imgamm2fun_3(x)+Imgamm2fun_4(x))/((Imgamm2fun_3(x)+Imgamm2fun_4(x))^2+Regamm2fun(x)^2),1.0,400.0,yaxis=:log)

rho(x)=Imgamm2fun_All(x) / (Imgamm2fun_All(x)^2 + Regamm2fun(x)^2)

mv1=collect(1.0:0.1:400.0)

mv2=collect(202.9:0.1:400.0)

plot(
    mv1,rho.(mv1),
    yaxis = :log,
)

plot!(
    mv2,sign.(rho.(mv2)).*log.(abs.(rho.(mv2)))
    # yaxis = :log,
)

sign(x)*log(1+abs(x)/10^C)

plot(v1,outvRe/145^2,legend=:topleft)


text:
# Test wrapping a DiGraph to specialize behavior.

struct MyGraph
    graph::DiGraph

    MyGraph() = new(DiGraph())
end

# Delegate:

NahaGraphs.edges(g::MyGraph) = edges(g.graph)
NahaGraphs.nodes(g::MyGraph) = nodes(g.graph)
NahaGraphs.add_edge!(g::MyGraph, edge) = add_edge!(g.graph, edge)
NahaGraphs.remove_edge!(g::MyGraph, edge) = remove_edge!(g.graph, edge)
NahaGraphs.query(g::MyGraph, from, to) = query(g.graph, from, to)

# Delegation for Dict support:
Base.keys(g::MyGraph) = keys(g.graph)
Base.values(g::MyGraph) = values(g.graph)
Base.haskey(g::MyGraph, key) = haskey(g.graph, key)
Base.getindex(g::MyGraph, key) = getindex(g.graph, key)


text:
using Compat
using Distributions
using PDMats

using LinearAlgebra

struct NormalWishart <: ContinuousMultivariateDistribution
    dim::Int
    zeromean::Bool
    mu::Vector{Float64}
    kappa::Float64
    #Tchol::Cholesky{Float64}  # Precision matrix (well, sqrt of one)
    T::Matrix{Float64}
    nu::Float64

    function NormalWishart(mu::Vector{Float64}, kappa::Real,
                                  T::Matrix{Float64}, nu::Real)
        # Probably should put some error checking in here
        d = length(mu)
        zmean::Bool = true
        for i = 1:d
            if mu[i] != 0.
                zmean = false
                break
            end
        end
        @compat new(d, zmean, mu, Float64(kappa), collect(Symmetric(T)), Float64(nu))
    end
end

function NormalWishart(mu::Vector{Float64}, kappa::Real,
                       T::AbstractMatrix{Float64}, nu::Real)
    NormalWishart(mu, kappa, collect(T), nu)
end

import Distributions.rand

function rand(nw::NormalWishart)
    Lam = rand(Wishart(nw.nu, nw.T))
    mu = rand(MvNormal(nw.mu, collect(inv(Symmetric(Lam))) ./ nw.kappa))
    return (mu, Lam)
end


text:
using ForwardDiff

export  FermiOccupation,
        ∂ϵFermiOccupation,
        ∂μFermiOccupation,
        ∂TFermiOccupation

function FermiOccupation(ϵ,μ,T)
    kb = 8.6*10^(-5) #eV / K
    1/(exp( (ϵ-μ)/(kb*T)) +1 )
end

∂ϵFermiOccupation(ϵ,μ,T) = ForwardDiff.derivative( x -> FermiOccupation(x,μ,T))
∂μFermiOccupation(ϵ,μ,T) = ForwardDiff.derivative( x -> FermiOccupation(ϵ,x,T))
∂TFermiOccupation(ϵ,μ,T) = ForwardDiff.derivative( x -> FermiOccupation(ϵ,μ,x))

function Dx(kx,ky,Ham!,Matrix1,Matrix2,diff)
    dϵ = 1<<15
    ϵ = 1/dϵ
    Ham!(kx+ϵ,ky,Matrix1) ; Ham!(kx-ϵ,ky,Matrix2)
    @. diff = (Matrix1-Matrix2)*dϵ/2
end

function Dy(kx,ky,Ham!,Matrix1,Matrix2,diff)
    dϵ = 1<<15
    ϵ = 1/dϵ
    Ham!(kx,ky+ϵ,Matrix1) ; Ham!(kx,ky-ϵ,Matrix2)
    @. diff = (Matrix1-Matrix2)*dϵ/2
end


text:
bin2int(s::String) = parse(Int, s, base = 2)

function hex_to_bin(c::Char)
    i = parse(Int, string(c), base = 16)
    string(i, base = 2, pad = 4)
end

slurp_string(s) = string(hex_to_bin.(collect(s))...)
slurp_file(filename) = slurp_string(readlines(filename)[1])

function gobble_literal(s)
    literal = ""
    while true
        f = s[1]
        next = s[2:5]
        literal *= next
        if f == '0'
            return bin2int(literal), s[6:end]
        end
        s = s[6:end]
    end
end

gobble_literal("101111111000101000") == (2021, "000")

struct Literal
    version::Int
    value::Int
end

struct Operator
    version::Int
    type_id::Int
    subpackets::Vector{Union{Operator,Literal}}
end

function parse_header(s)
    version = bin2int(s[1:3])
    type_id = bin2int(s[4:6])
    remaining = s[7:end]
    version, type_id, remaining
end

function gobble_packet(s)
    version, type_id, remaining = parse_header(s)
    if type_id == 4
        value, remaining = gobble_literal(remaining)
        return Literal(version, value), remaining
    else
        length_type_id = remaining[1]
        if length_type_id == '1'
            num_packets = bin2int(remaining[2:12])
            subpackets, remaining = parse_n_subpackets(num_packets, remaining[13:end])
            return Operator(version, type_id, subpackets), remaining
        end
        if length_type_id == '0'
            l = bin2int(remaining[2:16])
            subpackets, remaining = parse_subpackets_l(l, remaining[17:end])
            return Operator(version, type_id, subpackets), remaining
        end
    end
end

function parse_subpackets_l(l, s)
    result = []
    to_parse = s[1:l]
    remaining = s[l+1:end]
    while length(to_parse) > 7
        next, to_parse = gobble_packet(to_parse)
        push!(result, next)
    end
    result, remaining
end

function parse_n_subpackets(n, s)
    result = []
    for _ = 1:n
        packet, s = gobble_packet(s)
        push!(result, packet)
    end
    result, s
end

parse_subpackets_l(27, "1101000101001010010001001000000000")

eg, remaining = gobble_packet("11101110000000001101010000001100100000100011000001100000")

gobble_packet("00111000000000000110111101000101001010010001001000000000")

gobble_packet(slurp_file("16.txt"))

get_versions(o::Operator) = get_versions(o, [])

function get_versions(o::Operator, so_far)
    push!(so_far, o.version)
    for item in o.subpackets
        for version in get_versions(item)
            push!(so_far, version)
        end
    end
    so_far
end

get_versions(l::Literal) = [l.version]

get_versions(eg)

function solve_s(s)
    parsed, _ = gobble_packet(slurp_string(s))
    sum(+, get_versions(parsed))
end

function solve(filename)
    parsed, _ = gobble_packet(slurp_file(filename))
    sum(+, get_versions(parsed))
end

eval(l::Literal) = l.value

f_lookup = Dict(
    0 => sum,
    1 => values -> reduce(*, values),
    2 => values -> min(values...),
    3 => values -> max(values...),
    5 => values -> first(values) > last(values) ? 1 : 0,
    6 => values -> first(values) < last(values) ? 1 : 0,
    7 => values -> all(==(first(values)), values),
)

function eval(p::Operator)
    values = map(eval, p.subpackets)
    f = f_lookup[p.type_id]
    f(values)
end

function eval(s::AbstractString)
    exp, _ = gobble_packet(slurp_string(s))
    eval(exp)
end

eval_file(filename) = eval(readlines(filename)[1])

eval("C200B40A82")
eval("9C0141080250320F1802104A08") == 1

eval_file("16.txt")


text:
"""
DefaultState: track iterations and elapsed time

@author : Spencer Lyon <spencer.lyon@stern.nyu.edu>
@date : 2015-04-13 11:29:12

"""

mutable struct DefaultState{T} <: IterationState{T}
    n::Int              # number of iterations
    change::Float64     # change between previous and current state
    elapsed::Float64    # total time elapsed
    prev::T             # previous value
    prev_time::Float64  # previous absolute time
end

DefaultState(v) = DefaultState(0, Inf, 0.0, v, time())

function Base.show(io::IO, ds::DefaultState)
    n, d, t = ds.n, round(ds.change, 4), round(ds.elapsed, 4)
    m = "DefaultState: $n iterations, $d current change $t seconds elapsed"
    print(io, m)
end

num_iter(ds::DefaultState) = ds.n
LinearAlgebra.norm(istate::DefaultState) = abs(istate.change)

function display_iter(io::IO,  ds::DefaultState, prefix="")
    if ds.n == 0  # print banner
        @printf "%s%-13s%-15s%-17s\n" prefix "Iteration" "Distance" "Elapsed (seconds)"
        println(io, prefix, repeat("-", 45))
    else
        @printf "%s%-13i%-15.5e%-18.5f\n" prefix ds.n ds.change ds.elapsed
    end
end

update_prev!(istate::DefaultState{T}, v::T) where T <: AbstractArray =
    copy!(istate.prev, v)

function update_prev!(istate::DefaultState{T}, v::T) where T
    istate.prev = copy(v)
end

function update!(istate::DefaultState{T}, v::T; by::Function=default_by) where T
    istate.n += 1
    istate.change = by(istate.prev, v)
    update_prev!(istate, v)

    new_time = time()
    istate.elapsed += new_time - istate.prev_time
    istate.prev_time = new_time
    nothing
end


text:
"""
Brute force search for the optimal solution over all
index sets with m segemets. The costs are evaluated using least squares.

(Intended for verifying the dynamic programming algorithms)
"""
function brute_force_search(ℓ::AbstractTransitionCost{T}, V_N::QuadraticPolynomial{T}, m::Integer) where {T}
    cost_best = Inf
    I_best = Vector{Int64}(undef, m+1)
    Y_best = Vector{T}(undef, m+1)

    N = size(ℓ, 2)

    I = zeros(Int64, m+1)
    I[1] = 1
    I[end] = N

    for I_inner=IterTools.subsets(2:N-1, m-1)

        I[2:m] .= I_inner

        P = SymTridiagonal(zeros(T, m+1), zeros(T, m+1))
        q = zeros(T, m+1)
        r = 0

        # Add cost at right endpoint
        P[end,end] = V_N.a
        q[end]     = V_N.b
        r          = V_N.c

        # Form quadratic cost function Y'*P*Y + q'*Y + r
        # corresponding to the y-values in the vector Y
        for j=1:m
            P.dv[j] += ℓ[I[j], I[j+1]].P[1,1]
            P.dv[j+1] += ℓ[I[j], I[j+1]].P[2,2]
            P.ev[j] += ℓ[I[j], I[j+1]].P[1,2]
            q[j:j+1] .+= ℓ[I[j], I[j+1]].q
            r += ℓ[I[j], I[j+1]].r
        end

        # find the optimal Y-vector, and compute the correspinding error
        Yopt = -(P \ q) / 2
        cost = dot(Yopt, P*Yopt) + q' * Yopt + r

        if cost < cost_best
            cost_best = cost
            Y_best .= Yopt
            I_best .= I
        end
    end
    return I_best, Y_best, cost_best
end


text:
T = 3
n = 2
m = 1

A = [1.0 1.0; 0.0 1.0]
B = [0.0; 1.0]

Q = [1.0 0.0; 0.0 1.0]
R = [1.0]

x_init = [1.0; 0.0]

nz = T * n + (T - 1) * m + T * n

function unpack(z)
	x3 = z[1:n]
	u2 = z[n .+ (1:m)]
	x2 = z[n + m .+ (1:n)]
	u1 = z[n + m + n .+ (1:m)]
	x1 = z[n + m + n + m .+ (1:n)]
	y3 = z[n + m + n + m + n .+ (1:n)]
	y2 = z[n + m + n + m + n + n .+ (1:n)]
	y1 = z[n + m + n + m + n + n + n .+ (1:n)]

	return x3, u2, x2, u1, x1, y3, y2, y1
end

z0 = rand(nz)

function lagrangian(z, x_init)
	x3, u2, x2, u1, x1, y3, y2, y1 = unpack(z)

	L = 0.0

	# objective
	L += transpose(x3) * Q * x3
	L += (transpose(u2) * R * u2)[1]
	L += transpose(x2) * Q * x2
	L += (transpose(u1) * R * u1)[1]
	L += transpose(x1) * Q * x1

	# constraints
	L += transpose(y3) * (A * x2 + B * u2[1] - x3)
	L += transpose(y2) * (A * x1 + B * u1[1] - x2)
	L += transpose(y1) * (x_init - x1)

	return L
end

@variables z_sym[1:nz]
@variables x_sym[1:n]

L = lagrangian(z_sym, x_sym)
L = simplify.(L)

dL = Symbolics.gradient(L, z_sym)
ddL = Symbolics.hessian(L, z_sym)
ddLx = Symbolics.jacobian(dL, x_sym)

L_grad = eval(Symbolics.build_function(dL, z_sym, x_sym)[1])
L_hess = eval(Symbolics.build_function(ddL, z_sym, x_sym)[1])
Lx_hess = eval(Symbolics.build_function(ddLx, z_sym, x_sym)[1])

function solve(z0, x_init)
	z = copy(z0)

	∇L = L_grad(z, x_init)

	println()
	for i = 1:10
		r_norm = norm(∇L)
		println("iter ($i)")
		println("	grad norm: $(r_norm)")
		r_norm < 1.0e-8 && break

		∇²L = L_hess(z, x_init)

		Δ = ∇²L \ ∇L

		α = 1.0
		ẑ = z - α * Δ

		∇L̂ = L_grad(ẑ, x_init)

		iter = 1
		while norm(∇L̂) > norm(∇L)
			α *= 0.5
			ẑ = z - α * Δ
			∇L̂ = L_grad(ẑ, x_init)
			iter += 1

			if iter > 100
				break
			end
		end

		z = ẑ
		∇L = ∇L̂
	end

	return z, ∇L, L_hess(z, x_init)
end

z_sol, ∇L_sol, ∇²L_sol = solve(z0, x_init)
K_direct = (∇²L_sol \ Lx_hess(z_sol, x_init))[n + m + n .+ (1:m), :]

# LQR solution
K, P = tvlqr([A for t = 1:T-1], [B for t = 1:T-1], [Q for t = 1:T], [R[1,1] for t = 1:T-1])

@assert norm(K[1] - K_direct) < 1.0e-8

max_issues_repo_name: thowell/DirectMotionPlanning


text:
"""
    MatDeforModule

Module to support general operations for deformation material models.
"""
module MatDeforModule

__precompile__(true)

using FinEtools.FTypesModule: FInt, FFlt, FCplxFlt, FFltVec, FIntVec, FFltMat, FIntMat, FMat, FVec, FDataDict
import FinEtoolsDeforLinear.DeforModelRedModule: DeforModelRed3D, DeforModelRed2DStrain, DeforModelRed2DStress, DeforModelRed2DAxisymm, DeforModelRed1D
import FinEtools.MatModule: AbstractMat
using LinearAlgebra

_RotationMatrix = Union{Array{T, 2}, Transpose{T, Array{T, 2}}} where {T}

"""
    AbstractMatDefor

Abstract type that represents deformable materials.
"""
abstract type AbstractMatDefor <: AbstractMat end

include("genconv.jl")
include("genrot.jl")
include("tens4impl.jl")

end # module


text:
# using Revise
# using Distributed
# addprocs(4)

using RoME
# @everywhere using RoME

#  Do some plotting
using RoMEPlotting
# for drawing Bayes tree with images (see debug tricks below)
# using Cairo, Fontconfig
# using Gadfly

0

function driveHex(fgl, posecount::Int; steps::Int=5)
    # Drive around in a hexagon
    for i in (posecount-1):(posecount-1+steps)
        psym = Symbol("x$i")
        posecount += 1
        nsym = Symbol("x$(i+1)")
        addVariable!(fgl, nsym, Pose2)
        pp = Pose2Pose2(MvNormal([10.0;0;pi/3], Matrix(Diagonal([0.1;0.1;0.1].^2))))
        addFactor!(fgl, [psym;nsym], pp, autoinit=false )
    end

    return posecount
end

function offsetHexLeg(dfg::G, posecount::Int; direction=:right) where G <: AbstractDFG
    psym = Symbol("x$(posecount-1)")
    nsym = Symbol("x$(posecount)")
    posecount += 1
    addVariable!(dfg, nsym, Pose2)
    pp = nothing
    if direction == :right
        pp = Pose2Pose2(MvNormal([10.0;0;-pi/3], Matrix(Diagonal([0.1;0.1;0.1].^2))))
    elseif direction == :left
        pp = Pose2Pose2(MvNormal([10.0;0;pi/3], Matrix(Diagonal([0.1;0.1;0.1].^2))))
    end
    addFactor!(dfg, [psym; nsym], pp, autoinit=false )
    return posecount
end

## start with an empty factor graph object
fg = initfg()

# fg.solverParams
# fg.solverParams.isfixedlag = true
# fg.solverParams.qfl = 20
posecount = 0

# Add the first pose :x0
addVariable!(fg, :x0, Pose2)
posecount += 1

# Add at a fixed location PriorPose2 to pin :x0 to a starting location (10,10, pi/4)
addFactor!(fg, [:x0], PriorPose2( MvNormal([0.0; 0.0; 0.0],
                                           Matrix(Diagonal([0.1;0.1;0.05].^2))) ), autoinit=false )

# Add landmarks with Bearing range measurements
addVariable!(fg, :l1, Point2, labels=[:LANDMARK;])
p2br = Pose2Point2BearingRange(Normal(0,0.03),Normal(20.0,0.5))
addFactor!(fg, [:x0; :l1], p2br, autoinit=false )

## hex 1

posecount = driveHex(fg, posecount)

# Add landmarks with Bearing range measurements
p2br2 = Pose2Point2BearingRange(Normal(0,0.03),Normal(20.0,0.5))
addFactor!(fg, [Symbol("x$(posecount-1)"); :l1], p2br2, autoinit=false )

## hex 2

posecount = offsetHexLeg(fg, posecount, direction=:right)

# Add landmarks with Bearing range measurements
addVariable!(fg, :l2, Point2, labels=[:LANDMARK])
p2br = Pose2Point2BearingRange(Normal(0,0.03),Normal(20.0,0.5))
addFactor!(fg, [Symbol("x$(posecount-1)"); :l2], p2br, autoinit=false )

posecount = driveHex(fg, posecount, steps=5)

# Add landmarks with Bearing range measurements
p2br2 = Pose2Point2BearingRange(Normal(0,0.03),Normal(20.0,0.5))
addFactor!(fg, [Symbol("x$(posecount-1)"); :l2], p2br2, autoinit=false )

# writeGraphPdf(fg,engine="neato")

getSolverParams(fg).drawtree = true
getSolverParams(fg).showtree = true
getSolverParams(fg).downsolve = false
fg.solverParams.async = true

tree, smt, chi = solveTree!(fg, recordcliqs=ls(fg))

# hist = getCliqSolveHistory(tree, :x1)

Gadfly.set_default_plot_size(35cm,25cm)
drawPosesLandms(fg, meanmax=:max) |> PDF("/tmp/test.pdf");  @async run(`evince /tmp/test.pdf`)
# tree = wipeBuildNewTree!(fg)
# drawTree(tree, imgs=true)

#
# using Dates
#
# dd = [true;]
# mkpath("/tmp/caesar/anitree")
# @async begin
#   global dd
#   i =0
#   while dd[1]
#     i += 1
#     tt = split(string(now()), 'T')[end]
#     drawTree(tree, filepath="/tmp/caesar/anitree/bt_$i.pdf", show=false)
#     sleep(0.1)
#   end
# end
#
# # dd[1] = false

#


text:
using SparseArrays

const SHOW_EQUATIONS = Ref(false)
struct SystemStructure
    sys
    dxvar_offset::Int
    fullvars::Vector # [xvar; dxvars; algvars]
    varassoc::Vector{Int}
    graph::BipartiteGraph{Int}
    solvable_graph::BipartiteGraph{Int}
end
function SystemStructure(sys)
    sys = ModelingToolkit.flatten(sys)
    sys, dxvar_offset, fullvars, varassoc, graph, solvable_graph = init_graph(sys)
    SystemStructure(sys, dxvar_offset, fullvars, varassoc, graph, solvable_graph)
end

ModelingToolkit.equations(s::SystemStructure) = equations(s.sys)

function Base.show(io::IO, s::SystemStructure)
    @unpack fullvars, dxvar_offset, solvable_graph, graph = s
    algvar_offset = 2dxvar_offset
    print(io, "xvars: ")
    print(io, fullvars[1:dxvar_offset])
    print(io, "\ndxvars: ")
    print(io, fullvars[dxvar_offset+1:algvar_offset])
    print(io, "\nalgvars: ")
    print(io, fullvars[algvar_offset+1:end], '\n')

    if SHOW_EQUATIONS[]
        println(io, "Edges:")
        eqs = equations(s)
        for ev in 𝑠vertices(graph)
            print(io, "  $(eqs[ev])\n    -> ")
            vars = 𝑠neighbors(graph, ev)
            solvars = 𝑠neighbors(solvable_graph, ev)
            solvable = intersect(vars, solvars)
            notsolvable = setdiff(vars, solvars)

            print(io, join(string.(fullvars[notsolvable]), ", "))
            for ii in solvable
                print(io, ", ")
                var = fullvars[ii]
                Base.printstyled(io, string(fullvars[ii]), color=:cyan)
            end
            println(io)
        end
    end

    S = incidence_matrix(graph, Num(Sym{Real}(:×)))
    print(io, "Incidence matrix:")
    show(io, S)
end

# V-nodes `[x_1, x_2, x_3, ..., dx_1, dx_2, ..., y_1, y_2, ...]` where `x`s are
# differential variables and `y`s are algebraic variables.
function get_vnodes(sys)
    dxvars = []
    eqs = equations(sys)
    for (i, eq) in enumerate(eqs)
        if eq.lhs isa Symbolic
            # Make sure that the LHS is a first order derivative of a var.
            @assert operation(eq.lhs) isa Differential "The equation $eq is not in the form of `D(...) ~ ...`"
            @assert !(arguments(eq.lhs)[1] isa Differential) "The equation $eq is not first order"

            push!(dxvars, eq.lhs)
        end
    end

    xvars = (first ∘ var_from_nested_derivative).(dxvars)
    algvars  = setdiff(states(sys), xvars)
    return xvars, dxvars, algvars
end

function init_graph(sys)
    xvars, dxvars, algvars = get_vnodes(sys)
    dxvar_offset = length(xvars)
    algvar_offset = 2dxvar_offset

    fullvars = [xvars; dxvars; algvars]
    sys = reordersys(sys, dxvar_offset, fullvars)
    eqs = equations(sys)
    idxmap = Dict(fullvars .=> 1:length(fullvars))
    graph = BipartiteGraph(length(eqs), length(fullvars))
    solvable_graph = BipartiteGraph(length(eqs), length(fullvars))

    for (i, eq) in enumerate(eqs)
        if isdiffeq(eq)
            v = eq.lhs
            haskey(idxmap, v) && add_edge!(graph, i, idxmap[v])
        end
        # TODO: custom vars that handles D(x)
        vs = vars(eq.rhs)
        for v in vs
            haskey(idxmap, v) && add_edge!(graph, i, idxmap[v])
        end
    end

    varassoc = Int[(1:dxvar_offset) .+ dxvar_offset; zeros(Int, length(fullvars) - dxvar_offset)] # variable association list
    sys, dxvar_offset, fullvars, varassoc, graph, solvable_graph
end

function reordersys(sys, dxvar_offset, fullvars)
    eqs = equations(sys)
    neweqs = Vector{Equation}(undef, length(eqs))
    eqidxmap = Dict(@view(fullvars[dxvar_offset+1:2dxvar_offset]) .=> (1:dxvar_offset))
    varidxmap = Dict([@view(fullvars[1:dxvar_offset]); @view(fullvars[2dxvar_offset+1:end])] .=> (1:length(fullvars)-dxvar_offset))
    algidx = dxvar_offset
    for eq in eqs
        if isdiffeq(eq)
            neweqs[eqidxmap[eq.lhs]] = eq
        else
            neweqs[algidx+=1] = eq
        end
    end
    sts = states(sys)
    @set! sys.eqs = neweqs
    @set! sys.states = sts[map(s->varidxmap[s], sts)]
end