text:
"""
    SurfacePanel{TF}

Lifting surface panel with attached vortex ring

**Fields**
 - `rtl`: position of the left side of the top bound vortex
 - `rtc`: position of the center of the top bound vortex
 - `rtr`: position of the right side of the top bound vortex
 - `rbl`: position of the left side of the bottom bound vortex
 - `rbc`: position of the center of the bottom bound vortex
 - `rbr`: position of the right side of the bottom bound vortex
 - `rcp`: position of the panel control point
 - `ncp`: normal vector at the panel control point
 - `core_size`: finite core size (for use when the finite core smoothing model is enabled)
 - `chord`: panel chord length (for determining unsteady forces)
"""
struct SurfacePanel{TF}
    rtl::SVector{3, TF}
    rtc::SVector{3, TF}
    rtr::SVector{3, TF}
    rbl::SVector{3, TF}
    rbc::SVector{3, TF}
    rbr::SVector{3, TF}
    rcp::SVector{3, TF}
    ncp::SVector{3, TF}
    core_size::TF
    chord::TF
end

"""
    SurfacePanel(rtl, rtr, rbl, rbr, rcp, ncp, core_size, chord; kwargs...)

Construct and return a vortex ring panel.

# Arguments
 - `rtl`: position of the left side of the top bound vortex
 - `rtr`: position of the right side of the top bound vortex
 - `rbl`: position of the left side of the bottom bound vortex
 - `rbr`: position of the right side of the bottom bound vortex
 - `rcp`: position of the panel control point
 - `ncp`: normal vector at the panel control point
 - `core_size`: finite core size (for use when the finite core smoothing model is enabled)
 - `chord`: panel chord length (for determining unsteady forces)

# Keyword Arguments
 - `rtc`: position of the center of the top bound vortex, defaults to `(rtl+rtr)/2`
 - `rbc`: position of the center of the bottom bound vortex, defaults to `(rbl+rbr)/2`
"""
SurfacePanel(args...; kwargs...)

function SurfacePanel(rtl, rtr, rbl, rbr, rcp, ncp, core_size, chord;
    rtc = (rtl+rtr)/2, # default to average of corners
    rbc = (rbl+rbr)/2) # default to average of corners

    return SurfacePanel(rtl, rtc, rtr, rbl, rbc, rbr, rcp, ncp, core_size, chord)
end

function SurfacePanel(rtl, rtc, rtr, rbl, rbc, rbr, rcp, ncp, core_size, chord)

    TF = promote_type(eltype(rtl), eltype(rtc), eltype(rtr), eltype(rbl), eltype(rbc),
        eltype(rbr), eltype(rcp), eltype(ncp), typeof(core_size), typeof(chord))

    return SurfacePanel{TF}(rtl, rtc, rtr, rbl, rbc, rbr, rcp, ncp, core_size, chord)
end

@inline Base.eltype(::Type{SurfacePanel{TF}}) where TF = TF
@inline Base.eltype(::SurfacePanel{TF}) where TF = TF

"""
    top_left(panel::SurfacePanel)

Return the top left vertex of the vortex ring associated with `panel`
"""
@inline top_left(panel::SurfacePanel) = panel.rtl

"""
    top_center(panel::SurfacePanel)

Return the top center vertex of the vortex ring associated with `panel`
"""
@inline top_center(panel::SurfacePanel) = panel.rtc

"""
    top_right(panel::SurfacePanel)

Return the top right vertex of the vortex ring associated with `panel`
"""
@inline top_right(panel::SurfacePanel) = panel.rtr

"""
    bottom_left(panel::SurfacePanel)

Return the bottom left vertex of the vortex ring associated with `panel`
"""
@inline bottom_left(panel::SurfacePanel) = panel.rbl

"""
    bottom_center(panel::SurfacePanel)

Return the bottom center vertex of the vortex ring associated with `panel`
"""
@inline bottom_center(panel::SurfacePanel) = panel.rbc

"""
    bottom_right(panel::SurfacePanel)

Return the bottom right vertex of the vortex ring associated with `panel`
"""
@inline bottom_right(panel::SurfacePanel) = panel.rbr

"""
    controlpoint(panel::SurfacePanel)

Return the control point of `panel` (typically located at the 3/4 chord)
"""
@inline controlpoint(panel::SurfacePanel) = panel.rcp

"""
    normal(panel::SurfacePanel)

Return the normal vector of `panel` at the panel control point
"""
@inline normal(panel::SurfacePanel) = panel.ncp

"""
    get_core_size(panel::SurfacePanel)

Return the core size (smoothing parameter) corresponding to the vortex ring
associated with `panel`
"""
@inline get_core_size(panel::SurfacePanel) = panel.core_size

"""
    translate(panel::SurfacePanel, r)

Return a copy of `panel` translated the distance specified by vector `r`
"""
@inline function translate(panel::SurfacePanel, r)

    rtl = panel.rtl + r
    rtc = panel.rtc + r
    rtr = panel.rtr + r
    rbl = panel.rbl + r
    rbc = panel.rbc + r
    rbr = panel.rbr + r
    rcp = panel.rcp + r
    ncp = panel.ncp
    core_size = panel.core_size
    chord = panel.chord

    return SurfacePanel(rtl, rtc, rtr, rbl, rbc, rbr, rcp, ncp, core_size, chord)
end

"""
    translate(surface, r)

Return a copy of the panels in `surface` translated the distance specified by vector `r`
"""
translate(surface::AbstractMatrix, r) = translate.(surface, Ref(r))

"""
    translate!(surface, r)

Translate the panels in `surface` the distance specified by vector `r`
"""
function translate!(surface::AbstractMatrix, r)

    for i in eachindex(surface)
        surface[i] = translate(surface[i], r)
    end

    return surface
end

"""
    translate(surfaces, r)

Return a copy of the surfaces in `surfaces` translated the distance specified by vector `r`
"""
translate(surfaces::AbstractVector{<:AbstractMatrix}, r) = translate.(surfaces, Ref(r))

"""
    translate!(surfaces, r)

Translate the surfaces in `surfaces` the distance specified by vector `r`
"""
function translate!(surfaces::AbstractVector{AbstractMatrix}, r)

    for i in eachindex(surfaces)
        surfaces[i] = translate!(surfaces[i], r)
    end

    return surfaces
end

"""
    rotate(panel::SurfacePanel, R, r = [0,0,0])

Return a copy of `panel` rotated about point `r` using the rotation matrix `R`
"""
@inline function rotate(panel::SurfacePanel, R, r = (@SVector zeros(3)))

    rtl = R*(panel.rtl - r) + r
    rtc = R*(panel.rtc - r) + r
    rtr = R*(panel.rtr - r) + r
    rbl = R*(panel.rbl - r) + r
    rbc = R*(panel.rbc - r) + r
    rbr = R*(panel.rbr - r) + r
    rcp = R*(panel.rcp - r) + r
    ncp = R*panel.ncp
    core_size = panel.core_size
    chord = panel.chord

    return SurfacePanel(rtl, rtc, rtr, rbl, rbc, rbr, rcp, ncp, core_size, chord)
end

"""
    rotate(surface, R, r = [0,0,0])

Return a copy of the panels in `surface` rotated about point `r` using the rotation matrix `R`
"""
rotate(surface::AbstractMatrix, R, r = (@SVector zeros(3))) = rotate.(surface, Ref(R), Ref(r))

"""
    rotate!(surface, R, r = [0,0,0])

Rotate the panels in `surface` about point `r` using the rotation matrix `R`
"""
function rotate!(surface::AbstractMatrix, R, r = (@SVector zeros(3)))

    for i in eachindex(surface)
        surface[i] = rotate(surface[i], R, r)
    end

    return surface
end

"""
    rotate(surfaces, R, r = [0,0,0])

Return a copy of the surfaces in `surfaces` rotated about point `r` using the
rotation matrix `R`
"""
rotate(surfaces::AbstractVector{<:AbstractMatrix}, R, r = (@SVector zeros(3))) =
    rotate.(surfaces, Ref(R), Ref(r))

"""
    rotate!(surfaces, R, r = [0,0,0])

Rotate the surfaces in `surfaces` about point `r` using the rotation matrix `R`
"""
function rotate!(surfaces::AbstractVector{AbstractMatrix}, R, r = (@SVector zeros(3)))

    for i in eachindex(surfaces)
        surfaces[i] = rotate!(surfaces[i], R, r)
    end

    return surfaces
end

"""
    reflect(panel::SurfacePanel)

Reflect `panel` across the X-Z plane.
"""
@inline function reflect(panel::SurfacePanel)

    rtl = flipy(panel.rtr)
    rtc = flipy(panel.rtc)
    rtr = flipy(panel.rtl)
    rbl = flipy(panel.rbr)
    rbc = flipy(panel.rbc)
    rbr = flipy(panel.rbl)
    rcp = flipy(panel.rcp)
    ncp = flipy(panel.ncp)
    core_size = panel.core_size
    chord = panel.chord

    return SurfacePanel(rtl, rtc, rtr, rbl, rbc, rbr, rcp, ncp, core_size, chord)
end

"""
    reflect(surface)

Reflects the panels in `surface` across the X-Z plane
"""
reflect(surface::AbstractMatrix) = reflect.(reverse(surface, dims=2))

"""
    set_normal(panel::SurfacePanel, ncp)

Return a copy of `panel` with the new normal vector `ncp`
"""
@inline set_normal(panel::SurfacePanel, ncp) = SurfacePanel(panel.rtl, panel.rtc,
    panel.rtr, panel.rbl, panel.rbc, panel.rbr, panel.rcp, ncp, panel.core_size,
    panel.chord)

"""
    left_center(panel::SurfacePanel)

Return the center of the left bound vortex of the vortex ring associated with `panel`
"""
@inline left_center(panel::SurfacePanel) = (top_left(panel) + bottom_left(panel))/2

"""
    right_center(panel::SurfacePanel)

Return the center of the right bound vortex of the vortex ring associated with `panel`
"""
@inline right_center(panel::SurfacePanel) = (top_right(panel) + bottom_right(panel))/2

"""
    top_vector(panel::SurfacePanel)

Return the path of the top bound vortex of the vortex ring associated with `panel`
"""
@inline top_vector(panel::SurfacePanel) = top_right(panel) - top_left(panel)

"""
    left_vector(panel)

Return the path of the left bound vortex of the vortex ring associated with `panel`
"""
@inline left_vector(panel::SurfacePanel) = top_left(panel) - bottom_left(panel)

"""
    right_vector(panel)

Return the path of the right bound vortex of the vortex ring associated with `panel`
"""
@inline right_vector(panel::SurfacePanel) = bottom_right(panel) - top_right(panel)

"""
    bottom_vector(panel)

Return the path of the bottom bound vortex of the vortex ring associated with `panel`
"""
@inline bottom_vector(panel::SurfacePanel) = bottom_left(panel) - bottom_right(panel)

"""
    WakePanel{TF}

SurfacePanel used for modeling wakes.

**Fields**
 - `rtl`: position of the left side of the top bound vortex
 - `rtr`: position of the right side of the top bound vortex
 - `rbl`: position of the left side of the bottom bound vortex
 - `rbr`: position of the right side of the bottom bound vortex
 - `core_size`: finite core size (for use when the finite core smoothing model is enabled)
 - `gamma`: circulation strength of the panel
"""
struct WakePanel{TF}
    rtl::SVector{3, TF}
    rtr::SVector{3, TF}
    rbl::SVector{3, TF}
    rbr::SVector{3, TF}
    core_size::TF
    gamma::TF
end

"""
    WakePanel(rtl, rtr, rbl, rbr, core_size, gamma)

Construct and return a wake panel.

**Arguments**
 - `rtl`: position of the left side of the top bound vortex
 - `rtr`: position of the right side of the top bound vortex
 - `rbl`: position of the left side of the bottom bound vortex
 - `rbr`: position of the right side of the bottom bound vortex
 - `core_size`: finite core size
 - `gamma`: circulation strength of the panel
"""
WakePanel(args...; kwargs...)

function WakePanel(rtl, rtr, rbl, rbr, core_size, gamma)

    TF = promote_type(eltype(rtl), eltype(rtr), eltype(rbl), eltype(rbr), typeof(core_size), typeof(gamma))

    return WakePanel{TF}(rtl, rtr, rbl, rbr, core_size, gamma)
end

@inline Base.eltype(::Type{WakePanel{TF}}) where TF = TF
@inline Base.eltype(::WakePanel{TF}) where TF = TF

@inline top_left(panel::WakePanel) = panel.rtl

@inline top_right(panel::WakePanel) = panel.rtr

@inline bottom_left(panel::WakePanel) = panel.rbl

@inline bottom_right(panel::WakePanel) = panel.rbr

@inline get_core_size(panel::WakePanel) = panel.core_size

