text:
@testset "isgeomagdatavalid.jl" begin
    nms = (:year, :month, :day, :hour, :minute, :second, :data)
    @test !isgeomagdatavalid(NamedTuple{nms}([2022, 12, 12, 3, 51, 46, 99999.99]))
    @test !isgeomagdatavalid(NamedTuple{nms}([2022, 6, 31, 3, 51, 46, 99.23]))
    @test isgeomagdatavalid(NamedTuple{nms}([2022, 1, 29, 23, 59, 32, 99.23]))
    @test !isgeomagdatavalid(NamedTuple{nms}([2022, 1, 29, 23, 59, 32, 100000]))
end


text:
function ftau!(champu :: Array{ComplexF64,2},
               champv :: Array{ComplexF64,2},
               m      :: MicMac,
               t      :: Float64,
               fft_u  :: Vector{ComplexF64},
               fft_v  :: Vector{ComplexF64})

    m.v  .= exp.(1im * t * m.A1)
    m.u .= -1im * m.llambda * m.A2 .* conj.(m.v)

    m.v .*= fft_u
    ifft!(m.v,1)
    m.ut .= m.v .* m.matr

    m.v  .= exp.(1im * t * m.A1)
    m.v .*= fft_v
    ifft!(m.v,1)
    m.vt .= m.v .* m.matr

    m.ut .+= conj.(m.vt)
    m.ut .*= 0.5
    m.vt = abs.(m.ut) .^ (2 * m.sigma) .* m.ut

    m.ut .= m.conjmatr .* m.vt
    fft!(m.ut,1)
    m.ut .*= m.u

    m.vt .= m.conjmatr .* conj.(m.vt)
    fft!(m.vt,1)
    m.vt .*= m.u

    transpose!(champu, m.ut)
    transpose!(champv, m.vt)

end

function ftau!(m     :: MicMac,
               t     :: Float64,
               fft_u :: Array{ComplexF64,2},
               fft_v :: Array{ComplexF64,2})

    m.v .= exp.(1im * t * m.A1)

    transpose!(m.ut, fft_u)
    transpose!(m.vt, fft_v)

    m.ut .*= m.v
    m.vt .*= m.v

    ifft!(m.ut,1)
    ifft!(m.vt,1)

    m.ut .*= m.matr
    m.vt .*= m.matr

    m.ut .+= conj.(m.vt)
    m.ut .*= 0.5
    m.vt .= abs.(m.ut) .^ (2 * m.sigma) .* m.ut

    m.u .= -1im * m.llambda .* m.A2 .* conj.(m.v)

    m.ut .= m.conjmatr .* m.vt
    m.vt .= m.conjmatr .* conj.(m.vt)

    fft!(m.ut,1)
    fft!(m.vt,1)

    m.ut .*= m.u
    m.vt .*= m.u

    transpose!(fft_u, m.ut)
    transpose!(fft_v, m.vt)

end

ext: jl


text:
######### StanSample Bernoulli example  ###########

using StanSample

ProjDir = @__DIR__

bernoulli_model = "
data {
  int<lower=1> N;
  int<lower=0,upper=1> y[N];
}
parameters {
  real<lower=0,upper=1> theta;
}
model {
  theta ~ beta(1,1);
  y ~ bernoulli(theta);
}
";

data = Dict("N" => 10, "y" => [0, 1, 0, 1, 0, 0, 0, 0, 0, 1])

# Keep tmpdir across multiple runs to prevent re-compilation
tmpdir = joinpath(@__DIR__, "tmp")

sm = SampleModel("kw_bern", bernoulli_model, tmpdir);

sm |> display

rc = stan_sample(sm; data, num_threads=4, num_cpp_chains=4, num_chains=2, seed=12);

if success(rc)
  st = read_samples(sm)
  display(st)
  println()
  display(read_samples(sm, :dataframe))
end

max_stars_repo_name: StanJulia/StanSample

max_issues_repo_name: StanJulia/StanSample

max_forks_repo_name: StanJulia/StanSample


text:
function value(obj::RealFieldUnivariateMatrixFunction, arg0::RealFieldElement)
    return jcall(obj, "value", Vector{Vector{RealFieldElement}}, (RealFieldElement,), arg0)
end


text:
function get_eop_history(obj::TODProvider)
    return jcall(obj, "getEOPHistory", EOPHistory, ())
end

function get_non_interpolating_provider(obj::TODProvider)
    return jcall(obj, "getNonInterpolatingProvider", TODProvider, ())
end

function get_transform(obj::TODProvider, arg0::AbsoluteDate)
    return jcall(obj, "getTransform", Transform, (AbsoluteDate,), arg0)
end

function get_transform(obj::TODProvider, arg0::FieldAbsoluteDate)
    return jcall(obj, "getTransform", FieldTransform, (FieldAbsoluteDate,), arg0)
end


text:
using BenchmarkTools
using Trixi

const SUITE = BenchmarkGroup()

let dimension = "2d"
  SUITE[dimension] = BenchmarkGroup()
  EXAMPLES_DIR = joinpath(examples_dir(), dimension)
  for elixir in ["elixir_advection_extended.jl", "elixir_advection_amr_nonperiodic.jl",
                  "elixir_euler_ec.jl", "elixir_euler_vortex_mortar.jl", "elixir_euler_vortex_mortar_shockcapturing.jl"]
    SUITE[dimension][elixir] = BenchmarkGroup()
    for polydeg in [3, 7]
      trixi_include(joinpath(EXAMPLES_DIR, elixir), tspan=(0.0, 1.0e-10); polydeg)
      SUITE[dimension][elixir]["p$(polydeg)_rhs!"] = @benchmarkable Trixi.rhs!(
        $(similar(sol.u[end])), $(copy(sol.u[end])), $(semi), $(first(tspan)))
      SUITE[dimension][elixir]["p$(polydeg)_analysis"] = @benchmarkable ($analysis_callback)($sol)
    end
  end
end

let dimension = "3d"
  SUITE[dimension] = BenchmarkGroup()
  EXAMPLES_DIR = joinpath(examples_dir(), dimension)
  for elixir in ["elixir_advection_extended.jl",
                 "elixir_euler_ec.jl", "elixir_euler_mortar.jl"]
    SUITE[dimension][elixir] = BenchmarkGroup()
    for polydeg in [3, 7]
      trixi_include(joinpath(EXAMPLES_DIR, elixir), tspan=(0.0, 1.0e-10); polydeg)
      SUITE[dimension][elixir]["p$(polydeg)_rhs!"] = @benchmarkable Trixi.rhs!(
        $(similar(sol.u[end])), $(copy(sol.u[end])), $(semi), $(first(tspan)))
      SUITE[dimension][elixir]["p$(polydeg)_analysis"] = @benchmarkable ($analysis_callback)($sol)
    end
  end
end

let
  SUITE["latency"] = BenchmarkGroup()
  SUITE["latency"]["default_example"] = @benchmarkable run(
    `$(Base.julia_cmd()) -e 'using Trixi; trixi_include(default_example())'`) seconds=60
  for polydeg in [3, 7]
    command = "using Trixi; trixi_include(joinpath(examples_dir(), \"2d\", \"elixir_advection_extended.jl\"), tspan=(0.0, 1.0e-10), polydeg=$(polydeg), save_restart=TrivialCallback(), save_solution=TrivialCallback())"
    SUITE["latency"]["polydeg_$polydeg"] = @benchmarkable run($`$(Base.julia_cmd()) -e $command`) seconds=60
  end
end


text:
addprocs(2)

using Base.Test
using Dagger
using DagScheduler

@everywhere include("customsch.jl")

include("domain.jl")
include("array.jl")
Dagger.cleanup()
#include("cache.jl")


text:
# TODO: doc is missing

function variance(a)
    avg = a / length(a)
    x = sum(a - avg)
    s_sq = x / (length(a) - 1)
    return(s_sq)
end

max_stars_repo_name: ashwani-rathee/Julia

max_issues_repo_name: ashwani-rathee/Julia

max_forks_repo_name: ashwani-rathee/Julia


text:
function read_groups(lines)
    groups = []
    current_group = []
    for line in lines
        if isempty(line)
            push!(groups, current_group)
            current_group = []
        else
            push!(current_group, line)
        end
    end
    push!(groups, current_group)
    groups
end

function calc_group_sum(groups)
    s = 0
    for group in groups
        people = map(item -> Set(item), group)
        common = reduce(intersect, people)
        s += length(common)
    end
    s
end

function main()
    lines = readlines("../resources/input_06.txt")
    groups = read_groups(lines)
    result = calc_group_sum(groups)
    println(result)
end

main()


text:
using Gnuplot, Random
Random.seed!(1234)
matrix = randn(50,60)
@gp matrix "w image notit" "set auto fix"
save(term="pngcairo size 800,800", output="./code/plt3d_ex5.png") # hide

max_stars_repo_name: lazarusA/gnuplot-examples

max_issues_repo_name: lazarusA/gnuplot-examples

max_forks_repo_name: lazarusA/gnuplot-examples


text:
"""
Blackbody radiation.

* [`ephoton`](@ref)
* [`energy2photon`](@ref)
* [`planck`](@ref)
* [`stefan_boltzmann`](@ref)
* [`blackbody_temp`](@ref)
"""
module Blackbody

include("../docstring_style.jl")

using Canopy.Constants: c, h, N_A, c_1L, c_2, sigma_SB

export ephoton,
    energy2photon,
    planck,
    stefan_boltzmann,
    blackbody_temp

"""
Calculate the energy carried by a photon [J] from its wavelength [m].

# Examples

```jldoctest
julia> ephoton(0.5e-6)
3.972891714297857e-19
```
"""
function ephoton(wl)
    return h * c / wl
end

"""
Convert the energy flux density `flux` [W m^-2] at the wavelength `wl` [m] to
photon flux density [mol m^-2 s^-1].

# Examples

```jldoctest
julia> energy2photon(300., 0.5e-6)
0.0012539020843667665
```
"""
function energy2photon(flux, wl)
    return flux * wl / (N_A * h * c)
end

"""
Calculate the spectral radiance [W sr^-1 m^-3] from wavelength `wl` [m] and
blackbody temperature `temp` [K] accroding to Planck's law.

# Examples

```jldoctest
julia> planck(0.5e-6, 6000)
3.175690669942503e13

julia> planck(10e-6, 288)
8.114245536536237e6
```
"""
function planck(wl, temp)
    return c_1L * wl^(-5.) / (exp(c_2 / (wl * temp)) - 1.)
end

"""
Calculate the energy flux density [W m^-2] of a blackbody of temperature `temp`
[K] from the Stefan--Boltzmann law.

# Examples

```jldoctest
julia> stefan_boltzmann(298.15)
448.0752867066491
```
"""
function stefan_boltzmann(temp)
    temp_sq = temp * temp
    return sigma_SB * temp_sq * temp_sq
end

"""
Calculate the blackbody temperature [K] from energy flux density [W m^-2].

# Examples

```jldoctest
julia> blackbody_temp(240.)
255.06441705528474
```
"""
function blackbody_temp(flux)
    return (flux / sigma_SB)^0.25
end

end  # module


text:
import Base: copy, +

#*********************************#
#       COST FUNCTION CLASS       #
#*********************************#

"""
Abstract type that represents a scalar-valued function that accepts a state and control
at a single knot point.
"""
abstract type CostFunction end

"""
An abstract type that represents any [`CostFunction`](@ref) of the form
```math
\\frac{1}{2} x^T Q x + \\frac{1}{2} u^T R u + u^T H x + q^T x + r^T u + c
```

These types all support the following methods
* [`is_diag`](@ref)
* [`is_blockdiag`](@ref)
* [`invert!`](@ref)

As well as standard addition.
"""
abstract type QuadraticCostFunction{n,m,T} <: CostFunction end

"""
    is_diag(::QuadraticCostFunction)

Determines if the hessian of a quadratic cost function is strictly diagonal.
"""
is_diag(cost::QuadraticCostFunction) = is_blockdiag(cost) && cost.Q isa Diagonal && cost.R isa Diagonal

"""
    is_diag(::QuadraticCostFunction)

Determines if the hessian of a quadratic cost function is block diagonal (i.e. ``\\norm(H) = 0``).
"""
is_blockdiag(::QuadraticCostFunction) = false
state_dim(::QuadraticCostFunction{n}) where n = n
control_dim(::QuadraticCostFunction{<:Any,m}) where m = m

