text:
export fit_sellmeier

const GUESS_MARGIN = 0.05  # how much λres guesses are separated from range of measured λ

function has_similar(λres::AbsVecReal; rtol=1e-3)  # assume λres is sorted
    for i = 1:length(λres)-1
        isapprox(λres[i], λres[i+1]; rtol) && return true  # if poles are different by less than 0.1%, return true
    end

    return false
end

"""
    fit_sellmeier(λ::AbstractVector{<:Real}, ε::AbstractVector{<:Real}; Nmax::Integer=10, rtol_λres::Real=1e-3)

Fit the dielectric constant data `ε` sampled at wavelengths `λ` to the Sellmeier equation
with up to `Nmax` terms.  If `λres` of different terms have relative difference less than
`rtol_λres`, they are considered the same term and merged.
"""
fit_sellmeier(λ::AbsVecReal, ε::AbsVecReal; kwargs...) = fit_sellmeier(float(λ), float(ε); kwargs...)
function fit_sellmeier(λ::AbsVecFloat,  # wavelengths where ε was measured
                       ε::AbsVecFloat;  # measured relative permittivities (= squared refractive indices)
                       Nmax::Integer=10,  # maximum number of terms in Sellmeier equation
                       rtol_λres::Real=1e-3)  # relative tolerance to distinguish λres
    me = fit_sellmeier(λ, ε, Val(1))  # me: model and error
    for N = 2:Nmax
        me′ = fit_sellmeier(λ, ε, Val(N))
        if me′.err > me.err  # if more terms than necessary are used, error starts to increase
            break
        elseif has_similar(me′.mdl.λres; rtol=rtol_λres)  # terms with similar poles are unnecessary
            break
        else
            me = me′
        end
    end

    return me
end

"""
    fit_sellmeier(λ::AbstractVector{<:Real}, ε::AbstractVector{<:Real}, N::Integer)

Fit the dielectric constant data `ε` sampled at wavelengths `λ` to the Sellmeier equation
with exactly `N` terms.  Among the `N` terms, some are below and the other are above the
range of `λ`.
"""
fit_sellmeier(λ::AbsVecReal, ε::AbsVecReal, N::Integer) = fit_sellmeier(λ, ε, Val(N))
fit_sellmeier(λ::AbsVecReal, ε::AbsVecReal, ::Val{N}) where {N} = fit_sellmeier(float(λ), float(ε), Val(N))
function fit_sellmeier(λ::AbsVecFloat,  # wavelengths where ε was measured
                       ε::AbsVecFloat,  # measured relative permittivities (= squared refractive indices)
                       ::Val{N}  # number of terms in Sellmeier model
                       ) where {N}
    local me
    err_min = Inf
    for Nₙ = 0:N
        Nₚ = N - Nₙ
        me′ = fit_sellmeier(λ, ε, Nₙ, Nₚ)
        if me′.err < err_min
            err_min = me′.err
            me = me′
        end
    end

    return me
end

"""
    fit_sellmeier(λ::AbstractVector{<:Real}, ε::AbstractVector{<:Real}, Nₙ::Integer, Nₚ::Integer)

Fit the dielectric constant data `ε` sampled at wavelengths `λ` to the Sellmeier equation
with exactly `Nₙ` terms below and `Nₚ` terms above the range of `λ`.
"""
fit_sellmeier(λ::AbsVecReal, ε::AbsVecReal, Nₙ::Integer, Nₚ::Integer) = fit_sellmeier(float(λ), float(ε), Nₙ, Nₚ)
function fit_sellmeier(λ::AbsVecFloat,
                       ε::AbsVecFloat,
                       Nₙ::Integer,  # number of terms with poles below λ range
                       Nₚ::Integer)  # number of terms with poles above λ range
    issorted(λ) || @error "λ = $λ should be sorted."

    N = Nₙ + Nₚ

    # Normaize λ.
    λₙ, λₚ = λ[1], λ[end]
    λ₀ = √(λₙ*λₚ)
    λnorm = λ ./ λ₀
    Nλ = length(λ)

    # Create guess parameters.
    cₙ, cₚ = (1-GUESS_MARGIN)λnorm[1], (1+GUESS_MARGIN)λnorm[end]

    p₀mat = ones(N,2)
    p₀mat[1:Nₙ,2] .= cₙ
    p₀mat[Nₙ+1:end,2] .= cₚ
    p₀ = reshape(p₀mat, :)

    # Perform fitting.
    fit = curve_fit((λ,p)->(BC = reshape(p,:,2);
                            str = Ref(@view(BC[:,1]));
                            λres = Ref(@view(BC[:,2]));
                            sellmeier_model.(λ, str, λres)),
                    λnorm, ε, p₀)

    str = @view(fit.param[1:N])  # entries of str are dimensionless, so no need to unnormalize them
    λres = @view(fit.param[N+1:2N]) .* λ₀  # entries of λres are wavelengths, so unnormalize them
    err = sqrt(sum(abs2, fit.resid) / Nλ)  # RMS error

    mdl = SellmeierModel(str, λres)

    return (; mdl, err)
end


text:
const __bodyfunction__ = Dict{Method,Any}()

# Find keyword "body functions" (the function that contains the body
# as written by the developer, called after all missing keyword-arguments
# have been assigned values), in a manner that doesn't depend on
# gensymmed names.
# `mnokw` is the method that gets called when you invoke it without
# supplying any keywords.
function __lookup_kwbody__(mnokw::Method)
    function getsym(arg)
        isa(arg, Symbol) && return arg
        @assert isa(arg, GlobalRef)
        return arg.name
    end

    f = get(__bodyfunction__, mnokw, nothing)
    if f === nothing
        fmod = mnokw.module
        # The lowered code for `mnokw` should look like
        #   %1 = mkw(kwvalues..., #self#, args...)
        #        return %1
        # where `mkw` is the name of the "active" keyword body-function.
        ast = Base.uncompressed_ast(mnokw)
        if isa(ast, Core.CodeInfo) && length(ast.code) >= 2
            callexpr = ast.code[end-1]
            if isa(callexpr, Expr) && callexpr.head == :call
                fsym = callexpr.args[1]
                if isa(fsym, Symbol)
                    f = getfield(fmod, fsym)
                elseif isa(fsym, GlobalRef)
                    if fsym.mod === Core && fsym.name === :_apply
                        f = getfield(mnokw.module, getsym(callexpr.args[2]))
                    elseif fsym.mod === Core && fsym.name === :_apply_iterate
                        f = getfield(mnokw.module, getsym(callexpr.args[3]))
                    else
                        f = getfield(fsym.mod, fsym.name)
                    end
                else
                    f = missing
                end
            else
                f = missing
            end
        else
            f = missing
        end
        __bodyfunction__[mnokw] = f
    end
    return f
end

function _precompile_()
    ccall(:jl_generating_output, Cint, ()) == 1 || return nothing
    # widgets
    ## slider
    @assert precompile(Tuple{typeof(slider),UnitRange{Int}})
    @assert precompile(Tuple{Core.kwftype(typeof(slider)),NamedTuple{(:signal, :orientation),Tuple{Signal{Int},Char}},typeof(slider),UnitRange{Int}})
    @assert precompile(Tuple{typeof(push!),Slider{Int},UnitRange{Int},Int})
    @assert precompile(Tuple{typeof(push!),Slider{Float32},StepRangeLen{Float32,Float64,Float64},Float32})
    @assert precompile(Tuple{typeof(push!),Slider{Float64},StepRangeLen{Float64,Base.TwicePrecision{Float64},Base.TwicePrecision{Float64}},Float64})
    @assert precompile(Tuple{Type{Slider},Signal{Int},Gtk.GtkScaleLeaf,UInt,Array{Any,1}})
    @assert precompile(Tuple{Type{Slider},Signal{Float64},GtkScaleLeaf,UInt,Array{Any,1}})
    @assert precompile(Tuple{Type{Slider},Signal{Float32},GtkScaleLeaf,UInt,Array{Any,1}})
    let fbody = try __lookup_kwbody__(which(slider, (StepRangeLen{Float32,Float64,Float64},))) catch missing end
        if !ismissing(fbody)
            precompile(fbody, (Nothing,Nothing,Signal{Float32},String,Bool,Nothing,typeof(slider),StepRangeLen{Float32,Float64,Float64},))
            precompile(fbody, (Nothing,Nothing,Signal{Float64},String,Bool,Nothing,typeof(slider),StepRangeLen{Float64,Base.TwicePrecision{Float64},Base.TwicePrecision{Float64}},))
        end
    end
    ## checkbox
    @assert precompile(Tuple{typeof(checkbox)})
    @assert precompile(Tuple{typeof(checkbox),Bool})
    @assert precompile(Tuple{Core.kwftype(typeof(checkbox)),NamedTuple{(:label,),Tuple{String}},typeof(checkbox)})
    @assert precompile(Tuple{typeof(push!),Checkbox,Bool})
    ## togglebutton
    @assert precompile(Tuple{typeof(togglebutton)})
    @assert precompile(Tuple{typeof(togglebutton),Bool})
    @assert precompile(Tuple{Core.kwftype(typeof(togglebutton)),NamedTuple{(:label,),Tuple{String}},typeof(togglebutton)})
    ## button
    @assert precompile(Tuple{typeof(button)})
    @assert precompile(Tuple{typeof(button),String})
    @assert precompile(Tuple{Core.kwftype(typeof(button)),NamedTuple{(:widget,),Tuple{Gtk.GtkToolButtonLeaf}},typeof(button)})
    @assert precompile(Tuple{Core.kwftype(typeof(button)),NamedTuple{(:widget,),Tuple{Gtk.GtkButtonLeaf}},typeof(button)})
    ## spinbutton
    @assert precompile(Tuple{typeof(spinbutton),UnitRange{Int}})
    @assert precompile(Tuple{Type{SpinButton},Signal{Int},Gtk.GtkSpinButtonLeaf,UInt,Array{Any,1}})
    @assert precompile(Tuple{Core.kwftype(typeof(spinbutton)),NamedTuple{(:widget, :signal),Tuple{Gtk.GtkSpinButtonLeaf,Signal{Int}}},typeof(spinbutton),UnitRange{Int}})
    @assert precompile(Tuple{typeof(push!),SpinButton{Int},UnitRange{Int},Int})
    ## cyclicspinbutton
    @assert precompile(Tuple{typeof(cyclicspinbutton),UnitRange{Int},Signal{Bool}})
    @assert precompile(Tuple{Type{CyclicSpinButton},Signal{Int},Gtk.GtkSpinButtonLeaf,UInt,Array{Any,1}})
    ## textbox
    @assert precompile(Tuple{typeof(textbox),String})
    @assert precompile(Tuple{Type{Textbox},Signal{String},Gtk.GtkEntryLeaf,UInt,Array{Any,1},Nothing})
    @assert precompile(Tuple{Type{Textbox},Signal{Int},Gtk.GtkEntryLeaf,UInt,Array{Any,1},UnitRange{Int}})
    @assert precompile(Tuple{Type{Textbox},Signal{Float64},GtkEntryLeaf,UInt,Array{Any,1},Nothing})
    @assert precompile(Tuple{Type{Textbox},Signal{Float32},GtkEntryLeaf,UInt,Array{Any,1},Nothing})
    @assert precompile(Tuple{Core.kwftype(typeof(textbox)),NamedTuple{(:signal,),Tuple{Signal{Int}}},typeof(textbox),Type{Int}})
    @assert precompile(Tuple{Core.kwftype(typeof(textbox)),NamedTuple{(:range,),Tuple{UnitRange{Int}}},typeof(textbox),Int})
    @assert precompile(Tuple{Core.kwftype(typeof(textbox)),NamedTuple{(:widget, :signal, :range),Tuple{GtkEntryLeaf,Signal{Int},Base.OneTo{Int}}},typeof(textbox),Int})
    let fbody = try __lookup_kwbody__(which(textbox, (Type{Float32},))) catch missing end
        if !ismissing(fbody)
            precompile(fbody, (GtkEntryLeaf,Nothing,Nothing,Signal{Float32},Bool,Nothing,Symbol,typeof(textbox),Type{Float32},))
        end
    end
    ## textarea
    @assert precompile(Tuple{typeof(textarea),String})
    ## label
    @assert precompile(Tuple{typeof(label),String})
    ## dropdown
    @assert precompile(Tuple{typeof(dropdown),Array{Pair{String,Function},1}})
    @assert precompile(Tuple{Core.kwftype(typeof(dropdown)),NamedTuple{(:label,),Tuple{String}},typeof(dropdown),Array{Pair{String,Function},1}})
    ## progressbar
    @assert precompile(Tuple{typeof(progressbar),Interval{:closed,:closed,Int}})
    ## player
    @assert precompile(Tuple{typeof(player),UnitRange{Int}})
    @assert precompile(Tuple{typeof(player),Signal{Int},UnitRange{Int}})
    let fbody = try __lookup_kwbody__(which(player, (Signal{Int},Base.OneTo{Int},))) catch missing end
        if !ismissing(fbody)
            @assert precompile(fbody, (String,Int,typeof(player),Signal{Int},Base.OneTo{Int},))
        end
    end
    ## timewidget
    @assert precompile(Tuple{Core.kwftype(typeof(timewidget)),NamedTuple{(:signal,),Tuple{Signal{Time}}},typeof(timewidget),Time})
    @assert precompile(Tuple{typeof(push!),TimeWidget{Time},Time})
    @assert precompile(Tuple{typeof(push!),TimeWidget{DateTime},DateTime})
    ## datetimewidget
    @assert precompile(Tuple{Core.kwftype(typeof(datetimewidget)),NamedTuple{(:signal,),Tuple{Signal{DateTime}}},typeof(datetimewidget),DateTime})

    # canvas
    for unit in (DeviceUnit, UserUnit)
        @assert precompile(Tuple{typeof(canvas),Type{unit}})
        @assert precompile(Tuple{typeof(mousedown_cb),  Ptr{GObject},Ptr{Gtk.GdkEventButton},MouseHandler{unit}})
        @assert precompile(Tuple{typeof(mouseup_cb),    Ptr{GObject},Ptr{Gtk.GdkEventButton},MouseHandler{unit}})
        @assert precompile(Tuple{typeof(mousemove_cb),  Ptr{GObject},Ptr{Gtk.GdkEventScroll},MouseHandler{unit}})
        @assert precompile(Tuple{typeof(mousescroll_cb),Ptr{GObject},Ptr{Gtk.GdkEventScroll},MouseHandler{unit}})
        @assert precompile(Tuple{typeof(init_pan_drag),Canvas{unit},Signal{ZoomRegion{RInt}}})
        @assert precompile(Tuple{typeof(init_zoom_rubberband),Canvas{unit},Signal{ZoomRegion{RInt}}})
        @assert precompile(Tuple{typeof(init_zoom_scroll),Canvas{unit},Signal{ZoomRegion{RInt}}})
        @assert precompile(Tuple{typeof(init_pan_scroll),Canvas{unit},Signal{ZoomRegion{RInt}}})
        @assert precompile(Tuple{typeof(Reactive.run_push),Signal{MouseScroll{unit}},MouseScroll{unit},Function})
        @assert precompile(Tuple{typeof(draw),Function,Canvas{unit},Signal{Int},Vararg{Signal{Int},N} where N})
    end
    precompile(Tuple{Type{ZoomRegion},Tuple{Base.OneTo{Int},Base.OneTo{Int}}})
    precompile(Tuple{typeof(push!),Signal{ZoomRegion{RInt}},Tuple{UnitRange{Int},UnitRange{Int}}})
    # image_surface
    @assert precompile(Tuple{typeof(image_surface),Matrix{RGB{N0f8}}})
    @assert precompile(Tuple{typeof(image_surface),Matrix{RGBA{N0f8}}})
    @assert precompile(Tuple{typeof(image_surface),Matrix{Gray{N0f8}}})
    # map
    # Note: @asserts on some of these cause problems on earlier versions of Julia
    precompile(Tuple{typeof(map),Function,Signal{ZoomRegion{RInt}}})
    @assert precompile(Tuple{typeof(map),Function,Textbox{String},Textbox{Int}})
    precompile(Tuple{typeof(map),Function,Button})
    for unit in (DeviceUnit, UserUnit)
        precompile(Tuple{typeof(map),Function,Signal{MouseButton{unit}}})
        precompile(Tuple{typeof(map),Function,Signal{MouseScroll{unit}}})
    end
    # gc_preserve
    @assert precompile(Tuple{typeof(gc_preserve),GtkWindowLeaf,Tuple{Dict{String,Any},Signal{Nothing},Signal{Nothing}}})
    @assert precompile(Tuple{typeof(gc_preserve),GtkWindowLeaf,Dict{String,Dict{String,Any}}})
    @assert precompile(Tuple{typeof(gc_preserve),GtkFrameLeaf,Player{PlayerWithTextbox}})
    # init_signal2widget
    @assert precompile(Tuple{typeof(init_signal2widget),GtkScaleLeaf,UInt,Signal{Int}})
    @assert precompile(Tuple{typeof(init_signal2widget),GtkScaleLeaf,UInt,Signal{Float64}})
    @assert precompile(Tuple{typeof(init_signal2widget),GtkScaleLeaf,UInt,Signal{Float32}})

    # anonymous functions
    # (many of these won't work, but worth a shot)
    isdefined(GtkReactive, Symbol("#104#113")) && precompile(Tuple{getfield(GtkReactive, Symbol("#104#113")),Int})
    isdefined(GtkReactive, Symbol("#126#144")) && precompile(Tuple{getfield(GtkReactive, Symbol("#126#144")),Int})
    isdefined(GtkReactive, Symbol("#106#115")) && precompile(Tuple{getfield(GtkReactive, Symbol("#106#115")),Int})
    isdefined(GtkReactive, Symbol("#108#117")) && precompile(Tuple{getfield(GtkReactive, Symbol("#108#117")),Int})
    isdefined(GtkReactive, Symbol("#132#150")) && precompile(Tuple{getfield(GtkReactive, Symbol("#132#150")),Int})
    isdefined(GtkReactive, Symbol("#136#154")) && precompile(Tuple{getfield(GtkReactive, Symbol("#136#154")),Int})
    isdefined(GtkReactive, Symbol("#130#148")) && precompile(Tuple{getfield(GtkReactive, Symbol("#130#148")),Int})
    isdefined(GtkReactive, Symbol("#134#152")) && precompile(Tuple{getfield(GtkReactive, Symbol("#134#152")),Int})
    isdefined(GtkReactive, Symbol("#128#146")) && precompile(Tuple{getfield(GtkReactive, Symbol("#128#146")),Int})
    isdefined(GtkReactive, Symbol("#159#161")) && precompile(Tuple{getfield(GtkReactive, Symbol("#159#161")),Gtk.GtkCanvas})
