text:
"""
    getARXregressor(y::AbstractVector,u::AbstractVecOrMat, na, nb)
Returns a shortened output signal `y` and a regressor matrix `A` such that the least-squares ARX model estimate of order `na,nb` is `y\\A`
Return a regressor matrix used to fit an ARX model on, e.g., the form
`A(z)y = B(z)f(u)`
with output `y` and input `u` where the order of autoregression is `na` and
the order of input moving average is `nb`
# Example
Here we test the model with the Function `f(u) = √(|u|)`
```julia
A     = [1,2*0.7*1,1] # A(z) coeffs
B     = [10,5] # B(z) coeffs
u     = randn(100) # Simulate 100 time steps with Gaussian input
y     = filt(B,A,u)
yr,A  = getARXregressor(y,u,3,2) # We assume that we know the system order 3,2
x     = A\\yr # Estimate model polynomials
plot([yr A*x], lab=["Signal" "Prediction"])
```
For nonlinear ARX-models, see [BasisFunctionExpansions.jl](https://github.com/baggepinnen/BasisFunctionExpansions.jl/). See also `arx`
"""
function getARXregressor(y::AbstractVector,u::AbstractVecOrMat, na, nb)
    length(nb) == size(u,2) || throw(ArgumentError("Length of nb must equal number of input signals"))
    m    = max(na,maximum(nb))+1 # Start of yr
    @assert m >= 1
    n    = length(y) - m + 1 # Final length of yr
    @assert n <= length(y)
    A    = toeplitz(y[m:m+n-1],y[m:-1:m-na])
    @assert size(A,2) == na+1
    y    = A[:,1] # extract yr
    A    = A[:,2:end]
    for i = 1:length(nb)
        nb[i] <= 0 && continue
        s = m-1
        A = [A toeplitz(u[s:s+n-1,i],u[s:-1:s-nb[i]+1,i])]
    end
    return y,A
end

function getARregressor(y::AbstractVector, na)
    m    = na+1 # Start of yr
    n    = length(y) - m + 1 # Final length of yr
    A    = toeplitz(y[m:m+n-1],y[m:-1:m-na])
    @assert size(A,2) == na+1
    y    = A[:,1] # extract yr
    A    = A[:,2:end]
    return y,A
end

@userplot Find_na
"""
    find_na(y::AbstractVector,n::Int)
Plots the RMSE and AIC For model orders up to `n`. Useful for model selection
"""
find_na
@recipe function find_na(p::Find_na)
    y,n = p.args[1:2]
    error = zeros(n,2)
    for i = 1:n
        yt,A = getARXregressor(y,0y,i,0)
        e = yt-A*(A\yt)
        error[i,1] = rms(e)
        error[i,2] = aic(e,i)
    end
    layout --> 2
    title --> ["RMS error" "AIC"]
    seriestype --> :scatter
    @series begin
        error
    end
end

@userplot Find_nanb
"""
    find_nanb(y::AbstractVector,u,na,nb)
Plots the RMSE and AIC For model orders up to `n`. Useful for model selection
"""
find_nanb
@recipe function find_nanb(p::Find_nanb; logrms=false)
    y,u,na,nb = p.args[1:4]
    error = zeros(na,nb, 2)
    for i = 1:na, j = 1:nb
        yt,A = getARXregressor(y,u,i,j)
        e = yt-A*(A\yt)
        error[i,j,1] = logrms ? log10.(rms(e)) : rms(e)
        error[i,j,2] = aic(e,i+j)
    end
    layout --> 2
    seriestype --> :heatmap
    xticks := (1:nb, 1:nb)
    yticks := (1:na, 1:na)
    ylabel := "na"
    xlabel := "nb"
    @series begin
        title := "RMS error"
        subplot := 1
        error[:,:,1]
    end
    @series begin
        title := "AIC"
        subplot := 2
        error[:,:,2]
    end
end

"""
    Gtf, Σ = arx(h,y, u, na, nb; λ = 0, estimator=\\)

Fit a transfer Function to data using an ARX model and equation error minimization.
`nb` and `na` are the length of the numerator and denominator polynomials. `h` is the sample time of the data. `λ > 0` can be provided for L₂ regularization. `estimator` defaults to \\ (least squares), alternatives are `estimator = tls` for total least-squares estimation. `arx(Δt,yn,u,na,nb, estimator=wtls_estimator(y,na,nb)` is potentially more robust in the presence of heavy measurement noise.
The number of free parameters is `na-1+nb`
`Σ` is the covariance matrix of the parameter estimate. See `bodeconfidence` for visualiztion of uncertainty.

Supports MISO estimation by supplying a matrix `u` where times is first dim, with nb = [nb₁, nb₂...]
"""
function arx(h,y::AbstractVector, u::AbstractVecOrMat, na, nb; λ = 0, estimator=\)
    all(nb .<= na) || throw(DomainError(nb,"nb must be <= na"))
    na >= 1 || throw(ArgumentError("na must be positive"))
    y_train, A = getARXregressor(y,u, na, nb)
    w = ls(A,y_train,λ,estimator)
    a,b = params2poly(w,na,nb)
    model = tf(b,a,h)
    Σ = parameter_covariance(y_train, A, w, λ)
    return model, Σ
end

function ar(h,y::AbstractVector, na; λ = 0, estimator=\)
    na >= 1 || throw(ArgumentError("na must be positive"))
    y_train, A = getARregressor(y, na)
    w = ls(A,y_train,λ,estimator)
    a,b = params2poly(w,na)
    model = tf(b,a,h)
    Σ = parameter_covariance(y_train, A, w, λ)
    return model, Σ
end

function ls(A,y,λ=0,estimator=\)
    if λ == 0
        w = estimator(A,y)
    else
        w = estimator([A; λ*I], [y;zeros(size(A,2))])
    end
    w
end

"""
G, Gn = plr(h,y,u,na,nb,nc; initial_order = 20)

Perform pseudo-linear regression to estimate a model on the form
`Ay = Bu + Cw`
The residual sequence is estimated by first estimating a high-order arx model, whereafter the estimated residual sequence is included in a second estimation problem. The return values are the estimated system model, and the estimated noise model. `G` and `Gn` will always have the same denominator polynomial.
"""
function plr(h,y,u,na,nb,nc; initial_order = 20)
    all(nb .<= na) || throw(DomainError(nb,"nb must be <= na"))
    na >= 1 || throw(ArgumentError("na must be positive"))
    na -= 1
    y_train, A = getARXregressor(y,u,initial_order,initial_order)
    w1 = A\y_train
    yhat = A*w1
    ehat = yhat - y_train
    ΔN = length(y)-length(ehat)
    y_train, A = getARXregressor(y[ΔN+1:end-1],[u[ΔN+1:end-1,:] ehat[1:end-1]],na,[nb;nc])
    w = A\y_train
    a,b = params2poly(w,na,nb)
    model = tf(b,a,h)
    c = w[na+sum(nb)+1:end]
    noise_model = tf(c,a,h)
    model, noise_model
end
# Helper constructor to make a MISO system after MISO arx estimation
function ControlSystems.tf(b::AbstractVector{<:AbstractVector{<:Number}}, a::AbstractVector{<:Number}, h)
    tfs = map(b) do b
        tf(b,a,h)
    end
    hcat(tfs...)
end

"""
wtls_estimator(y,na,nb)
Create an estimator function for estimation of arx models in the presence of measurement noise.
"""
function wtls_estimator(y,na,nb)
    rowQ = [diagm(0=>[ones(na);zeros(nb);1]) for i = 1:length(y)-(na)]
    Qaa,Qay,Qyy = rowcovariance(rowQ)
    (A,y)->wtls(A,y,Qaa,Qay,Qyy)
end

"""
    a,b = params2poly(params,na,nb)
Used to get numerator and denominator polynomials after arx fitting
"""
function params2poly(w,na,nb)
    a = [1; -w[1:na]]
    w = w[na+1:end]
    b = map(nb) do nb
        b = w[1:nb]
        w = w[nb+1:end]
        b
    end
    a,b
end
function params2poly(w,na)
    a = [1; -w]
    b = 1
    a,b