"""
    circulation_strength(panel::WakePanel)

Return the circulation strength of the wake panel.
"""
@inline circulation_strength(panel::WakePanel) = panel.gamma

"""
    set_circulation_strength(panel::WakePanel, gamma)

Return a copy of `panel` with the circulation strength `gamma`
"""
@inline function set_circulation_strength(panel::WakePanel, gamma)

    rtl = panel.rtl
    rtr = panel.rtr
    rbl = panel.rbl
    rbr = panel.rbr
    core_size = panel.core_size

    return WakePanel(rtl, rtr, rbl, rbr, core_size, gamma)
end

@inline function translate(panel::WakePanel, r)

    rtl = panel.rtl + r
    rtr = panel.rtr + r
    rbl = panel.rbl + r
    rbr = panel.rbr + r
    core_size = panel.core_size
    gamma = panel.gamma

    return WakePanel(rtl, rtr, rbl, rbr, core_size, gamma)
end

@inline function reflect(panel::WakePanel)

    rtl = flipy(panel.rtr)
    rtr = flipy(panel.rtl)
    rbl = flipy(panel.rbr)
    rbr = flipy(panel.rbl)
    core_size = panel.core_size
    gamma = panel.gamma

    return Ring(rtl, rtr, rbl, rbr, core_size, gamma)
end

"""
    TrefftzPanel{TF}

Panel in the Trefftz plane.
"""
struct TrefftzPanel{TF}
    rl::SVector{3,TF}
    rc::SVector{3,TF}
    rr::SVector{3,TF}
    Γ::TF
end

@inline Base.eltype(::Type{TrefftzPanel{TF}}) where TF = TF
@inline Base.eltype(::TrefftzPanel{TF}) where TF = TF

"""
    normal(panel::TrefftzPanel)

Return the normal vector of `panel`, including magnitude
"""
@inline function normal(panel::TrefftzPanel)

    rl = panel.rl
    rr = panel.rr

    dy = rr[2] - rl[2]
    dz = rr[3] - rl[3]

    nhat = SVector(0, -dz, dy)

    return nhat
end

"""
    trefftz_panels(surface[s], freestream, Γ)

Constructs a set of panels for Trefftz plane calculations
"""
trefftz_panels

function trefftz_panels(surface::AbstractMatrix, fs, Γ)

    TF = eltype(eltype(panels))

    panels = Vector{TrefftzPanel{TF}}(undef, size(surface, 2))

    return trefftz_panels!(panels, surface, fs, Γ)
end

function trefftz_panels(surfaces::AbstractVector{<:AbstractMatrix}, fs, Γ)

    TF = eltype(eltype(eltype(surfaces)))

    panels = [Vector{TrefftzPanel{TF}}(undef, size(surfaces[i], 2)) for i = 1:length(surfaces)]

    return trefftz_panels!(panels, surfaces, fs, Γ)
end

# single surface
function trefftz_panels!(panels, surface::AbstractMatrix, fs, Γ)

    n1, n2 = size(surface)

    rΓ = reshape(Γ, n1, n2)

    R = body_to_wind(fs)

    for j = 1:n2
        # use trailing edge panel for geometry
        rl = bottom_left(surface[end,j])
        rc = bottom_center(surface[end,j])
        rr = bottom_right(surface[end,j])

        # panel circulation is trailing edge panel circulation
        Γt = rΓ[end,j]

        # rotate into Trefftz plane frame of reference
        rl = R*rl
        rc = R*rc
        rr = R*rr

        # zero out x-components
        rl = SVector(0, rl[2], rl[3])
        rc = SVector(0, rc[2], rc[3])
        rr = SVector(0, rr[2], rr[3])

        panels[j] = TrefftzPanel(rl, rc, rr, Γt)
    end

    return panels
end

# multiple surfaces
function trefftz_panels!(panels, surfaces::AbstractVector{<:AbstractMatrix}, fs, Γ)

    nsurf = length(surfaces)

    # indices for keeping track of circulation
    iΓ = 0

    # loop through surfaces
    for i = 1:nsurf

        # number of panels on this surface
        n = length(surfaces[i])

        # create view into `panels` and `Γ` vector
        vΓ = view(Γ, iΓ+1:iΓ+n)

        # populate vector of Trefftz panels
        trefftz_panels!(panels[i], surfaces[i], fs, vΓ)

        # increment position in circulation index
        iΓ += n
    end

    return panels
end

"""
    trefftz_panel_induced_drag(receiving::TrefftzPanel, sending::TrefftzPanel; kwargs...)

Induced drag on `receiving` panel induced by `sending` panel.

# Keyword Arguments
 - `symmetric`: Flag indicating whether a mirror image of `sending` should be
    used when calculating the induced drag
"""
@inline function trefftz_panel_induced_drag(receiving::TrefftzPanel,
    sending::TrefftzPanel; symmetric)

    rl = sending.rl
    rr = sending.rr
    Γs = sending.Γ

    rc = receiving.rc
    nc = normal(receiving)
    Γr = receiving.Γ

    Di = vortex_induced_drag(rl, -Γs, rc, Γr, nc)
    Di += vortex_induced_drag(rr, Γs, rc, Γr, nc)
    if symmetric
        Di += vortex_induced_drag(flipy(rr), -Γs, rc, Γr, nc)
        Di += vortex_induced_drag(flipy(rl), Γs, rc, Γr, nc)
    end

    return Di
end

"""
    vortex_induced_drag(rj, Γj, ri, Γi, ni)

Return induced drag from vortex `j` induced on panel `i`
"""
@inline function vortex_induced_drag(rj, Γj, ri, Γi, ni)

    rij = ri - rj
    Vthetai = SVector(0, -Γj*rij[3], Γj*rij[2]) / (2*pi*(rij[2]^2 + rij[3]^2))
    Vn = -dot(Vthetai, ni)

    Di = RHO/2.0*Γi*Vn

    return Di
end


text:
using DataFrames
# using Ipaper

DataFrames.nrow(x::AbstractArray) = size(x, 1)
DataFrames.ncol(x::AbstractArray) = size(x, 2)

# rbind(args...) = cat(args..., dims = 1)
# cbind(args...) = cat(args..., dims = 2)
abind(args...; along=3) = cat(args..., dims=along)

# rbind = vcat
rbind(args...; kw...) = vcat(args...; kw...)
rbind(x) = x
rbind(x::DataFrame,
    y::Union{DataFrame,AbstractVecOrMat}; kw...) = begin
    # @assert (ncol(x) == ncol(y))
    x = as_dataframe(x)
    y = as_dataframe(y, names(x))
    vcat(x, y; kw...)
end

rbind(x::AbstractVecOrMat, y::DataFrame; kw...) = rbind(as_dataframe(x, names(y)), y; kw...)

function rbind(x::DataFrame, args...; kw...)
    x = as_dataframe(x)
    if length(args) == 0
        x
    elseif length(args) == 1
        rbind(x, args[1])
    else
        rbind(rbind(x, args[1]), args[2:end]...)
    end
end

# cbind = hcat # not work
cbind(args...; kw...) = hcat(args...; kw...)
cbind(x) = x
cbind(x::DataFrame, y::Union{DataFrame,AbstractVecOrMat}) =
    hcat(as_dataframe(x), as_dataframe(y); makeunique=true)
cbind(x::AbstractVecOrMat, y::DataFrame; kw...) = cbind(as_dataframe(x), y; kw...)

# by reference
function cbind(x::DataFrame, args...; kw...)
    x = as_dataframe(x)
    n = length(kw)
    if n > 0
        vars = keys(kw)
        for i = 1:n
            key = vars[i]
            val = kw[i]
            # @show key
            # @show val
            if !isa(val, AbstractArray) || length(val) == 1
                x[:, key] .= val
            else
                x[:, key] = val
            end
        end
    end
    if length(args) == 0
        x
    elseif length(args) == 1
        cbind(x, args[1])
    else
        cbind(cbind(x, args[1]), args[2:end]...)
    end
end

macro as_df(x)
    name = string(x)
    expr = :(DataFrame($name => $x))
    esc(expr)
end
export @as_df;

as_matrix(x::DataFrame) = Matrix(x)

as_dataframe(x::DataFrame, args...) = x
as_dataframe(x::AbstractVector) = @as_df(x)

as_dataframe(x::AbstractMatrix) = DataFrame(x, :auto)
function as_dataframe(x::AbstractVecOrMat, names::AbstractVector; kw...)
    DataFrame(x, names; kw...)
end

is_dataframe(d) = d isa DataFrame

# for data.frame by reference operation
function melt_list(list; kw...)
    if length(kw) > 0
        by = keys(kw)[1]
        vals = kw[1]
    else
        by = :I
        vals = 1:length(list)
    end

    for i = 1:length(list)
        d = list[i]
        if (d isa DataFrame)
            d[:, by] .= vals[i]
        end
    end
    ind = map(is_dataframe, list)
    rbind(list[ind]...)
end

# seealso: leftjoin, rightjoin, innerjoin, outerjoin
function dt_merge(x::DataFrame, y::DataFrame; by=nothing,
    all=false, all_x=all, all_y=all, makeunique=true, kw...)

    if by === nothing
        by = intersect(names(x), names(y))
    end
    if !all
        if all_x
            leftjoin(x, y; on=by, makeunique=true, kw...)
        elseif all_y
            rightjoin(x, y; on=by, makeunique=true, kw...)
        else
            # all_x = f && all_y = f
            innerjoin(x, y; on=by, makeunique=true, kw...)
        end
    else
        outerjoin(x, y; on=by, makeunique=true, kw...)
    end
end

fread(f) = DataFrame(CSV.File(f))
fwrite(df, file) = begin
    dirname(file) |> check_dir
    CSV.write(file, df)
end

# for data.frame by reference operation
function dataframe(; kw...)
    DataFrame(pairs(kw))
end
datatable = dataframe

function list(; kw...)
    Dict(pairs(kw))
end

#! This version not work
# function datatable(args...; kw...)
#     params = args..., kw...
#     datatable(; params...)
# end
const DF = dataframe;

export rbind, cbind, abind, melt_list,
    fread, fwrite, dt_merge,
    is_dataframe,
    as_dataframe,
    as_matrix, nrow, ncol,
    DataFrame, DF, names,
    list,
    datatable, dataframe

ext: jl


text:
using CartesianSphericalHarmonics
using Test

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


text:
# The Computer Language Benchmarks Game
# https://salsa.debian.org/benchmarksgame-team/benchmarksgame/
#
# Contributed by Adam Beckmeyer

module NBodyAB

import Printf: @printf
import Base: @propagate_inbounds

const DAYS_PER_YEAR = 365.24
const SOLAR_MASS = 4 * π * π
const DEFAULT_STEP = 0.01

# Use a tuple rather than struct to hint compiler towards simd code
# 4 floats in tuple instead of 3 generates better SIMD instructions
V3d(x=0.0, y=0.0, z=0.0) = (Float64(x), Float64(y), Float64(z), 0.0)

# Rust implementation #7 uses pos and v of length 3 with 1 filler
# Float64 field for a total struct width of 64 bytes instead of 72
# bytes as shown currently. This might be possible for Julia in the
# future using https://github.com/eschnett/SIMD.jl. Not sure if it
# would actually be advantageous though.
struct Body
    pos::NTuple{4,Float64}
    v::NTuple{4,Float64}
    m::Float64
end#struct

# don't bother adding the empty 4th element
Base.sum(v::NTuple{4,Float64}) = @inbounds +(v[1], v[2], v[3])

function init_sun(bodies)
    p = V3d()
    for b in bodies
        p = muladd.(b.v, b.m, p)
    end#for
    Body(V3d(), p .* (-inv(SOLAR_MASS)), SOLAR_MASS)
end#function

@propagate_inbounds function update_velocity(b1::Body, b2::Body, Δt::Float64)
    Δpos = b1.pos .- b2.pos
    d² = sum(Δpos .* Δpos)

    # Fastest implementations use an intrinsic to do a
    # single-precision sqrt approximation followed by two iterations
    # of the Newton-Raphson method. But the sqrt isn't currently this
    # implementations bottleneck, so I've left it unoptimized.

    # The fastest Rust/C implementations store this value instead of
    # immediately modifying the body objects. Possible performance
    # gains there due to cache-locality? I've been unable to replicate
    # this in Julia by calculating mag as an NTuple{10,Float64}.
    vmag = Δpos .* (Δt / (d² * √d²))

    (Body(b1.pos, muladd.(-b2.m, vmag, b1.v), b1.m),
     Body(b2.pos, muladd.( b1.m, vmag, b2.v), b2.m))
end#function

@propagate_inbounds update_pos(b::Body, Δt::Float64) =
    Body(muladd.(b.v, Δt, b.pos), b.v, b.m)