end

ext: jl


text:
using TrajectoryOptimization
using StaticArrays
using LinearAlgebra
using Random
using Test
using Distributions
const TO = TrajectoryOptimization

# Solve Quadrotor zig-zag problem
prob,opts = Problems.Quadrotor(:zigzag, costfun=:Quadratic)
solver = iLQRSolver(prob,opts)
solve!(solver)
iterations(solver)
@test norm(state(solver.Z[end])[1:3] - solver.xf[1:3]) < 0.1

# Convert trajectory to general RBState
Xref = Controllers.RBState(solver.model, solver.Z)
Uref = controls(solver)
push!(Uref, Dynamics.trim_controls(solver.model))
dt_ref = solver.Z[1].dt

# Track with different model
Rot = MRP{Float64}
model = Dynamics.Quadrotor2{Rot}(use_rot=true, mass=.5)
dt = 1e-4
Nsim = Int(solver.tf/dt) + 1
inds = Int.(range(1,Nsim,length=solver.N))
Q = Diagonal(Dynamics.fill_error_state(model, 200, 200, 50, 50))
R = Diagonal(@SVector fill(1.0, 4))

mlqr = Controllers.TVLQR(model, Q, R, Xref, Uref, dt_ref)

noise = Distributions.MvNormal(Diagonal(fill(5.,6)))
model_ = Dynamics.NoisyRB(model, noise)
X = Controllers.simulate(model_, mlqr, Xref[1], solver.tf, w=0.0)
res_lqr = [Controllers.RBState(model, x) for x in X[inds]]
e1 = norm(res_lqr .⊖ Xref)

# Track with SE3 controller

# get reference trajectories
b_ = map(Xref) do x
    # prefer pointing torwards the goal
    x.r - Xref[end].r + 0.01*@SVector rand(3)
end
X_ = map(Xref) do x
    Dynamics.build_state(model, x)
end
Xd_ = map(1:solver.N) do k
    dynamics(model, X_[k], Uref[k])
end
t = collect(0:dt_ref:solver.tf)

# cntrl = SE3Tracking(model, X_, Xd_, b_, t)
kx = 2.71 + 10.0
kv = 1.01 + 3.0
kR = 3.0*Diagonal(@SVector [1,1,1.])
kO = 0.5*Diagonal(@SVector [1,1,1.])
hfca = Controllers.HFCA(model, X_, Xd_, b_, t, kx=kx, kv=kv, kR=kR, kO=kO)

X2 = Controllers.simulate(model, hfca, Xref[1], solver.tf, w=0.0)
res_so3 = [Controllers.RBState(model, x) for x in X2[inds]]
e2 = norm(res_so3 .⊖ Xref)
@test e1 < e2

# Test MLQR
model = Dynamics.Quadrotor2()
Q_ = (200.,200,50,50)
R = Diagonal(@SVector fill(1.0,4))
xref = zeros(model)[1]
uref = Dynamics.trim_controls(model)
dt_ctrl = 0.01
x0 = Controllers.RBState([0,0,0.2], expm(deg2rad(80)*@SVector [1,0,0]), [0,0,1.], [0.5,0.5,0])
tf = 10.0

Q = Diagonal(Dynamics.fill_error_state(model, Q_...))
cntrl = Controllers.MLQR(model, dt_ctrl, Q, R, xref, uref)
xinit = Dynamics.build_state(model, x0)
res = Controllers.simulate(model, cntrl, xinit, tf, dt=dt, w=0.0)
err = TrajectoryOptimization.state_diff(model, res[end], xref)
@test norm(err) < 0.01

# Test LQR
model = Dynamics.Quadrotor2{RodriguesParam{Float64}}(use_rot=false)
Q = Diagonal(Dynamics.fill_error_state(model, Q_...))
xref = zeros(model)[1]
cntrl = Controllers.LQR(model, dt_ctrl, Q, R, xref, uref)
xinit = Dynamics.build_state(model, x0)
res = Controllers.simulate(model, cntrl, xinit, tf, dt=dt, w=0.0)

# Test HFCA
Rot = UnitQuaternion{Float64,IdentityMap}
model = Dynamics.Quadrotor2{Rot}(use_rot=false)

times = range(0,tf,step=dt)
xref = zeros(model)[1]
Xref = [copy(xref) for k = 1:length(times)]
bref = [@SVector [1,0,0.] for k = 1:length(times)]
Xdref = [copy(xref) for k = 1:length(times)]

cntrl = Controllers.HFCA(model, Xref, Xdref, bref, collect(times))
xinit = Dynamics.build_state(model, x0)
res = Controllers.simulate(model, cntrl, xinit, tf, dt=dt, w=0.0)
err = TrajectoryOptimization.state_diff(model, res[end], xref)
@test norm(err) < 0.01

# Test SE3
cntrl = Controllers.SE3Tracking(model, Xref, Xdref, bref, collect(times))
xinit = Dynamics.build_state(model, x0)
res = Controllers.simulate(model, cntrl, xinit, tf, dt=dt, w=0.0)
err = TrajectoryOptimization.state_diff(model, res[end], xref)
@test norm(err) < 0.01


text:
prepend!(LOAD_PATH, ["Project.toml"])

using Dates
import HTTP

years = collect(1992:2015)
months = collect(1:12)

println("years=", years)
println("months=", months)

make_url(starttime, endtime) = "https://earthquake.usgs.gov/fdsnws/event/1/query?format=csv&starttime=$(starttime)&endtime=$(endtime)&minmagnitude=2&minlongitude=-170&minlatitude=15&maxlongitude=-60&maxlatitude=72"

base = "./data/usgs_earthquake"
mkpath(base)

for year in years
    path = joinpath(base, "$(year).csv")
    open(path, "w") do file
        for month in months
            dt1 = DateTime(year, month, 1)
            dt2 = lastdayofmonth(dt1)
            format = DateFormat("yyyy-m-d")
            starttime = Dates.format(dt1, format)
            endtime = Dates.format(dt2, format)
            println("starttime=$(starttime)&endtime=$(endtime)")
            url =  make_url(starttime, endtime)
            # println(url)
            response = HTTP.request("GET", url)
            # println(String(response.body))
            content = String(response.body)
            if month != 1
                content = content[findfirst("\n", content).start+1:lastindex(content)]
            end
            write(file, content)
        end
        run(`wc -l $(path)`)
    end
end


text:
# # One-Dimensional Cases
#
# The 1D case in Optimal Transport is a special case where one can
# easily obtain closed form solutions efficiently
# when the cost function is convex. In this situation,
# one does no need to use Linear Programming solvers
# to obtain the exact solution to the problem.

# ## Packages
#
# We load the following packages into our environment:

using OptimalTransport
using Distances
using Distributions
using StatsPlots

using LinearAlgebra
using Random

Random.seed!(1234);

# ## Continuous Distribution
#
# In the 1D case, when the source measure $\mu$ is continuous and the cost function
# has the form ``c(x, y) = h(|x - y|)`` where ``h`` is a convex function,
# the optimal transport plan is the Monge map
# ```math
# T = F_\nu^{-1} \circ F_\mu
# ```
# where ``F_\mu`` is the cumulative distribution function of `μ` and ``F_\nu^{-1}`` is the
# quantile function of `ν`.
# In this setting, the optimal transport cost can be computed as
# ```math
# \int_0^1 c(F_\mu^{-1}(x), F_\nu^{-1}(x)) \mathrm{d}x
# ```
# where ``F_\mu^{-1}`` and ``F_\nu^{-1}`` are the quantile functions of `μ` and `ν`,
# respectively.
#
# We start by defining the distributions.

μ = Normal(0, 1)

N = 10
ν = Poisson(N);

# Nest, we define a cost function.
c(x, y) = (abs(x - y))^2 # could have used `sqeuclidean` from `Distances.jl`

T = ot_plan(c, μ, ν);

# `T` is the Monge Map. Let's visualize it.
p1 = plot(μ; label='μ')
p1 = plot!(ν; marker=:circle, label='ν')
p2 = plot(-2:0.1:2, T(-2:0.1:2); label="Monge map", color=:green, legend=:topleft)
plot(p1, p2)

# The optimal transport cost can be computed with
ot_cost(c, μ, ν)