function (::Type{QC})(Q::AbstractArray, R::AbstractArray;
        H::AbstractArray=SizedMatrix{size(R,1),size(Q,1)}(zeros(eltype(Q), size(R,1), size(Q,1))),
        q::AbstractVector=(@SVector zeros(eltype(Q), size(Q,1))),
        r::AbstractVector=(@SVector zeros(eltype(R), size(R,1))),
        c::Real=zero(eltype(Q)), kwargs...) where QC <: QuadraticCostFunction
    QC(Q, R, H, q, r, c; kwargs...)
end

function QuadraticCostFunction(Q::AbstractArray, R::AbstractArray,
        H::AbstractArray, q::AbstractVector, r::AbstractVector, c::Real;
        kwargs...)
    if LinearAlgebra.isdiag(Q) && LinearAlgebra.isdiag(R) && norm(H) ≈ 0
        DiagonalCost(diag(Q), diag(R), q, r, c)
    else
        QuadraticCost(Q, R, H, q, r, c; kwargs...)
    end
end

"""
    stage_cost(costfun::CostFunction, x, u)
    stage_cost(costfun::CostFunction, x)

Calculate the scalar cost using `costfun` given state `x` and control `u`. If only the
state is provided, it is assumed it is a terminal cost.
"""
function stage_cost(cost::QuadraticCostFunction, x::AbstractVector, u::AbstractVector)
    J = 0.5*u'cost.R*u + dot(cost.r,u) + stage_cost(cost, x)
    if !is_blockdiag(cost)
        J += u'cost.H*x
    end
    return J
end

function stage_cost(cost::QuadraticCostFunction, x::AbstractVector{T}) where T
    0.5*x'cost.Q*x + dot(cost.q,x) + cost.c
end

"""
    gradient!(E::QuadraticCostFunction, costfun::CostFunction, x, u)
    gradient!(E::QuadraticCostFunction, costfun::CostFunction, x)

Evaluate the gradient of the cost function `costfun` at state `x` and control `u`, storing
    the result in `E.q` and `E.r`. Return a `true` if the gradient is constant, and `false`
    otherwise.
"""
function gradient!(E::QuadraticCostFunction, cost::QuadraticCostFunction, x)
    E.q .= cost.Q*x .+ cost.q
    return false
end

function gradient!(E::QuadraticCostFunction, cost::QuadraticCostFunction, x, u)
    gradient!(E, cost, x)
    E.r .= cost.R*u .+ cost.r
    if !is_blockdiag(cost)
        E.q .+= cost.H'u
        E.r .+= cost.H*x
    end
    return false
end

"""
    hessian!(E::QuadraticCostFunction, costfun::CostFunction, x, u)
    hessian!(E::QuadraticCostFunction, costfun::CostFunction, x)

Evaluate the hessian of the cost function `costfun` at state `x` and control `u`, storing
    the result in `E.Q`, `E.R`, and `E.H`. Return a `true` if the hessian is constant, and `false`
    otherwise.
"""
function hessian!(E::QuadraticCostFunction, cost::QuadraticCostFunction, x)
    if is_diag(cost)
        for i = 1:length(x); E.Q[i,i] = cost.Q[i,i] end
    else
        E.Q .= cost.Q
    end
    return true
end

function hessian!(E::QuadraticCostFunction, cost::QuadraticCostFunction, x, u)
    hessian!(E, cost, x)
    if is_diag(cost)
        for i = 1:length(u); E.R[i,i] = cost.R[i,i]; end
    else
        E.R .= cost.R
    end
    if !is_blockdiag(cost)
        E.H .= cost.H
    end
    return true
end

function Base.copy(c::QC) where QC<:QuadraticCostFunction
    QC(copy(c.Q), copy(c.R), H=copy(c.H), q=copy(c.q), r=copy(c.r), c=c.c,
        terminal=c.terminal, checks=false)
end

# Other methods
"""
    set_LQR_goal!(cost::QuadraticCostFunction, xf)
    set_LQR_goal!(cost::QuadraticCostFunction, xf, uf)

Change the reference state and control for an LQR tracking cost.
Only changes `q` and `r`, and not the constant term `c`. If `uf` is
not passed in, it isn't changed.
"""
function set_LQR_goal!(cost::QuadraticCostFunction, xf)
    cost.q .= -cost.Q * xf
    return nothing
end

function set_LQR_goal!(cost::QuadraticCostFunction, xf, uf)
    set_LQR_goal!(cost, xf)
    cost.r .= -cost.R * uf
    return nothing
end

function +(c1::QuadraticCostFunction, c2::QuadraticCostFunction)
    @assert state_dim(c1) == state_dim(c2)
    @assert control_dim(c1) == control_dim(c2)
    n,m = state_dim(c1), control_dim(c1)
    H1 = is_diag(c1) ? zeros(m,n) : c1.H
    H2 = is_diag(c2) ? zeros(m,n) : c2.H
    QC = promote_type(typeof(c1), typeof(c2))
    QC(c1.Q + c2.Q, c1.R + c2.R, H1 + H2,
       c1.q + c2.q, c1.r + c2.r, c1.c + c2.c,
       checks=false, terminal=c1.terminal && c2.terminal)
end