function next!(bodies, Δt)
    # Using special iteration tools like eachindex or Iterators.drop
    # is actually measurably slower than for loop around a UnitRange
    @inbounds for i in 1:length(bodies), j in i+1:length(bodies)
        # Script spends roughly 90% of it's time in this loop. Of that
        # 90%, 36% is spent on float multiplication, 34% is spent on
        # the call to muladd, 12% is spent on vector manipulation and
        # iteration overhead, 11% is spent on division, 2% is spent on
        # sqrt, 2% is spent on summing. The remaining 3% is
        # unaccounted for. Overall, multiplication (including muladd)
        # dominates this benchmark.
        bodies[i], bodies[j] = update_velocity(bodies[i], bodies[j], Δt)
    end#for

    # this could be broadcast or mapped, but for-loop is faster
    @inbounds for i in 1:length(bodies)
        bodies[i] = update_pos(bodies[i], Δt)
    end#for
end#function

function energy(bodies)
    e = 0.0
    @inbounds for i in 1:length(bodies)
        bi = bodies[i]
        e += 0.5 * bi.m * sum(bi.v .* bi.v)
        for j in i+1:length(bodies)
            bj = bodies[j]
            Δpos = bi.pos .- bj.pos
            d = √sum(Δpos .* Δpos)
            e = muladd(bi.m * bj.m, -inv(d), e)
        end#for
    end#for
    e
end#function

function main(io, n, Δt=DEFAULT_STEP)
    jupyter = Body(V3d(4.84143144246472090e+00,
                       -1.16032004402742839e+00,
                       -1.03622044471123109e-01),
                   V3d(1.66007664274403694e-03 * DAYS_PER_YEAR,
                       7.69901118419740425e-03 * DAYS_PER_YEAR,
                       -6.90460016972063023e-05 * DAYS_PER_YEAR),
                   9.54791938424326609e-04 * SOLAR_MASS)
    saturn = Body(V3d(8.34336671824457987e+00,
                      4.12479856412430479e+00,
                      -4.03523417114321381e-01),
                  V3d(-2.76742510726862411e-03 * DAYS_PER_YEAR,
                      4.99852801234917238e-03 * DAYS_PER_YEAR,
                      2.30417297573763929e-05 * DAYS_PER_YEAR),
                  2.85885980666130812e-04 * SOLAR_MASS)
    uranus = Body(V3d(1.28943695621391310e+01,
                      -1.51111514016986312e+01,
                      -2.23307578892655734e-01),
                  V3d(2.96460137564761618e-03 * DAYS_PER_YEAR,
                      2.37847173959480950e-03 * DAYS_PER_YEAR,
                      -2.96589568540237556e-05 * DAYS_PER_YEAR),
                  4.36624404335156298e-05 * SOLAR_MASS)
    neptune = Body(V3d(1.53796971148509165e+01,
                       -2.59193146099879641e+01,
                       1.79258772950371181e-01),
                   V3d(2.68067772490389322e-03 * DAYS_PER_YEAR,
                       1.62824170038242295e-03 * DAYS_PER_YEAR,
                       -9.51592254519715870e-05 * DAYS_PER_YEAR),
                   5.15138902046611451e-05 * SOLAR_MASS)
    sun = init_sun((jupyter, saturn, uranus, neptune))
    # Although bodies is fixed size, vector instead of tuple for mutability
    bodies = [sun, jupyter, saturn, uranus, neptune]

    @printf(io, "%.9f\n", energy(bodies))
    for i in 1:n
        next!(bodies, Δt)
    end#for
    @printf(io, "%.9f\n", energy(bodies))
end#function

end#module

isinteractive() || NBodyAB.main(stdout, parse(Int, ARGS[1]))


text:
using PyCall
unshift!(PyVector(pyimport("sys")["path"]), "");
@pyimport matplotlib.pyplot as plt
@pyimport numpy as np
@pyimport json
@pyimport os
@pyimport pickle
@pyimport pandas as pd
@pyimport collections

@pyimport parameters_julia
@pyimport utils_julia

@pyimport OD_functions
@pyimport runGLS_f

out_dir = parameters_julia.out_dir;
files_ID = parameters_julia.files_ID;
month_w = parameters_julia.month_w;
year = parameters_julia.year;
time_instances = parameters_julia.instances_ID;
deg_grid = parameters_julia.deg_grid;
c_grid = parameters_julia.c_grid;
lamb_grid = parameters_julia.lamb_grid;
week_day_list = parameters_julia.week_day_list;
average_over_time_p = parameters_julia.average_over_time_p;
od_pairs = parameters_julia.od_pairs;

using JuMP, Ipopt
function GLS_p2(xi_list, P)
	mGLSJulia = Model(solver=IpoptSolver(max_iter=5000))

	@variable(mGLSJulia, lam[1:size(P,1)] >= 0)
	@variable(mGLSJulia, p[1:size(P,1), 1:size(P,2)] >= 0)

	for i = 1:size(P,1)
	    for j = 1:size(P,2)
	        if P[i,j] == 0
	            @constraint(mGLSJulia, p[i,j] == 0)
	        end
	    end
	end

	for i = 1:size(P,1)
	    @constraint(mGLSJulia, sum{p[i,j], j = 1:size(P,2)} == 1)
	end

	for l = 1:size(P,2)
	    @NLconstraint(mGLSJulia, sum{p[i,l] * lam[i], i = 1:size(P,1)} == xi_list[l])
	end

	@objective(mGLSJulia, Min, sum{lam[j], j = 1:size(P,1)})  # play no actual role, but could not use zero objective
	#@NLobjective(mGLSJulia, Min, obj)
	#@NLobjective(mGLSJulia, Min, sum{ (sum{p[i,l] * lam[i], i = 1:size(P,1)} - xi_list[l])^2, l=1:size(P,1)})
	solve(mGLSJulia)

	return getvalue(lam), getvalue(p), getobjectivevalue(mGLSJulia)

end

function saveDemandVec(od_pairs, file_name, lam_list)
	lam = Dict()
	idx = 0
	open(file_name, "w") do file
		for od  = 1: size(od_pairs,1)
			i = od_pairs[od,1]
			j = od_pairs[od,2]
			idx += 1
			lam[idx] = lam_list[idx]
			dem = lam_list[idx]
			write(file, "$i, $j, $dem")#,"\n")
			write(file, "\n")
			#for ln =1:size(lam,1)
			#	println("$(i),$(j),$(lam_list[idx])")
			#end
		end
	close(file)
	end
end