# If instead you want the 2-Wasserstein distance (which is the square root
# of the optimal transport with the Square Euclidean distatce, then use
wasserstein(μ, ν; p=2)

# ## Finite Discrete Distributions
#
# If the source and target measures are 1D finite discrete distributions
# (sometimes referred as empirical distributions, or as sample distributions),
# and if the cost function is convex, then the optimal
# transport plan can be written as a sorting algorithm,
# where the utmost left probability mass of the source is transported
# to the closest probability mass of the target, until everything is transported.
#
# Define your measures as DiscreteNonParametric, which is a type in Distributions.jl.
# Also, let's assume both point masses with equal weights and let's
# use the `sqeuclidean` function instead of creating our own cost function.

M = 15
μ = DiscreteNonParametric(1.5rand(M), fill(1 / M, M))

N = 10
ν = DiscreteNonParametric(1.5rand(N) .+ 2, fill(1 / N, N))

γ = ot_plan(sqeuclidean, μ, ν);

# This time γ is a sparse matrix containing the transport plan. Let's visualize the results.
# We create a function `curve` just as a helper to draw the transport plan.

function curve(x1, x2, y1, y2)
    a = min(y1, y2)
    b = (y1 - y2 + a * (x1^2 - x2^2)) / (x1 - x2)
    c = y1 + a * x1^2 - b * x1
    f(x) = -a * x^2 + b * x + c
    return f
end

p = plot(μ; marker=:circle, label='μ')
p = plot!(ν; marker=:circle, label='ν', ylims=(0, 0.2))
for i in 1:M, j in 1:N
    if γ[i, j] > 0
        transport = curve(μ.support[i], ν.support[j], 1 / M, 1 / N)
        x = range(μ.support[i], ν.support[j]; length=100)
        p = plot!(x, transport.(x); color=:green, label=nothing, alpha=0.5)
    end
end
p

# Again, the optimal transport cost can be calculated with

ot_cost(sqeuclidean, μ, ν)


text:
module Czoo

const LIB="src/libczoo"

export free_ptr, add, concat, cons, CPstring, Pstring, pstring, ptr_pstring, CLinkedPstring, LinkedPstring, linked_pstrings
export print_list_int, print_Pstruct, Pstruct, print_Pstring_as_Cstring, capture_c_stdout_to_file, capture_c_stdout_to_IOBuffer

"""
    free_ptr(ptr; free=true)
Free the given C ptr in C.
# Arguments
- `ptr` the C pointer to be freed
- `free` a conditional flag to prevent an `if` at the call site
"""
free_ptr(ptr; free=true) = if free ccall(:free, Cvoid, (Ptr{Cvoid},), Ptr{Cvoid}(ptr)) end

"""
    string_from_ptr_uchar(ptr_uchar; free=true)
Create a Julia String from a memory location a byte at a time.
# Arguments
- `ptr_uchar` pointer to use
- `free` flag to determine if the memory should be freed once copied
"""
function string_from_ptr_uchar(ptr_uchar; free=true)
    if ptr_uchar == C_NULL
        return ""
    end
    uchars = Vector{UInt8}()
    i = 1
    while (c = unsafe_load(ptr_uchar, i)) != 0
        push!(uchars, c)
        i += 1
    end
    free_ptr(ptr_uchar; free)
    String(uchars)
end

"""
    vec_string_from_ptr_ptr_uchar(ptr_ptr_uchar; free=true)
Create a Vector{String} from the supplied list of uchar*
# Arguments
- `ptr_ptr_uchar` the pointer to the list
- `free` flag to determine if the memory of both the strings and the list should be freed once copied
"""
function vec_string_from_ptr_ptr_uchar(ptr_ptr_uchar; free=true)
    if ptr_ptr_uchar == C_NULL
        return String[]
    end
    strings = Vector{String}()
    i = 1
    while (stringp = unsafe_load(ptr_ptr_uchar,i)) != C_NULL
        push!(strings, string_from_ptr_uchar(stringp; free))
        i += 1
    end
    free_ptr(ptr_ptr_uchar; free)
    strings
end

"""
    CPstring
A "Pascal" style string with a length and the uchars
# Properties
- `length::Cint` length of the uchars
- `uchars::Ptr{Cuchar}` pointer to the uchars
"""
struct CPstring
    length::Cint
    uchars::Ptr{Cuchar}
end

"""
    Pstring
The Julia version of the `CPstring`
# Properties
- `length::Int` length of the string
- `uchars::String` the string
# Constructor
- `Pstring(cp::CPstring; free=true)` - create a Pstring from a CPstring and optionally free the memory
"""
struct Pstring
    length::Int
    uchars::String
    Pstring(s::String) = new(length(s), s)
    Pstring(cp::CPstring; free=true) = new(Int(cp.length), string_from_ptr_uchar(cp.uchars; free))
    function Pstring(pcp::Ptr{CPstring}; free=true)
        if pcp == C_NULL
            new(0, "")
        else
            cp = unsafe_load(pcp, 1)
            new(Int(cp.length), string_from_ptr_uchar(cp.uchars; free))
        end
    end
end

"""
    CLinkedPstring
A "Pascal" style string which also links to the "next" in the list
# Properties
- `next::Ptr{CLinkedPstring}`
- `length::Cint`
- `uchars::Ptr{Cuchar}`
"""
struct CLinkedPstring
    next::Ptr{CLinkedPstring}
    length::Cint
    uchars::Ptr{Cuchar}
end

"""
    LinkedPstring
The Julia version of the CLinkedPstring
# Constructor
- `LinkedPstring(clp::CLinkedPstring; free=true)` - create a LinkedPstring (following all the links) from a CLinkedPstring and optionally free the memory
- `LinkedPstring(pclp::Ptr{CLinkedPstring}; free=true)` - create a LinkedPstring (following all the links) from a Ptr{CLinkedPstring} and optionally free the memory

"""
struct LinkedPstring
    next
    length
    uchars
    LinkedPstring(clp::CLinkedPstring; free=true) = new(clp.next == C_NULL ? nothing : LinkedPstring(clp.next), clp.length, string_from_ptr_uchar(clp.uchars; free))
    function LinkedPstring(pclp::Ptr{CLinkedPstring}; free=true)
        if pclp == C_NULL
            return new(nothing, 0, "")
        end
        lp = LinkedPstring(unsafe_load(pclp,1); free)
        free_ptr(pclp; free)
        lp
    end
    LinkedPstring(s::AbstractString) = new(nothing, length(s), s)
    LinkedPstring(lp, s::AbstractString) = new(lp, length(s), s)
end

Base.isequal(lp1::LinkedPstring, lp2::LinkedPstring) = isequal(lp1.next, lp2.next) && lp1.length == lp2.length && lp1.uchars == lp2.uchars

"""
    add(a, b)
Calls C to adds two Ints and returns the result
"""
add(a, b) = @ccall LIB.add(a::Int, b::Int)::Int

"""
    concat(a, b)
Calls C to concat two Strings and returns the result (and frees the memory)
"""
concat(a, b) = (@ccall LIB.concat(a::Cstring, b::Cstring)::Ptr{Cuchar}) |> string_from_ptr_uchar

"""
    cons(a, b)
Calls C to create a list of two Strings and returns the result as a Vector{String} (and frees the memory)
"""
cons(a, b) = (@ccall LIB.cons(a::Cstring, b::Cstring)::Ptr{Ptr{Cuchar}}) |> vec_string_from_ptr_ptr_uchar

"""
    pstring(a)
Calls C to create a CPstring struct and returns a PString (and frees the memory used for the string)
"""
pstring(a) = Pstring(@ccall LIB.pstring(a::Cstring)::CPstring)

"""
    ptr_pstring(a)
Calls C to create a Ptr{CPstring} struct and returns a PString (and frees the memory of the string and the malloc'd struct)
"""
ptr_pstring(a) = Pstring((@ccall LIB.ptr_pstring(a::Cstring)::Ptr{CPstring}))

"""
    linked_pstrings(a, b, c)
Calls C to create a Ptr{CLinkedPstring} struct and returns a LinkedPString (and frees the memory of the string and the malloc'd structs)
"""
linked_pstrings(a, b, c) = LinkedPstring((@ccall LIB.linked_pstrings(a::Cstring, b::Cstring, c::Cstring)::Ptr{CLinkedPstring}))

"""
    print_list_ints(ints::Vector{Cint})
Calls C to print a list of Ints
"""
print_list_int(ints::Vector{Cint}) = @ccall LIB.print_list_int(ints::Ptr{Cint}, length(ints)::Cint)::Cint

"""
    print_Pstring_as_Cstring(ps::Pstring)
Send a Pstring Julia struct as a converted Cstring struct. I imagine there is a way to do this without the conversion.
"""
function print_Pstring_as_Cstring(ps::Pstring)
    cp = CPstring(ps.length, pointer(ps.uchars))
    @ccall LIB.print_Pstring(cp::CPstring)::Cint
end

"""
    capture_c_stdout_to_file()
Redirect stdout of the C program to a file and print it in Julia.
"""
function capture_c_stdout_to_file(c::Cint)
    (path, io) = mktemp() # Julia will delete this on exit
    prev_stdout = stdout
    redirect_stdout(io) do
        ccall(:putchar, Cint, (Cint,), c)
        Base.Libc.flush_cstdio() # required or it sits in the buffer - even between invocations
    end
    seekstart(io)
    redirect_stdout(prev_stdout)
    captured = "Captured: $(String(read(io)))"
    println(captured)
    captured
end

"""
    capture_c_stdout_to_IOBuffer(c::Cint)
Call a C function and capture its stdout into an IOBuffer
Thanks to Carsten Bauer https://github.com/carstenbauer for this one
"""
function capture_c_stdout_to_IOBuffer(c::Cint)

    function call()
        ccall(:putchar, Cint, (Cint,), c)
        Base.Libc.flush_cstdio() # required or it sits in the buffer - even between invocations
        flush(stdout) # also needed
    end

    io = IOBuffer()
    prev_stdout = stdout
    rd, = redirect_stdout()
    task = @async write(io, rd)
    try
        call()
    finally
        close(rd)
        redirect_stdout(prev_stdout)
        wait(task)
    end

    seekstart(io)
    captured = "Captured: $(String(read(io)))"
    println(captured)
    captured
end

#==
for adding to the collection

puts(s::Union{String,SubString{String}}) = ccall(:puts, Cint, (Ptr{Cchar},), s)
jl_malloc(n::Integer) = ccall(:jl_malloc, Ptr{Cvoid}, (Csize_t,), n)

==#

#== These don't work yet / are untested

"""
    arity0_func()
A function without arguments or return type to send to C to call
"""
arity0_func() = println("J: Airty0 function")

"""
    send_arity0_func()
Create a pointer to a Julia function which takes no arguments and returns nothing.
Call a function C which will, in turn, call that Julia function.
"""
function send_arity0_func()
    jfunc = @cfunction(arity0_func, Ptr{Cvoid}, (Cvoid,))
    @ccall LIB.call_arity0_julia_func(jfunc::Ptr{Cvoid})::Ptr{Cvoid}
end

"""
    airity1_func(i::Cint)
An airity1 function for C to call with a C Int, it just prints it.
"""
airity1_func(i::Cint) = prinlnt("J: C called with $i")

"""
    send_airity1_int(i::Cint)
Create a pointer to a Julia function which takes a C Integer but returns nothing.
Call a function C which will, in turn, call that Julia function.
"""
function send_airity1_int(i::Cint)
    jfunc = @cfunction(airity1_func, Cvoid, (Cvoid,))
    @ccall LIB.call_airity1_julia_func(jfunc::Ptr{Cvoid})::Ptr{Cvoid}
end

"""
    airity2_func(a::Cint, b::Cint)::Cint
Call this function with two Cint arguments and get a Cint in return
"""
    airity2_func(a::Cint, b::Cint)::Cint = begin println("J: called with $a $b, returning $(a+b)"); a+b; end

"""
    send_airity2_func()
"""
function send_airity2_func()
    jfunc = @cfunction(airity2_func, Cint, (Cint, Cint))
    i = @ccall LIB.call_airity2_julia_func(jfunc::Ptr{Cvoid})::Cint
    println("J: Called C, got back $i")
end

==#

###
end


text:
module Rendering
	using Images

	struct PlotOpts
		width::Int
		height::Int
		maxValue::Float64
		pixelsPerUnit::Int
	end

	function GenPlaneImage(values::Array{Float64}, opts::PlotOpts)
		img = Array{RGB}(undef, opts.height, opts.width)

		for c in range(1, opts.width)
			for r in range(1, opts.height)
				val = clamp01nan(abs(values[c,r]) / opts.maxValue)

				img[r, c] = RGB(values[c,r] > 0 ? val : 0, 0, values[c,r] < 0 ? val : 0)
			end
		end

		save("output.png", img)
	end

	function GenFuncImage(func::Function, opts::PlotOpts)
		img = Array{RGB}(undef, opts.height, opts.width)

		for c in range(1, opts.width)
			for r in range(1, opts.height)
				funcVal = func(c / opts.pixelsPerUnit, r / opts.pixelsPerUnit)
				val = clamp01nan(abs(funcVal) / opts.maxValue)

				img[r, c] = RGB(funcVal > 0 ? val : 0, 0, funcVal < 0 ? val : 0)
			end
		end

		save("output.png", img)
	end

	export GenPlaneImage, GenFuncImage
	export PlotOpts
end


text:
module WSMP

## Functions to interface with the Watson Sparse Matrix Package, WSMP

