text:
using LinearAlgebra
using UnPack, StaticArrays, IntervalSets

import ClimaCore: Fields, Domains, Meshes, Topologies, Spaces
import ClimaCore: slab
import ClimaCore.Operators
using ClimaCore.Geometry
import ClimaCore.Geometry: Abstract2DPoint

using ClimaCore.RecursiveApply

const parameters = (
    ϵ = 0.1,  # perturbation size for initial condition
    l = 0.5, # Gaussian width
    k = 0.5, # Sinusoidal wavenumber
    ρ₀ = 1.0, # reference density
    c = 2,
    g = 10,
)

domain = Domains.RectangleDomain(
    Geometry.XPoint(-2π) .. Geometry.XPoint(2π),
    Geometry.YPoint(-2π) .. Geometry.YPoint(2π),
    x1periodic = true,
    x2periodic = true,
)

function init_state(coord, p)
    @unpack x, y = coord
    # set initial state
    ρ = p.ρ₀

    # set initial velocity
    U₁ = cosh(y)^(-2)

    # Ψ′ = exp(-(y + p.l / 10)^2 / 2p.l^2) * cos(p.k * x) * cos(p.k * y)
    # Vortical velocity fields (u₁′, u₂′) = (-∂²Ψ′, ∂¹Ψ′)
    gaussian = exp(-(y + p.l / 10)^2 / 2p.l^2)
    u₁′ = gaussian * (y + p.l / 10) / p.l^2 * cos(p.k * x) * cos(p.k * y)
    u₁′ += p.k * gaussian * cos(p.k * x) * sin(p.k * y)
    u₂′ = -p.k * gaussian * sin(p.k * x) * cos(p.k * y)

    u = Geometry.UVVector(U₁ + p.ϵ * u₁′, p.ϵ * u₂′)
    # set initial tracer
    θ = sin(p.k * y)

    return (ρ = ρ, ρu = ρ * u, ρθ = ρ * θ)
end

function flux(state, p)
    @unpack ρ, ρu, ρθ = state
    u = ρu / ρ
    return (ρ = ρu, ρu = ((ρu ⊗ u) + (p.g * ρ^2 / 2) * I), ρθ = ρθ * u)
end

roe_average(ρ⁻, ρ⁺, var⁻, var⁺) =
    (sqrt(ρ⁻) * var⁻ + sqrt(ρ⁺) * var⁺) / (sqrt(ρ⁻) + sqrt(ρ⁺))