#=
# Solve P1 and P2 of GLS for every day and every instance
for instance in time_instances
	for day in week_day_list
		xi_list, x, A, P, gls_cost = runGLS_f.GLS_p1(instance, day, average_over_time_p);
		P[P.>0] = 1;
		lam_, Pr, obj  = GLS_p2(xi_list, P);
		file_name = out_dir * "OD_demands/OD_demand_matrix_" * month_w * "_" * string(day) *  "_weekday_" * instance * files_ID * ".txt"
		saveDemandVec(od_pairs, file_name, lam_)
		P_file_name = out_dir * "OD_demands/OD_route_prob_matrix_" * month_w * "_" * string(day) *  "_weekday_" * instance * files_ID * ".txt"
		writedlm(P_file_name, Pr)
		# Calculate flows APg
		flow_file = out_dir * "OD_demands/flow_APg_" * month_w * "_" * string(day) *  "_weekday_" * instance * files_ID * ".txt"
		writedlm(flow_file, A*Pr'*lam_ )
		# Store data flows flows
		flow_file = out_dir * "OD_demands/flow_data_" * month_w * "_" * string(day) *  "_weekday_" * instance * files_ID * ".txt"
		writedlm(flow_file, mean(x,2) )
		#gls_name = out_dir * "OD_demands/gls_cost_vec_" * month_w * "_weekday_" * instance * files_ID * ".json"
		#outfile = open(gls_name, "w")
	    #JSON.print(outfile, gls_cost)
	    #close(outfile)
	end
end
=#

# Solve P1 and P2 of GLS for every instance using ALL observations
for instance in time_instances
	xi_list, x, A, P, gls_cost = runGLS_f.GLS_p1_all(instance, week_day_list, average_over_time_p, "all");
	P[P.>0] = 1;
	lam_, Pr, obj  = GLS_p2(xi_list, P);

	for i in ["all", "full"]
		file_name = out_dir * "OD_demands/OD_demand_matrix_" * month_w * "_" * i *  "_weekday_" * instance * files_ID * ".txt"
		saveDemandVec(od_pairs, file_name, lam_)

		P_file_name = out_dir * "OD_demands/OD_route_prob_matrix_" * month_w * "_" * i *  "_weekday_" * instance * files_ID * ".txt"
		writedlm(P_file_name, Pr)

		# Calculate flows APg
		flow_file = out_dir * "OD_demands/flow_APg_" * month_w * "_" * i *  "_weekday_" * instance * files_ID * ".txt"
		writedlm(flow_file, A*Pr'*lam_ )

		# Store data flows flows
		flow_file = out_dir * "OD_demands/flow_data_" * month_w * "_" * i *  "_weekday_" * instance * files_ID * ".txt"
		writedlm(flow_file, mean(x,2) )

		gls_name = out_dir * "OD_demands/gls_cost_vec_" * month_w * "_weekday_" * instance * files_ID * ".json"
		outfile = open(gls_name, "w")
	    JSON.print(outfile, gls_cost)
	    close(outfile)
	end
end
#=
# Solve P1 and P2 of GLS for every instance using daily averages
for instance in time_instances
	xi_list, x, A, P = runGLS_f.GLS_p1_all(instance, week_day_list, average_over_time_p, "all")
	P[P.>0] = 1;
	lam_, Pr, obj  = GLS_p2(xi_list, P);

	file_name = out_dir * "OD_demands/OD_demand_matrix_" * month_w * "_" * "full" *  "_weekday_" * instance * files_ID * ".txt"
	saveDemandVec(od_pairs, file_name, lam_)

	file_name = out_dir * "OD_demands/OD_demand_matrix_" * month_w * "_" * "all" *  "_weekday_" * instance * files_ID * ".txt"
	saveDemandVec(od_pairs, file_name, lam_)

end
=#


text:
module MagicaVoxel # VoxelPlots

include("MagicaVoxel/parser.jl")

end # module VoxelPlots.MagicaVoxel


text:
# Copyright 2021 The CLVR Authors
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
#     http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.

"""Data structure for formalizing a standard-form linear program."""
struct StandardLinearProgram

    A_T::SparseMatrixCSC{Float64, Int64}
    b::Vector{Float64}
    c::Vector{Float64}
    prox

    function StandardLinearProgram(A_T, b, c)
        prox(x, τ) = max.(x, 0.0)
        new(sparse(A_T), Vector{Float64}(b), Vector{Float64}(c), prox)
    end
end

max_stars_repo_name: ericlincc/Efficient-GLP

max_issues_repo_name: ericlincc/Efficient-GLP

max_forks_repo_name: ericlincc/Efficient-GLP


text:
using DandelionWebSockets: FrameWriter, CLOSE_STATUS_NO_STATUS, CLOSE_STATUS_GOING_AWAY, closesocket
import DandelionWebSockets: sendcloseframe

@testset "Frame writer           " begin
    @testset "Write close frame with CLOSE_STATUS_NO_STATUS; No body in the frame" begin
        frameio = FrameIOStub()
        rng = FakeRNG{UInt8}(b"\x01\x02\x03\x04")
        framewriter = FrameWriter(frameio, rng)

        sendcloseframe(framewriter, CLOSE_STATUS_NO_STATUS)

        frame = getframe(frameio, 1)
        @test isempty(frame.payload)
    end

    @testset "Write close frame with CLOSE_STATUS_GOING_AWAY; Body starts with that code" begin
        frameio = FrameIOStub()
        mask = b"\x01\x02\x03\x04"
        rng = FakeRNG{UInt8}(mask)
        framewriter = FrameWriter(frameio, rng)

        sendcloseframe(framewriter, CLOSE_STATUS_GOING_AWAY)

        frame = getframeunmasked(frameio, 1, mask)
        payload = IOBuffer(frame.payload)
        @test read(payload, UInt16) == hton(CLOSE_STATUS_GOING_AWAY.code)
    end

    @testset "Write close frame with reason; Body has text included" begin
        frameio = FrameIOStub()
        mask = b"\x01\x02\x03\x04"
        rng = FakeRNG{UInt8}(mask)
        framewriter = FrameWriter(frameio, rng)
        reason = "Some reason text"

        sendcloseframe(framewriter, CLOSE_STATUS_GOING_AWAY; reason=reason)

        frame = getframeunmasked(frameio, 1, mask)

        # Check first that length matches the status code plus reason
        @test length(frame.payload) == 2 + length(reason) # Status code and reason

        # Skip the 16 bit status code and read the reason
        payload = IOBuffer(frame.payload)
        skip(payload, 2)
        @test read(payload, length(reason)) == unsafe_wrap(Vector{UInt8}, (reason))
    end

    @testset "Write close frame with no status, but reason; No reason is included" begin
        frameio = FrameIOStub()
        rng = FakeRNG{UInt8}(b"\x01\x02\x03\x04")
        framewriter = FrameWriter(frameio, rng)

        sendcloseframe(framewriter, CLOSE_STATUS_NO_STATUS; reason = "Some reason")

        frame = getframe(frameio, 1)
        @test isempty(frame.payload)
    end

    @testset "Closing the socket; Underlying stream is not open" begin
        frameio = FrameIOStub()
        rng = FakeRNG{UInt8}(b"")
        framewriter = FrameWriter(frameio, rng)

        closesocket(framewriter)

        @test frameio.isopen == false
    end
end


text:
import Base.dump

export dump, load_data, load_data!, load_init

"""
    dump(file_name :: String, data :: Data)

Save `data` to the file with name `file_name` (and extension ".h5").

"""
function dump( h5file :: String, data :: Data )
    h5open(joinpath(h5file * ".h5"), "cw") do file  #"w" for write ("r" for read)
        write(file,  "/data/size", collect(size(first(data.U))) )
        write(file, "/data/length", length(data.U) )

        for (i,u) in enumerate(data.U)
            write(file, "/data/U$i", u )
        end
    end

end

"""
    load_data(file_name :: String)

Load data from the file with name `file_name` (and extension ".h5").

"""
function load_data( h5file :: String )

    nsteps = h5read(h5file * ".h5", "/data/length")
    np, nv = h5read(h5file * ".h5", "/data/size")

    v = zeros(ComplexF64, np, nv)

    data = Data(v)

    @assert length( data ) == 1
    @assert np == size(first(data.U))[1]
    @assert nv == size(first(data.U))[2]

    data.U[1] .= h5read(h5file * ".h5", "/data/U1")
    for i in 2:nsteps
        push!(data.U, h5read(h5file * ".h5", "/data/U$i"))
    end

    return data

end

"""
    dump(file_name  :: String, problem :: Problem)

Save the data of `problem` to the file with name `file_name` (and extension ".h5").

"""
function dump( h5file :: String, pb :: Problem)
    dump( h5file, pb.data )
end

"""
    load_data!(file_name :: String, problem :: Problem)

Fills `problem` with raw data extracted from the file with name `file_name` (and extension ".h5").

"""
function load_data!( h5file :: String, pb::Problem )
    data = load_data( h5file )
    pb.data = data
end

"""
    dump(file_name :: String, x::Vector, init::InitialData )

Save the values of the initial data `init` at collocation points `x` to the file with name `file_name` (and extension ".h5").

"""
function dump( h5file :: String, x::Vector, init::InitialData)
    h5open(joinpath(h5file * ".h5"), "cw") do file  #"w" for write ("r" for read)

        write(file, "/init/x", x )
        write(file, "/init/η", init.η(x) )
        write(file, "/init/v", init.v(x) )
        write(file, "/init/label", init.label)

    end
end

"""
    load_init(file_name :: String; fast = false)

Load initial data from  the file with name `file_name` (and extension ".h5").

Keyword argument `fast` is optional (default is `false`), and corresponds the the keyword argument of the function `Init`.

Return an object of type `InitialData`.

"""
function load_init( h5file :: String ; fast = false)

    x = h5read(joinpath( h5file * ".h5"), "/init/x")
    η = h5read(joinpath( h5file * ".h5"), "/init/η")
    v = h5read(joinpath( h5file * ".h5"), "/init/v")
    label = h5read(joinpath(h5file * ".h5"), "/init/label")

    return Init(x,η,v ; fast = fast, label = label)
end


text:
## Preamble
using Base.Test

using CompScienceMeshes
using BEAST

## The actual tests

κ = ω = 1.0

Γ = meshsegment(1.0, 0.5)
X = lagrangec0d1(Γ)
@test numvertices(Γ)-2 == numfunctions(X)

hypersingular = HyperSingular(κ)
identityop    = Identity()
doublelayer   = DoubleLayer(κ)

@time N = assemble(hypersingular, X, X)
@time I = assemble(identityop, X, X)

@test size(N) == (numfunctions(X), numfunctions(X))
@test size(I) == (numfunctions(X), numfunctions(X))
@test rank(I) == numfunctions(X)

@time e = assemble(PlaneWaveNeumann(κ, point(0.0, 1.0)), X)
@test length(e) == numfunctions(X)

x1 = N \ e;

# Testing duallagrangec0d1
Γ1 = meshcircle(1.0, 2.5,3)          # creating a triangle
Γ2 = barycentric_refinement(Γ1)      # creating the refined mesh

X1 =duallagrangec0d1(Γ1,Γ2)          # creating the basis functions
X1 = duallagrangec0d1(Γ1, Γ2, x->false, Val{2})
@test numcells(Γ1) == numfunctions(X1)  # making sure it is assigned according to the coarse mesh segments
@test length(X1.fns[1])== 6           # making sure each segment represent 6 shapes inside it
@test length(X1.fns[numfunctions(X1)])== 6  # making sure the last segment functions contains 6 shapes as well

## Test linear Lagrange elements on triangles and the computation of their curl
T = Float64
Degr = 1
Dim1 = 3
NumF = 3
f = BEAST.LagrangeRefSpace{T,Degr,Dim1,NumF}()
sphere = readmesh(joinpath(dirname(@__FILE__),"assets","sphere5.in"))
s = chart(sphere, first(cells(sphere)))
t = neighborhood(s, [1,1]/3)
v = f(t, Val{:withcurl})

A = volume(s)
@test v[1][2] == (s[3]-s[2])/2A
@test v[2][2] == (s[1]-s[3])/2A
@test v[3][2] == (s[2]-s[1])/2A

@test v[1][1] ≈ 1/3
@test v[2][1] ≈ 1/3
@test v[3][1] ≈ 1/3

## Test the construction of continuous linear Lagrange elements on 2D surfaces
using CompScienceMeshes
using BEAST
using Base.Test

m = meshrectangle(1.0, 1.0, 0.5, 3)
X = lagrangec0d1(m)
x = refspace(X)

@test numfunctions(x) == 3

@test numfunctions(X) == 1
@test length(X.fns[1]) == 6

## test the scalar trace for Lagrange functions
using CompScienceMeshes
using BEAST
using Base.Test

p1 = point(0,0,0)
p2 = point(1,0,0)
p3 = point(0,1,0)
p4 = point(1,1,1)
c1 = index(1,2,3)

m = Mesh([p1,p2,p3,p4],[c1])
b = Mesh([p1,p2], [index(1,2)])
X = lagrangec0d1(m, boundary(m))
@test numfunctions(X) == 3

Y = BEAST.strace(X, b)

@test numfunctions(X) == 3
@test numfunctions(Y) == 3

@test length(Y.fns[1]) == 1
@test length(Y.fns[2]) == 1
@test length(Y.fns[3]) == 0

sh = Y.fns[1][1]; @test (sh.cellid, sh.refid, sh.coeff) == (1, 1, 1.0)
sh = Y.fns[2][1]; @test (sh.cellid, sh.refid, sh.coeff) == (1, 2, 1.0)

x = refspace(X)

cell = chart(m, first(cells(m)))
face = chart(b, first(cells(b)))
Q = BEAST.strace(x, cell, 3, face)
@test Q == [1 0 0; 0 1 0]

## test Lagrange construction on Junctions
using CompScienceMeshes
using BEAST
using Base.Test

m1 = meshrectangle(1.0, 0.5, 0.5)
m2 = rotate(m1, 0.5π*[1,0,0])
m3 = rotate(m1, 1.0π*[1,0,0])
m = weld(m1, m2, m3)

X = lagrangec0d1(m)
x = refspace(X)
@test numfunctions(X) == 1
@test length(X.fns[1]) == 9

p = point(0.5, 0.0, 0.0)
for s in X.fns[1]
    cell = m.faces[s.cellid]
    patch = chart(m, cell)
    bary = carttobary(patch, p)
    mp = neighborhood(patch, bary)

end

## Test the dual pieweise constant lagrange elemetns
using CompScienceMeshes
using BEAST
using Base.Test

width, height = 1.0, 1.0
h = 0.5

m = meshrectangle(width, height, h)
b = meshsegment(width, width, 3)

X = duallagrangecxd0(m, b)
@test numfunctions(X) == 2
@test isa(refspace(X), BEAST.LagrangeRefSpace)

@test length(X.fns[1]) == 6
@test length(X.fns[2]) == 12

fine = geometry(X)
for fn in X.fns
    n = length(fn)
    for sh in fn
        @test sh.refid == 1
        cellid = sh.cellid
        cell = fine.faces[cellid]
        ptch = chart(fine, cell)
        @test sh.coeff * volume(ptch) ≈ 1/n
    end
end

## Test the construction of dual piecewise linear, globally continuous elements
using CompScienceMeshes
using BEAST
using Base.Test

m = meshrectangle(1.0, 1.0, 0.25)
j = meshsegment(1.0, 1.0, 3)

X = duallagrangec0d1(m,j)
x = refspace(X)

@test numfunctions(X) == numcells(m)

isonjunction = inclosure_gpredicate(j)
els, ad = BEAST.assemblydata(X)
for p in 1:numcells(m)
    el = els[p]
    for r in 1:numfunctions(x)
        vert = el[r]
        isonjunction(vert) || continue
        for (i,w) in ad[p,r]
            @test w == 0
        end
    end
end

## end of file


text:
#=
This file is auto-generated. Do not edit.
=#

#! format: off

"""
    mutable struct KauraPLL <: FrequencyEstimator
        ω_lp::Float64
        kp_pll::Float64
        ki_pll::Float64
        ext::Dict{String, Any}
        states::Vector{Symbol}
        n_states::Int
    end

Parameters of a Phase-Locked Loop (PLL) based on Kaura, Vikram, and Vladimir Blasko.
"Operation of a phase locked loop system under distorted utility conditions."
IEEE Transactions on Industry applications 33.1 (1997): 58-63.

# Arguments
- `ω_lp::Float64`: PLL low-pass filter frequency (rad/sec), validation range: `(0, nothing)`
- `kp_pll::Float64`: PLL proportional gain, validation range: `(0, nothing)`
- `ki_pll::Float64`: PLL integral gain, validation range: `(0, nothing)`
- `ext::Dict{String, Any}`
- `states::Vector{Symbol}`: The states of the KauraPLL model are:
	vd_pll: d-axis of the measured voltage in the PLL synchronous reference frame (SRF),
	vq_pll: q-axis of the measured voltage in the PLL SRF,
	ε_pll: Integrator state of the PI controller,
	θ_pll: Phase angle displacement in the PLL SRF
- `n_states::Int`: KauraPLL has 4 states
"""
mutable struct KauraPLL <: FrequencyEstimator
    "PLL low-pass filter frequency (rad/sec)"
    ω_lp::Float64
    "PLL proportional gain"
    kp_pll::Float64
    "PLL integral gain"
    ki_pll::Float64
    ext::Dict{String, Any}
    "The states of the KauraPLL model are:
	vd_pll: d-axis of the measured voltage in the PLL synchronous reference frame (SRF),
	vq_pll: q-axis of the measured voltage in the PLL SRF,
	ε_pll: Integrator state of the PI controller,
	θ_pll: Phase angle displacement in the PLL SRF"
    states::Vector{Symbol}
    "KauraPLL has 4 states"
    n_states::Int
end

function KauraPLL(ω_lp, kp_pll, ki_pll, ext=Dict{String, Any}(), )
    KauraPLL(ω_lp, kp_pll, ki_pll, ext, [:vd_pll, :vq_pll, :ε_pll, :θ_pll], 4, )
end

function KauraPLL(; ω_lp, kp_pll, ki_pll, ext=Dict{String, Any}(), states=[:vd_pll, :vq_pll, :ε_pll, :θ_pll], n_states=4, )
    KauraPLL(ω_lp, kp_pll, ki_pll, ext, states, n_states, )
end

# Constructor for demo purposes; non-functional.
function KauraPLL(::Nothing)
    KauraPLL(;
        ω_lp=0,
        kp_pll=0,
        ki_pll=0,
        ext=Dict{String, Any}(),
    )
end

"""Get [`KauraPLL`](@ref) `ω_lp`."""
get_ω_lp(value::KauraPLL) = value.ω_lp
"""Get [`KauraPLL`](@ref) `kp_pll`."""
get_kp_pll(value::KauraPLL) = value.kp_pll
"""Get [`KauraPLL`](@ref) `ki_pll`."""
get_ki_pll(value::KauraPLL) = value.ki_pll
"""Get [`KauraPLL`](@ref) `ext`."""
get_ext(value::KauraPLL) = value.ext
"""Get [`KauraPLL`](@ref) `states`."""
get_states(value::KauraPLL) = value.states
"""Get [`KauraPLL`](@ref) `n_states`."""
get_n_states(value::KauraPLL) = value.n_states

"""Set [`KauraPLL`](@ref) `ω_lp`."""
set_ω_lp!(value::KauraPLL, val) = value.ω_lp = val
"""Set [`KauraPLL`](@ref) `kp_pll`."""
set_kp_pll!(value::KauraPLL, val) = value.kp_pll = val
"""Set [`KauraPLL`](@ref) `ki_pll`."""
set_ki_pll!(value::KauraPLL, val) = value.ki_pll = val
"""Set [`KauraPLL`](@ref) `ext`."""
set_ext!(value::KauraPLL, val) = value.ext = val


text:
# # How similar do two strings look?

# ## VisualStringDistances.jl

#

# ---

# ## Let's compare strings

#

using StringDistances

# How many single-character edits are needed to turn "Julia" into "JuIia"?

StringDistances.Levenshtein()("Julia", "JuIia")

# What about "Julia" into "JuQia"?

StringDistances.Levenshtein()("Julia", "JuQia")

# We can also compare based on how many times consecutive pairs of letters appear in each string...

StringDistances.QGram(2)("Julia", "JuIia"), StringDistances.QGram(2)("Julia", "JuQia")

# ---
# ## Visual distances

#

# But none of these take into account that "Julia" and "JuIia" look pretty similar, while "Julia" and "JuQia" look pretty different.

using VisualStringDistances: VisualStringDistances
const VSD = VisualStringDistances
VSD.visual_distance("Julia",  "JuIia"), VSD.visual_distance("Julia", "JuQia")

#

# That seems better! But how do we know it does something reasonable in other cases too? And how does it work?

#

# Just need two tools:
# 1.
 A way to translate strings into images

# 2.
 A way to compare images

# ---

# ## 1. A way to translate strings into images: GNU Unifont

VSD.printglyph("GNU Unifont"; symbols=("#", "-"))

# A bitmap font!

# ---

# Unifont stores characters as bitmaps, making things quite easy for us:

VSD.Glyph("Julia")

#

# (see also FreeTypeAbstraction.jl to render bitmaps from many fonts!)

# ---

# It is low resolution, but simple and comprehensive, with 57086 supported characters, including...

chars = [VSD.get_char(k) for k in rand(collect(keys(VSD.UNIFONT_LOOKUP)), 5)];
permutedims(chars)

# Which render as:

VSD.printglyph(join(chars, " "))

# ---

VSD.printglyph("Julia vs JuIia"); VSD.printglyph("Julia vs JuQia") # hide

# ---

# ## 2. A way to compare images: Optimal transport

#

# *
 you have $a(x_1)$ amount of stuff at site $x_1$, $a(x_2)$ amount of stuff at $x_2$, ..., $a(x_n)$ stuff at $x_n$.

# *
 you want to move it around until you have $b(y_1)$ stuff at site $y_1$, $b(y_2)$ stuff at $y_2$, ..., $b(y_m)$ stuff at $y_m$

# *
 it costs $c(x_i, y_j)$ to move one unit of mass from $x_i$ to $y_j$

# ```math
# \begin{aligned}
# \operatorname{OT}(a,b) := \text{minimize} \quad & \sum\_{x,y} π(x,y)\, c(x,y)\\\\
# \text{such that} \quad & a(x) = \sum\_{y} \pi(x,y)\\\\
# & b(y) = \sum\_{x} \pi(x,y) \\\\
# & \pi(x,y) \geq 0
# \end{aligned}
# ```

# *
 We optimize to find the variables $\pi(x,y)$ (how much stuff to move from $x$ to $y$)

# ---

# ## How does optimal transport relate to our problem?

#

# -
 If we have a black pixel in the 3rd column and 2nd row of the bitmap, we can see that as $a(1) = 1$ unit of mass at site $x_1 = (2,3)$.

# -
 In this way, we can translate the bitmap representation of the string into the language of optimal transport.

# -
 $c(x,y)$ is just the distance between those points

# -
 Note: we do two modifications to this

#   -
 we solve an approximate version for speed ("entropic regularization")

#   -
 add penalties for creating/destroying stuff for the case $\sum_x a(x) \neq  \sum_y b(y)$     [1].

#

#

#

# [1]: Séjourné, T., Feydy, J., Vialard, F.-X., Trouvé, A., Peyré, G., 2019. *Sinkhorn Divergences for Unbalanced Optimal Transport*. https://arxiv.org/abs/1910.12958.

# ---

# ## What use does this have?

# Making gifs!

# ---

# ## What use does this have?

# Adding a check for new packages being added to the General registry to try to prevent the malicious impersonation another package.

# Two main concerns:

# 1. Possibly, one will make a typo, and end up at the wrong package ("typosquatting") $\leadsto$ edit distance check
# 2. Possibly, one will copy a malicious tutorial that has mimicked the appearance of the name of a popular package $\leadsto$ visual distance

#

# ---

# ## Is this the right visual distance for an automated registry check?

# I'm not sure.

# *
 Human perception is actually a bit different

#   * e.g. we mix up "p" vs "q" more than "a" vs "e"  [2], but `visual_distance` says "p" and "q" are further apart than "a" and "e"
# *
 optimal transport is a bit slow (though not prohibitively so, with entropic regularization and the low resolution font)

# *
 there are several parameters and cutoffs to tune

#
#
# Possibly a perceptually-weighted edit distance is more sensible.

#

#

# [2]: Courrieu, Pierre, Fernand Farioli, and Jonathan Grainger. *Inverse Discrimination Time as a Perceptual Distance for Alphabetic Characters*. Visual Cognition 11, no. 7 (October 2004): 901–19. https://doi.org/10.1080/13506280444000049.

# ---

# ## References & Notes

# * Package for `visual_distance`, `printglyph`, etc: VisualStringDistances.jl
# *
 Package with the underlying algorithm optimal transport algorithm: UnbalancedOptimalTransport.jl

#

#

# References:
#

# [1]: Séjourné, T., Feydy, J., Vialard, F.-X., Trouvé, A., Peyré, G., 2019. *Sinkhorn Divergences for Unbalanced Optimal Transport*. https://arxiv.org/abs/1910.12958.
#
# [2]: Courrieu, Pierre, Fernand Farioli, and Jonathan Grainger. *Inverse Discrimination Time as a Perceptual Distance for Alphabetic Characters*. Visual Cognition 11, no. 7 (October 2004): 901–19. https://doi.org/10.1080/13506280444000049.

#

# Slides made with the help of Remark.jl, Literate.jl, and Documenter.jl; gifs made with Makie.jl.

# Thanks to Stefan Karpinski for suggesting GNU Unifont.


text:
module CanonicalizeRequest

import ..Layer, ..request
using ..Messages
using ..Strings: tocameldash

"""
    request(CanonicalizeLayer, method, ::URI, headers, body) -> HTTP.Response

Rewrite request and response headers in Canonical-Camel-Dash-Format.
"""
abstract type CanonicalizeLayer{Next <: Layer} <: Layer{Next} end
export CanonicalizeLayer

function request(::Type{CanonicalizeLayer{Next}},
                 method::String, url, headers, body; kw...) where Next

    headers = canonicalizeheaders(headers)

    res = request(Next, method, url, headers, body; kw...)

    res.headers = canonicalizeheaders(res.headers)

    return res
end

canonicalizeheaders(h::T) where {T} = T([tocameldash(k) => v for (k,v) in h])

end # module CanonicalizeRequest


text:
using Flux

# Build a VGG block
#  ifilters: number of input filters
#  ofilters: number of output filters
#  batchnorm: add batchnorm
function vgg_block(ifilters, ofilters, depth, batchnorm)
  k = (3,3)
  p = (1,1)
  layers = []
  for l in 1:depth
    if batchnorm
      w = convfilter(k, ifilters=>ofilters)
      b = Zeros()
      push!(layers, Conv(weight=w, bias=b, pad=p))
      push!(layers, BatchNorm(ofilters, relu))
    else
      push!(layers, Conv(k, ifilters=>ofilters, relu, pad=p))
    end
    ifilters = ofilters
  end
  return layers
end

# Build convolutionnal layers
#  config: :A (vgg11) :B (vgg13) :D (vgg16) :E (vgg19)
#  inchannels: number of channels in input image (3 for RGB)
function convolutional_layers(config, batchnorm, inchannels)
  layers = []
  ifilters = inchannels
  for c in configs[config]
    push!(layers, vgg_block(ifilters, c..., batchnorm)...)
    push!(layers, MaxPool((2,2)))
    ifilters, _ = c
  end
  return layers
end

# Build classification layers
#  imsize: image size
#  nclasses: number of classes
#  fcsize: size of fully connected layers (usefull for smaller nclasses than ImageNet)
#  dropout: dropout obviously
function classifier_layers(imsize, nclasses, fcsize, dropout)
  layers = []
  push!(layers, flatten)
  push!(layers, Dense(Int(prod(imsize) / 2), fcsize, relu))
  push!(layers, Dropout(dropout))
  push!(layers, Dense(fcsize, fcsize, relu))
  push!(layers, Dropout(dropout))
  push!(layers, Dense(fcsize, nclasses))
  push!(layers, softmax)
  return layers
end

function vgg(imsize; config, batchnorm=false, nclasses)
  conv = convolutional_layers(config, batchnorm, 3)
  class = classifier_layers(imsize, nclasses, 4096, 0.5)
  return Chain(conv..., class...)
end

const configs = Dict(:A => [(64,1) (128,1) (256,2) (512,2) (512,2)],
                     :B => [(64,2) (128,2) (256,2) (512,2) (512,2)],
                     :D => [(64,2) (128,2) (256,3) (512,3) (512,3)],
                     :E => [(64,2) (128,2) (256,4) (512,4) (512,4)])

vgg11(imsize; nclasses) =
  vgg(imsize, config=:A, nclasses=nclasses)

vgg11bn(imsize; nclasses) =
  vgg(imsize, config=:A, batchnorm=true, nclasses=nclasses)

vgg13(imsize; nclasses) =
  vgg(imsize, config=:B, nclasses=nclasses)

vgg13bn(imsize; nclasses) =
  vgg(imsize, config=:B, batchnorm=true, nclasses=nclasses)

vgg16(imsize; nclasses) =
  vgg(imsize, config=:D, nclasses=nclasses)

vgg16bn(imsize; nclasses) =
  vgg(imsize, config=:D, batchnorm=true, nclasses=nclasses)

vgg19(imsize; nclasses) =
  vgg(imsize, config=:E, nclasses=nclasses)

vgg19bn(imsize; nclasses) =
  vgg(imsize, config=:E, batchnorm=true, nclasses=nclasses)


text:
using Cairo, GeometryTypes, Visualize
using Visualize: CairoWindow, AbstractCairoWindow, Canvas, Scene, ColorBuffer, DepthBuffer, clip2pixel_space, CairoContext
using Visualize.JLRasterization: JLRasterizer, edge_function, interpolate,standard_transparency, broadcastmin, broadcastmax
using FileIO, GeometryTypes, Colors
using Visualize: normalmesh, Light, Shading, MeshUniforms, Camera, visualize, Vertices, AbstractCairoWindow, FRect

window = CairoWindow(Camera => BasicCamera())

canvas = window[Canvas]
x = canvas[Visualize.CairoSurface]
y = canvas[Visualize.CairoContext]

resolution = (1024, 1024)

x = LineSegments(Position => rand(Vec2f0, 32), Color => Vec4f0(1, 1, 1, 1))

empty!(window[Visualize.Renderlist])

visualize(window, x)
Visualize.draw_window!(window)
Cairo.write_to_png(canvas[Visualize.CairoSurface], homedir()*"/Desktop/test.png")


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

using Markdown
using InteractiveUtils

# ╔═╡ 5b928827-1226-44f8-a295-0a5d072f3d9e
begin
	using Pkg
	Pkg.add(file="/home/graham_s/julia/vw/ScottishTaxBenefitModel" )

end

# ╔═╡ 6aace8e2-4efc-11ec-1460-d15ea3d01b48
begin
	using SurveyDataWeighting
	using BudgetConstraints
	using PovertyAndInequalityMeasures
	using Plots,DataFrames,CSV
end

# ╔═╡ a0fa5a6f-16c7-4141-af1a-def84e31546e
begin
	using ScottishTaxBenefitModel

end

# ╔═╡ Cell order:
# ╠═6aace8e2-4efc-11ec-1460-d15ea3d01b48
# ╠═5b928827-1226-44f8-a295-0a5d072f3d9e
# ╠═a0fa5a6f-16c7-4141-af1a-def84e31546e


text:
module Gen

using ..ProtoBuf
using ..ProtoBuf.GoogleProtoBuf
using ..ProtoBuf.GoogleProtoBufCompiler
using Logging

import ProtoBuf: meta, DEF_REQ, DEF_FNUM, DEF_VAL, DEF_PACK

export gen

const JTYPES              = [Float64, Float32, Int64, UInt64, Int32, UInt64,  UInt32,  Bool,  AbstractString, Any,     Any,     Vector{UInt8}, UInt32, Int32, Int32, Int64, Int32, Int64]
const JTYPE_DEFAULTS      = [0,       0,       0,     0,      0,     0,       0,       false, "",             nothing, nothing, UInt8[],       0,      0,     0,     0,     0,     0]

isprimitive(fldtype) = (1 <= fldtype <= 8) || (13 <= fldtype <= 18)

# maps protofile name to package name
const _packages = Dict{String,String}()

# maps protofile name to array of imported protofiles
const protofile_imports = Dict()

# maps resolved names in the "protofile module name" -> "type name" -> ["nested type name",...] form
const name_maps = Dict()

function add_name_map(name, fullname)
    parts = rsplit(fullname, "."; limit=2, keepempty=false)

    if length(parts) == 1
        pkgname = ""
        typname = parts[1]
    else
        pkgname = parts[1]
        typname = parts[2]
    end

    if name == typname
        nested_level = ""
    else
        nested_level = typname[1:end-length(name)-1]
        typname = name
    end

    names_in_module = get!(name_maps, pkgname) do
        Dict{String,Set{String}}()
    end
    nested_type_names = get!(names_in_module, typname) do
        Set{String}()
    end
    push!(nested_type_names, nested_level)
    nothing
end

# Treat Google Proto3 extensions specially as they are built into ProtoBuf.jl (for issue #77)
const GOOGLE_PROTO3_EXTENSIONS = "google.protobuf"

_module_postfix = false
_map_as_array = false

#--------------------------------------------------------------------
# begin keyword handling
#--------------------------------------------------------------------
const _keywords = [
    "if", "else", "elseif", "while", "for", "begin", "end", "quote",
    "try", "catch", "return", "local", "abstract", "function", "macro",
    "ccall", "finally", "typealias", "break", "continue", "type",
    "global", "module", "using", "import", "export", "const", "let",
    "bitstype", "do", "baremodule", "importall", "immutable",
    "Type", "Enum", "Any", "DataType", "Base", "Set"
]

chk_keyword(name) = (name in _keywords) ? string('_', name) : String(name)
#--------------------------------------------------------------------
# end keyword handling
#--------------------------------------------------------------------

#--------------------------------------------------------------------
# begin name resolution utilities
#--------------------------------------------------------------------

mutable struct Scope
    name::String
    syms::Vector{String}
    files::Vector{String}
    is_module::Bool
    children::Vector{Scope}
    parent::Scope

    Scope(name::String; is_module::Bool=false) = new(name, String[], String[], is_module, Scope[])
    function Scope(name::String, parent::Scope; is_module::Bool=false)
        s = new(name, String[], String[], is_module, Scope[], parent)
        push!(parent.children, s)
        s
    end
end

const top_scope = Scope("")

"""
Full Julia name of a scope, that can be resolved from Julia code.
Scope can be a module or a struct within a module.

Modules are separated with `.` and structs with `_`.

When an additional name is passed, it is prefixed with the scope according to the same rules.
"""
function fullname(s::Scope)
    if isempty(s.name) || isempty(s.parent.name)
        return s.name
    end
    string(fullname(s.parent), s.parent.is_module ? "." : "_", s.name)
end
function fullname(scope::Scope, name::String)
    if isempty(scope.name)
        name
    else
        sep = scope.is_module ? "." : "_"
        string(fullname(scope), sep, name)
    end
end

scope_has_file(s::Scope, f::String) = (f in s.files || any(x->scope_has_file(x,f), s.children))

"""
Given the module name, return the corresponding scope object.
Inserts a new scope if required.
"""
function get_module_scope(parent::Scope, module_name::String)
    if _module_postfix
        module_name = module_name * "_pb"
    end
    # return existing scope if possible
    for s in parent.children
        if s.name == module_name
            return s
        end
    end
    # create new scope otherwise
    Scope(module_name, parent; is_module=true)
end

"""
Search for `name` in the scope hierarchy and qualify it with full scope.
"""
function qualify_in_hierarchy(name::String, scope::Scope)
    if startswith(name, "ProtoBuf.google.protobuf.")
        return name
    elseif name in scope.syms
        return fullname(scope, name)
    elseif isdefined(scope, :parent)
        return qualify_in_hierarchy(name, scope.parent)
    else
        error("unresolved name $name at scope $(scope.name)")
    end
end

splitmodule(name::String) = rsplit(name, '.'; limit=2)
function splitmodule_chkkeyword(name::String)
    sm = splitmodule(name)
    if length(sm) > 1
        modul = String(sm[1])
        name = String(sm[2])
    else
        modul = ""
    end
    name = chk_keyword(name)
    modul,name
end

function findmodule(name::String)
    mlen = 0
    mpkg = ""
    for pkg in values(_packages)
        orig = _module_postfix ? join(map(m->m[1:end-3], split(pkg, ".")), ".") : pkg
        if (length(pkg) > mlen) && startswith(name, orig)
            mlen = length(orig)
            mpkg = pkg
        end
    end
    (mpkg, replace((0 == mlen) ? name : name[(mlen+2):end], '.'=>'_'))
end

function field_type_name(full_type_name::String)
    comps = split(full_type_name, '.'; keepempty=false)
    if isempty(comps)
        type_name = full_type_name
    else
        type_name = join(comps, '.')
        for level in (length(comps)-1):-1:1
            package_name = join(comps[1:level], '.')
            if package_name == GOOGLE_PROTO3_EXTENSIONS
                comps[end] = chk_keyword(comps[end])
                type_name = "ProtoBuf.$(join(comps, '.'))"
                break
            elseif package_name in keys(name_maps)
                type_maps = name_maps[package_name]
                if last(comps) in keys(type_maps)
                    nested_level = join(comps[(level+1):(end-1)], '_')
                    if nested_level in type_maps[last(comps)]
                        type_name = isempty(nested_level) ? last(comps) : string(nested_level, "_", last(comps))
                        isempty(package_name) || (type_name = string(package_name, ".", type_name))
                        break
                    end
                end
            end
        end
    end
    @debug("usable type name for $full_type_name is $type_name")
    type_name
end

#--------------------------------------------------------------------
# end name resolution utilities
#--------------------------------------------------------------------

#--------------------------------------------------------------------
# begin deferred writing (correct ordering) of generated entities
#--------------------------------------------------------------------
struct DeferredWrite
    iob::IOBuffer
    depends::Vector{String}
end
const _deferred = Dict{String,DeferredWrite}()
# set of fully-qualified names we've resolved, to handle dependencies in other files
const _all_resolved = Set{String}()

function defer(name::String, iob::IOBuffer, depends::String)
    @debug("defer $name due to $depends")
    if isdeferred(name)
        depsnow = _deferred[name].depends
        !(depends in depsnow) && push!(depsnow, depends)
        return
    end
    _deferred[name] = DeferredWrite(iob, String[depends])
    nothing
end

isdeferred(name::String) = haskey(_deferred, name)
function isresolved(dtypename::String, referenced_name::String, full_referenced_name::String, exports::Vector{String})
    (dtypename == referenced_name) && return true
    for jtype in JTYPES
        (referenced_name == string(jtype)) && return true
    end
    if '.' in referenced_name
        return !isdeferred(referenced_name)
    elseif referenced_name in exports
        return !(isdeferred(referenced_name) || isdeferred(full_referenced_name))
    end
    false
end

function resolve(iob::IOBuffer, name::String)
    fully_resolved = String[]
    for (typ,dw) in _deferred
        idx = something(findfirst(isequal(name), dw.depends), 0)
        (idx == 0) && continue
        splice!(dw.depends, idx)
        isempty(dw.depends) && push!(fully_resolved, typ)
    end

    # write all fully resolved entities
    for typ in fully_resolved
        @debug("fully resolved $typ")
        print(iob, String(take!(_deferred[typ].iob)))
        delete!(_deferred, typ)
        push!(_all_resolved, typ)
    end

    # mark them resolved as well
    for typ in fully_resolved
        resolve(iob, typ)
    end
end
#--------------------------------------------------------------------
# end deferred writing (correct ordering) of generated entities
#--------------------------------------------------------------------

#--------------------------------------------------------------------
# start utility methods
#--------------------------------------------------------------------
function protofile_name_to_module_name(n::String)
    name = splitext(basename(n))[1]
    name = replace(name, '.'=>'_')
    name = string(name, "_pb")
    return name
end

function append_response(resp::CodeGeneratorResponse, protofile::FileDescriptorProto, io::IOBuffer)
    outdir = dirname(protofile.name)
    filename = protofile_name_to_module_name(protofile.name)
    filename = string(filename, ".jl")

    append_response(resp, filename, io)
end

function append_response(resp::CodeGeneratorResponse, filename::String, io::IOBuffer)
    jfile = ProtoBuf.instantiate(CodeGeneratorResponse_File)

    jfile.name = filename
    jfile.content = String(take!(io))
    !hasproperty(resp, :file) && (resp.file = CodeGeneratorResponse_File[])
    push!(resp.file, jfile)
    resp
end

function err_response(errio::IOBuffer)
    resp = ProtoBuf.instantiate(CodeGeneratorResponse)
    resp.error = String(take!(errio))
    resp
end
#--------------------------------------------------------------------
# end utility methods
#--------------------------------------------------------------------

#--------------------------------------------------------------------
# begin code generation
#--------------------------------------------------------------------
function generate_enum(io::IO, errio::IO, enumtype::EnumDescriptorProto, scope::Scope, exports::Vector{String})
    _modul,enumname = splitmodule_chkkeyword(fullname(scope, enumtype.name))
    push!(scope.syms, enumname)

    @debug("begin enum $(enumname)")
    println(io, "const ", enumname, " = (;[")
    for value::EnumValueDescriptorProto in enumtype.value
        # If we find that the field name is a keyword prepend it with `_`
        fldname = chk_keyword(value.name)
        println(io, "    Symbol(\"", fldname, "\") => Int32(", value.number, "),")
    end
    println(io, "]...)")
    println(io, "")
    push!(exports, enumname)
    @debug("end enum $(enumname)")
    nothing
end

function generate_msgtype(outio::IO, errio::IO, dtype::DescriptorProto, scope::Scope, syntax::String, exports::Vector{String}, depends::Vector{String}, mapentries::Dict{String,Tuple{String,String}}, deferedmode::Bool)
    full_dtypename = fullname(scope, dtype.name)

    io = IOBuffer()
    modul,dtypename = splitmodule_chkkeyword(full_dtypename)
    full_dtypename = (modul=="") ? dtypename : "$(modul).$(dtypename)"
    deferedmode && !isdeferred(full_dtypename) && return
    @debug("begin type $(full_dtypename)")
    add_name_map(dtype.name, full_dtypename)

    scope = Scope(dtype.name, scope)

    # check oneof
    oneof_names = Vector{String}()
    if hasproperty(dtype, :oneof_decl)
        for oneof_decl in dtype.oneof_decl
            if hasproperty(oneof_decl, :name)
                push!(oneof_names,  "Symbol(\"$(oneof_decl.name)\")")
            end
        end
    end

    # generate enums
    if hasproperty(dtype, :enum_type)
        for enum_type in dtype.enum_type
            generate_enum(io, errio, enum_type, scope, exports)
            (errio.size > 0) && return
        end
    end

    # generate nested types
    if hasproperty(dtype, :nested_type)
        for nested_type::DescriptorProto in dtype.nested_type
            generate_msgtype(io, errio, nested_type, scope, syntax, exports, depends, mapentries, deferedmode)
            (errio.size > 0) && return
        end
    end

    # generate this type
    reqflds = String[]
    allflds = String[]
    getpropertylist = Pair[]
    packedflds = String[]
    fldnums = Int[]
    defvals = String[]
    wtypes = String[]
    realfldtypes = String[]
    oneofs = isempty(oneof_names) ? Int[] : zeros(Int, length(dtype.field))
    ismapentry = hasproperty(dtype, :options) && hasproperty(dtype.options, :map_entry) && dtype.options.map_entry
    map_key_type = ""
    map_val_type = ""

    if hasproperty(dtype, :field)
        for fld_idx in 1:length(dtype.field)
            field = dtype.field[fld_idx]
            if hasproperty(field, :oneof_index) && !isempty(oneofs) && !ProtoBuf.isdefaultproperty(field, :oneof_index)
                oneof_idx = field.oneof_index + 1
                oneofs[fld_idx] = oneof_idx
            end

            # If we find that the field name is a keyword prepend it with _type
            fldname = chk_keyword(field.name)
            if field._type == FieldDescriptorProto_Type.TYPE_GROUP
                println(errio, "Groups are not supported")
                return
            end

            full_typ_name = ""
            if (field._type == FieldDescriptorProto_Type.TYPE_MESSAGE) || (field._type == FieldDescriptorProto_Type.TYPE_ENUM)
                typ_name = field.type_name
                if startswith(typ_name, '.')
                    (m,t) = findmodule(typ_name[2:end])
                    t = chk_keyword(t)
                    full_typ_name = m=="" ? t : "$(m).$(t)"
                    typ_name = (m == modul) ? t : full_typ_name
                else
                    full_typ_name = qualify_in_hierarchy(typ_name, scope)
                    typ_name = full_typ_name
                    m,t = splitmodule_chkkeyword(typ_name)
                    (m == modul) && (typ_name = t)
                    full_typ_name = m=="" ? t : "$(m).$(t)"
                end
            elseif field._type == FieldDescriptorProto_Type.TYPE_SINT32
                push!(wtypes, ":$fldname => :sint32")
            elseif field._type == FieldDescriptorProto_Type.TYPE_SINT64
                push!(wtypes, ":$fldname => :sint64")
            elseif field._type == FieldDescriptorProto_Type.TYPE_FIXED32
                push!(wtypes, ":$fldname => :fixed32")
            elseif field._type == FieldDescriptorProto_Type.TYPE_SFIXED32
                push!(wtypes, ":$fldname => :sfixed32")
            elseif field._type == FieldDescriptorProto_Type.TYPE_FIXED64
                push!(wtypes, ":$fldname => :fixed64")
            elseif field._type == FieldDescriptorProto_Type.TYPE_SFIXED64
                push!(wtypes, ":$fldname => :sfixed64")
            end
            enum_typ_name = (field._type == FieldDescriptorProto_Type.TYPE_ENUM) ? typ_name : ""
            (field._type != FieldDescriptorProto_Type.TYPE_MESSAGE) && (typ_name = "$(JTYPES[field._type])")

            push!(fldnums, field.number)
            (FieldDescriptorProto_Label.LABEL_REQUIRED == field.label) && push!(reqflds, ":"*fldname)

            if hasproperty(field, :default_value) && !isempty(field.default_value)
                if field._type == FieldDescriptorProto_Type.TYPE_STRING
                    push!(defvals, ":$fldname => \"$(escape_string(field.default_value))\"")
                elseif field._type == FieldDescriptorProto_Type.TYPE_MESSAGE
                    println(errio, "Default values for message types are not supported. Field: $(dtypename).$(fldname) has default value [$(field.default_value)]")
                    return
                elseif field._type == FieldDescriptorProto_Type.TYPE_BYTES
                    println(errio, "Default values for byte array types are not supported. Field: $(dtypename).$(fldname) has default value [$(field.default_value)]")
                    return
                elseif field._type == FieldDescriptorProto_Type.TYPE_ENUM
                    push!(defvals, ":$fldname => $(field_type_name(enum_typ_name)).$(field.default_value)")
                elseif field._type == FieldDescriptorProto_Type.TYPE_DOUBLE
                    defval = replace(replace(field.default_value, "inf" => "Inf"), "nan" => "NaN")
                    push!(defvals, ":$fldname => $defval")
                elseif field._type == FieldDescriptorProto_Type.TYPE_FLOAT
                    defval = replace(replace(field.default_value, "inf" => "Inf32"), "nan" => "NaN32")
                    push!(defvals, ":$fldname => $defval")
                else
                    push!(defvals, ":$fldname => $(field.default_value)")
                end
            end

            packed = (hasproperty(field, :options) && field.options.packed) ||
                     ((syntax == "proto3") && (FieldDescriptorProto_Label.LABEL_REPEATED == field.label) && isprimitive(field._type))
            packed && push!(packedflds, ":"*fldname)

            gen_typ_name = ""
            if !(isresolved(dtypename, typ_name, full_typ_name, exports) || full_typ_name in _all_resolved)
                if deferedmode
                    gen_typ_name = "Base.Any"
                else
                    defer(full_dtypename, io, full_typ_name)
                end
            end

            typ_name = field_type_name(typ_name)
            is_typ_mapentry = typ_name in keys(mapentries)
            if is_typ_mapentry && !_map_as_array
                k,v = mapentries[typ_name]
                typ_name = "Base.Dict{$k,$v}"
                if deferedmode
                    # because we do not know if Dict key and value types are resolved yet
                    gen_typ_name = "Base.Dict"
                end
            elseif FieldDescriptorProto_Label.LABEL_REPEATED == field.label
                typ_name = "Base.Vector{$typ_name}"
            end

            if isempty(gen_typ_name)
                gen_typ_name = typ_name
                push!(allflds, ":$fldname => $(typ_name)")
                push!(getpropertylist, ":$fldname" => "$(typ_name)")
            else
                push!(realfldtypes, ":$fldname => \"$(typ_name)\"")
                push!(allflds, ":$fldname => \"$(typ_name)\"")
                push!(getpropertylist, ":$fldname" => "Any")
            end

            if ismapentry
                if field.number == 1
                    map_key_type = typ_name
                elseif field.number == 2
                    map_val_type = typ_name
                end
            end
        end
    end

    # generate struct body
    println(io, """
        mutable struct $(dtypename) <: ProtoType
            __protobuf_jl_internal_meta::ProtoMeta
            __protobuf_jl_internal_values::Dict{Symbol,Any}
            __protobuf_jl_internal_defaultset::Set{Symbol}

            function $(dtypename)(; kwargs...)
                obj = new(meta($(dtypename)), Dict{Symbol,Any}(), Set{Symbol}())
                values = obj.__protobuf_jl_internal_values
                symdict = obj.__protobuf_jl_internal_meta.symdict
                for nv in kwargs
                    fldname, fldval = nv
                    fldtype = symdict[fldname].jtyp
                    (fldname in keys(symdict)) || error(string(typeof(obj), " has no field with name ", fldname))
                    if fldval !== nothing
                        values[fldname] = isa(fldval, fldtype) ? fldval : convert(fldtype, fldval)
                    end
                end
                obj
            end""")
    println(io, "end # mutable struct $(dtypename)", ismapentry ? " (mapentry)" : "", deferedmode ? " (has cyclic type dependency)" : "")

    # generate the meta for this type
    @debug("generating meta", dtypename)
    _d_fldnums = [1:length(fldnums);]
    println(io, """const __meta_$(dtypename) = Ref{ProtoMeta}()
        function meta(::Type{$dtypename})
            ProtoBuf.metalock() do
                if !isassigned(__meta_$dtypename)
                    __meta_$(dtypename)[] = target = ProtoMeta($dtypename)""")
    !isempty(reqflds) && println(io, "            req = Symbol[$(join(reqflds, ','))]")
    !isempty(defvals) && println(io, "            val = Dict{Symbol,Any}($(join(defvals, ", ")))")
    (fldnums != _d_fldnums) && println(io, "            fnum = Int[$(join(fldnums, ','))]")
    !isempty(packedflds) && println(io, "            pack = Symbol[$(join(packedflds, ','))]")
    !isempty(wtypes) && println(io, "            wtype = Dict($(join(wtypes, ", ")))")
    println(io, "            allflds = Pair{Symbol,Union{Type,String}}[$(join(allflds, ", "))]")
    if !isempty(oneofs)
        println(io, "            oneofs = Int[$(join(oneofs, ','))]")
        println(io, "            oneof_names = Symbol[$(join(oneof_names, ','))]")
    end
    print(io, "            meta(target, $(dtypename), allflds, ")
    print(io, isempty(reqflds) ? "ProtoBuf.DEF_REQ, " : "req, ")
    print(io, (fldnums == _d_fldnums) ? "ProtoBuf.DEF_FNUM, " : "fnum, ")
    print(io, isempty(defvals) ? "ProtoBuf.DEF_VAL, " : "val, ")
    print(io, isempty(packedflds) ? "ProtoBuf.DEF_PACK, " : "pack, ")
    print(io, isempty(wtypes) ? "ProtoBuf.DEF_WTYPES, " : "wtype, ")
    print(io, isempty(oneofs) ? "ProtoBuf.DEF_ONEOFS, " : "oneofs, ")
    print(io, isempty(oneofs) ? "ProtoBuf.DEF_ONEOF_NAMES" : "oneof_names")
    println(io, ")")
    println(io, "        end")
    println(io, "        __meta_$(dtypename)[]")
    println(io, "    end")
    println(io, "end")

    # generate new getproperty method
    if !isempty(getpropertylist)
        @debug("generating getproperty", dtypename)
        println(io, "function Base.getproperty(obj::$(dtypename), name::Symbol)")
        pfx = "if"
        for prop in getpropertylist
            println(io, "    $(pfx) name === $(first(prop))")
            println(io, "        return (obj.__protobuf_jl_internal_values[name])::$(last(prop))")
            if pfx == "if"
                pfx = "elseif"
            end
        end
        println(io, "    else")
        println(io, "        getfield(obj, name)")
        println(io, "    end")
        println(io, "end")
    end

    println(io, "")
    push!(exports, dtypename)
    ismapentry && (mapentries[dtypename] = (map_key_type, map_val_type))

    deferedmode && (full_dtypename in keys(_deferred)) && delete!(_deferred, full_dtypename)

    if !isdeferred(full_dtypename)
        @debug("resolved (!deferred)", full_dtypename)
        print(outio, String(take!(io)))
        deferedmode || resolve(outio, full_dtypename)
        push!(_all_resolved, full_dtypename)
    end

    @debug("end type", full_dtypename)
    nothing
end

function scoped_svc_type_name(scope::Scope, typ_name)
    modul = fullname(scope)

    if isempty(modul)
        return typ_name
    elseif startswith(typ_name, string(modul, "."))
        return typ_name[(length(modul)+2):end]
    elseif startswith(typ_name, '.')
        (m,t) = findmodule(typ_name[2:end])
        t = chk_keyword(t)
        full_typ_name = m=="" ? t : "$(m).$(t)"
        return (m == modul) ? t : full_typ_name
    else
        full_typ_name = qualify_in_hierarchy(typ_name, scope)
        typ_name = full_typ_name
        m,t = splitmodule_chkkeyword(typ_name)
        (m == modul) && (typ_name = t)
        return typ_name
    end
end

function generate_svc(io::IO, errio::IO, stype::ServiceDescriptorProto, scope::Scope, svcidx::Int, exports::Vector{String})
    nmethods = hasproperty(stype, :method) ? length(stype.method) : 0

    # generate method and service descriptors
    println(io, "# service methods for $(stype.name)")
    println(io, "const _$(stype.name)_methods = MethodDescriptor[")
    for idx in 1:nmethods
        method = stype.method[idx]
        in_typ_name = scoped_svc_type_name(scope, field_type_name(method.input_type))
        out_typ_name = scoped_svc_type_name(scope, field_type_name(method.output_type))
        elem_sep = (idx < nmethods) ? "," : ""

        method.client_streaming && (in_typ_name = "Channel{" * in_typ_name * "}")
        method.server_streaming && (out_typ_name = "Channel{" * out_typ_name * "}")

        println(io, "        MethodDescriptor(\"$(method.name)\", $(idx), $(in_typ_name), $(out_typ_name))$(elem_sep)")
    end
    println(io, "    ] # const _$(stype.name)_methods")
    fullservicename = scope.is_module ? fullname(scope, stype.name) : stype.name
    println(io, "const _$(stype.name)_desc = ServiceDescriptor(\"$(fullservicename)\", $(svcidx), _$(stype.name)_methods)")
    println(io, "")

    # generate service
    println(io, "$(stype.name)(impl::Module) = ProtoService(_$(stype.name)_desc, impl)")
    println(io, "")
    push!(exports, stype.name)

    # generate stubs
    stub = "$(stype.name)Stub"
    println(io, "mutable struct $(stub) <: AbstractProtoServiceStub{false}")
    println(io, "    impl::ProtoServiceStub")
    println(io, "    $(stub)(channel::ProtoRpcChannel) = new(ProtoServiceStub(_$(stype.name)_desc, channel))")
    println(io, "end # mutable struct $(stub)")
    println(io, "")
    push!(exports, stub)

    nbstub = "$(stype.name)BlockingStub"
    println(io, "mutable struct $(nbstub) <: AbstractProtoServiceStub{true}")
    println(io, "    impl::ProtoServiceBlockingStub")
    println(io, "    $(nbstub)(channel::ProtoRpcChannel) = new(ProtoServiceBlockingStub(_$(stype.name)_desc, channel))")
    println(io, "end # mutable struct $(nbstub)")
    println(io, "")
    push!(exports, nbstub)

    for idx in 1:nmethods
        method = stype.method[idx]
        in_typ_name = scoped_svc_type_name(scope, field_type_name(method.input_type))
        method.client_streaming && (in_typ_name = "Channel{" * in_typ_name * "}")
        println(io, "$(method.name)(stub::$(stub), controller::ProtoRpcController, inp::$(in_typ_name), done::Function) = call_method(stub.impl, _$(stype.name)_methods[$(idx)], controller, inp, done)")
        println(io, "$(method.name)(stub::$(nbstub), controller::ProtoRpcController, inp::$(in_typ_name)) = call_method(stub.impl, _$(stype.name)_methods[$(idx)], controller, inp)")
        println(io, "")
        push!(exports, method.name)
    end
    nothing
end

function generate_file(io::IO, errio::IO, protofile::FileDescriptorProto)
    @debug("generate begin for $(protofile.name), package $(protofile.package)")

    scope = top_scope
    if !isempty(protofile.package)
        for pkgname in split(protofile.package, '.')
            scope = get_module_scope(scope, String(pkgname))
        end
    end

    push!(scope.files, protofile.name)
    @debug("generated scope for $(protofile.name), package $(protofile.package)")

    # generate module begin
    if !isempty(scope.name)
        scope.is_module = true
        _packages[protofile.name] = fullname(scope)
    end

    # generate syntax version
    syntax = (hasproperty(protofile, :syntax) && !isempty(protofile.syntax)) ? protofile.syntax : "proto2"
    println(io, "# syntax: $(syntax)")

    depends = Vector{String}()
    @debug("generating imports")
    # generate imports
    println(io, "using ProtoBuf")
    dep_imports = Vector{String}()
    add_import = (imp) -> begin
        (imp in dep_imports) || push!(dep_imports, imp)
        nothing
    end
    add_import("ProtoBuf.meta")
    if hasproperty(protofile, :dependency)
        protofile_imports[protofile.name] = protofile.dependency
        using_pkgs = Set{String}()
        for dependency in protofile.dependency
            if haskey(_packages, dependency)
                push!(using_pkgs, _packages[dependency])
                push!(depends, _packages[dependency])
            end
        end

        fullscopename = scope.is_module ? fullname(scope) : ""
        parentscope = (isdefined(scope, :parent) && scope.parent.is_module) ? fullname(scope.parent) : ""
        for dependency in using_pkgs
            # current scope is available by default
            (fullscopename == dependency) && continue
            # google extenstions are included with ProtoBuf
            if dependency == GOOGLE_PROTO3_EXTENSIONS
                dependency = "ProtoBuf." * dependency
                add_import(dependency)
            else
                comps = split(dependency, '.'; keepempty=false)
                if startswith(dependency, parentscope*".")
                    comps[1] = ".." * comps[1]
                elseif !isempty(fullscopename)
                    comps[1] = "._ProtoBuf_Top_." * comps[1]
                end
                add_import(comps[1])
                #for idx in 1:length(comps)
                #    add_import(join(comps[1:idx], '.'))
                #end
            end
        end
    end
    for imp in dep_imports
        println(io, "import ", imp)
    end
    println(io, "")

    exports = Vector{String}()
    mapentries = Dict{String, Tuple{String,String}}()

    # generate top level enums
    @debug("generating enums")
    if hasproperty(protofile, :enum_type)
        for enum_type in protofile.enum_type
            generate_enum(io, errio, enum_type, scope, exports)
            (errio.size > 0) && return
        end
    end

    # generate message types
    @debug("generating types")
    if hasproperty(protofile, :message_type)
        for message_type in protofile.message_type
            generate_msgtype(io, errio, message_type, scope, syntax, exports, depends, mapentries, false)
            (errio.size > 0) && return
        end
    end

    # generate service stubs
    @debug("generating services")
    if hasproperty(protofile, :service)
        nservices = length(protofile.service)
        for idx in 1:nservices
            service = protofile.service[idx]
            generate_svc(io, errio, service, scope, idx, exports)
            (errio.size > 0) && return
        end
    end

    # generate deferred message types
    @debug("generating deferred types")
    if hasproperty(protofile, :message_type)
        for message_type in protofile.message_type
            generate_msgtype(io, errio, message_type, scope, syntax, exports, depends, mapentries, true)
            (errio.size > 0) && return
        end
    end

    # check if everything is resolved
    if !isempty(_deferred)
        println(errio, "All types could not be fully resolved:")
        for (k,v) in _deferred
            println(errio, "$k depends on $(join(v.depends, ','))")
        end
        return
    end

    # generate exports
    !isempty(exports) && println(io, "export ", join(exports, ", "))

    # mention mapentries
    !isempty(mapentries) && println(io, "# mapentries: ", join(mapentries, ", "))

    @debug("generate end for $(protofile.name)")
    nothing
end

function print_package(io::IO, s::Scope, indent="")
    s.is_module || return
    println(io, "$(indent)module $(s.name)")
    nested = indent*"  "
    println(io, "$(nested)const _ProtoBuf_Top_ = @static isdefined(parentmodule(@__MODULE__), :_ProtoBuf_Top_) ? (parentmodule(@__MODULE__))._ProtoBuf_Top_ : parentmodule(@__MODULE__)")
    children = Set(s.children)
    for f in s.files
        if haskey(protofile_imports,f)
            deps = protofile_imports[f]
            # if the file we're about to `include` depends on any of our child modules,
            # generate those first.
            for d in deps
                for c in filter(c->scope_has_file(c,d), children)
                    delete!(children, c)
                    print_package(io, c, nested)
                end
            end
        end

        fname = string(protofile_name_to_module_name(f), ".jl")
        println(io, "$(nested)include(\"$fname\")")
    end
    for c in s.children
        # check if already included
        (c in children) || continue
        print_package(io, c, nested)
    end
    println(io, "$(indent)end")
end

function readreq(srcio::IO)
    req = ProtoBuf.instantiate(CodeGeneratorRequest)
    readproto(srcio, req)
    req
end

function codegen(srcio::IO)
    errio = IOBuffer()
    resp = ProtoBuf.instantiate(CodeGeneratorResponse)
    @debug("generate begin")
    while !eof(srcio)
        req = readreq(srcio)

        if !hasproperty(req, :file_to_generate)
            @debug("no files to generate!!")
            continue
        end

        @debug("generate request for $(length(req.file_to_generate)) proto files")
        @debug("$(req.file_to_generate)")

        #hasproperty(req, :parameter) && @debug("parameter $(req.parameter)")

        for protofile in req.proto_file
            io = IOBuffer()
            generate_file(io, errio, protofile)
            (errio.size > 0) && return err_response(errio)
            append_response(resp, protofile, io)
        end
    end
    if !isempty(top_scope.children)
        for pkg in top_scope.children
            if pkg.is_module
                io = IOBuffer()
                print_package(io, pkg)
                append_response(resp, "$(pkg.name).jl", io)
            end
        end
    end
    @debug("generate end")
    resp
end
#--------------------------------------------------------------------
# end code generation
#--------------------------------------------------------------------

function debug_log_filename(filename=nothing)
    for arg in ARGS
        if startswith(arg, "--debug-gen")
            parts = split(arg, '=')
            filename = (length(parts) == 2) ? strip(last(parts)) : "juliaprotoc.log"
            break
        end
    end
    filename
end

function with_debug(f)
    debug_log_file = debug_log_filename()
    if debug_log_file === nothing
        f()
    else
        open(debug_log_file, "a") do logstream
            with_logger(SimpleLogger(logstream, Logging.Debug)) do
                f()
            end
        end
    end
end

##
# the main read - write method
function gen()
    global _module_postfix = in("--module-postfix-enabled", ARGS)
    global _map_as_array = in("--map-as-array", ARGS)
    with_debug() do
        try
            writeproto(stdout, codegen(stdin))
        catch ex
            println(stderr, "Exception while generating Julia code")
            rethrow()
        end
    end
end

end # module Gen


text:
function bfgs_Stop(nlp :: AbstractNLPModel;
                   x :: Vector{T}=copy(nlp.meta.x0),
                   stp :: NLPStopping = NLPStopping(nlp,
                                                    NLPAtX(nlp.meta.x0)),
                   scaling :: Bool = true,
                   B₀ :: Union{AbstractLinearOperator,
                          AbstractMatrix,
                          UniformScaling{T},
                          Nothing}              = nothing
                   ) where T

    @info log_header([:iter, :f, :dual, :step, :slope], [Int, T, T, T, T],
                     hdr_override=Dict(:f=>"f(x)", :dual=>"‖∇f‖", :slope=>"∇fᵀd"))
    f = obj(nlp,x)
    ∇f = grad(nlp, x)

    xt = similar(x)
    ∇ft = similar(∇f)

    B = one(eltype(x))*I
    if B₀ != nothing
        B = B₀
    end

    τ₀ = 0.0005
    τ₁ = 0.9999

    OK = update_and_start!(stp, x = x, fx = f, gx = ∇f)
    @info log_row(Any[0, f, norm(∇f)])

    while !OK
        d = - B*∇f

        #------------------------------------------
        # Hard coded line search
        hp0 = ∇f'*d
        t=1.0
        # Simple Wolfe forward tracking
        xt = x + t*d
        ft = obj(nlp, xt)
        ∇ft = grad(nlp,xt)
        hp = ∇ft'*d
        #  while  ~wolfe & armijo
        nbW = 0
        while (hp <= τ₁ * hp0) && (ft <= ( f + τ₀*t*hp0)) && (nbW < 10)
            t *= 5
            xt = x + t*d
            ∇ft = grad(nlp,xt)
            hp = ∇ft'*d
            ft = obj(nlp, xt)
            nbW += 1
            @debug "W", ft
        end
        tw = t

        # Simple Armijo backtracking
        nbk = 0
        while ft > ( f + τ₀*t*hp0) && (nbk < 20)
            t *= 0.4
            xt = x + t*d
            ft = obj(nlp, xt)
            nbk += 1
            @debug "A", ft
        end
        #------------------------------------------

        if t!=tw   ∇ft = grad(nlp, xt) end

        # Update BFGS approximation.
        B = push!(B, t * d, ∇ft - ∇f)
        #B = push!(B, t * d, ∇ft - ∇f, scaling)

        #move on
        x .= xt
        f = ft
        ∇f .= ∇ft

        OK = update_and_stop!(stp, x = x, gx = ∇f, fx=f)

        @info log_row(Any[stp.meta.nb_of_stop, f, norm(∇f), t, hp0])
    end

    if !stp.meta.optimal
        #@warn "Optimalité non atteinte"
        @warn status(stp,list=true)
    end

    return stp, B
    #return f, norm(∇f), stp, B
end

""" L-BFGS wrapper of the gereral BFGS implementation.
"""
function L_bfgs_Stop(nlp     :: AbstractNLPModel;
                     x       :: Vector{T}=copy(nlp.meta.x0),
                     stp :: NLPStopping = NLPStopping(nlp,
                                                      NLPAtX(nlp.meta.x0)),
                     mem     :: Int = 5,
                     scaling :: Bool = true,
                     kwargs...
                     ) where T

    @debug "U_Solver = L_bfgs"
    n = nlp.meta.nvar
    B₀ =  InverseLBFGSOperator(Float64, n, mem=mem, scaling=scaling)

    return bfgs_Stop(nlp; x=x, stp = stp, scaling=scaling, B₀=B₀, kwargs...)
end

max_stars_repo_name: tmigot/LSDescent

max_issues_repo_name: tmigot/LSDescent

max_forks_repo_name: tmigot/LSDescent


text:
using Test
using WebSockets
import Dates.now

"Time since start"
tss() = " $(Int(round((now() - T0).value / 1000))) s "
"Server side 'websocket handler'"
function coroutine(ws)
    sttim = time_ns()
    push!(WSLIST, sttim => ws)
    data, success = readguarded(ws)
    s = ""
    if success
        s = String(data)
        @wslog "This websocket is called ", s, " at ", tss()
        push!(WSIDLIST, sttim => s)
    else
        @wslog "Received a closing message instead of a first message", tss()
        return
    end
    try
        data = read(ws)
    catch err
        @wslog s, " was closed at ", tss(), " with message \n\t", err
        push!(CLOSINGMESSAGES, sttim => string(err))
        return
    end
    @wslog "Unexpectedly received a second message from websocket ", s, ". Now closing."
end

function gatekeeper(req, ws)
    #@info "\nOrigin: $(WebSockets.origin(req))  subprotocol: $(subprotocol(req))"
    coroutine(ws)
end

handle(req) = read("limited life websockets.html") |> WebSockets.HTTP.Response

"Returns a function intended for taking a client side websocket,
a 'websocket handler'"
function clientwsh(timeout::Int)
    function wsclienthandler(ws)
        s = "$(timeout)s_timeout_Julia"
        if !writeguarded(ws, s)
            @wslog "Could not write from ", ws, " called ", s, " at",  tss()
            return
        end
        sleep(timeout)
        WebSockets.close(ws, statusnumber = 1000, freereason = "$timeout seconds are up!")
    end
end
function checktasks()
    count = 0
    for clita in CLIENTTASKS
        count +=1
        if clita[2].state == :failed
            @error "Client websocket task ", clita[1], " => " clita[2] , " failed"
        end
    end
    count
end