## Note that WSMP is not free software.  An evaluation version is available, see
##  http://researcher.watson.ibm.com/projects/wsmp
## Production use requires the purchase of a license.

if isfile(joinpath(Pkg.dir("WSMP"),"deps","deps.jl"))
    include("../deps/deps.jl")
else
    error("WSMP not properly installed. Please run Pkg.build(\"WSMP\")")
end

ENV["MALLOC_TRIM_THRESHOLD_"] = -1
ENV["MALLOC_MMAP_MAX_"] = 0
ENV["WSMPLICPATH"] = Pkg.dir("WSMP","deps","usr","lib")
blas_set_num_threads(1)

export Wssmp,
    wsafree,                   # release the internal copy of the system matrix
    wsclear,                   # release all internal storage
    wsffree,                   # release the internal numeric factor storage
    wsrecallmat,               # recall saved matrix k as active storage
    wssetlf,                   # set the load factor (0. < lf < 1.)
    wssetmaxstack,             # set the max stack size (relative to 1 Mb stack)
    wssetmaxthreads,           # set maximum number of threads to spawn
    wssfree,                   # release the internal symbolic factor storage
    wssmp,                     # call wssmp_ on matrix A for stage k
    wssortindicesd,            # sort row indices carrying along the non-zero values
    wssortindicesi,            # sort row indices
    wsstoremat,                # store the current active matrix in slot k (0 ≤ k ≤ 63)
    wsversion                  # report the version of WSMP in use (returns 3 Cint's)

include("misc.jl")
include("wssmp.jl")

end # module


text:
using LabeledArrays
using Test

@testset "LabeledArrays.jl" begin
    # Write your own tests here.
    @testset "NBijection" begin
        x = NBijection('a':'z')
        @test x[13] == 'm'
        @test eltype(x) == Char
        @test get!(x,'m') == 13
        @test get!(x,'ö') == 27
        @test collect(x) == ['a':'z'...,'ö']
        @test keytype(x) == Char
        @test keys(x) == vcat('a':'z',['ö'])
        @test valtype(x) == Int
        @test values(x) == 1:27
        @test indexin('m', x) == indexin('m', collect(x))
    end
end


text:
## efficient value-based hashing of integers ##

function hash_integer(n::Integer, h::Uint)
    h = hash_uint(uint(n & typemax(Uint)) $ h) $ h
    n = ifelse(n < 0, oftype(n,-n), n)
    n >>>= sizeof(Uint) << 3
    while n != 0
        h = hash_uint(uint(n & typemax(Uint)) $ h) $ h
        n >>>= sizeof(Uint) << 3
    end
    return h
end

function hash_integer(n::BigInt, h::Uint)
    s = n.size
    s == 0 && return hash_integer(0, h)
    p = convert(Ptr{Uint}, n.d)
    b = unsafe_load(p)
    h = hash_uint(ifelse(s < 0, -b, b) $ h) $ h
    for k = 2:abs(s)
        h = hash_uint(unsafe_load(p, k) $ h) $ h
    end
    return h
end

## generic hashing for rational values ##

function hash(x::Real, h::Uint)
    # decompose x as num*2^pow/den
    num, pow, den = decompose(x)

    # handle special values
    num == 0 && den == 0 && return hash(NaN, h)
    num == 0 && return hash(ifelse(den > 0, 0.0, -0.0), h)
    den == 0 && return hash(ifelse(num > 0, Inf, -Inf), h)

    # normalize decomposition
    if den < 0
        num = -num
        den = -den
    end
    z = trailing_zeros(num)
    if z != 0
        num >>= z
        pow += z
    end
    z = trailing_zeros(den)
    if z != 0
        den >>= z
        pow -= z
    end

    # handle values representable as Int64, Uint64, Float64
    if den == 1
        left = ndigits0z(num,2) + pow
        right = trailing_zeros(num) + pow
        if -1074 <= right
            if 0 <= right && left <= 64
                left <= 63                     && return hash(int64(num) << int(pow), h)
                signbit(num) == signbit(den)   && return hash(uint64(num) << int(pow), h)
            end # typemin(Int64) handled by Float64 case
            left <= 1024 && left - right <= 53 && return hash(float64(num) * 2.0^pow, h)
        end
    end

    # handle generic rational values
    h = hash_integer(den, h)
    h = hash_integer(pow, h)
    h = hash_integer(num, h)
    return h
end

#=
`decompose(x)`: non-canonical decomposition of rational values as `num*2^pow/den`.

The decompose function is the point where rational-valued numeric types that support
hashing hook into the hashing protocol. `decompose(x)` should return three integer
values `num, pow, den`, such that the value of `x` is mathematically equal to

    num*2^pow/den

The decomposition need not be canonical in the sense that it just needs to be *some*
way to express `x` in this form, not any particular way – with the restriction that
`num` and `den` may not share any odd common factors. They may, however, have powers
of two in common – the generic hashing code will normalize those as necessary.

Special values:

 - `x` is zero: `num` should be zero and `den` should have the same sign as `x`
 - `x` is infinite: `den` should be zero and `num` should have the same sign as `x`
 - `x` is not a number: `num` and `den` should both be zero
=#

decompose(x::Integer) = x, 0, 1
decompose(x::Rational) = num(x), 0, den(x)

function decompose(x::Float32)
    isnan(x) && return 0, 0, 0
    isinf(x) && return ifelse(x < 0, -1, 1), 0, 0
    n = reinterpret(Int32, x)
    s = int32(n & 0x007fffff)
    e = int32(n & 0x7f800000 >> 23)
    s |= int32(e != 0) << 23
    d = ifelse(signbit(n) == 1, -1, 1)
    int(s), int(e - 150 + (e == 0)), d
end

function decompose(x::Float64)
    isnan(x) && return 0, 0, 0
    isinf(x) && return ifelse(x < 0, -1, 1), 0, 0
    n = reinterpret(Int64, x)
    s = int64(n & 0x000fffffffffffff)
    e = int64(n & 0x7ff0000000000000 >> 52)
    s |= int64(e != 0) << 52
    d = ifelse(signbit(n) == 1, -1, 1)
    int(s), int(e - 1075 + (e == 0)), d
end

function decompose(x::BigFloat)
    isnan(x) && return big(0), 0, 0
    isinf(x) && return big(x.sign), 0, 0
    x == 0 && return big(0), 0, int(x.sign)
    s = BigInt()
    ccall((:__gmpz_realloc2, :libgmp), Void, (Ptr{BigInt}, Culong), &s, x.prec)
    s.size = -fld(-x.prec,(sizeof(Culong)<<3))
    ccall(:memcpy, Ptr{Void}, (Ptr{Void}, Ptr{Void}, Csize_t), s.d, x.d, s.size*sizeof(Culong))
    s, int(x.exp - x.prec), int(x.sign)
end

## streamlined hashing for smallish rational types ##

function hash{T<:Integer64}(x::Rational{T}, h::Uint)
    num, den = Base.num(x), Base.den(x)
    den == 1 && return hash(num, h)
    den == 0 && return hash(ifelse(num > 0, Inf, -Inf), h)
    if isodd(den)
        pow = trailing_zeros(num)
        num >>= pow
    else
        pow = trailing_zeros(den)
        den >>= pow
        pow = -pow
        if den == 1 && abs(num) < 9007199254740992
            return hash(float64(num) * 2.0^pow)
        end
    end
    h = hash_integer(den, h)
    h = hash_integer(pow, h)
    h = hash_integer(num, h)
    return h
end

## hashing Float16s ##

hash(x::Float16, h::Uint) = hash(float64(x), h)

## hashing strings ##

const memhash = Uint == Uint64 ? :memhash_seed : :memhash32_seed

function hash{T<:ByteString}(s::Union(T,SubString{T}), h::Uint)
    h += uint(0x71e729fd56419c81)
    ccall(memhash, Uint, (Ptr{Uint8}, Csize_t, Uint32), pointer(s), sizeof(s), h) + h
end
hash(s::String, h::Uint) = hash(bytestring(s), h)

## hashing collections ##

function hash(a::AbstractArray, h::Uint)
    h += uint(0x7f53e68ceb575e76)
    h += hash(size(a))
    for x in a
        h = hash(x, h)
    end
    return h
end

function hash(a::Associative, h::Uint)
    h += uint(0x6d35bb51952d5539)
    for (k,v) in a
        h $= hash(k, hash(v))
    end
    return h
end

function hash(s::Set, h::Uint)
    h += uint(0x852ada37cfe8e0ce)
    for x in s
        h $= hash(x)
    end
    return h
end

hash(::(), h::Uint) = h + uint(0x77cfa1eef01bca90)
hash(x::(Any,), h::Uint)    = hash(x[1], hash((), h))
hash(x::(Any,Any), h::Uint) = hash(x[1], hash(x[2], hash((), h)))
hash(x::Tuple, h::Uint)     = hash(x[1], hash(x[2], hash(tupletail(x), h)))

hash(r::Range{Bool}, h::Uint) = invoke(hash, (Range, Uint), r, h)
hash(B::BitArray, h::Uint) = hash((size(B),B.chunks), h)
hash(a::AbstractArray{Bool}, h::Uint) = hash(bitpack(a), h)

# hashing ranges by component at worst leads to collisions for very similar ranges
function hash{T<:Range}(r::T, h::Uint)
    h += uint(0x80707b6821b70087)
    h = hash(first(r), h)
    h = hash(step(r), h)
    h = hash(last(r), h)
end

## hashing general objects ##

hash(x::ANY,  h::Uint) = hash(object_id(x), h)

max_stars_repo_name: aplund/julia

max_issues_repo_name: aplund/julia

max_forks_repo_name: aplund/julia


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder, Pkg

name = "PMIx"
version = v"4.1.0"

# Collection of sources required to complete build
sources = [
    ArchiveSource("https://github.com/openpmix/openpmix/releases/download/v$(version)/pmix-$(version).tar.bz2", "145f05a6c621bfb3fc434776b615d7e6d53260cc9ba340a01f55b383e07c842e")
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir
cd pmix-*
./configure --prefix=${prefix} --build=${MACHTYPE} --host=${target} \
    --enable-shared \
    --with-libevent=${prefix} \
    --with-hwloc=${prefix} \
    --without-tests-examples \
    --disable-man-pages
make -j${nproc}
make install
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = [
    Platform("i686", "linux"; libc = "glibc"),
    Platform("x86_64", "linux"; libc = "glibc"),
    Platform("aarch64", "linux"; libc = "glibc"),
    Platform("armv7l", "linux"; call_abi = "eabihf", libc = "glibc"),
    Platform("powerpc64le", "linux"; libc = "glibc"),
    Platform("i686", "linux"; libc = "musl"),
    Platform("x86_64", "linux"; libc = "musl"),
    Platform("aarch64", "linux"; libc = "musl"),
    Platform("armv7l", "linux"; call_abi = "eabihf", libc = "musl"),
    Platform("x86_64", "macos"; ),
    Platform("aarch64", "macos"; )
]

# TODO: Configure fails on Windows with:
```
checking for library containing event_config_new... no
checking for event_getcode4name in -levent... no
checking will libevent support be built... no
configure: WARNING: Either libevent or libev support is required, but neither
configure: WARNING: was found. Please use the configure options to point us
configure: WARNING: to where we can find one or the other library
configure: error: Cannot continue
```

# The products that we will ensure are always built
products = [
    LibraryProduct("libpmix", :libpmix)
    ExecutableProduct("pmix_info", :pmix_info)
]

# Dependencies that must be installed before this package can be built
dependencies = [
    Dependency(PackageSpec(name="libevent_jll", uuid="1080aeaf-3a6a-583e-a51c-c537b09f60ec")),
    Dependency(PackageSpec(name="Hwloc_jll", uuid="e33a78d0-f292-5ffc-b300-72abe9b543c8")),
    Dependency("Zlib_jll"),
]

init_block = raw"""
ENV["PMIX_PREFIX"] = artifact_dir
"""

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies;
               julia_compat="1.6", init_block=init_block)

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil


text:
"""
    NWFn(g0,m=0)

Calculates covariance matrix of sqrt(T)*sample average.

# Usage
S = NWFn(g0,m)

# Input
- `g0::Array`: Txq array of q moment conditions
- `m:int`: scalar, number of lags to use

# Output
- `S::Array`: qxq covariance matrix

"""
function NWFn(g0,m=0)

  T = size(g0,1)                    #g0 is Txq
  m = min(m,T-1)                    #number of lags

  g = g0 .- mean(g0,dims=1)         #normalizing to Eg=0

  S = g'g/T                         #(qxT)*(Txq)
  for s = 1:m
    Γ_s = g[s+1:T,:]'g[1:T-s,:]/T   #same as Sum[g(t)*g(t-s)',t=s+1,T]
    S   = S  +  (1 - s/(m+1))*(Γ_s + Γ_s')
  end

  return S

end

max_issues_repo_name: PaulSoderlind/EmpiricalFinancePhD

max_forks_repo_name: PaulSoderlind/EmpiricalFinancePhD


text:
module AddInit
using JSON

export @add_init

"""
    macro add_init(type::Expr)

    Automatically add a constructor for building objects with dict and json to DataType

# Examples:
```julia
    @add_init struct Test
        field:AbstractString
    end

    Test("{"field":"a"}") == Test("a")
    Test(Dict("field"=>"a")) == Test("a")
```
"""
macro add_init(expr::Expr)
    struct_expr = expr.args[2]
    local typename = esc(
        # 可能含有继承
        typeof(struct_expr) == Symbol ? struct_expr : struct_expr.args[1],
    ) # 消除卫生宏
    quote
        $(esc(expr))

        # 添加外部构造函数，利用Dict生成对象
        @add_dictinit $typename
        @add_jsoninit $typename
    end
end

"""
    macro add_init(symbol::Symbol)

    Automatically add a construction method for building objects with dict and json
    to the DataType, and adapt to the DataType modified by base.@kwdef

# Examples:
```julia
    @Base.kwdef struct Test
        a::Int
        b::Int=2
        c::Int
    end
    @add_init Test

    Test(Dict("a"=>1, "c"=>3)) == Test(1,2,3)
```
"""
macro add_init(symbol::Symbol)
    local typename = esc(symbol)
    quote
        if hasmethod($typename, Tuple{}, fieldnames($typename))
            # 如果有kwargs构造函数，例如使用了@kwdef
            $typename(dict::Dict) = $typename(;
                map([
                    k for k in keys(dict) if hasfield($typename, Symbol(k))
                ]) do key
                    constructor = fieldtype($typename, Symbol(key))
                    Symbol(key) =>
                        length(methods(constructor)) > 0 ?
                        constructor(dict[key]) : dict[key]
                end...
            )
        else
            # 添加外部构造函数，利用Dict生成对象
            esc(@add_dictinit $typename)
        end
        esc(@add_jsoninit $typename)
    end
end

macro add_jsoninit(typename)
    typename = esc(typename)
    :($typename(json::AbstractString) = $typename(JSON.parse(json)))
end

macro add_dictinit(typename)
    typename = esc(typename)
    quote
        $typename(dict::Dict{String,<:Any}) = $typename(
            map(fieldnames($typename)) do field
                constructor = fieldtype($typename, field)
                # 如果类型为Any或者其他抽象类型，直接返回传入值
                length(methods(constructor)) > 0 ?
                constructor(dict[String(field)]) : dict[String(field)]
            end...
        )
    end
end

end


text:
module DebugAdapter

import JSON
import JSONRPC
import JuliaInterpreter
import JSONRPC: @dict_readable, Outbound

include("packagedef.jl")

end


text:
"""
MODEL DESCRIPTION

Z-section cantilever under torsional loading.

Linear elastic analysis, Young's modulus = 210 GPa, Poisson's ratio = 0.3.

All displacements are fixed at X=0.

Torque of 1.2 MN-m applied at X=10. The torque is applied by two
uniformly distributed shear loads of 0.6 MN at each flange surface.

Objective of the analysis is to compute the axial stress at X = 2.5 from fixed end.

NAFEMS REFERENCE SOLUTION

Axial stress at X = 2.5 from fixed end (point A) at the midsurface is -108 MPa.
"""
module LE5_Z_cantilever_examples

using LinearAlgebra
using FinEtools
using FinEtoolsDeforLinear
using FinEtoolsFlexStructures.FESetShellT3Module: FESetShellT3
using FinEtoolsFlexStructures.FESetShellQ4Module: FESetShellQ4
using FinEtoolsFlexStructures.FEMMShellT3FFModule
using FinEtoolsFlexStructures.RotUtilModule: initial_Rfield, update_rotation_field!
using VisualStructures: plot_nodes, plot_midline, render, plot_space_box, plot_midsurface, space_aspectratio, save_to_json
using FinEtools.MeshExportModule.VTKWrite: vtkwrite

function zcant!(csmatout::FFltMat, XYZ::FFltMat, tangents::FFltMat, fe_label::FInt)
    r = vec(XYZ);
    cross3!(r, view(tangents, :, 1), view(tangents, :, 2))
    csmatout[:, 3] .= vec(r)/norm(vec(r))
    csmatout[:, 1] .= (1.0, 0.0, 0.0)
    cross3!(view(csmatout, :, 2), view(csmatout, :, 3), view(csmatout, :, 1))
    return csmatout
end

function _execute_model(formul, input = "nle5xf3c.inp", nrefs = 0, visualize = true)
    E = 210e9;
    nu = 0.3;
    L = 10.0;
    thickness = 0.1

    tolerance = thickness/1000
    output = import_ABAQUS(joinpath(dirname(@__FILE__()), input))
    fens = output["fens"]
    fes = output["fesets"][1]

    connected = findunconnnodes(fens, fes);
    fens, new_numbering = compactnodes(fens, connected);
    fes = renumberconn!(fes, new_numbering);

    for r in 1:nrefs
        fens, fes = T3refine(fens, fes)
    end

    mater = MatDeforElastIso(DeforModelRed3D, E, nu)
    ocsys = CSys(3, 3, zcant!)

    # Report
    @info "Mesh: $input, nrefs = $nrefs"

    sfes = FESetShellT3()
    accepttodelegate(fes, sfes)
    femm = formul.make(IntegDomain(fes, TriRule(1), thickness), mater)
    stiffness = formul.stiffness
    associategeometry! = formul.associategeometry!

    # Construct the requisite fields, geometry and displacement
    # Initialize configuration variables
    geom0 = NodalField(fens.xyz)
    u0 = NodalField(zeros(size(fens.xyz,1), 3))
    Rfield0 = initial_Rfield(fens)
    dchi = NodalField(zeros(size(fens.xyz,1), 6))

    # plots = cat(plot_space_box([[0 0 -L/2]; [L/2 L/2 L/2]]),
    #     plot_nodes(fens),
    #     plot_midsurface(fens, fes);
    # dims = 1)
    # pl = render(plots)

    # Apply EBC's
    # plane of symmetry perpendicular to X
    l1 = selectnode(fens; box = Float64[0 0 -Inf Inf -Inf Inf], inflate = tolerance)
    for i in [1,2,3,]
        setebc!(dchi, l1, true, i)
    end

    applyebc!(dchi)
    numberdofs!(dchi);
    @show dchi.nfreedofs

    # Assemble the system matrix
    associategeometry!(femm, geom0)
    K = stiffness(femm, geom0, u0, Rfield0, dchi);

    # Load
    nl = selectnode(fens; box = Float64[10.0 10.0 1.0 1.0 0 0], tolerance = tolerance)
    loadbdry1 = FESetP1(reshape(nl, 1, 1))
    lfemm1 = FEMMBase(IntegDomain(loadbdry1, PointRule()))
    fi1 = ForceIntensity(FFlt[0, 0, +0.6e6, 0, 0, 0]);
    nl = selectnode(fens; box = Float64[10.0 10.0 -1.0 -1.0 0 0], tolerance = tolerance)
    loadbdry2 = FESetP1(reshape(nl, 1, 1))
    lfemm2 = FEMMBase(IntegDomain(loadbdry2, PointRule()))
    fi2 = ForceIntensity(FFlt[0, 0, -0.6e6, 0, 0, 0]);
    F = distribloads(lfemm1, geom0, dchi, fi1, 3) + distribloads(lfemm2, geom0, dchi, fi2, 3);

    # @infiltrate
    # Solve
    U = K\F
    scattersysvec!(dchi, U[:])
    targetu =  minimum(dchi.values[:, 3]), maximum(dchi.values[:, 3])
    @info "Target: $(round.(targetu, digits=8))"

    # Generate a graphical display of displacements and rotations
    scalars = []
    for nc in 1:6
        push!(scalars, ("dchi$nc", deepcopy(dchi.values[:, nc])))
    end
    vectors = []
    push!(vectors, ("U", deepcopy(dchi.values[:, 1:3])))
    push!(vectors, ("UR", deepcopy(dchi.values[:, 4:6])))
    vtkwrite("z_cant-$input-$nrefs-dchi.vtu", fens, fes; scalars = scalars, vectors = vectors)

    # Generate a graphical display of resultants
    nl = selectnode(fens; nearestto = Float64[2.5 1.0 1.0])

    scalars = []
    for nc in 1:3
        fld = fieldfromintegpoints(femm, geom0, dchi, :moment, nc, outputcsys = ocsys)
            # fld = elemfieldfromintegpoints(femm, geom0, dchi, :moment, nc)
        push!(scalars, ("m$nc", fld.values))
    end
    vtkwrite("z_cant-$input-$nrefs-m.vtu", fens, fes; scalars = scalars, vectors = [("u", dchi.values[:, 1:3])])
    pointAstresses = Float64[]
    scalars = []
    for nc in 1:3
        fld = fieldfromintegpoints(femm, geom0, dchi, :membrane, nc, outputcsys = ocsys)
            # fld = elemfieldfromintegpoints(femm, geom0, dchi, :moment, nc)
        push!(pointAstresses, fld.values[nl[1], 1]/thickness)
        push!(scalars, ("n$nc", fld.values))
    end
    vtkwrite("z_cant-$input-$nrefs-n.vtu", fens, fes; scalars = scalars, vectors = [("u", dchi.values[:, 1:3])])
    scalars = []
    for nc in 1:2
        fld = fieldfromintegpoints(femm, geom0, dchi, :shear, nc, outputcsys = ocsys)
            # fld = elemfieldfromintegpoints(femm, geom0, dchi, :moment, nc)
        push!(scalars, ("q$nc", fld.values))
    end
    vtkwrite("z_cant-$input-$nrefs-q.vtu", fens, fes; scalars = scalars, vectors = [("u", dchi.values[:, 1:3])])

    # Visualization
    if !visualize
        return true
    end
    scattersysvec!(dchi, (L/8)/maximum(abs.(U)).*U)
    update_rotation_field!(Rfield0, dchi)
    plots = cat(plot_space_box([[0 0 -L/2]; [L/2 L/2 L/2]]),
        #plot_nodes(fens),
        plot_midsurface(fens, fes; x = geom0.values, u = dchi.values[:, 1:3], R = Rfield0.values);
    dims = 1)
    pl = render(plots)
    return pointAstresses
end

function test_convergence()
    formul = FEMMShellT3FFModule
    @info "LE5 Z-cantilever, formulation=$(formul)"
    for n in [0, 1, 2, 3, 4, ]
        _execute_model(formul, "nle5xf3c.inp", n, false)
    end
    return true
end

function allrun()
    println("#####################################################")
    println("# test_convergence ")
    test_convergence()
    return true
end # function allrun

@info "All examples may be executed with "
println("using .$(@__MODULE__); $(@__MODULE__).allrun()")

end # module
nothing


text:
"""
    AbstractDiscount

An abstract type for discount functions in GVFs.
"""
abstract type AbstractDiscount end

function Base.get(γ::AbstractDiscount, state_t, action_t, state_tp1, action_tp1, preds_tp1)
    throw(DomainError("get(DiscountType, args...) not defined!"))
end

"""
    ConstantDiscount{T} <: AbstractDiscount

Always returns a constant value γ
"""
struct ConstantDiscount{T} <: AbstractDiscount
    γ::T
end

# Base.get(cd::ConstantDiscount, state_t, action_t, state_tp1, action_tp1, preds_tp1) = cd.γ
Base.get(cd::ConstantDiscount; kwargs...) = cd.γ

"""
    StateTerminationDiscount{T<:Number, F} <: AbstractDiscount

Returns 0 if the condition is true and γ otherwise.
"""
struct StateTerminationDiscount{T<:Number, F} <: AbstractDiscount
    γ::T
    condition::F
    terminal::T

end

StateTerminationDiscount(γ, condition) =
    StateTerminationDiscount(γ, condition, zero(γ))

# Base.get(td::StateTerminationDiscount, state_t, action_t, state_tp1, action_tp1, preds_tp1) =
#     td.condition(state_tp1) ? td.terminal : td.γ

Base.get(td::StateTerminationDiscount; kwargs...) =
    td.condition(;kwargs...) ? td.terminal : td.γ

max_stars_repo_name: MatthewMcLeod/curiosity

max_issues_repo_name: MatthewMcLeod/curiosity

max_forks_repo_name: MatthewMcLeod/curiosity


text:
struct TextDocument
    _uri::URI
    _content::String
    _version::Int
    _line_offsets::Union{Nothing,Vector{Int}} # TODO Legacy, remove eventually
    _line_indices::Union{Nothing,Vector{Int}}

    function TextDocument(uri::URI, text::AbstractString, version::Int)
        # TODO Remove this check eventually
        occursin('\0', text) && throw(LSInvalidFile("Tried to set a text with an embedded NULL as the document content."))

        line_offsets = _compute_line_offsets(text)
        line_indices = _compute_line_indices(text)

        return new(uri, text, version, line_offsets, line_indices)
    end
end

Base.display(doc::TextDocument) = println("TextDocument: $(basename(doc._uri)) ")

get_text(doc::TextDocument) = doc._content

get_uri(doc::TextDocument) = doc._uri

get_version(doc::TextDocument) = doc._version

get_line_indices(doc::TextDocument) = doc._line_indices

"""
    index_at(doc::TextDocument, p::Position, forgiving_mode=false)

Converts a 0-based `Position` that is UTF-16 encoded to a 1-based UTF-8
encoded Julia string index.
"""
function index_at(doc::TextDocument, p::Position, forgiving_mode=false)
    line = p.line
    character = p.character

    line_indices = get_line_indices(doc)
    text = get_text(doc)

    if line >= length(line_indices)
        forgiving_mode || throw(LSOffsetError("index_at crashed. More diagnostics:\nline=$line\nline_indices='$line_indices'"))

        return nextind(text, lastindex(text))
    elseif line < 0
        throw(LSOffsetError("index_at crashed. More diagnostics:\nline=$line\nline_indices='$line_indices'"))
    end

    line_index = line_indices[line + 1]

    next_line_index = line + 1 < length(line_indices) ? line_indices[line + 2] : nextind(text, lastindex(text))

    pos = line_index

    while character > 0
        if pos >= next_line_index
            pos = next_line_index
            break
        end

        if UInt32(text[pos]) >= 0x010000
            character -= 2
        else
            character -= 1
        end

        pos = nextind(text, pos)
    end

    return pos
end

function apply_text_edits(doc::TextDocument, edits, new_version)
    content = doc._content

    for edit in edits
        if ismissing(edit.range) && ismissing(edit.rangeLength)
            # No range given, replace all text
            content = edit.text
        else
            editrange = _convert_lsrange_to_jlrange(doc, edit.range)
            content = string(content[1:prevind(content, editrange.start)], edit.text, content[nextind(content, editrange.stop):lastindex(content)])
        end
    end

    return TextDocument(doc._uri, content, new_version)
end

# =====================
# Internal functions
# =====================

function _convert_lsrange_to_jlrange(doc::TextDocument, range::Range)
    start_index_ls = index_at(doc, range.start)
    stop_index = index_at(doc, range.stop)

    content = doc._content

    # we use prevind for the stop value here because Julia stop values in
    # a range are inclusive, while the stop value is exclusive in a LS
    # range
    return start_index_ls:prevind(content, stop_index)
end

function _compute_line_indices(text)
    line_indices = Int[1]

    ind = firstindex(text)
    while ind <= lastindex(text)
        c = text[ind]
        if c == '\n' || c == '\r'
            if c == '\r' && ind + 1 <= lastindex(text) && text[ind + 1] == '\n'
                ind += 1
            end
            push!(line_indices, ind + 1)
        end

        ind = nextind(text, ind)
    end
    return line_indices
end

# Note: to be removed
function _obscure_text(s)
    i = 1
    io = IOBuffer()
    while i <= sizeof(s) # AUDIT: OK, i is generated by nextind
        di = nextind(s, i) - i
        if di == 1
            if s[i] in ('\n', '\r')
                write(io, s[i])
            else
                write(io, "a")
            end
        elseif di == 2
            write(io, "α")
        elseif di == 3
            write(io, "—")
        else
            write(io, s[i])
        end
        i += di
    end
    String(take!(io))
end

function _get_line_of(doc::TextDocument, offset::Integer)
    line_offsets = doc._line_offsets
    nlines = length(line_offsets)
    if offset > last(line_offsets)
        line = nlines
    else
        line = 1
        while line < nlines
            if line_offsets[line] <= offset < line_offsets[line + 1]
                break
            end
            line += 1
        end
end
    return line, line_offsets[line]
end

# ======================
# TODO The following functions are legacy offset versions that should eventually be removed
# ======================

function _compute_line_offsets(text)
    line_offsets = Int[0]

    ind = firstindex(text)
    while ind <= lastindex(text)
        c = text[ind]
        nl = c == '\n' || c == '\r'
        if c == '\r' && ind + 1 <= lastindex(text) && text[ind + 1] == '\n'
            ind += 1
        end
        nl && push!(line_offsets, ind)
        ind = nextind(text, ind)
    end

    return line_offsets
end

"""
    get_line_offsets(doc::Document)

Updates the doc._line_offsets field, an n length Array each entry of which
gives the byte offset position of the start of each line. This always starts
with 0 for the first line (even if empty).
"""
get_line_offsets(doc::TextDocument) = doc._line_offsets

"""
    get_offset(doc, line, char)

Returns the 0 based byte offset position corresponding to a line/character position.
This takes 0 based line/char inputs. Corresponding functions are available for
Position and Range arguments, the latter returning a UnitRange{Int}.
"""
function get_offset(doc::TextDocument, line::Integer, character::Integer)
    c = ' '
    line_offsets = get_line_offsets(doc)
    io = IOBuffer(get_text(doc))
    try
        seek(io, line_offsets[line + 1])
        while character > 0
            c = read(io, Char)
            character -= 1
            if UInt32(c) >= 0x010000
                character -= 1
            end
        end
        if UInt32(c) < 0x0080
            return position(io)
        elseif UInt32(c) < 0x0800
            return position(io) - 1
        elseif UInt32(c) < 0x010000
            return position(io) - 2
        else
            return position(io) - 3
        end
    catch err
        throw(LSOffsetError("get_offset crashed. More diagnostics:\nline=$line\ncharacter=$character\nposition(io)=$(position(io))\nline_offsets='$line_offsets'\ntext='$(_obscure_text(get_text(doc)))'\n\noriginal_error=$(sprint(Base.display_error, err, catch_backtrace()))"))
    end
end
get_offset(doc::TextDocument, p::Position) = get_offset(doc, p.line, p.character)
get_offset(doc::TextDocument, r::Range) = get_offset(doc, r.start):get_offset(doc, r.stop)

"""
    get_position_from_offset(doc, offset)

Returns the 0-based line and character position within a document of a given
byte offset.
"""
function get_position_from_offset(doc::TextDocument, offset::Integer)
    offset > sizeof(get_text(doc)) && throw(LSPositionToOffsetException("offset[$offset] > sizeof(content)[$(sizeof(get_text(doc)))]")) # OK, offset comes from EXPR spans
    line_offsets = get_line_offsets(doc)
    line, _ = _get_line_of(doc, offset)
    io = IOBuffer(get_text(doc))
    seek(io, line_offsets[line])
    character = 0
    while offset > position(io)
        c = read(io, Char)
        character += 1
        if UInt32(c) >= 0x010000
    character += 1
        end
    end
    close(io)
    return line - 1, character
end

"""
    Range(Doc, rng)
Converts a byte offset range to a LSP Range.
"""
function Range(doc::TextDocument, rng::UnitRange)
    start_l, start_c = get_position_from_offset(doc, first(rng))
    end_l, end_c = get_position_from_offset(doc, last(rng))
    rng = Range(start_l, start_c, end_l, end_c)
end


text:
using PyPlot,
  JLD2,
  Statistics,
  Printf,
  Random,
  FourierFlows

using FFTW: irfft
using Random: seed!
using Statistics: mean

import GeophysicalFlows.BarotropicQG
import GeophysicalFlows.BarotropicQG: energy, enstrophy

# Numerical parameters and time-stepping parameters
nx = 256       # 2D resolution = nx^2
stepper = "FilteredRK4"   # timestepper
dt  = 0.01     # timestep
nsteps = 40000 # total number of time-steps
nsubs  = 500   # number of time-steps for plotting
               # (nsteps must be multiple of nsubs)

# Physical parameters
 Lx  = 2π      # domain size
 nu  = 0e-05   # viscosity
 nnu = 1       # viscosity order
beta = 10.0    # planetary PV gradient
mu   = 0.01    # bottom drag

# Forcing
kf, dkf = 14.0, 1.5     # forcing wavenumber and width of
                        # forcing ring in wavenumber space
ε = 0.001               # energy input rate by the forcing

gr  = TwoDGrid(nx, Lx)

x, y = gridpoints(gr)
Kr = [ gr.kr[i] for i=1:gr.nkr, j=1:gr.nl]

force2k = @. exp(-(sqrt(gr.Krsq)-kf)^2/(2*dkf^2))
@. force2k[gr.Krsq < 2.0^2 ] .= 0
@. force2k[gr.Krsq > 20.0^2 ] .= 0
force2k[Kr .< 2π/Lx] .= 0
ε0 = FourierFlows.parsevalsum(force2k.*gr.invKrsq/2.0, gr)/(gr.Lx*gr.Ly)
force2k .= ε/ε0 * force2k  # normalization so that forcing injects
                           # energy ε per domain area per unit time

# reset of the random number generator for reproducibility
seed!(1234)

# the function that updates the forcing realization
function calcFq!(Fh, sol, t, cl, v, p, g)
  ξ = exp.(2π*im*rand(Float64, size(sol)))/sqrt(cl.dt)
  ξ[1, 1] = 0
  @. Fh = ξ*sqrt(force2k)
  Fh[abs.(Kr).==0] .= 0
  nothing
end

# Initialize problem
prob = BarotropicQG.ForcedProblem(nx=nx, Lx=Lx, beta=beta, nu=nu, nnu=nnu,
                                  mu=mu, dt=dt, stepper=stepper, calcFq=calcFq!, stochastic=true)
sol, cl, v, p, g = prob.sol, prob.clock, prob.vars, prob.params, prob.grid

# Files
filepath = "."
plotpath = "./plots_forcedbetaturb"
plotname = "snapshots"
filename = joinpath(filepath, "forcedbetaturb.jld2")

# File management
if isfile(filename); rm(filename); end
if !isdir(plotpath); mkdir(plotpath); end

# Zero initial condition
BarotropicQG.set_zeta!(prob, 0*x)

# Create Diagnostic -- "energy" and "enstrophy" are functions imported at the top.
E = Diagnostic(energy, prob; nsteps=nsteps)
Z = Diagnostic(enstrophy, prob; nsteps=nsteps)
diags = [E, Z] # A list of Diagnostics types passed to "stepforward!" will
# be updated every timestep.

# Create Output
get_sol(prob) = sol # extracts the Fourier-transformed solution
get_u(prob) = irfft(im*g.l.*g.invKrsq.*sol, g.nx)
out = Output(prob, filename, (:sol, get_sol), (:u, get_u))

function plot_output(prob, fig, axs; drawcolorbar=false)
  # Plot the vorticity and streamfunction fields as well as the zonal mean
  # vorticity and the zonal mean zonal velocity.

  sol, v, p, g = prob.sol, prob.vars, prob.params, prob.grid
  BarotropicQG.updatevars!(prob)

  sca(axs[1])
  cla()
  pcolormesh(x, y, v.q)
  axis("square")
  xticks(-2:2:2)
  yticks(-2:2:2)
  title(L"vorticity $\zeta = \partial_x v - \partial_y u$")
  if drawcolorbar==true
    colorbar()
  end

  sca(axs[2])
  cla()
  contourf(x, y, v.psi)
  if maximum(abs.(v.psi))>0
    contour(x, y, v.psi, colors="k")
  end
  axis("square")
  xticks(-2:2:2)
  yticks(-2:2:2)
  title(L"streamfunction $\psi$")
  if drawcolorbar==true
    colorbar()
  end

  sca(axs[3])
  cla()
  plot(Array(transpose(mean(v.zeta, dims=1))), y[1,:])
  plot(0*y[1,:], y[1,:], "k--")
  ylim(-Lx/2, Lx/2)
  xlim(-4, 4)
  title(L"zonal mean $\zeta$")

  sca(axs[4])
  cla()
  plot(Array(transpose(mean(v.u, dims=1))), y[1,:])
  plot(0*y[1,:], y[1,:], "k--")
  ylim(-Lx/2, Lx/2)
  xlim(-0.7, 0.7)
  title(L"zonal mean $u$")

  sca(axs[5])
  cla()
  plot(mu*E.t[1:E.i], E.data[1:E.i], label="energy")
  xlabel(L"\mu t")
  legend()

  sca(axs[6])
  cla()
  plot(mu*Z.t[1:Z.i], Z.data[1:E.i], label="enstrophy")
  xlabel(L"\mu t")
  legend()

  pause(0.001)
end

fig, axs = subplots(ncols=3, nrows=2, figsize=(14, 8))
plot_output(prob, fig, axs; drawcolorbar=false)

# Step forward
startwalltime = time()

while cl.step < nsteps
  stepforward!(prob, diags, nsubs)

  BarotropicQG.updatevars!(prob)

  # Message
  cfl = cl.dt*maximum([maximum(v.u)/g.dx, maximum(v.v)/g.dy])
  log = @sprintf("step: %04d, t: %d, cfl: %.2f, E: %.4f, Q: %.4f, τ: %.2f min",
    cl.step, cl.t, cfl, E.data[E.i], Z.data[Z.i],
    (time()-startwalltime)/60)

  println(log)

  plot_output(prob, fig, axs; drawcolorbar=false)
end

# how long did it take?
println((time()-startwalltime))

plot_output(prob, fig, axs; drawcolorbar=false)

# save the figure as png
savename = @sprintf("%s_%09d.png", joinpath(plotpath, plotname), cl.step)
savefig(savename)


text:
using Test, GuessworkQuantumSideInfo, LinearAlgebra, Statistics, Random, UnPack, StableRNGs
using GenericLinearAlgebra
using SCS

TEST_MATLAB = false # requires MATLAB and MATLAB.jl installed
TEST_MISDP = false # requires Pajarito or another MISDP solver
TEST_BB84_MISDP = false # takes ~100 seconds; requires TEST_MISDP
TEST_MOI = true # incompatible with the current version of Pajarito

default_sdp_solver() = TEST_MOI ? SCS.Optimizer(verbose = 0, eps = 1e-6) : SCSSolver(verbose = 0, eps = 1e-6)

if TEST_MATLAB
    include("test_matlab.jl")
end

include("test_problems.jl")


text:
""" Categories of graphs and other categorical and algebraic aspects of graphs.
"""
module GraphCategories

using ..FinSets, ..CSets, ..Limits
using ...Graphs.BasicGraphs
import ...Graphs.GraphAlgorithms: connected_component_projection

function connected_component_projection(g::AbstractACSet)::FinFunction
  proj(coequalizer(FinFunction(src(g), nv(g)),
                   FinFunction(tgt(g), nv(g))))
end

end


text:
module LIPO

using LinearAlgebra

function lipo(f; k=0.0, n = 100, lower=0.0, upper=2.0, max_ndraws=1e4)
    Xdom = [0.0, upper]

    X = [rand()*(upper-lower)+lower]
    fX = [f(X[1])]

    U(Xtp1, X, fX, k) = minimum(fXt+k*norm(Xtp1 - Xt, 2) for (Xt, fXt) in zip(X, fX))
    t = 1
    ndraws = 0
    while t < n
       Xtp1 = rand()*(upper-lower)+lower
       ndraws += 1
       if U(Xtp1, X, fX, k) >= maximum(fX)
           X = push!(X, Xtp1)
           fX = push!(fX, f(X[end]))
           t += 1
       end
       if ndraws >= max_ndraws
           break
       end
    end
    return X, fX, ndraws
end

function adalipo(f; k_α=0.5, n = 100, lower=0.0, upper=2.0, max_ndraws=1e4, choose_max=true)
    Xdom = [lower, upper]

    # first, check that function isn't constant, run for max
    # safeguarding number of iterations (used to ensure a proper)

    X = [rand()*(upper-lower)+lower]
    fX = [f(X[1])]

    U(Xtp1, X, fX, k) = minimum(fXt+k*norm(Xtp1-Xt, 2) for (Xt, fXt) in zip(X, fX))
    t = 1
    ndraws = 0
    nexploit = 0
    nexplore = 0
    p = 0.4
    khat_i = 0.0
    while t < n
        if rand() < p
            Xtp1 = rand()*(upper-lower)+lower
            nexplore += 1
        else
            Xtp1 = exploit(fX, X, khat_i, choose_max, lower, upper)
            nexploit += 1
        end

       # a draw can come from
       ndraws += 1
       if U(Xtp1, X, fX, khat_i) >= maximum(fX)
           X = push!(X, Xtp1)
           fX = push!(fX, f(X[end]))
           t += 1
       end

       khat_i = khat(k_α, fX, X)

       if ndraws >= max_ndraws
           break
       end

    end
    return X, fX, ndraws
end

function exploit(fX, X, khat_i, choose_max, lower, upper)
    # pick either the largest or from the set of
    # tops
    val_at_tops, x_at_tops = find_tops(fX, X, khat_i, l)
    if chose_max
        i_m = argmax(val_at_tops)
    else
        i_m = rand(val_at_tops)
    end
    Xtp1 = x_at_tops[i_m]
end

function find_tops(fX, X, khat_i, lower, upper)
    xperm = sortperm(X)
    x = X[xperm]
    f = fX[perm]
    if f[1]+norm(x) < f[end]
        x_top = f[end]
        f_top = f[end]
    else
        x_top = x[1]
        f_top = f[1]
    end

    x_top = 0.0
    val_top = NaN
    for 1:length(x)-1
        if norm(x[i]-x[i+1])
    end
end

function khat(k_α, fX::AbstractVector{T}, X; max_depth = 100) where T
    #find estimates

    argmax_val = -T(1)
    for j = 1:length(X)
        for i = 1:length(X)
            if i == j
                continue
            end
            num = abs(fX[i] - fX[j])
            den = norm(X[i]- X[j], 2)
            candidate = num/den
            if candidate > argmax_val
                argmax_val = candidate
            end
        end
    end

    # Using the cloesd form solution in the paper
    k_i = ceil(log(argmax_val)/log(1+k_α))
    k_i
end

end # module


text:
type VcovCluster  <: AbstractVcovMethod
    clusters::Vector{Symbol}
end
VcovCluster(x::Symbol) = VcovCluster([x])
allvars(x::VcovCluster) = x.clusters

type VcovClusterData <: AbstractVcovMethodData
    clusters::DataFrame
    size::Dict{Symbol, Int}
end

function VcovMethodData(v::VcovCluster, df::AbstractDataFrame)
    vclusters = DataFrame(Vector, size(df, 1), length(v.clusters))
    names!(vclusters, v.clusters)
    vsize = Dict{Symbol, Int}()
    for c in v.clusters
        p = df[c]
        typeof(p) <: PooledDataVector || error("Cluster variable $(c) is of type $(typeof(p)), but should be a PooledDataVector.")
        vclusters[c] = p
        # may be subset / NA
        vsize[c] = length(unique(p.refs))
    end
    return VcovClusterData(vclusters, vsize)
end

df_FStat(v::VcovClusterData, ::VcovData, ::Bool) = minimum(values(v.size)) - 1

function vcov!(v::VcovClusterData, x::VcovData)
    S = shat!(v, x)
    out = sandwich(x.crossmatrix, S)
    # Cameron, Gelbach, & Miller (2011)
    pinvertible(out)
    return out
end
function shat!{T}(v::VcovClusterData, x::VcovData{T, 1})
    # Cameron, Gelbach, & Miller (2011).
    clusternames = names(v.clusters)
    X = x.regressors
    broadcast!(*, X, X, x.residuals)
    S = fill(zero(Float64), (size(X, 2), size(X, 2)))
    for i in 1:length(clusternames)
        for c in combinations(clusternames, i)
            if length(c) == 1
                f = (v.clusters)[c[1]]
                # no need to group in this case
                fsize = (v.size)[c[1]]
            else
                df = v.clusters[c]
                f = group(df)
                fsize = length(f.pool)
            end
            if rem(length(c), 2) == 1
                S += helper_cluster(X, f, fsize)
            else
                S -= helper_cluster(X, f, fsize)
            end
        end
    end
    scale!(S, (size(X, 1) - 1) / x.df_residual)
    return S
end

function helper_cluster(Xu::Matrix{Float64}, f::PooledDataVector, fsize::Int)
    if fsize == size(Xu, 1)
        # if only one obs by pool, use White, as in Petersen (2009) & Thomson (2011)
        return At_mul_B(Xu, Xu)
    else
        # otherwise
        X2 = fill(zero(Float64), (fsize, size(Xu, 2)))
        for j in 1:size(Xu, 2)
             @inbounds @simd for i in 1:size(Xu, 1)
                X2[f.refs[i], j] += Xu[i, j]
            end
        end
        out = At_mul_B(X2, X2)
        scale!(out, fsize / (fsize- 1))
        return out
    end
end

function shat!{T}(v::VcovClusterData, x::VcovData{T, 2})
    # Cameron, Gelbach, & Miller (2011).
    clusternames = names(v.clusters)
    X = x.regressors
    res = x.residuals
    dim = (size(X, 2) *size(res, 2))
    S = fill(zero(Float64), (dim, dim))
    for i in 1:length(clusternames)
        for c in combinations(clusternames, i)
            if length(c) == 1
                f = (v.clusters)[c[1]]
                # no need to group in this case
                fsize = (v.size)[c[1]]
            else
                df = v.clusters[c]
                f = group(df)
                fsize = length(f.pool)
            end
            if rem(length(c), 2) == 1
                S += helper_cluster(X, res, f, fsize)
            else
                S -= helper_cluster(X, res, f, fsize)
            end
        end
    end
    scale!(S, (size(X, 1) - 1) / x.df_residual)
    return S
end

# S_{(l-1) * K + k, (l'-1)*K + k'} = \sum_g (\sum_{i in g} X[i, k] res[i, l]) (\sum_{i in g} X[i, k'] res[i, l'])
function helper_cluster(X::Matrix{Float64}, res::Matrix{Float64}, f::PooledDataVector, fsize::Int)
    dim = size(X, 2) * size(res, 2)
    nobs = size(X, 1)
    S = fill(zero(Float64), (dim, dim))
    temp = fill(zero(Float64), fsize, dim)
    if fsize == nobs
        index = zero(Int)
        @inbounds for l in 1:size(res, 2), k in 1:size(X, 2)
            index += 1
            @simd for i in 1:nobs
                temp[i, index] = X[i, k] * res[i, l]
            end
        end
        S = At_mul_B(temp, temp)
    else
        index = zero(Int)
        @inbounds for l in 1:size(res, 2), k in 1:size(X, 2)
            index += 1
            @simd for i in 1:nobs
                temp[f.refs[i], index] += X[i, k] * res[i, l]
            end
        end
        S = At_mul_B(temp, temp)
        scale!(S, fsize / (fsize - 1))
    end
    return S
end

function pinvertible(A::Matrix, tol = eps(real(float(one(eltype(A))))))
    SVD         = svdfact(A, thin=true)
    Stype       = eltype(SVD.S)
    small = SVD.S .<= tol
    if any(small)
        warn("estimated covariance matrix of moment conditions not of full rank.
                 model tests should be interpreted with caution.")
        SVD.S[small] = 0
        return  SVD.U * diagm(SVD.S) * SVD.Vt
    else
        return A
    end
end


text:
isreal(::Type{<:Real}) = true
isreal(::Type{<:Complex}) = false
isreal(::Type{T}) where {T} = isreal(eltype(T))

const StaticTypes = Union{Number,<:StaticVector{N} where N,<:NTuple{N,Any} where N}

"What is the euclidean dimension of the given type (if applicable)?"
euclideandimension(::Type{T}) where {T <: Number} = 1
euclideandimension(::Type{T}) where {N,T <: StaticVector{N}} = N
euclideandimension(::Type{T}) where {N,T <: NTuple{N,Any}} = N
# Does not apply to Vector{T}: we don't know its dimension

unitvector(d::Domain{T}) where {N,S,T<:SVector{N,S}} = SVector{N,S}(ntuple(i -> i==1, N))
function unitvector(d::Domain{T}) where {T<:AbstractVector}
    p = zeros(eltype(T), dimension(d))
    p[1] = 1
    p
end
unitvector(d::Domain{T}) where {T<:Number} = one(T)

#################
# Precision type
#################

"The floating point precision type associated with the argument."
prectype(x) = prectype(typeof(x))
prectype(::Type{<:Complex{T}}) where {T} = prectype(T)
prectype(::Type{<:AbstractArray{T}}) where {T} = prectype(T)
prectype(::Type{NTuple{N,T}}) where {N,T} = prectype(T)
prectype(::Type{Tuple{A}}) where {A} = prectype(A)
prectype(::Type{Tuple{A,B}}) where {A,B} = prectype(A,B)
prectype(::Type{Tuple{A,B,C}}) where {A,B,C} = prectype(A,B,C)
@generated function prectype(T::Type{<:Tuple{Vararg}})
    quote $(promote_type(map(prectype, T.parameters[1].parameters)...)) end
end
prectype(::Type{T}) where {T<:AbstractFloat} = T
prectype(::Type{T}) where {T<:Number} = prectype(float(T))

prectype(a...) = promote_type(map(prectype, a)...)

"Convert `x` such that its `prectype` equals `U`."
convert_prectype(x, ::Type{U}) where {U} = convert(to_prectype(typeof(x),U), x)

# function convert_prectype(x, ::Type{U}) where {U}
#     @warn "The order of arguments of convert_prectype has changed."
#     prectype(U, x)
# end

"Return the type to which `U` can be converted, such that the `prectype` becomes `T`."
to_prectype(::Type{T}, ::Type{U}) where {T,U} = error("Don't know how to convert the numtype of $(T) to $(U).")
to_prectype(::Type{T}, ::Type{U}) where {T <: Real,U <: Real} = U
to_prectype(::Type{Complex{T}}, ::Type{U}) where {T <: Real,U <: Real} = Complex{U}
to_prectype(::Type{SVector{N,T}}, ::Type{U}) where {N,T,U} = SVector{N,to_prectype(T,U)}
to_prectype(::Type{Vector{T}}, ::Type{U}) where {T,U} = Vector{to_prectype(T,U)}

"Promote the precision types of the arguments to a joined supertype."
promote_prectype(a) = a
promote_prectype(a, b) = _promote_prectype(prectype(a,b), a, b)
promote_prectype(a, b, c...) = _promote_prectype(prectype(a,b,c...), a, b, c...)
_promote_prectype(U, a) = convert_prectype(a, U)
_promote_prectype(U, a, b) = convert_prectype(a, U), convert_prectype(b, U)
_promote_prectype(U, a, b, c...) =
    (convert_prectype(a, U), convert_prectype(b, U), _promote_prectype(U, c...)...)

#################
# Numeric type
#################

"The numeric element type of x in a Euclidean space."
numtype(x) = numtype(typeof(x))
numtype(::Type{T}) where {T<:Number} = T
numtype(::Type{T}) where {T} = eltype(T)
numtype(::Type{NTuple{N,T}}) where {N,T} = T
numtype(::Type{Tuple{A,B}}) where {A,B} = promote_type(numtype(A), numtype(B))
numtype(::Type{Tuple{A,B,C}}) where {A,B,C} = promote_type(numtype(A), numtype(B), numtype(C))
numtype(::Type{Tuple{A,B,C,D}}) where {A,B,C,D} = promote_type(numtype(A), numtype(B), numtype(C), numtype(D))
@generated function numtype(T::Type{<:Tuple{Vararg}})
    quote $(promote_type(T.parameters[1].parameters...)) end
end

numtype(a...) = promote_type(map(numtype, a)...)

"Convert `x` such that its `numtype` equals `U`."
convert_numtype(x, ::Type{U}) where {U} = convert(to_numtype(typeof(x),U), x)

to_numtype(::Type{T}, ::Type{U}) where {T,U} = error("Don't know how to convert the numtype of $(T) to $(U).")
to_numtype(::Type{T}, ::Type{U}) where {T <: Number,U <: Number} = U
to_numtype(::Type{SVector{N,T}}, ::Type{U}) where {N,T,U} = SVector{N,to_numtype(T,U)}
to_numtype(::Type{Vector{T}}, ::Type{U}) where {T,U} = Vector{to_numtype(T,U)}

"Promote the numeric types of the arguments to a joined supertype."
promote_numtype(a) = a
promote_numtype(a, b) = _promote_numtype(numtype(a,b), a, b)
promote_numtype(a, b, c...) = _promote_numtype(numtype(a,b,c...), a, b, c...)
_promote_numtype(U, a) = convert_numtype(a, U)
_promote_numtype(U, a, b) = convert_numtype(a, U), convert_numtype(b, U)
_promote_numtype(U, a, b, c...) =
    (convert_numtype(a, U), convert_numtype(b, U), _promote_numtype(U, c...)...)

## Conversion from nested vectors to flat vectors and back

"""
Convert a vector from a cartesian format to a nested tuple according to the
given dimensions.

For example:
`convert_fromcartesian([1,2,3,4,5], Val{(2,2,1)}()) -> ([1,2],[3,4],5)`
"""
@generated function convert_fromcartesian(x::AbstractVector, ::Val{DIM}) where {DIM}
	dimsum = [0; cumsum([d for d in DIM])]
	E = Expr(:tuple, [ (dimsum[i+1]-dimsum[i] > 1 ? Expr(:call, :SVector, [:(x[$j]) for j = dimsum[i]+1:dimsum[i+1]]...) : :(x[$(dimsum[i+1])])) for i in 1:length(DIM)]...)
	return quote $(E) end
end

"The inverse function of `convert_fromcartesian`."
@generated function convert_tocartesian(x, ::Val{DIM}) where {DIM}
    dimsum = [0; cumsum([d for d in DIM])]
    E = vcat([[:(x[$i][$j]) for j in 1:DIM[i]] for i in 1:length(DIM)]...)
    quote SVector($(E...)) end
end


text:
function sphericalImpedance(
    f::Real,
    r::Real,
    ρ::Real,
    c::Real
    )

    kr = 2π * f * r / c

    return im * ρ * c * kr / (1.0 + im * kr)

end

function sphericalWave(
    U::Real,
    a::Real,
    f::Real,
    t::Real,
    r::Real,
    ρ::Real,
    c::Real
    )

    Z = sphericalImpedance(f, a, ρ, c)

    ω = 2π * f
    k = ω / c

    ϕ = mod2pi(ω * t - k * (r - a))

    return Z * a * U * exp(im * ϕ) / r

end

function decibel(x::T, ref::Real) where {T<:Number}
    return 20.0 * log10(abs(x) / ref)
end

max_stars_repo_name: CrocoDuckoDucks/Acoustics-Models

max_issues_repo_name: CrocoDuckoDucks/Acoustics-Models

max_forks_repo_name: CrocoDuckoDucks/Acoustics-Models


text:
eval(Expr(
    :struct, true, :PLEXOSSolutionDatasetSummary, Expr(:block,
        [:($(t.fieldname)::Tuple{Int,Int}) for t in plexostables]...
    )
))

PLEXOSSolutionDatasetSummary() =
    PLEXOSSolutionDatasetSummary(((0,0) for _ in 1:length(plexostables))...)

function PLEXOSSolutionDatasetSummary(
    zippath::String, xmlname::String=defaultxml(zippath)
)

    resultsarchive = _open_plexoszip(zippath)
    xml = parsexml(resultsarchive[xmlname])
    return PLEXOSSolutionDatasetSummary(xml)

end

function PLEXOSSolutionDatasetSummary(xml::Document)

    summary = PLEXOSSolutionDatasetSummary()

    for element in eachelement(xml.root)

        # Ignore the band table
        element.name == "t_band" && continue

        table = plexostables_lookup[element.name]
        count, maxidx = getfield(summary, table.fieldname)

        if isnothing(table.identifier)
            setfield!(summary, table.fieldname, (count + 1, count + 1))
        else
            idx = getchildint(table.identifier, element) + table.indexoffset
            setfield!(summary, table.fieldname, (count + 1, max(maxidx, idx)))
        end

    end

    return summary

end


text:
module MaxHelpingHandMadelinePonytailTrigger

using ..Ahorn, Maple

@mapdef Trigger "MaxHelpingHand/MadelinePonytailTrigger" MadelinePonytailTrigger(x::Integer, y::Integer, width::Integer=Maple.defaultTriggerWidth, height::Integer=Maple.defaultTriggerHeight, enable::Bool=true)

const placements = Ahorn.PlacementDict(
    "Madeline Ponytail Trigger (max480's Helping Hand)" => Ahorn.EntityPlacement(
        MadelinePonytailTrigger,
        "rectangle"
    )
)

end

max_stars_repo_name: Kitty-Cats/MaxHelpingHand

max_issues_repo_name: Kitty-Cats/MaxHelpingHand

max_forks_repo_name: Kitty-Cats/MaxHelpingHand


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct IoK8sApiFlowcontrolV1alpha1NonResourcePolicyRule <: SwaggerModel
    nonResourceURLs::Any # spec type: Union{ Nothing, Vector{String} } # spec name: nonResourceURLs
    verbs::Any # spec type: Union{ Nothing, Vector{String} } # spec name: verbs

    function IoK8sApiFlowcontrolV1alpha1NonResourcePolicyRule(;nonResourceURLs=nothing, verbs=nothing)
        o = new()
        validate_property(IoK8sApiFlowcontrolV1alpha1NonResourcePolicyRule, Symbol("nonResourceURLs"), nonResourceURLs)
        setfield!(o, Symbol("nonResourceURLs"), nonResourceURLs)
        validate_property(IoK8sApiFlowcontrolV1alpha1NonResourcePolicyRule, Symbol("verbs"), verbs)
        setfield!(o, Symbol("verbs"), verbs)
        o
    end
end # type IoK8sApiFlowcontrolV1alpha1NonResourcePolicyRule

const _property_map_IoK8sApiFlowcontrolV1alpha1NonResourcePolicyRule = Dict{Symbol,Symbol}(Symbol("nonResourceURLs")=>Symbol("nonResourceURLs"), Symbol("verbs")=>Symbol("verbs"))
const _property_types_IoK8sApiFlowcontrolV1alpha1NonResourcePolicyRule = Dict{Symbol,String}(Symbol("nonResourceURLs")=>"Vector{String}", Symbol("verbs")=>"Vector{String}")
Base.propertynames(::Type{ IoK8sApiFlowcontrolV1alpha1NonResourcePolicyRule }) = collect(keys(_property_map_IoK8sApiFlowcontrolV1alpha1NonResourcePolicyRule))
Swagger.property_type(::Type{ IoK8sApiFlowcontrolV1alpha1NonResourcePolicyRule }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_IoK8sApiFlowcontrolV1alpha1NonResourcePolicyRule[name]))}
Swagger.field_name(::Type{ IoK8sApiFlowcontrolV1alpha1NonResourcePolicyRule }, property_name::Symbol) =  _property_map_IoK8sApiFlowcontrolV1alpha1NonResourcePolicyRule[property_name]

