text:
const _UNARY_STUFF = [
    (-, b1),
    (abs, b1),
    (exp, b1),
    (log, b1),
    (log10, b1),
    (log2, b1),
    (sqrt, b1),
    (cbrt, b1),
    (sign, b1),
    (!, b_boolean),
    (inv, b1),
]

@testset "unary" begin
    for (f, block) in _UNARY_STUFF
        @testset "$f" begin
            _test_unary_op(f, block)
        end
    end

    @testset "inverse identities" begin
        @test -(-n1) === n1
        @test inv(inv(n1)) === n1
        @test !(!n_boolean) === n_boolean
    end
end

# FIXME: really need extended tests (especially for dot) where we use vector inputs too.
const _BINARY_FUNCTIONS = [+, -, *, /, ^, min, max, >, <, >=, <=, dot]

@testset "binary" begin
    for f in _BINARY_FUNCTIONS
        @testset "$f" begin
            # Test evaluations with different alignments.
            _test_binary_op(f)
        end
    end
end


text:
module Nystrom

using DocStringExtensions
using StaticArrays
using LinearAlgebra
using SparseArrays
using IterativeSolvers
using Base.Threads
using SpecialFunctions
using OrderedCollections

using WaveProp.Utils
using WaveProp.Geometry
using WaveProp.Integration
using WaveProp.Integration
using WaveProp.Mesh
using WaveProp.PhysicalProblem

import WaveProp.Geometry: geometric_dimension, ambient_dimension, domain, range
import WaveProp.Integration: singular_weights, qnodes, qweights, qnormals
import WaveProp.Mesh: etypes, elements

export
    Laplace,
    Helmholtz,
    Elastostatic,
    Stokes,
    Maxwell,
    SingleLayerKernel,
    DoubleLayerKernel,
    IntegralOperator,
    NystromMesh,
    Density,
    γ₀,
    γ₁,
    SingleLayerPotential,
    DoubleLayerPotential,
    SingleLayerOperator,
    DoubleLayerOperator,
    AdjointDoubleLayerOperator,
    HyperSingularOperator,
    GreensCorrection,
    singular_weights_dim,
    single_double_layer,
    isinside,
    NystromMesh,
    dom2dof,
    assemble

include("nystrommesh.jl")
include("greensformulae.jl")
include("lebedev.jl")
include("kernels.jl")
include("laplace.jl")
include("helmholtz.jl")
include("elastostatic.jl")
include("maxwell.jl")
include("density.jl")
include("potential.jl")
include("integraloperators.jl")
include("dim.jl")
include("ldim.jl")
include("assemble.jl")
include("martensen_kussmaul.jl")

end # module

ext: jl

max_stars_repo_name: WaveProp/WaveProp

max_issues_repo_name: WaveProp/WaveProp

max_forks_repo_name: WaveProp/WaveProp


text:
export buildPyNNModel_01_FromElements, buildPyNNModel_01_FromWeights, loadTfModelIntoFlux
export buildPose2OdoNN_01_FromElements, buildPose2OdoNN_01_FromWeights, loadPose2OdoNNModelIntoFlux

function buildPose2OdoNN_01_FromElements( W1::AbstractMatrix{<:Real}=zeros(4,8),
                                          b1::AbstractVector{<:Real}=zeros(8),
                                          W2::AbstractMatrix{<:Real}=zeros(8,48),
                                          b2::AbstractVector{<:Real}=zeros(8),
                                          W3::AbstractMatrix{<:Real}=zeros(2,8),
                                          b3::AbstractVector{<:Real}=zeros(2),
                                          pad::Bool=false)
  #
  # W1 = randn(Float32, 4,8)
  # b1 = randn(Float32,8)
  # outPad = pad ? x->[x;0] : x->x
  modjl = Chain(
    x -> (x*W1)' .+ b1 .|> relu,
    x -> reshape(x', 25,8,1),
    x -> maxpool(x, PoolDims(x, 4)),
    x -> reshape(x[:,:,1]',:),
    Dense(48,8,relu),
    Dense(8,2),
    x->[x;0]
  )

  modjl[5].W .= W2
  modjl[5].b .= b2

  modjl[6].W .= W3
  modjl[6].b .= b3

  return modjl
end

# As loaded from tensorflow get_weights
# Super specialized function
function buildPose2OdoNN_01_FromWeights(pywe; pad::Bool=false)
  buildPose2OdoNN_01_FromElements(pywe[1], pywe[2][:], pywe[3]', pywe[4][:], pywe[5]', pywe[6][:], pad)
end

# function buildPose2OdoNN_01_FromWeightsGPU(pywe)
#   buildPose2OdoNN_01_FromElements(collect(pywe[1]) |> gpu, collect(pywe[2][:]) |> gpu, collect(pywe[3]') |> gpu, collect(pywe[4][:]) |> gpu, collect(pywe[5]') |> gpu, collect(pywe[6][:]) |> gpu) |> gpu
# end

@deprecate buildPyNNModel_01_FromElements(args...;kw...) buildPose2OdoNN_01_FromElements(args...;kw...)
@deprecate buildPyNNModel_01_FromWeights(args...;kw...) buildPose2OdoNN_01_FromWeights(args...;kw...)

#


text:
using Test

@testset "Van-Wijngaarden" begin include("test_van_wijngaarden.jl") end
@testset "Weibull Count" begin include("test_weibull_count.jl") end
@testset "Distribution" begin include("test_distribution.jl") end


text:
mv("input.txt", "output.txt")
mv("docs", "mydocs")
mv("/input.txt", "/output.txt")
mv("/docs", "/mydocs")

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
#--------------------------------------------#
#    Environmental conditions rH and temperature    #
#--------------------------------------------#

using PlotlyJS
using DelimitedFiles
using LinearAlgebra
using Distributed
using DistributedArrays
using SparseArrays
using Krylov
# using ArrayFire
# using CUDAdrv, CUDAnative, CuArrays
Env=readdlm("inputs/environmental_year.csv",',')
Tem=Vector(undef,2919)	#hourly values from 3rd to 6th month (location Zagreb, year 2018... -9 h in final month)
r_H=Vector(undef,2919)	#hourly values from 3rd to 6th month (location Zagreb, year 2018... -9 h in final month)
j=1
for i=1:size(Env)[1]
	if (mod(Env[i,1],60)<=5)
		global j=j+1
	end
end
Env2=Array{Float64}(undef,j-1,3)
j=1
for i=1:size(Env)[1]
	if (mod(Env[i,1],60)<=5)
		Env2[j,1]=Env[i,1]
		Env2[j,2]=Env[i,2]
		Env2[j,3]=Env[i,3]
		global j=j+1
	end
end
j=1
for i=2:size(Env2)[1]
	if ~((Env2[i,1]-Env2[i-1,1])<=5)
		global j=j+1
	end
end
Env=Array{Float64}(undef,j-1,3)
j=1
for i=2:size(Env2)[1]
	if ~((Env2[i,1]-Env2[i-1,1])<=5)
		Env[j,1]=Env2[i,1]
		Env[j,2]=Env2[i,2]
		Env[j,3]=Env2[i,3]
		global j=j+1
	end
end
j=1
for i=2:size(Env)[1]
	if ~((Env[i,1]-Env[i-1,1])==60)
		global j=j+1
	end
end
Tem=Env[:,2]
r_H=Env[:,3]
function soil3(H)
	h=1.23-H
	if h<0.02
		7.8
	elseif h>=0.02&&h<0.05
		7.8-0.1/0.03*(h-0.02)
	elseif h>=0.05&&h<=0.1
		7.7+0.6/0.05*(h-0.05)
	elseif h>=0.1&&h<=0.2
		8.3-0.8/0.1*(h-0.1)
	elseif h>=0.2&&h<=0.3
		7.5-0.3/0.1*(h-0.2)
	elseif h>=0.3&&h<=0.5
		7.2-0.1/0.2*(h-0.3)
	elseif h>=0.5
		7.1
	end
end
function soil4(H)
	h=1.23-H
	if h<0.02
		13.1
	elseif h>=0.02&&h<0.05
		13.1-0.1/0.03*(h-0.02)
	elseif h>=0.05&&h<=0.1
		13.0+0.5/0.05*(h-0.05)
	elseif h>=0.1&&h<=0.2
		13.5-1.3/0.1*(h-0.1)
	elseif h>=0.2&&h<=0.3
		12.2-0.6/0.1*(h-0.2)
	elseif h>=0.3&&h<=0.5
		11.6-0.5/0.2*(h-0.3)
	elseif h>=0.5
		11.1-1.5/0.5*(h-0.5)
	end
end
function soil5(H)
	h=1.23-H
	if h<0.02
		20.5
	elseif h>=0.02&&h<0.05
		20.5-0.2/0.03*(h-0.02)
	elseif h>=0.05&&h<=0.1
		20.3+0.3/0.05*(h-0.05)
	elseif h>=0.1&&h<=0.2
		20.6-2/0.1*(h-0.1)
	elseif h>=0.2&&h<=0.3
		18.6-0.7/0.1*(h-0.2)
	elseif h>=0.3&&h<=0.5
		17.9-0.9/0.2*(h-0.3)
	elseif h>=0.5
		16.8-3/0.5*(h-0.5)
	end
end
function soil6(H)
	h=1.23-H
	if h<0.02
		24.8
	elseif h>=0.02&&h<0.05
		24.8-0.3/0.03*(h-0.02)
	elseif h>=0.05&&h<=0.1
		24.5+0.2/0.05*(h-0.05)
	elseif h>=0.1&&h<=0.2
		24.7-1.9/0.1*(h-0.1)
	elseif h>=0.2&&h<=0.3
		22.8-0.8/0.1*(h-0.2)
	elseif h>=0.3&&h<=0.5
		22-1.3/0.2*(h-0.3)
	elseif h>=0.5
		20.7-3.4/0.5*(h-0.5)
	end
end


text:
# exports
export create, Sphere, show, Shell

mutable struct Sphere <: Shape
    radius::Float64
end

mutable struct Shell <: Shape
    outer_radius::Float64
    inner_radius::Float64
end

function Base.show(io::IO, x::Sphere)
    println(io, "Sphere")
    println(io, "Radius: $(x.radius)")
end

function create(c::Sphere; resolution=nothing, rand_=0.0, type::Int64=1)
    radius = c.radius
    bounds = [-radius radius; -radius radius; -radius radius]
    x, v, y, vol, type_ = unpack(create(Cuboid(bounds), resolution=resolution, rand_=rand_))
    mask = vec(sum(x.^2, dims=1) .<= radius^2)
    mesh = x[:, mask]
    return Dict(
        :x => mesh,
        :v => zeros(size(mesh)),
        :y => copy(mesh),
        :volume => vol,
        :type => id*ones(Int64, length(vol))
    )
end

max_stars_repo_name: sajidmannan/PDMesh

max_issues_repo_name: sajidmannan/PDMesh

max_forks_repo_name: sajidmannan/PDMesh


text:
__precompile__(true)

module YAML

import Base: isempty, length, show

if VERSION < v"1.0.0-rc1.0"
    import Base: start, next, done
else
    import Base: iterate
end

import Codecs
using Compat
using Compat.Dates
using Compat.Printf

if VERSION < v"0.7.0-DEV.2915"
    isnumeric(c::Char) = isnumber(c)
end

if VERSION < v"0.7.0-DEV.3526"
    Base.parse(T::Type{Int}, s::AbstractString; base = 10) = parse(T, s, base)
end

include("scanner.jl")
include("parser.jl")
include("composer.jl")
include("constructor.jl")

const _constructor = Union{Nothing, Dict}

function load(ts::TokenStream, more_constructors::_constructor=nothing)
    events = EventStream(ts)
    node = compose(events)
    return construct_document(Constructor(more_constructors), node)
end

function load(input::IO, more_constructors::_constructor=nothing)
    load(TokenStream(input), more_constructors)
end

mutable struct YAMLDocIterator
    input::IO
    ts::TokenStream
    more_constructors::_constructor
    next_doc

    function YAMLDocIterator(input::IO, more_constructors::_constructor=nothing)
        it = new(input, TokenStream(input), more_constructors, nothing)
        it.next_doc = eof(it.input) ? nothing : load(it.ts, it.more_constructors)
        return it
    end
end

# Old iteration protocol:

start(it::YAMLDocIterator) = nothing

function next(it::YAMLDocIterator, state)
    doc = it.next_doc
    if eof(it.input)
        it.next_doc = nothing
    else
        reset!(it.ts)
        it.next_doc = load(it.ts, it.more_constructors)
    end
    return doc, nothing
end

done(it::YAMLDocIterator, state) = it.next_doc === nothing

# 0.7 iteration protocol:

iterate(it::YAMLDocIterator) = next(it, start(it))
iterate(it::YAMLDocIterator, s) = done(it, s) ? nothing : next(it, s)

function load_all(input::IO, more_constructors::_constructor=nothing)
    YAMLDocIterator(input, more_constructors)
end

function load(input::AbstractString, more_constructors::_constructor=nothing)
    load(IOBuffer(input), more_constructors)
end

function load_all(input::AbstractString, more_constructors::_constructor=nothing)
    load_all(IOBuffer(input), more_constructors)
end

function load_file(filename::AbstractString, more_constructors::_constructor=nothing)
    open(filename, "r") do input
        load(input, more_constructors)
    end
end

function load_all_file(filename::AbstractString, more_constructors::_constructor=nothing)
    open(filename, "r") do input
        load_all(input, more_constructors)
    end
end

end  # module


text:
using LocalSearchSolvers
using CBLS
using ConstraintModels
using ConstraintDomains
using CompositionalNetworks
using Constraints
using Documenter

@info "Makeing documentation..."
makedocs(;
    modules=[LocalSearchSolvers, CBLS, ConstraintModels,
        ConstraintDomains, CompositionalNetworks, Constraints],
    expandfirst = ["variables.md", "constraints.md", "objectives.md", "solving.md"],
    authors="Jean-François Baffier",
    repo="https://github.com/JuliaConstraints/LocalSearchSolvers.jl/blob/{commit}{path}#L{line}",
    sitename="LocalSearchSolvers.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", nothing) == "true",
        canonical="https://JuliaConstraints.github.io/LocalSearchSolvers.jl",
        assets = ["assets/favicon.ico"; "assets/github_buttons.js"; "assets/custom.css"],
    ),
    pages=[
        "Home" => "index.md",
        "Manual" => [
            "Quick Start Guide" => "quickstart.md",
            "Variables" => "variables.md",
            "Constraints" => "constraints.md",
            "Objectives" => "objectives.md",
            "Solving" => "solving.md",
        ],
        "Examples" => [
            "Sudoku" => "sudoku.md",
            "Golomb ruler" => "golomb.md",
            "Mincut" => "mincut.md",
        ],
        "Related" => [
            "ConstraintDomains.jl" => "domain.md",
            "Constraints.jl" => "d_constraint.md",
            "CompositionalNetworks.jl" => "icn.md",
            "CBLS.jl" => "cbls.md",
            "ConstraintModels.jl" => "models.md",
        ],
        "Library" => [
            "Public" => "public.md",
            "Internals" => "internals.md",
        ],
        "Constributing" => "contributing.md",
    ],
)