"""
    invert!(Ginv, cost::QuadraticCostFunction)

Invert the hessian of the cost function, storing the result in `Ginv`. Performs the inversion
    efficiently, depending on the structure of the Hessian (diagonal or block diagonal).
"""
function invert!(Ginv, cost::QuadraticCostFunction{n,m}) where {n,m}
    ix = 1:n
    iu = n .+ (1:m)
    if is_diag(cost)
        for i = 1:n; Ginv[i,i] = inv(cost.Q[i,i]); end
        if !cost.terminal
            for i = 1:m; Ginv[i+n,i+n] = inv(cost.R[i,i]); end
        end
    elseif is_blockdiag(cost)
        Ginv[ix,ix] .= inv(SizedMatrix{n,n}(cost.Q))
        if !cost.terminal
            Ginv[iu,iu] .= inv(SizedMatrix{m,m}(cost.R))
        end
    else
        G1 = [cost.Q cost.H']
        G2 = [cost.H cost.R ]
        G = [G1; G2]
        Ginv .= inv(G)
    end
end

#######################################################
#              COST FUNCTION INTERFACE                #
#######################################################

"""
    DiagonalCost{n,m,T}

Cost function of the form
```math
\\frac{1}{2} x^T Q x + \\frac{1}{2} u^T R u + q^T x + r^T u + c
```
where ``Q`` and ``R`` are positive semi-definite and positive definite diagonal matrices,
respectively, and ``x`` is `n`-dimensional and ``u`` is `m`-dimensional.

# Constructors
    DiagonalCost(Qd, Rd, q, r, c; kwargs...)
    DiagonalCost(Q, R, q, r, c; kwargs...)
    DiagonalCost(Qd, Rd; [q, r, c, kwargs...])
    DiagonalCost(Q, R; [q, r, c, kwargs...])

where `Qd` and `Rd` are the diagonal vectors, and `Q` and `R` are matrices.

Any optional or omitted values will be set to zero(s). The keyword arguments are
* `terminal` - A `Bool` specifying if the cost function is terminal cost or not.
* `checks` - A `Bool` specifying if `Q` and `R` will be checked for the required definiteness.
"""
struct DiagonalCost{n,m,T} <: QuadraticCostFunction{n,m,T}
    Q::Diagonal{T,SVector{n,T}}
    R::Diagonal{T,SVector{m,T}}
    q::MVector{n,T}
    r::MVector{m,T}
    c::T
    terminal::Bool
    function DiagonalCost(Qd::StaticVector{n}, Rd::StaticVector{m},
                          q::StaticVector{n},  r::StaticVector{m},
                          c::Real; terminal::Bool=false, checks::Bool=true, kwargs...) where {n,m}
        T = promote_type(typeof(c), eltype(Qd), eltype(Rd), eltype(q), eltype(r))
        if checks
            if any(x->x<0, Qd)
                @warn "Q needs to be positive semi-definite."
            elseif any(x->x<=0, Rd) && !terminal
                @warn "R needs to be positive definite."
            end
        end
        new{n,m,T}(Diagonal(SVector(Qd)), Diagonal(SVector(Rd)), q, r, T(c), terminal)
    end
end

Base.copy(c::DiagonalCost) =
    DiagonalCost(c.Q.diag, c.R.diag, copy(c.q), copy(c.r), c.c,
        checks=false, terminal=c.terminal)

# strip away the H
@inline function (::Type{<:DiagonalCost})(Q::AbstractArray, R::AbstractArray, H::AbstractArray,
        q::AbstractVector, r::AbstractVector, c; kwargs...)
    DiagonalCost(Q, R, q, r, c; kwargs...)
end

# convert to form for inner constructor
function DiagonalCost(Q::AbstractArray, R::AbstractArray,
        q::AbstractVector, r::AbstractVector, c; kwargs...)
    n,m = length(q), length(r)
    Qd = SVector{n}(diag(Q))
    Rd = SVector{m}(diag(R))
    DiagonalCost(Qd, Rd, convert(MVector{n},q), convert(MVector{m},r), c; kwargs...)
end

# Pass in vectors for Q and R
function DiagonalCost(Q::AbstractVector, R::AbstractVector,
        q::AbstractVector, r::AbstractVector, c::Real; kwargs...)
    DiagonalCost(Diagonal(Q), Diagonal(R), q, r, c; kwargs...)
end

function DiagonalCost(Q::AbstractVector, R::AbstractVector; kwargs...)
    DiagonalCost(Diagonal(Q), Diagonal(R); kwargs...)
end

is_blockdiag(::DiagonalCost) = true
is_diag(::DiagonalCost) = true

function LinearAlgebra.inv(cost::DiagonalCost)
    return DiagonalCost(inv(cost.Q), inv(cost.R), cost.q, cost.r, cost.c, terminal=cost.terminal)
end

function Base.:\(cost::DiagonalCost, z::AbstractKnotPoint)
    x = cost.Q\state(z)
    u = cost.R\control(z)
    return StaticKnotPoint([x;u], z._x, z._u, z.dt, z.t)
end

function change_dimension(cost::DiagonalCost, n::Int, m::Int, ix, iu)
    Qd = zeros(n)
    Rd = zeros(m)
    q = zeros(n)
    r = zeros(m)
    Qd[ix] = diag(cost.Q)
    Rd[iu] = diag(cost.R)
    q[ix] = cost.q
    r[iu] = cost.r
    DiagonalCost(Qd, Rd, q, r, cost.c, terminal=cost.terminal, checks=false)
end

"""
    QuadraticCost{n,m,T,TQ,TR}

Cost function of the form
```math
\\frac{1}{2} x^T Q x + \\frac{1}{2} u^T R u + u^T H x + q^T x + r^T u + c
```
where ``R`` must be positive definite, ``Q`` and ``Q_f`` must be positive semidefinite.

The type parameters `TQ` and `TR` specify the type of ``Q`` and ``R``.

# Constructor
    QuadraticCost(Q, R, H, q, r, c; kwargs...)
    QuadraticCost(Q, R; H, q, r, c, kwargs...)

Any optional or omitted values will be set to zero(s). The keyword arguments are
* `terminal` - A `Bool` specifying if the cost function is terminal cost or not.
* `checks` - A `Bool` specifying if `Q` and `R` will be checked for the required definiteness.
"""
mutable struct QuadraticCost{n,m,T,TQ,TR} <: QuadraticCostFunction{n,m,T}
    Q::TQ                     # Quadratic stage cost for states (n,n)
    R::TR                     # Quadratic stage cost for controls (m,m)
    H::SizedMatrix{m,n,T,2,Matrix{T}}   # Quadratic Cross-coupling for state and controls (m,n)
    q::MVector{n,T}           # Linear term on states (n,)
    r::MVector{m,T}           # Linear term on controls (m,)
    c::T                      # constant term
    terminal::Bool
    zeroH::Bool
    function (::Type{QC})(Q::TQ, R::TR, H::TH,
            q::Tq, r::Tr, c::Real; checks=true, terminal=false, kwargs...) where {TQ,TR,TH,Tq,Tr,QC<:QuadraticCost}
        @assert size(Q,1) == length(q)
        @assert size(R,1) == length(r)
        @assert size(H) == (length(r), length(q))
        if checks
            if !isposdef(Array(R)) && !terminal
                @warn "R is not positive definite"
            end
            if !ispossemidef(Array(Q))
                @warn "Q is not positive semidefinite"
            end
        end
        zeroH = norm(H,Inf) ≈ 0
        m,n = size(H)
        T = promote_type(eltype(Q), eltype(R), eltype(H), eltype(q), eltype(r), typeof(c))
        new{n,m,T,TQ,TR}(Q,R,H,q,r,c,terminal,zeroH)
    end
    function QuadraticCost{n,m,T,TQ,TR}(qcost::QuadraticCost) where {n,m,T,TQ,TR}
        new{n,m,T,TQ,TR}(qcost.Q, qcost.R, qcost.H, qcost.q, qcost.r, qcost.c, qcost.terminal, qcost.zeroH)
    end
end

state_dim(cost::QuadraticCost) = length(cost.q)
control_dim(cost::QuadraticCost) = length(cost.r)
is_blockdiag(cost::QuadraticCost) = cost.zeroH

function QuadraticCost{T}(n::Int,m::Int; terminal=false) where T
    Q = SizedMatrix{n,n}(Matrix(one(T)*I,n,n))
    R = SizedMatrix{m,m}(Matrix(one(T)*I,m,m))
    H = SizedMatrix{m,n}(zeros(T,m,n))
    q = SizedVector{n}(zeros(T,n))
    r = SizedVector{m}(zeros(T,m))
    c = zero(T)
    QuadraticCost(Q,R,H,q,r,c, checks=false, terminal=terminal)
end

QuadraticCost(cost::QuadraticCost) = cost

function LinearAlgebra.inv(cost::QuadraticCost)
    if norm(cost.H,Inf) ≈ 0
        QuadraticCost(inv(cost.Q), inv(cost.R), cost.H, cost.q, cost.r, cost.c,
            checks=false, terminal=cost.terminal)
    else
        m,n = size(cost.H)
        H1 = [cost.Q cost.H']
        H2 = [cost.H cost.R ]
        H = [H1; H2]
        Hinv = inv(H)
        ix = 1:n
        iu = n .+ (1:m)
        Q = SizedMatrix{n,n}(Hinv[ix, ix])
        R = SizedMatrix{m,m}(Hinv[iu, iu])
        H = SizedMatrix{m,n}(Hinv[iu, ix])
        QuadraticCost(Q,R,H, cost.q, cost.r, cost.c,
            checks=false, terminal=cost.terminal)
    end
end

@inline function (::Type{<:QuadraticCost})(dcost::DiagonalCost)
    QuadraticCost(dcost,Q, dcost.R, q=dcost.q, r=dcost.r, c=dcost.c, terminal=dcost.terminal)
end

function Base.convert(::Type{<:QuadraticCost}, dcost::DiagonalCost)
    QuadraticCost(dcost.Q, dcost.R, q=dcost.q, r=dcost.r, c=dcost.c, terminal=dcost.terminal)
end

Base.promote_rule(::Type{<:QuadraticCostFunction}, ::Type{<:QuadraticCostFunction}) = QuadraticCost

function Base.promote_rule(::Type{<:QuadraticCost{n,m,T1,Q1,R1}},
                           ::Type{<:QuadraticCost{n,m,T2,Q2,R2}}) where {n,m,T1,T2,Q1,Q2,R1,R2}
    T = promote_type(T1,T2)
    function diag_type(T1,T2,n)
        elT = promote_type(eltype(T1), eltype(T2))
        if T1 == T2
            return T1
        elseif (T1 <: Diagonal) && (T2 <: Diagonal)
            return Diagonal{elT,MVector{n,elT}}
        else
            return SizedMatrix{n,n,elT,2}
        end
    end
    Q = diag_type(Q1,Q2,n)
    R = diag_type(R1,R2,m)
    QuadraticCost{n,m,T, Q, R}
end

@inline Base.convert(::Type{QC}, cost::QuadraticCost) where QC <: QuadraticCost = QC(cost)

"""
    LQRCost(Q, R, xf, [uf; kwargs...])

Convenience constructor for a `QuadraticCostFunction` of the form:
```math
\\frac{1}{2} (x-x_f)^T Q (x-xf) + \\frac{1}{2} (u-u_f)^T R (u-u_f)
```

If ``Q`` and ``R`` are diagonal, the output will be a `DiagonalCost`, otherwise it will
be a `QuadraticCost`.
"""
function LQRCost(Q::AbstractArray, R::AbstractArray,
        xf::AbstractVector, uf=(@SVector zeros(size(R,1))); kwargs...)
    H = @SMatrix zeros(size(R,1),size(Q,1))
    q = -Q*xf
    r = -R*uf
    c = 0.5*xf'*Q*xf + 0.5*uf'R*uf
    return QuadraticCostFunction(Q, R, H, q, r, c, kwargs...)
end

function LQRCost(Q::Diagonal{<:Any,<:SVector{n}},R::Diagonal{<:Any,<:SVector{m}},
        xf::AbstractVector, uf=(@SVector zeros(m)); kwargs...) where {n,m}
    q = -Q*xf
    r = -R*uf
    c = 0.5*xf'*Q*xf + 0.5*uf'R*uf
    return DiagonalCost(Q, R, q, r, c; kwargs...)
end


text:
export ## Types,
       Dipole,
       FlipFlop,
       FluxTransmon,
       SemiClassicalDipole,
       ## Methods
       update_params

#Iteractions linearly add new Hamiltonians to the system

type Dipole <: Interaction
    system1::QSystem
    system2::QSystem
    strength::Float64
end

function hamiltonian(f::Dipole)
    return f.strength*kron(X(f.system1), X(f.system2))
end

type FlipFlop <: Interaction
    system1::QSystem
    system2::QSystem
    strength::Float64
end
function hamiltonian(f::FlipFlop)
    return f.strength*(kron(raising(f.system1), lowering(f.system2)) + kron(lowering(f.system1), raising(f.system2)))
end

type SemiClassicalDipole <: Interaction
    system1::Field
    system2::QSystem
    strength::Float64
end
function hamiltonian(scd::SemiClassicalDipole, t::Float64)
    return scd.strength*amplitude(scd.system1, t)*X(scd.system2)
end

# ParametricInteractions are interactions which modify
# parameters of systems rather than linearly add new Hamiltonians

#Flux control
type FluxTransmon <: ParametricInteraction
    flux::Field
    transmon::String #refer to transmon by label in parent CompositeQSystem to avoid broken reference on copy
    strength::Float64
end
FluxTransmon(flux::Field, transmon::TunableTransmon, strength::Float64) = FluxTransmon(flux, label(transmon), strength)

#Set the flux in the associated transmon
function update_params(sys::CompositeQSystem, ft::FluxTransmon, t::Float64)
    myTransmon = sys[ft.transmon]
    myTransmon.flux = myTransmon.fluxBias + ft.strength*amplitude(ft.flux,t)
end


text:
"""
See documentation at https://github.com/baggepinnen/Robotlib.jl
"""
module Robotlib
using LinearAlgebra, Statistics, StaticArrays, SparseArrays
using TotalLeastSquares
using Quaternions
import Quaternions: Quaternion, rotationmatrix
using Optim
export rotationmatrix
const I4 = SMatrix{4,4,Float64,16}(Matrix{Float64}(I, 4, 4))
const I3 = SMatrix{3,3,Float64,9}(Matrix{Float64}(I, 3, 3))
include("DH.jl")
include("utils.jl")
include("kinematics.jl")
include("Frames.jl") # Must come before robotplot
include("robotplot.jl")
include("read_log.jl")
include("posDepFric.jl")

include("calibLPOE.jl")
include("calib_force.jl")
include("calibNAXP.jl")
include("calibAXYB.jl")

export calibLPOE,
    calibLPOEdual, calib_force, calib_force_iterative, calib_force_eigen, calibNAXP, calibAXYB, leds_to_tf

using Requires
function __init__()
    @require Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80" begin
        function plotLines(points, lines)
            P = size(points, 2)
            for i = 1:P
                p = points[1:3, i]
                l = lines[1:3, i]
                Plots.plot3d!(
                    eachcol([p' + 0.05 * l'; p' - 0.05 * l'])...,
                    c = :red,
                    lab = "",
                )
            end
            display(Plots.current())
        end

        function plotPlanes(normals)
            P = size(normals, 1)
            p = zeros(P, 4)
            Plots.plot()
            for i = 1:P
                n = normals[i, :]
                xdir = n + [1, 0, 0]
                ydir = normalize(cross(n, xdir))
                xdir = normalize(cross(ydir, n))
                p[:, 1] = n + xdir + ydir
                p[:, 2] = n + xdir - ydir
                p[:, 3] = n - xdir - ydir
                p[:, 4] = n - xdir + ydir

                Plots.surface!(p[1, :], p[2, :], p[3, :])
            end
            Plots.plot!(alpha = 0.5, zlims = (-0.3, 1))
        end
    end
end

export skewcoords,
    twistcoords,
    skewify,
    skew,
    skew4,
    expω,
    expξ,
    expξ2,
    expξ!,
    logT,
    logR,
    trinv,
    isrot,
    isse3,
    Rangle,
    conformize,
    DH2twistsPOE,
    DH2twistsLPOE,
    dh2Tn,
    orthonormal,
    orthonormal!,
    rpy2R,
    Quaternion,
    xyθ,
    smartdiff,
    R2rpy,
    traj2quat,
    centraldiff
export fkinePOE,
    fkineLPOE,
    ikinePOE,
    dh2Tn,
    jacobianPOE,
    jacobianPOEikine,
    jacobian,
    get_kinematic_functions
export trajplot, trajplot3, plot3smart
export DH,
    DH2400,
    DHYuMi,
    DH7600,
    DHtest,
    abb2logical!,
    logical2abb!,
    abb2logical,
    logical2abb
export csv2dict, csv2mat, orcalog2mat, getdata
export frictionRBFN, getcenters

dir(x...) = joinpath(dirname(pathof(Robotlib)), x...)

end

max_forks_repo_name: baggepinnen/Robotlib


text:
"""
The Mathenize module exists mainly to make math calculations in
strings.
"""
module Mathenize
export calculate
include("computation.jl")

"""
    calculate(math::String, print_info::Bool=false)

Return calculated value(s) in a string.

# Arguments
- `math::String` Math operation.
- `print_info::Bool` Shows log of what happened during the calculation.

# Examples
```julia-repl
julia> calculate("sqrt(complex(-90)) + 10im")
0.0 + 19.486832980505138im
```

"""
function calculate(math::String, print_info::Bool=false)

    r = nothing
    math = Meta.parse(math)
    hasproperty(math, :args) ? tasks = length(math.args) : tasks = 0

    #Adding to log info.
    LOG_INFO = []
    push!(LOG_INFO, "-> $(math) <- \n └Tasks: $(tasks)\n └$(hasproperty(math, :args) ? math.args : "Empty")")
    if math !== nothing
    if tasks >= 1
        push!(LOG_INFO, "   └ Performing subtasks:")
        r = subtasking(math, tasks, math.args, LOG_INFO, print_info)
    elseif math isa Number
        push!(LOG_INFO, "   └$(math) is a number")
        r = math
    elseif math in sym
        push!(LOG_INFO, "   └$(math) is a valid symbol")
        r = Core.eval(Base.Math, math)
    elseif tasks == 1 && ispermitted(math, LOG_INFO)
        push!(LOG_INFO, "   └$(math) is a valid a permitted number or expression")
        r = Core.eval(Base.Math, math)
    else
        r = nothing
        unknownmath(math, LOG_INFO, print_info)
    end
    end

    if print_info == true
        @info join(LOG_INFO, "\n")
    end

    return r #Return result

end
end #End module


text:
using BlackBoxOptim: makefitness, NewFitness, SingleObjectiveFitness, SingleNumberFitness
using BlackBoxOptim: MultiObjectiveFitness, VectorFitness
using BlackBoxOptim: fitnessvalues, characteristics, simulations, focusversion, goalversion, evaluator

describe("Fitness types") do
  @repeat test("SingleObjectiveFitness") do
    fval = randn()
    f1 = makefitness(fval)
    @check typeof(f1) <: NewFitness
    @check typeof(f1) <: SingleObjectiveFitness
    @check typeof(f1) == SingleNumberFitness{Float64}

    @check fitnessvalues(f1) == [fval]
    @check characteristics(f1) == [fval]
    @check simulations(f1) == [fval]
    @check evaluator(f1) == nothing
    @check focusversion(f1) == 0
    @check goalversion(f1) == 0
  end

  @repeat test("MultiObjectiveFitness") do
    size = rand(2:117)
    vals = randn(size)
    f2 = makefitness(vals)
    @check typeof(f2) <: NewFitness
    @check typeof(f2) <: MultiObjectiveFitness
    @check typeof(f2) == VectorFitness{Float64}

    @check fitnessvalues(f2) == vals
    @check characteristics(f2) == vals
    @check simulations(f2) == vals
    @check evaluator(f2) == nothing
    @check focusversion(f2) == 0
    @check goalversion(f2) == 0
  end
end


text:
@testmodule DataReaderTest begin
  using MLLibrary

  function test_DataSetMemory()

    resources=["abalone.data","3percent-miscategorization.csv","10percent-miscatergorization.csv","Seperable.csv"]
    for res in resources
      delimiter=","
      path = joinpath(Pkg.dir("MLLibrary"),"resources","datasets",res)
      count = MLLibrary.DataParser.DataStream.getDataSetCount(path)

      featureCount = MLLibrary.DataParser.DataStream.getFeatureCount(path)
      threshold=.1
      srand(123)
      trainingsetindices = MLLibrary.DataSampling.ThresholdSampling.getThresholdTrainingSetIndicies(path,count,threshold)
      srand(123)
      testsetindices = MLLibrary.DataSampling.ThresholdSampling.getThresholdTestSetIndices(path,count,threshold)

      trainingset = MLLibrary.DataSampling.DataReader.getDataSetInMemoryFromIndices(path,delimiter,trainingsetindices)

      testset = MLLibrary.DataSampling.DataReader.getDataSetInMemoryFromIndices(path,delimiter,testsetindices)
      @test size(trainingset)==(convert(Int,round(count*(1-threshold))),featureCount)
      @test size(testset)==(convert(Int,round(count*threshold)),featureCount)

    end

  end

end

max_stars_repo_name: ajermaky/CustomMLLibrary

max_issues_repo_name: ajermaky/CustomMLLibrary

max_forks_repo_name: ajermaky/CustomMLLibrary


text:
#!/usr/bin/env julia
# Pattern Matching and Regular Expressions
# https://teaching.healthtech.dtu.dk/22110/index.php/Pattern_Matching_and_Regular_Expressions

# Exercise 1
println("Check if is number:")
num = chomp(readline())
if match(r"^-?\d+(\.\d+)?$", num) !== nothing
    println("Yes, it is a number")
else
    println("No, it is not a number")
end

# Exercise 2
struct InvalidLengthError <: Exception
    val1::Int64
    val2::Int64
end
Base.showerror(io::IO, e::InvalidLengthError) = print(io, "The parsed length (", e.val1, ") does not equal the calculated length (",e.val2,")")

function parse_prot(f)
    seq = ""
    sp_id = ""
    sp_ac = ""
    aa_len = 0
    FLAG = false

    for line in eachline(f)
        key, value = line[begin:2], rstrip(line[6:end])
        if key == "ID"
            id_reg = match(r"^(\w+)\s+", value)
            sp_id = string(id_reg[1])
            println("The ID is: ", sp_id)
        elseif key == "AC"
            ac_reg = match(r"^(\w+);", value)
            sp_ac = string(ac_reg[1])
            println("The accession number is: ", sp_ac)
        end

        if key == "//"
            FLAG = false
        end
        if FLAG
            seq_parts_reg = findall(r"(\w+)", value)
            seq_parts = [value[x] for x in seq_parts_reg]
            seq *= join(seq_parts)
        end
        if key == "SQ"
            FLAG = true
            aa_end = aa_start = 12
            while value[aa_end] != ' '
                aa_end += 1
            end
            aa_len = parse(Int64, value[aa_start:aa_end])
        end
    end
    println("Sequence: ",seq)

    if aa_len != length(seq)
        throw(InvalidLengthError(aa_len,length(seq)))
    else
        println("AA length: ", aa_len)
    end
    seq, aa_len, sp_ac, sp_id
end

read_prot(filename) = try
    open(parse_prot, filename)
catch e
    msg = sprint(showerror, e)
    println(msg)
end

println("Enter file containing SwissProt database entry to analyze: ")

seq, aa_len, sp_ac, sp_id = read_prot(chomp(readline()))

function to_outfile(seq, aa_len, sp_ac, sp_id)
    outfile = open("sprot.fsa", "w")
    println(outfile,">sp|", sp_ac, "|", sp_id)

    for i in 1:60:length(seq)
        if i+59 > length(seq)
            println(outfile, seq[i:end])
        else
            println(outfile, seq[i:(i+59)])
        end
    end

    close(outfile)
end

to_outfile(seq, aa_len, sp_ac, sp_id)

# Exercise 3
function read_fasta(f)
    headlines = String[]
    seqs = String[]
    seq = ""
    for line = eachline(f)
        ln = chomp(line)
        if ln[1] == '>'
            if seq != ""
                push!(seqs, seq)
                seq = ""
            end
            push!(headlines, ln)
        else
            seq *= ln
        end
    end
    push!(seqs, seq)
    headlines, seqs
end

open("data/revdna.fsa", "w") do outfile
    try
        headlines, seqs = open(read_fasta, "data/dna.fsa")
        for (headline, seq) in zip(headlines, seqs)
            if length(findfirst(r"[ATCG]+", seq)) != length(seq)
                throw(ArgumentError("The sequence does not entirely contain ATCG."))
            end
            println(outfile, "$headline ComplementStrand")
            revseq = reverse_complement(seq)
            for i in 1:60:length(revseq)
                if i + 59 > length(revseq)
                    println(outfile, revseq[i:end])
                else
                    println(outfile, revseq[i:(i + 59)])
                end
            end
        end
    catch err
        msg = sprint(showerror, err)
        print(msg)
        exit(1)
    end
end

# Exercise 4-9
read_data() = try
    if length(ARGS) == 1
        filename = ARGS[1]
    else
        filename = "data/data1.gb"
    end
    open(_read_data, filename)
catch e
    msg = sprint(showerror, e)
    println(msg)
    exit(1)
end

function _read_data(f)
    MED = String[]
    SEQ_flag = AA_flag = CDS_flag = false
    SEQ = AA = CDS = ""
    CDS_range = ""

    for line in eachline(f)
        if startswith(line, "//")
            SEQ_flag = false
        end
        if SEQ_flag
            SEQ_re = findall(r"([atcg]+)", line)
            SEQ_parts = [line[x] for x in SEQ_re]
            SEQ *= join(SEQ_parts)
        else
            ACC_re = match(r"^ACCESSION\s+(\S+)$", line)
            DEF_re = match(r"^DEFINITION\s+(.+)$", line)
            ORG_re = match(r"^\s+ORGANISM\s+(.+)$", line)
            MED_re = match(r"^\s+MEDLINE\s+(\S+)$", line)

            if ACC_re !== nothing
                println("Accession number: ", ACC_re[1])
            elseif DEF_re !== nothing
                println("Definition: ", DEF_re[1])
            elseif ORG_re !== nothing
                println("Organism: ", ORG_re[1])
            elseif MED_re !== nothing
                push!(MED, MED_re[1])
            end
        end

        if startswith(line, "ORIGIN")
            SEQ_flag = true
        end

        CDS_re = match(r"^\s+CDS\s+(\S+)$", line)
        AA_re = match(r"^\s+/translation=(\S+)$", line)

        if CDS_flag
            CDS_data = strip(line)
            if startswith(CDS_data, "/")
                CDS_flag = false
            else
                CDS_range *= CDS_data
            end
        end
        if CDS_re !== nothing
            CDS_range *= CDS_re[1]
            CDS_flag = true
        end

        if AA_flag
            AA_data = match(r"^\s+(\S+)$", line)
            if AA_data !== nothing
                AA *= AA_data[1]
            else
                AA_flag = false
            end
        end
        if AA_re !== nothing
            AA *= AA_re[1]
            AA_flag = true
        end
    end

    CDS_list = [split(CDS_range[x], "..") for x in findall(r"(\d+\.\.\d+)", CDS_range)]

    for x in CDS_list
        CDS *= SEQ[parse(Int64, x[1]):parse(Int64, x[2])]
    end

    println("MEDLINE article number(s): ", join(MED, ", "))
    println()
    println("Sequence: ", SEQ)
    println()
    println("Coding sequence (CDS): ", CDS)
    println()
    println("Amino acid sequence: ", AA)
end

read_data()


text:
# This is a common snippet to be inserted into other modules to enforce a
# docstring style. Do not use it as a module.

using DocStringExtensions

@template DEFAULT =
    """
    $(SIGNATURES)

    $(DOCSTRING)
    """


text:
################################################################################
#                          PARAMETER FUNCTION METHODS
################################################################################
# Extend dispatch_variable_ref
function dispatch_variable_ref(
    model::InfiniteModel,
    index::ParameterFunctionIndex
    )::ParameterFunctionRef
    return ParameterFunctionRef(model, index)
end

# Extend _add_data_object
function _add_data_object(
    model::InfiniteModel,
    object::ParameterFunctionData
    )::ParameterFunctionIndex
    return MOIUC.add_item(model.param_functions, object)
end

# Extend _data_dictionary (reference based)
function _data_dictionary(
    fref::ParameterFunctionRef
    )::MOIUC.CleverDict{ParameterFunctionIndex, ParameterFunctionData{<:ParameterFunction}}
    return JuMP.owner_model(fref).param_functions
end

# Extend _data_object
function _data_object(fref::ParameterFunctionRef)
  object = get(_data_dictionary(fref), JuMP.index(fref), nothing)
  object === nothing && error("Invalid parameter function reference, cannot find " *
                              "corresponding object in the model. This is likely " *
                              "caused by using the reference of a deleted function.")
  return object
end

# Extend _core_variable_object
function _core_variable_object(fref::ParameterFunctionRef)
    return _data_object(fref).func
end

# Extend _object_numbers
function _object_numbers(fref::ParameterFunctionRef)::Vector{Int}
    return _core_variable_object(fref).object_nums
end

# Extend _parameter_numbers
function _parameter_numbers(fref::ParameterFunctionRef)::Vector{Int}
    return _core_variable_object(fref).parameter_nums
end

"""
    build_parameter_function(
        _error::Function,
        func::Function,
        parameter_refs::Union{GeneralVariableRef, AbstractArray{<:GeneralVariableRef}, Tuple}
        )::ParameterFunction

Build an [`ParameterFunction`](@ref) object that employs a parameter function
`func` that takes instances of the infinite parameter(s) as input. This can
ultimately by incorporated into expressions to enable nonlinear infinite parameter
behavior and/or incorporate data over infinite domains.

Here `func` should be of the form:
```
func(paramvals...)::Float64
```
where the formatting of `paramvals` is analagous to point variables (and will be
based on the tuple of infinite parameter references given in `parameter_refs`).

Errors if the infinite parameter tuple is formatted incorrectly. The allowed
format follows that of infinite variables. Also errors if the function doesn't
accept a support realization of the `parameter_refs` as input.

**Example**
```julia-repl
julia> f = build_parameter_function(error, sin, t);
```
"""
function build_parameter_function(
    _error::Function,
    func::Function,
    parameter_refs::Union{GeneralVariableRef, AbstractArray{<:GeneralVariableRef}, Tuple};
    extra_kwargs...
    )
    # check for unneeded keywords
    for (kwarg, _) in extra_kwargs
        _error("Keyword argument $kwarg is not for use with parameter functions.")
    end
    # process the parameter reference inputs and check
    prefs = Collections.VectorTuple(parameter_refs)
    _check_parameter_tuple(_error, prefs)
    # check the function
    supp_type = typeof(Tuple(Vector{Float64}(undef, length(prefs)), prefs))
    if !hasmethod(func, supp_type)
        _error("Parameter function method `$(func)` is not defined for `$(func)(",
               join(Tuple(prefs), ", ") * ")`. Note that the infinite parameter ",
               "arguments `" * join(Tuple(prefs), ", ") * "` are checked via a ",
               "numeric support (each parameter is a `Float64`) of the same format.")
    end
    # get the parameter object numbers
    object_nums = Int[]
    for pref in prefs
        union!(object_nums, _object_number(pref))
    end
    # get the parameter numbers
    param_nums = [_parameter_number(pref) for pref in prefs]
    # make the variable and return
    return ParameterFunction(func, prefs, object_nums, param_nums)
end

# Fallback for weird macro inputs
function build_parameter_function(_error::Function, func, prefs; kwargs...)
    _error("Invalid syntax causing unexpected parsing of the function and " *
           "infinite parameter arguments.")
end

# Used to update parameter-parameter function mappings
function _update_param_var_mapping(
    fref::ParameterFunctionRef,
    prefs::Collections.VectorTuple
    )::Nothing
    for pref in prefs
        dependency_list = _parameter_function_dependencies(pref)
        if !(JuMP.index(fref) in dependency_list)
            push!(dependency_list, JuMP.index(fref))
        end
    end
    return
end

"""
    add_parameter_function(model::InfiniteModel, pfunc::ParameterFunction,
                           [name::String])::GeneralVariableRef

Add an [`ParameterFunction`](@ref) `pfunc` to the `model` using `name` for
printing and return a `GeneralVariableRef` such that it can be embedded in
expressions. Errors if the parameter function `pfunc` points to do not belong to
`model`. Note that `pfunc` should be created using [`build_parameter_function`](@ref).

**Example**
```julia-repl
julia> f = build_parameter_function(error, sin, t);

julia> fref = add_parameter_function(model, f)
sin(t)
```
"""
function add_parameter_function(
    model::InfiniteModel,
    pfunc::ParameterFunction,
    name::String = String(nameof(pfunc.func))
    )::GeneralVariableRef
    _check_parameters_valid(model, pfunc.parameter_refs)
    data_object = ParameterFunctionData(pfunc, name)
    findex = _add_data_object(model, data_object)
    fref = ParameterFunctionRef(model, findex)
    _update_param_var_mapping(fref, pfunc.parameter_refs)
    return _make_variable_ref(model, findex)
end

"""
    JuMP.name(fref::ParameterFunctionRef)::String

Extend `JuMP.name` to return the base name of `fref`.

**Example**
```julia-repl
julia> name(fref)
"func_name"
```
"""
function JuMP.name(fref::ParameterFunctionRef)::String
    object = get(_data_dictionary(fref), JuMP.index(fref), nothing)
    return object === nothing ? "" : object.name
end

"""
    JuMP.set_name(fref::ParameterFunctionRef, name::String)::Nothing

Extend `JuMP.set_name` to set the name of a parameter function.

**Example**
```julia-repl
julia> set_name(fref, "func_name")

julia> name(fref)
"func_name"
```
"""
function JuMP.set_name(fref::ParameterFunctionRef, name::String)::Nothing
    _data_object(fref).name = name
    return
end

"""
    parameter_function(func::Function,
                       pref_inputs::Union{GeneralVariableRef, AbstractArray{GeneralVariableRef}, Tuple};
                       [name::String = [the name of `func`]]
                       )::GeneralVariableRef

Make a parameter function and return a `GeneralVariableRef` that can be
embedded in InfiniteOpt expressions. This serves as a convenient wrapper for
[`build_parameter_function`](@ref) and [`add_parameter_function`](@ref). For an
even more convenient definition method see [`@parameter_function`](@ref).

Here `func` denotes the function that will take a support of infinite parameters as
input (formatted like `pref_inputs`) and will return a scalar value. Specifically,
`func` should be of the form:
```
func(paramvals...)::Float64
```
where the formatting of `paramvals` is analagous to point variables (and will be
based on the tuple of infinite parameter references given in `parameter_refs`).
Moreover, `func` must be a function that returns a scalar numeric value.

Errors if `func` will not take a support formatted like `pref_inputs` in
combination with the `fargs` and `fkwargs` specified. Also errors if `pref_inputs`
follow an invalid input format.

**Example**
```julia-repl
julia> p_func = parameter_function(sin, t)
sin(t)

julia> p_func3 = parameter_function((t_supp) -> 2 * sin(2 * t_supp), t, name = "mysin")
mysin(t)

julia> p_func4 = parameter_function(t, name = "mysin") do t_supp
                    if t_supp <= 5
                        return sin(t_supp)
                    else
                        return 2 * sin(2 * t_supp)
                    end
                 end

mysin(t)
```
"""
function parameter_function(
    func::Function,
    pref_inputs;
    name::String = String(nameof(func))
    )::GeneralVariableRef
    f = build_parameter_function(error, func, pref_inputs)
    model = JuMP.owner_model(first(f.parameter_refs))
    return add_parameter_function(model, f, name)
end

"""
    raw_parameter_refs(fref::ParameterFunctionRef)::VectorTuple

Return the raw [`VectorTuple`](@ref InfiniteOpt.Collections.VectorTuple) of the
parameter references that `fref` depends on. This is primarily an internal method
where [`parameter_refs`](@ref parameter_refs(::ParameterFunctionRef))
is intended as the preferred user function.
"""
function raw_parameter_refs(fref::ParameterFunctionRef)
    return _core_variable_object(fref).parameter_refs
end

"""
    parameter_refs(fref::ParameterFunctionRef)::Tuple

Return the parameter references associated with `fref`. This
is formatted as a Tuple of containing the parameter references as they inputted
to define `fref`.

**Example**
```julia-repl
julia> parameter_refs(p_func)
(t,)
```
"""
function parameter_refs(fref::ParameterFunctionRef)
    return Tuple(raw_parameter_refs(fref))
end

"""
    parameter_list(fref::ParameterFunctionRef)::Vector{GeneralVariableRef}

Return a vector of the parameter references that `fref` depends on. This is
primarily an internal method where [`parameter_refs`](@ref parameter_refs(::ParameterFunctionRef))
is intended as the preferred user function.
"""
function parameter_list(fref::ParameterFunctionRef)::Vector{GeneralVariableRef}
    return raw_parameter_refs(fref).values
end

"""
    raw_function(fref::ParameterFunctionRef)::Function

Returns the raw function behind `fref` that takes a particular support of `fref`'s
infinite parameters as input.
"""
function raw_function(fref::ParameterFunctionRef)
    return _core_variable_object(fref).func
end

"""
    call_function(fref::ParameterFunctionRef, support...)::Float64

Safely evaluates the [`raw_function`](@ref) of `fref` at a particular support `support`
point that matches the format of the infinite parameter tuple given when the `fref`
was defined. This is essentially equivalent to `raw_function(fref)(supps...)`.
"""
function call_function(fref::ParameterFunctionRef, supps...)::Float64
    pfunc = _core_variable_object(fref)
    return pfunc.func(supps...)
end

# Extend _semi_infinite_variable_dependencies
function _semi_infinite_variable_dependencies(
    fref::ParameterFunctionRef
     )::Vector{SemiInfiniteVariableIndex}
    return _data_object(fref).semi_infinite_var_indices
end

# Extend _derivative_dependencies
function _derivative_dependencies(
    fref::ParameterFunctionRef
     )::Vector{DerivativeIndex}
    return _data_object(fref).derivative_indices
end

# Extend _measure_dependencies
function _measure_dependencies(fref::ParameterFunctionRef)::Vector{MeasureIndex}
    return _data_object(fref).measure_indices
end

# Extend _constraint_dependencies
function _constraint_dependencies(
    fref::ParameterFunctionRef
    )::Vector{InfOptConstraintIndex}
    return _data_object(fref).constraint_indices
end

"""
    used_by_semi_infinite_variable(fref::ParameterFunctionRef)::Bool

Return a `Bool` indicating if `fref` is used by a semi-infinite infinite variable.

**Example**
```julia-repl
julia> used_by_semi_infinite_variable(fref)
false
```
"""
function used_by_semi_infinite_variable(fref::ParameterFunctionRef)::Bool
    return !isempty(_semi_infinite_variable_dependencies(fref))
end

"""
    used_by_derivative(fref::ParameterFunctionRef)::Bool

Return a `Bool` indicating if `fref` is used by a derivative.

**Example**
```julia-repl
julia> used_by_derivative(vref)
true
```
"""
function used_by_derivative(fref::ParameterFunctionRef)::Bool
    return !isempty(_derivative_dependencies(fref))
end

"""
    used_by_measure(fref::ParameterFunctionRef)::Bool

Return a `Bool` indicating if `fref` is used by a measure.

**Example**
```julia-repl
julia> used_by_measure(fref)
true
```
"""
function used_by_measure(fref::ParameterFunctionRef)::Bool
    return !isempty(_measure_dependencies(fref))
end

"""
    used_by_constraint(fref::ParameterFunctionRef)::Bool

Return a `Bool` indicating if `fref` is used by a constraint.

**Example**
```julia-repl
julia> used_by_constraint(fref)
false
```
"""
function used_by_constraint(fref::ParameterFunctionRef)::Bool
    return !isempty(_constraint_dependencies(fref))
end

"""
    is_used(fref::ParameterFunctionRef)::Bool

Return a `Bool` indicating if `fref` is used in the model.

**Example**
```julia-repl
julia> is_used(fref)
true
```
"""
function is_used(fref::ParameterFunctionRef)::Bool
    return used_by_measure(fref) || used_by_constraint(fref) ||
           used_by_semi_infinite_variable(fref) || used_by_derivative(fref)
end

"""
    JuMP.delete(model::InfiniteModel, fref::ParameterFunctionRef)::Nothing

Extend `JuMP.delete` to delete parameter functions and their dependencies. Errors
if `fref` is invalid, meaning it has already been deleted or it belongs to
another model.
"""
function JuMP.delete(model::InfiniteModel, fref::ParameterFunctionRef)::Nothing
    @assert JuMP.is_valid(model, fref) "Parameter function is invalid."
    # update the optimizer model status
    if is_used(fref)
        set_optimizer_model_ready(model, false)
    end
    # update parameter mapping
    all_prefs = parameter_list(fref)
    for pref in all_prefs
        filter!(e -> e != JuMP.index(fref), _parameter_function_dependencies(pref))
    end
    gvref = _make_variable_ref(model, JuMP.index(fref))
    # delete associated semi-infinite variables and mapping
    for index in _semi_infinite_variable_dependencies(fref)
        JuMP.delete(model, dispatch_variable_ref(model, index))
    end
    # delete associated derivative variables and mapping
    for index in _derivative_dependencies(fref)
        JuMP.delete(model, dispatch_variable_ref(model, index))
    end
    # remove from measures if used
    for mindex in _measure_dependencies(fref)
        mref = dispatch_variable_ref(model, mindex)
        func = measure_function(mref)
        data = measure_data(mref)
        if func isa GeneralVariableRef
            new_func = zero(JuMP.GenericAffExpr{Float64, GeneralVariableRef})
            new_meas = Measure(new_func, data, Int[], Int[], true)
        else
            _remove_variable(func, gvref)
            new_meas = build_measure(func, data)
        end
        _set_core_variable_object(mref, new_meas)
    end
    # remove from constraints if used
    for cindex in copy(_constraint_dependencies(fref))
        cref = _make_constraint_ref(model, cindex)
        func = JuMP.jump_function(JuMP.constraint_object(cref))
        if func isa GeneralVariableRef
            set = JuMP.moi_set(JuMP.constraint_object(cref))
            new_func = zero(JuMP.GenericAffExpr{Float64, GeneralVariableRef})
            new_constr = JuMP.ScalarConstraint(new_func, set)
            _set_core_constraint_object(cref, new_constr)
            empty!(_object_numbers(cref))
        elseif func isa AbstractArray && any(isequal(gvref), func)
            JuMP.delete(model, cref)
        else
            _remove_variable(func, gvref)
            _data_object(cref).object_nums = sort(_object_numbers(func))
        end
    end
    # delete the data object
    _delete_data_object(fref)
    return
end

################################################################################
#                               BASIC EXTENSIONS
################################################################################
# Convert to NLPExpr
function Base.convert(::Type{NLPExpr}, expr)
    return NLPExpr(_LCRST.Node(_process_child_input(expr)))
end
function Base.convert(::Type{NLPExpr}, expr::NLPExpr)
    return expr
end

# Redefine Base.isequal for UnorderedPair{GeneralVariableRef}
# This avoids symbolic conflicts with ==
function Base.isequal(
    p1::P,
    p2::P
    ) where {P <: JuMP.UnorderedPair{GeneralVariableRef}}
    return (isequal(p1.a, p2.a) && isequal(p1.b, p2.b)) ||
           (isequal(p1.a, p2.b) && isequal(p1.b, p2.a))
end

# Define Base.isequal to avoid default to ==
function Base.isequal(
    v::Union{GeneralVariableRef, JuMP.GenericAffExpr, JuMP.GenericQuadExpr, NLPExpr},
    w
    )
    return false
end
function Base.isequal(
    w,
    v::Union{GeneralVariableRef, JuMP.GenericAffExpr, JuMP.GenericQuadExpr, NLPExpr}
    )
    return false
end
function Base.isequal(
    w::Union{GeneralVariableRef, JuMP.GenericAffExpr, JuMP.GenericQuadExpr, NLPExpr},
    v::Union{GeneralVariableRef, JuMP.GenericAffExpr, JuMP.GenericQuadExpr, NLPExpr}
    )
    return false # relies on underlying extension
end

################################################################################
#                              VARIABLE ITERATION
################################################################################
## Create helper methods to interrogate the variables of an expr w/ a function
# Constant
function _interrogate_variables(interrogator::Function, c)
    return
end

# GeneralVariableRef
function _interrogate_variables(
    interrogator::Function,
    v::JuMP.AbstractVariableRef
    )
    interrogator(v)
    return
end

# AffExpr
function _interrogate_variables(
    interrogator::Function,
    aff::JuMP.GenericAffExpr
    )
    for (v, _) in aff.terms
        interrogator(v)
    end
    return
end

# QuadExpr
function _interrogate_variables(
    interrogator::Function,
    quad::JuMP.GenericQuadExpr
    )
    for (p, _) in quad.terms
        interrogator(p.a)
        interrogator(p.b)
    end
    _interrogate_variables(interrogator, quad.aff)
    return
end

# NLPExpr
function _interrogate_variables(interrogator::Function, nlp::NLPExpr)
    for n in AbstractTrees.Leaves(nlp.tree_root)
        _interrogate_variables(interrogator, _node_value(n.data))
    end
    return
end

# AbstractArray
function _interrogate_variables(interrogator::Function, arr::AbstractArray)
    for ex in arr
        _interrogate_variables(interrogator, ex)
    end
    return
end

################################################################################
#                            VARIABLE LIST MAKING
################################################################################
## Determine which variables are present in a function
# GeneralVariableRef
function _all_function_variables(f::GeneralVariableRef)
    return [f]
end

# GenericAffExpr
function _all_function_variables(f::JuMP.GenericAffExpr)
    return collect(keys(f.terms))
end

# GenericQuadExpr
function _all_function_variables(f::JuMP.GenericQuadExpr)
    vref_set = Set(keys(f.aff.terms))
    for (pair, _) in f.terms
        push!(vref_set, pair.a, pair.b)
    end
    return collect(vref_set)
end

# NLPExpr or array of expressions
function _all_function_variables(f::Union{NLPExpr, AbstractArray})
    vref_set = Set{GeneralVariableRef}()
    _interrogate_variables(v -> push!(vref_set, v), f)
    return collect(vref_set)
end

# Fallback
function _all_function_variables(f)
    error("`_all_function_variables` not defined for expression of type $(typeof(f)).")
end

################################################################################
#                        OBJECT/PARAMETER NUMBER METHODS
################################################################################
## Return the unique set of object numbers in an expression
# Dispatch fallback (--> should be defined for each non-empty variable type)
_object_numbers(v::DispatchVariableRef) = Int[]

# GeneralVariableRef
function _object_numbers(v::GeneralVariableRef)
    return _object_numbers(dispatch_variable_ref(v))
end

# Other
function _object_numbers(expr)
    obj_nums = Set{Int}()
    _interrogate_variables(v -> union!(obj_nums, _object_numbers(v)), expr)
    return collect(obj_nums)
end

## Return the unique set of parameter numbers in an expression
# Dispatch fallback (--> should be defined for each non-empty variable type)
_parameter_numbers(v::DispatchVariableRef) = Int[]

# GeneralVariableRef
function _parameter_numbers(v::GeneralVariableRef)
    return _parameter_numbers(dispatch_variable_ref(v))
end

# Other
function _parameter_numbers(expr)
    param_nums = Set{Int}()
    _interrogate_variables(v -> union!(param_nums, _parameter_numbers(v)), expr)
    return collect(param_nums)
end

################################################################################
#                             MODEL EXTRACTION METHODS
################################################################################
## Get the model from an expression
# Constant
function _model_from_expr(::Union{Number, Bool})
    return
end

# GeneralVariableRef
function _model_from_expr(expr::GeneralVariableRef)
    return JuMP.owner_model(expr)
end

# AffExpr
function _model_from_expr(expr::JuMP.GenericAffExpr)
    if isempty(expr.terms)
        return
    else
        return JuMP.owner_model(first(keys(expr.terms)))
    end
end

# QuadExpr
function _model_from_expr(expr::JuMP.GenericQuadExpr)
    result = _model_from_expr(expr.aff)
    if result !== nothing
        return result
    elseif isempty(expr.terms)
        return
    else
        return JuMP.owner_model(first(keys(expr.terms)).a)
    end
end

# NLPExpr
function _model_from_expr(expr::NLPExpr)
    for node in AbstractTrees.Leaves(expr.tree_root)
        result = _model_from_expr(_node_value(node.data))
        if result !== nothing
            return result
        end
    end
    return
end

# Vector{GeneralVariableRef}
function _model_from_expr(vrefs::Vector{GeneralVariableRef})
    if isempty(vrefs)
        return
    else
        return JuMP.owner_model(first(vrefs))
    end
end

# Fallback
function _model_from_expr(expr)
    error("`_model_from_expr` not defined for expr of type $(typeof(expr)).")
end

################################################################################
#                            VARIABLE REMOVAL BOUNDS
################################################################################
## Delete variables from an expression
# GenericAffExpr
function _remove_variable(f::JuMP.GenericAffExpr, vref::GeneralVariableRef)
    if haskey(f.terms, vref)
        delete!(f.terms, vref)
    end
    return
end

# GenericQuadExpr
function _remove_variable(f::JuMP.GenericQuadExpr, vref::GeneralVariableRef)
    _remove_variable(f.aff, vref)
    for (pair, _) in f.terms
        if isequal(pair.a, vref) || isequal(pair.b, vref)
            delete!(f.terms, pair)
        end
    end
    return
end

# Helper functions for NLP variable deletion
function _remove_variable_from_node(node, c, vref)
    return
end
function _remove_variable_from_node(node, n_vref::GeneralVariableRef, vref)
    if isequal(n_vref, vref)
        node.data = NodeData(0.0)
    end
    return
end
function _remove_variable_from_node(
    node,
    ex::Union{JuMP.GenericAffExpr, JuMP.GenericQuadExpr},
    vref
    )::Nothing
    _remove_variable(ex, vref)
    return
end

# NLPExpr
function _remove_variable(f::NLPExpr, vref::GeneralVariableRef)
    for node in AbstractTrees.Leaves(f.tree_root)
        _remove_variable_from_node(node, _node_value(node.data), vref)
    end
    return
end

# AbstractArray
function _remove_variable(arr::AbstractArray, vref::GeneralVariableRef)
    for f in arr
        _remove_variable(f, vref)
    end
    return
end

################################################################################
#                                MAPPING METHODS
################################################################################
"""
    map_expression(transform::Function,
                   expr::JuMP.AbstractJuMPScalar)::JuMP.AbstractJuMPScalar

Map and return a new expression of `expr` where each variable is transformed
via `transform`. This can be helpful for writing user extensions.
"""
function map_expression(transform::Function, v::JuMP.AbstractVariableRef)
    return transform(v)
end

# AffExpr
function map_expression(transform::Function, aff::JuMP.GenericAffExpr)
    return _MA.@rewrite(sum(c * transform(v)
                        for (c, v) in JuMP.linear_terms(aff)) +
                        JuMP.constant(aff))
end

# QuadExpr
function map_expression(transform::Function, quad::JuMP.GenericQuadExpr)
    return _MA.@rewrite(sum(c * transform(v1) * transform(v2)
                        for (c, v1, v2) in JuMP.quad_terms(quad)) +
                        map_expression(transform, quad.aff))
end

# Helper methods for mapping NLP trees
function _process_node(expr::NLPExpr)
    return expr.tree_root
end
function _process_node(expr)
    return _LCRST.Node(NodeData(expr))
end
function _map_expr_node(transform, data::JuMP.AbstractJuMPScalar)
    return _process_node(map_expression(transform, data))
end
function _map_expr_node(transform, data)
    return _LCRST.Node(NodeData(data))
end

# NLPExpr
function map_expression(transform::Function, nlp::NLPExpr)
    return NLPExpr(_map_tree(n -> _map_expr_node(transform, _node_value(n.data)),
                   nlp.tree_root))
end

################################################################################
#                             COEFFICIENT METHODS
################################################################################
## Modify linear coefficient of variable in expression
# GeneralVariableRef
function _set_variable_coefficient!(expr::GeneralVariableRef,
    var::GeneralVariableRef,
    coeff::Real
    )::JuMP.GenericAffExpr{Float64, GeneralVariableRef}
    # Determine if variable is that of the expression and change accordingly
    if isequal(expr, var)
        return Float64(coeff) * var
    else
        return expr + Float64(coeff) * var
    end
end

# GenericAffExpr
function _set_variable_coefficient!(expr::JuMP.GenericAffExpr{C, V},
    var::V,
    coeff::Real
    )::JuMP.GenericAffExpr{C, V} where {C, V <: GeneralVariableRef}
    # Determine if variable is in the expression and change accordingly
    if haskey(expr.terms, var)
        expr.terms[var] = coeff
        return expr
    else
        return expr + coeff * var
    end
end

# GenericQuadExpr
function _set_variable_coefficient!(expr::JuMP.GenericQuadExpr{C, V},
    var::V,
    coeff::Real
    )::JuMP.GenericQuadExpr{C, V} where {C, V <: GeneralVariableRef}
    # Determine if variable is in the expression and change accordingly
    if haskey(expr.aff.terms, var)
        expr.aff.terms[var] = coeff
        return expr
    else
        return expr + coeff * var
    end
end

# Fallback
function _set_variable_coefficient!(expr, var::GeneralVariableRef, coeff::Real)
    error("Unsupported function type for coefficient modification.")
end

## Query the coefficient of a variable
# GeneralVariableRef
function _affine_coefficient(
    func::GeneralVariableRef,
    var::GeneralVariableRef
    )::Float64
    return isequal(func, var) ? 1.0 : 0.0
end

# GenericAffExpr
function _affine_coefficient(
    func::GenericAffExpr,
    var::GeneralVariableRef
    )::Float64
    return get(func.terms, var, 0.0)
end

# GenericQuadExpr
function _affine_coefficient(
    func::GenericQuadExpr,
    var::GeneralVariableRef
    )::Float64
    return get(func.aff.terms, var, 0.0)
end

# Fallback
function _affine_coefficient(func, var::GeneralVariableRef)
    error("Unsupported function type for coefficient queries.")
end

################################################################################
#                         PARAMETER REFERENCE METHODS
################################################################################
## Return an element of a parameter reference tuple given the model and index
# IndependentParameterIndex
function _make_param_tuple_element(model::InfiniteModel,
    idx::IndependentParameterIndex,
    )::GeneralVariableRef
    return _make_parameter_ref(model, idx)
end

# DependentParametersIndex
function _make_param_tuple_element(model::InfiniteModel,
    idx::DependentParametersIndex,
    )::Vector{GeneralVariableRef}
    dpref = DependentParameterRef(model, DependentParameterIndex(idx, 1))
    num_params = _num_parameters(dpref)
    return [GeneralVariableRef(model, idx.value, DependentParameterIndex, i)
            for i in 1:num_params]
end

"""
    parameter_refs(expr)::Tuple

Return the tuple of parameter references that determine the infinite
dependencies of `expr`.

**Example**
```julia-repl
julia> parameter_refs(my_expr)
(t,)
```
"""
function parameter_refs(
    expr::Union{JuMP.GenericAffExpr, JuMP.GenericQuadExpr, NLPExpr}
    )::Tuple
    model = _model_from_expr(expr)
    if model === nothing
        return ()
    else
        obj_nums = sort!(_object_numbers(expr))
        obj_indices = _param_object_indices(model)[obj_nums]
        return Tuple(_make_param_tuple_element(model, idx) for idx in obj_indices)
    end
end

################################################################################
#                        VARIABLE ITERATION (IN PROGRESS)
################################################################################
# struct Variables{T}
#     expr::T
# end

# function Base.iterate(itr::Variables{<:JuMP.AbstractVariableRef})
#     return itr.tree_root, nothing
# end

# function Base.iterate(itr::Variables{<:JuMP.AbstractVariableRef}, ::Nothing)
#     return
# end

# Base.length(itr::Variables{<:JuMP.AbstractVariableRef}) = 1
# Base.eltype(::Variables{V}) where {V <: JuMP.AbstractVariableRef} = V

# function Base.iterate(itr::Variables{<:JuMP.GenericAffExpr})
#     out = iterate(itr.expr.terms)
#     return out === nothing ? out : (out[1][1], out[2])
# end

# function Base.iterate(itr::Variables{<:JuMP.GenericAffExpr}, state)
#     out = iterate(itr.expr.terms, state)
#     return out === nothing ? out : (out[1][1], out[2])
# end

# Base.length(itr::Variables{<:JuMP.GenericAffExpr}) = length(itr.expr.terms)
# Base.eltype(::Variables{JuMP.GenericAffExpr{C, V}}) where {C, V} = V

# mutable struct _QuadItrData{V <: JuMP.AbstractVariableRef}
#     state::Int
#     use_aff::Bool
#     has_prev_term::Bool
#     next_term::V
# end

# function Base.iterate(itr::Variables{JuMP.GenericQuadExpr{C, V}}) where {C, V}
#     out = iterate(itr.expr.terms)
#     if out === nothing
#         aff_out = iterate(itr.expr.aff.terms)
#         if aff_out === nothing
#             return
#         else
#             return aff_out[1][1], _QuadItrData{V}(aff_out[2], true, false, aff_out[1][1])
#         end
#     end
#     return out[1][1].a, _QuadItrData{V}(out[2], false, true, out[1][1].b)
# end

# function Base.iterate(itr::Variables{<:JuMP.GenericQuadExpr}, state)
#     if state.has_prev_term
#         state.has_prev_term = false
#         return state.next_term, state
#     elseif !state.use_aff
#         out = iterate(itr.expr.terms, state.state)
#         if out === nothing
#             aff_out = iterate(itr.expr.aff.terms)
#             if aff_out === nothing
#                 return
#             else
#                 state.use_aff = true
#                 state.state = aff_out[2]
#                 return aff_out[1][1], state
#             end
#         end
#         state.state = out[2]
#         state.has_prev_term = true
#         state.next_term = out[1][1].b
#         return out[1][1].a, state
#     else
#         out = iterate(itr.expr.aff.terms, state.state)
#         if out === nothing
#             return
#         else
#             state.state = out[2]
#             return out[1][1], state
#         end
#     end
# end

# function Base.length(itr::Variables{<:JuMP.GenericQuadExpr})
#     return 2 * length(itr.expr.terms) + length(itr.expr.aff.terms)
# end
# Base.eltype(::Variables{JuMP.GenericQuadExpr{C, V}}) where {C, V} = V

# mutable struct _NLPItrData
#     leaf_itr::AbstractTrees.Leaves{_LCRST.Node{NodeData}}
#     state::_LCRST.Node{NodeData}
#     has_internal::Bool
#     internal_itr::Any
#     internal_state::Any
# end

# function _process_itr(raw::GeneralVariableRef, state)
#     state.has_internal = false
#     return raw
# end

# function _process_itr(::Real, state)
#     out = iterate(state.leaf_itr, state.state)
#     out === nothing && return
#     state.state = out[2]
#     raw = _node_value(out[1].data)
#     return _process_itr(raw, state)
# end

# function _process_itr(raw, state)
#     itr = Variables(raw)
#     out = iterate(itr)
#     if out === nothing
#         new_out = iterate(state.leaf_itr, state.state)
#         new_out === nothing && return
#         state.state = new_out[2]
#         raw = _node_value(new_out[1].data)
#         return _process_itr(raw, state)
#     else
#         state.internal_itr = itr
#         state.has_internal = true
#         state.internal_state = out[2]
#         return out[1]
#     end
# end

# function Base.iterate(itr::Variables{NLPExpr})
#     leaf_itr = AbstractTrees.Leaves(itr.expr.tree_root)
#     out = iterate(leaf_itr)
#     out === nothing && return
#     state = _NLPItrData(leaf_itr, out[2], false, nothing, nothing)
#     raw = _process_itr(_node_value(out[1].data), state)
#     return raw === nothing ? raw : (raw, state)
# end

# function Base.iterate(itr::Variables{NLPExpr}, state)
#     if state.has_internal
#         int_out = iterate(state.internal_itr, state.internal_state)
#         if int_out !== nothing
#             state.internal_state = int_out[2]
#             return int_out[1], state
#         end
#     end
#     out = iterate(state.leaf_itr, state.state)
#     out === nothing && return
#     state.state = out[2]
#     raw = _process_itr(_node_value(out[1].data), state)
#     return raw === nothing ? raw : (raw, state)
# end

# Base.IteratorSize(::Variables{NLPExpr}) = Base.SizeUnknown()
# Base.eltype(::Variables{NLPExpr}) = GeneralVariableRef


text:
using RegERMs, Base.Test

# Seed random generator to avoid random effects in testing
srand(1)

X = [1 1; 2 2;  1 -1]
y = [-1; -1; 1]

@test_approx_eq_eps optimize(SVM(X, y), optimizer=:l_bfgs).theta [-0.1499984861841396,-0.3499984861803246] 1e-5
@test_approx_eq_eps optimize(SVM(X, y; λ=1.0), optimizer=:l_bfgs).theta [1.34394e-16,-1.0] 1e-5
@test_approx_eq_eps optimize(SVM(X, y; λ=10.1), optimizer=:l_bfgs).theta [1.34394e-16,-1.0] 1e-5
@test_approx_eq_eps optimize(SVM(X, y), optimizer=:sgd).theta [-0.1499984861841396,-0.3499984861803246] 5e-2
@test_approx_eq_eps optimize(SVM(X, y; λ=1.0), optimizer=:sgd).theta [1.34394e-16,-1.0] 5e-2
@test_approx_eq_eps optimize(SVM(X, y; λ=10.0), optimizer=:sgd).theta [1.34394e-16,-1.0] 5e-2
@test_approx_eq_eps optimize(SVM(X, y, L1RDAParameters(0.1,1.0,1.0)), optimizer=:l1_rda).theta [-14.811388300841895,-37.39005079444413] 5e-2
@test_approx_eq_eps optimize(SVM(X, y, L1RDAParameters(1.0,1.0,1.0)), optimizer=:l1_rda).theta [0.0,-8.866306299725341] 5e-2
@test_approx_eq_eps optimize(SVM(X, y, L1RDAParameters(10.0,1.0,1.0)), optimizer=:l1_rda).theta [0.0,0.0] 5e-2
model1 = optimize(SVM(X, y, L1RDAParameters(0.1,1.0,1.0)), optimizer=:l1_rda)
model2 = optimize(SVM(X, y, λ=10.0), optimizer=:sgd)

@test predict(model1, X) == y
@test predict(model2, X) == y
show(IOBuffer(), SVM(X, y))

@test_throws DimensionMismatch SVM(X', y)
@test_throws ArgumentError SVM(X, [3; 3; 2])
@test_throws ArgumentError optimize(SVM(X, y, L1RDAParameters(-0.1,1.0,1.0)), optimizer=:l1_rda)
@test_throws ArgumentError optimize(SVM(X, y, L1RDAParameters(0.1,-1.0,1.0)), optimizer=:l1_rda)
@test_throws ArgumentError optimize(SVM(X, y, L1RDAParameters(0.1,1.0,-1.0)), optimizer=:l1_rda)


text:
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Description
# ==============================================================================
#
#   Conversions between representations using Julia built-in system.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #

Base.convert(::Type{KeplerianElements{T}}, k::KeplerianElements) where T =
    KeplerianElements{T}(k.t, k.a, k.e, k.i, k.Ω, k.ω, k.f)

Base.convert(::Type{OrbitStateVector{T}}, sv::OrbitStateVector) where T =
    OrbitStateVector{T}(sv.t, sv.r, sv.v, sv.a)

Base.convert(::Type{KeplerianElements}, sv::OrbitStateVector{T}) where T =
    Base.convert(KeplerianElements{T}, sv)

function Base.convert(KT::Type{KeplerianElements{T}}, sv::OrbitStateVector) where T
    k = sv_to_kepler(sv)
    return convert(KT, k)
end

Base.convert(::Type{OrbitStateVector}, k::KeplerianElements{T}) where T =
    Base.convert(OrbitStateVector{T}, k)

function Base.convert(ST::Type{OrbitStateVector{T}}, k::KeplerianElements) where T
    sv = kepler_to_sv(k)
    return convert(ST, sv)
end


text:
module JDD

using Debugger
using HTTP
using HTTP.WebSockets
using JSON2
using Rematch

const INSTANCE = Ref{Tuple{String, Int}}()

include("rpcs.jl")
include("endpoints.jl")

function server(;host=HTTP.Sockets.localhost, port=9229)
  @info "Starting server on" host port
  INSTANCE[] = (string(host), port)

  ROUTER = HTTP.Router()
  HTTP.@register(ROUTER, "GET", "/json/version", version)
  HTTP.@register(ROUTER, "GET", "/json", main)
  HTTP.@register(ROUTER, "GET", "/json/list", main)
  HTTP.@register(ROUTER, "GET", "/ws", handle_rpcs)

  println("""
      Debugger started
        to begin debugging open: chrome-devtools://devtools/bundled/js_app.html?experiments=true&ws=$host:$port/ws
      """)
  HTTP.serve(ROUTER, host, port)
end

end # module


text:
module ImGuiGLFWBackend

using Libdl
using GLFW_jll
using LibCImGui

include("LibGLFW.jl")
using .LibGLFW

const GLFW_GET_CLIPBOARD_TEXT_FUNCPTR = Ref{Ptr{Cvoid}}(C_NULL)
const GLFW_SET_CLIPBOARD_TEXT_FUNCPTR = Ref{Ptr{Cvoid}}(C_NULL)

function __init__()
    glfwInit() == GLFW_TRUE && atexit(glfwTerminate)
    GLFW_GET_CLIPBOARD_TEXT_FUNCPTR[] = dlsym(dlopen(libglfw), :glfwGetClipboardString)
    GLFW_SET_CLIPBOARD_TEXT_FUNCPTR[] = dlsym(dlopen(libglfw), :glfwSetClipboardString)
end

function c_get(x::Ptr{NTuple{N,T}}, i) where {N,T}
    unsafe_load(Ptr{T}(x), Integer(i)+1)
end

function c_set!(x::Ptr{NTuple{N,T}}, i, v) where {N,T}
    unsafe_store!(Ptr{T}(x), v, Integer(i)+1)
end

include("macros.jl")

include("constants.jl")
export GlfwClientApi, GlfwClientApi_Unknown, GlfwClientApi_OpenGL, GlfwClientApi_Vulkan

include("callbacks.jl")

include("context.jl")
export set_default_glfw_window_size, set_default_glfw_window_title

include("platform.jl")
include("interface.jl")

end


text:
"""
RIntegratorNextStateReward: Estimate Rsas'
REstimateDummy: Do nothing - used for original SmallBackups implementation of Prioritized Sweeping
RIntegratorStateActionReward: Estimate Rsa (when reward is a function of current s-a pair and NOT s')
"""
struct RIntegratorNextStateReward
    nr_states::Int
    init_R::Float64
    Ns::Array{Int, 1}
    Rsum::Array{Float64, 1}
    R::Array{Float64, 1}
end
function RIntegratorNextStateReward(; nr_states = 10, init_R = 0.)
    Ns = zeros(Int, nr_states)
    Rsum = zeros(nr_states)
    R = zeros(nr_states) .+ init_R
    RIntegratorNextStateReward(nr_states, init_R, Ns, Rsum, R)
end
export RIntegratorNextStateReward
function updater!(learnerR::RIntegratorNextStateReward, s0, a0, s1, r)
    learnerR.Ns[s1] += 1
    learnerR.Rsum[s1] += r
    learnerR.R[s1] = learnerR.Rsum[s1] / learnerR.Ns[s1]
end

struct REstimateDummy end
REstimateDummy()
updater!(::REstimateDummy, s0, a0, s1, r) = nothing
export REstimateDummy

struct RIntegratorStateActionReward
    nr_states::Int
    nr_actions::Int
    init_R::Float64
    Nsa::Array{Int, 2}
    Rsum::Array{Float64, 2}
    R::Array{Float64, 2}
end
function RIntegratorStateActionReward(; nr_states = 10, nr_actions = 4, init_R = 0.)
    Nsa = zeros(Int, nr_actions, nr_states)
    Rsum = zeros(nr_actions, nr_states)
    R = zeros(nr_actions, nr_states) .+ init_R
    RIntegratorStateActionReward(nr_states, nr_actions, init_R, Nsa, Rsum, R)
end
export RIntegratorStateActionReward
function updater!(learnerR::RIntegratorStateActionReward, s0, a0, s1, r)
    learnerR.Nsa[a0, s0] += 1
    learnerR.Rsum[a0, s0] += r
    learnerR.R[a0, s0] = learnerR.Rsum[a0, s0] / learnerR.Nsa[a0, s0]
    # learnerR.R[a0, s0] -= learnerR.R[a0, s0] / learnerR.Nsa[a0, s0]
    # learnerR.R[a0, s0] += r / learnerR.Nsa[a0, s0]
end

export updater!


text:
"""
    AbstractMetadata{X,T}

Abstract supertype for all metadata wrappers.

Metadata wrappers allow tracking the contents and origin of metadata. This can
facilitate conversion between metadata types (for saving a file to a differenet format)
or simply saving data back to the same file type with identical metadata.

Using a wrapper instead of `Dict` or `NamedTuple` also lets us pass metadata
objects to [`set`](@ref) without ambiguity about where to put them.
"""
abstract type AbstractMetadata{X,T} end

const _MetadataContents =Union{AbstractDict,NamedTuple}
const AllMetadata = Union{AbstractMetadata,AbstractDict}

Base.get(m::AbstractMetadata, args...) = get(val(m), args...)
Base.getindex(m::AbstractMetadata, key) = getindex(val(m), Symbol(key))
Base.setindex!(m::AbstractMetadata, x, key) = setindex!(val(m), x, Symbol(key))
Base.haskey(m::AbstractMetadata, key) = haskey(val(m), Symbol(key))
Base.keys(m::AbstractMetadata) = keys(val(m))
Base.iterate(m::AbstractMetadata, args...) = iterate(val(m), args...)
Base.IteratorSize(m::AbstractMetadata) = Base.IteratorSize(val(m))
Base.IteratorEltype(m::AbstractMetadata) = Base.IteratorEltype(val(m))
Base.eltype(m::AbstractMetadata) = eltype(val(m))
Base.length(m::AbstractMetadata) = length(val(m))
Base.:(==)(m1::AbstractMetadata, m2::AbstractMetadata) = m1 isa typeof(m2) && val(m1) == val(m2)

"""
    Metadata <: AbstractMetadata

    Metadata{X}(val::Union{Dict,NamedTuple})
    Metadata{X}(pairs::Pair...) => Metadata{Dict}
    Metadata{X}(; kw...) => Metadata{NamedTuple}

General [`Metadata`](@ref) object. The `X` type parameter
categorises the metadata for method dispatch, if required.
"""
struct Metadata{X,T<:_MetadataContents} <: AbstractMetadata{X,T}
    val::T
end
Metadata(val::T) where {T<:_MetadataContents} = Metadata{Nothing,T}(val)
Metadata{X}(val::T) where {X,T<:_MetadataContents} = Metadata{X,T}(val)

# NamedTuple/Dict constructor
# We have to combine these because the no-arg method is overwritten by empty kw.
function (::Type{M})(ps...; kw...) where M <: Metadata
    if length(ps) > 0 && length(kw) > 0
        throw(ArgumentError("Metadata can be constructed with args of Pair to make a Dict, or kw for a NamedTuple. But not both."))
    end
    length(kw) > 0 ? M((; kw...)) : M(Dict(ps...))
end

ConstructionBase.constructorof(::Type{<:Metadata{X}}) where {X} = Metadata{X}

val(m::Metadata) = m.val

# Metadata nearly always contains strings, which break GPU compat.
# For now just remove everything, but we could strip all strings
# and replace Dict with NamedTuple.
# This might also be a problem with non-GPU uses of Adapt.jl where keeping
# the metadata is fine.
Adapt.adapt_structure(to, m::Metadata) = NoMetadata()

"""
    NoMetadata <: AbstractMetadata

    NoMetadata()

Indicates an object has no metadata. But unlike using `nothing`,
`get`, `keys` and `haskey` will still work on it, `get` always
returning the fallback argument. `keys` returns `()` while `haskey`
always returns `false`.
"""
struct NoMetadata <: AbstractMetadata{Nothing,NamedTuple{(),Tuple{}}} end

val(m::NoMetadata) = NamedTuple()

Base.keys(::NoMetadata) = ()
Base.haskey(::NoMetadata, args...) = false
Base.get(::NoMetadata, key, fallback) = fallback
Base.length(::NoMetadata) = 0

function Base.show(io::IO, mime::MIME"text/plain", metadata::Metadata{N}) where N
    print(io, "Metadata")
    if N !== Nothing
        print(io, "{")
        show(io, N)
        print(io, "}")
    end
    printstyled(io, " of "; color=:light_black)
    show(io, mime, val(metadata))
end

# Metadata utils

function metadatadict(dict)
    symboldict = Dict{Symbol,Any}()
    for (k, v) in dict
        symboldict[Symbol(k)] = v
    end
    symboldict
end

metadata(x) = NoMetadata()

units(x) = units(metadata(x))
units(m::NoMetadata) = nothing
units(m::Nothing) = nothing
units(m::Metadata) = get(m, :units, nothing)
units(m::AbstractDict) = get(m, :units, nothing)


text:
using Printf

include("gen_rmesh_exp.jl")
include("gen_drmesh_exp.jl")
include("integrate.jl")

function my_func(r)
    return exp(-0.2*r)
    #return exp(-0.2*r*r)
end

function main()
    r_min = 0.0
    r_max = 200.0
    a = 1e9
    N = 500

    rmesh = gen_rmesh_exp(r_min, r_max, a, N)
    drmesh = gen_drmesh_exp(r_min, r_max, a, N)

    f = my_func.(rmesh)
    s = integrate_trapz_7(f, drmesh)

    println("s = ", s)

end

main()


text:
# # clip-01.jl

# Load Julia packages (libraries) needed  for the snippets in chapter 0

using StatisticalRethinking

# ### snippet 3.1

PrPV = 0.95
PrPM = 0.01
PrV = 0.001
PrP = PrPV*PrV + PrPM*(1-PrV)
PrVP = PrPV*PrV / PrP

# End of `03/clip-01.jl`


text:
Pkg.add("StatsBase")
Pkg.add("Mocha")
Pkg.add("Compat")


text:
############
# BuildRef #
############

mutable struct BuildRef
    repo::String  # the build repo
    sha::String   # the build + status SHA
    vinfo::String # versioninfo() taken during the build
end

BuildRef(repo, sha) = BuildRef(repo, sha, "retrieving versioninfo() failed")

Base.copy(x::BuildRef) = BuildRef(x.repo, x.sha, x.vinfo)

function Base.:(==)(a::BuildRef, b::BuildRef)
    return (a.repo == b.repo &&
            a.sha == b.sha &&
            a.vinfo == b.vinfo)
end

Base.summary(build::BuildRef) = string(build.repo, SHA_SEPARATOR, snipsha(build.sha))

# reduce download bandwidth and time by keeping copies after the build
# TODO: intercept JLDOWNLOAD instead?
function sync_srcs!(fromdir, todir, link::Bool)
    for n in readdir(fromdir)
        endswith(n, ".tmp") && continue
        src = abspath(fromdir, n)
        dst = abspath(todir, n)
        if isfile(src) && !islink(src) && !ispath(dst)
            if link
                symlink(src, dst)
            else
                dsttmp = dst * ".tmp"
                cp(src, dsttmp, force=true)
                chmod(dsttmp, 0o444)
                Base.Filesystem.rename(dsttmp, dst, force=true) # mv is implemented badly in Base, so avoid it
            end
        end
    end
end

# if a PR number is included, attempt to build from the PR's merge commit
# TODO: re-use PkgEval's BinaryBuilder-based build?
function build_julia!(config::Config, build::BuildRef, logpath, prnumber::Union{Int,Nothing}=nothing)
    # make a temporary workdir for our build
    gid = parse(Int, readchomp(`id -g $(config.user)`))
    tmpdir = mktempdir(workdir(config))
    chown(tmpdir, -1, gid)
    chmod(tmpdir, 0o775)
    srcdir = joinpath(tmpdir, "julia")
    run(`sudo -n -u $(config.user) -- mkdir -m 775 $srcdir`)
    chmod(tmpdir, 0o555)

    mirrordir = joinpath(workdir(config), "mirrors", split(config.trackrepo, "/")...)
    mkpath(dirname(mirrordir), mode=0o755)
    mkpidlock(mirrordir * ".lock") do
        if ispath(mirrordir)
            run(setenv(`git fetch --quiet --all`; dir=mirrordir))
        else
            mkpath(mirrordir)
            gitclone!(config.trackrepo, mirrordir, `--mirror`)
        end
    end

    # clone/fetch the appropriate Julia version
    if prnumber !== nothing
        # clone from `trackrepo`, not `build.repo`, since that's where the merge commit is
        gitclone!(config.trackrepo, srcdir, `-c core.sharedRepository=group --reference $mirrordir --dissociate`; user=config.user)
        try
            run(setenv(`sudo -n -u $(config.user) -- git fetch --quiet origin +refs/pull/$(prnumber)/merge:`; dir=srcdir))
        catch
            # if there's not a merge commit on the remote (likely due to
            # merge conflicts) then fetch the head commit instead.
            run(setenv(`sudo -n -u $(config.user) -- git fetch --quiet origin +refs/pull/$(prnumber)/head:`; dir=srcdir))
        end
        run(setenv(`sudo -n -u $(config.user) -- git checkout --quiet --force FETCH_HEAD`; dir=srcdir))
        build.sha = readchomp(setenv(`sudo -n -u $(config.user) -- git rev-parse HEAD`; dir=srcdir))
    else
        gitclone!(build.repo, srcdir, `-c core.sharedRepository=group --reference $mirrordir --dissociate`; user=config.user)
        run(setenv(`sudo -n -u $(config.user) -- git checkout --quiet $(build.sha)`; dir=srcdir))
    end

    # set up logs for STDOUT and STDERR
    logname = string(build.sha, "_build")
    outfile = joinpath(logpath, string(logname, ".out"))
    errfile = joinpath(logpath, string(logname, ".err"))

    mirrordir1 = joinpath(workdir(config), "srccache", "deps")
    mirrordir2 = joinpath(workdir(config), "srccache", "stdlib")
    mkpath(mirrordir1, mode=0o755)
    mkpath(mirrordir2, mode=0o755)
    srccache1 = joinpath(srcdir, "deps", "srccache")
    srccache2 = joinpath(srcdir, "stdlib", "srccache")
    run(`sudo -n -u $(config.user) -- mkdir -m 775 $srccache2 $srccache1`)

    # TODO: support user build flags (like PkgEval)
    buildflags = ["JULIA_PRECOMPILE=0"]

    # run the build
    cpus = mycpus(config)
    sync_srcs!(mirrordir1, srccache1, true)
    sync_srcs!(mirrordir2, srccache2, true)
    run(pipeline(setenv(`sudo -n -u $(config.user) -- make -j$(length(cpus)) --output-sync=target $buildflags`; dir=srcdir), stdout=outfile, stderr=errfile))
    sync_srcs!(srccache1, mirrordir1, false)
    sync_srcs!(srccache2, mirrordir2, false)
    run(`sudo -n -u $(config.user) -- rm -rf $srccache2 $srccache1`) # erase the cloned files (which might have difficult permissions)
    run(`sudo -n -u $(config.user) -- chmod -R a-w $srcdir`) # make it r-x to all
    # TODO: symlink("bin/julia", joinpath(tmpdir, "julia"))
    return tmpdir
end


text:
using CPLEX
using Test

@testset "MathOptInterface Tests" begin
    for file in readdir("MathOptInterface")
        include(joinpath("MathOptInterface", file))
    end
end

@testset "Deprecated functions" begin
    err = ErrorException(CPLEX._DEPRECATED_ERROR_MESSAGE)
    @test_throws err newlongannotation()
    @test_throws err CPLEX.get_status()
    model = CPLEX.Optimizer()
    @test_throws err model.inner
    @test model.lp isa Ptr{Cvoid}
end


text:
module Feather

include("arrow/Arrow.jl")
using .Arrow
using Mmap, Dates
using FlatBuffers, CategoricalArrays, DataFrames, Tables

# in many cases we can load from a general IO, but if we get a string we assume it's a file
const FilenameOrIO = Union{<:IO,<:AbstractString,<:AbstractVector{UInt8}}

const FEATHER_VERSION = 2
# wesm/feather/cpp/src/common.h
const FEATHER_MAGIC_BYTES = Vector{UInt8}(codeunits("FEA1"))
const MIN_FILE_LENGTH = 12
const SHOULD_USE_MMAP = !Sys.iswindows()

include("metadata.jl")  # flatbuffer defintions
include("loaddata.jl")
include("source.jl")
include("sink.jl")

end # module


text:
using Test

@testset "SSPs Component" begin
    include("test_SSPs.jl")
end

@testset "RegionAggregatorSum Component" begin
    include("test_RegionAggregatorSum.jl")
end

@testset "Coupled" begin
    include("test_Coupled.jl")
end