function check_required(o::IoK8sApiFlowcontrolV1alpha1NonResourcePolicyRule)
    (getproperty(o, Symbol("nonResourceURLs")) === nothing) && (return false)
    (getproperty(o, Symbol("verbs")) === nothing) && (return false)
    true
end

function validate_property(::Type{ IoK8sApiFlowcontrolV1alpha1NonResourcePolicyRule }, name::Symbol, val)
end


text:
module SceneComposer

using StaticArrays
using CoordinateTransformations
using Rotations

export

    SceneNode,
    Vector3D,
    Quat

const Vector3D{T} = SVector{3,T}

mutable struct SceneNode{T}
    parent::Nullable{SceneNode{T}}
    children::Vector{SceneNode{T}}
    local_position::Vector3D{T}
    local_orientation::Quat{T}
    local_scale::Vector3D{T}
    global_position::Vector3D{T}
    global_orientation::Quat{T}
    global_scale::Vector3D{T}

    function SceneNode(position::Vector3D{T} = Vector3D(0.,0,0),
                       orientation::Quat{T} = Quat(1.,0,0,0),
                       scale::Vector3D{T} = Vector3D(1.,1,1)) where T
        new{T}(
            Nullable{SceneNode{T}}(),
            SceneNode{T}[],
            position,
            orientation,
            scale,
            position,
            orientation,
            scale,
        )
    end