deploydocs(;
    repo="github.com/JuliaConstraints/LocalSearchSolvers.jl.git",
    devbranch="main",
)


text:
# (GEM) iJR904 (Reed et al., 2003) (download link: https://darwin.di.uminho.pt/models)
# TODO: change to alternative link!!!
# (GEM) Alternative link http://bigg.ucsd.edu/static/models/iJR904.mat
module Chemostat_Heerden2013

    import Chemostat
    const Ch = Chemostat

    import UtilsJL
    const UJL = UtilsJL
    UJL.gen_top_proj(@__MODULE__)

    include("Utils/Utils.jl")
    include("BegData/BegData.jl")
    include("HeerdenData/HeerdenData.jl")
    include("iJR904/iJR904.jl")
    # include("EColiCore/EColiCore.jl")

    function __init__()
        UJL.create_proj_dirs(@__MODULE__)
    end

end # module


text:
# License for this file: MIT (expat)
# Copyright 2017-2018, DLR Institute of System Dynamics and Control
#
# This file is part of module
#   ModiaMath.Result (ModiaMath/Result/_module.jl)
#

"""
    ModiaMath.plot(result, names; heading="", grid=true, xAxis= :time,
                   figure=1, prefix="", reuse=false, maxLegend=10)

Plot time series of the result defined by the names keys (Symbol or String).
The keys (and their units, if available in the result) are automatically used as legend.
Units can be either added by using package `Unitful` if result is just a dictionary, or it
can be added by using package `ModiaMath.Result`, where units are defined as elements
of the variable definition.

# Arguments
Argument `result` maybe one of the following:
- A dictionary `Dict{AbstractString,Any}`
- An instance of struct [`ModiaMath.Result`](@ref)
- An object for which function [`ModiaMath.resultTimeSeries`](@ref) is defined.

Argument `names` defines the diagrams to be drawn and the time series to be included in the respective diagram:
- If names is a **Symbol** or **String**, generate one diagram with one time series.
- If names is a **Tuple** of Symbols/Strings, generate one diagram with the time series of the given keys
- If names is a **Vector** or **Matrix** of Symbols/Strings/Tuples, generate a vector or matrix of diagrams.

Remaining arguments:
- `heading::AbstractString`: Optional heading above the diagram.
- `grid::Bool`: Optional grid.
- `xAxis`: Name of x-axis (Symbol or AbstractString).
- `figure::Int`: Integer identifier of the window in which the diagrams shall be drawn.
- `prefix::AbstractString`: String that is appended in front of every legend label (useful especially if reuse=true)
- `reuse::Bool`: If figure already exists and reuse=false, clear the figure before adding the plot.
- `maxLegend::Int`: If the number of legend entries in one plot command > maxLegend, the legend is suppressed.
  All curves have still their names as labels. The curves can be inspected by their names by clicking in the
  toolbar of the plot on button `Edit axis, curve ..` and then on `Curves`.

# Examples
```julia
import ModiaMath
using Unitful

t = linspace(0.0, 10.0, 100)
result = Dict{AbstractString,Any}(
            "time" => t*u"s", "phi1" => sin.(t)u"rad"  , "phi2" => 0.5*sin.(t),
                              "w1"   => cos.(t)u"rad/s", "w2"   => 0.6*cos.(t))

# 1 signal in one diagram
#   (legend = "phi1 [rad]")
ModiaMath.plot(result, :phi1)

# 3 signals in one diagram
ModiaMath.plot(result, ("phi1", :phi2, :w1), figure=2)

# 3 diagrams in form of a vector (every diagram has one signal)
ModiaMath.plot(result, [:phi1, :phi2, :w1], figure=3)

# 4 diagrams in form of a matrix (every diagram has one signal)
ModiaMath.plot(result, ["phi1" "phi2";
                        "w1"   "w2"   ], figure=4)

# 2 diagrams in form of a vector
ModiaMath.plot(result, [ (:phi1,:phi2), (:w1) ], figure=5)

# 4 diagrams in form of a matrix
ModiaMath.plot(result, [ (:phi1,)           (:phi2,:w1);
                         (:phi1,:phi2,:w1)  (:w2,)     ],figure=6)

# Plot w1=f(phi1) in one diagram
ModiaMath.plot(result, :w1, xAxis=:phi1, figure=7)

# Append signal of the next simulation run to figure=1
# (legend = "Sim 2: phi1 [rad]")
result[:phi1] = 0.5*result[:phi1]
ModiaMath.plot(result, :phi1, prefix="Sim 2: ", reuse=true)
```

The 5th example above (2 diagrams in form of a vector) give the following plot:

![Figure 5](../../resources/images/plot_figure5.svg)
"""
plot(result, names::Symbol; heading::AbstractString="", grid::Bool=true, xAxis="time", figure::Int=1, prefix::AbstractString="", reuse::Bool=false, maxLegend::Integer=10) =
    plot(result, string(names), heading=heading, grid=grid, xAxis=string(xAxis), figure=figure, prefix=prefix, reuse=reuse, maxLegend=maxLegend)

plot(result, names::AbstractVector; heading::AbstractString="", grid::Bool=true, xAxis="time", figure::Int=1, prefix::AbstractString="", reuse::Bool=false, maxLegend::Integer=10) =
    plot(result, reshape(names, length(names), 1); heading=heading, grid=grid, xAxis=string(xAxis), figure=figure, prefix=prefix, reuse=reuse, maxLegend=maxLegend)


text:
## From https://github.com/JuliaGNI/SimpleSolvers.jl/blob/master/src/linear/lu_solver.jl

mutable struct LUSolver{T}
    n::Int
    A::Matrix{T}
    b::Vector{T}
    x::Vector{T}
    pivots::Vector{Int}
    perms::Vector{Int}
    info::Int

    LUSolver{T}(n) where {T} = new(n, zeros(T, n, n), zeros(T, n), zeros(T, n), zeros(Int, n), zeros(Int, n), 0)
end

function LUSolver(A::Matrix{T}) where {T}
    n = LinearAlgebra.checksquare(A)
    lu = LUSolver{eltype(A)}(n)
    lu.A .= A
    lu
end

function LUSolver(A::Matrix{T}, b::Vector{T}) where {T}
    n = LinearAlgebra.checksquare(A)
    @assert n == length(b)
    lu = LUSolver{eltype(A)}(n)
    lu.A .= A
    lu.b .= b
    lu
end

function simplelu_factorize!(lu::LUSolver{T}, pivot=true) where {T}
    A = lu.A

    begin
        @inbounds for i in eachindex(lu.perms)
            lu.perms[i] = i
        end

        @inbounds for k = 1:lu.n
            # find index max
            kp = k
            if pivot
                amax = real(zero(T))
                for i = k:lu.n
                    absi = abs(A[i,k])
                    if absi > amax
                        kp = i
                        amax = absi
                    end
                end
            end
            lu.pivots[k] = kp
            lu.perms[k], lu.perms[kp] = lu.perms[kp], lu.perms[k]

            if A[kp,k] != 0
                if k != kp
                    # Interchange
                    for i = 1:lu.n
                        tmp = A[k,i]
                        A[k,i] = A[kp,i]
                        A[kp,i] = tmp
                    end
                end
                # Scale first column
                Akkinv = inv(A[k,k])
                for i = k+1:lu.n
                    A[i,k] *= Akkinv
                end
            elseif lu.info == 0
                lu.info = k
            end
            # Update the rest
            for j = k+1:lu.n
                for i = k+1:lu.n
                    A[i,j] -= A[i,k]*A[k,j]
                end
            end
        end

        lu.info
    end
end

function simplelu_solve!(lu::LUSolver{T}) where {T}
    @inbounds for i = 1:lu.n
        lu.x[i] = lu.b[lu.perms[i]]
    end

    @inbounds for i = 2:lu.n
        s = zero(T)
        for j = 1:i-1
            s += lu.A[i,j] * lu.x[j]
        end
        lu.x[i] -= s
    end

    lu.x[lu.n] /= lu.A[lu.n,lu.n]
    @inbounds for i = lu.n-1:-1:1
        s = zero(T)
        for j = i+1:lu.n
            s += lu.A[i,j] * lu.x[j]
        end
        lu.x[i] -= s
        lu.x[i] /= lu.A[i,i]
    end

    copyto!(lu.b,lu.x)

    lu.x
end

### Wrapper

struct SimpleLUFactorization <: AbstractFactorization
    pivot::Bool
    SimpleLUFactorization(pivot=true) = new(pivot)
end

function SciMLBase.solve(cache::LinearCache, alg::SimpleLUFactorization; kwargs...)
    if cache.isfresh
        cache.cacheval.A = cache.A
        simplelu_factorize!(cache.cacheval, alg.pivot)
    end
    cache.cacheval.b = cache.b
    cache.cacheval.x = cache.u
    y = simplelu_solve!(cache.cacheval)
    SciMLBase.build_linear_solution(alg,y,nothing,cache)
end

init_cacheval(alg::SimpleLUFactorization, A, b, u, Pl, Pr, maxiters, abstol, reltol, verbose) = LUSolver(convert(AbstractMatrix,A))


text:
export  =>ˡ,  =>ʳ, =>ˡʳ, HigherIntervention, hi

struct HigherIntervention{T} <: AbstractIntervention
  f::T
end

"""
`hi(f::RV{RV})`

Higher order intervention.

If `f` is a variable over interventions then `hi` allows you to intervene
on the variable it points to.

```julia
using OmegaCore, Distributions
μ = 1 ~ Normal(0, 1)
σ = 2 ~ Uniform(1, 3)
y = 3 ~ Normal(μ, σ)
choice(ω) = ifelse((4 ~ Bernoulli(0.5))(ω), μ, σ)
int_dist(ω) = Intervention(choice(ω) => 5.0)

joint = @joint(y_, μ, σ, choice)
joint_ = joint |ᵈ hi(int_dist)
randsample(joint_)
```

"""
hi(f) = HigherIntervention(f)

(hi::HigherIntervention)(ω) = hi.f(ω)

"Random variable over interventions"
struct LR{M, X, V}
  x::X
  v::V
end

LR{M}(x::X, v::V) where M where {X, V} = LR{M, X, V}(x, v)

x =>ˡ v = HigherIntervention(LR{:L}(x, v))

"`x ²¦ v`: `ω -> x(ω) => v`"
x =>ʳ v = HigherIntervention(LR{:R}(x, v))

"`x ²¦² v`: `ω -> x(ω) => v(ω)`"
x =>ˡʳ v = HigherIntervention(LR{:LR}(x, v))

(i::LR{:L})(ω) = autointervention(i.x(ω) => i.v)
(i::LR{:R})(ω) = autointervention(i.x => i.v(ω))
(i::LR{:LR})(ω) = autointervention(i.x(ω) => i.v(ω))

max_stars_repo_name: zenna/expect

max_issues_repo_name: zenna/expect

max_forks_repo_name: zenna/expect


text:
module Uploaders

using Genie, Stipple, StippleUI, StippleUI.API
import Genie.Renderer.Html: HTMLString, normal_element, register_normal_element

export uploader

register_normal_element("q__uploader", context = @__MODULE__)

function uploader(args...; kwargs...)
  q__uploader(args...; kw(kwargs)...)
end

end


text:
import KernelDensity: UnivariateKDE, default_bandwidth, kde
import Distributions: Distribution

"""
    bin(x::AbstractUncertainIndexValueDataset, binning::BinnedWeightedResampling{UncertainScalarKDE}) -> UncertainIndexValueDataset
    bin(x::AbstractUncertainIndexValueDataset, binning::BinnedWeightedResampling{UncertainScalarPopulation}) -> UncertainIndexValueDataset

Resample every element of `x` a number of times. After resampling, distribute the
values according to their indices, into the `N` bins given by the `N-1`-element
grid defined by `binning.left_bin_edges`. In total, `length(x)*binning.n` draws
are distributed among the bins. The precise number of times `x[i]` is resampled is
given by `binning.weights[i]` (probability weights are always normalised to 1).

## Returns

Returns an `UncertainIndexValueDataset`. Indices are assumed to be uniformly distributed within each
bin, and are represented as `CertainValue`s at the bin centers. Values of the dataset have different
representations depending on what `binning` is:

- If `binning isa BinnedWeightedResampling{UncertainScalarKDE}`, then values in each bin are
    represented by a kernel density estimate to the distribution of the resampled values whose
    resampled indices fall in that bin.
- If `binning isa BinnedWeightedResampling{UncertainScalarPopulation}`, then values in each bin are
    represented by equiprobable populations consisting of the resampled values whose resampled
    indices fall in the bins.
"""
function bin(x::AbstractUncertainIndexValueDataset, binning::BinnedWeightedResampling); end