function roeflux(n, (y⁻, parameters⁻), (y⁺, parameters⁺))
    Favg = RecursiveApply.rdiv(flux(y⁻, parameters⁻) ⊞ flux(y⁺, parameters⁺), 2)

    λ = sqrt(parameters⁻.g)

    ρ⁻, ρu⁻, ρθ⁻ = y⁻.ρ, y⁻.ρu, y⁻.ρθ
    ρ⁺, ρu⁺, ρθ⁺ = y⁺.ρ, y⁺.ρu, y⁺.ρθ

    u⁻ = ρu⁻ / ρ⁻
    θ⁻ = ρθ⁻ / ρ⁻
    uₙ⁻ = u⁻' * n

    u⁺ = ρu⁺ / ρ⁺
    θ⁺ = ρθ⁺ / ρ⁺
    uₙ⁺ = u⁺' * n

    # in general thermodynamics, (pressure, soundspeed)
    p⁻ = (λ * ρ⁻)^2 * 0.5
    c⁻ = λ * sqrt(ρ⁻)

    p⁺ = (λ * ρ⁺)^2 * 0.5
    c⁺ = λ * sqrt(ρ⁺)

    # construct roe averges
    ρ = sqrt(ρ⁻ * ρ⁺)
    u = roe_average(ρ⁻, ρ⁺, u⁻, u⁺)
    θ = roe_average(ρ⁻, ρ⁺, θ⁻, θ⁺)
    c = roe_average(ρ⁻, ρ⁺, c⁻, c⁺)

    # construct normal velocity
    uₙ = u' * n

    # differences
    Δρ = ρ⁺ - ρ⁻
    Δp = p⁺ - p⁻
    Δu = u⁺ - u⁻
    Δρθ = ρθ⁺ - ρθ⁻
    Δuₙ = Δu' * n

    # constructed values
    c⁻² = 1 / c^2
    w1 = abs(uₙ - c) * (Δp - ρ * c * Δuₙ) * 0.5 * c⁻²
    w2 = abs(uₙ + c) * (Δp + ρ * c * Δuₙ) * 0.5 * c⁻²
    w3 = abs(uₙ) * (Δρ - Δp * c⁻²)
    w4 = abs(uₙ) * ρ
    w5 = abs(uₙ) * (Δρθ - θ * Δp * c⁻²)

    # fluxes!!!

    fluxᵀn_ρ = (w1 + w2 + w3) * 0.5
    fluxᵀn_ρu =
        (w1 * (u - c * n) + w2 * (u + c * n) + w3 * u + w4 * (Δu - Δuₙ * n)) *
        0.5
    fluxᵀn_ρθ = ((w1 + w2) * θ + w5) * 0.5

    Δf = (ρ = -fluxᵀn_ρ, ρu = -fluxᵀn_ρu, ρθ = -fluxᵀn_ρθ)
    RecursiveApply.rmap(f -> f' * n, Favg) ⊞ Δf
end

function volume!(dydt, y, (parameters,), t)
    div = Operators.WeakDivergence()
    rparameters = Ref(parameters)
    @. dydt = -(div(flux(y, rparameters)))
    return dydt
end

function add_face!(dydt, y, (parameters,), t)
    Operators.add_numerical_flux_internal!(roeflux, dydt, y, parameters)
    return dydt
end


text:
using Distributions
using ForwardDiff
dof=1
n=1000
chidist=Chisq(dof)
X1=rand(chidist,n)
print("First Moment Conditions is E[X1-$dof]=0")

mhat(d)=(1/n)sum(X1-d)

#First Step
W=1
f(d)=mhat(d)*W*mhat(d)
g(d)=ForwardDiff(f(d),d)

pr

include("sgdfunc.jl")
W,valid = SGDBabysitter(g, nn_predict::Function, maxIter, nHidden,
                        nParams, xtrain, ytrain, xvalid, yvalid )

ext: jl


text:
"""
A submodule that provides a selection of themes that can be used to colourise source code.

$(EXPORTS)
"""
module Themes

using DocStringExtensions

import ..Highlights: AbstractTheme, AbstractLexer

# Public interface.

export AbstractTheme, @S_str, @theme

# Colour, Style, Theme types.

"""
$(TYPEDEF)

Represents a single RGB colour value that can be 'active' or 'inactive'.
"""
struct RGB
    r::UInt8
    g::UInt8
    b::UInt8
    active::Bool

    function RGB(str::AbstractString)
        n = length(str)
        r, g, b = n == 3 ? rgb3(str) : n == 6 ? rgb6(str) :
            error("invalid colour code. Must be a 3 or 6 digit hex value.")
        return new(r, g, b, true)
    end

    RGB() = new(0x00, 0x00, 0x00, false)
end

"Convert a three digit hex string to a 3-tuple of `UInt8`s."
function rgb3(str)
    r = parse(UInt8, str[1], base=16)
    g = parse(UInt8, str[2], base=16)
    b = parse(UInt8, str[3], base=16)
    return (r << 4 + r, g << 4 + g, b << 4 + b)
end

"Convert a six digit hex string to a 3-tuple of `UInt8`s."
function rgb6(str)
    r1 = parse(UInt8, str[1], base=16)
    r2 = parse(UInt8, str[2], base=16)
    g1 = parse(UInt8, str[3], base=16)
    g2 = parse(UInt8, str[4], base=16)
    b1 = parse(UInt8, str[5], base=16)
    b2 = parse(UInt8, str[6], base=16)
    return (r1 << 4 + r2, g1 << 4 + g2, b1 << 4 + b2)
end

"""
$(TYPEDEF)

An internal type used to track colour scheme definition information such as foreground and
background colours as well as bold, italic, and underlining.
"""
struct Style
    fg::RGB
    bg::RGB
    bold::Bool
    italic::Bool
    underline::Bool

    function Style(spec::AbstractString)
        fg = bg = RGB()
        bold = italics = underline = false
        for part in split(spec, r"\s*;\s*")
            startswith(part, "fg:") && (fg = RGB(strip(part[4:end])))
            startswith(part, "bg:") && (bg = RGB(strip(part[4:end])))
            part == "bold" && (bold = true)
            part == "italic" && (italics = true)
            part == "underline" && (underline = true)
        end
        return new(fg, bg, bold, italics, underline)
    end
end

"""
$(SIGNATURES)

Construct a [`Style`](@ref) object.
"""
macro S_str(spec)
    Style(spec)
end

import ..Highlights.Tokens

"""
$(TYPEDEF)

Represents a "compiled" colour scheme.
"""
struct Theme
    base::Style
    styles::Vector{Style}
    Theme(base::Style, n::Integer) = new(base, Vector{Style}(undef, n))
end

function metadata end
function theme end

function maketheme(T)
    dict = metadata(T)
    n = length(Tokens.__TOKENS__)
    theme = Theme(get(dict, :style, S""), n)
    styles = get(dict, :tokens, Dict{Tokens.TokenValue, Style}())
    get!(styles, Tokens.TEXT, S"") # Set default TEXT if not already done.
    for (nth, t) in enumerate(Tokens.__TOKENS__)
        theme.styles[nth] = fallback(Tokens.__FALLBACKS__, styles, Tokens.TokenValue(t))
    end
    return theme
end

fallback(mapping, styles, token) = haskey(styles, token) ? styles[token] :
    fallback(mapping, styles, Tokens.TokenValue(mapping[token.value]))

"""
$(SIGNATURES)

Declare the theme definition for theme `T` based on `dict`. `dict` must be a `Dict` and `T`
must be a subtype of `AbstractTheme`.

# Examples

```jldoctest
julia> using Highlights.Themes

julia> abstract type CustomTheme <: AbstractTheme end

julia> @theme CustomTheme Dict(
           :name => "Custom",
           :tokens => Dict(
               # ...
           )
       );
```
"""
macro theme(T, dict)
    tx, dx = map(esc, (T, dict))
    quote
        $(Themes).metadata(::Type{$tx}) = $dx
        let data = $(Themes).maketheme($tx)
            $(Themes).theme(::Type{$tx}) = data::Theme
        end
    end
end

# Theme names.

export
    DefaultTheme,

    EmacsTheme,
    GitHubTheme,
    MonokaiTheme,
    MonokaiMiniTheme,
    PygmentsTheme,
    TangoTheme,
    TracTheme,
    VimTheme,
    VisualStudioTheme,
    XcodeTheme

"The default colour scheme with colours based on the Julia logo."
abstract type DefaultTheme <: AbstractTheme end

"A theme based on the Emacs colour scheme."
abstract type EmacsTheme <: AbstractTheme end
"A GitHub inspired colour scheme."
abstract type GitHubTheme <: AbstractTheme end
"A colour scheme similar to the Monokai theme."
abstract type MonokaiTheme <: AbstractTheme end
"A colour scheme similar to the Monokai theme that works on both light and dark backgrounds."
abstract type MonokaiMiniTheme <: AbstractTheme end
"Based on the default colour scheme used by the Pygments highlighter."
abstract type PygmentsTheme <: AbstractTheme end
"A Tango-inspired colour scheme."
abstract type TangoTheme <: AbstractTheme end
"Based on the default trac highlighter."
abstract type TracTheme <: AbstractTheme end
"A Vim 7.0 based colour scheme."
abstract type VimTheme <: AbstractTheme end
"A theme based on the default Visual Studio colours."
abstract type VisualStudioTheme <: AbstractTheme end
"A theme based on the default Xcode colour scheme."
abstract type XcodeTheme <: AbstractTheme end

# Theme definitions.

using ..Highlights.Tokens

include("themes/default.jl")

include("themes/pygments.jl")
include("themes/emacs.jl")
include("themes/github.jl")
include("themes/monokai.jl")
include("themes/monokaimini.jl")
include("themes/tango.jl")
include("themes/trac.jl")
include("themes/vim.jl")
include("themes/vs.jl")
include("themes/xcode.jl")

end # module


text:
# Form: a thing that is visible once drawn.

export polygon, rectangle, lines, curve, arc, circle, ellipse, text,
       hcenter, hleft, hright, vcenter, vtop, vbottom, empty_form

# A primitive form: typcially a shape or sequence of drawing operations.
abstract FormPrimitive

# Form is a (possibly empty) binary tree. Each node is a form primitive
# with an associated property. When the form is rendered, the graph is
# traversed from the root. The property at a node is applied to all downstream
# nodes (children, grandchildren, etc).
abstract Form

# The empty form, which forms the identity element of the Form monoid.
type EmptyForm <: Form end
const empty_form = EmptyForm()

contents(io, f::EmptyForm, n::Int, indent) = nothing
length(f::EmptyForm) = 0

# A non-empty sequence of form primitives each with an associated (possibly
# empty) property.
type FormTree <: Form
    primitive::Union(Nothing, FormPrimitive)
    property::Property
    children::List{Form}

    function FormTree(primitive::Union(Nothing, FormPrimitive),
                      property::Property,
                      children::List{Form})
        new(primitive, property, convert(List{Form}, children))
    end

    function FormTree(primitive::Union(Nothing, FormPrimitive),
                      property::Property)
        new(primitive, property, ListNil{Form}())
    end

    function FormTree(primitive::FormPrimitive)
        new(primitive, empty_property, ListNil{Form}())
    end

    # shallow copy constructor
    function FormTree(a::FormTree)
        new(a.primitive, a.property, a.children)
    end
end
length(f::FormTree) = length(f.children)

function contents(io, f::FormTree, n::Int, indent)
    if !(f.primitive === nothing)
        contents(io, f.primitive, 1, indent)
    end
    i = 1
    for fc in f.children
        contents(io, fc, n - 1, string(indent, "  "))
        if i > 10
            println(io, indent, "  ...")
            break
        end
        i += 1
    end
end

copy(a::FormTree) = FormTree(a)

children(a::FormTree) = a.children

function removable(a::FormTree)
    a.primitive === nothing && a.property === empty_property
end

# Conceptually, combining forms is tree joining by introducing a new root.
#
#   a      b               c
#  / \    / \    --->     / \
# X   Y  U   V           a   b
#                       / \ / \
#                      X  Y U  V
#
# There is a trick here to avoid an exceess of nop or "removable" nodes.
#
function combine(forms::Form...)
    children = ListNil{Form}()
    for form in forms
        if form === empty_form
            continue
        end
        if removable(form)
            for child in form.children
                children = cons(child, children)
            end
        else
            children = cons(form, children)
        end
    end

    FormTree(nothing, empty_property, children)
end

# Composition of properties into forms
function compose(a::FormTree, b::Property)
    FormTree(a.primitive, combine(a.property, b), a.children)
end

function compose(a::EmptyForm, b::Property)
    a
end

# Nop drawing of empty forms.
function draw(backend::Backend,
              t::NativeTransform,
              unit_box::BoundingBox,
              box::NativeBoundingBox,
              root_form::EmptyForm)
end

# Does a property in a node apply to it's siblings? No!
function draw(backend::Backend,
              t::NativeTransform,
              unit_box::BoundingBox,
              box::NativeBoundingBox,
              root_form::FormTree)

    S = {root_form}

    while !isempty(S)
        form = pop!(S)

        if form === :POP_PROPERTY
            pop_property(backend)
        elseif form === empty_form
            continue
        else
            if !is(form.property, empty_property)
                push!(S, :POP_PROPERTY)
                push_property(backend,
                              native_measure(form.property, t, unit_box,
                                             box, backend))
            end

            for child in children(form)
                push!(S, child)
            end

            if typeof(form) === FormTree
                if !(form.primitive === nothing)
                    draw(backend, t, unit_box, box, form.primitive)
                end
            else
                draw(backend, t, unit_box, box, form)
            end
        end
    end
end

type Lines <: FormPrimitive
    points::Vector{Point}
end

function contents(io, f::Lines, n::Int, indent)
    println(io, indent, "Line with ", length(f.points), " points")
end

function lines(points::XYTupleOrPoint...)
    FormTree(Lines([convert(Point, point) for point in points]))
end

function draw(backend::Backend, t::NativeTransform, unit_box::BoundingBox,
              box::NativeBoundingBox, form::Lines)
    native_form = Lines([native_measure(point, t, unit_box, box, backend)
                         for point in form.points])
    draw(backend, native_form)
end

type Curve <: FormPrimitive
    anchor0::Point
    ctrl0::Point
    ctrl1::Point
    anchor1::Point
end

function contents(io, f::Curve, n::Int, indent)
    println(io, indent, "Curve between ", f.anchor0, " and ", f.anchor1)
end

function curve(anchor0::XYTupleOrPoint, ctrl0::XYTupleOrPoint,
               ctrl1::XYTupleOrPoint, anchor1::XYTupleOrPoint)
    FormTree(Curve(convert(Point, anchor0), convert(Point, ctrl0),
                   convert(Point, ctrl1), convert(Point, anchor1)),
             fill(nothing))
end

function draw(backend::Backend, t::NativeTransform, unit_box::BoundingBox,
              box::NativeBoundingBox, form::Curve)
    native_form = Curve(native_measure(form.anchor0, t, unit_box, box, backend),
                        native_measure(form.ctrl0, t, unit_box, box, backend),
                        native_measure(form.ctrl1, t, unit_box, box, backend),
                        native_measure(form.anchor1, t, unit_box, box, backend))
    draw(backend, native_form)
end

type Polygon <: FormPrimitive
    points::Vector{Point}

end

function contents(io, f::Polygon, n::Int, indent)
    println(io, indent, "Polygon with ", length(f.points), " points")
end

function polygon(points::XYTupleOrPoint...)
    FormTree(Polygon([convert(Point, point) for point in points]))
end

function draw(backend::Backend, t::NativeTransform, unit_box::BoundingBox,
              box::NativeBoundingBox, form::Polygon)
    native_form = Polygon([native_measure(point, t, unit_box, box, backend)
                           for point in form.points])
    draw(backend, native_form)
end

function rectangle(x0::MeasureOrNumber, y0::MeasureOrNumber,
                   width::MeasureOrNumber, height::MeasureOrNumber)
    x0 = x_measure(x0)
    y0 = y_measure(y0)
    width = x_measure(width)
    height = y_measure(height)
    xy0 = Point(x0, y0)
    xy1 = Point(x0 + width, y0 + height)
    polygon(xy0, Point(xy1.x, xy0.y), xy1, Point(xy0.x, xy1.y))
end

function rectangle()
    rectangle(0.0w, 0.0h, 1.0w, 1.0h)
end

type Ellipse <: FormPrimitive
    center::Point
    x_point::Point
    y_point::Point

    function Ellipse(center::Point, x_point::Point, y_point::Point)
        new(center, x_point, y_point)
    end

    function Ellipse(x::MeasureOrNumber, y::MeasureOrNumber,
                     rx::MeasureOrNumber, ry::MeasureOrNumber)
        x = x_measure(x)
        y = y_measure(y)
        new(Point(x, y),
            Point(x + x_measure(rx), y),
            Point(x, y + y_measure(ry)))
    end
end
function contents(io, f::Ellipse, n::Int, indent)
    println(io, indent, "Ellipse centered at ", f.center)
end

function ellipse(x::MeasureOrNumber, y::MeasureOrNumber,
                 x_radius::MeasureOrNumber, y_radius::MeasureOrNumber)
    x = x_measure(x)
    y = y_measure(y)
    FormTree(Ellipse(Point(x, y),
                     Point(x + x_measure(x_radius), y),
                     Point(x, y + y_measure(y_radius))))
end

function ellipse()
    ellipse(0.5w, 0.5h, 0.5w, 0.5h)
end

# Despite the innacuracy, this is convenient is a few places.
typealias Circle Ellipse

function circle(x::MeasureOrNumber, y::MeasureOrNumber, radius::MeasureOrNumber)
    ellipse(x, y, radius, radius)
end

function draw(backend::Backend, t::NativeTransform, unit_box::BoundingBox,
              box::NativeBoundingBox, form::Ellipse)
    native_form = Ellipse(
        native_measure(form.center, t, unit_box, box, backend),
        native_measure(form.x_point, t, unit_box, box, backend),
        native_measure(form.y_point, t, unit_box, box, backend))
    draw(backend, native_form)
end

abstract HAlignment
type HLeft   <: HAlignment end
type HCenter <: HAlignment end
type HRight  <: HAlignment end

const hleft   = HLeft()
const hcenter = HCenter()
const hright  = HRight()

abstract VAlignment
type VTop    <: VAlignment end
type VCenter <: VAlignment end
type VBottom <: VAlignment end

const vtop    = VTop()
const vcenter = VCenter()
const vbottom = VBottom()

type Text <: FormPrimitive
    pos::Point
    value::String
    halign::HAlignment
    valign::VAlignment

    # Text forms need their own rotation field unfortunately, since there is no
    # way to give orientation with just a position point.
    t::NativeTransform
end

function text(x::MeasureOrNumber, y::MeasureOrNumber, value::String,
              halign::HAlignment, valign::VAlignment)
    FormTree(Text(Point(x_measure(x), y_measure(y)),
                  value, halign, valign, NativeTransform()))
end

function text(x::MeasureOrNumber, y::MeasureOrNumber, value::String)
    text(x, y, value, hleft, vbottom)
end

function draw(backend::Backend, t::NativeTransform, unit_box::BoundingBox,
              box::NativeBoundingBox, form::Text)
    form = Text(native_measure(form.pos, t, unit_box, box, backend),
                form.value, form.halign, form.valign, t)

    draw(backend, form)
end

function contents(io, f::Text, n::Int, indent)
    println(io, indent, "Text: ", f.value)
end


text:
struct Pendulum{T} <: Model
    n::Int
    m::Int
    d::Int

    mass::T    # mass
    b::T    # friction
    lc::T   # length to center of mass
    I::T    # inertia
    g::T    # gravity
end

function f(model::Pendulum, x, u, w)
    @SVector [x[2],
              (u[1] / ((model.mass * model.lc * model.lc))
                - model.g * sin(x[1]) / model.lc
                - model.b * x[2] / (model.mass * model.lc * model.lc))]
end

function k(model::Pendulum, x)
    @SVector [model.lc * sin(x[1]),
              -1.0 * model.lc * cos(x[1])]
end

n, m, d = 2, 1, 0
model = Pendulum(n, m, d, 1.0, 0.1, 0.5, 0.25, 9.81)

max_stars_repo_name: zacmanchester/DirectMotionPlanning

max_issues_repo_name: zacmanchester/DirectMotionPlanning

max_forks_repo_name: zacmanchester/DirectMotionPlanning


text:
export stop_PARSDMM

"""
implement stopping conditions for PARSDMM.jl
if activated, we set the flag stop==true and pass it back to PARSDMM.jl
"""
function stop_PARSDMM(
                    log_PARSDMM,
                    i                        ::Integer,
                    evol_rel_tol             ::TF,
                    feas_tol                 ::TF,
                    obj_tol                  ::TF,
                    adjust_rho               ::Bool,
                    adjust_gamma             ::Bool,
                    adjust_feasibility_rho   ::Bool,
                    ind_ref                  ::Integer,
                    counter                  ::Integer
                    ) where {TF<:Real}

    stop = false;

    #stop if objective value does not change and x is sufficiently feasible for all sets
    if i>6 && maximum(log_PARSDMM.set_feasibility[counter-1,:])<feas_tol && maximum(abs.( (log_PARSDMM.obj[i-5:i]-log_PARSDMM.obj[i-1-5:i-1])./log_PARSDMM.obj[i-1-5:i-1] )) < obj_tol
        println("stationary objective and reached feasibility, exiting PARSDMM (iteration ",i,")")
        stop=true;
    end

    #stop if x doesn't change significantly anyjore
    if i>5 && maximum(log_PARSDMM.evol_x[i-5:i])<evol_rel_tol
      println("relative evolution to small, exiting PARSDMM (iteration ",i,")")
      stop=true;
    end

    # fix rho to ensure regular ADMM convergence if primal residual does not decrease over a 20 iteration window
    if i>20 && adjust_rho==true && log_PARSDMM.r_pri_total[i]>maximum(log_PARSDMM.r_pri_total[(i-1):-1:max((i-50),1)])
      println("no primal residual reduction, fixing PARSDMM rho & gamma (iteration ",i,")")
      adjust_rho = false;
      adjust_feasibility_rho = false;
      adjust_gamma = false;
      if nprocs()>1
        [ @spawnat pid adjust_gamma for pid in workers() ]
        [ @spawnat pid adjust_rho for pid in workers() ]
        [ @spawnat pid adjust_feasibility_rho for pid in workers() ]
      end
      ind_ref = i;
    end

    #if rho is fixed and still no decrease in primal residual is observed over a window, we give up
    if adjust_rho==false && i>(ind_ref+25) && log_PARSDMM.r_pri_total[i]>maximum(log_PARSDMM.r_pri_total[(i-1):-1:max(ind_ref,max((i-50),1))])
      println("no primal residual reduction, exiting PARSDMM (iteration ",i,")")
      stop = true;
    end
    return stop,adjust_rho,adjust_gamma,adjust_feasibility_rho,ind_ref
end #stop_psdmm

max_stars_repo_name: slimgroup/SetIntersectionProjection

max_issues_repo_name: slimgroup/SetIntersectionProjection

max_forks_repo_name: slimgroup/SetIntersectionProjection


text:
# NOTE: This script is only supposed to be executed by our own build script.

import Pkg

# https://docs.travis-ci.com/user/languages/julia/#dependency-management
@static if VERSION < v"1.4"
    Pkg.add([Pkg.PackageSpec("Distributions", v"0.22"),
             Pkg.PackageSpec("Documenter", v"0.25"),
             Pkg.PackageSpec("GLPK", v"0.13"),
             Pkg.PackageSpec("GLPKMathProgInterface", v"0.5"),
             Pkg.PackageSpec("GR", v"0.44"),
             Pkg.PackageSpec("IntervalMatrices", v"0.6"),
             Pkg.PackageSpec("ModelingToolkit", v"0.8.0"),
             Pkg.PackageSpec("Optim", v"0.21"),
             Pkg.PackageSpec("RecipesBase", v"0.7"),
             Pkg.PackageSpec("TaylorModels", v"0.1")])
end


text:
"""
    center(A)

Centers columns of a 2d array

# Arguments

- A = 2d array of floats

# Value

2d array of floats

"""
function center(A::AbstractArray{Float64,2})

    # Row means of A
    m = mean(A, dims=1)
    # Initialize centered matrix
    W = zeros(size(A))

    # Subtract column means
    W .= A .- m

    return W
end

"""
    cov_est(resid)

Estimates error variance and its variance/covariance

# Arguments

- resid = 2d array of floats consisting of the residuals

# Value

Tuple
- est: 2d array of floats; estimate
- varest: 2d array of floats; variance/covariance estimate

2d array of floats

"""
function cov_est(resid::AbstractArray{Float64,2})

    # Dimensions of residuals
    n = size(resid, 1)
    p = size(resid, 2)

    # Centered residual matrix
    W = center(resid)

    # Allocate space for the estimates and their variances
    est = zeros(p, p)
    varEst = zeros(p, p)

    # Loop through the possible entries
    for i = 1:p
        for j = i:p
            # Multiply the ith and jth columns
            ww = W[:,i].*W[:,j]

            if i==j # Diagonal elements
                est[i,i] = (n/(n-1)) * mean(ww)
                varEst[i,i] = (n/(n-1)^2) * var(ww)
            else # Non-diagonal elements
                est[i,j] = est[j,i] = (n/(n-1)) * mean(ww)
                varEst[i,j] = varEst[j,i] = (n/(n-1)^2) * var(ww)
            end

        end
    end

    return est, varEst
end

"""
    shrink_sigma(resid, targetType)

Estimates variance of errors and the shrinkage coefficient

# Arguments

- resid = 2d array of floats consisting of the residuals
- targetType = string indicating the target type toward which to shrink the
  variance. Acceptable inputs are "A", "B", "C", and "D".
    - "A": Target is identity matrix
    - "B": Target is diagonal matrix with constant diagonal
    - "C": Target is has same diagonal element, and same off-diagonal element
    - "D": Target is diagonal matrix with unequal entries

# Value

Tuple
- sigma: 2d array of floats; shrunk estimated variance of errors
- lambda: floating scalar; estimated shrinkage coefficient
  (0 = no shrinkage, 1 = complete shrinkage)

# Reference

Ledoit, O., & Wolf, M. (2003). Improved estimation of the covariance matrix
    of stock returns with an application to portfolio selection. Journal of
    empirical finance, 10(5), 603-621.

"""
function shrink_sigma(resid::AbstractArray{Float64,2}, targetType::String)

    # Dimensions of resid
    (n, p) = size(resid)

    # Estimates and the variance of the error variance
    (est, varEst) = cov_est(resid)

    if targetType=="A"  # Shrink to identity
        # Create identity target matrix
        T = Matrix{Float64}(I, p, p)
        # Estimate optimal lambda
        lambda = sum(varEst) / sum((est-T).^2)

    elseif targetType=="B"  # Shrink to common variance
        # Create target matrix
        T = Matrix{Float64}(I, p, p) * mean(diag(est))
        # Estimate optimal lambda
        lambda = sum(varEst) / sum((est-T).^2)

    elseif targetType=="C"  # Shrink to equal variance and covariance
        v = mean(diag(est))
        c = (sum(est) - sum(diag(est))) / (n*(n-1))
        # Create target matrix
        T = fill(c,(p,p)) + (v-c) * Matrix{Float64}(I, p, p)
        # Estimate optimal lambda
        lambda = sum(varEst) / sum((est-T).^2)

    elseif targetType=="D"  # Shrink to zero correlation
        v = diag(est)
        # Create target matrix
        T = diagm(0 => v)
        # Estimate optimal lambda
        lambda = (sum(varEst) - sum(diag(varEst))) /
                 (sum(est.^2) - sum(diag(est).^2))
    end

    return lambda*T + (1-lambda)*est, lambda
end


text:
# ---
# notebook: nothing
# ---


text:
# This file was generated, do not modify it. # hide
r = range(atom,
          :min_samples_split,
          lower=2,
          upper=100, scale=:log)

mach = machine(atom, X, y)

figure()
curve = learning_curve!(mach,
                        range=r,
                        measure=mav,
                        resampling=CV(nfolds=9),
                        verbosity=0)

plot(curve.parameter_values, curve.measurements)
xlabel(curve.parameter_name)

savefig(joinpath(@OUTPUT, "e1.svg")) # hide

max_issues_repo_name: alan-turing-institute/MLJTutorials

max_forks_repo_name: alan-turing-institute/MLJTutorials


text:
#
# While a problem consists of all things one has to know about the Optimization
# problem itself (independent of the solver), the options collect parameters,
# that steer the solver (indipendent of the problem at hand)
#
import Base: copy

export StoppingCriterion, Stepsize
export EvalOrder, LinearEvalOrder, RandomEvalOrder, FixedRandomEvalOrder
export Options, getOptions, getReason
export IsOptionsDecorator

export Action, StoreOptionsAction
export hasStorage, getStorage, updateStorage!

"""
    IsOptionsDecorator{O}

A trait to specify that a certain `Option` decorates, i.e. internally
stores the original [`Options`](@ref) under consideration.
"""
@traitdef IsOptionsDecorator{O}

"""
    Options

A general super type for all options.

# Fields
The following fields are assumed to be default. If you use different ones,
provide the access functions accordingly
* `x` an [`MPoint`](@ref) with the current iterate
* `stop` a [`StoppingCriterion`](@ref).

"""
abstract type Options end
#
# StoppingCriterion meta
#
@doc doc"""
    StoppingCriterion

An abstract type for the functors representing stoping criteria, i.e. they are
callable structures. The naming Scheme follows functions, see for
example [`stopAfterIteration`](@ref).

Every StoppingCriterion has to provide a constructor and its function has to have
the interface `(p,o,i)` where a [`Problem`](@ref) as well as [`Options`](@ref)
and the current number of iterations are the arguments and returns a Bool whether
to stop or not.

By default each `StoppingCriterion` should provide a fiels `reason` to provide
details when a criteion is met (and that is empty otherwise).
"""
abstract type StoppingCriterion end
#
#
# StepsizeOptions
#
#
"""
    Stepsize

An abstract type for the functors representing step sizes, i.e. they are callable
structurs. The naming scheme is `TypeOfStepSize`, e.g. `ConstantStepsize`.

Every Stepsize has to provide a constructor and its function has to have
the interface `(p,o,i)` where a [`Problem`](@ref) as well as [`Options`](@ref)
and the current number of iterations are the arguments
and returns a number, namely the stepsize to use.

# See also
[`Linesearch`](@ref)
"""
abstract type Stepsize end
#
#
# Evalualtion Orders
#
#
"""
    EvalOrder
type for specifying an evaluation order for any cyclicly evaluated algorithms
"""
abstract type EvalOrder end
"""
    LinearEvalOrder <: EvalOrder
evaluate in a linear order, i.e. for each cycle of length l evaluate in the
order 1,2,...,l.
"""
mutable struct LinearEvalOrder <: EvalOrder end
"""
    RandomEvalOrder <: EvalOrder
choose a random order for each evaluation of the l functionals.
"""
mutable struct RandomEvalOrder <: EvalOrder end
"""
    FixedRandomEvalOrder <: EvalOrder
Choose a random order once and evaluate always in this order, i.e. for
l elements there is one chosen permutation used for each iteration cycle.
"""
mutable struct FixedRandomEvalOrder <: EvalOrder end

@doc doc"""
    getOptions(O)

return the undecorated [`Options`](@ref) of the (possibly) decorated `O`.
As long as your decorated options stores the options within `o.options` and
implements the `SimpleTrait` `IsOptionsDecorator`, this is behaviour is optained
automatically.
"""
getOptions(O) = error("Not implemented for types that are not `Options`")
# this might seem like a trick/fallback just for documentation reasons
@traitfn getOptions(o::O) where {O <: Options; !IsOptionsDecorator{O}} = o
@traitfn getOptions(o::O) where {O <: Options; IsOptionsDecorator{O}} = getOptions(o.options)

@doc doc"""
    getReason(o)

return the current reason stored within the [`StoppingCriterion`](@ref) from
within the [`Options`](@ref) This reason is empty if the criterion has never
been met.
"""
getReason(o::O) where O <: Options = getReason( getOptions(o).stop )

#
# Common Actions for decorated Options
#
@doc doc"""
    Action

a common `Type` for `Actions` that might be triggered in decoraters,
for example [`DebugOptions`](@ref) or [`RecordOptions`](@ref).
"""
abstract type Action end

@doc doc"""
    StoreTupleAction <: Action

internal storage for [`Action`](@ref)s to store a tuple of fields from an
[`Options`](@ref)s

This functor posesses the usual interface of functions called during an
iteration, i.e. acts on `(p,o,i)`, where `p` is a [`Problem`](@ref),
`o` is an [`Options`](@ref) and `i` is the current iteration.

# Fields
* `values` – a dictionary to store interims values based on certain `Symbols`
* `keys` – an `NTuple` of `Symbols` to refer to fields of `Options`
* `once` – whether to update the internal values only once per iteration
* `lastStored` – last iterate, where this `Action` was called (to determine `once`

# Constructiors

    StoreOptionsAction([keys=(), once=true])

Initialize the Functor to an (empty) set of keys, where `once` determines
whether more that one update per iteration are effective

    StoreOptionsAction(keys, once=true])

Initialize the Functor to a set of keys, where the dictionary is initialized to
be empty. Further, `once` determines whether more that one update per iteration
are effective, otherwise only the first update is stored, all others are ignored.
"""
mutable struct StoreOptionsAction <: Action
    values::Dict{Symbol,<:Any}
    keys::NTuple{N,Symbol} where N
    once::Bool
    lastStored::Int
    StoreOptionsAction(keys::NTuple{N,Symbol} where N = NTuple{0,Symbol}(),once=true) = new(Dict{Symbol,Any}(), keys, once,-1 )
end
function (a::StoreOptionsAction)(p::P,o::O,i::Int) where {P <: Problem, O <: Options}
    #update values (maybe only once)
    if !a.once || a.lastStored != i
        merge!(a.values, Dict( key => getproperty(o,key) for key in a.keys) )
    end
    a.lastStored = i
end
"""
    getStorage(a,key)

return the internal value of the [`StoreOptionsAction`](@ref) `a` at the
`Symbol` `key`.
"""
getStorage(a::StoreOptionsAction,key) = a.values[key]
"""
    getStorage(a,key)

return whether the [`StoreOptionsAction`](@ref) `a` has a value stored at the
`Symbol` `key`.
"""
hasStorage(a::StoreOptionsAction,key) = haskey(a.values,key)
"""
    updateStorage!(a,o)

update the [`StoreOptionsAction`](@ref) `a` internal values to the ones given on
the [`Options`](@ref) `o`.
"""
updateStorage!(a::StoreOptionsAction,o::O) where {O <: Options} = updateStorage!(a, Dict( key => getproperty(o, key) for key in a.keys) )
"""
    updateStorage!(a,o)

update the [`StoreOptionsAction`](@ref) `a` internal values to the ones given in
the dictionary `d`. The values are merged, where the values from `d` are preferred.
"""
function updateStorage!(a::StoreOptionsAction,d::Dict{Symbol,<:Any}) where {O <: Options}
    merge!(a.values, d)
    # update keys
    a.keys = Tuple( keys(a.values) )
end


text:
function FreeConvectionNDE(NN, ds; iterations=nothing)
    weights, reconstruct = Flux.destructure(NN)

    T = ds["T"]
    wT = ds["wT"]
    FT = eltype(T)
    Nz = size(T, 3)
    Nt = size(T, 4)
    zc = znodes(T)
    zf = znodes(wT)
    times = T.times

    H = abs(zf[1]) # Domain depth/height
    τ = times[end]  # Simulation length

    Δẑ = diff(zc)[1] / H  # Non-dimensional grid spacing
    Dzᶜ = Dᶜ(Nz, Δẑ) # Differentiation matrix operator
    Dzᶜ = convert(Array{FT}, Dzᶜ)

    if isnothing(iterations)
        iterations = 1:length(times)
    end

    """
    Non-dimensional PDE is

        ∂T/∂t = - σ_wT/σ_T * τ/H * ∂/∂z(wT)
    """
    function ∂T∂t(T, p, t)
        weights = p[1:end-6]
        bottom_flux, top_flux, σ_T, σ_wT, H, τ = p[end-5:end]

        NN = reconstruct(weights)
        wT_interior = NN(T)
        wT = [bottom_flux; wT_interior; top_flux]
        ∂z_wT = Dzᶜ * σ_wT/σ_T * τ/H * wT
        return -∂z_wT
    end

    tspan = FT.( (0.0, maximum(iterations) / Nt) )
    saveat = range(tspan[1], tspan[2], length=length(iterations))

    # See: https://github.com/SciML/DiffEqFlux.jl/blob/449efcecfc11f1eab65d0e467cf57db9f5a5dbec/src/neural_de.jl#L66-L67
    # We set the initial condition to `nothing`. Then we will set it to some actual initial condition when calling `solve`.
    ff = ODEFunction{false}(∂T∂t, tgrad=DiffEqFlux.basic_tgrad)
    return ODEProblem{false}(ff, nothing, tspan, saveat=saveat)
end

function FreeConvectionNDEParameters(ds, T_scaling, wT_scaling)
    wT = ds["wT"]
    zf = znodes(wT)
    times = wT.times

    H = abs(zf[1]) # Domain depth/height
    τ = times[end]  # Simulation length

    bottom_flux = wT_scaling(interior(wT)[1, 1, 1, 1])
    top_flux = wT_scaling(ds.metadata["temperature_flux"])

    FT = eltype(wT)
    return FT.([bottom_flux, top_flux, T_scaling.σ, wT_scaling.σ, H, τ])
end


text:
module Network
    using CBinding
    c`-std=c99 -fparse-all-comments`
    const c"__u8" = UInt8
    const c"__u16" = UInt16
    const c"__u32" = UInt32
    const c"__be8" = UInt8
    const c"__be16" = UInt16
    const c"__be32" = UInt32
    const c"__sum16" = UInt16
    c"#include <linux/if_ether.h>"j
    c"#include <linux/ip.h>"j
    c"#include <linux/in.h>"j
    c"#include <linux/tcp.h>"j
end


text:
function _call(re)
    # Make it dumb so that it's easier in debugger:
    f = re.f
    args = re.args
    kwargs = re.kwargs
    JuliaInterpreter.@bp
    # NOTE: Hit `s` (maybe twice) to step into the recorded function, if you are in a debugger.
    return f(args...; kwargs...)
end


text:
"""
    f(x)=scalaransatz(x)

Evaluate the simple vortex core ansatz at radial point `x`.

The ansatz is the rational function approximation to the vortex core
```math
f(x)=\\sqrt{\\frac{x^2}{1+x^2}}
```
"""
scalaransatz(x) = sqrt(x^2/(1+x^2))

"""
    core = Ansatz(ψ,ξ,Λ)

Construct a fast interpolation for the vortex core ansatz.

Returns a callable type that can be evaluated at position `x`.

# Arguments
- `ψ::Interpolation=ψa`: ansatz wavefunction.
- `ξ::Float64=1.0`: healing length.
- `Λ::Float64=0.8249`: slope of wavefunction at vortex core.

# Examples
```jldoctest
julia> f = Ansatz()
julia> x = LinRange(0,10,100)
julia> y = f.(x)
```

See also: [`Exact`](@ref)
"""
Ansatz() = Ansatz(ψa,1.0,Λ)

function (core::Ansatz)(x)
    @unpack f,ξ,Λ = core
    return f(Λ*x/ξ)
end
(core::Ansatz)(x,y) = core(hypot(x,y))

"""
    core = Exact(ψ,ξ)

Construct fast interpolation for exact vortex core.

Returns a callable type that can be evaluated at position `x`.

# Examples
```jldoctest
julia> f = Exact()
julia> x = LinRange(0,10,100)
julia> y = f.(x)
```

See also: [`Ansatz`](@ref)
"""
Exact(ξ::Float64) = Exact(VortexDistributions.ψi,ξ::Float64)
Exact() = Exact(1.0)
function (core::Exact)(x)
    @unpack ξ,f = core
    return f(x/ξ)
end
(core::Exact)(x,y) = core(hypot(x,y))

"""
    vort = ScalarVortex(ξ::Float64=1.0;pv::PointVortex)
    vort = ScalarVortex(pv::Array{PointVortex,1})
    vort = ScalarVortex(pv::PointVortex)
    vort = ScalarVortex(ξ::Float64,pv::Array{PointVortex,1})
    vort = ScalarVortex(ξ::Array{Float64,1},pv::Array{PointVortex,1})

Construct a scalar vortex with healing length `ξ` and point vortex coordinates `pv`.

# Arguments
- `ξ::Float64=1.0`: healing length.
- `pv::PointVortex`: scalar or vector of point vortices.

# Examples
```jldoctest
julia> pv = PointVortex(-2.,3.,1)
julia> v1 = ScalarVortex(pv)
```
returns a scalarvortex `v1` suitable for vortex construction in a wavefunction.

See also: [`vortex!`](@ref), [`randScalarVortex`](@ref), [`randVortex`](@ref), [`PointVortex`](@ref), [`randPointVortex`](@ref)
"""
ScalarVortex(vort::PointVortex) = ScalarVortex(Exact(),vort)
function (s::ScalarVortex{T})(x,y) where T <: CoreShape
    @unpack xv,yv,qv = s.vort
    return s.core(x - xv, y - yv)*exp(im*qv*atan(y - yv,x - xv))
end
ScalarVortex(ξ::Float64,pv::Array{PointVortex,1}) = ScalarVortex.([Exact(ξ)],pv::Array{PointVortex,1})
ScalarVortex(ξ::Array{Float64,1},pv::Array{PointVortex,1}) = @. ScalarVortex(Exact(ξ),pv::Array{PointVortex,1})
ScalarVortex(ξ::Float64,pv::PointVortex) = ScalarVortex(ξ,[pv])
ScalarVortex(pv::PointVortex) = ScalarVortex(1.0,pv::PointVortex)
ScalarVortex(pv::Array{PointVortex,1}) = ScalarVortex(1.0,pv)

"""
    rv = randScalarVortex()
    rv = randScalarVortex(n::Int)
    rv = randScalarVortex(ψ::Field)
    rv = randScalarVortex(n::Int,ψ::Field)

Sample `n` random scalar vortices, using the sptial domain of the field `ψ`.

See also: [`Field`](@ref), [`randVortex`](@ref), [`ScalarVortex`](@ref), [`randPointVortex`](@ref)
"""
randScalarVortex() = ScalarVortex(randPointVortex())
randScalarVortex(n) = ScalarVortex.(randPointVortex(n))
randScalarVortex(psi::Field) = ScalarVortex(Exact(),randPointVortex(psi))
randScalarVortex(n,psi::Field) = ScalarVortex.([Exact()],randPointVortex(n,psi))

randVortex() = randScalarVortex()
randVortex(n) = randScalarVortex(n)
randVortex(n,psi::Field) = randScalarVortex(n,psi)

"""
    vortex!(ψ<:Field, vort<:Vortex)

Density and phase imprint a vortex onto the field `ψ`, writing in place.

See also: [`Field`](@ref), [`ScalarVortex`](@ref), [`PointVortex`](@ref), [`randPointVortex`](@ref), [`randScalarVortex`](@ref)
"""
function vortex!(psi::F,vort::ScalarVortex{T}) where {T <: CoreShape, F<:Field}
    @unpack ψ,x,y = psi
    @. ψ *= vort.(x,y')
    @pack! psi = ψ
end
function vortex!(psi::F,vort::S) where {F <: Field, S <: Vortex}
    vortex!(psi,ScalarVortex(vort))
end
function vortex!(psi::F,vort::Array{S}) where {F <: Field, S <: Vortex}
    for j in eachindex(vort)
        vortex!(psi,vort[j])
    end
end

function randVortexField(n)
    Nx = 400; Ny = 400
    Lx = 200; Ly = 200
    x = LinRange(-Lx / 2, Ly / 2, Nx)
    y = LinRange(-Ly / 2, Ly / 2, Ny)

    psi0 = one.(x*y') |> complex; psi = Torus(psi0,x,y)
    vort = randVortex(n,psi)
    vortex!(psi,vort)
    return psi,PointVortex(vort)
end

#--- periodic dipole phase
# Billam et al, PRL 112, 145301 (2014), Supplemental
H(x) = x > 0. ? 1.0 : 0.0
shift(x,xi) = x - xi
tans(x,xk) = tan((shift(x,xk) - π)*0.5)
tanhs(x,xk,j) = tanh((shift(x,xk) + 2*π*j)*0.5)

function kernel(x,y,xp,yp,xn,yn,j)
    return atan(tanhs(y,yn,j)*tans(x,xn)) -
    atan(tanhs(y,yp,j)*tans(x,xp))
end

# arbitrary domains and dipole sizes:
function thetad(x,y,xp,yp,xn,yn)
    s = 0.0
    for j = -5:5
        s += kernel(x,y,xp,yp,xn,yn,j)
    end
    s += π*(H(shift(x,xp)) - H(shift(x,xn))) - y*(xp - xn)/(2*π)
    return s - x*H(abs(yp - yn) - π) + y*H(abs(xp - xn) - π)
end

function Thetad(x,y,xp,yp,xn,yn)
    Lx = x[end]-x[1]
    Ly = y[end]-y[1]
    return @. thetad(x*2*pi/Lx,y'*2*pi/Ly,xp*2*pi/Lx,yp*2*pi/Ly,xn*2*pi/Lx,yn*2*pi/Ly)
end

function periodic_dipole!(psi::F,dip::Array{ScalarVortex{T},1}) where {T <: CoreShape, F<:Field}
    @assert length(dip) == 2
    @assert dip[1].vort.qv + dip[2].vort.qv == 0
    @assert hypot(dip[1].vort.xv-dip[2].vort.xv,dip[1].vort.yv-dip[2].vort.yv) >= 2*pi
    @unpack ψ,x,y = psi
    (dip[1].vort.qv > 0) ? (jp = 1;jn = 2) : (jp = 2;jn = 1)
    vp = rawData(dip[jp].vort)[1:2]
    vn = rawData(dip[jn].vort)[1:2]
    @. ψ *= abs(dip[jn](x,y')*dip[jp](x,y'))
    ψ .*= exp.(im*Thetad(x,y,vp...,vn...))
    @pack! psi = ψ
end

#Chebyshev methods

function gpecore_exact(K,L=2,N=100,R = K)
    #currently r does  nothing!
    #N = 100
    #L = 2
    #Κ = 1
    #R = 2 # Stretching coordinate.  R ~ kappa seems to be a good ballpark

# Convert Chebyshev grid to [0,L] from [1 -1];
z,Dz = getChebDMatrix(N)
blank,D2z = getChebD2Matrix(N)
z = vec(reverse(z,dims=2))
z = (z)*L./2
Dz = -2*Dz./L
D2z = 4*D2z./L.^2

# y is the physical coordinate, range [0 infinity] (i.e y = r)
y =  R*(1 .+z)./(1 .-z)

#initial guess based on ansatz
ψ = @. y/sqrt(y^2 + Λ^-2)
ψ[1] = 0
ψ[end] = 1
ψ0 = vec(ψ)

Q = vec(z.^2 .- 2*z .+ 1)
Qmat = repeat(Q,1,N)
Zmat = repeat(2*(z .-1),1,N)
Ymat = repeat(1 ./y,1,N);

# Second Derivative
residuals = -0.5*( (Q.*(D2z*ψ) + 2*(z .-1).*(Dz*ψ) ).*Q/(4*R^2)
        + (Q/(2*R) ./y).*(Dz*ψ) )+ 0.5*K^2 *ψ ./y.^2 + ψ.^3 .- ψ
residuals[1] = 0; residuals[end] =0

while sum(abs.(residuals).^2) > 1e-12
    # Second Derivative
    residuals = -0.5*( (Q.*(D2z*ψ) + 2*(z .-1).*(Dz*ψ) ).*Q/(4*R^2)
            + (Q/(2*R) ./y).*(Dz*ψ) )+ 0.5*K^2 *ψ ./y.^2 + ψ.^3 .- ψ
    residuals[1] = 0; residuals[end] =0

    Jacobi = -0.5*( (Qmat.*(D2z) + Zmat.*(Dz) ).*(Qmat/(4*R^2))
    + (Qmat/(2*R).*Ymat).*(Dz) ) + diagm(0 => 0.5*K^2 ./y.^2)+ diagm(0 => 3*ψ.^2 .- 1)

    Jacobi[1,:] = [1 zeros(1,N-1)]
    Jacobi[N,:] = [zeros(1,N-1) 1]

    Δ = vec(-4/7*(Jacobi\residuals))

    ψ = ψ + Δ
    ψ[1] = 0
    ψ[end] =1
end
    res = sum(abs.(residuals).^2)
    return y,ψ,res
end

"""
`x, DM =  chebdif(N,M)`

computes the differentiation matrices `D1, D2, ..., DM` on Chebyshev nodes.

Input:

`N`         Size of differentiation matrix.

`M`         Number of derivatives required (integer).

Note      `0 < M <= N-1`.

Output:

DM        `DM(1:N,1:N,ell)` contains ell-th derivative matrix, ell=1..M.

The code implements two strategies for enhanced
accuracy suggested by W. Don and S. Solomonoff in
SIAM J. Sci. Comp. Vol. 6, pp. 1253--1268 (1994).
The two strategies are (a) the use of trigonometric
identities to avoid the computation of differences
x(k)-x(j) and (b) the use of the "flipping trick"
which is necessary since sin t can be computed to high
relative precision when t is small whereas sin (pi-t) cannot.
Note added May 2003:  It may, in fact, be slightly better not to
implement the strategies (a) and (b).   Please consult the following
paper for details:   "Spectral Differencing with a Twist", by
R. Baltensperger and M.R. Trummer, to appear in SIAM J. Sci. Comp.

J.A.C. Weideman, S.C. Reddy 1998.  Help notes modified by
JACW, May 2003."""
function chebdif(N, M)
   @assert 0<M<=N-1

   n1,n2 = (floor(N/2)|> Int,ceil(N/2)|> Int)  # Indices used for flipping trick.

   k = collect(0:N-1)' # Compute theta vector.
   th = k*π/(N-1)

   x = sin.(π*collect(N-1:-2:1-N)'./(2*(N-1))) # Compute Chebyshev points.

   T = repeat(th/2,N,1)'
   DX = 2*sin.(T'.+T).*sin.(T'.-T) # Trigonometric identity.
   DX = [DX[1:n1,:]; -reverse(reverse(DX[1:n2,:],dims=2),dims=1)] # Flipping trick.
   DX[diagind(DX)] .= 1 # Put 1's on the main diagonal of DX.

   C = Matrix(Toeplitz((-1).^k',(-1).^k'))  # C is the matrix with
   C[1,:] = C[1,:]*2; C[N,:] = C[N,:]*2  # entries c(k)/c(j)
   C[:,1] = C[:,1]/2; C[:,N] = C[:,N]/2

   Z = 1 ./DX  # Z contains entries 1/(x(k)-x(j))
   Z[diagind(Z)] .= 0  # with zeros on the diagonal.

   D = Matrix{Float64}(I, N, N)  # D contains diff. matrices.

   DM = zeros(N,N,M)

   for ell = 1:M
       D = ell*Z.*(C.*repeat(diag(D),1,N) - D) # Off-diagonals
       D[diagind(D)] = -sum(D',dims=1)  # Correct main diagonal of D
       DM[:,:,ell] = D  # Store current D in DM
   end

   return x, DM
   end

function getChebDMatrix(n)
   z, M = chebdif(n, 1)
   Dz = M[:,:,1]
   return z, Dz
end

function getChebD2Matrix(n)
   z, M = chebdif(n, 2)
   D2z = M[:,:,2]
   return z, D2z
end


text:
function f(x::Nothing) end
f(42)

max_stars_repo_name: gabriel-doriath-dohler/pjulia

max_issues_repo_name: gabriel-doriath-dohler/pjulia

max_forks_repo_name: gabriel-doriath-dohler/pjulia


text:
## Diversity Measures

"""
    shannon(v::Union{AbstractVector, AbstractSparseMatrix})
    shannon(abt::AbstractAbundanceTable)

Computes the Shannon alpha diversity metric for a vector.
When called on an `AbstractAbundanceTable`,
returns a 1 x nsamples matrix with 1 entry per sample.
See also [`shannon!`](@ref).
"""
function shannon(v::Union{AbstractVector{T}, AbstractSparseMatrix{T}}) where T<:Real
    total = sum(v)
    relab = map(x-> x/total, v)
    return -sum([log(x^x) for x in relab])
end

function shannon(abt::AbstractAbundanceTable)
    return hcat((shannon(abundances(abt[:, s])) for s in samplenames(abt))...)
end

"""
    shannon!(abt::AbstractAbundanceTable; overwrite=false)

Adds a `:shannon` entry to the metadata for each sample in `abt`
with the Shannon alpha diversity of that sample (see [`shannon`](@ref)).
If `overwrite=false` (the default), uses `insert!` to perform this operation,
so an error will be thrown if any sample already contains a `:shannon` entry.
Otherwise, uses `set!`.
"""
function shannon!(abt::AbstractAbundanceTable; overwrite=false)
    func! = overwrite ? set! : insert!
    for s in samples(abt)
        col = abt[:, s]
        sh = shannon(abundances(col))
        func!(s, :shannon, sh)
    end
    return abt
end

"""
    ginisimpson(v::Union{AbstractVector, AbstractSparseMatrix})
    ginisimpson(abt::AbstractAbundanceTable, overwrite=false)

Computes the Gini-Simpson alpha diversity metric for a vector.
When called on an `AbstractAbundanceTable`,
returns a 1 x nsamples matrix with 1 entry per sample.
See also [`ginisimpson!`](@ref).
"""
function ginisimpson(v::Union{AbstractVector{T}, AbstractSparseMatrix{T}}) where T<:Real
    total = sum(v)
    relab = map(x-> x/total, v)
    return 1 - sum([x^2 for x in relab])
end

function ginisimpson(abt::AbstractAbundanceTable)
    return hcat((ginisimpson(abundances(abt[:, s])) for s in samplenames(abt))...)
end

"""
    ginisimpson!(abt::AbstractAbundanceTable; overwrite=false)

Adds a `:ginisimpson` entry to the metadata for each sample in `abt`
with the Gini-Simpson alpha diversity of that sample (see [`ginisimpson`](@ref)).
If `overwrite=false` (the default), uses `insert!` to perform this operation,
so an error will be thrown if any sample already contains a `:ginisimpson` entry.
Otherwise, uses `set!`.
"""
function ginisimpson!(abt::AbstractAbundanceTable; overwrite=false)
    func! = overwrite ? set! : insert!
    for s in samples(abt)
        col = abt[:, s]
        sh = ginisimpson(abundances(col))
        func!(s, :ginisimpson, sh)
    end
    return abt
end

"""
    braycurtis(abt::AbstractAbundanceTable)

Returns a pairwise Bray-Curtis dissimilarity matrix.
"""
braycurtis(abt::AbstractAbundanceTable) = pairwise(BrayCurtis(), collect(abundances(abt)), dims=2)

"""
    jaccard(abt::AbstractAbundanceTable)

Returns a pairwise Jaccard distance matrix.
"""
jaccard(abt::AbstractAbundanceTable) = pairwise(Jaccard(), collect(abundances(abt)), dims=2)

"""
    hellinger(abt::AbstractAbundanceTable)

Returns a pairwise Hellinger distance matrix.
"""
hellinger(abt::AbstractAbundanceTable) = pairwise(HellingerDist(), collect(abundances(abt)), dims=2)

"""
    pcoa(abt::AbstractAbundanceTable, f=braycurtis)

Returns eigenvectors from fitting `MDS` to a distance metric generated by `f`,
by default `braycurtis`.
"""
pcoa(abt::AbstractAbundanceTable, f=braycurtis) = fit(MDS, f(abt), distances = true)


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

name = "JSON_C"
version = v"0.15.0"

# Collection of sources required to complete build
sources = [
    ArchiveSource("https://s3.amazonaws.com/json-c_releases/releases/json-c-0.15.tar.gz", "b8d80a1ddb718b3ba7492916237bbf86609e9709fb007e7f7d4322f02341a4c6")
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir
cd json-c-0.15
cmake -DCMAKE_INSTALL_PREFIX=$prefix -DCMAKE_TOOLCHAIN_FILE=${CMAKE_TARGET_TOOLCHAIN} -DCMAKE_BUILD_TYPE=Release
make
make install
exit
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = supported_platforms()

# The products that we will ensure are always built
products = [
    LibraryProduct("libjson-c", :libjson_c)
]

# Dependencies that must be installed before this package can be built
dependencies = Dependency[
]

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

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil


text:
struct CPUEvent <: Event
    task::Core.Task
end
isdone(ev::CPUEvent) = Base.istaskdone(ev.task)
failed(ev::CPUEvent) = Base.istaskfailed(ev.task)

function Event(::CPU)
    return NoneEvent()
end

"""
    Event(f, args...; dependencies, progress, sticky=true)

Run function `f` with `args` in a Julia task. If `sticky` is `true` the task
is run on the thread that launched it.
"""
function Event(f, args...; dependencies=nothing, progress=nothing, sticky=true)
    T = Task() do
        wait(MultiEvent(dependencies), progress)
        f(args...)
    end
    T.sticky = sticky
    Base.schedule(T)
    return CPUEvent(T)
end

wait(ev::Union{CPUEvent, NoneEvent, MultiEvent}, progress=nothing) = wait(CPU(), ev, progress)
wait(::CPU, ev::NoneEvent, progress=nothing) = nothing

function wait(cpu::CPU, ev::MultiEvent, progress=nothing)
    events = ev.events
    N = length(events)
    alldone = ntuple(i->false, N)

    while !all(alldone)
        alldone = ntuple(N) do i
            if alldone[i]
                true
            else
                isdone(events[i])
            end
        end
        if progress === nothing
            yield()
        else
            progress()
        end
    end

    if any(failed, events)
        ex = CompositeException()
        for event in events
            if failed(event) && event isa CPUEvent
                push!(ex, Base.TaskFailedException(event.task))
            end
        end
        throw(ex)
    end
end

function wait(::CPU, ev::CPUEvent, progress=nothing)
    if progress === nothing
        wait(ev.task)
    else
        while !Base.istaskdone(ev.task)
            progress()
        end
        if Base.istaskfailed(ev.task)
            throw(Base.TaskFailedException(ev.task))
        end
    end
end

function async_copy!(::CPU, A, B; dependencies=nothing, progress=nothing)
    Event(copyto!, A, B, dependencies=dependencies, progress=progress)
end

function (obj::Kernel{CPU})(args...; ndrange=nothing, workgroupsize=nothing, dependencies=nothing, progress=nothing)
    ndrange, workgroupsize, iterspace, dynamic = launch_config(obj, ndrange, workgroupsize)

    if dependencies isa Event
        dependencies = (dependencies,)
    end

    if length(blocks(iterspace)) == 0
        return MultiEvent(dependencies)
    end

    Event(__run, obj, ndrange, iterspace, args, dynamic,
          dependencies=dependencies, progress=progress)
end

function launch_config(kernel::Kernel{CPU}, ndrange, workgroupsize)
    if ndrange isa Integer
        ndrange = (ndrange,)
    end
    if workgroupsize isa Integer
        workgroupsize = (workgroupsize, )
    end

    if KernelAbstractions.workgroupsize(kernel) <: DynamicSize && workgroupsize === nothing
        workgroupsize = (1024,) # Vectorization, 4x unrolling, minimal grain size
    end
    iterspace, dynamic = partition(kernel, ndrange, workgroupsize)
    # partition checked that the ndrange's agreed
    if KernelAbstractions.ndrange(kernel) <: StaticSize
        ndrange = nothing
    end

    return ndrange, workgroupsize, iterspace, dynamic
end

# Inference barriers
function __run(obj, ndrange, iterspace, args, dynamic)
    N = length(iterspace)
    Nthreads = Threads.nthreads()
    if Nthreads == 1
        len, rem = N, 0
    else
        len, rem = divrem(N, Nthreads)
    end
    # not enough iterations for all the threads?
    if len == 0
        Nthreads = N
        len, rem = 1, 0
    end
    if Nthreads == 1
        __thread_run(1, len, rem, obj, ndrange, iterspace, args, dynamic)
    else
        @sync for tid in 1:Nthreads
            Threads.@spawn __thread_run(tid, len, rem, obj, ndrange, iterspace, args, dynamic)
        end
    end
    return nothing
end

function __thread_run(tid, len, rem, obj, ndrange, iterspace, args, dynamic)
    # compute this thread's iterations
    f = 1 + ((tid-1) * len)
    l = f + len - 1
    # distribute remaining iterations evenly
    if rem > 0
        if tid <= rem
            f = f + (tid-1)
            l = l + tid
        else
            f = f + rem
            l = l + rem
        end
    end
    # run this thread's iterations
    for i = f:l
        block = @inbounds blocks(iterspace)[i]
        ctx = mkcontext(obj, block, ndrange, iterspace, dynamic)
        obj.f(ctx, args...)
    end
    return nothing
end

function mkcontext(kernel::Kernel{CPU}, I, _ndrange, iterspace, ::Dynamic) where Dynamic
    return CompilerMetadata{ndrange(kernel), Dynamic}(I, _ndrange, iterspace)
end

@inline function __index_Local_Linear(ctx, idx::CartesianIndex)
    indices = workitems(__iterspace(ctx))
    return @inbounds LinearIndices(indices)[idx]
end

@inline function __index_Group_Linear(ctx, idx::CartesianIndex)
    indices = blocks(__iterspace(ctx))
    return @inbounds LinearIndices(indices)[__groupindex(ctx)]
end

@inline function __index_Global_Linear(ctx, idx::CartesianIndex)
    I = @inbounds expand(__iterspace(ctx), __groupindex(ctx), idx)
    @inbounds LinearIndices(__ndrange(ctx))[I]
end

@inline function __index_Local_Cartesian(_, idx::CartesianIndex)
    return idx
end

@inline function __index_Group_Cartesian(ctx, ::CartesianIndex)
    __groupindex(ctx)
end

@inline function __index_Global_Cartesian(ctx, idx::CartesianIndex)
    return @inbounds expand(__iterspace(ctx), __groupindex(ctx), idx)
end

@inline function __validindex(ctx, idx::CartesianIndex)
    # Turns this into a noop for code where we can turn of checkbounds of
    if __dynamic_checkbounds(ctx)
        I = @inbounds expand(__iterspace(ctx), __groupindex(ctx), idx)
        return I in __ndrange(ctx)
    else
        return true
    end
end

###
# CPU implementation of shared memory
###
@inline function SharedMemory(::Type{T}, ::Val{Dims}, ::Val) where {T, Dims}
    MArray{__size(Dims), T}(undef)
end

###
# CPU implementation of scratch memory
# - memory allocated as a MArray with size `Dims`
###

struct ScratchArray{N, D}
    data::D
    ScratchArray{N}(data::D) where {N, D} = new{N, D}(data)
end

@inline function Scratchpad(ctx, ::Type{T}, ::Val{Dims}) where {T, Dims}
    return ScratchArray{length(Dims)}(MArray{__size((Dims..., prod(__groupsize(ctx)))), T}(undef))
end

# Base.view creates a boundscheck which captures A
# https://github.com/JuliaLang/julia/issues/39308
@inline function aview(A, I::Vararg{Any, N}) where N
     J = Base.to_indices(A, I)
     Base.unsafe_view(Base._maybe_reshape_parent(A, Base.index_ndims(J...)), J...)
end

@inline function Base.getindex(A::ScratchArray{N}, idx) where N
    return @inbounds aview(A.data, ntuple(_->:, Val(N))..., idx)
end

# Argument conversion
KernelAbstractions.argconvert(k::Kernel{CPU}, arg) = arg


text:
# wrapper for API requests.
function wmata_request(url::String)
    subscription_key = Dict("api_key" => wmata.api_key)
    api_response = parse(String(request("GET", url, subscription_key).body))

    return api_response
end

# verify that a linecode is one that makes sense
function verify_line_input(line_input)
    line_colors = ["RD", "BL", "YL", "OR", "GR", "SV", "All"]
    if !(line_input in line_colors)
        error("LineCode must be one of: RD, BL, YL, OR, GR, SV, or All")
    else
        return line_input
    end
end

# verify that a station code has a match in the WMATA endpoint.
function verify_station_input(station_input)
    # this function relies on the same endpoint that the station_list function does.
    r = wmata_request(wmata.station_list_url)

    valid_station_codes = push!([station["Code"] for station in r["Stations"]], "All")

    if !(station_input in valid_station_codes)
        error("$station_input is not a valid station code.\nTry using station_list to find and verify your station code.")
    else
        return station_input
    end
end

#=
support optional argument in functions that involve pulling details
 based on a station code - enables a user to use a station name if they
 don't know the station code.
=#
function get_station_code(StationName::String)
    r = wmata_request(wmata.station_list_url)

    stations = Dict(
        [station["Name"] for station in r["Stations"]] .=> [station["Code"] for station in r["Stations"]]
    )

    if !(StationName in keys(stations))
        error("$StationName is not a valid station name.")
    else
        return stations[StationName]
    end
end

#=
it makes sense that we'd want to parse arrival times to a better format (aka one that we can do arithmetic on),
  there is probably a more Juli-onic way to do this but this is my best attempt at the moment.
=#
function convert_arrival_times(arrival_times::Vector{String})
    converted_times = []
    for time in arrival_times
        if time == "ARR" || time == "BRD"
            push!(converted_times, 0)
        elseif time == "---"
            push!(converted_times, missing)
        else
            push!(converted_times, Base.parse(Int64, time))
        end
    end
    return convert(Vector{Union{Missing, Int64}}, converted_times)
end


text:
# Conjugates for normal distribution

using Base.Test
using Distributions
using ConjugatePriors

import ConjugatePriors: NormalGamma, NormalInverseGamma
import ConjugatePriors: posterior, posterior_rand, posterior_mode, posterior_randmodel, fit_map

n = 100
w = rand(100)

# Νormal - Νormal (known sigma)

pri = Normal(1.0, 5.0)

x = rand(Normal(2.0, 3.0), n)
p = posterior((pri, 3.0), Normal, x)
@test isa(p, Normal)
@test_approx_eq mean(p)  (mean(pri) / var(pri) + sum(x) / 9.0) / (1.0 / var(pri) + n / 9.0)
@test_approx_eq var(p) inv(1.0 / var(pri) + n / 9.0)

r = posterior_mode((pri, 3.0), Normal, x)
@test_approx_eq r mode(p)

f = fit_map((pri, 3.0), Normal, x)
@test isa(f, Normal)
@test f.μ == r
@test f.σ == 3.0

p = posterior((pri, 3.0), Normal, x, w)
@test isa(p, Normal)
@test_approx_eq mean(p)  (mean(pri) / var(pri) + dot(x, w) / 9.0) / (1.0 / var(pri) + sum(w) / 9.0)
@test_approx_eq var(p) inv(1.0 / var(pri) + sum(w) / 9.0)

r = posterior_mode((pri, 3.0), Normal, x, w)
@test_approx_eq r mode(p)

f = fit_map((pri, 3.0), Normal, x, w)
@test isa(f, Normal)
@test f.μ == r
@test f.σ == 3.0

# ΙnverseGamma - Νormal (known mu)

pri = InverseGamma(1.5, 0.5)

x = rand(Normal(2.0, 3.0), n)
p = posterior((2.0, pri), Normal, x)
@test isa(p, InverseGamma)
@test_approx_eq shape(p) shape(pri) + n / 2
@test_approx_eq scale(p) scale(pri) + sum(abs2(x .- 2.0)) / 2

r = posterior_mode((2.0, pri), Normal, x)
@test_approx_eq r mode(p)

f = fit_map((2.0, pri), Normal, x)
@test isa(f, Normal)
@test f.μ == 2.0
@test_approx_eq abs2(f.σ) r

p = posterior((2.0, pri), Normal, x, w)
@test isa(p, InverseGamma)
@test_approx_eq shape(p) shape(pri) + sum(w) / 2
@test_approx_eq scale(p) scale(pri) + dot(w, abs2(x .- 2.0)) / 2

r = posterior_mode((2.0, pri), Normal, x, w)
@test_approx_eq r mode(p)

f = fit_map((2.0, pri), Normal, x, w)
@test isa(f, Normal)
@test f.μ == 2.0
@test_approx_eq abs2(f.σ) r

# Gamma - Νormal (known mu)

pri = Gamma(1.5, 2.0)

x = rand(Normal(2.0, 3.0), n)
p = posterior((2.0, pri), Normal, x)
@test isa(p, Gamma)
@test_approx_eq shape(p) shape(pri) + n / 2
@test_approx_eq scale(p) scale(pri) + sum(abs2(x .- 2.0)) / 2

r = posterior_mode((2.0, pri), Normal, x)
@test_approx_eq r mode(p)

f = fit_map((2.0, pri), Normal, x)
@test isa(f, Normal)
@test f.μ == 2.0
@test_approx_eq abs2(f.σ) inv(r)

p = posterior((2.0, pri), Normal, x, w)
@test isa(p, Gamma)
@test_approx_eq shape(p) shape(pri) + sum(w) / 2
@test_approx_eq scale(p) scale(pri) + dot(w, abs2(x .- 2.0)) / 2

r = posterior_mode((2.0, pri), Normal, x, w)
@test_approx_eq r mode(p)

f = fit_map((2.0, pri), Normal, x, w)
@test isa(f, Normal)
@test f.μ == 2.0
@test_approx_eq abs2(f.σ) inv(r)

# NormalInverseGamma - Normal

mu_true = 2.
sig2_true = 3.
x = rand(Normal(mu_true, sig2_true), n)

mu0 = 2.
v0 = 3.
shape0 = 5.
scale0 = 2.
pri = NormalInverseGamma(mu0, v0, shape0, scale0)

post = posterior(pri, Normal, x)
@test isa(post, NormalInverseGamma)

@test_approx_eq post.mu (mu0/v0 + n*mean(x))/(1./v0 + n)
@test_approx_eq post.v0 1./(1./v0 + n)
@test_approx_eq post.shape shape0 + 0.5*n
@test_approx_eq post.scale scale0 + 0.5*(n-1)*var(x) + n./v0./(n + 1./v0)*0.5*(mean(x)-mu0).^2

ps = posterior_randmodel(pri, Normal, x)

@test isa(ps, Normal)
@test insupport(ps,ps.μ) && ps.σ > zero(ps.σ)

# NormalGamma - Normal

mu_true = 2.
tau2_true = 3.
x = rand(Normal(mu_true, 1./tau2_true), n)

mu0 = 2.
nu0 = 3.
shape0 = 5.
rate0 = 2.
pri = NormalGamma(mu0, nu0, shape0, rate0)

post = posterior(pri, Normal, x)
@test isa(post, NormalGamma)

@test_approx_eq post.mu (nu0*mu0 + n*mean(x))./(nu0 + n)
@test_approx_eq post.nu nu0 + n
@test_approx_eq post.shape shape0 + 0.5*n
@test_approx_eq post.rate rate0 + 0.5*(n-1)*var(x) + n*nu0/(n + nu0)*0.5*(mean(x)-mu0).^2

ps = posterior_randmodel(pri, Normal, x)

@test isa(ps, Normal)
@test insupport(ps, ps.μ) && ps.σ > zero(ps.σ)


text:
module RNN
export beer_95

using Parameters: @unpack
using ..ExampleBase: LEDemo, ContinuousExample

σ(x) = 1 / (1 + exp(-x))
σ′(x) = σ(x) * (1 - σ(x))

struct ContinuousRNN{W,TH,TA}
    w::W
    θ::TH
    τ::TA
end

@inline function phase_dynamics!(du, u, p, t)
    @unpack w, θ, τ = p
    du .= (- u .+ w * σ.(u .+ θ)) ./ τ
end

@inline @views function tangent_dynamics!(du, u, p, t)
    @unpack w, θ, τ = p
    phase_dynamics!(du[:, 1], u[:, 1], p, t)
    Y = u[:, 2:end]
    du[:, 2:end] .= (- Y .+ w * Diagonal(σ′.(u[:, 1] .+ θ)) * Y) ./ τ
end

"""
Return a [`LEDemo`](@ref) for a low-dimensional chaotic
continuous-time recurrent neural networks by Beer (1995).

* Beer, R. D. (1995). On the dynamics of small continuous-time recurrent
  neural networks. Adapt. Behav., 3(4), 469–509.
  <https://doi.org/10.1177/105971239500300405>.
  (Figure 9D)
"""
function beer_95(;
        u0=[0.1, 0.1, 0.1],
        t_renorm=10.0,
        t_attr = t_renorm * 4000,
        atol=1e-2, rtol=5e-1,
        terminator_options = [:atol => atol, :rtol => rtol,
                              :max_tail_corr => 0.2,
                              # Use max_mle smaller than known_exponents[1]:
                              :max_mle => 0.001],
        kwargs...)
    param = ContinuousRNN(
        # w
        [  5.422  -0.018  2.75
          -0.24    4.59   1.21
           0.535  -2.25   3.885 ],
        # θ
        [-4.108, -2.787, -1.114],
        # τ
        [1.0, 2.5, 1.0],
    )
    LEDemo(ContinuousExample(
        "Beer (1995)",
        phase_dynamics!, u0, t_renorm, param,
        tangent_dynamics!,
        t_attr,
        [0.010, 0],   # known_exponents
        atol, rtol,
        terminator_options,
    ); kwargs...)
end

end


text:
using RHEOS; using PyPlot

data = importcsv("../examples/example1_data.csv", t_col=1, ϵ_col=2, σ_col=3)

maxwell_model = modelfit(data, Maxwell, strain_imposed)
maxwell_predict = modelpredict(extract(data, strain_only), maxwell_model)

fig, ax = subplots()

ax.plot(data.t, data.ϵ, "--", label="Strain (Data)")
ax.plot(data.t, data.σ, ".", label="Stress (Data)")
ax.plot(maxwell_predict.t, maxwell_predict.σ, label="Stress (Predicted)")

ax.set_xlabel("Time (s)")
ax.set_ylabel("Stress, Strain")
ax.legend()


text:
struct IdealManifold{N, T <: SMatrix{3,N}} <: AbstractManifold{N}
    scaled_antenna_positions::T
end

"""
Calculates the manifold based on antenna position and signal frequency.

# Examples
```julia-repl
julia> manifold = IdealManifold(1575420e3, 0.1904 / 4 * SVector(SVector(1, 1, 0), SVector(-1, 1, 0), SVector(1, -1, 0), SVector(-1, -1, 0)))
julia> get_steer_vec(manifold, SVector(0.0,0.0,1.0))
julia> get_steer_vec(manifold, SVector(0.0,0.0,1.0), RotXYZ(0.0,0.0,0.0))
```
"""
function IdealManifold(f_0::Real, antenna_positions::SMatrix{3,N,T}; c0 = 299_792_458) where {N, T <: Real}
    λ = c0 / f_0
    IdealManifold(2π / λ * antenna_positions)
end

function IdealManifold(f_0::Real, antenna_positions::SVector{N,SVector{3,T}}; c0 = 299_792_458) where {N, T <: Real}
    IdealManifold(f_0, reduce(hcat, antenna_positions), c0 = c0)
end

function IdealManifold(f_0::Real, antenna_positions::SVector{3,T}...; c0 = 299_792_458) where T <: Real
    IdealManifold(f_0, SVector(antenna_positions), c0 = c0)
end

function IdealManifold(f_0::Real, num_ant_x, num_ant_y, dist_x, dist_y; kwargs...)
    max_x = (num_ant_x - 1) / 2
    max_y = (num_ant_y - 1) / 2
    pos_x = (-max_x:max_x) .* dist_x
    pos_y = (-max_y:max_y) .* dist_y
    ant_pos = SMatrix{3, num_ant_x * num_ant_y}(
        vcat(
            repeat(pos_x, outer = num_ant_y)',
            repeat(pos_y, inner = num_ant_x)',
            zeros(1, num_ant_x * num_ant_y)
        )
    )
    IdealManifold(f_0, ant_pos; kwargs...)
end

function IdealManifold(f_0::Real, num_ant_xy, dist_xy; kwargs...)
    IdealManifold(f_0, num_ant_xy, num_ant_xy, dist_xy, dist_xy; kwargs...)
end

function get_steer_vec(manifold::IdealManifold, doa)
    cis.(transpose(manifold.scaled_antenna_positions) * doa)
end

function IdealManifold()
    IdealManifold(zeros(SMatrix{3, 1, Float64}))
end

function get_steer_vec(manifold::IdealManifold{1}, doa)
    1.0
end

function get_steer_vec(manifold::IdealManifold{1}, doa, attitude)
    1.0
end

function get_steer_vec(manifold::IdealManifold, doa::Spherical)
    get_steer_vec(manifold, sph2cart(doa))
end


text:
using ProximalOperators
using ProximalAlgorithms
using LinearAlgebra

function assert_lp_solution(c, A, b, x, y, tol)
    # Check and print solution quality measures (for some reason the
    # returned dual iterate is the negative of the dual LP variable y)

    nonneg = -minimum(min.(0.0, x))
    @test nonneg <= tol

    primal_feasibility = norm(A * x - b)
    @test primal_feasibility <= tol

    dual_feasibility = maximum(max.(0.0, -A' * y - c))
    @test dual_feasibility <= tol

    complementarity = abs(dot(c + A' * y, x))
    @test complementarity <= tol
end

@testset "Linear programs ($T)" for T in [Float32, Float64]

    # Solving LP with AFBA
    #
    #   minimize    c'x             -> f = <c,.>
    #   subject to  Ax = b          -> h = ind{b}, L = A
    #               x >= 0          -> g = ind{>=0}
    #
    # Dual LP
    #
    #   maximize    b'y
    #   subject to  A'y <= c
    #
    # Optimality conditions
    #
    #   x >= 0              [primal feasibility 1]
    #   Ax = b              [primal feasibility 2]
    #   A'y <= c            [dual feasibility]
    #   X'(c - A'y) = 0     [complementarity slackness]
    #

    n = 10 # primal dimension
    m = 8 # dual dimension (i.e. number of linear equalities)
    k = 5 # number of active dual constraints (must be 0 <= k <= n)
    # primal optimal point
    x_star = T[0.03606099647643202, 0.6641306619990367, 0.14689326835593086, 0.8616058527226432, 0.6518888386753204, 0.0, 0.0, 0.0, 0.0, 0.0]
    # dual optimal slack variable
    s_star = T[0.0, 0.0, 0.0, 0.0, 0.0, 0.7652088547602296, 0.5474850582278736, 0.6291635905927829, 0.5505791259743739, 0.2787959059899776]
    # dual optimal point
    y_star = T[-0.18584225497010673, 1.0289179383971105, -0.5464276767930338, 0.6916826379378143, 0.0052105118567898154, -1.3874339670318891, 2.6448648322602337, 0.4436510755789434]
    # equality constraints coefficients
    A = T[
        0.9670566934980286 0.3533472582831891 -0.2819205403663176 -0.3071751292615649 -0.008197151577028876 1.9968628487924958 -1.3978028828102382 -0.10109443461339453 0.3540170332321368 -0.5983697948494805;
        -0.08643002247762153 1.0822796020885652 -0.8133412010128359 0.5073579509595517 -0.2793798574957181 -0.3437059320143032 -0.10596751644319548 0.8296541102523696 -0.28406001552384064 0.33566140783820164;
        -0.6474810505335304 1.2524424221527595 0.7261766049639965 0.6930548839308885 0.9064585656428545 0.5197152355759463 1.3746077918877961 1.5289983684262054 0.9258506062644877 -1.3774011601531342;
        0.8725430405646543 1.1911309824177332 -0.7221185305116811 0.44324697589700257 -0.15551892138880116 -0.5167033349322372 -1.4774045165687548 -0.8670756465083638 -1.4024841888738206 0.296570674868075;
        0.5252662593900623 -0.6817132148621747 0.37446015899900237 0.4180282147408251 -0.8405436435394317 -1.742941478155391 -0.06419427972299957 0.43910205564784205 1.0643883425210827 2.3063869854427335;
        -0.7148913270640012 -0.769028546306448 -1.059257097999333 -0.6795170119545777 -0.0498976886779524 -1.2392107698826862 -0.4415384005606088 -1.058758868936871 -0.108504245219676 -0.5576550366602419;
        -1.208388774142618 -0.15206129387542855 2.311520055340236 0.8043266793420988 -0.5692874893454578 1.1246423711381501 0.5335942753441769 2.6595405998250876 0.09162292399585106 0.3749905036072034;
        -0.2300660921924555 -0.7014271654627467 -0.20170532145095504 -0.02503019691724233 -2.2191605023268512 0.9110653907470295 1.8238644805628141 -1.1908921287611471 0.12168786553115268 0.17399181994853638
    ]
    b = A * x_star
    c = A' * y_star + s_star

    tol = 100 * eps(T)
    maxit = 100_000

    @testset "AFBA" begin

        f = Linear(c)
        g = IndNonnegative()
        h = IndPoint(b)

        x0 = zeros(T, n)
        x0_backup = copy(x0)
        y0 = zeros(T, m)
        y0_backup = copy(y0)

        solver = ProximalAlgorithms.AFBA(tol = tol, maxit = maxit)
        (x, y), it = solver(x0 = x0, y0 = y0, f = f, g = g, h = h, L = A, beta_f = 0)

        @test eltype(x) == T
        @test eltype(y) == T

        @test it <= maxit

        assert_lp_solution(c, A, b, x, y, 1000 * tol)

        @test x0 == x0_backup
        @test y0 == y0_backup

    end

    @testset "VuCondat" begin

        f = Linear(c)
        g = IndNonnegative()
        h = IndPoint(b)

        x0 = zeros(T, n)
        x0_backup = copy(x0)
        y0 = zeros(T, m)
        y0_backup = copy(y0)

        solver = ProximalAlgorithms.VuCondat(tol = tol, maxit = maxit)
        (x, y), it = solver(x0 = x0, y0 = y0, f = f, g = g, h = h, L = A, beta_f = 0)

        @test eltype(x) == T
        @test eltype(y) == T

        @test it <= maxit

        assert_lp_solution(c, A, b, x, y, 1000 * tol)

        @test x0 == x0_backup
        @test y0 == y0_backup

    end

    @testset "ChambollePock" begin
        g = Linear(c)
        h = SlicedSeparableSum((IndPoint(b), IndNonnegative()), ((1:m,), (m+1:m+n,)))

        x0 = zeros(T, n)
        x0_backup = copy(x0)
        y0 = zeros(T, m + n)
        y0_backup = copy(y0)

        solver = ProximalAlgorithms.ChambollePock(tol = tol, maxit = maxit)
        (x, y), it = solver(x0 = x0, y0 = y0, g = g, h = h, L = vcat(A, Matrix{T}(I, n, n)))

        @test eltype(x) == T
        @test eltype(y) == T

        @test it <= maxit

        assert_lp_solution(c, A, b, x, y[1:m], 1000 * tol)

        @test x0 == x0_backup
        @test y0 == y0_backup
    end

    @testset "DavisYin" begin

        f = Linear(c)
        g = IndNonnegative()
        h = IndAffine(A, b)

        x0 = zeros(T, n)
        x0_backup = copy(x0)

        solver = ProximalAlgorithms.DavisYin(gamma = T(1), tol = tol, maxit = maxit)
        xf, it = solver(x0 = x0, f = f, g = g, h = h)

        @test eltype(xf) == T

        @test it <= maxit

        @assert norm(xf - x_star) <= 1e2 * tol

        @test x0 == x0_backup

    end

end


text:
function get_tspn05()

m = Model()

# ----- Variables ----- #
@variable(m, objvar)
x_Idx = Any[1, 2, 3, 4, 5, 6, 7, 8, 9, 10]
@variable(m, x[x_Idx])
b_Idx = Any[11, 12, 13, 14, 15, 16, 17, 18, 19, 20]
@variable(m, b[b_Idx], Bin)
JuMP.set_lower_bound(x[1], 68.0)
JuMP.set_upper_bound(x[1], 71.0)
JuMP.set_lower_bound(x[2], 65.0)
JuMP.set_upper_bound(x[2], 87.0)
JuMP.set_lower_bound(x[3], 107.0)
JuMP.set_upper_bound(x[3], 126.0)
JuMP.set_lower_bound(x[4], 38.0)
JuMP.set_upper_bound(x[4], 49.0)
JuMP.set_lower_bound(x[5], 40.0)
JuMP.set_upper_bound(x[5], 55.0)
JuMP.set_lower_bound(x[6], 54.0)
JuMP.set_upper_bound(x[6], 68.0)
JuMP.set_lower_bound(x[7], 92.0)
JuMP.set_upper_bound(x[7], 106.0)
JuMP.set_lower_bound(x[8], 113.0)
JuMP.set_upper_bound(x[8], 117.0)
JuMP.set_lower_bound(x[9], 82.0)
JuMP.set_upper_bound(x[9], 87.0)
JuMP.set_lower_bound(x[10], 76.0)
JuMP.set_upper_bound(x[10], 85.0)

# ----- Constraints ----- #
@NLconstraint(m, e1, sqrt( (x[1]-x[3])^2+ (x[2]-x[4])^2)*b[11]+sqrt( (x[1]-x[5])^2+ (x[2]-x[6])^2)*b[12]+sqrt( (x[1]-x[7])^2+ (x[2]-x[8])^2)*b[13]+sqrt( (x[1]-x[9])^2+ (x[2]-x[10])^2)*b[14]+sqrt( (x[3]-x[5])^2+ (x[4]-x[6])^2)*b[15]+sqrt( (x[3]-x[7])^2+ (x[4]-x[8])^2)*b[16]+sqrt( (x[3]-x[9])^2+ (x[4]-x[10])^2)*b[17]+sqrt( (x[5]-x[7])^2+ (x[6]-x[8])^2)*b[18]+sqrt( (x[5]-x[9])^2+ (x[6]-x[10])^2)*b[19]+sqrt( (x[7]-x[9])^2+ (x[8]-x[10])^2)*b[20]-objvar == 0.0)
@NLconstraint(m, e2, 0.444444444444444* (x[1])^2-61.7777777777778*x[1]+0.00826446280991736* (x[2])^2-1.25619834710744*x[2] <= -2193.51331496786)
@NLconstraint(m, e3, 0.0110803324099723* (x[3])^2-2.58171745152355*x[3]+0.0330578512396694* (x[4])^2-2.87603305785124*x[4] <= -211.938760559511)
@NLconstraint(m, e4, 0.0177777777777778* (x[5])^2-1.68888888888889*x[5]+0.0204081632653061* (x[6])^2-2.48979591836735*x[6] <= -115.049886621315)
@NLconstraint(m, e5, 0.0204081632653061* (x[7])^2-4.04081632653061*x[7]+0.25* (x[8])^2-57.5*x[8] <= -3505.27040816327)
@NLconstraint(m, e6, 0.16* (x[9])^2-27.04*x[9]+0.0493827160493827* (x[10])^2-7.95061728395062*x[10] <= -1461.45234567901)
@constraint(m, e7, b[11]+b[12]+b[13]+b[14] == 2.0)
@constraint(m, e8, b[11]+b[15]+b[16]+b[17] == 2.0)
@constraint(m, e9, b[12]+b[15]+b[18]+b[19] == 2.0)
@constraint(m, e10, b[13]+b[16]+b[18]+b[20] == 2.0)
@constraint(m, e11, b[14]+b[17]+b[19]+b[20] == 2.0)

# ----- Objective ----- #
@objective(m, Min, objvar)

return m
end


text:
using BufferedStreams
using Test

struct InfiniteStream <: IO
    byte::UInt8
end

function Base.eof(stream::InfiniteStream)
    return false
end

function Base.read(stream::InfiniteStream, ::Type{UInt8})
    return stream.byte
end

if isdefined(Base, :unsafe_read)
    function Base.unsafe_read(stream::InfiniteStream, pointer::Ptr{UInt8}, n::UInt)
        ccall(:memset, Cvoid, (Ptr{Cvoid}, Cint, Csize_t), pointer, stream.byte, n)
        return nothing
    end
end

# A few things that might not be obvious:
#   * In a few places we wrap an array in an IOBuffer before wrapping it in a
#     BufferedInputStream. This is to force the BufferedInputStream to read
#     incrementally to expose possible bugs in buffer refilling logic.
#   * Similar, we manually set the buffer size to be smaller than the default to
#     force more buffer refills.

@testset "BufferedInputStream" begin
    @testset "read" begin
        data = rand(UInt8, 1000000)
        stream = BufferedInputStream(IOBuffer(data), 1024)
        read_data = UInt8[]
        while !eof(stream)
            push!(read_data, read(stream, UInt8))
        end
        @test data == read_data
        @test data == read(BufferedInputStream(IOBuffer(data), 1024))

        halfn = div(length(data), 2)
        @test data[1:halfn] == read(BufferedInputStream(IOBuffer(data), 1024), halfn)

        # read multi-byte data
        buffer = IOBuffer()
        write(buffer, UInt8(1))
        write(buffer, UInt16(2))
        write(buffer, UInt32(3))
        write(buffer, UInt64(4))
        write(buffer, UInt128(5))
        write(buffer, Float16(6))
        write(buffer, Float32(7))
        write(buffer, Float64(8))
        seekstart(buffer)
        stream = BufferedInputStream(buffer)
        @test read(stream, UInt8) === UInt8(1)
        @test read(stream, UInt16) === UInt16(2)
        @test read(stream, UInt32) === UInt32(3)
        @test read(stream, UInt64) === UInt64(4)
        @test read(stream, UInt128) === UInt128(5)
        @test read(stream, Float16) === Float16(6)
        @test read(stream, Float32) === Float32(7)
        @test read(stream, Float64) === Float64(8)

        # EOFError
        stream = BufferedInputStream(IOBuffer())
        @test_throws EOFError read(stream, UInt8)
        @test_throws EOFError read(stream, UInt16)
        @test_throws EOFError read(stream, UInt32)
        @test_throws EOFError read(stream, UInt64)
        @test_throws EOFError read(stream, UInt128)
        @test_throws EOFError read(stream, Float16)
        @test_throws EOFError read(stream, Float32)
        @test_throws EOFError read(stream, Float64)
    end

    if isdefined(Base, :unsafe_read)
        @testset "unsafe_read" begin
            stream = BufferedInputStream(IOBuffer("abcdefg"), 3)
            data = Vector{UInt8}(undef, 7)
            unsafe_read(stream, pointer(data, 1), 1)
            @test data[1] == UInt8('a')
            unsafe_read(stream, pointer(data, 2), 2)
            unsafe_read(stream, pointer(data, 4), 4)
            @test data == b"abcdefg"
            @test_throws EOFError unsafe_read(stream, pointer(data), 1)
        end
    end

    @testset "peek" begin
        stream = BufferedInputStream(IOBuffer([0x01, 0x02]))
        @test peek(stream) === 0x01
        @test peek(stream) === 0x01
        read(stream, UInt8)
        @test peek(stream) === 0x02
        @test peek(stream) === 0x02
    end

    @testset "bytesavailable" begin
        stream = BufferedInputStream(IOBuffer([0x01, 0x02]))
        @test bytesavailable(stream) == 2
        read(stream, 1)
        @test bytesavailable(stream) == 1
        read(stream, 1)
        @test bytesavailable(stream) == 0
    end

    @testset "peekbytes!" begin
        data = rand(UInt8, 1000000)
        stream = BufferedInputStream(IOBuffer(data), 1024)

        read_data = Array{UInt8}(undef, 1000)
        @test peekbytes!(stream, read_data, 1000) == 1000
        @test data[1:1000] == read_data
        # Check that we read the bytes we just peeked, i.e. that the position
        # wasn't advanced on peekbytes!()
        readbytes!(stream, read_data, 1000)
        @test data[1:1000] == read_data
        # Check that we now peek the next 24 bytes, as we don't go past the end
        # of the buffer
        @test peekbytes!(stream, read_data, 1000) == 24
        @test data[1001:1024] == read_data[1:24]

        # Reset the data
        data = rand(UInt8, 1000000)
        stream = BufferedInputStream(IOBuffer(data), 1024)

        read_data = Array{UInt8}(undef, 2000)
        @test peekbytes!(stream, read_data, 2000) == 1024
        # Note that we truncate at buffer size, as
        @test data[1:1024] == read_data[1:1024]

        # Check that we only read up to the buffer size
        read_data = Array{UInt8}(undef, 5)
        @test peekbytes!(stream, read_data) == 5
        @test data[1:5] == read_data

        close(stream)
        @test_throws Exception peekbytes!(stream, read_data)
    end

    @testset "readbytes!" begin
        stream = BufferedInputStream(IOBuffer([0x01:0xff;]), 4)
        @test !eof(stream)
        out = zeros(UInt8, 2)
        @test BufferedStreams.readbytes!(stream, out) === 2
        @test out == [0x01, 0x02]
        out = zeros(UInt8, 3)
        @test BufferedStreams.readbytes!(stream, out) === 3
        @test out == [0x03, 0x04, 0x05]
        out = zeros(UInt8, 5)
        @test BufferedStreams.readbytes!(stream, out) === 5
        @test out == [0x06, 0x07, 0x08, 0x09, 0x0a]
        @test !eof(stream)
        @test read(stream) == [0x0b:0xff;]
        @test eof(stream)
    end

    @testset "readuntil" begin
        data = rand(UInt8, 1000000)
        stream = BufferedInputStream(IOBuffer(data), 1024)

        true_num_zeros = 0
        zero_positions = Int[]
        for (i, b) in enumerate(data)
            if b == 0x00
                push!(zero_positions, i)
                true_num_zeros += 1
            end
        end

        num_zeros = 0
        chunk_results = Bool[]
        while true
            # are we extracting the right chunk?
            chunk = readuntil(stream, 0x00)
            first = num_zeros == 0 ? 1 : zero_positions[num_zeros]+1
            last = num_zeros < length(zero_positions) ?  zero_positions[num_zeros+1] : length(data)
            true_chunk = data[first:last]
            push!(chunk_results, true_chunk == chunk)

            if !eof(stream)
                num_zeros += 1
            else
                break
            end
        end

        @test all(chunk_results)
        @test num_zeros == true_num_zeros
    end

    @testset "arrays" begin
        data = rand(UInt8, 1000000)
        stream = BufferedInputStream(data)
        read_data = UInt8[]
        while !eof(stream)
            push!(read_data, read(stream, UInt8))
        end
        @test data == read_data
        @test data == read(BufferedInputStream(data))
    end

    @testset "marks" begin
        # very small buffer
        stream = BufferedInputStream(IOBuffer([0x01:0xff;]), 2)
        @test !ismarked(stream)
        mark(stream)
        @test ismarked(stream)
        a = read(stream, UInt8)
        b = read(stream, UInt8)
        c = read(stream, UInt8)
        reset(stream)
        @test !ismarked(stream)
        a′ = read(stream, UInt8)
        b′ = read(stream, UInt8)
        c′ = read(stream, UInt8)
        @test (a, b, c) == (a′, b′, c′) == (0x01, 0x02, 0x03)
        mark(stream)
        @test unmark(stream)
        @test !unmark(stream)
        @test_throws ErrorException reset(stream)
    end

    @testset "anchors" begin
        data = rand(UInt8, 100000)

        function random_range(n)
            a = rand(1:n)
            b = rand(a:n)
            return a:b
        end

        # test that anchors work correctly by extracting random intervals from a
        # buffered stream.
        function test_anchor()
            r = random_range(length(data))
            i = 1
            stream = BufferedInputStream(IOBuffer(data), 1024)
            while !eof(stream)
                if i == r.start
                    anchor!(stream)
                end
                if i == r.stop + 1
                    return takeanchored!(stream) == data[r]
                end
                read(stream, UInt8)
                i += 1
            end
            if i == r.stop + 1
                return takeanchored!(stream) == data[r]
            end
            error("nothing extracted (range: $r)")
        end

        @test all(Bool[test_anchor() for _ in 1:100])
    end

    @testset "seek" begin
        n = 100000
        data = rand(UInt8, n)
        positions = rand(0:n-1, 1000)

        function test_seek(stream, p)
            seek(stream, p)
            return position(stream) == p && read(stream, UInt8) == data[p + 1]
        end

        stream = BufferedInputStream(data)
        @test all(Bool[test_seek(stream, p) for p in positions])

        stream = BufferedInputStream(IOBuffer(data), 1024)
        @test all(Bool[test_seek(stream, p) for p in positions])
    end

    @testset "skip" begin
        n = 100000
        data = rand(UInt8, n)
        positions = rand(0:n-1, 1000)
        sort!(positions)

        last = 1
        offsets = Int[]
        for p in positions
            push!(offsets, p - last)
            last = p
        end

        function test_skip(stream, p, offset)
            skip(stream, offset)
            peek(stream) == data[p]
        end

        stream = BufferedInputStream(IOBuffer(data), 1024)
        @test_throws Exception skip(stream, n + 1)
        @test_throws Exception skip(stream, -1)

        stream = BufferedInputStream(IOBuffer(data), 1024)
        @test all(Bool[test_skip(stream, p, offset)
                       for (p, offset) in zip(positions, offsets)])
    end

    @testset "close" begin
        iobuffer = IOBuffer([0x00, 0x01])
        stream = BufferedInputStream(iobuffer)
        @test isopen(stream)
        @test isopen(iobuffer)
        read(stream, UInt8)
        @test close(stream) === nothing
        @test !isopen(stream)
        @test !isopen(iobuffer)
        @test_throws Exception read(stream, UInt8)
        @test close(stream) === nothing
    end

    @testset "iostream" begin
        mktemp() do path, input
            write(input, [0x01, 0x02, 0x03, 0x04, 0x05])
            flush(input)
            seekstart(input)

            stream = BufferedInputStream(input, 2)
            @test !eof(stream)
            @test read(stream, UInt8) === 0x01
            @test !eof(stream)
            @test read(stream, UInt8) === 0x02
            @test !eof(stream)
            @test read(stream) == [0x03, 0x04, 0x05]
            @test eof(stream)

            @test isopen(stream)
            @test isopen(input)
            close(stream)
            @test !isopen(stream)
            @test !isopen(input)
        end
    end

    @testset "immobilized buffer" begin
        stream = BufferedInputStream(IOBuffer("abcdefg"), 2)
        stream.immobilized = false
        @assert read(stream, UInt8) == UInt8('a')
        mark(stream)
        @assert read(stream, UInt8) == UInt8('b')
        BufferedStreams.fillbuffer!(stream)
        @test stream.buffer[1] == UInt8('b')

        stream = BufferedInputStream(IOBuffer("abcdefg"), 2)
        stream.immobilized = true
        @assert read(stream, UInt8) == UInt8('a')
        mark(stream)
        @assert read(stream, UInt8) == UInt8('b')
        BufferedStreams.fillbuffer!(stream)
        @test stream.buffer[2] == UInt8('b')

        stream = BufferedInputStream(IOBuffer("abcdefg"), 6)
        stream.immobilized = true
        data = Vector{UInt8}(undef, 7)
        BufferedStreams.readbytes!(stream, data, 1, 3)
        @test data[1:3] == b"abc"
        BufferedStreams.readbytes!(stream, data, 4, 7)
        @test data[4:7] == b"defg"
    end

    @testset "shiftdata!" begin
        stream = BufferedInputStream(IOBuffer("abcdefg"), 2)
        read(stream, 1)
        @test BufferedStreams.shiftdata!(stream) > 0
        read(stream, 2)
        @test BufferedStreams.shiftdata!(stream) > 0
    end

    @testset "misc." begin
        stream = BufferedInputStream(IOBuffer("foobar"), 10)
        @test !BufferedStreams.ensurebuffered!(stream, 10)
        repr_regex = r"^BufferedInputStream{.*}\(<.* \d+% filled>\)$"
        @test occursin(repr_regex, repr(stream))

        stream = BufferedInputStream(IOBuffer("foobar"), 4 * 2^10)
        @test occursin(repr_regex, repr(stream))
        @test occursin(r"KiB", repr(stream))

        close(stream)
        @test occursin(r"^BufferedInputStream{.*}\(<closed>\)$", repr(stream))
        @test_throws ArgumentError BufferedInputStream(IOBuffer("foo"), 0)
    end

    @testset "massive read" begin
        byte = 0x34
        bufsize = 1024
        stream = BufferedInputStream(InfiniteStream(byte), bufsize)

        # read byte by byte
        ok = true
        for i in 1:2^30
            ok = ok && (read(stream, UInt8) == byte)
        end
        @test ok
        @test length(stream.buffer) == 1024

        # read in one shot
        @test all(read(stream, 5 * bufsize) .== byte)
        @test length(stream.buffer) == 1024
    end

    @testset "readavailable" begin
        stream = BufferedInputStream(IOBuffer("some data"))
        @test readavailable(stream) == b"some data"
    end
end

@testset "BufferedOutputStream" begin
    @testset "write" begin
        data = rand(UInt8, 1000000)
        stream1 = BufferedOutputStream()
        sink = IOBuffer()
        stream2 = BufferedOutputStream(sink, 1024)
        for c in data
            write(stream1, c)
            write(stream2, c)
        end
        flush(stream1)
        flush(stream2)
        @test take!(stream1) == data
        @test take!(sink) == data
        close(stream1)
        close(stream2)
        @test !isopen(sink)
    end

    @testset "arrays" begin
        expected = UInt8[]
        stream1 = BufferedOutputStream()
        sink = IOBuffer()
        stream2 = BufferedOutputStream(sink, 1024)
        for _ in 1:1000
            data = rand(UInt8, rand(1:1000))
            append!(expected, data)
            write(stream1, data)
            write(stream2, data)
        end
        flush(stream1)
        flush(stream2)
        @test take!(stream1) == expected
        @test take!(sink) == expected
        close(stream1)
        close(stream2)
    end

    @testset "takebuf_string" begin
        data = rand('A':'z', 1000000)
        iobuf = IOBuffer()
        stream = BufferedOutputStream()
        for c in data
            write(stream, c)
            write(iobuf, c)
        end
        @test String(take!((stream))) == String(take!((iobuf)))
    end

    @testset "write_result" begin
        sink = IOBuffer()
        stream = BufferedOutputStream(sink, 16)
        for len in 0:10:100
            result = write(stream, repeat("x", len))
            @test result == len
        end
    end

    @testset "close" begin
        iobuffer = IOBuffer()
        stream = BufferedOutputStream(iobuffer)
        @test isopen(stream)
        @test isopen(iobuffer)
        write(stream, 0x00)
        close(stream)
        @test !isopen(stream)
        @test !isopen(iobuffer)
        @test_throws Exception write(stream, 0x00)
    end

    @testset "vector sink" begin
        sink = UInt8[]
        stream = BufferedOutputStream(sink)
        write(stream, 0x00)
        write(stream, 0x01)
        write(stream, 0x02)
        flush(stream)
        @test take!(stream) == [0x00, 0x01, 0x02]
        @test isopen(stream)
        close(stream)
        @test !isopen(stream)
        @test_throws Exception write(stream, 0x00)
    end

    @testset "iostream" begin
        mktemp() do path, out
            stream = BufferedOutputStream(out, 10)
            write(stream, "hello")
            @test stat(path).size == 0
            write(stream, "world")
            @test stat(path).size == 0
            write(stream, "!")
            # BufferedOutputStream buffer has run out of space,
            # but IOStream buffer has not
            @test stat(path).size == 0
            flush(stream)
            @test stat(path).size == 11
            write(stream, "...")
            @test stat(path).size == 11
            close(stream)
            @test !isopen(out)
            @test stat(path).size == 14
        end
    end

    @testset "misc." begin
        stream = BufferedOutputStream(IOBuffer(), 5)
        @test eof(stream)
        @test pointer(stream) == pointer(stream.buffer)
        @test occursin(r"^BufferedOutputStream{.*}\(<.* \d+% filled>\)$", string(stream))
        close(stream)
        @test occursin(r"^BufferedOutputStream{.*}\(<closed>\)$", string(stream))
        @test_throws ArgumentError BufferedOutputStream(IOBuffer(), 0)
    end
end