text:
"""
    drawpointorbitR2(f, x0, y0 [; preiterations, iterations])

Return the drawing of the orbit of a point \$(x_0,y_0)\$ under a function
    \$f:\\mathbb{R}^2\\rightarrow\\mathbb{R}^2\$,
    in a rectangular region in \$\\mathbb{R}^2\$.

The orbit of \$(x_0,y_0)\$ under \$f\$ is defined as
    \$o(x_0,y_0,f)=\\{(x_0,y_0),f(x_0,y_0),\\dots,f^n(x_0,y_0),\\dots\\}\$

#### Arguments
- `f::Function`: A function \$f:\\mathbb{R}^2\\rightarrow\\mathbb{R}^2\$.
- `x0::Real`: X coordinate of the point.
- `y0::Real`: Y coordinate of the point.
- `preiterations::Integer`: Number of first iterations to calculate but not to draw.
- `iterations::Integer`: Number of iterations to calculate (after `preiterations`).
"""
function drawpointorbitR2(f::Function, x0::Real, y0::Real;
    preiterations::Integer=0, iterations::Integer=100)

    # Verifying functions
    @assert typeof(f(1., 1.)) <: Tuple{Real,Real}

    SDDGraphics.newdrawing()
    SDDGraphics.updatecolorarray(iterations)

    xn,yn = x0,y0

    for n in 1:preiterations
        xn,yn = f(xn,yn)
    end # for n preiterations

    for n in 1:iterations
        if SDDGraphics.insiderectregion(xn,yn)
            SDDGraphics.color(n)
            SDDGraphics.drawpoint(xn,yn)
        end # if
        xn,yn = f(xn,yn)
    end # for n iterations

    SDDGraphics.drawing()

end # function drawpointorbitR2

"""
    drawpointssetorbitR2(f, ps [; preiterations, iterations])

Return the drawing of the orbit of a points set under a function
    \$f:\\mathbb{R}^2\\rightarrow\\mathbb{R}^2\$,
    in a rectangular region in \$\\mathbb{R}^2\$.

The orbit of the set \$P\$ under \$f\$ is defined as
    \$o(P,f)=\\bigcup_{n\\geq 0}f^n(P)\$.

#### Arguments
- `f::Function`: A function \$f:\\mathbb{R}^2\\rightarrow\\mathbb{R}^2\$.
- `ps::Array{Tuple{Real,Real},1}`: The set of points (array of pairs of real numbers).
- `preiterations::Integer`: Number of first iterations to calculate but not to draw.
- `iterations::Integer`: Number of iterations to calculate (after `preiterations`).
"""
function drawpointssetorbitR2(f::Function, ps::Array{Tuple{T,T},1};
    preiterations::Integer=0, iterations::Integer=100) where T <: Real

    # Verifying functions
    @assert typeof(f(1., 1.)) <: Tuple{Real,Real}

    SDDGraphics.newdrawing()
    SDDGraphics.updatecolorarray(iterations)

    pns = deepcopy(ps)

    for n in 1:preiterations
        for k in 1:length(pns)
            pns[k] = f(pns[k]...)
        end # for k pns
    end # for n preiterations

    for n in 1:iterations
        SDDGraphics.color(n)
        for k in 1:length(pns)
            if SDDGraphics.insiderectregion(pns[k]...)
                SDDGraphics.drawpoint(pns[k]...)
            end # if
            pns[k] = f(pns[k]...)
        end # for k pns
    end # for n iterations

    SDDGraphics.drawing()

end # function drawpointorbitR2

"""
    drawpointorbitC(f, z0 [; preiterations, iterations])

Return the drawing of the orbit of a point \$z_0\$ under a function
    \$f:\\mathbb{C}\\rightarrow\\mathbb{C}\$,
    in a rectangular region in \$\\mathbb{C}\$.

The orbit of \$z_0\$ under \$f\$ is defined as
    \$o(z_0,f)=\\{z_0,f(z_0),\\dots,f^n(z_0),\\dots\\}\$.

#### Arguments
- `f::Function`: A function \$f:\\mathbb{C}\\rightarrow\\mathbb{C}\$.
- `z0::Number`: The point.
- `preiterations::Integer`: Number of first iterations to calculate but not to draw.
- `iterations::Integer`: Number of iterations to calculate (after `preiterations`).
"""
function drawpointorbitC(f::Function, z0::Number;
    preiterations::Integer=0, iterations::Integer=100)

    # Verifying functions
    @assert typeof(f(1.0im)) <: Number

    SDDGraphics.newdrawing()
    SDDGraphics.updatecolorarray(iterations)

    zn = z0

    for n in 1:preiterations
        zn = f(zn)
    end # for n preiterations

    for n in 1:iterations
        if SDDGraphics.insiderectregion(zn)
            SDDGraphics.color(n)
            SDDGraphics.drawpoint(zn)
        end #if
        zn = f(zn)
    end # for n iterations

    SDDGraphics.drawing()

end # function drawpointorbitC

"""
    drawpointssetorbitC(f, ps [; preiterations, iterations])

Return the drawing of the orbit of a points set under a function
    \$f:\\mathbb{C}\\rightarrow\\mathbb{C}\$,
    in a rectangular region in \$\\mathbb{C}\$.

The orbit of the set \$P\$ under \$f\$ is defined as
    \$o(P,f)=\\bigcup_{n\\geq 0}f^n(P)\$.

#### Arguments
- `f::Function`: A function \$f:\\mathbb{C}\\rightarrow\\mathbb{C}\$.
- `ps::Array{Number,1}`: The set of points (array of numbers).
- `preiterations::Integer`: Number of first iterations to calculate but not to draw.
- `iterations::Integer`: Number of iterations to calculate (after `preiterations`).
"""
function drawpointssetorbitC(f::Function, ps::Array{T,1};
    preiterations::Integer=0, iterations::Integer=100) where T <: Number

    # Verifying functions
    @assert typeof(f(1.0im)) <: Number

    SDDGraphics.newdrawing()
    SDDGraphics.updatecolorarray(iterations)

    pns = deepcopy(ps)

    for n in 1:preiterations
        pns = f.(pns)
    end # for n preiterations

    for n in 1:iterations
        SDDGraphics.color(n)
        for k in 1:length(pns)
            if SDDGraphics.insiderectregion(pns[k])
                SDDGraphics.drawpoint(pns[k])
            end # if
            pns[k] = f(pns[k])
        end # for k pns
    end # for n iterations

    SDDGraphics.drawing()

end # function drawpointorbitR2

"""
    drawpointorbitpathR2(f, x0, y0 [; preiterations, iterations, kind])

Return the drawing of the orbit's path of a point \$(x_0,y_0)\$ under a function
\$f:\\mathbb{R}^2\\rightarrow\\mathbb{R}^2\$,
in a rectangular region in \$\\mathbb{R}^2\$.

The orbit of \$(x_0,y_0)\$ under \$f\$ is defined as
\$o(x_0,y_0,f)=\\{(x_0,y_0),f(x_0,y_0),\\dots,f^n(x_0,y_0),\\dots\\}\$.

#### Arguments
- `f::Function`: A function \$f:\\mathbb{R}^2\\rightarrow\\mathbb{R}^2\$.
- `x0::Real`: X coordinate of the point.
- `y0::Real`: Y coordinate of the point.
- `preiterations::Integer`: Number of first iterations to calculate but not to draw.
- `iterations::Integer`: Number of iterations to calculate (after `preiterations`).
- `kind::Symbol`: Kind of path: `:line` or `:arc`.
"""
function drawpointorbitpathR2(f::Function, x0::Real, y0::Real;
    preiterations::Integer=0, iterations::Integer=100, kind::Symbol=:line)

    drawsegment = SDDGraphics.drawlinesegment

    if kind == :arc
        SDDGraphics.supported(:drawarc)
        drawsegment = SDDGraphics.drawarc
    else
        SDDGraphics.supported(:drawlinesegment)
    end

    # Verifying functions
    @assert typeof(f(1., 1.)) <: Tuple{Real,Real}

    SDDGraphics.newdrawing()
    SDDGraphics.updatecolorarray(iterations)

    xn0,yn0 = x0,y0

    for n in 1:preiterations
        xn0,yn0 = f(xn0,yn0)
    end # for n preiterations

    for n in 1:iterations
        xn1,yn1 = f(xn0,yn0)
        SDDGraphics.color(n)
        drawsegment(xn0,yn0,xn1,yn1)
        xn0,yn0 = xn1,yn1
    end # for n iterations

    SDDGraphics.drawing()

end # function drawpointorbitpathR2

"""
    drawpointorbitpathC(f, z0 [; preiterations, iterations, kind])

Return the drawing of the orbit's path of a point \$z_0\$ under a function
\$f:\\mathbb{C}\\rightarrow\\mathbb{C}\$,
in a rectangular region in \$\\mathbb{C}\$.

The orbit of \$z_0\$ under \$f\$ is defined as
\$o(z_0,f)=\\{z_0,f(z_0),\\dots,f^n(z_0),\\dots\\}\$.

#### Arguments
- `f::Function`: A function \$f:\\mathbb{C}\\rightarrow\\mathbb{C}\$.
- `z0::Number`: The point.
- `preiterations::Integer`: Number of first iterations to calculate but not to draw.
- `iterations::Integer`: Number of iterations to calculate (after `preiterations`).
- `kind::Symbol`: Kind of path: `:line` or `:arc`.
"""
function drawpointorbitpathC(f::Function, z0::Number;
    preiterations::Integer=0, iterations::Integer=100, kind::Symbol=:line)

    drawsegment = SDDGraphics.drawlinesegment

    if kind == :arc
        SDDGraphics.supported(:drawarc)
        drawsegment = SDDGraphics.drawarc
    else
        SDDGraphics.supported(:drawlinesegment)
    end

    # Verifying functions
    @assert typeof(f(1.0im)) <: Number

    SDDGraphics.newdrawing()
    SDDGraphics.updatecolorarray(iterations)

    zn0 = z0

    for n in 1:preiterations
        zn0 = f(zn0)
    end # for n preiterations

    for n in 1:iterations
        zn1 = f(zn0)
        SDDGraphics.color(n)
        drawsegment(zn0,zn1)
        zn0 = zn1
    end # for n iterations

    SDDGraphics.drawing()

end # function drawpointorbitpathC

"""
    drawpointorbitpathR(f, x0, y0 [; preiterations, iterations])

Return the drawing of the orbit's arc-path of a point \$x_0\$ under a function
\$f:\\mathbb{R}\\rightarrow\\mathbb{R}\$.

The orbit of \$x_0\$ under \$f\$ is defined as
\$o(x_0,f)=\\{x_0,f(x_0,\\dots,f^n(x_0,\\dots\\}\$.

#### Arguments
- `f::Function`: A function \$f:\\mathbb{R}^2\\rightarrow\\mathbb{R}^2\$.
- `x0::Real`: X coordinate of the point.
- `y0::Real`: Y coordinate of the point.
- `preiterations::Integer`: Number of first iterations to calculate but not to draw.
- `iterations::Integer`: Number of iterations to calculate (after `preiterations`).
"""
function drawpointorbitpathR(f::Function, x0::Real;
    preiterations::Integer=0, iterations::Integer=100)

    SDDGraphics.supported(:drawarc)

    # Verifying functions
    @assert typeof(f(1.)) <: Real

    SDDGraphics.newdrawing()
    SDDGraphics.updatecolorarray(iterations)

    xn0 = x0

    for n in 1:preiterations
        xn0 = f(xn0)
    end # for n preiterations

    for n in 1:iterations
        xn1 = f(xn0)
        SDDGraphics.color(n)
        SDDGraphics.drawarc(xn0,0,xn1,0)
        xn0 = xn1
    end # for n iterations

    SDDGraphics.drawing()

end # function drawpointorbitpathR2


text:
module MangalTestKnownBugs

    using Mangal
    using EcologicalNetworks
    using Test

    _problematic_ids = [28, 83, 86, 100, 101, 102, 103, 1280, 2463, 2464, 3258, 922, 924, 889, 890, 907, 906, 1516]

    for (idx,netid) in enumerate(_problematic_ids)
        thisnet = Mangal.network(netid)
        @test typeof(thisnet) <: MangalNetwork
        U = convert(UnipartiteNetwork, thisnet)
        @test typeof(U) <: UnipartiteNetwork
    end

end

ext: jl


text:
"""
    ExplicitRungeKuttaSolver <: AbstractRungeKuttaSolver

A composite type for an explicit [`AbstractRungeKuttaSolver`](@ref).

# Constructors
```julia
ExplicitRungeKuttaSolver(tableau, h[, adaptive])
ERK(args...; kwargs...)
```

# Arguments
- `tableau  :: ButcherTableau`     : Butcher tableau.
- `stepsize :: StepSize`           : step-size.
- `adaptive :: AdaptiveParameters` : embedded method's parameters.

# Functions
- [`show`   ](@ref) : shows name and contents.
- [`summary`](@ref) : shows name.
"""
mutable struct ExplicitRungeKuttaSolver{tableau_T, stepsize_T, adaptive_T} <: AbstractRungeKuttaSolver
    tableau::tableau_T
    stepsize::stepsize_T
    adaptive::adaptive_T
end

ExplicitRungeKuttaSolver(tableau, h::Real, adaptive) = ExplicitRungeKuttaSolver(tableau, StepSize(h), adaptive)
ExplicitRungeKuttaSolver(tableau, stepsize) = ExplicitRungeKuttaSolver(tableau, stepsize, nothing)
@doc (@doc ExplicitRungeKuttaSolver) ERK(args...; kwargs...) = ExplicitRungeKuttaSolver(args...; kwargs...)

############################################################################################
#                                         PRINTING                                         #
############################################################################################

"""
    show(io::IO, solver::ExplicitRungeKuttaSolver)

prints a full description of `solver` and its contents to a stream `io`.
"""
Base.show(io::IO, solver::ExplicitRungeKuttaSolver) = NSDEBase._show(io, solver)

"""
    summary(io::IO, solver::ExplicitRungeKuttaSolver)

prints a brief description of `solver` to a stream `io`.
"""
Base.summary(io::IO, solver::ExplicitRungeKuttaSolver) = NSDEBase._summary(io, solver)

############################################################################################
#                                          METHODS                                         #
############################################################################################