"""
    bin(x::AbstractUncertainIndexValueDataset,
        binning::BinnedWeightedResampling{RawValues}) -> Tuple(Vector, Vector{Vector})

Resample every element of `x` a number of times. After resampling, distribute the
values according to their indices, into the `N` bins given by the `N-1`-element
grid defined by `binning.left_bin_edges`. In total, `length(x)*binning.n` draws
are distributed among the bins. The precise number of times `x[i]` is resampled is
given by the `binning.weights[i]` (probability weights are always normalised to 1).

## Returns

Return a tuple containing the `N` different bin centers and a `N`-length vector of
resampled values whose resampled indices fall in the `N` different bins.

## Example

```julia
using Plots, UncertainData
# Some example data with unevenly spaced time indices
function ar1(n::Int, x0 = 0.5, p = 0.3)
    vals = zeros(n)
    [vals[i] = vals[i - 1]*p + rand()*0.5 for i = 2:n]
    return vals
end

npts = 50
time, vals = sort(rand(1:1000, npts)), ar1(npts)

# Add uncertainties to indices and values, and represent as
# UncertainIndexValueDataset
utime = [UncertainValue(Normal, t, 5) for t in time]
uvals = [UncertainValue(Normal, v, 0.03) for v in vals]

udata = UncertainIndexValueDataset(utime, uvals)

# Bin data into fall in 25 time step wide time bins ranging
# from time indices 100 to 900 and return a vector of raw
# values for each bin. Do this by resampling each uncertain
# data point on average 10000 times and distributing those
# draws among the bins.
time_grid = 100:40:900
n_draws = 5000
# Let odd-indexed values be three times as likely to be
# sampled compared to even-indexed values.
wts = Weights([i % 2 == 0 ? 1 : 3 for i = 1:length(udata)])
binning = BinnedWeightedResampling(RawValues, time_grid, wts, n_draws)

bin_centers, bin_draws = bin(udata, binning);
```
"""
function bin(x::AbstractUncertainIndexValueDataset, binning::BinnedWeightedResampling{RawValues};
        nan_threshold = 0.0)

    # Determine how many times each element of `x` should be resampled
    # based on the provided weights.
    wts = binning.weights ./ sum(binning.weights)

    n_total_draws = binning.n * length(x)
    Ns = ceil.(Int, n_total_draws .* binning.weights)

    # Separately convert indices and values to weighted populations
    pop_inds = UncertainValue(x.indices.indices, wts)
    pop_vals = UncertainValue(x.values.values, wts)

    bin_edges = binning.left_bin_edges
    n_bins = length(bin_edges) - 1

    # Used to compute the index of the bin into which a draw belongs
    mini = minimum(binning.left_bin_edges)
    s = step(binning.left_bin_edges)

    # Empty vectors that will contain draws.
    binvecs = [Vector{Float64}(undef, 0) for i = 1:n_bins]
    #[sizehint!(bv, n_total_draws) for bv in binvecs]

    for (i, N) in enumerate(Ns)
        # Pre-allocate some arrays into which we resample the values of the
        # index and value populations.
        idxs = fill(NaN, N)
        vals = fill(NaN, N)

        perminds = zeros(Int, N)
        sorted_idxs = fill(NaN, N)
        sorted_vals = fill(NaN, N)

        # Resample the i-th index and i-th value
        resample!(idxs, pop_inds[i])
        resample!(vals, pop_vals[i])

        # Get the vector that sorts the index vector, and use that to
        # sort the draws.
        sortperm!(perminds, idxs)
        sorted_idxs .= idxs[perminds]
        sorted_vals .= vals[perminds]

        # The vectors above are sorted sorted, so this can be done faster
        for i in 1:n_bins
            inbin = findall(bin_edges[i] .<= sorted_idxs .<= bin_edges[i+1])
            if length(inbin) > nan_threshold
                append!(binvecs[i], sorted_vals[inbin])
            end
        end
    end

    bin_centers = bin_edges[1:end-1] .+ step(bin_edges)/2

    return bin_centers, binvecs
end

function bin(x::AbstractUncertainIndexValueDataset, binning::BinnedWeightedResampling{UncertainScalarKDE};
        nan_threshold = 0.0)

    bin_edges = binning.left_bin_edges
    n = binning.n
    wts = binning.weights
    binning = BinnedWeightedResampling(RawValues, bin_edges, wts, n)

    bin_centers, binvecs = bin(x, binning)

    # Estimate distributions in each bin by kernel density estimation
    n_bins = length(binning.left_bin_edges) - 1
    estimated_value_dists = Vector{Union{CertainValue, UncertainScalarKDE}}(undef, n_bins)

    for i in 1:n_bins
        if length(binvecs[i]) > nan_threshold
            estimated_value_dists[i] = UncertainValue(binvecs[i])
        else
            estimated_value_dists[i] = UncertainValue(NaN)
        end
    end

    new_inds = UncertainIndexDataset(UncertainValue.(bin_edges[1:end-1] .+ step(bin_edges)/2))
    new_vals = UncertainValueDataset(estimated_value_dists)

    UncertainIndexValueDataset(new_inds, new_vals)
end

function bin(x::AbstractUncertainIndexValueDataset, binning::BinnedWeightedResampling{UncertainScalarPopulation};
        nan_threshold = 0.0)

    bin_edges = binning.left_bin_edges
    n = binning.n
    wts = binning.weights
    binning = BinnedWeightedResampling(RawValues, bin_edges, wts, n)

    bin_centers, binvecs = bin(x, binning)

    # Estimate distributions in each bin by kernel density estimation
    n_bins = length(binning.left_bin_edges) - 1
    estimated_value_dists = Vector{Union{CertainValue, UncertainScalarPopulation}}(undef, n_bins)

    for i in 1:n_bins
        if length(binvecs[i]) > nan_threshold
            L = length(binvecs[i])
            estimated_value_dists[i] = UncertainValue(binvecs[i], repeat([1 / L], L))
        else
            estimated_value_dists[i] = UncertainValue(NaN)
        end
    end

    new_inds = UncertainIndexDataset(UncertainValue.(bin_edges[1:end-1] .+ step(bin_edges)/2))
    new_vals = UncertainValueDataset(estimated_value_dists)

    UncertainIndexValueDataset(new_inds, new_vals)
end

import KernelDensity: UnivariateKDE, default_bandwidth, kde
import Distributions: Distribution


text:
"""
长方形有关的功能
"""

struct AbstractRectangle{T} <: AbstractPolygon{T}
    vertex :: Array{Point2D, 1}
    edges :: Array{Segment, 1}
    center :: Point2D
    width :: Float64
    height :: Float64
end

"""
构造一个平行坐标轴的长方形
"""
function Rectangle(center::Point2D, width, height)
    verts = [
        Point2D(center.x - width/2, center.y + height/2),
        Point2D(center.x - width/2, center.y - height/2),
        Point2D(center.x + width/2, center.y - height/2),
        Point2D(center.x + width/2, center.y + height/2)
    ]
    edges = [
        Segment(verts[1], verts[2]),
        Segment(verts[2], verts[3]),
        Segment(verts[3], verts[4]),
        Segment(verts[4], verts[1])
    ]
    return AbstractRectangle{:AXISRECT}(verts, edges, center, width, height)
end

function Square(center::Point2D, width)
    verts = [
        Point2D(center.x - width/2, center.y + width/2),
        Point2D(center.x - width/2, center.y - width/2),
        Point2D(center.x + width/2, center.y - width/2),
        Point2D(center.x + width/2, center.y + width/2)
    ]
    edges = [
        Segment(verts[1], verts[2]),
        Segment(verts[2], verts[3]),
        Segment(verts[3], verts[4]),
        Segment(verts[4], verts[1])
    ]
    return AbstractRectangle{:AXISSQUARE}(verts, edges, center, width, width)
end

"""
长方型的面积
"""
function area(rect::AbstractRectangle{N}) where N
    return rect.height * rect.width
end

#display函数在中间加::MIME"text/plain"
function Base.show(io::IO, tri::AbstractRectangle{N}) where N
    name = Dict(
        :AXISRECT => "Rectangle",
        :AXISSQUARE => "Square"
    )
    @printf(io, "%s\n center: %s\n width: %.6f\n height: %.6f",
    name[N], tri.center, tri.width, tri.height)
end


text:
module TestEcology
using Test

# Looking at relations to historical measures, updated with similarity
# and partitioning
using Diversity
using Diversity.Ecology
using Diversity.ShortNames

numspecies = 100;
numcommunities = 8;
communities = rand(numspecies, numcommunities);
communities /= sum(communities);
eco = Metacommunity(communities)
weights = rand(numspecies);
weights /= sum(weights);
Z1 = ones(typeof(weights[1]), (length(weights), length(weights)));