end

function Base.show(io::IO, node::SceneNode)
    println(io, summary(node))
    println(io, " hierachy")
    println(io, "  - parent: ", !isnull(node.parent))
    println(io, "  - children: ", length(node.children))
    if !isnull(node.parent)
        println(io, " local")
        println(io, "  - position: ", node.local_position)
        println(io, "  - orientation: ", node.local_orientation)
        println(io, "  - scale: ", node.local_scale)
    end
    println(io, " global")
    println(io, "  - position: ", node.global_position)
    println(io, "  - orientation: ", node.global_orientation)
    print(io, "  - scale: ", node.global_scale)
end

function update_global!(node::SceneNode)
    if isnull(node.parent)
        # no parent means that local == global
        node.global_position = node.local_position
        node.global_orientation = node.local_orientation
        node.global_scale = node.local_scale
    else
        # cache global values using the local relative to parent
        parent = get(node.parent)
        node.global_scale = node.local_scale .* parent.global_scale
        node.global_orientation = parent.global_orientation * node.local_orientation
        node.global_position = parent.global_position .+ (parent.global_orientation * (node.local_position .* parent.global_scale))
    end
    foreach(update_global!, node.children)
    node
end

function detach_parent!(node::SceneNode{T}, update=true) where T
    if !isnull(node.parent)
        # remove from parents children list
        parent_children = get(node.parent).children
        for i in 1:length(parent_children)
            if parent_children[i] == node
                deleteat!(parent_children, i)
                break
            end
        end
        # set new parent as null
        node.parent = Nullable{SceneNode{T}}()
    end
    update && update_global!(child)
    node