function (solver::ExplicitRungeKuttaSolver)(problem::AbstractInitialValueProblem; save_stages::Bool = false)
    solve(problem, solver; save_stages=save_stages)
end


text:
#=
function offline_tracking(wt,max_whiskers=10)

    for i=1:size(wt.vid,3)
        image_preprocessing(wt.vid,i)
    end

    for i=1:size(wt.vid,3)
        #Need to transpose becuase row major in C vs column major in julia
        wt.all_whiskers[i]=WT_trace(i,wt.vid[:,:,i]',wt.min_length,wt.pad_pos,wt.mask)
        println(string(i,"/",size(wt.vid,3)))
    end

    reorder_whiskers(wt)

    nothing
end

function offline_tracking_parallel(wt,max_whiskers=10)

    pmap(t->image_preprocessing(wt.vid,t),1:size(wt.vid,3),batch_size=ceil(Int,size(wt.vid,3)/nworkers()))
    #pmap(t->image_preprocessing(wt.vid,t),1:10)
    println("Preprocessing complete")

    wt.all_whiskers=pmap(t->WT_trace(t,wt.vid[:,:,t]',wt.min_length,wt.pad_pos,wt.mask),1:size(wt.vid,3),
    batch_size=ceil(Int,size(wt.vid,3)/nworkers()))
    #wt.all_whiskers=pmap(t->WT_trace(t,wt.vid[:,:,t]',wt.min_length,wt.pad_pos,wt.mask),1:10)

    reorder_whiskers(wt)

    nothing
end
=#
function image_preprocessing(vid,i)

    temp_img = zeros(Float64,size(vid,1),size(vid,2))
    temp_img2=zeros(Float64,size(vid,1),size(vid,2))

    #Adjust contrast

    for j=1:size(vid,1)
        for k=1:size(vid,2)
            temp_img[j,k] = convert(Float64,vid[j,k,i])
        end
    end

    local_contrast_enhance!(temp_img,temp_img)

    #anisotropic diffusion to smooth while perserving edges
    anisodiff!(temp_img, 20,20.0,0.05,1,temp_img2)

    for j=1:size(vid,1)
        for k=1:size(vid,2)
            vid[j,k,i] = round(UInt8,temp_img[j,k])
        end
    end

    nothing
end

function mean_image_uint8(han::Tracker_Handles)
    #round.(UInt8,squeeze(mean(han.wt.vid,3),3)) #These should use the normalized from teh whole video
end

function mean_image(han::Tracker_Handles)
    #squeeze(mean(han.wt.vid,3),3) #These should use the normalized from the whole video
end


text:
@recipe function f(fSTL::STL)

    # subplots configuration
    legend:= false
    grid:= false
    layout := @layout [Data
                       Seasonal
                       Trend
                       Remainder]

    xtickfont:= font(3, "Courier")
    ytickfont:= font(3, "Courier")
    yguidefont:= font(5, "Courier")

    # reference bar
    a,b = extrema(skipmissing(values(fSTL.decomposition[!,:Seasonal].+
                                     fSTL.decomposition[!,:Trend].+
                                     fSTL.decomposition[!,:Remainder]))); hd = b-a
    a,b = extrema(skipmissing(values(fSTL.decomposition[!,:Seasonal]))); hs = b-a
    a,b = extrema(skipmissing(values(fSTL.decomposition[!,:Trend]))); ht = b-a
    a,b = extrema(skipmissing(values(fSTL.decomposition[!,:Remainder]))); hr = b-a
    mh = min(hd,hs,ht,hr)

    inset_subplots := [(1, bbox(-0.012, 0, 0.01, 1.0, :left)),
                       (2, bbox(-0.012, 0, 0.01, 1.0, :right)),
                       (3, bbox(-0.012, 0, 0.01, 1.0, :left)),
                       (4, bbox(-0.012, 0, 0.01, 1.0, :right))]

    @series begin
        subplot := 1
        yguide := "Data"
        xaxis := nothing
        bottom_margin := -5Plots.mm
        fSTL.decomposition[!,:Seasonal] .+
            fSTL.decomposition[!,:Trend] .+
            fSTL.decomposition[!,:Remainder]
    end

    @series begin
        subplot := 2
        yguide := "Trend"
        xaxis := nothing
        ymirror:=true
        guide_position:=:left
        bottom_margin := -5Plots.mm
        fSTL.decomposition[!,:Trend]
    end

    @series begin
        subplot := 3
        yguide := "Seasonal"
        xaxis := nothing
        seriestype := :sticks
        bottom_margin := -5Plots.mm
        fSTL.decomposition[!,:Seasonal]
    end

    @series begin
        subplot := 4
        yguide := "Remainder"
        seriestype := :sticks
        #bottom_margin:=0Plots.mm
        ymirror:=true
        guide_position:=:left

        fSTL.decomposition[!,:Timestamp], fSTL.decomposition[!,:Remainder]
    end

    @series begin
        subplot := 5
        background_color_inside := nothing
        framestyle := :none
        seriestype := :bar
        ylims := (0,hd)
        mh:mh
    end

    @series begin
        subplot := 6
        background_color_inside := nothing
        framestyle := :none
        seriestype := :bar
        ylims := (0,ht)
        mh:mh
    end

    @series begin
        subplot := 7
        framestyle := :none
        seriestype := :bar
        ylims := (0,hs)
        mh:mh
    end

    @series begin
        subplot := 8
        framestyle := :none
        seriestype := :bar
        ylims := (0,hr)
        mh:mh
    end

end


text:
export ZabbixHistory, get_history

struct ZabbixHistory
  id::Int
  name::String
end
ZabbixHistory(dict) =
  ZabbixHistory(
    parse(Int, dict["historyid"]),
    dict["name"]
  )

showres(h::ZabbixHistory; offset=0) =
  println(repeat(" ", offset), "ZabbixHistory: ", h.id, " => ", h.name)

function get_history(z, item::ZabbixItem)
  history = ZabbixHistory[]
  params = Dict("itemids"=>item.id)
  _history = @checkerr make_request(z, "history.get", params)
  for h in _history
    push!(history, ZabbixHistory(h))
  end
  history
end


text:
mutable struct ADM{U, O} <: AbstractSubspaceOptimiser
    λ::U
    R::O
end

get_threshold(opt::ADM) = opt.λ
set_threshold!(opt::ADM, λ) = (opt.λ = λ; return)

"""
    ADM(λ = 0.1)

Optimiser for finding a sparse basis vector in a subspace based on [this paper](https://arxiv.org/pdf/1412.4659.pdf).
`λ` is the weight for the soft-thresholding operation.
"""
function ADM(λ::T = 0.1) where T <: Real
    f = NormL1(λ)
    return ADM(λ, f)
end

function fit!(q::AbstractArray{T, 1}, Y::AbstractArray, opt::ADM; maxiter::Int64= 10) where T <: Real
    normalize!(q)
    x = Y*q
    for k in 1:maxiter
        prox!(x, opt.R, Y*q)
        mul!(q, Y', x/norm(Y'*x, 2))
    end

    q[abs.(q) .< get_threshold(opt)] .= zero(eltype(q))
    normalize!(q, 2)
    return
end

function fit!(q::AbstractArray{T, 2}, Y::AbstractArray, opt::ADM; maxiter::Int64= 10) where T <: Real
    @inbounds for i in 1:size(q, 2)
        fit!(view(q, :, i), Y, opt, maxiter = maxiter)
    end

    return
end


text:
using JuMP, Ipopt
model = Model(with_optimizer(Ipopt.Optimizer));

@variable(model, x, start = 4)
@variable(model, y, start = -9.66);

@NLparameter(model, p == 0.003); # Providing a starting value is necessary for parameters
@NLparameter(model, l[i = 1:10] == 4 - i); # A collection of parameters

value(l[1])

set_value(l[1], -4)
value(l[1])

@NLexpression(model, expr_1, sin(x))
@NLexpression(model, expr_2, asin(expr_1)); # Inserting one expression into another

@NLconstraint(model, exp(x) + y^4 <= 0)
@NLobjective(model, Min, tan(x) + log(y))

my_function(a,b) = (a * b)^-6 + (b / a)^3
register(model, :my_function, 2, my_function, autodiff = true)

using Random, Statistics

n = 1_000
data = randn(n)

mle = Model(with_optimizer(Ipopt.Optimizer, print_level = 0))
@NLparameter(mle, problem_data[i = 1:n] == data[i])
@variable(mle, μ, start = 0.0)
@variable(mle, σ >= 0.0, start = 1.0)
@NLexpression(mle, likelihood,
(2 * π * σ^2)^(-n / 2) * exp(-(sum((problem_data[i] - μ)^2 for i in 1:n) / (2 * σ^2)))
)

@NLobjective(mle, Max, log(likelihood))

optimize!(mle)

println("μ = ", value(μ))
println("mean(data) = ", mean(data))
println("σ^2 = ", value(σ)^2)
println("var(data) = ", var(data))
println("MLE objective: ", objective_value(mle))

# Changing the data

data = randn(n)
optimize!(mle)

println("μ = ", value(μ))
println("mean(data) = ", mean(data))
println("σ^2 = ", value(σ)^2)
println("var(data) = ", var(data))
println("MLE objective: ", objective_value(mle))


text:
module Models

export
    NonhydrostaticModel, ShallowWaterModel,
    HydrostaticFreeSurfaceModel, VectorInvariant,
    ExplicitFreeSurface, ImplicitFreeSurface, SplitExplicitFreeSurface,
    HydrostaticSphericalCoriolis, VectorInvariantEnstrophyConserving,
    PrescribedVelocityFields, PressureField

using Oceananigans: AbstractModel

import Oceananigans.Architectures: device_event

device_event(model::AbstractModel) = device_event(model.architecture)

abstract type AbstractNonhydrostaticModel{TS} <: AbstractModel{TS} end

include("NonhydrostaticModels/NonhydrostaticModels.jl")
include("HydrostaticFreeSurfaceModels/HydrostaticFreeSurfaceModels.jl")
include("ShallowWaterModels/ShallowWaterModels.jl")

using .NonhydrostaticModels: NonhydrostaticModel, PressureField

using .HydrostaticFreeSurfaceModels:
    HydrostaticFreeSurfaceModel, VectorInvariant,
    ExplicitFreeSurface, ImplicitFreeSurface, SplitExplicitFreeSurface,
    HydrostaticSphericalCoriolis,
    VectorInvariantEnstrophyConserving, VectorInvariantEnergyConserving,
    PrescribedVelocityFields

using .ShallowWaterModels: ShallowWaterModel

end


text:
# # [Implicit Nonlinear Dynamics : Michaelis Menten](@id michaelis_menten)
#
# What if you want to estimate an implicitly defined system of the form ``f(u_t, u, p, t) = 0``?
# The solution : Implicit Sparse Identification. This method was originally described in [this paper](http://ieeexplore.ieee.org/document/7809160/), and currently there exist [robust algorithms](https://royalsocietypublishing.org/doi/10.1098/rspa.2020.0279) to identify these systems.
# We will focus on [Michaelis Menten Kinetics](https://en.wikipedia.org/wiki/Michaelis%E2%80%93Menten_kinetics). As before, we will define the [`DataDrivenProblem`](@ref) and the [`Basis`](@ref) containing possible candidate functions for our [`sparse_regression!`](@ref).
# Lets generate some data! We will use two experiments starting from different initial conditions.

using DataDrivenDiffEq
using LinearAlgebra
using ModelingToolkit
using OrdinaryDiffEq
#md using Plots

function michaelis_menten(u, p, t)
    [0.6 - 1.5u[1]/(0.3+u[1])]
end

u0 = [0.5]

problem_1 = ODEProblem(michaelis_menten, u0, (0.0, 4.0));
solution_1 = solve(problem_1, Tsit5(), saveat = 0.1);
problem_2 = ODEProblem(michaelis_menten, 2*u0, (4.0, 8.0));
solution_2 = solve(problem_2, Tsit5(), saveat = 0.1);

# Since we have multiple trajectories at hand, we define a [`DataDrivenDataset`](@ref), which collects multiple problems but handles them as a unit
# for the processing.

function michaelis_menten(X::AbstractMatrix, p, t::AbstractVector)
    reduce(hcat, map((x,ti)->michaelis_menten(x, p, ti), eachcol(X), t))
end

data = (
    Experiment_1 = (X = Array(solution_1), t = solution_1.t, DX = michaelis_menten(Array(solution_1),[], solution_1.t) ),
    Experiment_2 = (X = Array(solution_2), t = solution_2.t, DX = michaelis_menten(Array(solution_2),[], solution_2.t))
)

prob = DataDrivenDiffEq.ContinuousDataset(data);
#md plot(prob)

# Next, we define our [`Basis`](@ref). Since we want to identify an implicit system, we have to include
# some candidate terms which use these as an argument and inform our constructor about the meaning of these variables.

@parameters t
D = Differential(t)
@variables u[1:1](t)
h = [monomial_basis(u[1:1], 4)...]
basis = Basis([h; h .* D(u[1])], u, implicits = D.(u), iv = t)
println(basis) # hide

# Next, we define the [`ImplicitOptimizer`](@ref) and `solve` the problem. It wraps a standard optimizer, by default [`STLSQ`](@ref), and performs
# implicit sparse regression upon the selected basis.
# To improve our result, we batch the data by using a [`DataSampler`](@ref). Here, we use a train-test split of 0.8 and
# divide the training data into 10 batches. Since we are using a batching process, we can also use a different

sampler = DataSampler(
    Split(ratio = 0.8), Batcher(n = 10)
)

opt = ImplicitOptimizer(1e-1:1e-1:5e-1)
res = solve(prob, basis, opt,  normalize = false, denoise = false, by = :min, sampler = sampler, maxiter = 1000);
println(res) # hide

# As we can see, the [`DataDrivenSolution`](@ref) has good metrics. Furthermore, inspection of the underlying system shows that the original equations have been recovered correctly:

system = result(res)
println(system) # hide

#md plot(
#md     plot(prob), plot(res), layout = (1,2)
#md )

#md # ## [Copy-Pasteable Code](@id michaelis_menten_copy_paste)
#md #
#md # ```julia
#md # @__CODE__
#md # ```

## Test #src
for r_ in [res] #src
    @test all(l2error(r_) .< 0.4) #src
    @test all(aic(r_) .> 5e2) #src
    @test all(determination(r_) .>= 0.7) #src
end #src


text:
println("\nTesting MCDTS complete tree, coupled Logistic, CCM:")
@time begin
@testset "MCDTS single rollout on univariate Logistic map data" begin

    Lval = 3500
    x = zeros(Lval)
    y = zeros(Lval)
    r = 3.8
    r2 = 3.5
    βxy = 0.02
    βyx = 0.1
    x[1]=0.4
    y[1]=0.2

    for i = 2:Lval
        x[i]=x[i-1]*(r-r*x[i-1]-βxy*y[i-1])
        y[i]=y[i-1]*(r2-r2*y[i-1]-βyx*x[i-1])
    end

    w1 = DelayEmbeddings.estimate_delay(x, "mi_min")
    w2 = DelayEmbeddings.estimate_delay(y, "mi_min")

    test1 = DelayEmbeddings.standardize(x)
    test2 = DelayEmbeddings.standardize(y)
    # try MCDTS with CCM
    taus1 = 0:10 # the possible delay vals
    trials = 20 # the sampling of the tree

    Random.seed!(1234)
    optmodel = MCDTS.MCDTSOptimGoal(MCDTS.CCM_ρ(test2), MCDTS.Range_function())
    tree = mcdts_embedding(Dataset(test1), optmodel, w1, taus1, trials)
    best_node = MCDTS.best_embedding(tree)
    τ_mcdts = best_node.τs
    ts_mcdts = best_node.ts
    Lval = L(best_node)

    # less fid points
    Random.seed!(1234)
    optmodel2 = MCDTS.MCDTSOptimGoal(MCDTS.CCM_ρ(test2,1,0.5), MCDTS.Range_function())
    tree2 = mcdts_embedding(Dataset(test1), optmodel2, w1, taus1, trials)
    best_node2 = MCDTS.best_embedding(tree2)
    τ_mcdts2 = best_node2.τs
    ts_mcdts2 = best_node2.ts
    L2 = L(best_node2)

    @test Lval < -0.95
    @test Lval - .01 < L2 < Lval + .01
    @test length(τ_mcdts) ==  3 == length(τ_mcdts2)
    @test τ_mcdts[1] == 0 == τ_mcdts2[1]
    @test τ_mcdts[2] == 2 == τ_mcdts2[2]
    @test τ_mcdts[3] == 1 == τ_mcdts2[3]

    Random.seed!(1234)
    optmodel2 = MCDTS.MCDTSOptimGoal(MCDTS.CCM_ρ(test1), MCDTS.Range_function())
    tree2 = mcdts_embedding(Dataset(test2), optmodel2, w2, taus1, trials)
    best_node = MCDTS.best_embedding(tree2)
    τ_mcdts = best_node.τs
    ts_mcdts = best_node.ts
    Lval = L(best_node)

    @test Lval < -0.997
    @test length(τ_mcdts) ==  3
    @test τ_mcdts[1] == 0
    @test τ_mcdts[2] == 1
    @test τ_mcdts[3] == 2

end
end
true


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

using Markdown
using InteractiveUtils

# ╔═╡ 3a8ef900-d5d9-11eb-1e22-edc33d208d0e
import Pkg; Pkg.add(path="c:/Users/4vektor/.julia/dev/SSpline")

# ╔═╡ 1e7150a7-40ce-44e7-a719-fa995c10ac86
let
   # import Pkg
	Pkg.activate(".")
    #Pkg.add("SSpline")
end

# ╔═╡ 1ff97a4d-92bf-4f4f-ab60-c5f311b6b2b4
using Plots, SSpline

# ╔═╡ fcc81fdb-a8c0-489c-b36e-4065a91dca6e
pyplot()

# ╔═╡ c56bd9d8-6829-47c8-b610-e0bb0c1f88ed
html"""<style> main { max-width: 900px;} """

# ╔═╡ 101c094a-b0eb-46fb-a06d-1415de4018da
begin
	x=[10,20,35,50,55,57,60];      y=[5,30,20,22,18,30,12];
	#x=[0,10,30,50,70,90,100];     y=[30,130,150,150,170,220,320];
	#x=[0.0,25.0,50.0,75.0,100.0]; y=[20.0,50.0,10.0,60.0,35.0];
	#x=[0,10,30,45,70,90,100];     y=[20,50,70,-70,60,50,20];
	#x=[0,10,30,50,70,80,82];      y=[150,200,200,200,180,100,0];
end

# ╔═╡ ee8e5634-b8c3-4a99-8dad-2ee228221609
begin
	ra = x[1]:0.2:x[end] # interpolation range
end

# ╔═╡ a8eb6119-8ba8-4a60-9be4-7787731da630
begin
	nn = spline0(x,y)
	ls = spline1(x,y)
	ni = interp(nn, ra)
	li = interp(ls, ra)
end

# ╔═╡ 2cf7565e-6d87-4219-bd8b-64a0a3e7ec41
begin
	plot(size = (850,500), legend = :bottom, title = "nearest neighbor and linear interpolations")
	scatter!(x, y, label = "knots")
	plot!(ni[1], ni[2], color = :green, label = "nearest neighbor")
	plot!(li[1], li[2], color = :blue, style = :dashdot, label = "linear")
end

# ╔═╡ 8081af65-98bc-4160-8b73-915f671e21a9
begin
	qs = spline2(x,y)
	qs2 = spline2(x,y, 6.1)
	qi = interp(qs, ra)
	qi2 = interp(qs2, ra)
	qsl = slope(qs2, ra)
end

# ╔═╡ 4f01a157-9a1a-4bfc-b367-df142313e7e6
begin
	plot(size = (850,500), legend = :topleft, title = "linear and quadratic interpolations")
	scatter!(x, y, label = "knots")
	plot!(li[1], li[2], color = :grey, label = "linear")
	plot!(qi[1], qi[2], color = :orange, style = :dashdot, label = "quadratic, f'(x₁) = 0.0 (default)")
	plot!(qi2[1], qi2[2], color = :blue, thickness = 2.0, style = :dash, label = "quadratic, f'(x₁) = 6.1")
end

# ╔═╡ 0bcea8ae-002b-4fb6-be6a-89a4b1abadb1
begin
	csn = spline3(x,y, :natural)
	csc = spline3(x,y, :clamped, -5.0, 5.0)
	csp = spline3(x,y,:notaknot)
	ci = interp(csn,ra)
	cic = interp(csc, ra)
	csl = slope(csn, ra)
	cci = interp(csp,ra)
end

# ╔═╡ f72de351-b2ae-4a96-b6ef-d8627a45c70a
begin
	plot(size = (850,500), legend = :topleft, title = "quadratic vs. cubic natural")
	scatter!(x, y, label = "knots", color = :white, markersize = 10, markershape = :hexagon)
	plot!(qi2[1], qi2[2], color = :green, style = :dashdot, thickness = 2.0, label = "quadratic, slope 6.1 at start")
	plot!(ci[1], ci[2], color = :purple, style = :dash, thickness = 2.0, label = "cubic natural")
	plot!(cic[1], cic[2], color = :orange2, thickness = 2.0, label = "cubic clamped (-5.0, 5.0)")
	plot!(cci[1], cci[2], color = :skyblue, label = "notaknot f")
end

# ╔═╡ 03e94833-ab83-46c8-9eb9-4612727d072d
begin
	xh=[0,10,30,50,70,90,100]
	yh=[30,130,150,150,170,220,320]

	rah = xh[1]:0.2:xh[end]

	hk = hspline(xh, yh, :kruger)
	hki = interp(hk, rah)

	hf = hspline(xh, yh, :free)
	hfi = interp(hf, rah)

	#hg = hspline(xh, yh, :geom)
	#hgi = interp(hg, rah)
end

# ╔═╡ 16be7b7b-c784-4519-8ec8-13a8c1ac2d26
begin
	plot(size = (850,500), legend = :bottom, title = "hermite splines")
	scatter!(xh, yh, label = "knots", color = :white, markersize = 10, markershape = :hexagon)
	plot!(hki[1], hki[2], color = :orange2, label = "hermite kruger (harmonic mean)")
	#plot!(hgi[1], hgi[2], color = :blue, style = :dashdot, label = "hermite geometric mean")
	plot!(hfi[1], hfi[2], color = :green, style = :dash, label = "hermite free")
end

# ╔═╡ 5cefbebb-fe49-4ef1-b3fb-5e6577341383
begin
	cel = extrapolate(hk, -10:110, :linear)
	ces = extrapolate(hk, -10:110, :smooth)
	ceil = interp(cel, -10:0.2:110)
	ceis = interp(ces, -10:0.2:110)
end

# ╔═╡ df66e2da-d6a0-4afe-8fc6-dc36398a5bc2
begin
	plot(size = (850,400), legend = :bottom, title = "hermite extrapolations")
	scatter!(xh, yh, label = "knots", color = :white, markersize = 10, markershape = :hexagon)
	plot!(hki[1], hki[2], color = :orange, width = 2.0, label = "hermite kruger")
	plot!(ceil[1], ceil[2], color = :skyblue, style = :dash, label = "herm. kruger linear extr.")
	plot!(ceis[1], ceis[2], color = :green, style = :dashdot, label = "herm. kruger smooth extr.")
end

# ╔═╡ 7fc4737e-9b75-4d53-bae1-f977a1466d8e
begin
	cnn = spline2(xh,yh)
	ceq = extrapolate(cnn, -10:110, :zero)
	ceiq = interp(ceq, -10:0.2:110)
	plot(size = (850,400), legend = :bottom, title = "cubic extrapolations")
	scatter!(xh, yh, label = "knots", color = :white, markersize = 10, markershape = :hexagon)
	plot!(ceiq[1], ceiq[2], color = :red, style = :dot, label = "cub. nat. lin.")
	scatter!(ceq.x, ceq.a)
end

# ╔═╡ a5fb28d3-4b1b-4bda-8d20-1b5796f1374f
ceiq[2][end]

# ╔═╡ 680e2580-f4bf-47d1-bca8-19086f017b04
ceq

# ╔═╡ d7189e79-63d5-46bd-bf8c-d37f5ed8d50c
ceq.a[end-1]

# ╔═╡ 014045fc-e2b2-48cf-aeb6-0e71cbe18426
begin
	xx=[0,10,30,50,70,90,100];	yy=[30,130,150,150,170,220,320]
	#xx=[-30,-15, -10,-7, -5, 0, 5, 7, 10, 15, 30]; yy=[0, 20, 15, 20, 8, 8, 8, 20, 15, 20, 0];
	#xx=[10,18,20,35,50,55,57,60]; yy=[0,20,20,20,22,18,30,12];
	#xx=[10,18,20,35,50,55,57,60]; yy=[12,30,18,22,20,20,20,0];
	#x2=[10,13,15,20,35,50,52,60]; y2=[0,20,20,20,22,18,30,12];

	cn = cspline(xx,yy,:natural)
	cp = cspline(xx,yy,:parabolic)
	cc = cspline(xx,yy,:notaknot)
	cr = nakspline(xx,yy)
	#natural
	nx,ny = interp(cn, xx[1]:.1:xx[end])
	n1x,n1y = deriv1(cn, xx[1]:1:xx[end])
	n2x,n2y = deriv2(cn, xx[1]:1:xx[end])
	n3x,n3y = deriv3(cn, xx[1]:1:xx[end])
	nkx,nky = deriv2(cn, xx)
	#parabolic
	px,py = interp(cp, xx[1]:.1:xx[end])
	p1x,p1y = deriv1(cp, xx[1]:1:xx[end])
	p2x,p2y = deriv2(cp, xx[1]:1:xx[end])
	p3x,p3y = deriv3(cp, xx[1]:1:xx[end])
	pkx,pky = deriv2(cp, xx)
	#notaknot
	cx,cy = interp(cc, xx[1]:.1:xx[end])
	c1x,c1y = deriv1(cc, xx[1]:1:xx[end])
	c2x,c2y = deriv2(cc, xx[1]:1:xx[end])
	c3x,c3y = deriv3(cc, xx[1]:1:xx[end])
	ckx,cky = deriv2(cc, xx)

	#rx,ry = interp(cr, x2[1]:0.1:x2[end])

	plot(size = (950,600), xlimits = (xx[1], xx[end]), legend = :outertopright, title = "extrapolations", fg_color_grid = :black, gridalpha =  0.5)
	scatter!(xx, yy, label = "knots", color = :white, markersize = 10, markershape = :hexagon)
	#plot!(e,f, color = :blue)
	#plot!(hx,hy, color = :red, label = "natural")
	plot!(nx,ny, color = :blue, style = :solid, label = "natural f")
	#plot!(n1x,10*n1y, color = :blue, style = :dash, label = "natural f'")
	#plot!(n2x,10*n2y, color = :blue, style = :dashdot, label = "natural f''")
	#plot!(n3x,10*n3y, color = :blue, style = :dot, label = "natural f'''")
	#scatter!(nkx, nky, color = :blue)
	#plot!(px,py, color = :red, style = :dash, label = "parabolic f")
	#plot!(p1x,10*p1y, color = :red, style = :solid, label = "parabolic f'")
	#plot!(p2x,10*p2y, color = :red, style = :dashdot, label = "parabolic f''")
	plot!(p3x,100*p3y, color = :red, style = :dot, label = "parabolic f'''")
	#scatter!(pkx, pky, color = :red)
	plot!(cx, cy,color = :orange, style = :dash, label = "notaknot f")
	#plot!(c1x,10*c1y, color = :green, style = :solid, label = "notaknot f'")
	#plot!(c2x,100*c2y, color = :orange, style = :dashdot, label = "notaknot f''")
	plot!(c3x,100*c3y, color = :orange, style = :dot, label = "notaknot f'''")
	#scatter!(nkx, nky, color = :green)
	#plot!(rx, ry,color = :green, style = :dash, label = "notaknot reverse")
end

# ╔═╡ b2889364-f087-434a-a137-9ca6dd85bb19

# ╔═╡ 2ef87041-69f6-4e37-8a46-22e480a2b417
begin
        n = length(x)-1
        Δx= Vector{Float64}(undef, n-1)
        Δy= Vector{Float64}(undef, n-1)
        λ = Vector{Float64}(undef, n)    # upper diagonal
        δ = ones(n) * 2 # main diagonal
        μ = Vector{Float64}(undef, n)    # lower diagonal
        r = Vector{Float64}(undef, n)    # derivatives
        M = Vector{Float64}(undef, n)

        @inbounds for i in 1:n-1
            Δx[i] = x[i+1] - x[i]
            Δy[i] = y[i+1] - y[i]
        end
            A = zeros(n, n)
            @inbounds for i in 2:n-2
                A[i, i-1] = Δx[i-1]
                A[i, i] = 2.0 * (Δx[i-1] + Δx[i])
                #A[i, i] = 2.0 * (Δx[i-1] + Δx[i])
                A[i, i+1] = Δx[i+1]
            end
            #A[1, 1] = -Δx[2]
            #A[1, 2] = Δx[1] + Δx[2]
            #A[1, 3] = -Δx[1]
            A[begin, begin] = -Δx[begin + 2]
            A[begin, begin + 1] = Δx[begin + 1] + Δx[begin + 2]
            A[begin, begin + 2] = -Δx[begin + 1]

            A[end, end-2] = -Δx[end - 1]
            A[end, end-1] = Δx[end - 2] + Δx[end - 1]
            A[end, end]   = -Δx[end - 2]
            for i in 2:n-1
            #    r[i] = 6.0 / (Δx[i-1] + Δx[i]) * (Δy[i] / Δx[i] - Δy[i-1] / Δx[i-1])
                r[i] = 6.0 * (Δy[i] / Δx[i] - Δy[i-1] / Δx[i-1])
            end
end

# ╔═╡ b4a2a3c5-7f10-4330-b7ad-adbb54addf01
A

# ╔═╡ Cell order:
# ╠═3a8ef900-d5d9-11eb-1e22-edc33d208d0e
# ╠═1e7150a7-40ce-44e7-a719-fa995c10ac86
# ╠═1ff97a4d-92bf-4f4f-ab60-c5f311b6b2b4
# ╠═fcc81fdb-a8c0-489c-b36e-4065a91dca6e
# ╠═c56bd9d8-6829-47c8-b610-e0bb0c1f88ed
# ╠═101c094a-b0eb-46fb-a06d-1415de4018da
# ╠═ee8e5634-b8c3-4a99-8dad-2ee228221609
# ╠═a8eb6119-8ba8-4a60-9be4-7787731da630
# ╠═2cf7565e-6d87-4219-bd8b-64a0a3e7ec41
# ╠═8081af65-98bc-4160-8b73-915f671e21a9
# ╠═4f01a157-9a1a-4bfc-b367-df142313e7e6
# ╠═0bcea8ae-002b-4fb6-be6a-89a4b1abadb1
# ╠═f72de351-b2ae-4a96-b6ef-d8627a45c70a
# ╠═03e94833-ab83-46c8-9eb9-4612727d072d
# ╠═16be7b7b-c784-4519-8ec8-13a8c1ac2d26
# ╠═5cefbebb-fe49-4ef1-b3fb-5e6577341383
# ╠═df66e2da-d6a0-4afe-8fc6-dc36398a5bc2
# ╠═7fc4737e-9b75-4d53-bae1-f977a1466d8e
# ╠═a5fb28d3-4b1b-4bda-8d20-1b5796f1374f
# ╠═680e2580-f4bf-47d1-bca8-19086f017b04
# ╠═d7189e79-63d5-46bd-bf8c-d37f5ed8d50c
# ╠═014045fc-e2b2-48cf-aeb6-0e71cbe18426
# ╠═b2889364-f087-434a-a137-9ca6dd85bb19
# ╠═2ef87041-69f6-4e37-8a46-22e480a2b417
# ╠═b4a2a3c5-7f10-4330-b7ad-adbb54addf01

max_stars_repo_path: pluto/Revolutionary science.jl

max_issues_repo_path: pluto/Revolutionary science.jl

max_forks_repo_path: pluto/Revolutionary science.jl


text:
using Test
using QuakeML

"For each field in `fields`, test that x.field is `missing`"
function test_all_missing(x, fields...)
    for f in fields
        getfield(x, f) === missing || (@show f)
        @test getfield(x, f) === missing
    end
end

"For each field in `fields`, test that x.field is empty"
function test_all_empty(x, fields...)
    for f in fields
        @test isempty(getfield(x, f))
    end
end

@testset "All tests" begin
    include("util.jl")
    include("construction.jl")
    include("comparison.jl")
    include("parsing.jl")
    include("io.jl")
    include("accessors.jl")
    include("field_access.jl")
    include("types_module.jl")
end

max_stars_repo_name: anowacki/QuakeML


text:
include("building_tree.jl")
include("utilities.jl")
include("merge.jl")

function main_merge()
    for dataSetName in ["iris", "seeds", "wine"]

        print("=== Dataset ", dataSetName)

        # Préparation des données
        include("../data/" * dataSetName * ".txt")
        train, test = train_test_indexes(length(Y))
        X_train = X[train,:]
        Y_train = Y[train]
        X_test = X[test,:]
        Y_test = Y[test]

        println(" (train size ", size(X_train, 1), ", test size ", size(X_test, 1), ", ", size(X_train, 2), ", features count: ", size(X_train, 2), ")")

        # Temps limite de la méthode de résolution en secondes
        time_limit = 10

        for D in 2:4
            println("\tD = ", D)
            println("\t\tUnivarié")
            testMerge(X_train, Y_train, X_test, Y_test, D, time_limit = time_limit, isMultivariate = false)
            println("\t\tMultivarié")
            testMerge(X_train, Y_train, X_test, Y_test, D, time_limit = time_limit, isMultivariate = true)
        end
    end
end

function testMerge(X_train, Y_train, X_test, Y_test, D; time_limit::Int=-1, isMultivariate::Bool = false)

    # Pour tout pourcentage de regroupement considéré
    println("\t\t\tGamma\t\t# clusters\tGap")
    for gamma in 0:0.2:1
        print("\t\t\t", gamma * 100, "%\t\t")
        clusters = simpleMerge(X_train, Y_train, gamma)
        print(length(clusters), " clusters\t")
        T, obj, resolution_time, gap = build_tree(clusters, D, multivariate = isMultivariate, time_limit = time_limit)
        print(round(gap, digits = 1), "%\t")
        print("Erreurs train/test : ", prediction_errors(T,X_train,Y_train))
        print("/", prediction_errors(T,X_test,Y_test), "\t")
        println(round(resolution_time, digits=1), "s")
    end
    println()
end

max_stars_repo_name: boumpteryx/RODM

max_issues_repo_name: boumpteryx/RODM

max_forks_repo_name: boumpteryx/RODM


text:
# ------------------------------------------------------------------
# Licensed under the MIT License. See LICENSE in the project root.
# ------------------------------------------------------------------

"""
    LCV(ranges, nfolds=10)

Likelihood cross-validation with parameter ranges
`ranges` and number of folds `nfolds`.

## References

* Sugiyama et al. 2008. Direct importance estimation for
  covariate shift adaptation.

### Author

* Júlio Hoffimann (julio.hoffimann@gmail.com)
"""
struct LCV <: EstimatorFitter
  ranges::NamedTuple
  nfolds::Int
end

LCV(ranges::NamedTuple) = LCV(ranges, 10)

function _fit(::Type{<:KLIEP}, x_nu, x_de,
              fitter::EstimatorFitter,
              optlib::Type{<:OptimizationLibrary})
  # retrieve parameters
  ranges = fitter.ranges
  nfolds = fitter.nfolds
  npts   = length(x_nu)

  @assert nfolds ≤ npts "number of folds must be smaller than number of numerator samples"

  # partition numerator samples into folds
  folds = collect(Iterators.partition(1:npts, npts ÷ nfolds))

  # initialize loss and optimal estimator
  Ĵmax, dre⭐ = -Inf, KLIEP()

  # loop over hyperparameters
  for σ in ranges.σ, b in ranges.b
    # density ratio estimator
    dre = KLIEP(σ=σ, b=b)

    # estimate loss with cross-validation
    Ĵₖ = map(1:nfolds) do k
      # training and hold-out samples
      train = [ind for i in vcat(1:k-1, k+1:nfolds) for ind in folds[i]]
      hold  = folds[k]

      # perform estimation with training samples
      r = densratiofunc(x_nu[train], x_de, dre, optlib=optlib)

      # evaluate loss with hold-out samples
      mean(log(r(x_nu[j])) for j in hold)
    end

    # mean loss for hyperparameters
    Ĵ = mean(Ĵₖ)

    # update and continue
    if Ĵ > Ĵmax
      Ĵmax = Ĵ
      dre⭐ = dre
    end
  end

  # optimal estimator
  dre⭐
end


text:
using MultivariateStats: PCA, pcacov
import IterTools: takenth, product
export PCAEmbedding

"""
	PCAEmbedding(s, em::AbstractSpatialEmbedding; kwargs...) → embedding
A spatio temporal delay coordinates structure with
Principal Component Analysis as a means of dimension reduction,
`embedding` can be used as a functor:
```julia
embedding(rvec, s, t, α)
```
which operates inplace on `rvec` and reconstructs vector from spatial time series `s` at
timestep `t` and cartesian index `α`.

To instantiate this `embedding`, give the data to be reconstructed `s` as well as an
instance of [`SpatioTemporalEmbedding`](@ref) to `PCAEmbedding`.

## Keyword Arguments
* `pratio = 0.99` : Ratio of variances that needs to be preserved in low-dimensional
  PCA-reconstruction.
* `maxoutdim = 25`: Upper limit for output dimension. May break `pratio` criterion.
* `every_t = 1` : Speed up computation by only using every n-th point in time.
* `every_α = 1` : Speed up computation further by only using every n-th point in space
  (linear indexing).

To set the output dimension to a certain value `X`, pass `pratio=1, maxoutdim=X`.
"""
struct PCAEmbedding{T,Φ,BC,X,Y} <: AbstractSpatialEmbedding{Φ,BC,X}
	stem::AbstractSpatialEmbedding{Φ,BC,Y}
	meanv::T
	covmat::Matrix{T}
	drmodel::PCA{T}
	tmp::Vector{T} #length Y
end

function PCAEmbedding(
		s::AbstractArray{<:AbstractArray{T,Φ}},
		stem::AbstractSpatialEmbedding{Φ,BC,Y};
		pratio   = 0.99,
		maxoutdim= 25,
		every_t::Int    = 1,
		every_α::Int    = 1) where {T,Φ,BC,Y}
	meanv = mean(mean.(s))
	tsteps = (length(s) - get_τmax(stem))
	usable_idxs = get_usable_idxs(stem)
	num_pt = length(usable_idxs)
	subset = Iterators.product(takenth(usable_idxs, every_α), 1:every_t:tsteps)
	L      = length(subset)

	recv = zeros(T,Y)
	covmat = zeros(T,Y,Y)
	for (α, t) ∈ subset
		stem(recv, s,t,α)
		recv  .-= meanv
		@inbounds for j=1:Y, i=1:Y
			covmat[i,j] += recv[i]*recv[j] / L
		end
	end
	drmodel = compute_pca(covmat,pratio, maxoutdim)
	X = size(drmodel.proj,2)
	tmp = zeros(Y)
	return PCAEmbedding{T,Φ,BC,X,Y}(stem, meanv, covmat, drmodel,tmp)
end

function (r::PCAEmbedding{T,Φ,BC,X,Y})(data, s, t, α) where {T,Φ,BC,X,Y}
	r.stem(r.tmp, s,t,α)
	r.tmp .-= r.meanv
	mul!(data,transpose(r.drmodel.proj),r.tmp)
end

Base.show(io::IO, em::PCAEmbedding) = (show(io, em.stem); show(io, em.drmodel))

compute_pca(covmat::Matrix{T}, pratio, maxoutdim) where T=
	pcacov(covmat, T[]; maxoutdim=maxoutdim, pratio=pratio)

get_τmax(em::PCAEmbedding) = get_τmax(em.stem)::Int #somehow type inference doesn't figure this out
get_num_pt(em::PCAEmbedding) = get_num_pt(em.stem)::Int

outdim(em::PCAEmbedding{T,Φ,BC,X}) where {T,Φ,BC,X} = X


text:
module GeoIP

using ZipFile
using Geodesy
using CSV
using IPNets
import Sockets: IPv4
import Base: getindex

export
    # types
    Location,
    # methods
    geolocate,
    setlocale,
    load

include("data.jl")
include("geoip-module.jl")

end # module


text:
const E = [
    0 0 1 1 0 0 0 0 0 0 0 0 0 0 1 1 0 0;
    0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0;
    0 1 0 0 0 1 0 0 0 0 0 0 1 0 0 0 1 0;
    0 1 0 0 1 0 1 0 0 0 0 1 0 1 0 0 1 0;
    0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0;
    1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1;
    0 1 0 0 0 0 1 0 0 0 0 1 0 0 0 0 1 0;
    0 1 0 0 0 0 0 1 0 0 1 0 0 0 0 0 1 0;
    0 1 0 0 0 0 0 0 1 1 0 0 0 0 0 0 1 0;
    0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0;
    0 0 1 1 0 0 0 0 0 0 0 0 0 0 1 1 0 0;
]

function frown!(E)
    E[7, 7] = E[7, 12] = E[9, 9] = E[9, 10] = 0
    E[9, 7] = E[9, 12] = E[7, 9] = E[7, 10] = 1
    E
end

frown(E) = frown!(copy(E))

rot270(E) = E'

# TODO This should probably be written out since rotr90 is unlike to be the solution a student would come up with
rot90(E) = rotr90(E)

stickerwall(E) = vcat(
    hcat(E, zeros(Int, size(E, 1)), ones(Int, size(E, 1)), zeros(Int, size(E, 1)), vcat(frown(E))),
    hcat(zeros(Int, size(E, 2) + 1)', [1], zeros(Int, size(E, 2) + 1)'),
    ones(Int, 2size(E, 2) + 3)',
    hcat(zeros(Int, size(E, 2) + 1)', [1], zeros(Int, size(E, 2) + 1)'),
    hcat(vcat(frown(E)), zeros(Int, size(E, 1)), ones(Int, size(E, 1)), zeros(Int, size(E, 1)), E)
)


text:
__precompile__()

module PhyloPlots

# standard libraries
using Markdown
using Printf: @printf, @sprintf

# dependencies that need explicit declaration
using DataFrames
using ColorTypes # used by Gadfly already. To resolve data type names (Colorant)
using Gadfly
import Gadfly.plot
using RCall
import RCall.sexp
using PhyloNetworks

export plot
export sexp

include("phylonetworksPlots.jl")
include("plotGadfly.jl")
include("plotRCall.jl")
include("substitutionmodels.jl")
include("rexport.jl")

end # of module


text:
module SimplePartitions
using DataStructures, Permutations

import Base.show, Base.==, Base.join, Base.+, Base.*, Base.in
import Base.<, Base.<=, Base.>, Base.>=

import AbstractLattices: ∧, ∨

export Partition, set_element_type, num_elements, num_parts, parts
export ground_set, merge_parts!, PartitionBuilder
export in_same_part, find_part, meet, refines, ∧, ∨

"""
`set_element_type(A)` gives the element type of a set `A`.
"""
function set_element_type(A::Set{T}) where T
  return T
end
set_element_type(B::BitSet) = Int

"""
A `Partition` is a set of nonempty, pairwise disjoint sets.
A new `Partition` is created by specifying the ground set `A`
and calling `Partition(A)`. The set `A` may be either a `Set{T}`
for some type `T` or an `BitSet`.

The parameter `A` may also be a list (one-dimensional array).

In addition, `Partition(n)` for a nonnegative integer `n` creates
a partition of the set {1,2,...,n}.

The datatype `Partition` is, essentially, a wrapper around the
`DataStructures.DisjointSets` type.
"""
mutable struct Partition{T}
  elements::Set{T}
  parts::DisjointSets{T}
  SOS::Set{Set{T}}  # set of sets representation

  function Partition{T}() where T
    A = Set{T}()
    return Partition(A)
  end

  function Partition{T}(A::Set{T}) where T
    elts = deepcopy(A)
    pts  = DisjointSets{T}(elts)
    sos = Set{Set{T}}()
    for a in A
      S = Set{T}([a])
      push!(sos,S)
    end
    new(elts,pts,sos)
  end
end # end type definition

function Partition(A::Set)
  T = set_element_type(A)
  return Partition{T}(A)
end
Partition(B::BitSet) = Partition{Int}(Set{Int}(B))

# Also construct from a vector
Partition(list::Vector{T}) where T = Partition(Set(list))
Partition(n::Int) = Partition(Set(1:n))
Partition() = Partition{Any}()

"""
`Partition(p::Permutation)` constructs a new partition in which the
parts are the cycles of the permutation `p`.
"""
function Partition(p::Permutation)
  n = length(p)
  P = Partition(n)
  clist = cycles(p)
  for C in clist
    merge_parts!(P,C)
  end
  return P
end

"""
`Partition(d::Dict)` constructs a new partition
whose elements are the keys of `d`. Two elements
`a` and `b` are in the same part of the partition
iff `d[a]==d[b]`.
"""
function Partition(d::Dict)
  K = collect(keys(d))
  P = Partition(K)
  n = length(K)
  for i=1:n-1
    a = K[i]
    for j=i+1:n
      b = K[j]
      if d[a]==d[b]
        merge_parts!(P,a,b)
      end
    end
  end
  return P
end

"""
`kill_sos!(P)` wipes out the set-of-sets cache. Only used
internally.
"""
kill_sos!(P::Partition{T}) where T = P.SOS = Set{Set{T}}()

# Construct a Partition from a set of sets
"""
`PartitionBuilder(A,check=true)` takes a set of nonempty, pairwise disjoint
sets and creates the corresponding partition. It is the inverse operation
to `parts(P)`. The optional parameter `check` causes sanity checks to be return
on the input set of sets (throwing errors if it is invalid).
"""
function PartitionBuilder(A::Set{Set{T}}, check::Bool=true) where T
  parts_list = collect(A)
  np = length(parts_list)
  parts_sum = 0

  if check
    for p in parts_list
      np = length(p)
      @assert np>0 "The sets in the input must be nonempty."
      parts_sum += np
    end
  end

  ground = Set{T}()
  for p in parts_list
    for a in p
      push!(ground, a)
    end
  end

  if check
    @assert length(ground) == parts_sum "The sets in the input must be pairwise disjoint."
  end

  P = Partition(ground)

  for p in parts_list
    plist = collect(p)
    np = length(plist)
    for k=1:np-1
      merge_parts!(P,plist[k], plist[k+1])
    end
  end
  kill_sos!(P)
  return P
end

# Equality check
==(P::Partition, Q::Partition) = parts(P)==parts(Q)

# Hashing
import Base.hash
function hash(P::Partition{T},h::UInt64=UInt64(0)) where T
  if length(P.SOS)==0 && num_parts(P)>0
    build_sos!(P)
  end
  hash(P.SOS,h)
end

# Display
function show(io::IO, P::Partition)
  print(io, parts(P))
end

"""
`in(a,P)` checks if `a` is in the ground set of `P` and
`in(A,P)` checks if `A` is a part of `P`.
"""
function in(a::T,P::Partition{T}) where T
  return in(a,P.elements)
end

function in(A::Set{T}, P::Partition{T}) where T
  return in(A, parts(P))
end

"""
`merge_parts!(P,a,b)` updates `P` by merging the parts
that contain elements `a` and `b`.
"""
function merge_parts!(P::Partition{T},a::T,b::T) where T
  @assert in(a,P)&&in(b,P) "One or both of these elements is not in the partition."
  union!(P.parts,a,b)
  kill_sos!(P)
  nothing
end

"""
`merge_parts!(P,elts)` merges all the elements in `elts` into a
part. Thus `merge_parts!(P,[a,b])` is equivalent to
`merge_parts!(P,a,b)`.
"""
function merge_parts!(P::Partition{T}, elts::Vector{T}) where T
  ne = length(elts)
  for k=1:ne-1
    merge_parts!(P,elts[k],elts[k+1])
  end
end

"""
`num_parts(P)` gives the number of parts in the partition `P`.
"""
num_parts(P::Partition) = num_groups(P.parts)

"""
`num_elements(P)` gives the number of elements in the ground
set of the partition `P`. This equals the sum of the sizes of the
parts.
"""
num_elements(P::Partition) = length(P.elements)

"""
`elements(P)` returns (a copy of) the ground set of the
partition `P`.
"""
ground_set(P::Partition) = deepcopy(P.elements)

import Base.collect

"""
`collect(P::Partition)` returns a list of the parts of the
partition `P`.
"""
collect(P::Partition) = collect(parts(P))

"""
`parts(P)` returns a set containing the parts of the partition `P`.
That is, we return a set of sets.
"""
function parts(P::Partition{T}) where T
  if length(P.SOS)==0 && num_parts(P)>0
    build_sos!(P)
  end
  return deepcopy(P.SOS)
end

# using internally by parts(P); not exposed
function build_sos!(P::Partition{T}) where T
  n = num_parts(P)
  GS = P.elements

  # Find the indices of the roots of the parts
  root_set = Set{T}()
  for item in GS
    r = find_root(P.parts,item)
    push!(root_set,r)
  end
  roots = collect(root_set)  # make it an array

  # Create a mapping from root numbers to [1:n]
  rootmap = Dict{T,Int}()
  for k=1:n
    rootmap[roots[k]] = k
  end

  # Create an array of sets to hold the parts
  plist = Array{Set{T}}(undef,n)
  for k=1:n
    plist[k] = Set{T}()  # make sure they're empty sets
  end

  # For each element in the ground set, look up its root
  # in order to place it in its appropriate set.
  for item in GS
    idx = rootmap[find_root(P.parts,item)]
    push!(plist[idx], item)
  end

  # now take those parts and pack them into a set
  S = Set{Set{T}}()
  for item in plist
    push!(S,item)
  end
  P.SOS = S
  nothing
end

"""
`in_same_part(P,a,b)` returns `true` if `a` and `b` are in the same part
of the partition `P`. An error is thrown if either is not in the ground
set of `P`.
"""
function in_same_part(P::Partition{T},a::T,b::T) where T
  @assert in(a,P)&&in(b,P) "One or both of these elements is not in the partition."
  return find_root(P.parts,a) == find_root(P.parts,b)
end

"""
`find_part(P,a)` returns the part of `P` that contains `a` (or throws an
error if `a` is not in the ground set).
"""
function find_part(P::Partition{T},a::T) where T
  @assert in(a,P) "$a is not in the ground set of this partition."
  r = find_root(P.parts,a)
  A = Set{T}()
  for x in P.elements
    if find_root(P.parts,x) == r
      push!(A,x)
    end
  end
  return A
end

"""
`join(P,Q)` computes the join of the two partitions. This may also
be invoked as `P+Q`.
"""
function join(P::Partition{T}, Q::Partition{T}) where T
  @assert P.elements==Q.elements "The ground sets of the two partitions must be the same."

  R = Partition(P.elements)
  PP = union(parts(P), parts(Q))

  for p in PP
    items = collect(p)
    n = length(items)
    for i=1:n-1
      merge_parts!(R,items[i],items[i+1])
    end
  end
  return R
end

"""
For partitions `P` and `Q`, `P+Q` is their join.
"""
(+)(P::Partition, Q::Partition) = join(P,Q)

"""
For partitions `P` and `Q`, `P∨Q` is their join.
"""
(∨)(P::Partition, Q::Partition) = join(P,Q)

"""
For a partition `P` and element `x`, `P+x` builds a new partition
that adds `x` as a singleton element.
"""
function (+)(P::Partition{T}, x::T) where T
  @assert !(in(x,P.elements)) "This element is already in the ground set of this partition."
  sets = parts(P)
  singlet = Set{T}(x)
  push!(sets,singlet)
  return PartitionBuilder(sets)
end

"""
For a partition `P` and a set `A`, `P+A` creates a new partition
with `A` as a new part. Note that `A` and the ground set of `P`
must be disjoint.
"""
function (+)(P::Partition{T}, A::Set{T}) where T
  @assert length(intersect(P.elements,A))==0 "The added set must be dijsoint from the ground set of the partition."
  sets = parts(P)
  push!(sets,A)
  return PartitionBuilder(sets)
end

"""
`meet(P,Q)` computes the meet of the partitions. This may
also be invoked as `P*Q`.
"""
function meet(P::Partition{T}, Q::Partition{T}) where T
  @assert P.elements==Q.elements "The ground sets of the two partitions must be the same."

  R = Partition(P.elements)
  elist = collect(P.elements)
  n = length(elist)

  for i=1:n-1
    a = elist[i]
    for j=i+1:n
      b = elist[j]
      if in_same_part(P,a,b) && in_same_part(Q,a,b)
        merge_parts!(R,a,b)
      end
    end
  end
  return R
end

"""
For partitions `P` and `Q`, `P*Q` is their meet.
"""
(*)(P::Partition, Q::Partition) = meet(P,Q)

"""
For partitions `P` and `Q`, `P∧Q` is their meet.
"""
(∧)(P::Partition, Q::Partition) = meet(P,Q)

"""
`refines(P,Q)` determines if `P` is a refinement of `Q`. That is,
is every part of `P` a subset of a part of `Q`? The two partitions
must have the same ground set of else an error is thrown.

`refines(P,Q)` can be invoked as `P<=Q`. The variants
`P<Q`, `P>=Q`, and `P>Q` operate as expected. Note that partitions
are only partially ordered by refinement and one can easily construct
partitions `P` and `Q` for which both `P<=Q` and `Q<=P` are false.
"""
function refines(P::Partition{T}, Q::Partition{T}) where T
  @assert P.elements==Q.elements "The two partitions must have the same ground set."

  Pparts = collect(parts(P))
  for S in Pparts
    elts = collect(S)
    n = length(elts)
    for i=1:n-1
      if !in_same_part(Q,elts[i],elts[i+1])
        return false
      end
    end
  end
  return true
end

# Comparison operators based on "is finer". Only a partial order!!
(<=)(P::Partition,Q::Partition) = refines(P,Q)
(<)(P::Partition,Q::Partition) = refines(P,Q) && !(P==Q)
(>=)(P::Partition,Q::Partition) = (Q<=P)
(>)(P::Partition,Q::Partition) = (Q<P)

# include("random_partitions.jl")  # not ready for prime time
include("all_partitions.jl")
include("integer_partition.jl")

end  # end of module


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

JLLWrappers.@generate_main_file_header("CompilerSupportLibraries")
JLLWrappers.@generate_main_file("CompilerSupportLibraries", UUID("e66e0078-7015-5450-92f7-15fbd957f2ae"))
end  # module CompilerSupportLibraries_jll


text:
cd(@__DIR__)
using Pkg
Pkg.activate(".")
Pkg.instantiate()
using JSON, Statistics, GitHub, Base64, SHA, Downloads, Dates, CairoMakie

include("benchmark-library.jl")

ctx = github_context()

project = "current-pr-project/"
# It seems, that between julia versions, the manifest must be deleted to not get problems
isdir(project) && rm(project; force=true, recursive=true)
mkdir(project)

Pkg.activate(project)
Pkg.develop([(;path="../../MakieCore"), (;path="../../"), (;path="../../CairoMakie"), (;path="../../../GeometryBasics")])
this_pr = BenchInfo(
    project=project,
    branch="current-pr",
    commit=string(current_commit())
)

Pkg.activate(".")

get_benchmark_data(ctx, this_pr; n=10, force=true)

plot_url = run_benchmarks(ctx, [GitHub.branch(ctx.repo, "master"), this_pr])

i1 = BenchInfo(commit="444eae5ce174d23d53c181144b357382bd57afa8", branch="julia1.9+this-pr")
i2 = BenchInfo(commit="7ccf35e789fccdda2429b39b6b10f4e91adcc5fd", branch="julia1.9+no-static-arrays")
i3 = BenchInfo(julia="julia-1-7-1", commit="f04ea93cde4384bad3f008bddf291a15ac35774a", branch="julia1.7+more-precompiles")
i4 = BenchInfo(julia="julia-1-7-1", commit="4a127bfd8dee589bbc3c2da1c91573a73752f36c", branch="julia1.7+v0.16.5")
i5 = BenchInfo(julia="julia-1-9-0-DEV-5", commit="7841b37e84921c4be8ec028b6aeb6cd12c036650", branch="master-yesterday")

fig = run_benchmarks(ctx, [i4, i3, i5, GitHub.branch(ctx.repo, "master"), i2, i1])

resize!(fig.scene, (1500, 600))
fig


text:
#NOTE this might be an overkill and using vectors with indices might be
# sufficient

# assumes throughout that T is an Ntuple of Booleans
import Base: length

"""
    reformat_updt_coord(updt_coord::Nothing, θ)

Chosen not to update parameters, returned object is not important
"""
reformat_updt_coord(updt_coord::Nothing, θ) = (Val((true,)),)

IntContainer = Union{Number,NTuple{N,<:Integer},Vector{<:Integer}} where N
IntLongContainer = Union{NTuple{N,<:Integer},Vector{<:Integer}} where N
DoubleContainer = Union{Vector{<:IntLongContainer},NTuple{N,Tuple}} where N
"""
    reformat_updt_coord(updt_coord::S, θ) where S<:IntContainer

Single joint update of multiple parameters at once
"""
#function reformat_updt_coord(updt_coord::S, θ) where S<:IntContainer
#    (int_to_val(updt_coord, θ),)
#end
#WARNING the above might have been commented out whereas it is lines 128--131 in
#the file mcmc_components that should have been commented out...

"""
    reformat_updt_coord(updt_coord::S, θ) where S<:DoubleContainer

Multiple updates, reformat from indices to update to a tuple of Val{...}()
"""
function reformat_updt_coord(updt_coord::S, θ) where S<:DoubleContainer
    Tuple([int_to_val(uc, θ) for uc in updt_coord])
end

"""
    int_to_val(updt_coord::S, θ)

Transform from a container of indices to update to a
Val{tuple of true/false one-hot-encoding}()
"""
function int_to_val(updt_coord::S, θ) where S<:IntContainer
    @assert all([1 <= uc <= length(θ) for uc in updt_coord])
    Val{Tuple([i in updt_coord for i in 1:length(θ)])}()
end

"""
    reformat_updt_coord(updt_coord::Nothing, θ)

If the user does not use indices of coordinates to be updated it is assumed that
appropriate Val{(...)}() object is passed and nothing is done, use at your own risk
"""
reformat_updt_coord(updt_coord, θ) = updt_coord

@generated function thetaex(::Val{T}, θ) where T
    z = Expr(:tuple, 1.0, (:(θ[$i]) for i in 1:length(T) if  !T[i])...)
    return z
end

@generated function thetainc(::Val{T}, θ) where T
    z = Expr(:tuple, (:(θ[$i]) for i in 1:length(T) if T[i])...)
    return z
end

@generated function indices(::Val{T}) where T
    z = Expr(:tuple, (i for i in 1:length(T) if T[i])...)
    return z
end

function matinc(mat, coord_idx)
    idx = [indices(coord_idx)...]
    Matrix(view(mat, idx, idx))
end

"""
    moveToProperPlace(ϑ, θ, ::Val{T}) where {T}

Update parameter vector `θ` at indices specified by `Val{T}` with the
values collected in `ϑ`.

# Examples
```julia-repl
julia> moveToProperPlace([10,20,30], [1,2,3,4,5],
                         Val((true, false, true, false, true)))
5-element Array{Float64,1}:
 10.0
  2.0
 20.0
  4.0
 30.0
```
"""
function move_to_proper_place(ϑ, θ, coord_idx::Val{T}) where T
    θᵒ = copy(θ)
    θᵒ[[indices(coord_idx)...]] = wrap(ϑ)
    θᵒ
end

function move_to_proper_place!(ϑ, θ, coord_idx::Val{T}) where T
    θ[[indices(coord_idx)...]] = wrap(ϑ)
    θ
end

wrap(x::Any) = x
wrap(x::Number) = [x]
truelength(::Val{T}) where T = sum(T)
length(::Val{T}) where T = length(T)


text:
using LinearAlgebra
function vech(x)
    a = triu(x)
    a = a[findall(x->x!=0,a)]
end

max_stars_repo_name: manuelmontesinos/Econometrics

max_issues_repo_name: manuelmontesinos/Econometrics

max_forks_repo_name: manuelmontesinos/Econometrics


text:
using Distributions
using Base.Test

# Test for gradloglik on univariate distributions

@test_approx_eq_eps gradloglik(Beta(1.5, 3.0), 0.7) -5.9523809523809526 1.0e-8
@test_approx_eq_eps gradloglik(Chi(5.0), 5.5) -4.7727272727272725 1.0e-8
@test_approx_eq_eps gradloglik(Chisq(7.0), 12.0) -0.29166666666666663 1.0e-8
@test_approx_eq_eps gradloglik(Exponential(2.0), 7.0) -0.5 1.0e-8
@test_approx_eq_eps gradloglik(Gamma(9.0, 0.5), 11.0) -1.2727272727272727 1.0e-8
@test_approx_eq_eps gradloglik(Gumbel(3.5, 1.0), 4.0) -1.6065306597126334 1.0e-8
@test_approx_eq_eps gradloglik(Laplace(7.0), 34.0) -1.0 1.0e-8
@test_approx_eq_eps gradloglik(Logistic(-6.0), 1.0) -0.9981778976111987 1.0e-8
@test_approx_eq_eps gradloglik(LogNormal(5.5), 2.0) 1.9034264097200273 1.0e-8
@test_approx_eq_eps gradloglik(Normal(-4.5, 2.0), 1.6) -1.525 1.0e-8
@test_approx_eq_eps gradloglik(TDist(8.0), 9.1) -0.9018830525272548 1.0e-8
@test_approx_eq_eps gradloglik(Weibull(2.0), 3.5) -6.714285714285714 1.0e-8

# Test for gradloglik on multivariate distributions

@test_approx_eq_eps gradloglik(MvNormal([1., 2.], [1. 0.1; 0.1 1.]), [0.7, 0.9]) [0.191919191919192, 1.080808080808081] 1.0e-8
@test_approx_eq_eps gradloglik(MvTDist(5., [1., 2.], [1. 0.1; 0.1 1.]), [0.7, 0.9]) [0.2150711513583442, 1.2111901681759383] 1.0e-8


text:
applymany(fs, x) = map(xi->xi(x), fs)
@spec all([_res[i] = f[i](x) for f in fs])

"""
Transpose for nested list / tuples.
Useful for output of rand(::NTuple{RandVar})

```jldoctest
x = normal(0.0, 1.0)
y = normal(0.0, 1.0)
samples = rand((x, y), x == y)
x_, y_ = ntranspose(samples)
```
"""
ntranspose(xs) = [[x[i] for x in xs] for i = 1:length(xs[1])]
@spec :incomplete same([length(x) for x in xs])

"Counter"
mutable struct Counter{T}
  count::T
end
Counter(::Type{T}) where T = Counter(zero(T))
Counter() = Counter(UInt)

"Increment counter"
increment!(c::Counter) = x::Int = c.count += 1
@spec UTupleec c.count == _pre(c.count) += 1

reset!(c::Counter{T}) where T = c.count = zero(T)

UTuple{T} = Tuple{Vararg{T, N}} where N

# FIXME, dynamic dispatch
"merge using combine if key is shared.
Assumes result of combine same type as in `b`

```jldoctest
x = (a = 1, b = 2)
y = (b = 2, c = 3)
merge(*, x, y)
```
"
function Base.merge(combine::Function, a::NamedTuple{an}, b::NamedTuple{bn}) where {an, bn}
  names = merge_names(an, bn)
  types = merge_types(names, typeof(a), typeof(b))
  function resolve(n)
    if sym_in(n, bn)
      if sym_in(n, an)
        combine(Val{n}, getfield(a, n), getfield(b, n))
      else
        getfield(b, n)
      end
    else
      getfield(a, n)
    end
  end
  NamedTuple{names, types}(map(resolve, names))::NamedTuple{names, types} # Inference fails without this
end

"Symbol concatenation"
x::Symbol *ₛ y::Symbol = Symbol(string(x), string(y)) # Is there better way?


text:
using FastGaussQuadrature, Memoization, SpecialFunctions

const RTOL = eps()

abstract type AbstractShape end
struct GaussianShape <: AbstractShape
  σ::Float64
  function GaussianShape(σ)
    σ > 0 || throw(error(ArgumentError("GuassianShape σ must be > 0 but is $σ")))
    return new(σ)
  end
end
width(s::GaussianShape) = 12.5 .* s.σ
sigma(s::GaussianShape) = s.σ
(s::GaussianShape)(x, centre) = exp(-(x-centre)^2 / s.σ^2) / √π / s.σ
knots(s::GaussianShape) = 0:1

struct BSpline{N} <: AbstractShape
  Δ::Float64
  function BSpline{N}(Δ) where {N}
    Δ > 0 || throw(error(ArgumentError("BSpline{$N} Δ must be > 0 but is $Δ")))
    return new{N}(Δ)
  end
end
width(s::BSpline{N}) where N = (N + 1) * s.Δ

function (s::BSpline{0})(x, centre)
  return (-s.Δ/2 <= x - centre < s.Δ/2) / s.Δ
end
const TopHatShape = BSpline{0}

function (s::BSpline{1})(x, centre)
  z = 2 * (x - centre + 1 * s.Δ) / width(s) # z is between 0 and 2
  value = if 0 <= z < 1
    z
  elseif 1 <= z < 2
    (2 - z)
  else
    zero(z)
  end
  return value / s.Δ
end
const TentShape = BSpline{1}

function (s::BSpline{2})(x, centre)
  z = 3 * (x - centre + 1.5 * s.Δ) / width(s) # z is between 0 and 3
  value = if 0 <= z < 1
    z^2 / 2
  elseif 1 <= z < 2
    3/4 - (1.5 - z)^2
  elseif 2 <= z < 3
    (3 - z)^2 / 2
  else
    zero(z)
  end
  return value / s.Δ
end
const QuadraticBSplineShape = BSpline{2}

knots(b::BSpline{N}) where N = 0:1/(N + 1):1

struct DeltaFunctionShape <: AbstractShape end
DeltaFunctionShape(_) = DeltaFunctionShape() # so has same ctor as others
width(s::DeltaFunctionShape) = 0

mutable struct BasisFunction{S<:AbstractShape, T}
  shape::S
  centre::Float64
  weight::T
end
BasisFunction(s::AbstractShape, centre) = BasisFunction(s, centre, 0.0)

lower(b::BasisFunction) = b.centre - width(b.shape) / 2
upper(b::BasisFunction) = b.centre + width(b.shape) / 2
lower(a::BasisFunction, b::BasisFunction) = max(lower(a), lower(b))
upper(a::BasisFunction, b::BasisFunction) = min(upper(a), upper(b))
width(b::BasisFunction) = width(b.shape)
weight(b::BasisFunction) = b.weight
centre(b::BasisFunction) = b.centre
shape(b::BasisFunction) = b.shape
knots(b::BasisFunction) = knots(shape(b)) .* width(b) .+ lower(b)
function getkey(u::BasisFunction, v::BasisFunction)::UInt64
  hashtuple = (width(u), width(v), trunc((centre(u) - centre(v)), digits=14))
  return foldr(hash, hashtuple; init=hash(typeof(u), hash(typeof(v))))
end
function translate(b::BasisFunction, x::Number)
  @assert isfinite(x)
  translated = deepcopy(b)
  translated.centre += x
  return translated
end
function translate!(b::BasisFunction, x::Number, bc::AbstractBC)
  b.centre = bc(b.centre + x)
  return b
end
Base.:+(b::BasisFunction, x) = (@assert isfinite(x); b.weight += x; b)

zero!(b::BasisFunction) = (b.weight *= false)
function Base.in(x::Number, b::BasisFunction)
  return (b.centre - width(b)/2 <= x < b.centre + width(b)/2)
end

function knots(a::BasisFunction, b::BasisFunction)
  overlap(a, b) || return output
  ka = knots(a)
  kb = knots(b)
  lo = max(ka[1], kb[1])
  hi = min(ka[end], kb[end])
  output = filter(x->lo <= x <= hi, unique(vcat(ka, kb)))
  sort!(output)
  return output
end

function translate(a::BasisFunction{S1}, b::BasisFunction{S2},
    p::PeriodicGridBC) where {S1<:AbstractShape, S2<:AbstractShape}
  lower(p) <= lower(a) && upper(a) < upper(p) && lower(p) <= lower(b) && upper(b) < upper(p) && return (a, b)
  in(a, b) && return (a, b)
  t = translate(a, length(p)); in(t, b) && return (t, b)
  t = translate(a,-length(p)); in(t, b) && return (t, b)
  return (a, b)
end

(b::BasisFunction)(x::Number) = b.shape(x, b.centre)
function (b::BasisFunction)(x::Number, p::PeriodicGridBC)
  in(x, b) && return b(x)
  t = translate(b, length(p)); in(x, t) && return t(x)
  t = translate(b,-length(p)); in(x, t) && return t(x)
  return b(x)
end

function overlap(a::BasisFunction, b::BasisFunction)
  return lower(a) == lower(b) || (lower(a) < upper(b) && upper(a) > lower(b))
end
Base.in(a::BasisFunction, b::BasisFunction) = overlap(a, b)
function overlap(u::BasisFunction, v::BasisFunction, p::PeriodicGridBC)
  return in(translate(u, v, p)...)
end

function integral(b::BasisFunction{DeltaFunctionShape}, f::F, _::AbstractBC) where {F<:Function}
  return f(centre(b))
end

function integral(b::BasisFunction, f::F, _::AbstractBC) where {F}
  ks = knots(b)
  return mapreduce(i->
    QuadGK.quadgk(x->b(x) * f(x), ks[i], ks[i+1]; order=7, atol=eps(), rtol=RTOL)[1],
    +, 1:length(ks)-1)
end

function integral(u::BasisFunction{S1},
                  v::BasisFunction{S2},
                  p::PeriodicGridBC) where {S1<:AbstractShape, S2<:AbstractShape}
  u, v = translate(u, v, p)
  return in(u, v) ? integral(u, v) : 0.0
end

function integral(u::BasisFunction{BSpline{N1}}, v::BasisFunction{BSpline{N2}}
    ) where {N1, N2}
  output = 0.0
  kn = knots(u, v)
  x, w = FastGaussQuadrature.gausslegendre((N1 + 1) * (N2 + 1))
  for i ∈ 1:length(kn)-1
    a, b = kn[i], kn[i+1]
    for j ∈ eachindex(x)
      xj = (x[j] + 1) / 2 * (b - a) + a
      wj = w[j] / 2 * (b - a)
      output += wj * u(xj) * v(xj)
    end
  end
  return output
end
function integral(u::BasisFunction{GaussianShape, T1},
                  v::BasisFunction{BSpline{N}, T2}) where {T1, N, T2}
  ks = knots(v)
  return mapreduce(i->
    QuadGK.quadgk(x->u(x) * v(x), ks[i], ks[i+1]; order=7, atol=eps(), rtol=RTOL)[1],
    +, 1:length(ks)-1)
end
function integral(u::BasisFunction{BSpline{N}, T1},
                  v::BasisFunction{GaussianShape, T2}) where {N, T1, T2}
  return integral(v, u)
end

function integral(u::BasisFunction{TopHatShape}, v::BasisFunction{GaussianShape})
  return integral(v, u)
end
function integral(u::BasisFunction{GaussianShape}, v::BasisFunction{TopHatShape})
  return (erf(-(u.centre - upper(v)) / u.shape.σ) -
          erf(-(u.centre - lower(v)) / u.shape.σ)) / 2 / width(v)
end
function integral(u::BasisFunction{TopHatShape},
                  v::BasisFunction{TopHatShape})
  w = upper(u, v) - lower(u, v)
  h = 1 / (u.shape.Δ * v.shape.Δ)
  return w * h * overlap(u, v)
end
integral(a::BasisFunction{DeltaFunctionShape}, b::BasisFunction) = b(a.centre)
integral(a::BasisFunction, b::BasisFunction{DeltaFunctionShape}) = a(b.centre)

function integral(a::BasisFunction{GaussianShape},
                  b::BasisFunction{GaussianShape})
  σ₁, σ₂ = a.shape.σ, b.shape.σ
  μ₁, μ₂ = a.centre, b.centre
  σ² = σ₁^2 + σ₂^2
  return exp(-(μ₁ - μ₂)^2 / σ²) / √(π * σ²)
end

function integral(a::BasisFunction{GaussianShape}, b::BasisFunction{TentShape})
  σ = a.shape.σ
  μ = centre(a)
  Δ = b.shape.Δ
  c = centre(b)
  foo(l, m) = (σ * exp(-(m - l)^2/σ^2)/√π + m * erf((m - l)/σ)) / 2Δ^2
  return foo(0, μ - c + Δ) - foo(Δ, μ - c + Δ) +
    foo(0, μ - c - Δ) - foo(-Δ, μ - c - Δ)
end


text:
##
#reload("EMST.jl")
using EMST
@static if VERSION < v"0.7.0-DEV.2005"
    using Base.Test
else
    using Test
end

function test_emst_with_uniform(d,n;check_verification=false)
    x = rand(d,n)
    # x = rand(2,2000)

    x_emst    = EMST.compute_emst(x)
    good_emst = EMST.verify_emst(x,x_emst,size(x,2))
    @test good_emst

    if(check_verification)
        # and show that verify works..
        x_emst_kaputt = copy(x_emst)
        # change n edges..
        for zi in 1:min(n-1,5) ; x_emst_kaputt[ rand( 1:( size(x,2))-1) , rand(1:2) ] = rand(1:size(x,2)); end
        emst_verification_working =  ~verify_emst(x,x_emst_kaputt,size(x,2))
        @test emst_verification_working
    end
end

@testset "EMST Tests A" begin
    for zi=1:100
        test_emst_with_uniform(2,500)
    end
    for zi=1:100
        test_emst_with_uniform(17,500)
    end
    # and check that our verification actually works.. :)
    for zi=1:20
        test_emst_with_uniform(4,1000;check_verification=true)
    end
end

@testset "EMST Tests B" begin
    test_emst_with_uniform(2,1000)
    test_emst_with_uniform(123,1000)
    test_emst_with_uniform(10,4000)
    test_emst_with_uniform(10,8000)
    test_emst_with_uniform(2,16000)
    test_emst_with_uniform(44,16000)
end

#using Plots
#using Gadfly
#function plot_emst_2d(x,edges)
#    Plots.gr()
#    Plots.plot()
#    Plots.scatter!(x[1,:],x[2,:],ms=2.0,marker=(stroke(0,:gray)))
#    for zr in 1:size(edges,1)
#        Plots.plot!( x[1,edges[zr,:]] , x[2,edges[zr,:]] )
#    end
#    Plots.scatter!() # for some reason with this the plot shows im juno.. :)
#end
#x_test  = rand(2,2000)
#tic(); e_test  = EMST.compute_emst(x_test;nmin=64); toc();
#EMST.verify_emst(x_test,e_test,200)
#plot_emst_2d(x_test,e_test)

##


text:
import BSON
# import CSV
import CSVFiles
import DataFrames
import FileIO
import JLD2
import ProgressMeter

"""
"""
function save_model(
        filename::AbstractString,
        fittable_object_to_save::AbstractFittable,
        )
    if filename_extension(filename) == ".jld2"
        save_result = save_model_jld2(filename,fittable_object_to_save)
    elseif filename_extension(filename) == ".bson"
        save_result = save_model_bson(filename,fittable_object_to_save)
    else
        error("extension must be one of: .jld2, .bson")
    end
    return filename
end

function save_model_jld2(
        filename::AbstractString,
        fittable_object_to_save::AbstractFittable,
        )
    if filename_extension(filename) != ".jld2"
        error(
            string(
                "Filename \"",
                filename,
                "\" does not end in \".jld2\"")
            )
    end
    dict_of_objects_to_save = Dict(
        "jld2_saved_model" => fittable_object_to_save,
        )
    @debug("Attempting to save model...")
    # make sure the parent directory exists
    parent_directory = Base.Filesystem.dirname(filename)
    try
        Base.Filesystem.mkpath(parent_directory)
    catch
    end
    # save the .jld2 file
    FileIO.save(filename, dict_of_objects_to_save)
    @debug(string("Saved model to file \"", filename, "\""))
    return filename
end

function save_model_bson(
        filename::AbstractString,
        fittable_object_to_save::AbstractFittable,
        )
    if filename_extension(filename) != ".bson"
        error(
            string(
                "Filename \"",
                filename,
                "\" does not end in \".bson\"")
            )
    end
    dict_of_objects_to_save = Dict(
        :bson_saved_model => fittable_object_to_save,
        )
    @debug("Attempting to save model...")
    # make sure the parent directory exists
    parent_directory = Base.Filesystem.dirname(filename)
    try
        Base.Filesystem.mkpath(parent_directory)
    catch
    end
    # save the .bson file
    BSON.bson(filename, dict_of_objects_to_save)
    @debug(string("Saved model to file \"", filename, "\""))
    return filename
end

"""
"""
function load_model(filename::AbstractString)
    if filename_extension(filename) == ".jld2"
        load_result = load_model_jld2(filename)
        return load_result
    elseif filename_extension(filename) == ".bson"
        load_result = load_model_bson(filename)
        return load_result
    else
        error("extension must be one of: .jld2, .bson")
    end
end

function load_model_jld2(filename::AbstractString)
    if filename_extension(filename) != ".jld2"
        error(
            string(
                "Filename \"",
                filename,
                "\" does not end in \".jld2\"")
            )
    end
    @debug("Attempting to load model...")
    dict_of_loaded_objects = FileIO.load(filename)
    loaded_fittable_object = dict_of_loaded_objects["jld2_saved_model"]
    @debug(string("Loaded model from file \"", filename, "\""))
    return loaded_fittable_object
end

function load_model_bson(filename::AbstractString)
    if filename_extension(filename) != ".bson"
        error(
            string(
                "Filename \"",
                filename,
                "\" does not end in \".bson\"")
            )
    end
    @debug("Attempting to load model...")
    dict_of_loaded_objects = BSON.load(filename)
    loaded_fittable_object = dict_of_loaded_objects[:bson_saved_model]
    @debug(string("Loaded model from file \"", filename, "\""))
    return loaded_fittable_object
end


text:
"""
Module to test the PlayingCards module
"""
module PlayingCardsTest

using  Test

using  PlayingCards

@testset "Card Constructors" begin
    for value in 0:51
        c = Card(value)
        @test c isa Card
        @test Card(suitint(c), rankint(c)) == c
        @test Card(suit(c), rank(c)) == c
        @test Card(string(c)) == c
    end
    @test_throws AssertionError Card(52)

end;

@testset "Card Suits and Ranks" begin
    c = Card(0)
    @test suitint(c) == 1
    @test suit(c) == "♣"
    @test rankint(c) == 1
    @test rank(c) == "A"

    c = Card(51)
    @test suitint(c) == 4
    @test suit(c) == "♠"
    @test rankint(c) == 13
    @test rank(c) == "K"
end;

@testset "Deck creation" begin
    deck = createdeck()
    @test length(deck) == 52
    for (i, card) in enumerate(deck)
        @test i == card.value + 1
    end;
end;

@testset "Dealing cards" begin
    deck = createdeck()
    hands = deal!(deck, 13, 4)
    @test length(hands) == 4
    for h in hands
        @test length(h) == 13
        for c in h
            @test c ∉ deck
        end
    end

    deck = createdeck()
    hands = deal!(deck, 6, 2)
    @test length(hands) == 2
    for h in hands
        @test length(h) == 6
        for c in h
            @test c ∉ deck
        end
    end

    deck = createdeck()
    @test_throws AssertionError deal!(deck, 15, 5)
end;

@testset "Parsing Cards" begin
    for suit in PlayingCards.suits, rank in PlayingCards.ranks
        cardstr = join([rank, suit], "")
        @test string(parse(Card, cardstr)) == cardstr
    end

    hand = parse(Array{Card,1}, "8♣ 7♥ 7♠ 6♦")
    @test length(hand) == 4
    for card in hand
        @test card isa Card
    end
    @test string(hand) == "8♣ 7♥ 7♠ 6♦"
    handtuple = tuple(hand...)
    @test string(handtuple) == "8♣ 7♥ 7♠ 6♦"
end;

end


text:
import Base: isapprox, ctranspose
import StatsBase: sample, pweights

isapprox(x::Ket, y::Ket) = (x.basis == y.basis) && x.data ≈ y.data
isapprox(x::Bra, y::Bra) = (x.basis == y.basis) && x.data ≈ y.data
isapprox(x::SparseOperator, y::SparseOperator) = (x.basis_l == y.basis_l) && (x.basis_r == y.basis_r) && x.data ≈ y.data
isapprox(x::DenseOperator, y::DenseOperator) = (x.basis_l == y.basis_l) && (x.basis_r == y.basis_r) && x.data ≈ y.data
isapprox(x::SparseOperator, y::DenseOperator) = (x.basis_l == y.basis_l) && (x.basis_r == y.basis_r) && x.data ≈ y.data
isapprox(x::DenseOperator, y::SparseOperator) = (x.basis_l == y.basis_l) && (x.basis_r == y.basis_r) && x.data ≈ y.data

ctranspose(op::QuantumOptics.particle.FFTKets) = dagger(op::QuantumOptics.particle.FFTKets)
ctranspose(op::QuantumOptics.particle.FFTOperators) = dagger(op::QuantumOptics.particle.FFTOperators)
ctranspose(op::QuantumOptics.operators_lazytensor.LazyTensor) = dagger(op::QuantumOptics.operators_lazytensor.LazyTensor)
ctranspose(op::QuantumOptics.operators_lazyproduct.LazyProduct) = dagger(op::QuantumOptics.operators_lazyproduct.LazyProduct)
ctranspose(op::QuantumOptics.operators_lazysum.LazySum) = dagger(op::QuantumOptics.operators_lazysum.LazySum)
ctranspose(x::QuantumOptics.operators_sparse.SparseOperator) = dagger(x::QuantumOptics.operators_sparse.SparseOperator)
ctranspose(x::QuantumOptics.operators_dense.DenseOperator) = dagger(x::QuantumOptics.operators_dense.DenseOperator)
ctranspose(x::QuantumOptics.states.Ket) = dagger(x::QuantumOptics.states.Ket)
ctranspose(x::QuantumOptics.states.Bra) = dagger(x::QuantumOptics.states.Bra)
ctranspose(a::QuantumOptics.operators.Operator) = dagger(a::QuantumOptics.operators.Operator)

"""
    outer(x::Ket, y::Ket)

Construct outer product from two kets, ``|x⟩⟨y|``.
"""
outer(x::QuantumOptics.states.Ket, y::QuantumOptics.states.Ket) = x ⊗ y'

"""
    inner(x::Ket, y::Ket)

Calculate inner product from two kets, ``⟨x|y⟩``.
"""
inner(x::QuantumOptics.states.Ket, y::QuantumOptics.states.Ket) = x' * y

"""
    qubit(x)

Prepare n-qubit state from binary.

# Example
```julia
julia> qubit("010")
Ket(dim=8)
  basis: [Spin(1/2) ⊗ Spin(1/2) ⊗ Spin(1/2)]
 0.0+0.0im
 0.0+0.0im
 1.0+0.0im
 0.0+0.0im
 0.0+0.0im
 0.0+0.0im
 0.0+0.0im
 0.0+0.0im
```
"""
function qubit(x::String)
    nq = length(x)
    ar = parse.(Int, split(x, "")) + 1
    basis = SpinBasis(1//2)
    return tensor( basisstate.(basis, ar)... )
end

"""
    qudit(x, dim)

Prepare n-qudit state

# Example
```
julia> qudit("010", 3)
Ket(dim=27)
  basis: [NLevel(N=3) ⊗ NLevel(N=3) ⊗ NLevel(N=3)]
 0.0+0.0im
 0.0+0.0im
 0.0+0.0im
 1.0+0.0im
    ⋮
 0.0+0.0im
 0.0+0.0im
 0.0+0.0im

 julia> qudit("010", 3) |> tex
|ψ⟩ = |010⟩
```
"""
function qudit(x::String, dim::Int)
    nq = length(x)
    ar = parse.(Int, split(x, "")) + 1
    basis = NLevelBasis(dim)
    return tensor( nlevelstate.(basis, ar)... )
end

"""
    quxit(x::String, dims::Vector{Int})

Prepare a composite state which the dimension of i th system is `dims[i]`.

# Example
```
julia> quxit("021", [2,3,2])
Ket(dim=12)
  basis: [NLevel(N=2) ⊗ NLevel(N=3) ⊗ NLevel(N=2)]
 0.0+0.0im
 0.0+0.0im
 0.0+0.0im
 0.0+0.0im
    ⋮
 0.0+0.0im
 1.0+0.0im
 0.0+0.0im
```
"""
function quxit(x::String, dims::Vector{Int})
    @assert length(x) == length(dims)
    nq = length(x)
    ar = parse.(Int, split(x, "")) + 1
    return tensor(basisstate.(NLevelBasis.(dims), ar)...)
end

"""
    sparse_spinallup_dm(n)

n qubits all up state density matrix
```math
| 0 \\rangle \\langle 0 |^{\\otimes n}
```
"""
function sparse_spinallup_dm(n::Int)
       x = spzeros(Complex{Float64}, 2^n, 2^n)
       x[1,1] = 1. + 0im
       return SparseOperator(SpinBasis(1//2)^n, x)
end

"""
    sparse_spinalldown_dm(n)

n qubits all down state density matrix
```math
| 1 \\rangle \\langle 1 |^{\\otimes n}
```
"""
function sparse_spinalldown_dm(n::Int)
       x = spzeros(Complex{Float64}, 2^n, 2^n)
       x[2^n,2^n] = 1. + 0im
       return SparseOperator(SpinBasis(1//2)^n, x)
end

"""
       ghz(n)

n qubits GHZ(Greenberger–Horne–Zeilinger) state
```math
| \\mathrm{GHZ} \\rangle = \\frac{| 0 \\rangle^{\\otimes n} + | 1 \\rangle^{\\otimes n}}{\\sqrt{2}}
```

# Example
```julia
julia> n = 3;

julia> ghz(n)
Ket(dim=8)
  basis: [Spin(1/2) ⊗ Spin(1/2) ⊗ Spin(1/2)]
 0.707107+0.0im
      0.0+0.0im
      0.0+0.0im
      0.0+0.0im
      0.0+0.0im
      0.0+0.0im
      0.0+0.0im
 0.707107+0.0im
```
"""
function ghz(n::Int=3)
       b = SpinBasis(1//2)^n
       state = basisstate(b, 1) + basisstate(b, 2^n)
       normalize!(state)
       return state
end

"""
    measure(state, n)

Projective measurement with standard basis.
Measure n th qubit w.r.t standard basis.

Return outcome (0 or 1) and post-meamurement state
"""
function measure(state::Ket, n::Int=1)
    nqubit = state.basis.shape |> length
    mops = [identityoperator(SpinBasis(1//2)) for i in 1:nqubit]
    mops[n] = spinup(SpinBasis(1//2)) |> dm |> sparse
    mop = tensor(mops...)
    p = real( trace( mop * dm(state) ) )

    outcome = sample(0:1, pweights([p, 1-p]))
    if outcome == 0
        state = mop * state / sqrt(p)
    else
        mops[n] = spindown(SpinBasis(1//2)) |> dm |> sparse
        mop = tensor(mops...)
        state = mop * state / sqrt(1-p)
    end
    return outcome, state
end


text:
function jlty2llvmty(typ::Type, isboxed_ref=Ref{Bool}(false))
    ccall(:jl_type_to_llvm, LLVM.API.LLVMTypeRef, (Any, Ptr{Bool}), typ, isboxed_ref)
end


text:
module StructsOfArrays
export StructOfArrays

immutable StructOfArrays{T,N,U<:Tuple} <: AbstractArray{T,N}
    arrays::U
end

@generated function StructOfArrays{T}(::Type{T}, dims::Integer...)
    (!isleaftype(T) || T.mutable) && return :(throw(ArgumentError("can only create an StructOfArrays of leaf type immutables")))
    isempty(T.types) && return :(throw(ArgumentError("cannot create an StructOfArrays of an empty or bitstype")))
    N = length(dims)
    arrtuple = Tuple{[Array{T.types[i],N} for i = 1:length(T.types)]...}
    :(StructOfArrays{T,$N,$arrtuple}(($([:(Array($(T.types[i]), dims)) for i = 1:length(T.types)]...),)))
end
StructOfArrays(T::Type, dims::Tuple{Vararg{Integer}}) = StructOfArrays(T, dims...)

Base.linearindexing{T<:StructOfArrays}(::Type{T}) = Base.LinearFast()

@generated function Base.similar{T}(A::StructOfArrays, ::Type{T}, dims::Dims)
    if isbits(T) && length(T.types) > 1
        :(StructOfArrays(T, dims))
    else
        :(Array(T, dims))
    end
end

Base.convert{T,S,N}(::Type{StructOfArrays{T,N}}, A::AbstractArray{S,N}) =
    copy!(StructOfArrays(T, size(A)), A)
Base.convert{T,S,N}(::Type{StructOfArrays{T}}, A::AbstractArray{S,N}) =
    convert(StructOfArrays{T,N}, A)
Base.convert{T,N}(::Type{StructOfArrays}, A::AbstractArray{T,N}) =
    convert(StructOfArrays{T,N}, A)

Base.size(A::StructOfArrays) = size(A.arrays[1])
Base.size(A::StructOfArrays, d) = size(A.arrays[1], d)

@generated function Base.getindex{T}(A::StructOfArrays{T}, i::Integer...)
    Expr(:block, Expr(:meta, :inline),
         Expr(:new, T, [:(A.arrays[$j][i...]) for j = 1:length(T.types)]...))
end
@generated function Base.setindex!{T}(A::StructOfArrays{T}, x, i::Integer...)
    quote
        $(Expr(:meta, :inline))
        v = convert(T, x)
        $([:(A.arrays[$j][i...] = getfield(v, $j)) for j = 1:length(T.types)]...)
        x
    end
end
end # module


text:
@testset "$(rpad("Vandermonde Matrix Tests",80))" begin

    n = 11
    x = collect(range(0., stop=1., length=n))
    y = rand(n)

    V = vandermonde_matrix(x)
    Vinv = vandermonde_matrix_inverse(x)

    x = V\y
    a = *(Vinv, y)

    @test maximum(abs.((a .- x) ./ a)) < eps(Float32)

end


text:
module NurseData

using CSV
using Dates
using DataFrames

export n_nurses

basepath = normpath(@__DIR__, "../../")

function n_nurses(states::Array{String,1}; source::Symbol=:ahrf)
	if source == :ahrf
		return num_nurses_from_ahrf(states)
	elseif source == :employment
		return num_nurses_from_employment(states)
	elseif source == :beds
		return num_nurses_from_beds(states)
	else
		error("Invalid nurses data source")
	end
end

########################
##### Employment #######
########################

function load_nurse_employment_data()
	nurse_data_employment = DataFrame(CSV.File(joinpath(basepath, "data/nurses/deaggregated_by_hospital_beds.csv")), copycols=true)
	nurse_data_employment = combine(groupby(nurse_data_employment, :state), :weighted_emp_distribution => sum)
	sort!(nurse_data_employment, :state)
	return nurse_data_employment
end

function num_nurses_from_employment(states::Array{String,1})
	@assert states == sort(states)
	@assert !("HI" in states)

	nurse_data_employment = load_nurse_employment_data()
	nurse_data_local = filter(row -> row.state in states, nurse_data_employment)
	return Float32.(nurse_data_local.weighted_emp_distribution_sum .* (1_713_120 / 2_982_280))
end

########################
######## Beds ##########
########################

function load_nurse_beds_data()
	nurse_data_beds = DataFrame(CSV.File(joinpath(basepath, "data/hospitals/definitivehc.csv")), copycols=true)
	filter!(row -> !(row.HOSPITAL_TYPE in ["Psychiatric Hospital", "Rehabilitation Hospital"]), nurse_data_beds)
	filter!(row -> !(ismissing(row.NUM_STAFFED_BEDS) || ismissing(row.NUM_ICU_BEDS) || ismissing(row.HQ_STATE)), nurse_data_beds)
	filter!(row -> (row.NUM_STAFFED_BEDS > 0) && (row.NUM_ICU_BEDS > 0), nurse_data_beds)

	nurse_data_beds = combine(groupby(nurse_data_beds, :HQ_STATE), [:NUM_STAFFED_BEDS => sum, :NUM_ICU_BEDS => sum])
	nurse_data_beds.non_icu_beds_sum = nurse_data_beds.NUM_STAFFED_BEDS_sum - nurse_data_beds.NUM_ICU_BEDS_sum
	sort!(nurse_data_beds, :HQ_STATE)

	return nurse_data_beds
end

function num_nurses_from_beds(
		states::Array{String,1};
		nurse_hrs_per_week::Real=36,
		nurses_per_bed_regular::Real=0.2,
		nurses_per_bed_icu::Real=0.5
)
	@assert states == sort(states)
	nurse_data_beds = load_nurse_beds_data()
	nurse_data_local = filter(row -> row.HQ_STATE in states, nurse_data_beds)
	nurse_data_local.est_nurses_per_shift =
		  (nurses_per_bed_icu * nurse_data_local.NUM_ICU_BEDS_sum)
		+ (nurses_per_bed_regular * nurse_data_local.non_icu_beds_sum)
	n_shifts = (24*7 / nurse_hrs_per_week)
	return Float32.(nurse_data_local.est_nurses_per_shift * n_shifts)
end

########################
######## AHRF ##########
########################

function load_nurse_ahrf_data()
	nurse_data = DataFrame(CSV.File(joinpath(basepath, "data/nurses/nurses_per_county.csv")))
	nurse_data.registered_nurses = coalesce.(nurse_data.registered_nurses, 0.0)
	nurse_data.nurse_practitioners = coalesce.(nurse_data.nurse_practitioners, 0.0)
	nurse_data.tot_nurses = coalesce.(nurse_data.tot_nurses, 0.0)
	nurse_data_bystate = groupby(nurse_data, :state)
	nurse_data_ahrf = combine(nurse_data_bystate, [
		:registered_nurses => sum,
		:nurse_practitioners => sum,
		:tot_nurses => sum,
	])
	sort!(nurse_data_ahrf, :state)
	return nurse_data_ahrf
end

function num_nurses_from_ahrf(states::Array{String,1})
	@assert states == sort(states)
	nurse_data_ahrf = load_nurse_ahrf_data()
	nurses_local = filter(row -> row.state in states, nurse_data_ahrf)
	nurses = nurses_local.registered_nurses_sum
	replace!(nurses, missing => 0)
	return Float32.(nurses)
end

end;

max_stars_repo_name: flixpar/covid-resource-allocation

max_issues_repo_name: flixpar/covid-resource-allocation

max_forks_repo_name: flixpar/covid-resource-allocation


text:
using JLD2; using PyPlot

include("../../../src/CGS.jl")
using Main.CGS

@load "../../../data/T1_JWST_all_transits_mass_error.jld2" JWST_sims nplanet nsim mstar planet

clf()
for ip=1:nplanet
  mp_in = zeros(Float64,nsim); mp_out = zeros(Float64,nsim); sp_out = zeros(Float64,nsim);
  for i=1:nsim
    mp_in[i]  = JWST_sims[i].elements_in[ip+1,1]*mstar*MSUN/MEARTH
    mp_out[i] = JWST_sims[i].elements_best[ip+1,1]*mstar*MSUN/MEARTH
    sp_out[i] = sqrt(JWST_sims[i].cov[(ip-1)*5+1,(ip-1)*5+1])*mstar*MSUN/MEARTH
  end
  errorbar(mp_in,mp_out-mp_in,sp_out,fmt="o",label=string("T1",planet[ip]),linewidth=2)
end
plot([0,1.5],[0,0])
axis([0,1.5,-0.002,0.002])
legend(loc="upper left",fontsize=15,ncol=2)
xlabel(L"$M_{in}/M_\oplus$",fontsize=15)
ylabel(L"$(M_{out}-M_{in})/M_\oplus$",fontsize=15)
tight_layout()
savefig("../Recovered_masses_JWST_5yr_all_transits_NIRSPEC.pdf",bbox_inches="tight")


text:
using KernelDensityEstimate, RoME, Gadfly

XY = Float64[]
for i in 1:200  XY=[XY;pol2cart([100.0+randn();2pi*rand()-pi],ones(2))[1]]; end
xy1 = reshape(XY[1:200],2,100);
xy2 = reshape(XY[201:400],2,100);
xy1[1,:] -= 25.0
xy2[1,:] += 25.0

plot(layer(x=xy1[1,:],y=xy1[2,:],Geom.point),layer(x=xy2[1,:],y=xy2[2,:],Geom.point))

p1 = kde!(xy1);
p2 = kde!(xy2);
p12 = p1*p2;

plotKDE(marginal(p12,[1]))
plotKDE(marginal(p12,[2]))


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

## dummy stub for https://github.com/JuliaBinaryWrappers/libLLVM_jll.jl

baremodule libLLVM_jll
using Base, Libdl
Base.Experimental.@compiler_options compile=min optimize=0 infer=false

const PATH_list = String[]
const LIBPATH_list = String[]

export libLLVM

# These get calculated in __init__()
const PATH = Ref("")
const LIBPATH = Ref("")
artifact_dir = ""
libLLVM_handle = C_NULL
libLLVM_path = ""

if Sys.iswindows()
    const libLLVM = "libLLVM-14jl.dll"
elseif Sys.isapple()
    const libLLVM = "@rpath/libLLVM.dylib"
else
    const libLLVM = "libLLVM-14jl.so"
end

function __init__()
    global libLLVM_handle = dlopen(libLLVM)
    global libLLVM_path = dlpath(libLLVM_handle)
    global artifact_dir = dirname(Sys.BINDIR)
    LIBPATH[] = dirname(libLLVM_path)
    push!(LIBPATH_list, LIBPATH[])
end

# JLLWrappers API compatibility shims.  Note that not all of these will really make sense.
# For instance, `find_artifact_dir()` won't actually be the artifact directory, because
# there isn't one.  It instead returns the overall Julia prefix.
is_available() = true
find_artifact_dir() = artifact_dir
dev_jll() = error("stdlib JLLs cannot be dev'ed")
best_wrapper = nothing
get_libLLVM_path() = libLLVM_path

end  # module libLLVM_jll

max_stars_repo_name: TechPenguineer/julia

max_issues_repo_name: TechPenguineer/julia

max_forks_repo_name: TechPenguineer/julia