@testset "Standard ecological diversities" begin
    @test richness(communities)[!,:diversity] ≈ subdiv(ᾱ(eco), 0)[!,:diversity]
    @test richness(Metacommunity(communities))[!,:diversity] ==
        richness(communities)[!,:diversity]
    @test_throws ErrorException richness(Metacommunity(communities,
        GeneralTypes(rand(numspecies,numspecies))))

    @test shannon(communities)[!,:diversity] ≈ log.(subdiv(ᾱ(eco), 1)[!,:diversity])
    @test shannon(Metacommunity(communities))[!,:diversity] ==
        shannon(communities)[!,:diversity]
    @test_throws ErrorException shannon(Metacommunity(communities,
        GeneralTypes(rand(numspecies,numspecies))))

    @test simpson(communities)[!,:diversity] ≈ subdiv(ᾱ(eco), 2)[!,:diversity] .^ -1
    @test simpson(Metacommunity(communities))[!,:diversity] ==
        simpson(communities)[!,:diversity]
    @test_throws ErrorException simpson(Metacommunity(communities,
        GeneralTypes(rand(numspecies,numspecies))))

    @test jaccard([1 0 0; 0 1 1]'/3)[!,:diversity] .+ 1.0 ≈ [1.0]
    @test jaccard(Metacommunity([1 0 1; 0 1 1]'/4))[!,:diversity] ==
        jaccard([1 0 1; 0 1 1]'/4)[!,:diversity]
    @test_throws ErrorException jaccard(Metacommunity([1 0 1; 0 1 1]'/4,
                                                      GeneralTypes(rand(3, 3))))
    @test jaccard([1 0 1; 0 1 1]'/4)[!,:diversity] ≈ [1.0 / 3.0]
    @test_throws ErrorException jaccard([1 1 0; 0 1 1; 1 1 1]/7)
    @test_throws ErrorException jaccard(Metacommunity([1 1 0; 0 1 1; 1 1 1]/7))
end

@testset "Generalised ecological diversities" begin
    @test generalisedrichness(metacommunityDiversity,
                              communities, Z1)[!,:diversity] ≈ [1]
    @test generalisedrichness(subcommunityDiversity,
                              communities)[!,:diversity] ≈
          richness(communities)[!,:diversity]

    @test_throws ErrorException generalisedrichness(individualDiversity,
                                                    communities, Z1)

    @test generalisedshannon(metacommunityDiversity,
                             communities, Z1)[!,:diversity] .+ 1.0 ≈ [1.0]
    @test generalisedshannon(subcommunityDiversity,
                             communities)[!,:diversity] ≈
          shannon(communities)[!,:diversity]
    @test_throws ErrorException generalisedshannon(individualDiversity,
                                                   communities, Z1)

    @test generalisedsimpson(metacommunityDiversity,
                             communities, Z1)[!,:diversity] ≈ [1.0]
    @test generalisedsimpson(subcommunityDiversity,
                             communities)[!,:diversity] ≈
          simpson(communities)[!,:diversity]
    @test_throws ErrorException generalisedsimpson(individualDiversity,
                                                   communities, Z1)

    @test generalisedjaccard([1 0 1; 0 1 1]'/4, [0, Inf])[!,:diversity] ≈ [1.0/3.0, 1.0]
    @test generalisedjaccard([1 1 1; 1 1 1]'/6, [0, 1])[!,:diversity] ≈ [1.0, 1.0]
end

end


text:
struct GitHub <: DataRepo
end

base_url(::GitHub) = "https://github.com/"

git_repo_page_url(page) = base_url(GitHub()) * getattr(first(eachmatch(sel"strong[itemprop=\"name\"] a", page.root)), "href")

function description(::GitHub, mainpage)
    # Just load the readme -- it is all that we can do as this level of generic
    desc = get_docfile(GitHub(), mainpage, "README", 12)

    license_text = get_docfile(GitHub(), mainpage, "LICENSE", 4)
    if license_text == ""
        desc = "License: Unknown\n\n" * desc
    else
        desc *= "\n\nLICENSE\n --------\n" * license_text
    end
    desc
end

function get_docfile(::GitHub, page, docname, max_lines=typemax(Int))
    function inner(page)
        nodes = eachmatch(Selector(".content span a :contains($(docname))"), page.root)
        if length(nodes)>0
            node = first(nodes)
            url = "https://rawgit.com" * getattr(node.parent, "href")
            url = replace(url, "blob/" => "")
            text = getpage_raw(url) # It is plain-text/markdown probs.
            lines = split(text, "\n")
            if length(lines) > max_lines
                text = join(lines[1:max_lines], "\n")
                text *= "...\n (Read more at $(url))"
            end
            text
        else
            ""
        end
    end

    doc = inner(page)
    if doc==""
        doc = inner(getpage(git_repo_page_url(page))) # Try the base dir of the repository
    end
    doc
end

function get_cdn_url_converter(mainpage)
    commit = match(r"\b(?<=(commit\/))[0-9a-f]{5,40}\b", string(mainpage.root))
    function(urlsub)
        ret = "https://cdn.rawgit.com"*urlsub
        if commit === nothing
            @warn("Not able to retrieve commit hash. Switching to master.")
            rep_string = "master"
        else
            rep_string = commit.match
        end
        replace(ret, "blob/master" => rep_string)
    end
end

function get_urls(repo::GitHub, page, cdn_url_converter=get_cdn_url_converter(page))
    links = eachmatch(sel".content span a", page.root)
    urls = Any[]
    for link in links
        urlsub =  getattr(link, "href")
        github_url = base_url(repo) * urlsub
        linkclass = getattr(link.parent.parent.parent[1][1], "class")
        if occursin("directory", linkclass) # is a dirctory
            push!(urls, get_urls(repo, getpage(github_url), cdn_url_converter)) # making a list of lists
        else # it is a file
            push!(urls, cdn_url_converter(urlsub))
        end

    end
    urls
end

function data_fullname(::GitHub, mainpage)
    text_only(last(eachmatch(sel"h1", mainpage.root)))
end


text:
println("Dv) Test of the PairAnnihilation1Photon module with ASF from an internally generated initial- and final-state multiplet.")

@warn("\n\n !!! This example does not work properly at present !!! \n\n")
#
wa = Atomic.Computation("xx",  Nuclear.Model(26.), JAC.Radial.Grid("grid: exponential"), false, [JAC.NoProperty],
     Configuration[],
     Configuration[Configuration("1s^2")],
     Configuration[],
     Configuration[Configuration("1s") ],
     Einstein.Settings(),  Hfs.Settings(),  IsotopeShift.Settings(),   PlasmaShift.Settings(), LandeZeeman.Settings(),
     JAC.PairA1P,  PairAnnihilation1Photon.Settings([E1, M1], [JAC.UseCoulomb, JAC.UseBabushkin], [0.], true, false, Tuple{Int64,Int64}[])  )

wb = perform(wa)

#     [Configuration("[Ne] 3s^2 3p^5"), Configuration("[Ne] 3s 3p^6"), Configuration("[Ne] 3s^2 3p^4 3d")],


text:
using Plots

# Model
include_model("flamingo")
model = Flamingo{Discrete, FixedTime}(n, m, d,
			  g_world, 0.5,
			  l_torso, d_torso, m_torso, J_torso,
			  l_thigh, d_thigh, m_thigh, J_thigh,
			  l_calf, d_calf, m_calf, J_calf,
			  l_foot, d_foot, m_foot, J_foot,
			  l_thigh, d_thigh, m_thigh, J_thigh,
			  l_calf, d_calf, m_calf, J_calf,
			  l_foot, d_foot, m_foot, J_foot,
			  qL, qU,
			  uL, uU,
			  nq,
			  nu,
			  nc,
			  nf,
			  nb,
			  ns,
			  idx_u,
			  idx_λ,
			  idx_b,
			  idx_ψ,
			  idx_η,
			  idx_s,
			  0.0 * joint_friction)# model = free_time_model(model)

# Visualize
# - Pkg.add any external deps from visualize.jl
include(joinpath(pwd(), "models/visualize.jl"))
vis = Visualizer()
render(vis)
open(vis)
"""
Useful methods

joint-space inertia matrix:
	M_func(model, q)

position (px, pz) of foot 1 heel:
	kinematics_2(model,
		get_q⁺(x), body = :calf_1, mode = :ee)

position (px, pz) of foot 2 heel:
	kinematics_2(model,
		get_q⁺(x), body = :calf_2, mode = :ee)

position (px, pz) of foot 1 toe:
		kinematics_3(model,
	        get_q⁺(x), body = :foot_1, mode = :ee)

position (px, pz) of foot 2 toe:
		kinematics_3(model,
	        get_q⁺(x), body = :foot_2, mode = :ee)

jacobian to foot 1 heel
	jacobian_2(model, get_q⁺(x), body = :calf_1, mode = :ee)

jacobian to foot 2 heel
	jacobian_2(model, get_q⁺(x), body = :calf_2, mode = :ee)

jacobian to foot 1 toe
	jacobian_3(model, get_q⁺(x), body = :foot_1, mode = :ee)

jacobian to foot 2 toe
	jacobian_3(model, get_q⁺(x), body = :foot_2, mode = :ee)

gravity compensating torques
	initial_torque(model, q, h)[model.idx_u] # returns torques for system
	initial_torque(model, q, h) # returns torques for system and contact forces

	-NOTE: this method is iterative and is not differentiable (yet)
		-TODO: implicit function theorem to compute derivatives

contact models
	- maximum-dissipation principle:
		contact_constraints(model, T)
	- no slip:
		contact_no_slip_constraints(T)

loop constraint
	- constrain variables at two indices to be the same:
		loop_constraints(model, x_idx, t1_idx, t2_idx)
"""

function get_q⁺(x)
	view(x, model.nq .+ (1:model.nq))
end

function ellipse_traj(x_start, x_goal, z, T)
	dist = x_goal - x_start
	a = 0.5 * dist
	b = z

	z̄ = 0.0

	x = circular_projection_range(x_start, stop = x_goal, length = T)
	# x = range(x_start, stop = x_goal, length = T)

	z = sqrt.(max.(0.0, (b^2) * (1.0 .- ((x .- (x_start + a)).^2.0) / (a^2.0))))

	return x, z
end

function circular_projection_range(start; stop=1.0, length=10)
	dist = stop - start
	θr = range(π, stop=0, length=length)
	r = start .+ dist * ((1 .+ cos.(θr))./2)
	return r
end

# Horizon
T = 2
h = 0.01

# Permutation matrix
perm = @SMatrix [1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0;
                 0.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0;
				 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0;
				 0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0;
				 0.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0;
				 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0;
				 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0;
				 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0;
				 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0]

# Configurations
# 1: x pos
# 2: z pos
# 3: torso angle (rel. to downward vertical)
# 4: thigh 1 angle (rel. to downward vertical)
# 5: calf 1 (rel. to downward vertical)
# 6: thigh 2 (rel. to downward vertical)
# 7: calf 2 (rel. to downward vertical)
# 8: foot 1 (rel. to downward vertical)
# 9: foot 2 (rel. to downward vertical)
function initial_configuration(model, θ_torso, θ_thigh_1, θ_leg_1, θ_thigh_2)
    q1 = zeros(model.nq)
    q1[3] = θ_torso
    q1[4] = θ_thigh_1
    q1[5] = θ_leg_1
    z1 = model.l_thigh1 * cos(q1[4]) + model.l_calf1 * cos(q1[5])
    q1[6] = θ_thigh_2
    q1[7] = 1.0 * acos((z1 - model.l_thigh2 * cos(q1[6])) / model.l_calf2)
	q1[2] = z1
	q1[8] = pi / 2.0
	q1[9] = pi / 2.0

    return q1
end

# q1 = initial_configuration(model, -π / 50.0, -π / 6.5, π / 10.0 , -π / 10.0)
# q1 = initial_configuration(model, -0.00*π / 50.0, -π / 5.6, π / 25.0 , -π / 7.0)
q1 = initial_configuration(model, -0.00*π / 50.0, -π / 6.0, π / 12.0 , -π / 6.0)
pf1 = kinematics_3(model, q1, body = :foot_1, mode = :com)
pf2 = kinematics_3(model, q1, body = :foot_2, mode = :com)

ph1 = kinematics_3(model, q1, body = :foot_1, mode = :heel)
ph2 = kinematics_3(model, q1, body = :foot_2, mode = :heel)[1]
ph2 = kinematics_2(model, q1, body = :calf_2, mode = :ee)[1]
q1[1]

strd = 2 * (pf2 - pf1)[1]

qT = Array(perm) * copy(q1)
qT[1] += 0.5 * strd
q_ref = linear_interpolation(q1, qT, T+1)
x_ref = configuration_to_state(q_ref)
visualize!(vis, model, q_ref, Δt = h)

# Control

# ul <= u <= uu
_uu = Inf * ones(model.m)
_uu[model.idx_u] .= Inf
# _uu[end] = 2.0 * h
_ul = zeros(model.m)
_ul[model.idx_u] .= -Inf
# _ul[end] = 0.25 * h
ul, uu = control_bounds(model, T, _ul, _uu)

qL = [-Inf; -Inf; q1[3] - pi / 750.0; q1[4:end] .- pi / 6.0; -Inf; -Inf; q1[3] - pi / 50.0; q1[4:end] .- pi / 6.0]
qU = [Inf; q1[2] + 0.001; q1[3] + pi / 750.0; q1[4:end] .+ pi / 6.0; Inf; Inf; q1[3:end] .+ pi / 6.0]
qL[8] = q1[8] - pi / 10.0
qL[9] = q1[9] - pi / 10.0
qL[9 + 8] = q1[8] - pi / 10.0
qL[9 + 9] = q1[9] - pi / 10.0

qU[8] = q1[8] + pi / 10.0
qU[9] = q1[9] + pi / 10.0
qU[9 + 8] = q1[8] + pi / 10.0
qU[9 + 9] = q1[9] + pi / 10.0

xl, xu = state_bounds(model, T,
    qL, qU,
    x1 = [q1; q1],
	xT = [q1; q1])

# Objective
include_objective(["velocity", "nonlinear_stage", "control_velocity"])

x0 = configuration_to_state(q_ref)

# penalty on slack variable
obj_penalty = PenaltyObjective(1.0e5, model.m)

obj = MultiObjective([obj_penalty])

# Constraints
include_constraints(["contact"])
con_contact = contact_constraints(model, T)
con = multiple_constraints([con_contact])

# Problem
prob = trajectory_optimization_problem(model,
               obj,
               T,
			   h = h,
               xl = xl,
               xu = xu,
               ul = ul,
               uu = uu,
               con = con)

# trajectory initialization
u0 = [[1.0e-2 * randn(model.nu); 0.01 * randn(model.m - model.nu)] for t = 1:T-1] # random controls

# Pack trajectories into vector
z0 = pack(x0, u0, prob)

# Solve
# NOTE: run multiple times to get good trajectory
@time z̄, info = solve(prob, copy(z0),
    nlp = :ipopt,
	max_iter = 2000,
    tol = 1.0e-5, c_tol = 1.0e-5, mapl = 5,
    time_limit = 60 * 3)

@show check_slack(z̄, prob)
x̄, ū = unpack(z̄, prob)
# _tf, _t, h̄ = get_time(ū)

q = state_to_configuration(x̄)
u = [u[model.idx_u] for u in ū]
γ = [u[model.idx_λ] for u in ū]
b = [u[model.idx_b] for u in ū]
ψ = [u[model.idx_ψ] for u in ū]
η = [u[model.idx_η] for u in ū]

plot(u, linetype = :steppost)

μ = model.μ
@save joinpath(@__DIR__, "flamingo_stand_100hz.jld2") q u γ b ψ η μ h

vis = Visualizer()
render(vis)
# open(vis)
visualize!(vis, model, q, Δt = h)

max_issues_repo_name: thowell/DirectMotionPlanning


text:
function cma_es_animated(f,x0;tol=1e-5,maxiter=1e3, animate=true,
       xrange=[-2., 3.], yrange=[-2.,6.])
  if animate
     # make a contour plot of the function we're minimizing. This is for
     # illustrating; you wouldn't have this normally
     x = range(xrange[1],xrange[2], length=100)
     y = range(yrange[1],yrange[2], length=100)
     c = contour(x,y,(x,y) -> log(f([x,y])) )
     anim = Animation()
   else
     anim = nothing
   end
    N=length(x0)
    λ=4+floor(3*log(N))
    m=x0
      σ=0.3
      μ=floor(λ/2)
      w=log(μ+0.5).-log.(1:μ)
      w=w/sum(w)
      μ_w=(sum(w)^2)/sum(w.^2)
      c_c=(4+μ_w/N)/(N+4+2*μ_w/N)
      c_σ=(μ_w+2)/(N+μ_w+5)
      c_1=2/((N+1.3)^2+μ_w)
      c_μ=min(1-c_1,
              2*(μ_w-2+1/μ_w)/((N+2)^2+μ_w) )
      damp_σ=1+c_σ+ 2*max(0.0,
                  sqrt( (μ_w-1)/(N+1) )-1    )
      p_c=zeros(N,1)
      p_σ=p_c
      B=I(N)
      D=ones(N,1)
      C=B*Diagonal(D.^2)*B'
      invsqrtC=B*Diagonal( D.^-1 )*B'
      eigeneval=0

      χ_N=sqrt(N)*(1-1/(4*N) +1/(21*N^2) )

      count=0
      normdiff=Inf
      N=length(x0)

      fx=zeros(Integer(λ))

      C=Diagonal(ones(N))
      equis=zeros(Integer(λ),N)

      function loop(equis,p_c,p_σ,c_1,c_μ,c_c,σ,C,m,
        B,D,invsqrtC,eigeneval)

          for count in 1:maxiter

          oldxbest=equis[1,:]
          for i in 1:Integer(λ)
            equis[i,:]=(σ.*B*(D.*randn(N,1))+m)'
            fx[i]=f(equis[i,:])
            #count=count+1
          end
          if animate
          # plot the search points
          p = deepcopy(c)
          scatter!(p, equis[:,1], equis[:,2], markercolor=:black,
          markeralpha=0.5, legend=false, xlims=xrange, ylims=yrange)
          end
          equis=equis[sortperm(fx),:]
          if animate
          # plot the best point so far
          scatter!(p, equis[1:Integer(μ),1], equis[1:Integer(μ),2],
           markercolor=:red, legend=false)
          end
          frame(anim)
          dm=(equis[1:Integer(μ),:]'-repeat(m,1,Integer(μ)))*w
          m=m+dm
          #dm=-dm
          p_σ=(1-c_σ)*p_σ+
                  ((sqrt(c_σ*(2-c_σ)*μ_w)*invsqrtC)*dm)/σ

          hsig=norm(p_σ)/sqrt(1-(1-c_σ)^(2*count/λ))χ_N<1.4+2/(N+1)

          p_c=(1-c_c)*p_c+(hsig*sqrt(c_c*(2-c_c)*μ)/σ) *dm

          aux=(1/σ).*(equis[1:Integer(μ),:]'-repeat(m,1,Integer(μ)) )
          C=(1-c_1-c_μ)*C+
                  c_1*(p_c*p_c'+(1-hsig)*c_c*(2-c_c)*C)+
                  c_μ*aux*Diagonal(w)*aux'
          σ=σ*exp( (c_σ/damp_σ)*(norm(p_σ)/χ_N-1) )

          if count-eigeneval>λ/(c_1+c_μ)/N/10
            eigeneval=count
            C=UpperTriangular(C)+UpperTriangular(C)'-Diagonal(C)
            B=eigvecs(C)
            D=eigvals(C).^(0.5)
          end
          invsqrtC=B'*Diagonal( D.^-1 )*B
          count=count+1
          normdiff=norm(equis[1,:]-oldxbest)
          #normdiff=norm(dm)
          if normdiff< tol || maximum(D)>1e7*minimum(D)
                  break
          end
          end
          return(f(equis[1,:]),equis[1,:],anim)
      end
      loop(equis,p_c,p_σ,c_1,c_μ,c_c,σ,C,m,B,D,invsqrtC,eigeneval)
end


text:
#=
This file is part of the replication code for: Hasenzagl, T., Pellegrino, F., Reichlin, L., & Ricco, G. (2020). A Model of the Fed's View on Inflation.
Please cite the paper if you are using any part of the code for academic work (including, but not limited to, conference and peer-reviewed papers).
=#

function get_random_disturbance(x, Σ)

     # Init random perturbation
     rnd_per     = zeros(x);
     rnd_per_ind = sum(Σ, dims=2)[:] .!= 0;

     if sum(rnd_per_ind) > 0

          # Get random perturbation
          rnd_per_cov                   = Σ[findall(rnd_per_ind), findall(rnd_per_ind)];
          rnd_per[findall(rnd_per_ind)] = rand(MvNormal(zeros(sum(rnd_per_ind)), rnd_per_cov));
     end

     return rnd_per;
end


text:
using Luxor, Test, Colors

function imagematrix()
    Drawing(2, 2, :image)

    # get image as matrix
    mat = image_as_matrix()
    @test mat[1] == Colors.ARGB32(0.0,0.0,0.0,0.0)
    @test mat[1] == mat[2] == mat[3] == mat[4]

    # make it red, then get image as matrix again
    setcolor("red")
    paint()
    mat = image_as_matrix()

    # make it green, then get image as matrix again
    setcolor(0, 1, 0)
    paint()
    mat = image_as_matrix()
    @test mat[1] == Colors.ARGB32(0.0, 1.0, 0.0, 1.0)
    @test mat[1] == mat[2] == mat[3] == mat[4]

    @test finish() == true

    # make 20 by 20 (400) grid
    Drawing(20, 20, :image)
    Luxor.origin()
    tiles = Tiler(20, 20, 20, 20, margin=0)
    background("white")
    sethue("blue")
    # fill odd ones with blue
    for (pos, n) in tiles
        isodd(n) && box(pos, 1, 1, :fill)
    end
    mat = image_as_matrix()

    # half the 400 are white
    @test count(i -> i == ARGB32(1, 1, 1, 1), mat) == 200

    # half are blue - 0r, 0g, 1b, 1a
    @test count(i -> i == ARGB32(0, 0, 1, 1), mat) == 200
    @test finish() == true
end

imagematrix()

println("...finished test")


text:
using moeum
using Base.Test

#default test
mo_a = moeum.MOEUM(dat_sori = ["dat1", "dat2"], hol_sori = [[1, 2], [3, 4]], name = "mo_a")
dat1 = mo_a.select("dat1")
dat1.describe()

mo_a = moeum.from_csv("mo_a.csv")
dat1 = mo_a.select("dat1")
dat1.describe()

mo_a.where("dat1 >= 2")


text:
using Pkg
Pkg.add("Clang")
Pkg.activate(joinpath(@__DIR__, ".."))

using Clang
using AMCL_jll

# LIBCLANG_HEADERS are those headers to be wrapped.
const LIBAMCL_INCLUDE = joinpath(dirname(AMCL_jll.libamcl_core_path), "..", "include", "amcl") |> normpath
const LIBAMCL_HEADERS = [joinpath(LIBAMCL_INCLUDE, header) for header in readdir(LIBAMCL_INCLUDE) if endswith(header, ".h")]
#const LIBAMCL_HEADERS = [joinpath(LIBAMCL_INCLUDE, "amcl.h")]

clang_args = ["-I$LIBAMCL_INCLUDE"]
for header in find_std_headers()
    push!(clang_args, "-I$header")
end

function libmap(include_path)
    lib = strip(splitext(basename(include_path))[1])

    # Core APIs
    if lib in ("amcl", "utils", "version", "randapi", "ecdh_support", "pbc_support", "rsa_support")
        return "libamcl_core"
    end

    # BLS24 and BLS48 specific
    lib != "pair192_BLS24" || return "libamcl_pairing_BLS24"
    lib != "bls192_BLS24" || return "libamcl_bls_BLS24"
    lib != "wcc192_BLS24" || return "libamcl_wcc_BLS24"
    lib != "mpin192_BLS24" || return "libamcl_mpin_BLS24"
    lib != "pair256_BLS48" || return "libamcl_pairing_BLS48"
    lib != "bls256_BLS48" || return "libamcl_bls_BLS48"
    lib != "wcc256_BLS48" || return "libamcl_wcc_BLS48"
    lib != "mpin192_BLS24" || return "libamcl_mpin_BLS48"

    # Curve APIs
    m = match(r"(?:ecdh_|ecp_|fp_)(.*)", lib)
    isnothing(m) || return "libamcl_curve_$(m.captures[1])"

    # Pairing APIs
    m = match(r"(?:ecp\d+_|fp\d+_|pair_)(.*)", lib)
    isnothing(m) || return "libamcl_pairing_$(m.captures[1])"

    # RSA APIs
    m = match(r"ff_(.*)", lib)
    isnothing(m) || return "libamcl_rsa_$(m.captures[1])"

    # big number APIs
    lib != "big_256_56" || return "libamcl_curve_BN254"
    lib != "big_336_60" || return "libamcl_curve_HIFIVE"
    lib != "big_384_56" || return "libamcl_curve_NIST384"
    lib != "big_384_58" || return "libamcl_curve_BLS381"
    lib != "big_416_60" || return "libamcl_curve_C41417"
    lib != "big_448_58" || return "libamcl_curve_GOLDILOCKS"
    lib != "big_464_60" || return "libamcl_curve_BLS461"
    lib != "big_480_56" || return "libamcl_curve_BLS24"
    lib != "big_512_56" || return "libamcl_curve_NUMS512E"
    lib != "big_512_60" || return "libamcl_rsa_4096"
    lib != "big_528_60" || return "libamcl_curve_NIST521"
    lib != "big_560_58" || return "libamcl_curve_BLS48"
    lib != "big_1024_58" || return "libamcl_rsa_2048"

    return "libamcl_$lib"
end

wc = init(
    headers = LIBAMCL_HEADERS,
    output_file = joinpath(@__DIR__, "libamcl_api.jl"),
    common_file = joinpath(@__DIR__, "libamcl_common.jl"),
    clang_includes = vcat(LIBAMCL_INCLUDE, CLANG_INCLUDE),
    clang_args = clang_args,
    header_wrapped = (root, current)->root == current,
    header_library = libmap,
)

#run(wc, false)

module AMCL
    const int8_t = Int8
    const int16_t = Int16
    const int32_t = Int32
    const int64_t = Int64
    const uint8_t = UInt8
    const uint16_t = UInt16
    const uint32_t = UInt32
    const uint64_t = UInt64
    const NLEN_1024_58 = 1
    const NLEN_256_56 = 1
    const NLEN_336_60 = 1
    const NLEN_384_56 = 1
    const NLEN_384_58 = 1
    const NLEN_416_60 = 1
    const NLEN_448_58 = 1
    const NLEN_464_60 = 1
    const NLEN_480_56 = 1
    const NLEN_512_56 = 1
    const NLEN_512_60 = 1
    const NLEN_528_60 = 1
    const NLEN_560_58 = 1
    include(joinpath(@__DIR__, "libamcl_common.jl"))
end

const BIG_TYPES = r"^BIG_\d+_\d+$"
const DBIG_TYPES = r"^DBIG_\d+_\d+$"
const FP_TYPES = r"^FP\d*_[A-Z0-9]+$"
const ECP_TYPES = r"^ECP\d*_[A-Z0-9]+$"
const RSA_TYPES = r"^rsa_(?:private|public)_key_\d+$"
const OTHER_TYPES = r"^(?:hash256|hash512|sha3|csprng|gcm|amcl_aes|pktype|octet)$"
const ALL_TYPES = (BIG_TYPES, DBIG_TYPES, FP_TYPES, ECP_TYPES, RSA_TYPES, OTHER_TYPES)

mathches_any(s::String, patterns::Regex...) = mathches_any(s, patterns)
mathches_any(s::String, patterns::NTuple{N,Regex}) where {N} = any(p->!isnothing(match(p, s)), patterns)

open(joinpath(@__DIR__, "libamcl_base_methods.jl"), "w") do f
    println(f, "# Automatically generated with generator.jl")

    for sym in names(AMCL, all=true)
        T = getproperty(AMCL, sym)
        T isa Type || continue
        n = String(sym)

        if mathches_any(n, ALL_TYPES)
            # pointer(::T)
            println(f, "Base.pointer(a::$n) = reinterpret(Ptr{$n}, pointer_from_objref(a))")
            # unsafe_convert(::Type{Ptr{T}}, ::T)
            println(f, "Base.unsafe_convert(::Type{Ptr{$n}}, a::$n) = pointer(a)")
        else
            @info "skipped: $n"
            continue
        end
        if mathches_any(n, BIG_TYPES)
            println(f, "Base.:(==)(a::$n, b::$n) = iszero($(n)_comp(a, b))")
            println(f, "Base.isless(a::$n, b::$n) = isone(-$(n)_comp(a, b))")
            println(f, "$n() = (a = $n(undef); $(n)_zero(a); a)")
        end
        if mathches_any(n, DBIG_TYPES)
            println(f, "Base.:(==)(a::$n, b::$n) = iszero($(n[2:end])_dcomp(a, b))")
            println(f, "Base.isless(a::$n, b::$n) = isone(-$(n[2:end])_dcomp(a, b))")
            println(f, "$n() = (a = $n(undef); $(n[2:end])_dzero(a); a)")
        end
        if mathches_any(n, FP_TYPES)
            println(f, "Base.:(==)(a::$n, b::$n) = isone($(n)_equals(a, b))")
            println(f, "$n() = (a = $n(undef); $(n)_zero(a); a)")
        end
        if mathches_any(n, ECP_TYPES)
            println(f, "Base.:(==)(a::$n, b::$n) = isone($(n)_equals(a, b))")
            println(f, "$n() = (a = $n(undef); $(n)_inf(a); a)")
        end
        if mathches_any(n, r"octet")
            println(f, "Base.:(==)(a::octet, b::octet) = isone(AMCL.OCT_comp(a, b))")
        end
    end
end


text:
# module Hangul

# 자음 모음
const Consonant = :Consonant
const Vowel     = :Vowel

# Hangul Compatibility Letter
primitive type Letter{T} <: AbstractChar 32 end

# Hangul Syllables
struct Syllable
    initial::Union{Nothing,Letter{Consonant}}
    medial::Union{Nothing,Letter{Vowel}}
    final::Union{Nothing,Letter{Consonant}}
end

struct ComposeError <: Exception
    message::String
end

struct DecomposeError <: Exception
    message::String
end

struct LetterError <: Exception
    message::String
end

Syllable(i::Union{Nothing,Letter{Consonant}}, m::Union{Nothing,Letter{Vowel}}) = Syllable(i, m, nothing)

Letter{T}(c::UInt32) where T = reinterpret(Letter{T}, c)
Base.codepoint(l::Letter{T}) where T = reinterpret(UInt32, l)

"""
    Jamo
"""
function Jamo(c::Char)::Union{Letter,YetJamo} # LetterError
    Jamo(codepoint(c))
end

function Jamo(c::Union{UInt16,UInt32})::Union{Letter,YetJamo} # LetterError
    # 0x3130 : undef
    # 0x318f : undef
    if 0x3131 <= c <= 0x314e || # ㄱ   ... ㅎ  # Consonant letters
       0x3165 <= c <= 0x3186    # ㄴㄴ ... ㆆ  # Old consonant letters
        Letter{Consonant}(c)
    elseif 0x314f <= c <= 0x3163 || # ㅏ ... ㅣ   # Vowel letters
           0x3187 <= c <= 0x318e    # ㆇ ... ㆎ   # Old vowel letters
        Letter{Vowel}(c)
    elseif is_yetjamo(c)
        YetJamo(c)
    else
        throw(LetterError(string(c)))
    end
end

macro _w_str(s::String) # Vector{Char}
    Expr(:call, ()->first.(split(s, ' ')))
end

const 모음채움 = Letter{Vowel}(0x3164)     # HANGUL FILLER
const 자음채움 = Letter{Consonant}(0x318f) # ?

const 초성표 = Letter{Consonant}.(_w"ㄱ ㄲ ㄴ ㄷ ㄸ ㄹ ㅁ ㅂ ㅃ ㅅ ㅆ ㅇ ㅈ ㅉ ㅊ ㅋ ㅌ ㅍ ㅎ")
const 중성표 = Letter{Vowel}.(_w"ㅏ ㅐ ㅑ ㅒ ㅓ ㅔ ㅕ ㅖ ㅗ ㅘ ㅙ ㅚ ㅛ ㅜ ㅝ ㅞ ㅟ ㅠ ㅡ ㅢ ㅣ")
const 종성표 = Letter{Consonant}.(vcat([자음채움], _w"ㄱ ㄲ ㄳ ㄴ ㄵ ㄶ ㄷ ㄹ ㄺ ㄻ ㄼ ㄽ ㄾ ㄿ ㅀ ㅁ ㅂ ㅄ ㅅ ㅆ ㅇ ㅈ ㅊ ㅋ ㅌ ㅍ ㅎ"))

const 중성오프셋 = 28  # length(종성표)
const 초성오프셋 = 588 # length(중성표) * 중성오프셋  # 21 * 28
const 가 = 0xac00 # codepoint('가')
const 힣 = 0xd7a3 # codepoint('힣')

indexof(순서::Vector{Letter{T}}, 소리::Letter{T}) where T = findall(in(소리), 순서)[1]

function throw_compose_error(syl::Syllable) # ComposeError
    throw(ComposeError(string(syl)))
end

function throw_decompose_error(c::Char) # DecomposeError
    throw(DecomposeError(string(c)))
end

function to_char(syl::Syllable)::Char # ComposeError
    초성, 중성, 종성 = syl.initial, syl.medial, syl.final
    if 초성 isa Nothing
        if 중성 isa Nothing
            if 종성 isa Nothing
                throw_compose_error(syl)
            else
                Char(종성)
            end
        else # 초성x, 중성o
            if 종성 isa Nothing
                Char(중성)
            else
                throw_compose_error(syl)
            end
        end
    else # 초성o
        if 중성 isa Nothing
            if 종성 isa Nothing
                Char(초성)
            else
                throw_compose_error(syl)
            end
        else # 초성o, 중성o
            초성인덱스 = (indexof(초성표, 초성) - 1) * 초성오프셋
            중성인덱스 = (indexof(중성표, 중성) - 1) * 중성오프셋
            value = 가 + 초성인덱스 + 중성인덱스
            if 종성 isa Nothing
            else
                종성인덱스 = indexof(종성표, 종성) - 1
                value += 종성인덱스
            end
            Char(value)
        end
    end
end

function to_syllable(c::Char)::Syllable # DecomposeError
    code = codepoint(c)
    if 가<= code <= 힣
        value::UInt32 = code - 가
        초성인덱스 = Int(round(value / 초성오프셋, RoundToZero)) + 1
        n = value % 초성오프셋
        중성인덱스 = Int(round(n / 중성오프셋, RoundToZero)) + 1
        종성인덱스 = Int(n % 중성오프셋) + 1
        Syllable(초성표[초성인덱스], 중성표[중성인덱스], 1 == 종성인덱스 ? nothing : 종성표[종성인덱스])
    else
        if c in 초성표
            Syllable(Letter{Consonant}(c), nothing, nothing)
        elseif c in 중성표
            Syllable(nothing, Letter{Vowel}(c), nothing)
        else
            throw_decompose_error(c)
        end
    end
end

# module Hangul


text:
# This file contains functions for the DEA optimizers
"""
    AbstractDEAOptimizer
An abstract type representing a DEA optimizer.
"""
abstract type AbstractDEAOptimizerModel end

"""
    DEAOptimizer(optimizer; time_limit, silent)
An data structure storing the configuration of a DEA optimizer.

# Optimizer specification:
- `LP`: linear programming default optimizer, GLPK.
- `NLP`: nonlinear programmin default optimizer, Ipopt.
- Any JuMP supported solver.

# Optional Arguments
- `time_limit=:60`: time limit in seconds.
- `silent=:true`: suppress optimizer output.

# Examples
```jldoctest
julia> myoptimizer = DEAOptimizer(GLPK.Optimizer, time_limit = 10, silent = true);
```
"""
struct DEAOptimizer <: AbstractDEAOptimizerModel
    optimizer
    time_limit::Float64
    silent::Bool

    function DEAOptimizer(optimizer; time_limit = 60, silent = true)

        if optimizer == :LP
            optimizer =  GLPK.Optimizer
        elseif optimizer == :NLP
            optimizer = Ipopt.Optimizer
        end

        new(optimizer, time_limit, silent)
    end

end

"""
    newdeamodel(DEAOptimizer)
Generate a new JuMP model for DEA with the specified optimizer.

This function is used internally and for packages that want to extend the functionality of this package.

# Examples
```jldoctest
julia> deamodel = newdeamodel(DEAOptimizer(GLPK.Optimizer));
```
"""
function newdeamodel(deaoptimizer::DEAOptimizer)::Model
    deamodel = Model(deaoptimizer.optimizer)

    set_time_limit_sec(deamodel, deaoptimizer.time_limit)

    if deaoptimizer.silent
        set_silent(deamodel)
    end

    return deamodel
end


text:
import Linda

import GLPK
import MathOptInterface
const MOI = MathOptInterface

import Tulip

# GLPK does not support range constraints, so we bridge
MOI.get(::GLPK.Optimizer, ::MOI.BarrierIterations) = 0
MOI.get(::GLPK.Optimizer, ::MOI.SimplexIterations) = 0

function test_colgen1()
    # Here we optimize a problem with one load in one house
    # Aggregator must enforce more restrictive total load bounds,
    # and coordinate for price.

    # Problem data
    T = 2
    ptou = [1.0, 1.0]
    pmkt = [1.0, 5.0]
    ptot_min = [0.0, 0.0]
    ptot_max = [1.0, 1.0]

    R = 1  # number of resources

    # Instantiate master problem
    rmp = MOI.Bridges.full_bridge_optimizer(GLPK.Optimizer(), Float64)
    MOI.set(rmp, MOI.Silent(), true)

    fobj = MOI.ScalarAffineFunction(MOI.ScalarAffineTerm{Float64}[], 0.0)
    MOI.set(rmp, MOI.ObjectiveFunction{typeof(fobj)}(), fobj)
    MOI.set(rmp, MOI.ObjectiveSense(), MOI.MIN_SENSE)

    pagg = MOI.add_variables(rmp, T)  # Aggregator net load

    # Bounds
    for t in 1:T
        MOI.add_constraint(rmp,
            MOI.SingleVariable(pagg[t]),
            MOI.Interval(ptot_min[t], ptot_max[t])
        )
    end

    # Set objective coefficient
    for t in 1:T
        MOI.modify(rmp,
            MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(),
            MOI.ScalarCoefficientChange(pagg[t], pmkt[t] - ptou[t])
        )
    end

    # Convexity constraints
    con_cvx = MOI.ConstraintIndex[]
    for r in 1:R
        cidx = MOI.add_constraint(rmp,
            MOI.ScalarAffineFunction(MOI.ScalarAffineTerm{Float64}[], 0.0),
            MOI.EqualTo(1.0)
        )
        push!(con_cvx, cidx)
    end

    # Linking constraints
    # This is simply -pagg[t] + p1[t] + ... + pR[t] + yp[t] - ym[t] == 0
    con_link = MOI.ConstraintIndex[]
    for t in 1:T
        cidx = MOI.add_constraint(rmp,
            MOI.ScalarAffineFunction([
                MOI.ScalarAffineTerm(-1.0, pagg[t])
            ], 0.0),
            MOI.EqualTo(0.0)
        )
        push!(con_link, cidx)
    end

    mp = Linda.Master(rmp, con_cvx, con_link, zeros(T), pagg, MOI.VariableIndex[])

    # Create the resource
    l1 = DR.DER.DeferrableLoad(
        index=0, T=T, ncycles=1,
        cycle_begin=[1], cycle_end=[T],
        cycle_energy_min=[1.5], cycle_energy_max=[2.0],
        pwr_min=zeros(T), pwr_max=2.0 * ones(T),
        binflag=false
    )

    l2 = DR.DER.ShiftableLoad(
        T=T,
        cycle_begin_min=1,
        cycle_begin_max=2,
        cycle_length=1,
        cycle_load=[1.0]
    )

    # Create sub-problem
    h1 = DR.DER.House(
        1, T,
        [0.0, 0.0], [2.0, 2.0], ptou,
        [l1, l2]
    )
    sp1 = MOI.Bridges.full_bridge_optimizer(GLPK.Optimizer(), Float64)
    o1 = DR.build_oracle(sp1, h1)

    # Add one artifical column
    # This ensures that the master is feasible
    col0 = Linda.Column(
        100.0,
        h1.index,
        [1, 2], [0.0, 0.0], true
    )
    Linda.add_column!(mp, col0)

    # Solve with column generation
    env = Linda.LindaEnv()
    env.verbose.val = 1
    env.num_cgiter_max.val = 10
    Linda.solve_colgen!(env, mp, [o1])

    # Check solution status
    pagg_ = [
        MOI.get(rmp, MOI.VariablePrimal(), pagg[t])
        for t in 1:T
    ]

    # Solution should be optimal
    @test mp.mp_status == MOI.OPTIMAL

    @test isapprox(pagg_[1], 1.0, atol = 1e-5)
    @test isapprox(pagg_[2], 0.5, atol = 1e-5)

    return nothing
end

function test_colgen2()
    # Here we optimize a model with one load in a first house, and a battery
    # in a second house.
    # The aggregator must coordinate so that the battery is charged at the first
    # period and discharged at the second.

    # Problem data
    T = 2
    # Prices
    # Houses do not see that second period is in fact more expensive
    ptou = [2.0, 1.0]
    pmkt = [1.0, 5.0]
    ptot_min = [0.0, 0.0]
    ptot_max = [10.0, 10.0]

    R = 2  # number of resources

    # Instantiate master problem
    rmp = MOI.Bridges.full_bridge_optimizer(GLPK.Optimizer(), Float64)
    MOI.set(rmp, MOI.Silent(), true)

    fobj = MOI.ScalarAffineFunction(MOI.ScalarAffineTerm{Float64}[], 0.0)
    MOI.set(rmp, MOI.ObjectiveFunction{typeof(fobj)}(), fobj)
    MOI.set(rmp, MOI.ObjectiveSense(), MOI.MIN_SENSE)

    pagg = MOI.add_variables(rmp, T)  # Aggregator net load
    # We add artificial variables to ensure RMP is feasible
    yp = MOI.add_variables(rmp, T)
    ym = MOI.add_variables(rmp, T)

    # Bounds
    for t in 1:T
        MOI.add_constraint(rmp,
            MOI.SingleVariable(pagg[t]),
            MOI.Interval(ptot_min[t], ptot_max[t])
        )

        MOI.add_constraint(rmp, MOI.SingleVariable(yp[t]), MOI.GreaterThan(0.0))
        MOI.add_constraint(rmp, MOI.SingleVariable(ym[t]), MOI.GreaterThan(0.0))
    end

    # Set objective coefficient
    for t in 1:T
        MOI.modify(rmp,
            MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(),
            MOI.ScalarCoefficientChange(pagg[t], pmkt[t] - ptou[t])
        )

        MOI.modify(rmp,
            MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(),
            MOI.ScalarCoefficientChange(yp[t], 1000.0)
        )
        MOI.modify(rmp,
            MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(),
            MOI.ScalarCoefficientChange(ym[t], 1000.0)
        )
    end

    # Convexity constraints
    con_cvx = MOI.ConstraintIndex[]
    for r in 1:R
        cidx = MOI.add_constraint(rmp,
            MOI.ScalarAffineFunction(MOI.ScalarAffineTerm{Float64}[], 0.0),
            MOI.EqualTo(1.0)
        )
        push!(con_cvx, cidx)
    end

    # Linking constraints
    # This is simply -pagg[t] + p1[t] + ... + pR[t] + yp[t] - ym[t] == 0
    con_link = MOI.ConstraintIndex[]
    for t in 1:T
        cidx = MOI.add_constraint(rmp,
            MOI.ScalarAffineFunction([
                MOI.ScalarAffineTerm(-1.0, pagg[t]),
                MOI.ScalarAffineTerm( 1.0, yp[t]),
                MOI.ScalarAffineTerm(-1.0, ym[t])
            ], 0.0),
            MOI.EqualTo(0.0)
        )
        push!(con_link, cidx)
    end

    mp = Linda.Master(rmp, con_cvx, con_link, zeros(T), pagg, MOI.VariableIndex[])

    # Create the resources
    # House 1: one deferrable load
    l1 = DR.DER.DeferrableLoad(
        index=0, T=T, ncycles=1,
        cycle_begin=[1], cycle_end=[T],
        cycle_energy_min=[3.0], cycle_energy_max=[4.0],
        pwr_min=[0.0, 0.0], pwr_max=[2.0, 2.0],
        binflag=false
    )
    h1 = DR.DER.House(
        1, T,
        [0.0, 0.0], [2.0, 2.0], ptou,
        [l1]
    )
    sp1 = MOI.Bridges.full_bridge_optimizer(GLPK.Optimizer(), Float64)
    o1 = DR.build_oracle(sp1, h1)

    # House 2: one fixed load and one battery
    # House 2 should discharge its battery at T=2
    l2 = DR.DER.FixedLoad(index=0, T=T, load=[1.0, 1.0])
    bat2 = DR.DER.Battery(
        index=0, T=T,
        soc_min=[0.0, 0.0], soc_max=[10.0, 10.0],
        charge_pwr_max=6.0, discharge_pwr_max=6.0,
        charge_eff=1.0,
        discharge_eff=1.0
    )

    # Create sub-problem
    # Second house can inject power into the grid
    h2 = DR.DER.House(
        2, T,
        [-10.0, -10.0], [10.0, 10.0], ptou,
        [l2, bat2]
    )
    sp2 = MOI.Bridges.full_bridge_optimizer(GLPK.Optimizer(), Float64)
    o2 = DR.build_oracle(sp2, h2)

    # Add one initial column per resource
    # This ensures that the master is feasible
    for o in [o1, o2]
        Linda.optimize!(o)
        col, rc = Linda.get_columns(o)[1]
        Linda.add_column!(mp, col)
    end

    # Solve with column generation
    env = Linda.LindaEnv()
    env.verbose.val = 0
    env.num_cgiter_max.val = 10
    Linda.solve_colgen!(env, mp, [o1, o2])

    # Check solution status
    pagg_ = [
        MOI.get(rmp, MOI.VariablePrimal(), pagg[t])
        for t in 1:T
    ]

    # Solution should be optimal
    @test mp.mp_status == MOI.OPTIMAL

    # Check aggregated load
    @test isapprox(pagg_[1], 5.0, atol = 1e-5)
    @test isapprox(pagg_[2], 0.0, atol = 1e-5)

    # TODO: check consumption of each house and battery levels

    return nothing
end

@testset "ColGen" begin
    test_colgen1()
    test_colgen2()
end


text:
"""
# Examples
```julia
f = x -> x == :a
@test ChainMap.negate(f)(:b)
```
"""
negate(f) = (args...; kwargs...) -> !(f(args...; kwargs...))

function dots_to_back(o::DataStructures.OrderedDict)
    is_dots = (k, v) -> MacroTools.isexpr(k, :...)
    to_back = filter(is_dots, o)
    if length(to_back) > 1
        error("Can splat no more than one positional argument")
    end
    @chain begin
        o
        filter(negate(is_dots), _)
        merge(_, to_back)
    end
end

function parameters_to_front(o::DataStructures.OrderedDict)
    is_parameters = (k, v) -> double_match(k, :parameters, :...)
    to_front = filter(is_parameters, o)
    if length(to_front) > 1
        error("Can splat no more than one keyword argument")
    end
    @chain begin
        o
        filter(negate(is_parameters), _)
        merge(to_front, _)
    end
end

"""
    split_anonymous(e::Expr)

Return the component anonymous function and arguments for [`unweave`](@ref).

```julia
e = :(~_ + 1)
ChainMap.split_anonymous(e)
```
"""
function split_anonymous(e::Expr)
    d = Dict()
    e_replace = replace_record!(e, d)

    d_reorder = @chain begin
        d
        DataStructures.OrderedDict(_)
        parameters_to_front
        dots_to_back
    end

    anonymous_function = @chain begin
        d_reorder
        values
        Expr(:tuple, _...)
        Expr(:->, _, e_replace)
    end

    (anonymous_function, keys(d_reorder))
end

export unweave

unweave(e) = e
unweave(f, e) = e

"""
    @unweave e

Interprets `e` as a function with its positional arguments wrapped in tildas and
interwoven into it.

Will return a [`LazyCall`](@ref) object containing both an anonymous function
and its arguments. No more than one splatted positional argument can be woven
in. No more than one splatted keyword argument can be woven in provided there is
a `;` visible both inside and outside the tilda. Make multi-line functions by
wrapping in a begin block. To use `~` as a function, use the alias
[`bitnot`](@ref). If there are no woven arguments in `e`, return `e`.

# Examples

```julia
@test 1 == @unweave 1
@test 1 == @unweave +(1)

A = [1, 2]
B = ( [5, 6], [7, 8] )

unweave_test = @chain begin
    @unweave vcat(~A, ~[3, 4], ~(B...) )
    run(_, map)
end

@test unweave_test ==
      map((a, c, b...) -> vcat(a, c, b...), A, [3, 4], B...)

keyword_test(; keyword_arguments...) = keyword_arguments

a = keyword_test(a = 1, b = 2)

unweave_keyword_test = @chain begin
    @unweave keyword_test(c = 3; ~(a...))
    run
end

@test unweave_keyword_test == keyword_test(c = 3; a... )

# Can splat no more than one positional argument
@test_throws ErrorException unweave(:( ~(a...) + ~(b...) ))

# Can splat no more than one keyword argument
@test_throws ErrorException unweave(:( ~(;a...) + ~(;b...) ))
```
"""
function unweave(e::Expr)
    anonymous_function, arguments = split_anonymous(e)

    if length(arguments) == 0
        return(e)
    end

    @chain begin
        arguments
        Expr(:call, :collect_arguments, _...)
        Expr(:call, :LazyCall, _, anonymous_function)
    end
end

"""
    @unweave f e

[`unweave`](@ref) `e` then run `f` on the component parts.

If there are no woven arguments in `e`, return `e`.

# Examples
```julia
@test 1 == @unweave 1 1
@test 1 == @unweave 1 +(1)

e = :(vcat(~a, ~b) )
f = :broadcast
unweave(f, e)

a = [1, 2]
b = [3, 4]

@test broadcast((a, b) -> vcat(a, b), a, b) ==
    @unweave broadcast vcat(~a, ~b)
```
"""
function unweave(f, e::Expr)

    anonymous_function, anonymous_arguments = split_anonymous(e)

    if length(anonymous_arguments) == 0
        return(e)
    end

    Expr(:call, f, anonymous_function, anonymous_arguments...)

end

"""
    @unweave f::Expr e

[`unweave`](@ref) `e` then insert the function as the first argument to `f` and
the woven arguments at the end of the arguments of `f`.

If there are no woven arguments in `e`, return `e`.

# Examples
```julia
@test 1 == @unweave +(1) 1
@test 1 == @unweave +(1) +(1)

e = :(~a + ~b)
f = :(NullableArrays.broadcast(lift = true))

unweave(f, e)

a = NullableArrays.NullableArray([1, 2])
b = NullableArrays.NullableArray([3, 4], [false, true])

result = @unweave broadcast(lift = true) ~a + ~b

@test result.values[1] == 4
@test result.isnull == [false, true]

# `f` must be a call
@test_throws ErrorException unweave(:(import ChainMap), :(~_ + 1) )
```
"""
function unweave(f::Expr, e::Expr)

    function_test = MacroTools.@capture f function_call_(arguments__)

    if !(function_test)
        error("`f` must be a call")
    end

    anonymous_function, anonymous_arguments = split_anonymous(e)

    if length(anonymous_arguments) == 0
        return(e)
    end

    Expr(:call, function_call, anonymous_function,
         arguments..., anonymous_arguments...)

end

@nonstandard unweave
export @unweave

"""
    @broadcast e

A convenience macro for [`unweave`](@ref)` where `f` = `broadcast`

# Examples

```julia
a = [1, 2]
b = [3, 4]

@test broadcast((a, b) -> vcat(a, b), a, b) ==
    @broadcast vcat(~a, ~b)
```
"""
macro broadcast(e)
    esc(unweave(:broadcast, e))
end
export @broadcast

export bitnot
"""
    bitnot

Alias for `~` for use within [`@unweave`](@ref)

# Examples
```julia
@test bitnot(1) == ~1
```
"""
bitnot = ~


text:
# A. De Marchi, Feb 2021.

using DataFrames
using Printf
using Plots

"""
    time_profile(stats, cost)
Produce a time profile comparing solvers in `stats` using the `cost` function.
Inputs:
- `stats::Dict{Symbol,DataFrame}`: pairs of `:solver => df`;
- `cost::Function`: cost function applyed to each `df`. Should return a vector with the cost of solving the problem at each row;
  - 0 cost is not allowed;
  - If the solver did not solve the problem, return Inf, a negative number, or NaN.
Examples of cost functions:
- `cost(df) = df.elapsed_time`: Simple `elapsed_time` cost. Assumes the solver solved the problem.
- `cost(df) = (df.status .!= :first_order) * Inf + df.elapsed_time`: (default) Takes into consideration the status of the solver.
"""
function time_profile(
    stats::Dict{Symbol,DataFrame},
    cost::Function = df -> (df.status .!= :first_order) * Inf + df.elapsed_time,
    args...;
    kwargs...,
)
    solvers = keys(stats)
    dfs = (stats[s] for s in solvers)
    T = hcat([cost(df) for df in dfs]...)
    time_profile(T, string.(solvers), args...; kwargs...)
end

"""Produce a time profile.
Each column of the matrix `T` defines the cost for a solver (smaller is better).
Failures on a given problem are represented by a negative value, an infinite value, or `NaN`.
The optional argument `logscale` is used to produce a logarithmic (base 10) time plot.
"""
function time_profile(
    T::Array{Float64,2},
    labels::Vector{AbstractString};
    logscale::Bool = true,
    title::AbstractString = "",
    kwargs...,
)
    kwargs = Dict{Symbol,Any}(kwargs)
    # compute time profiles
    (fmat, tvec) = time_fractions(T)
    (nt, ns) = size(fmat)
    # setup style
    length(labels) == 0 && (labels = [@sprintf("column %d", col) for col = 1:ns])
    xlabel = pop!(kwargs, :xlabel, "Runtime" * (logscale ? " (log scale)" : ""))
    ylabel = pop!(kwargs, :ylabel, "Fraction of problems")
    linestyles = pop!(kwargs, :linestyles, Symbol[])
    # initial plot
    profile = Plots.plot(
        xlabel = xlabel,
        ylabel = ylabel,
        title = title,
        xlims = ((logscale && (tvec[1] == 0.0)) ? tvec[2] : tvec[1], tvec[nt]),
        ylims = (0.0, 1.0),
        xscale = (logscale ? :log10 : :none),
    )
    # plot solvers
    for s = 1:ns
        if length(linestyles) > 0
            kwargs[:linestyle] = linestyles[s]
        end
        Plots.plot!(tvec, fmat[:, s], t = :steppost, label = labels[s]; kwargs...)
    end
    return profile
end

time_profile(
    T::Array{Tn,2},
    labels::Vector{S};
    kwargs...,
) where {Tn<:Number,S<:AbstractString} = time_profile(
    convert(Array{Float64,2}, T),
    convert(Vector{AbstractString}, labels);
    kwargs...,
)

"""
    Compute time fractions used to produce a time profile.
    There is normally no need to call this function directly.
"""
function time_fractions(T::Array{Float64,2})

    # number of problems, number of solvers
    (np, ns) = size(T)

    # manage NaN, Inf, and negative values
    T[isinf.(T)] .= NaN
    T[T.<0] .= NaN
    failures = isnan.(T)

    # create time vector (x-axis)
    tvec = sort(unique(T[.!failures]))
    if tvec[1] == 0.0
        tvec = tvec[2:end]
    end
    append!(tvec, 1.1 * tvec[end])
    nt = length(tvec)

    # Compute fractions for each solver (y-axis)
    fmat = zeros(nt, ns)
    for is = 1:ns
        for it = 1:nt
            fmat[it, is] = sum(T[:, is] .≤ tvec[it]) / np
        end
    end

    return (fmat, tvec)
end

time_fractions(T::Array{Td,2}) where {Td<:Number} =
    time_fractions(convert(Array{Float64,2}, T))


text:
module DiskCollections
using Serialization
using Dates

export DiskDict, LoggingDict, checkpoint!, update_key

abstract type AbstractDiskDict{K,V} <: AbstractDict{K,V} end
mutable struct DiskDict{K,V} <: AbstractDiskDict{K,V}
    cache::Dict{K,V}
    s::Serializer
    io::IO
    n::Int
    counter::Int
    buffered_checkpoint::Bool
end

@enum ActionTag setaction deleteaction

function init_dict(K,V,io)
    seek(io, 0)
    s = Serializer(io)
    dict::Dict{K,V} = if eof(io)
        d = Dict{K,V}()
        serialize(s, d)
        d
    else
        deserialize(s)
    end
    s, dict
end

"""
    DiskDict{K,V}(io::IO, n::Int; buffered_checkpoint=true)
    DiskDict{K,V}(filename::String, n::Int; buffered_checkpoint=true)
    DiskDict{K,V}(f::Function, filename::String, n::Int; buffered_checkpoint=true)

Constructs a persistent dictionary behaving like a `Dict{K,V}` that serializes
keys and values to the seekable IO `io`. If a filename is provided instead that
file is opened with `read`, `write`, and `create` set to `true`. Every `n`
updates the dictionary reserializes its current state. If `buffered_checkpoint`
is `true`, this reserialization is performed in memory before writing the data
to the io.

Because the constructor opens a file that should be closed a function can be
provided as the first argument, enabling the use of Do-Block Syntax just as for
`open`.

```jldoctest
julia> new_disk_dict = DiskDict(open("mydict.jlp","w+"), 1000)
DiskDict{Any,Any}()

julia> new_disk_dict["name"] = "DiskCollections"
"DiskCollections"

julia> close(new_disk_dict)

julia> plain_dict = DiskDict("mydict.jlp", 1000) do disk_dict
           Dict(disk_dict)
       end
Dict{Any,Any} with 1 entry:
  "name" => "DiskCollections"
```
"""
function DiskDict{K,V}(io::IO, n::Int; buffered_checkpoint=true) where {K,V}
    s, cache = init_dict(K,V,io)
    counter = 0
    while !eof(io)
        tag = read(io, UInt8)
        if tag == UInt8(deleteaction)
            key = deserialize(s)
            delete!(cache, key)
        elseif tag == UInt8(setaction)
            key = deserialize(s)
            value = deserialize(s)
            cache[key] = value
        end
        counter += 1
    end
    DiskDict(cache, s, io, n, counter, buffered_checkpoint)
end

function DiskDict{K,V}(filename::String, n::Int; kwargs...) where {K,V}
    io = open(filename, create=true, read=true, write=true)
    DiskDict{K,V}(io, n; kwargs...)
end

function DiskDict{K,V}(f::Function, filename::String, args...; kwargs...) where {K,V}
    open(filename, create=true, read=true, write=true) do io
        f(DiskDict{K,V}(io, args...; kwargs...))
    end
end

DiskDict(args...; kwargs...) = DiskDict{Any,Any}(args...; kwargs...)

Base.close(dict::DiskDict) = close(dict.io)
Base.isopen(dict::DiskDict) = isopen(dict.io)

Base.length(dict::AbstractDiskDict) = length(dict.cache)
Base.iterate(dict::AbstractDiskDict, state...) = iterate(dict.cache, state...)
Base.getindex(dict::AbstractDiskDict{K,V}, key::K) where {K,V} = dict.cache[key]

function Base.setindex!(h::DiskDict{K,V}, v0, key0) where {K,V}
    key = convert(K, key0)
    if !isequal(key, key0)
        throw(ArgumentError("$(Base.limitrepr(key0)) is not a valid key for type $K"))
    end
    setindex!(h, v0, key)
end

function Base.setindex!(dict::DiskDict{K,V}, value, key::K) where {K,V}
    v = convert(V, value)
    if dict.counter >= dict.n
        checkpoint!(dict)
    end

    write_action(dict, setaction, key, v)
    dict.counter += 1

    dict.cache[key] = v
end

update_key(h::AbstractDiskDict, key) = setindex!(h, key, h[key])

function write_action(coll, action::ActionTag, args...)
    write(coll.io, UInt8(action))
    for a in args
        serialize(coll.s, a)
    end
    flush(coll.io)
end

function checkpoint!(dict::DiskDict)
    seek(dict.io, 0)
    if dict.buffered_checkpoint
        buf = IOBuffer()
        serialize(Serializer(buf), dict.cache)
        write(dict.io, take!(buf))
    else
        serialize(Serializer(dict.io), dict.cache)
    end
    truncate(dict.io, position(dict.io))
    dict.s = Serializer(dict.io)
    flush(dict.io)
    dict.counter = 0
end

function Base.delete!(dict::DiskDict, key)
    if dict.counter >= dict.n
        checkpoint!(dict)
    end
    write_action(dict, deleteaction, key)
    dict.counter += 1

    delete!(dict.cache, key)
end

function Base.pop!(dict::AbstractDiskDict, key)
    v = dict.cache[key]
    delete!(dict, key)
    return v
end

mutable struct LoggingDict{K,V} <: AbstractDiskDict{K,V}
    cache::Dict{K,V}
    s::Serializer
    io::IO
    log_start::Int
end

"""
    LoggingDict{K,V}(io::IO; until::Union{DateTime,Nothing}=nothing)
    LoggingDict{K,V}(filename::String; until::Union{DateTime,Nothing}=nothing)
    LoggingDict{K,V}(f::Function, filename::String; until::Union{DateTime,Nothing}=nothing)

Constructs a persistent dictionary behaving like a `Dict{K,V}` that serializes
keys and values to the seekable IO `io`. If a filename is provided instead that
file is opened with `read`, `write`, and `create` set to `true`. For each update
the `LoggingDict` stores the time; by providing a `DateTime` for the `until`
parameter of the `LoggingDict` constructor the updates are loaded from the file
up until that point in time.

Because the constructor opens a file that should be closed a function can be
provided as the first argument, enabling the use of Do-Block Syntax just as for
[`open`](@ref).

```jldoctest
julia> using Dates

julia> LoggingDict{Int,String}("mylogdict.jlp") do logdict
           for i in 1:5
               logdict[i] = string(i)
           end
           sleep(10)
           for i in 6:10
               logdict[i] = string(i)
           end
       end

julia> LoggingDict{Int,String}("mylogdict.jlp", until=now()-Second(10))
LoggingDict{Int64,String} with 5 entries:
  4 => "4"
  2 => "2"
  3 => "3"
  5 => "5"
  1 => "1"
```
"""
function LoggingDict{K,V}(io::IO; until::Union{DateTime,Nothing}=nothing) where {K,V}
    s, cache = init_dict(K,V,io)
    log_start = position(io)
    while !eof(io)
        datetime = unix2datetime(read(io, Float64))
        if !isnothing(until) && datetime > until
            break
        end
        tag = read(io, UInt8)
        if tag == UInt8(deleteaction)
            key = deserialize(s)
            delete!(cache, key)
        elseif tag == UInt8(setaction)
            key = deserialize(s)
            value = deserialize(s)
            cache[key] = value
        end
    end
    LoggingDict(cache, s, io, log_start)
end

function Base.delete!(dict::LoggingDict, key)
    write_action(dict, deleteaction, key)
    delete!(dict.cache, key)
end

function Base.setindex!(dict::LoggingDict{K,V}, value, key::K) where {K,V}
    v = convert(V, value)

    write(dict.io, datetime2unix(now()))
    write_action(dict, setaction, key, v)

    dict.cache[key] = v
end

function LoggingDict{K,V}(filename::String; kwargs...) where {K,V}
    io = open(filename, create=true, read=true, write=true)
    LoggingDict{K,V}(io; kwargs...)
end

function LoggingDict{K,V}(f::Function, filename::String, args...; kwargs...) where {K,V}
    open(filename, create=true, read=true, write=true) do io
        f(LoggingDict{K,V}(io, args...; kwargs...))
    end
end

LoggingDict(args...; kwargs...) = LoggingDict{Any,Any}(args...; kwargs...)

end


text:
"""
    l1ca_chipping_rate

L1 C/A Code chipping rate = `1.023e6Hz`
"""
const l1ca_chipping_rate = 1.023e6  # Hz

"""
    l1ca_chip_length

L1 C/A code chip length in 1ms = `1023`
"""
const l1ca_code_length = 1023  # chips

"""
    l1ca_db_chipping_rate

The chipping rate of the data bits for the L1 C/A code = `50Hz`
"""
const l1ca_db_chipping_rate = 50.  # Hz

"""
    l1ca_taps

Dictionary containing the phase selector taps (columns 1 and 2) and the code
chip delay amount in chips (column 3). The dictionary keys represent the PRNs.
"""
const l1ca_taps = Dict( 1 => (2,  6,   5),
                        2 => (3,  7,   6),
                        3 => (4,  8,   7),
                        4 => (5,  9,   8),
                        5 => (1,  9,  17),
                        6 => (2, 10,  18),
                        7 => (1,  8, 139),
                        8 => (2,  9, 140),
                        9 => (3, 10, 141),
                       10 => (2,  3, 251),
                       11 => (3,  4, 252),
                       12 => (5,  6, 254),
                       13 => (6,  7, 255),
                       14 => (7,  8, 256),
                       15 => (8,  9, 257),
                       16 => (9, 10, 258),
                       17 => (1,  4, 469),
                       18 => (2,  5, 470),
                       19 => (3,  6, 471),
                       20 => (4,  7, 472),
                       21 => (5,  8, 473),
                       22 => (6,  9, 474),
                       23 => (1,  3, 509),
                       24 => (4,  6, 512),
                       25 => (5,  7, 513),
                       26 => (6,  8, 514),
                       27 => (7,  9, 515),
                       28 => (8, 10, 516),
                       29 => (1,  6, 859),
                       30 => (2,  7, 860),
                       31 => (3,  8, 861),
                       32 => (4,  9, 862))

"""
    ca_code_chips(prn::Int)

Generates the CA Code chip-by-chip sequence for a given GPS satellite PRN
number.

Required Arguments:

- `prn::Int`: the pseudorandom noise (PRN) number, which is also related to the
              satellite vehicle

Returns:

- `ca_code::Vector{Int}`: a 1023-element array of C/A code chips for the PRN
                          number argument, `prn`
"""
function ca_code_chips(prn::Int)
    # Initialize G1 and G2 registers to 1's
    G1 = ones(Int64, 10)
    G2 = ones(Int64, 10)
    # Get phase selector taps for specific PRN
    phase_sel_taps = l1ca_taps[prn][1:2]
    ca_code = Array{Int64}(undef, l1ca_code_length)
    @inbounds for chip in 1:l1ca_code_length
        G1_output = G1[10]
        G1_epoch = xor(G1_output, G1[3])
        # Calculate phase selector output
        phase_sel_out = xor(G2[phase_sel_taps[1]],
                            G2[phase_sel_taps[2]])
        # XOR phase selector and G1 outputs to yield
        # C/A code value
        ca_code[chip] = xor(G1_output, phase_sel_out)
        # Calculate G2 epoch by XORing the G2 taps
        # 2, 3, 6, 8, 9, and 10
        G2_epoch = xor(G2[2], G2[3], G2[6], G2[8],
                       G2[9], G2[10])
        # Shift registers forward by one chip
        @inbounds for i in 10:-1:2
            G1[i] = G1[i-1]
            G2[i] = G2[i-1]
        end
        # Save G1 and G2 epoch values to beginning of
        # shifted G1 and G2 arrays
        G1[1] = G1_epoch
        G2[1] = G2_epoch
    end
    return ca_code
end

"""
    gen_l1ca_codes()

Generates a dictionary containing all the L1 C/A codes for PRNS 1 through 32.
The PRN number is the dictionary key.

Arguments: `None`

Returns:

- `l1ca_codes::Dict`: dictionary containing a 1023-element array of C/A code
                      chips for each PRN
"""
function gen_l1ca_codes()
    l1ca_codes = Dict{Int64, Array{Int64,1}}()
    for prn in keys(l1ca_taps)
        l1ca_codes[prn] = ca_code_chips(prn)
    end
    return l1ca_codes
end

"""
    l1ca_codes

Constant dictionary containing the L1 C/A codes. The keys are the PRN numbers.
"""
const l1ca_codes = gen_l1ca_codes()

"""
    define_l1ca_code_type(t_length=missing; prns=collect(keys(l1ca_codes)),
                          sig_freq=L1_freq)

Generates a `SignalType` struct which defines the L1 C/A code that can be used
to define and generate L1 C/A signals. The databits generated are random.

Optional Arguments:

- `t_length:Float64`: the length of the planned signal in seconds
    * this argument is used only to generate a `Vector` of databits so they do
      not repeat throughout the duration of the signal, once it is generated
    * `(default = missing)`
    * if given, databits are included in signal type
- `prns`: codes to define for signal type `(default = keys(l1ca_codes))`
	* default is to define all PRN codes and assume that databits for each
	  are different
- `sig_freq`: carrier frequency in Hz `(default = L1_freq)`
- `databits`: optional databits `(default = missing)`, where random bits are used
- `B`: receiver bandwidth `(default = missing)`

Returns:

- `siganl_type::SignalType`: can be used in `definesignal` function to create
                             a signal, which can be used with `generatesignal!`
                             to generate the signal
"""
function define_l1ca_code_type(t_length=missing; prns=collect(keys(l1ca_codes)),
                               sig_freq=L1_freq, databits=missing, B=missing)
    if ~ismissing(t_length)
        # Generate random databits
        if ismissing(databits)
            databits = random_databits(l1ca_db_chipping_rate, t_length; prns=prns)
        end
        # Include databits in I channel codes
        I_codes = definecodetype(copy_dictionary(l1ca_codes, prns),
                                 l1ca_chipping_rate;
                                 databits=[databits, l1ca_db_chipping_rate])
    else
        # Define I channel codes without databits
        I_codes = definecodetype(copy_dictionary(l1ca_codes, prns),
                                 l1ca_chipping_rate)
    end
    signal_type = definesignaltype(I_codes, sig_freq, "I"; name="L1 C/A", B=B)
    return signal_type
end


text:
module ReactiveMPModelsHMMTest

using Test, InteractiveUtils
using Rocket, ReactiveMP, GraphPPL, Distributions
using BenchmarkTools, Random, Plots, Dates, LinearAlgebra, StableRNGs

## Model definition
## -------------------------------------------- ##
@model [ default_factorisation = MeanField() ] function transition_model(n)

    A ~ MatrixDirichlet(ones(3, 3))
    B ~ MatrixDirichlet([ 10.0 1.0 1.0; 1.0 10.0 1.0; 1.0 1.0 10.0 ])

    s_0 ~ Categorical(fill(1.0 / 3.0, 3))

    s = randomvar(n)
    x = datavar(Vector{Float64}, n)

    s_prev = s_0

    for t in 1:n
        s[t] ~ Transition(s_prev, A) where { q = q(out, in)q(a) }
        x[t] ~ Transition(s[t], B)
        s_prev = s[t]
    end

    return s, x, A, B
end
## -------------------------------------------- ##
## Inference definition
## -------------------------------------------- ##
function inference(data, vmp_iters)
    n = length(data)

    model, (s, x, A, B) = transition_model(n, options = (limit_stack_depth = 500, ))

    sbuffer = keep(Vector{Marginal})
    Abuffer = keep(Marginal)
    Bbuffer = keep(Marginal)
    fe      = ScoreActor(Float64)

    ssub  = subscribe!(getmarginals(s), sbuffer)
    Asub  = subscribe!(getmarginal(A), Abuffer)
    Bsub  = subscribe!(getmarginal(B), Bbuffer)
    fesub = subscribe!(score(Float64, BetheFreeEnergy(), model), fe)

    setmarginal!(A, vague(MatrixDirichlet, 3, 3))
    setmarginal!(B, vague(MatrixDirichlet, 3, 3))

    foreach(s) do svar
        setmarginal!(svar, vague(Categorical, 3))
    end

    for i in 1:vmp_iters
        update!(x, data)
    end

    unsubscribe!(ssub)
    unsubscribe!(Asub)
    unsubscribe!(Bsub)
    unsubscribe!(fesub)

    return map(getvalues, (sbuffer, Abuffer, Bbuffer, fe))
end

@testset "Hidden Markov Model" begin

    @testset "Full graph" begin
        ## -------------------------------------------- ##
        ## Data creation
        ## -------------------------------------------- ##
        function rand_vec(rng, distribution::Categorical)
            k = ncategories(distribution)
            s = zeros(k)
            s[ rand(rng, distribution) ] = 1.0
            s
        end

        function generate_data(rng, n_samples)

            # Transition probabilities (some transitions are impossible)
            A = [0.9 0.0 0.1; 0.1 0.9 0.0; 0.0 0.1 0.9]
            # Observation noise
            B = [0.9 0.05 0.05; 0.05 0.9 0.05; 0.05 0.05 0.9]
            # Initial state
            s_0 = [1.0, 0.0, 0.0]
            # Generate some data
            s = Vector{Vector{Float64}}(undef, n_samples) # one-hot encoding of the states
            x = Vector{Vector{Float64}}(undef, n_samples) # one-hot encoding of the observations

            s_prev = s_0

            for t = 1:n_samples
                a = A * s_prev
                s[t] = rand_vec(rng, Categorical(a ./ sum(a)))
                b = B * s[t]
                x[t] = rand_vec(rng, Categorical(b ./ sum(b)))
                s_prev = s[t]
            end

            return x, s
        end
        rng = StableRNG(123)
        x_data, s_data = generate_data(rng, 100);
        ## -------------------------------------------- ##
        ## Inference execution
        sbuffer, Abuffer, Bbuffer, fe = inference(x_data, 20);
        ## -------------------------------------------- ##
        ## Test inference results
        @test length(sbuffer) === 20
        @test length(Abuffer) === 20
        @test length(Bbuffer) === 20
        @test length(fe) === 20 && all(filter(e -> abs(e) > 1e-3, diff(fe)) .< 0)
        @test abs(last(fe) - 60.614480654) < 0.01

        ## -------------------------------------------- ##
        ## Form debug output
        base_output = joinpath(pwd(), "_output", "models")
        mkpath(base_output)
        timestamp        = Dates.format(now(), "dd-mm-yyyy-HH-MM")
        plot_output      = joinpath(base_output, "hmm_plot_$(timestamp)_v$(VERSION).png")
        benchmark_output = joinpath(base_output, "hmm_benchmark_$(timestamp)_v$(VERSION).txt")
        ## -------------------------------------------- ##
        ## Create output plots
        p1 = scatter(argmax.(s_data))
        p2 = scatter(argmax.(ReactiveMP.probvec.(sbuffer[end])))
        p3 = plot(fe)
        p = plot(p1, p2, p3, layout = @layout([ a b; c ]))
        savefig(p, plot_output)
        ## -------------------------------------------- ##
        ## Create output benchmarks
        benchmark = @benchmark inference($x_data, 20) seconds=15
        open(benchmark_output, "w") do io
            show(io, MIME("text/plain"), benchmark)
            versioninfo(io)
        end
        ## -------------------------------------------- ##
    end

end

end


text:
using StokesParticles, Test

d = 1.3
ϵ = 0.1
s = 0.3

@test ( StokesParticles.LJForceMagnitude((1+2*ϵ)*d,s,d,ϵ) == 0)
@test ( StokesParticles.LJForceMagnitude(d,s,d,ϵ) == 0)
@test ( StokesParticles.LJForceMagnitude((1-ϵ)*d,s,d,ϵ) < 0 )
@test ( StokesParticles.LJForceMagnitude((1+ϵ/2)*d,s,d,ϵ) > 0 )