end

function attach_parent!(child::SceneNode{T}, parent::SceneNode{T}) where T
    # detach child from its current parent (without re-caching)
    detach_parent!(child, false)
    # attach child to new parent
    push!(parent.children, child)
    child.parent = Nullable{SceneNode{T}}(parent)
    update_global!(child)
    child
end

# package code goes here

end # module


text:
struct Selector
    gid::UInt64
    tag::Symbol # :default, :invalid, :Gibbs, :HMC, etc.
    rerun::Bool
end
function Selector(tag::Symbol=:default, rerun=tag != :default)
    return Selector(time_ns(), tag, rerun)
end
function Selector(gid::Integer, tag::Symbol=:default)
    return Selector(gid, tag, tag != :default)
end
hash(s::Selector) = hash(s.gid)
==(s1::Selector, s2::Selector) = s1.gid == s2.gid

max_stars_repo_name: TuringLang/MicroPPL

max_issues_repo_name: TuringLang/MicroPPL

max_forks_repo_name: TuringLang/MicroPPL


text:
# This file is a part of AstroLib.jl. License is MIT "Expat".
# Copyright (C) 2016 Mosè Giordano.

# TODO: give sense to "-0" ("-0" is bitwise equal to "0").
ten(degrees::T, minutes::T, seconds::T) where {T<:AbstractFloat} =
    copysign(1, degrees)*(abs(degrees) + minutes/60 + seconds/3600)

ten(d::Real, m::Real=0, s::Real=0) =
    ten(promote(float(d), float(m), float(s))...)

# TODO: improve performance, if possible.  There are a couple of slow tests to
# make parsing of the string work.
function ten(strng::AbstractString)
    # Convert strings into numbers, empty strings into 0s.
    # Replace in the string multiple spaces or colons with a single space, strip leading
    # whitespaces, and split the resulting string using the space as separator.
    tmp = map(x-> x=="" ? 0.0 : parse(Float64, x),
              split(lstrip(replace(strng, r"(\:| )+" => s" ")), " "))
    # Concatenate "tmp" with 3 zeros, so that "angle" has at least 3 elements
    # also with an empty string in input.
    angle = vcat(tmp, zeros(Float64, 3))::Vector{Float64}
    ten(angle[1], angle[2], angle[3])
end

ten(itr) = ten(itr...)

"""
    ten(deg[, min, sec]) -> decimal
    ten("deg:min:sec") -> decimal

### Purpose ###

Converts a sexagesimal number or string to decimal.

### Explanation ###

`ten` is the inverse of the `sixty` function.

### Arguments ###

`ten` takes as argument either three scalars (`deg`, `min`, `sec`) or a string.
The string should have the form `"deg:min:sec"` or `"deg min sec"`.  Also any
iterable like `(deg, min, sec)` or `[deg, min, sec]` is accepted as argument.

If minutes and seconds are not specified they default to zero.

### Output ###

The decimal conversion of the sexagesimal numbers provided is returned.

### Method ###

The formula used for the conversion is

\$\$\\mathrm{sign}(\\mathrm{deg})·\\left(|\\mathrm{deg}| + \\frac{\\mathrm{min}}{60} + \\frac{\\mathrm{sec}}{3600}\\right)\$\$

### Example ###

```jldoctest
julia> using AstroLib

julia> ten(-0.0, 19, 47)
-0.3297222222222222

julia> ten("+5:14:58")
5.249444444444444

julia> ten("-10 26")
-10.433333333333334

julia> ten((-10, 26))
-10.433333333333334
```

### Notes ###

These functions cannot deal with `-0` (negative integer zero) in numeric input.
If it is important to give sense to negative zero, you can either make sure to
pass a floating point negative zero `-0.0` (this is the best option), or use
negative minutes and seconds, or non-integer negative degrees and minutes.
"""
ten