end
"""
w, a, b = params(G::TransferFunction)
w = [a;b]
"""
function params(G::TransferFunction)
    am, bm = -denpoly(G)[1].a[1:end-1], G.matrix[1].num.a
    wm     = [am; bm]
    wm, am, bm
end

"""
    Σ = parameter_covariance(y_train, A, w, λ=0)
"""
function parameter_covariance(y_train, A, w, λ=0)
    σ² = var(y_train .- A*w)
    iATA = if λ == 0
        inv(A'A)
    else
        ATA = A'A
        ATAλ = ATA + λ*I
        ATAλ\ATA/ATAλ
    end
    iATA = (iATA+iATA')/2
    Σ = σ²*iATA + sqrt(eps())*Matrix(LinearAlgebra.I,size(iATA))
end

"""
TransferFunction(T::Type{<:AbstractParticles}, G::TransferFunction, Σ, N=500)

Create a `TransferFunction` where the coefficients are `Particles` from [`MonteCarloMeasurements.jl`](https://github.com/baggepinnen/MonteCarloMeasurements.jl) that can represent uncertainty.
# Example
```julia
using MonteCarloMeasurements
Gls,Σls = arx(Δt,y,u,na,nb)
Glsp    = TransferFunction(Particles, Gls, Σls)
w       = exp10.(LinRange(-3,log10(π/Δt),100))
mag     = bode(Glsp,w)[1][:]
errorbarplot(w,mag,0.01; yscale=:log10, xscale=:log10, layout=3, subplot=1, lab="ls")

See full example [here](https://github.com/baggepinnen/MonteCarloMeasurements.jl/blob/master/examples/controlsystems.jl)
```
"""
function ControlSystems.TransferFunction(T::Type{<:MonteCarloMeasurements.AbstractParticles}, G::TransferFunction, Σ, N=500)
      wm, am, bm = ControlSystemIdentification.params(G)
      na,nb      = length(am), length(bm)
      p          = T(N, MvNormal(wm, Σ))
      a,b        = ControlSystemIdentification.params2poly(p,na,nb)
      arxtf      = tf(b,a,G.Ts)
end

function ControlSystems.TransferFunction(T::Type{<:MonteCarloMeasurements.AbstractParticles}, G::TransferFunction, p::AbstractMatrix)
      wm, am, bm = ControlSystemIdentification.params(G)
      na,nb      = length(am), length(bm)
      p          = T(p .+ wm')
      a,b        = ControlSystemIdentification.params2poly(p,na,nb)
      arxtf      = tf(b,a,G.Ts)
end

function DSP.filt(tf::ControlSystems.TransferFunction, y)
    b,a = numvec(tf)[], denvec(tf)[]
    filt(b,a,y)
end

"""
    bodeconfidence(arxtf::TransferFunction, Σ::Matrix, ω = logspace(0,3,200))
Plot a bode diagram of a transfer function estimated with [`arx`](@ref) with confidence bounds on magnitude and phase.
"""
bodeconfidence

@userplot BodeConfidence

@recipe function BodeConfidence(p::BodeConfidence)
    arxtfm = p.args[1]
    Σ      = p.args[2]
    ω      = length(p.args) >= 3 ? p.args[3] : exp10.(LinRange(-2,3,200))
    L      = cholesky(Hermitian(Σ)).L
    wm, am, bm = params(arxtfm)
    na,nb  = length(am), length(bm)
    mc     = 100
    res = map(1:mc) do _
        w             = L*randn(size(L,1)) .+ wm
        a,b           = params2poly(w,na,nb)
        arxtf         = tf(b,a,arxtfm.Ts)
        mag, phase, _ = bode(arxtf, ω)
        mag[:], phase[:]
    end
    magmc      = reduce(hcat, getindex.(res,1))
    phasemc    = reduce(hcat, getindex.(res,2))
    mag        = mean(magmc,dims=2)[:]
    phase      = mean(phasemc,dims=2)[:]
    # mag,phase,_ = bode(arxtfm, ω) .|> x->x[:]
    uppermag   = getpercentile(magmc,0.95)[:]
    lowermag   = getpercentile(magmc,0.05)[:]
    upperphase = getpercentile(phasemc,0.95)[:]
    lowerphase = getpercentile(phasemc,0.05)[:]
    layout := (2,1)

    @series begin
        subplot := 1
        title --> "ARX estimate"
        ylabel --> "Magnitude"
        fillrange := [lowermag, uppermag]
        yscale --> :log10
        xscale --> :log10
        alpha --> 0.3
        ω, mag
    end
    @series begin
        subplot := 2
        fillrange := [lowerphase, upperphase]
        ylabel --> "Phase [deg]"
        xlabel --> "Frequency [rad/s]"
        xscale --> :log10
        alpha --> 0.3
        ω, phase
    end
    nothing

end

"""
    getpercentile(x,p)

calculates the `p`th percentile along dim 2
"""
function getpercentile(mag,p)
    uppermag = mapslices(mag, dims=2) do magω
        sort(magω)[round(Int,length(magω)*p)]
    end
end


text:
function reorder_complex_schur{T<:FloatingPoint}(s::Array{T,2},q::Array{T,2})

  n = size(s,1)

  for ilast = n:-1:2

    mag = norm(s[ilast,ilast])
    icurrent = ilast

    for j = (ilast-1):-1:1

      if norm(s[j,j]) > mag
        mag = norm(s[j,j])
        icurrent = j
      end

    end

    for i = icurrent:(ilast-1)

      u = s[i:i+1,i+1]
      u[2] = u[2]-s[i,i]
      g = givens(u'[:])
      s[i:i+1,:] = g's[i:i+1,:]
      s[:,i:i+1] = s[:,i:i+1]*g
      q[i:i+1,:] = g'q[i:i+1,:]

    end

  end

  return s,q

end

function reorder_complex_schur{T<:FloatingPoint}(s::Array{Complex{T},2},q::Array{Complex{T},2})

  n = size(s,1)

  for ilast = n:-1:2

    mag = norm(s[ilast,ilast])
    icurrent = ilast

    for j = (ilast-1):-1:1

      if norm(s[j,j]) > mag
        mag = norm(s[j,j])
        icurrent = j
      end

    end

    for i = icurrent:(ilast-1)

      u = s[i:i+1,i+1]
      u[2] = u[2]-s[i,i]
      g = givens(u'[:])
      s[i:i+1,:] = g's[i:i+1,:]
      s[:,i:i+1] = s[:,i:i+1]*g
      q[i:i+1,:] = g'q[i:i+1,:]

    end

  end

  return s,q

end

function reorder_generalized_complex_schur{T<:FloatingPoint}(s::Array{T,2},t::Array{T,2},q::Array{T,2},z::Array{T,2})

  # This draws on draws on Kressner, D., (2006), "BLOCK ALGORITHMS FOR REORDERING STANDARD AND GENERALIZED SCHUR FORMS", LAPACK WORKING NOTE 171

  n = size(s,1)

  for ilast = n:-1:2

    if norm(t[ilast,ilast]) < 2*eps(T)*norm(t)
      mag = Inf
    else
      mag = norm(s[ilast,ilast]/t[ilast,ilast])
    end
    icurrent = ilast

    if mag != Inf
      for j = (ilast-1):-1:1

        if norm(t[j,j]) < 2*eps(T)*norm(t)
          mag = Inf
          icurrent = j
        elseif norm(s[j,j]/t[j,j]) > mag
          mag = norm(s[j,j]/t[j,j])
          icurrent = j
        end

      end
    end

    for i = icurrent:(ilast-1)

      s1 = s[i,i]
      s2 = s[i,i+1]
      s4 = s[i+1,i+1]
      t1 = t[i,i]
      t2 = t[i,i+1]
      t4 = t[i+1,i+1]

      if norm(t1) < 2*eps(T)*norm(t) && norm(s1) < 2*eps(T)*norm(s)
        error("Bz-A is not a regular matrix pencil: your model is not well posed")
      elseif norm(t4) < 2*eps(T)*norm(t) && norm(s4) < 2*eps(T)*norm(s)
        error("Bz-A is not a regular matrix pencil: your model is not well posed")
      elseif norm(t1) < 2*eps(T)*norm(t) && norm(t4) < 2*eps(T)*norm(s)
        # Both eigenvalues are infinite and do not need to be switched
        gq = eye(2)
        gz = eye(2)
      elseif norm(t1*s4-t4*s1) < 2*eps(T)*maximum([norm(s),norm(t)])
        # Both eigenvalues have equal modulus and do not need to be switched
        gq = eye(2)
        gz = eye(2)
      else
        h = [s1 -s4; t1 -t4]
        f = [s2; t2]
        p = inv(h)*f
        if norm(p[1]) < 2*eps(T)*norm(h) && norm(p[2]) < 2*eps(T)*norm(h)
          gq = [0.0 1.0; 1.0 0.0]
          gz = [0.0 1.0; 1.0 0.0]
        else
          u = [-p[2], complex(1.0)]
          (gq,r) = householder(u)
          u = [complex(1.0) p[1]]'
          (gz,r) = householder(u)
          gz = [gz[:,2] gz[:,1]]
        end
      end

      s[i:i+1,:] = gq's[i:i+1,:]
      t[i:i+1,:] = gq't[i:i+1,:]
      s[:,i:i+1] = s[:,i:i+1]*gz
      t[:,i:i+1] = t[:,i:i+1]*gz
      q[i:i+1,:] = gq'q[i:i+1,:]
      z[:,i:i+1] = z[:,i:i+1]*gz

    end

  end

  return s,t,q,z

end

function reorder_generalized_complex_schur{T<:FloatingPoint}(s::Array{Complex{T},2},t::Array{Complex{T},2},q::Array{Complex{T},2},z::Array{Complex{T},2})

  # This method draws on Kressner, D., (2006), "BLOCK ALGORITHMS FOR REORDERING STANDARD AND GENERALIZED SCHUR FORMS", LAPACK WORKING NOTE 171

  n = size(s,1)

  for ilast = n:-1:2

    if norm(t[ilast,ilast]) < 2*eps(T)*norm(t)
      mag = Inf
    else
      mag = norm(s[ilast,ilast]/t[ilast,ilast])
    end
    icurrent = ilast

    if mag != Inf
      for j = (ilast-1):-1:1

        if norm(t[j,j]) < 2*eps(T)*norm(t)
          mag = Inf
          icurrent = j
        elseif norm(s[j,j]/t[j,j]) > mag
          mag = norm(s[j,j]/t[j,j])
          icurrent = j
        end

      end
    end

    for i = icurrent:(ilast-1)

      s1 = s[i,i]
      s2 = s[i,i+1]
      s4 = s[i+1,i+1]
      t1 = t[i,i]
      t2 = t[i,i+1]
      t4 = t[i+1,i+1]

      if norm(t1) < 2*eps(T)*norm(t) && norm(s1) < 2*eps(T)*norm(s)
        error("Bz-A is not a regular matrix pencil: your model is not well posed")
      elseif norm(t4) < 2*eps(T)*norm(t) && norm(s4) < 2*eps(T)*norm(s)
        error("Bz-A is not a regular matrix pencil: your model is not well posed")
      elseif norm(t1) < 2*eps(T)*norm(t) && norm(t4) < 2*eps(T)*norm(s)
        # Both eigenvalues are infinite and do not need to be switched
        gq = eye(2)
        gz = eye(2)
      elseif norm(t1*s4-t4*s1) < 2*eps(T)*maximum([norm(s),norm(t)])
        # Both eigenvalues have equal modulus and do not need to be switched
        gq = eye(2)
        gz = eye(2)
      else
        h = [s1 -s4; t1 -t4]
        f = [s2; t2]
        p = inv(h)*f
        if norm(p[1]) < 2*eps(T)*norm(h) && norm(p[2]) < 2*eps(T)*norm(h)
          gq = [0.0 1.0; 1.0 0.0]
          gz = [0.0 1.0; 1.0 0.0]
        else
          u = [-p[2], complex(1.0)]
          (gq,r) = householder(u)
          u = [complex(1.0) p[1]]'
          (gz,r) = householder(u)
          gz = [gz[:,2] gz[:,1]]
        end
      end

      s[i:i+1,:] = gq's[i:i+1,:]
      t[i:i+1,:] = gq't[i:i+1,:]
      s[:,i:i+1] = s[:,i:i+1]*gz
      t[:,i:i+1] = t[:,i:i+1]*gz
      q[i:i+1,:] = gq'q[i:i+1,:]
      z[:,i:i+1] = z[:,i:i+1]*gz

    end

  end

  return s,t,q,z

end

#=

function reorder_complex_schur{T<:FloatingPoint}(s::Array{T,2},q::Array{T,2})

  # Alternative reordering using Householder transform.  This method draws on Kressner, D., (2006), "BLOCK ALGORITHMS FOR REORDERING
  # STANDARD AND GENERALIZED SCHUR FORMS", LAPACK WORKING NOTE 171

  n = size(s,1)

  for ilast = n:-1:2

    mag = norm(s[ilast,ilast])
    icurrent = ilast

    for j = (ilast-1):-1:1

      if norm(s[j,j]) > mag
        mag = norm(s[j,j])
        icurrent = j
      end

    end

    for i = icurrent:(ilast-1)

      s1 = s[i,i]
      s2 = s[i,i+1]
      s4 = s[i+1,i+1]

      if norm(s1-s4) < 2*eps(T)*norm(s)
        g = eye(2)
      else
        (g,r) = householder([-(s1-s4)^(-1)*s2,1.0])
      end

      s[i:i+1,:] = g's[i:i+1,:]
      s[:,i:i+1] = s[:,i:i+1]*g
      q[i:i+1,:] = g'q[i:i+1,:]

    end

  end

  return s,q

end

function reorder_complex_schur{T<:FloatingPoint}(s::Array{Complex{T},2},q::Array{Complex{T},2})

  # Alternative reordering using Householder transform.  This method draws on Kressner, D., (2006), "BLOCK ALGORITHMS FOR REORDERING
  # STANDARD AND GENERALIZED SCHUR FORMS", LAPACK WORKING NOTE 171

  n = size(s,1)

  for ilast = n:-1:2

    mag = norm(s[ilast,ilast])
    icurrent = ilast

    for j = (ilast-1):-1:1

      if norm(s[j,j]) > mag
        mag = norm(s[j,j])
        icurrent = j
      end

    end

    for i = icurrent:(ilast-1)

      s1 = s[i,i]
      s2 = s[i,i+1]
      s4 = s[i+1,i+1]

      if norm(s1-s4) < 2*eps(T)*norm(s)
        g = eye(2)
      else
        p = [-(s1-s4)^(-1)*s2,complex(1.0)]
        (g,r) = householder(p)
      end

      s[i:i+1,:] = g's[i:i+1,:]
      s[:,i:i+1] = s[:,i:i+1]*g
      q[i:i+1,:] = g'q[i:i+1,:]

    end

  end

  return s,q

end

=#

ext: jl

max_stars_repo_name: dingliumath/SolveDSGE

max_issues_repo_name: dingliumath/SolveDSGE

max_forks_repo_name: dingliumath/SolveDSGE


text:
list = readlines("day1_input.txt")

function get_two(list)
    years = parse.(Int, list)
    for i in years
        for j in years
            if i+j == 2020
                return i*j
            end
        end
    end
    return "Not found"
end

function get_two_ind(list)
    years = parse.(Int, list)
    for i in 1:length(years)
        for j in (i+1):length(years)
            if years[i]+years[j] == 2020
                return years[i]*years[j]
            end
        end
    end
    return "Not found"
end

function get_two_comprehension(list)
    years = parse.(Int, list)
    prod = [a*b for a in years for b in years if a+b==2020]
    if length(prod) > 1
        return prod[1]
    else
        return "Not found"
    end
end

function get_two_sorting(list)
    years = parse.(Int, list)
    list_inc = sort(years)
    list_dec = sort(years, rev = true)
    for a in list_inc
        for b in list_dec
            if a + b == 2020
                return a*b
            end
        end
    end
    return "Not found"
end

function findsol(list)
    years = parse.(Int, list)
    sorted=sort(years)
    sm,la=1,length(sorted)
    while sm<la
        cur=sorted[sm]+sorted[la]
        if cur<2020
            sm+=1
        elseif cur>2020
            la-=1
        else
            return sorted[sm]*sorted[la]
        end
    end
    error("dead end")
end

function get_three(list)
    years = parse.(Int, list)
    for i in years
        for j in years
            for k in years
                if i+j+k == 2020
                    return i*j*k
                end
            end
        end
    end
    return "Not found"
end

@show get_two(list) == 355875
@show get_two_ind(list) == 355875
@show get_two_comprehension(list) == 355875
@show get_two_sorting(list) == 355875
@show findsol(list) == 355875
@show get_three(list) == 140379120

#=
```julia
julia> @btime get_two(list)
  19.653 μs (2 allocations: 1.78 KiB)

julia> @btime get_two_ind(list);
  18.567 μs (2 allocations: 1.78 KiB)

julia> @btime get_two_comprehension(list);
  42.493 μs (5 allocations: 1.98 KiB)

julia> @btime get_two_sorting(list);
  14.270 μs (4 allocations: 5.31 KiB)

julia> @btime findsol(list);
  12.201 μs (3 allocations: 3.55 KiB)

julia> @btime get_three(list);
  437.935 μs (2 allocations: 1.78 KiB)
```
 =#

nothing


text:
function sierpinski(n::Integer, token::AbstractString="*")
    x = fill(token, 1, 1)
    for _ in 1:n
        t = fill(" ", size(x))
        x = [x x x; x t x; x x x]
    end
    return x
end

function printsierpinski(m::Matrix)
    for r in 1:size(m, 1)
        println(join(m[r, :]))
    end
end

sierpinski(2, "#") |> printsierpinski

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
using JuliaPerformanceVerification
using Base.Test

# write your own tests here
@test 1 == 1

iteration_cnt_range = [ 10^i for i in 1:7 ]
vector_size_range   = [ 10^i for i in 1:7 ]
vector_wins = trues(length(vector_size_range), length(iteration_cnt_range))
vector_faster = zeros(length(vector_size_range), length(iteration_cnt_range))
for (i, current_vector_size) in enumerate(vector_size_range), (j, iteration_cnt) in enumerate(iteration_cnt_range)
	if current_vector_size * iteration_cnt >= 10^10
	    continue
	end
	a = rand(current_vector_size)
	b = rand(current_vector_size)
	c_vec = zeros(current_vector_size)
	c_iter = zeros(current_vector_size)
	c_vec_oped = zeros(current_vector_size)
	c_iter_oped = zeros(current_vector_size)
	tic()
	vecaddbyvec!(a, b, c_vec, iteration_cnt)
	vector_time = toq();
	tic()
	vecaddbyiter!(a, b, c_iter, iteration_cnt)
	iteration_time = toq();
	@test c_vec == c_iter
	tic()
	vecaddbyvec_oped!(a, b, c_vec_oped, iteration_cnt)
	vector_oped_time = toq();
	tic()
	vecaddbyiter_oped!(a, b, c_iter_oped, iteration_cnt)
	iteration_oped_time = toq();
	@test c_vec == c_iter
	@test c_vec_oped == c_iter
	@test c_vec == c_iter_oped
	vector_wins[i, j] = iteration_time > vector_time ? true: false
	vector_faster[i, j] = iteration_time / vector_time
	if vector_wins[i, j]
		@printf "Two %d dim vectors adding for %d times, vector add wins %.03f faster!\n" current_vector_size iteration_cnt vector_faster[i, j]
	else
		@printf "Two %d dim vectors adding for %d times, iteration add wins %.03f faster!\n" current_vector_size iteration_cnt 1/vector_faster[i, j]
	end
	@printf "  @inbounds and @simd accelerate ratio Vector: %.03f; Iteration: %.03f\n\n" vector_time/vector_oped_time iteration_time/iteration_oped_time
end

@printf "The vector add wins:\n %s \n" vector_wins
@printf "The vector add is x times fast as iteration:\n %s \n" vector_faster


text:
using Reagents
using Reagents: Computed, Return, CAS, PostCommit

struct HMNode{T,R}
    next::R
    data::T
    HMNode{T,R}(next) where {T,R} = new{T,R}(next)
    HMNode{T,R}(next, data) where {T,R} = new{T,R}(next, data)
end

struct Deleted{T}
    value::T
end

const HMNodeRef{T} =
    Reagents.Ref{Union{Nothing,HMNode{T},Deleted{<:Union{Nothing,HMNode{T}}}}}

HMNode(next::HMNodeRef{T}) where {T} = HMNode{T,typeof(next)}(next)
HMNode(next::HMNodeRef{T}, data) where {T} = HMNode{T,typeof(next)}(next, data)

struct HMList{T,K,R<:HMNodeRef{T}}
    key::K
    head::R
    function HMList{T}(key = identity) where {T}
        head = HMNodeRef{T}(nothing)
        return new{T,typeof(key),typeof(head)}(key, head)
    end
end

Base.eltype(::Type{<:HMList{T}}) where {T} = T

function search(list::HMList{T}, x::T) where {T}
    needle = list.key(x)
    while true
        prev = list.head
        curr = prev[]
        while true
            curr === nothing && return false, prev, nothing, nothing
            next = curr.next[]
            prev[] === curr || break
            if next isa Deleted
                if Reagents.try!(CAS(prev, curr, next.value)) === nothing
                    break
                end
                curr = next.value
            else
                ckey = list.key(curr.data)
                if !isless(ckey, needle)
                    return isequal(ckey, needle), prev, curr, next
                end
                prev = curr.next
                curr = next
            end
        end
    end
end

pushing(list::HMList{T}) where {T} =
    Computed() do x
        found, prev, curr, _next = search(list, x)
        if found
            Return(false)
        else
            node = HMNode(HMNodeRef{T}(curr), x)
            CAS(prev, curr, node) ⨟ Return(true)
        end
    end

trydeleting(list::HMList) =
    Computed() do x
        found, prev, curr, next = search(list, x)
        if found
            deleted = Deleted(next)
            CAS(curr.next, next, deleted) ⨟ PostCommit() do _
                Reagents.try!(CAS(prev, curr, next))
            end
        else
            Return(false)
        end
    end

function Base.push!(list::HMList, x)
    pushing(list)(convert(eltype(list), x))
    return list
end

function Base.delete!(list::HMList, x)
    trydeleting(list)(convert(eltype(list), x))
    return list
end

Base.in(x, list::HMList) = search(list, convert(eltype(list), x))[1]


text:
@testset "Test basis" begin

using HDiscontinuousGalerkin
using Tensors

#Test Dubiner base
interpolation = Dubiner{2,RefTetrahedron,4}()
quad_rule = QuadratureRule{2,RefTetrahedron}(Strang(),5)
for j = 1:15
    @test integrate(x->(dubiner_basis(x[1],x[2],j) - value(interpolation, j, x))^2,quad_rule) < eps()
    @test integrate(x->sum((∇dubiner_basis(x[1],x[2],j) - gradient_value(interpolation, j, x)).^2),quad_rule) < eps()
end

# Test nodal points for Lagrange base
nodal_points, topology = get_nodal_points(RefTetrahedron, Val{2}, 1)
@test  nodal_points == [Vec{2, Float64}((0.0, 0.0)),
                        Vec{2, Float64}((1.0, 0.0)),
                        Vec{2, Float64}((0.0, 1.0))]
@test topology == Dict(0=>3,1=>0,2=>0)
nodal_points, topology = get_nodal_points(RefTetrahedron, Val{2}, 2)
@test nodal_points == [Vec{2, Float64}((0.0, 0.0)),
                        Vec{2, Float64}((1.0, 0.0)),
                        Vec{2, Float64}((0.0, 1.0)),
                        Vec{2, Float64}((0.5, 0.5)),
                        Vec{2, Float64}((0.0, 0.5)),
                        Vec{2, Float64}((0.5, 0.0))]
@test topology == Dict(0=>3,1=>3,2=>0)
nodal_points, topology = get_nodal_points(RefTetrahedron, Val{2}, 3)
@test length(nodal_points) == 10
@test topology == Dict(0=>3,1=>6,2=>1)
@test nodal_points[10] ≈ Vec{2, Float64}((1/3, 1/3))
# Test Lagrange base
function ref_value(i::Int, ξ::Vec{2})
    ξ_x = ξ[1]
    ξ_y = ξ[2]
    i == 2 && return ξ_x
    i == 3 && return ξ_y
    i == 1 && return 1. - ξ_x - ξ_y
    throw(ArgumentError("no shape function $i for interpolation $ip"))
end

interpolation = Lagrange{2,RefTetrahedron,1}()
ξ_ref = Vec{2}((0.5,0.5))
for j = 1:3
    @test abs(ref_value(j,ξ_ref) - value(interpolation,j,ξ_ref)) < eps(Float64)
    @test norm(gradient(ξ -> ref_value(j, ξ), ξ_ref) - gradient_value(interpolation,j,ξ_ref), Inf) < eps(Float64)
end

function ref_value1(i::Int, ξ::Vec{1})
    ξ_x = ξ[1]
    i == 1 && return 1-ξ_x
    i == 2 && return ξ_x
    throw(ArgumentError("no shape function $i for interpolation $ip"))
end
ξ_ref = Vec{1}((0.5,))
interpolation = Lagrange{1,RefTetrahedron,1}()
for j = 1:2
    @test abs(ref_value1(j,ξ_ref) - value(interpolation,j,ξ_ref)) < eps(Float64)
    @test norm(gradient(ξ -> ref_value1(j, ξ), ξ_ref) - gradient_value(interpolation,j,ξ_ref), Inf) < eps(Float64)
end

#Test Legendre basis
function ref_value(i::Int, x::Real)
    @assert i <= 3
    if i == 0
        return 1.0
    elseif i == 1
        return sqrt(3)*(2*x - 1)
    elseif i == 2
        return sqrt(5)*(6*x^2 - 6*x + 1)
    else
        return sqrt(7)*(2*x - 1)*(10*x^2 - 10*x + 1)
    end
end
function ref_dvalue(i::Int, x::Real)
    @assert i <= 3
    if i == 0
        return 0
    elseif i == 1
        return sqrt(3)*2
    elseif i == 2
        return sqrt(5)*(12*x - 6)
    else
        return sqrt(7)*12*(5*x^2 - 5*x + 1)
    end
end

interpolation = Legendre{1,RefTetrahedron,3}()
quad_rule = QuadratureRule{1,RefTetrahedron}(GaussLegendre(),3)
for i = 0:3
    @test integrate(x->(ref_value(i,x[1]) - value(interpolation, i+1, x))^2,quad_rule) < eps()
    @test integrate(x->(ref_dvalue(i,x[1]) - gradient_value(interpolation, i+1, x)[1])^2,quad_rule) < eps()
end

end


text:
#!jl # This is a title
# # [This is another title](@id custom_id)

# This is a description


text:
fname = tempname()

# known quirks

NCDatasets.Dataset(fname,"c") do ds

    ds.attrib["single_element"] = [1]

    ds.attrib["vector_of_chars"] = ['a','b','c']

    ds.attrib["single_char"] = 'a'

    # issue 12
    str = "some_string_with_null\0"
    ccall((:nc_put_att,NCDatasets.libnetcdf),Cint,
          (Cint,Cint,Cstring,NCDatasets.nc_type,Csize_t,Ptr{Nothing}),
          ds.ncid,NCDatasets.NC_GLOBAL,"attrib",NCDatasets.NC_CHAR,length(str),pointer.(str))

end

NCDatasets.Dataset(fname,"r") do ds
    # same behaviour in python netCDF4 1.3.1
    @test ds.attrib["single_element"] == 1

    # no chars in python and numpy.charray are not supported
    @test ds.attrib["vector_of_chars"] == "abc"

    # all chars are strings in python
    @test ds.attrib["single_char"] == "a"

    # issue 12
    @test ds.attrib["attrib"] == "some_string_with_null"
end


text:
export logrange

# like the log function, but if it is broadcasted over values that are close to
# a range (in log space), it converts the result to a range representation
logrange(x) = log(x)
PlotAxes.fn_prefix(x::typeof(logrange)) = "log"
function Base.Broadcast.broadcasted(f::typeof(logrange),vals)
  logvals = log.(ustrip.(vals))
  if std(diff(logvals)) < 1e-8
    range(first(logvals),last(logvals),length=length(logvals))
  else
    logvals
  end
end

function PlotAxes.asplotable(x::AuditorySpectrogram,args...;quantize=(100,128),
    kwds...)
  args = replace(collect(args),:freq => (:freq => logrange))

  asplotable(getcontents(x),args...;quantize=quantize,kwds...)
end

default_quantize(x) = (100,)
default_quantize(x,y) = (100,128,)
default_quantize(x,y,args...) = (100,128,fill(10,length(args))...)

function PlotAxes.asplotable(x::MetaAxisArray;kwds...)
  asplotable(x,axisnames(x)...;kwds...)
end
function PlotAxes.asplotable(x::MetaAxisArray,ax1,axes...;
  quantize=default_quantize(ax1,axes...),kwds...)
  # TODO: handle time or freq not being present
  args = (ax1,axes...)
  args = replace(collect(args),:freq => (:freq => logrange))
  infront = filter(!isnothing,indexin([:time,:freq => logrange],args))
  others = setdiff(1:length(args),infront)
  args = args[vcat(infront,others)]

  asplotable(getcontents(x),args...;quantize=quantize,kwds...)
end

max_stars_repo_name: haberdashPI/ShammaModel

max_issues_repo_name: haberdashPI/ShammaModel

max_forks_repo_name: haberdashPI/ShammaModel


text:
# Tests for Revenue DDF DEA Models
@testset "RevenueDDFDEAModel" begin

    # Test using Book data
    X = [1; 1; 1; 1; 1; 1; 1; 1];
    Y = [7 7; 4 8; 8 4; 3 5; 3 3; 8 2; 6 4; 1.5 5];
    P = [1 1; 1 1; 1 1; 1 1; 1 1; 1 1; 1 1; 1 1];

    # With Gy = :Ones
    dearevenueddfones = dearevenueddf(X, Y, P, Gy = :Ones)

    @test typeof(dearevenueddfones) == RevenueDDFDEAModel

    @test nobs(dearevenueddfones) == 8
    @test ninputs(dearevenueddfones) == 1
    @test noutputs(dearevenueddfones) == 2
    @test ismonetary(dearevenueddfones) == false

    @test efficiency(dearevenueddfones) == efficiency(dearevenueddfones, :Economic)
    @test efficiency(dearevenueddfones, :Economic)   ≈ [0.0; 1.0; 1.0; 3.0; 4.0; 2.0; 2.0; 3.75]
    @test efficiency(dearevenueddfones, :Technical)  ≈ [0.0; 0.0; 0.0; 2.5; 4.0; 0.0; 1.5; 2.875]
    @test efficiency(dearevenueddfones, :Allocative) ≈ [0.0; 1.0; 1.0; 0.5; 0.0; 2.0; 0.5; 0.875]
    @test normfactor(dearevenueddfones) ≈ [2.0; 2.0; 2.0; 2.0; 2.0; 2.0; 2.0; 2.0]

    dearevenueddfonesmon = dearevenueddf(X, Y, P, Gy = :Ones, monetary = true)

    @test ismonetary(dearevenueddfonesmon) == true

    @test efficiency(dearevenueddfonesmon) ≈ efficiency(dearevenueddfones) .* normfactor(dearevenueddfones)
    @test efficiency(dearevenueddfonesmon, :Technical) ≈ efficiency(dearevenueddfones, :Technical) .* normfactor(dearevenueddfones)
    @test efficiency(dearevenueddfonesmon, :Allocative) ≈ efficiency(dearevenueddfones, :Allocative) .* normfactor(dearevenueddfones)

    # With Gy = :Ones and CRS (Same result in this example)
    dearevenueddfonescrs = dearevenueddf(X, Y, P, Gy = :Ones, rts = :CRS)

    @test efficiency(dearevenueddfonescrs, :Economic)   ≈ [0.0; 1.0; 1.0; 3.0; 4.0; 2.0; 2.0; 3.75]
    @test efficiency(dearevenueddfonescrs, :Technical)  ≈ [0.0; 0.0; 0.0; 2.5; 4.0; 0.0; 1.5; 2.875]
    @test efficiency(dearevenueddfonescrs, :Allocative) ≈ [0.0; 1.0; 1.0; 0.5; 0.0; 2.0; 0.5; 0.875]

    # With Gx = :Observed
    dearevenueddfobs = dearevenueddf(X, Y, P, Gy = :Observed)

    @test efficiency(dearevenueddfobs, :Economic)   ≈ [0.0; 6/36; 6/36; 0.75; 4/3; 0.4; 0.4; 1.153846] atol = 1e-5
    @test efficiency(dearevenueddfobs, :Technical)  ≈ [0.0; 0.0; 0.0; 5/9; 4/3; 0.0; 3/11; 3/5] atol = 1e-5
    @test efficiency(dearevenueddfobs, :Allocative) ≈ [0.0; 1/6; 1/6; 7/36; 0.0; 0.4; 7/55; 0.553846] atol = 1e-5

    # With Gx = :Mean
    dearevenueddfmean = dearevenueddf(X, Y, P, Gy = :Mean)

    @test efficiency(dearevenueddfmean, :Economic)   ≈ [0.0; 0.203821; 0.203822; 0.611465; 0.815287; 0.407643; 0.407643; 0.764331] atol = 1e-5
    @test efficiency(dearevenueddfmean, :Technical)  ≈ [0.0; 0.0; 0.0; 0.517799; 0.802508; 0.0; 0.300940; 0.595469] atol = 1e-5
    @test efficiency(dearevenueddfmean, :Allocative) ≈ [0.0; 0.203822; 0.203822; 0.093666; 0.012779; 0.407643; 0.106703; 0.168862] atol = 1e-5

    # With Gx = :Monetary
    dearevenueddfmon = dearevenueddf(X, Y, P, Gy = :Monetary)

    @test efficiency(dearevenueddfmon, :Economic)   ≈ [0.0; 2.0; 2.0; 6.0; 8.0; 4.0; 4.0; 7.5]
    @test efficiency(dearevenueddfmon, :Technical)  ≈ [0.0; 0.0; 0.0; 5.0; 8.0; 0.0; 3.0; 5.75]
    @test efficiency(dearevenueddfmon, :Allocative) ≈ [0.0; 2.0; 2.0; 1.0; 0.0; 4.0; 1.0; 1.75]

    # With custom direction
    dearevenueddfcustom = dearevenueddf(X, Y, P, Gy = ones(size(Y)))
    @test efficiency(dearevenueddfcustom) ≈ efficiency(dearevenueddfones)

    # Print
    show(IOBuffer(), dearevenueddfones)

    # Test errors
    @test_throws DimensionMismatch dearevenueddf([1; 2 ; 3], [4 ; 5], [1; 1; 1] , Gy = :Ones) #  Different number of observations
    @test_throws DimensionMismatch dearevenueddf([1; 2; 3], [4; 5; 6], [1; 2; 3; 4], Gy = :Ones) # Different number of observation in prices
    @test_throws DimensionMismatch dearevenueddf([1 1; 2 2; 3 3 ], [4; 5; 6], [1 1 1; 2 2 2; 3 3 3], Gy = :Ones) # Different number of output prices and outputs
    @test_throws DimensionMismatch dearevenueddf([1; 2; 3] , [4; 5; 6], [1; 1; 1], Gy = [1 1; 2 2; 3 3]) # Different size of outputs direction
    @test_throws ArgumentError dearevenueddf([1; 2; 3], [4; 5; 6], [1; 2; 3], Gy = :Ones, rts = :Error) # Invalid returns to scale
    @test_throws ArgumentError dearevenueddf([1; 2; 3], [1; 2; 3], [1; 1; 1], Gy = :Error) # Invalid outputs direction

end


text:
function moi_01(optimizer)

    # original problem
    # min -4x -3y
    # s.t.
    #      2x + y <= 4
    #       x +2y <= 4
    #       x     >= 0
    #           y >= 0
    #
    # sol x = y = 4/3
    # obj = 28/3

    MOI.empty!(optimizer)
    @test MOI.is_empty(optimizer)

    # add 10 variables - only diagonal is relevant
    X = MOI.add_variables(optimizer, 2)

    c1 = MOIU.normalize_and_add_constraint(optimizer,
        MOI.ScalarAffineFunction([
            MOI.ScalarAffineTerm(2.0, X[1]),
            MOI.ScalarAffineTerm(1.0, X[2])
        ], 0.0), MOI.EqualTo(4.0))

    c2 = MOIU.normalize_and_add_constraint(optimizer,
        MOI.ScalarAffineFunction([
            MOI.ScalarAffineTerm(1.0, X[1]),
            MOI.ScalarAffineTerm(2.0, X[2])
        ], 0.0), MOI.EqualTo(4.0))

    b1 = MOIU.normalize_and_add_constraint(optimizer,
        MOI.ScalarAffineFunction([
            MOI.ScalarAffineTerm(1.0, X[1])
        ], 0.0), MOI.GreaterThan(0.0))

    b2 = MOIU.normalize_and_add_constraint(optimizer,
        MOI.ScalarAffineFunction([
            MOI.ScalarAffineTerm(1.0, X[2])
        ], 0.0), MOI.GreaterThan(0.0))

    MOI.set(optimizer,
        MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(),
        MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([-4.0, -3.0], [X[1], X[2]]), 0.0)
        )
    MOI.set(optimizer, MOI.ObjectiveSense(), MOI.MIN_SENSE)
    MOI.optimize!(optimizer)

    obj = MOI.get(optimizer, MOI.ObjectiveValue())

    @test obj ≈ -9.33333 atol = 1e-2

    Xr = MOI.get(optimizer, MOI.VariablePrimal(), X)

    @test Xr ≈ [1.3333, 1.3333] atol = 1e-2

end

function moi_02(optimizer, mode = BilevelJuMP.SOS1Mode())

    # min -4x -3y
    # s.t.
    # y = argmin_y y
    #      2x + y <= 4
    #       x +2y <= 4
    #       x     >= 0
    #           y >= 0
    #
    # sol: x = 2, y = 0
    # obj_upper = -8
    # obj_lower =  0

    ### Build upper level model

    upper = MOIU.Model{Float64}()
    @test MOI.is_empty(upper)

    X = MOI.add_variables(upper, 2)

    MOI.set(upper, MOI.VariableName(), X[1], "x_u")
    MOI.set(upper, MOI.VariableName(), X[2], "y_u")

    MOI.set(upper,
        MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(),
        MOI.ScalarAffineFunction(
            MOI.ScalarAffineTerm.([-4.0, -3.0], [X[1], X[2]]), 0.0)
        )
    MOI.set(upper, MOI.ObjectiveSense(), MOI.MIN_SENSE)

    ### Build lower level model

    lower = MOIU.Model{Float64}()
    @test MOI.is_empty(lower)

    X2 = MOI.add_variables(lower, 2)
    MOI.set(lower, MOI.VariableName(), X2[1], "x_l")
    MOI.set(lower, MOI.VariableName(), X2[2], "y_l")

    c1 = MOIU.normalize_and_add_constraint(lower,
        MOI.ScalarAffineFunction([
            MOI.ScalarAffineTerm(2.0, X2[1]),
            MOI.ScalarAffineTerm(1.0, X2[2])
        ], 0.0), MOI.LessThan(4.0))
    MOI.set(lower, MOI.ConstraintName(), c1, "l_c1")

    c2 = MOIU.normalize_and_add_constraint(lower,
        MOI.ScalarAffineFunction([
            MOI.ScalarAffineTerm(1.0, X2[1]),
            MOI.ScalarAffineTerm(2.0, X2[2])
        ], 0.0), MOI.LessThan(4.0))
    MOI.set(lower, MOI.ConstraintName(), c2, "l_c2")

    b1 = MOIU.normalize_and_add_constraint(lower,
        MOI.ScalarAffineFunction([
            MOI.ScalarAffineTerm(1.0, X2[1])
        ], 0.0), MOI.GreaterThan(0.0))
    MOI.set(lower, MOI.ConstraintName(), b1, "l_b1")

    b2 = MOIU.normalize_and_add_constraint(lower,
        MOI.ScalarAffineFunction([
            MOI.ScalarAffineTerm(1.0, X2[2])
        ], 0.0), MOI.GreaterThan(0.0))
    MOI.set(lower, MOI.ConstraintName(), b2, "l_b12")

    MOI.set(lower,
        MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(),
        MOI.ScalarAffineFunction(
            MOI.ScalarAffineTerm.([1.0], [X2[2]]), 0.0)
        )
    MOI.set(lower, MOI.ObjectiveSense(), MOI.MIN_SENSE)

    ### Build bilevel model combining both models

    links = Dict(X[1] => X2[1], X[2] => X2[2])

    parametric = [X2[1]]

    blp, _, _ = BilevelJuMP.build_bilevel(upper, lower, links, parametric, mode)

    # MOI.optimize!(bridged)

    # obj = MOI.get(bridged, MOI.ObjectiveValue())

    # @test obj ≈ -9.33333 atol = 1e-2

    # Xr = MOI.get(bridged, MOI.VariablePrimal(), X)

    # @test Xr ≈ [1.3333, 1.3333] atol = 1e-2

    # lp_model = MathOptFormat.MOF.Model()
    # MOI.copy_to(lp_model, blp)
    # MOI.write_to_file(lp_model, "my_model.LP")
    # @show pwd()

    MOI.empty!(optimizer)
    @test MOI.is_empty(optimizer)
    MOI.copy_to(optimizer, blp)

    MOI.optimize!(optimizer)

end

function moi_03(optimizer, mode = BilevelJuMP.SOS1Mode())

    upper = MOIU.Model{Float64}()
    @test MOI.is_empty(upper)

    x = MOI.add_variable(upper)
    y = MOI.add_variable(upper)

    MOI.set(upper, MOI.VariableName(), x, "x_u")
    MOI.set(upper, MOI.VariableName(), y, "y_u")

    MOI.set(upper,
        MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(),
        MOI.ScalarAffineFunction(
            MOI.ScalarAffineTerm.([1.0, 3.0], [x, y]), 0.0)
        )
    MOI.set(upper, MOI.ObjectiveSense(), MOI.MIN_SENSE)

    lower = MOIU.Model{Float64}()
    @test MOI.is_empty(lower)

    x2 = MOI.add_variable(lower)
    y2 = MOI.add_variable(lower)

    MOI.set(lower, MOI.VariableName(), x2, "x_l")
    MOI.set(lower, MOI.VariableName(), y2, "y_l")

    c1 = MOIU.normalize_and_add_constraint(lower,
        MOI.ScalarAffineFunction([
            MOI.ScalarAffineTerm(1.0, x2),
            MOI.ScalarAffineTerm(1.0, y2)
        ], 0.0), MOI.LessThan(8.0))
    MOI.set(lower, MOI.ConstraintName(), c1, "l_c1")

    c2 = MOIU.normalize_and_add_constraint(lower,
        MOI.ScalarAffineFunction([
            MOI.ScalarAffineTerm(1.0, x2),
            MOI.ScalarAffineTerm(4.0, y2)
        ], 0.0), MOI.GreaterThan(8.0))
    MOI.set(lower, MOI.ConstraintName(), c2, "l_c2")

    c3 = MOIU.normalize_and_add_constraint(lower,
    MOI.ScalarAffineFunction([
        MOI.ScalarAffineTerm(1.0, x2),
        MOI.ScalarAffineTerm(2.0, y2)
    ], 0.0), MOI.LessThan(13.0))
    MOI.set(lower, MOI.ConstraintName(), c3, "l_c3")

    b1 = MOIU.normalize_and_add_constraint(lower,
        MOI.ScalarAffineFunction([
            MOI.ScalarAffineTerm(1.0, x2)
        ], 0.0), MOI.GreaterThan(1.0))
    MOI.set(lower, MOI.ConstraintName(), b1, "l_b1")

    b2 = MOIU.normalize_and_add_constraint(lower,
        MOI.ScalarAffineFunction([
            MOI.ScalarAffineTerm(1.0, x2)
        ], 0.0), MOI.LessThan(6.0))
    MOI.set(lower, MOI.ConstraintName(), b2, "l_b12")

    MOI.set(lower,
        MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(),
        MOI.ScalarAffineFunction(
            MOI.ScalarAffineTerm.([-1.0], [y2]), 0.0)
        )
    MOI.set(lower, MOI.ObjectiveSense(), MOI.MIN_SENSE)

    links = Dict(x => x2, y => y2)

    parametric = [x2]

    blp, _, _ = BilevelJuMP.build_bilevel(upper, lower, links, parametric, mode)

    # MOI.optimize!(bridged)

    # obj = MOI.get(bridged, MOI.ObjectiveValue())

    # @test obj ≈ -9.33333 atol = 1e-2

    # Xr = MOI.get(bridged, MOI.VariablePrimal(), X)

    # @test Xr ≈ [1.3333, 1.3333] atol = 1e-2

    # lp_model = MathOptFormat.MOF.Model()
    # MOI.copy_to(lp_model, blp)
    # MOI.write_to_file(lp_model, "my_model.LP")
    # @show pwd()

    MOI.empty!(optimizer)
    @test MOI.is_empty(optimizer)
    MOI.copy_to(optimizer, blp)

    MOI.optimize!(optimizer)

end


text:
"""
    A = JopSlantStack(dom[; dz=10.0, dh=10.0, h0=-1000.0, ...])

where `A` is the 2D slant-stack operator mapping for `z-h` to `tau-p`.  The domain of the operator
is `nz` x `nh` with precision T, `dz` is the depth spacing, `dh` is the half-offset spacing, and `h0` is the
origin of the half-offset axis.  The additional named optional arguments along with their default
values are,

* `theta=-60:1.0:60` - range of opening angles.  The ray parameter is: p=sin(theta/2)/c
* `padz=0.0,padh=0.0` - padding in depth and offset to apply before applying the Fourier transfrom
* `taperz=(0,0)` - beginning and end taper in the z-direction before transforming from `z-h` to `kz-kh`
* `taperh=(0,0)` - beginning and end taper in the h-direction before transforming from `z-h` to `kz-kh`
* `taperkz=(0,0)` - beginning and end taper in the kz-direction before transforming from `kz-kh` to `z-h`
* `taperkh=(0,0)` - beginning and end taper in the kh-direction before transforming from `kz-kh` to `z-h`

# Notes
* It should be possible to extend this operator to 3D

* If your domain is t-h rather than z-h, you can still use this operator
"""
function JopSlantStack(
        dom::JetAbstractSpace{T};
        theta = collect(-60.0:1.0:60.0),
        dz = -1.0,
        dh = -1.0,
        h0 = NaN,
        padz = 0.0,
        padh = 0.0,
        taperz = (0,0),
        taperh = (0,0),
        taperkz = (0,0),
        taperkh = (0,0)) where {T}
    dz < 0.0 && error("expected dz>0.0, got dz=$(dz)")
    dh < 0.0 && error("expected dh>0.0, got dh=$(dh)")
    isnan(h0) && error("expected finite h0, got h0=$(h0)")

    nz,nh = size(dom)

    # kz
    nzfft = nextprod([2,3,5,7], round(Int, nz*(1 + padz)))
    kn = pi/dz
    dk = kn/nzfft
    kz = dk*[0:div(nzfft,2)+1;]

    # kh
    nhfft = nextprod([2,3,5,7], round(Int, nh*(1+padh)))
    kn = pi/dh
    dk = kn/nhfft
    local kh
    if rem(nhfft,2) == 0
        kh = 2*dk*[ [0:div(nhfft,2);] ; [-div(nhfft,2)+1:1:-1;] ] # factor of 2 is for kg+ks with kg=ks -- i.e. go from (kg,ks)->kh
    else
        kh = 2*dk*[ [0:div(nhfft,2);] ; [-div(nhfft,2)+0:1:-1;] ] # factor of 2 is for kg+ks with kg=ks -- i.e. go from (kg,ks)->kh
    end

    # c*p
    cp = sin.(.5*theta*pi/180) # factor of .5 is to make theta the opening angle -- i.e. go from (theta_s, theta_g)->theta

    # conversions
    kz,kh,cp = map(x->convert(Array{T,1}, x), (kz,kh,cp))
    nzfft,nhfft = map(x->convert(Int64, x), (nzfft,nhfft))
    h0 = T(h0)

    # tapers
    TX = JopTaper(dom, (1,2), (taperz[1],taperh[1]), (taperz[2],taperh[2]))
    TK = JopTaper(JetSpace(Complex{eltype(dom)},div(nzfft,2)+1,length(cp)), (1,2), (taperkz[1], taperkh[1]), (taperkz[2], taperkh[2]), mode=(:normal,:fftshift))

    JopLn(dom = dom, rng = JetSpace(T, nz, length(cp)), df! = JopSlantStack_df!, df′! = JopSlantStack_df′!,
        s = (nzfft=nzfft, nhfft=nhfft, kz=kz, kh=kh, cp=cp, h0=h0, TX=TX, TK=TK))
end
export JopSlantStack

function JopSlantStack_df!(d::AbstractArray{T,2}, m::AbstractArray{T,2}; nzfft, nhfft, kz, kh, cp, h0, TX, TK, kwargs...) where {T}
    nz, nh, np, dh = size(m)..., length(cp), abs(kh[2]-kh[1])

    mpad = zeros(T, nzfft, nhfft)
    mpad[1:nz,1:nh] = TX*m

    M = rfft(mpad)
    dtmp = similar(d)

    D = zeros(eltype(M), size(M,1), np)
    for ikz = 1:div(nzfft,2)+1, ip = 1:np
        ikh_m1, ikh_p1, _kh = slantstack_compute_kh(ikz, ip, cp, kz, kh, nhfft)

        ikh_m1 < 1 && continue

        if ikh_m1 == ikh_p1
            D[ikz,ip] = M[ikz,ikh_p1]*exp(-im*kh[ikh_p1]*h0)
            continue
        end

        local d_p1, a_p1
        if 1 <= ikh_p1 <= nhfft
            d_p1 = M[ikz,ikh_p1]*exp(-im*kh[ikh_p1]*h0)
            a_p1 = abs(kh[ikh_p1] - _kh)/dh
        else
            a_p1 = 0.0
        end

        local d_m1, a_m1
        if 1 <= ikh_m1 <= nhfft
            d_m1 = M[ikz,ikh_m1]*exp(-im*kh[ikh_m1]*h0)
            a_m1 = abs(_kh - kh[ikh_m1])/dh
        else
            a_m1 = 0.0
        end

        D[ikz,ip] = a_m1*d_m1 + a_p1*d_p1
    end

    d .= brfft(TK*D, nzfft, 1)[1:nz,1:np] ./ nzfft
end

function JopSlantStack_df′!(m::AbstractArray{T,2}, d::AbstractArray{T,2}; nzfft, nhfft, kz, kh, cp, h0, TX, TK, kwargs...) where {T}
    nz, nh, np, dh = size(m)..., length(cp), abs(kh[2]-kh[1])

    dpad = zeros(T, nzfft, np)
    dpad[1:nz,:] = d

    D = TK * (rfft(dpad, 1) ./ nzfft)

    M = zeros(Complex{T}, div(nzfft,2)+1, nhfft)
    for ikz = 1:div(nzfft,2)+1, ip = 1:np
        ikh_m1, ikh_p1, _kh = slantstack_compute_kh(ikz, ip, cp, kz, kh, nhfft)

        ikh_m1 < 1 && continue

        if ikh_m1 == ikh_p1
            M[ikz,ikh_p1] += D[ikz,ip]*exp(im*kh[ikh_p1]*h0)
            continue
        end

        if 1 <= ikh_p1 <= nhfft
            m_p1 = D[ikz,ip]*exp(im*kh[ikh_p1]*h0)
            a_p1 = (kh[ikh_p1] - _kh)/dh
            M[ikz,ikh_p1] += a_p1*m_p1
        end

        if 1 <= ikh_m1 <= nhfft
            m_m1 = D[ikz,ip]*exp(im*kh[ikh_m1]*h0)
            a_m1 = (_kh - kh[ikh_m1])/dh
            M[ikz,ikh_m1] += a_m1*m_m1
        end
    end
    m .= TX * (brfft(M, nzfft)[1:nz,1:nh])
end

@inline function slantstack_compute_kh(ikz::Int64, ip::Int64, cp, kz, kh, nhfft)
    _kh = 1/(1-cp[ip]^2) - 1
    _kh < 0 && return -1,-1,1.0
    _kh = 2*kz[ikz]*sqrt(_kh)

    ikh_m1 = floor(Int64, _kh/kh[2]) + 1
    ikh_p1 = ceil(Int64, _kh/kh[2]) + 1
    ikh_m1 = ikh_m1 < 1 ? nhfft + ikh_m1 : ikh_m1
    ikh_p1 = ikh_p1 < 1 ? nhfft + ikh_p1 : ikh_p1

    ikh_m1, ikh_p1, _kh
end


text:
"""
    ASB07{N, AM, RM, S, R} <: AbstractContinuousPost

Implementation of Althoff - Stursberg - Buss algorithm for reachability of
linear systems with uncertain parameters and inputs using zonotopes.

## Fields

- `δ`                -- step-size of the discretization
- `approx_model`     -- (optional, default: `Forward`) approximation model;
                        see `Notes` below for possible options
- `max_order`        -- (optional, default: `5`) maximum zonotope order
- `reduction_method` -- (optional, default: `GIR05()`) zonotope order reduction method used
- `static`           -- (optional, default: `false`) if `true`, convert the problem data
                        to statically sized arrays
- `recursive`        -- (optional default: `true`) if `true`, use the implementation that
                        recursively computes each reach-set; otherwise, use the implementation
                        that unwraps the sequence until the initial set

## Notes

The type fields are:

- `N`  -- number type of the step-size
- `AM` -- type of the approximation model
- `RM` -- type associated to the reduction method
- `S`  -- value type associated to the `static` option
- `R`  -- value type associated to the `recursive` option

The sole parameter which doesn't have a default value is the step-size,
associated to the type parameter `N`.

The default approximation model is

```julia
approx_model=CorrectionHull(order=10, exp=:base)
```
Here, `CorrectionHull` refers to an implementation of the interval matrix
approximation method described in [[ASB07]](@ref). For technicalities on
interval matrix operations, we refer to the package `IntervalMatrices.jl`.

## References

The main ideas behind this algorithm can be found in [[ASB07]](@ref).
These methods are discussed at length in the dissertation [[ALT10]](@ref).

Regarding the zonotope order reduction methods, we refer to [[COMB03]](@ref),
[[GIR05]](@ref) and the review article [[YS18]](@ref).
"""
struct ASB07{N, AM, RM, S, R, D, NG} <: AbstractContinuousPost
    δ::N
    approx_model::AM
    max_order::Int
    reduction_method::RM
    static::S
    recursive::R
    dim::D
    ngens::NG
end

# convenience constructor using symbols
function ASB07(; δ::N,
               approx_model::AM=CorrectionHull(order=10, exp=:interval),
               max_order::Int=5,
               reduction_method::RM=GIR05(),
               static::Bool=false,
               recursive::Bool=true,
               dim::Union{Int, Missing}=missing,
               ngens::Union{Int, Missing}=missing) where {N, AM, RM}

    n = ismissing(dim) ? missing : Val(dim)
    p = ismissing(ngens) ? missing : Val(ngens)
    return ASB07(δ, approx_model, max_order, reduction_method, Val(static), Val(recursive), n, p)
end

step_size(alg::ASB07) = alg.δ
numtype(::ASB07{N}) where {N} = N

function rsetrep(alg::ASB07{N, AM, RM, Val{false}}) where {N, AM, RM}
    RT = ReachSet{N, Zonotope{N, Vector{N}, Matrix{N}}}
end

function rsetrep(alg::ASB07{N, AM, RM, S, R, Val{n}, Val{p}}) where {N, AM, RM, S, R, n, p}
    VT = SVector{n, N}
    MT = SMatrix{n, p, N, n*p}
    ZT = Zonotope{N, VT, MT}
    RT = ReachSet{N, ZT}
end

include("post.jl")
include("reach_homog.jl")
include("reach_inhomog.jl")