text:
using SparseArrays
using LinearAlgebra
using UnicodePlots

include("diagonal_sbp.jl")

# flatten tuples to arrays
flatten_tuples(x) = reshape(collect(Iterators.flatten(x)), length(x[1]),
                            length(x))

⊗(A,B) = kron(A, B)

const BC_DIRICHLET        = 1
const BC_NEUMANN          = 2
const BC_LOCKED_INTERFACE = 0
const BC_JUMP_INTERFACE   = 7

#{{{ Transfinite Blend
function transfinite_blend(α1, α2, α3, α4, α1s, α2s, α3r, α4r, r, s)
  # +---4---+
  # |       |
  # 1       2
  # |       |
  # +---3---+
  @assert [α1(-1) α2(-1) α1( 1) α2( 1)] ≈ [α3(-1) α3( 1) α4(-1) α4( 1)]

  x = (1 .+ r) .* α2(s)/2 + (1 .- r) .* α1(s)/2 +
      (1 .+ s) .* α4(r)/2 + (1 .- s) .* α3(r)/2 -
     ((1 .+ r) .* (1 .+ s) .* α2( 1) +
      (1 .- r) .* (1 .+ s) .* α1( 1) +
      (1 .+ r) .* (1 .- s) .* α2(-1) +
      (1 .- r) .* (1 .- s) .* α1(-1)) / 4

  xr =  α2(s)/2 - α1(s)/2 +
        (1 .+ s) .* α4r(r)/2 + (1 .- s) .* α3r(r)/2 -
      (+(1 .+ s) .* α2( 1) +
       -(1 .+ s) .* α1( 1) +
       +(1 .- s) .* α2(-1) +
       -(1 .- s) .* α1(-1)) / 4

  xs = (1 .+ r) .* α2s(s)/2 + (1 .- r) .* α1s(s)/2 +
       α4(r)/2 - α3(r)/2 -
      (+(1 .+ r) .* α2( 1) +
       +(1 .- r) .* α1( 1) +
       -(1 .+ r) .* α2(-1) +
       -(1 .- r) .* α1(-1)) / 4

  return (x, xr, xs)
end

function transfinite_blend(α1, α2, α3, α4, r, s, p)
  (Nrp, Nsp) = size(r)
  (Dr, _, _, _) = diagonal_sbp_D1(p, Nrp-1; xc = (-1,1))
  (Ds, _, _, _) = diagonal_sbp_D1(p, Nsp-1; xc = (-1,1))

  α2s(s) = α2(s) * Ds'
  α1s(s) = α1(s) * Ds'
  α4r(s) = Dr * α4(r)
  α3r(s) = Dr * α3(r)

  transfinite_blend(α1, α2, α3, α4, α1s, α2s, α3r, α4r, r, s)
end

function transfinite_blend(v1::T1, v2::T2, v3::T3, v4::T4, r, s
                          ) where {T1 <: Number, T2 <: Number,
                                   T3 <: Number, T4 <: Number}
  e1(α) = v1 * (1 .- α) / 2 + v3 * (1 .+ α) / 2
  e2(α) = v2 * (1 .- α) / 2 + v4 * (1 .+ α) / 2
  e3(α) = v1 * (1 .- α) / 2 + v2 * (1 .+ α) / 2
  e4(α) = v3 * (1 .- α) / 2 + v4 * (1 .+ α) / 2
  e1α(α) = -v1 / 2 + v3 / 2
  e2α(α) = -v2 / 2 + v4 / 2
  e3α(α) = -v1 / 2 + v2 / 2
  e4α(α) = -v3 / 2 + v4 / 2
  transfinite_blend(e1, e2, e3, e4, e1α, e2α, e3α, e4α, r, s)
end
#}}}

#{{{ connectivityarrays
function connectivityarrays(EToV, EToF)
  # number of elements
  nelems = size(EToV, 2)
  nfaces = maximum(maximum(EToF))

  # Determine secondary arrays
  # FToE : Unique Global Face to Element Number
  # FToLF: Unique Global Face to Element local face number
  # EToO : Element to Unique Global Faces Orientation
  # EToS : Element to Unique Global Face Side

  FToE  = zeros(Int64, 2, nfaces)
  FToLF = zeros(Int64, 2, nfaces)
  EToO  = Array{Bool,2}(undef, 4, nelems)
  EToS  = zeros(Int64, 4, nelems)

  # Local Face to Local Vertex map
  LFToLV = flatten_tuples(((1,3), (2, 4), (1,2), (3,4)))
  for e = 1:nelems
    for lf = 1:4
      gf = EToF[lf, e]
      if FToE[1, gf] == 0
        @assert FToLF[1, gf] == 0
        FToE[1, gf] = e
        FToLF[1, gf] = lf
        EToO[lf, e] = true
        EToS[lf, e] = 1
      else
        @assert FToE[2, gf] == 0
        @assert FToLF[2, gf] == 0
        FToE[2, gf] = e
        FToLF[2, gf] = lf
        EToS[lf, e] = 2

        ne = FToE[1, gf]
        nf = FToLF[1, gf]

        nv = EToV[LFToLV[:,nf], ne]
        lv = EToV[LFToLV[:,lf], e]
        if nv == lv
          EToO[lf, e] = true
        elseif nv[end:-1:1] == lv
          EToO[lf, e] = false
        else
          error("problem with connectivity")
        end
      end
    end
  end
  (FToE, FToLF, EToO, EToS)
end
#}}}

#{{{ locoperator
function create_metrics(pm, Nr, Ns,
                        xf=(r,s)->(r, ones(size(r)), zeros(size(r))),
                        yf=(r,s)->(s, zeros(size(s)), ones(size(s))))
  Nrp = Nr + 1
  Nsp = Ns + 1
  Np = Nrp * Nsp

  # Derivative operators for the metric terms
  @assert pm <= 8
  pp = pm == 6 ? 8 : pm

  r = range(-1, stop=1, length=Nrp)
  s = range(-1, stop=1, length=Nsp)

  # Create the mesh
  r = ones(1, Nsp) ⊗ r
  s = s' ⊗ ones(Nrp)
  (x, xr, xs) = xf(r, s)
  (y, yr, ys) = yf(r, s)

  J = xr .* ys - xs .* yr
  @assert minimum(J) > 0

  rx =  ys ./ J
  sx = -yr ./ J
  ry = -xs ./ J
  sy =  xr ./ J

  # variable coefficient matrix components
  crr = J .* (rx .* rx + ry .* ry)
  crs = J .* (sx .* rx + sy .* ry)
  css = J .* (sx .* sx + sy .* sy)

  #
  # Block surface matrices
  #
  (xf1, yf1) = (view(x, 1, :), view(y, 1, :))
  nx1 = -ys[1, :]
  ny1 =  xs[1, :]
  sJ1 = hypot.(nx1, ny1)
  nx1 = nx1 ./ sJ1
  ny1 = ny1 ./ sJ1

  (xf2, yf2) = (view(x, Nrp, :), view(y, Nrp, :))
  nx2 =  ys[end, :]
  ny2 = -xs[end, :]
  sJ2 = hypot.(nx2, ny2)
  nx2 = nx2 ./ sJ2
  ny2 = ny2 ./ sJ2

  (xf3, yf3) = (view(x, :, 1), view(y, :, 1))
  nx3 =  yr[:, 1]
  ny3 = -xr[:, 1]
  sJ3 = hypot.(nx3, ny3)
  nx3 = nx3 ./ sJ3
  ny3 = ny3 ./ sJ3

  (xf4, yf4) = (view(x, :, Nsp), view(y, :, Nsp))
  nx4 = -yr[:, end]
  ny4 =  xr[:, end]
  sJ4 = hypot.(nx4, ny4)
  nx4 = nx4 ./ sJ4
  ny4 = ny4 ./ sJ4

  (coord = (x,y),
   facecoord = ((xf1, xf2, xf3, xf4), (yf1, yf2, yf3, yf4)),
   crr = crr, css = css, crs = crs,
   J=J,
   sJ = (sJ1, sJ2, sJ3, sJ4),
   nx = (nx1, nx2, nx3, nx4),
   ny = (ny1, ny2, ny3, ny4),
   rx = rx, ry = ry, sx = sx, sy = sy)
end

function locoperator(p, Nr, Ns, metrics=create_metrics(p,Nr,Ns),
                     LFToB = (BC_DIRICHLET, BC_DIRICHLET,
                              BC_DIRICHLET, BC_DIRICHLET);
                     τscale = 2,
                     crr = metrics.crr,
                     css = metrics.css,
                     crs = metrics.crs)
  csr = crs
  J = metrics.J

  Nrp = Nr + 1
  Nsp = Ns + 1
  Np = Nrp * Nsp

  # Derivative operators for the rest of the computation
  (Dr, HrI, Hr, r) = diagonal_sbp_D1(p, Nr; xc = (-1,1))
  Qr = Hr * Dr
  QrT = sparse(transpose(Qr))

  (Ds, HsI, Hs, s) = diagonal_sbp_D1(p, Ns; xc = (-1,1))
  Qs = Hs * Ds
  QsT = sparse(transpose(Qs))

  # Identity matrices for the comuptation
  Ir = sparse(I, Nrp, Nrp)
  Is = sparse(I, Nsp, Nsp)

  #{{{ Set up the rr derivative matrix
  ISr0 = Array{Int64,1}(undef,0)
  JSr0 = Array{Int64,1}(undef,0)
  VSr0 = Array{Float64,1}(undef,0)
  ISrN = Array{Int64,1}(undef,0)
  JSrN = Array{Int64,1}(undef,0)
  VSrN = Array{Float64,1}(undef,0)

  (_, S0e, SNe, _, _, Ae, _) = variable_diagonal_sbp_D2(p, Nr, rand(Nrp))
  IArr = Array{Int64,1}(undef,Nsp * length(Ae.nzval))
  JArr = Array{Int64,1}(undef,Nsp * length(Ae.nzval))
  VArr = Array{Float64,1}(undef,Nsp * length(Ae.nzval))
  stArr = 0

  ISr0 = Array{Int64,1}(undef,Nsp * length(S0e.nzval))
  JSr0 = Array{Int64,1}(undef,Nsp * length(S0e.nzval))
  VSr0 = Array{Float64,1}(undef,Nsp * length(S0e.nzval))
  stSr0 = 0

  ISrN = Array{Int64,1}(undef,Nsp * length(SNe.nzval))
  JSrN = Array{Int64,1}(undef,Nsp * length(SNe.nzval))
  VSrN = Array{Float64,1}(undef,Nsp * length(SNe.nzval))
  stSrN = 0
  for j = 1:Nsp
    rng = (j-1) * Nrp .+ (1:Nrp)
    (_, S0e, SNe, _, _, Ae, _) = variable_diagonal_sbp_D2(p, Nr, crr[rng])
    (Ie, Je, Ve) = findnz(Ae)
    IArr[stArr .+ (1:length(Ve))] = Ie .+ (j-1) * Nrp
    JArr[stArr .+ (1:length(Ve))] = Je .+ (j-1) * Nrp
    VArr[stArr .+ (1:length(Ve))] = Hs[j,j] * Ve
    stArr += length(Ve)

    (Ie, Je, Ve) = findnz(S0e)
    ISr0[stSr0 .+ (1:length(Ve))] = Ie .+ (j-1) * Nrp
    JSr0[stSr0 .+ (1:length(Ve))] = Je .+ (j-1) * Nrp
    VSr0[stSr0 .+ (1:length(Ve))] =  Hs[j,j] * Ve
    stSr0 += length(Ve)

    (Ie, Je, Ve) = findnz(SNe)
    ISrN[stSrN .+ (1:length(Ve))] = Ie .+ (j-1) * Nrp
    JSrN[stSrN .+ (1:length(Ve))] = Je .+ (j-1) * Nrp
    VSrN[stSrN .+ (1:length(Ve))] =  Hs[j,j] * Ve
    stSrN += length(Ve)
  end
  Ãrr = sparse(IArr[1:stArr], JArr[1:stArr], VArr[1:stArr], Np, Np)
  Sr0 = sparse(ISr0[1:stSr0], JSr0[1:stSr0], VSr0[1:stSr0], Np, Np)
  SrN = sparse(ISrN[1:stSrN], JSrN[1:stSrN], VSrN[1:stSrN], Np, Np)
  Sr0T = sparse(JSr0[1:stSr0], ISr0[1:stSr0], VSr0[1:stSr0], Np, Np)
  SrNT = sparse(JSrN[1:stSrN], ISrN[1:stSrN], VSrN[1:stSrN], Np, Np)
  #= affine mesh test
  # @assert Ãrr ≈ Ãrr'
  (D2, S0, SN, _, _, _) = diagonal_sbp_D2(p, Nr)
  Ar = SN - S0 - Hr * D2
  @assert Ãrr ≈ Hs ⊗ Ar
  =#
  # @assert Sr0 ≈ ((sparse(Diagonal(crr[1   .+ Nrp*(0:Ns)])) * Hs) ⊗ S0)
  # @assert SrN ≈ ((sparse(Diagonal(crr[Nrp .+ Nrp*(0:Ns)])) * Hs) ⊗ SN)
  #}}}

  #{{{ Set up the ss derivative matrix
  (_, S0e, SNe, _, _, Ae, _) = variable_diagonal_sbp_D2(p, Ns, rand(Nsp))
  IAss = Array{Int64,1}(undef,Nrp * length(Ae.nzval))
  JAss = Array{Int64,1}(undef,Nrp * length(Ae.nzval))
  VAss = Array{Float64,1}(undef,Nrp * length(Ae.nzval))
  stAss = 0

  ISs0 = Array{Int64,1}(undef,Nrp * length(S0e.nzval))
  JSs0 = Array{Int64,1}(undef,Nrp * length(S0e.nzval))
  VSs0 = Array{Float64,1}(undef,Nrp * length(S0e.nzval))
  stSs0 = 0

  ISsN = Array{Int64,1}(undef,Nrp * length(SNe.nzval))
  JSsN = Array{Int64,1}(undef,Nrp * length(SNe.nzval))
  VSsN = Array{Float64,1}(undef,Nrp * length(SNe.nzval))
  stSsN = 0
  for i = 1:Nrp
    rng = i .+ Nrp * (0:Ns)
    (_, S0e, SNe, _, _, Ae, _) = variable_diagonal_sbp_D2(p, Ns, css[rng])
    R = Ae - Dr' * Hr * Diagonal(css[rng]) * Dr

    (Ie, Je, Ve) = findnz(Ae)
    IAss[stAss .+ (1:length(Ve))] = i .+ Nrp * (Ie .- 1)
    JAss[stAss .+ (1:length(Ve))] = i .+ Nrp * (Je .- 1)
    VAss[stAss .+ (1:length(Ve))] = Hr[i,i] * Ve
    stAss += length(Ve)

    (Ie, Je, Ve) = findnz(S0e)
    ISs0[stSs0 .+ (1:length(Ve))] = i .+ Nrp * (Ie .- 1)
    JSs0[stSs0 .+ (1:length(Ve))] = i .+ Nrp * (Je .- 1)
    VSs0[stSs0 .+ (1:length(Ve))] = Hr[i,i] * Ve
    stSs0 += length(Ve)

    (Ie, Je, Ve) = findnz(SNe)
    ISsN[stSsN .+ (1:length(Ve))] = i .+ Nrp * (Ie .- 1)
    JSsN[stSsN .+ (1:length(Ve))] = i .+ Nrp * (Je .- 1)
    VSsN[stSsN .+ (1:length(Ve))] = Hr[i,i] * Ve
    stSsN += length(Ve)
  end
  Ãss = sparse(IAss[1:stAss], JAss[1:stAss], VAss[1:stAss], Np, Np)
  Ss0 = sparse(ISs0[1:stSs0], JSs0[1:stSs0], VSs0[1:stSs0], Np, Np)
  SsN = sparse(ISsN[1:stSsN], JSsN[1:stSsN], VSsN[1:stSsN], Np, Np)
  Ss0T = sparse(JSs0[1:stSs0], ISs0[1:stSs0], VSs0[1:stSs0], Np, Np)
  SsNT = sparse(JSsN[1:stSsN], ISsN[1:stSsN], VSsN[1:stSsN], Np, Np)
  # @assert Ãss ≈ Ãss'
  #= affine mesh test
  (D2, S0, SN, _, _, _) = diagonal_sbp_D2(p, Ns)
  As = SN - S0 - Hs * D2
  @assert Ãss ≈ As ⊗ Hr
  =#
  # @assert Ss0 ≈ (S0 ⊗ (Hr * sparse(Diagonal(css[1:Nrp]))))
  # @assert SsN ≈ (SN ⊗ (Hr * sparse(Diagonal(css[Nrp*Ns .+ (1:Nrp)]))))
  #}}}

  #{{{ Set up the sr and rs derivative matrices
  Ãsr = (QsT ⊗ Ir) * sparse(1:length(crs), 1:length(crs), view(crs, :)) * (Is ⊗ Qr)
  Ãrs = (Is ⊗ QrT) * sparse(1:length(csr), 1:length(csr), view(csr, :)) * (Qs ⊗ Ir)
  #}}}

  Ã = Ãrr + Ãss + Ãrs + Ãsr

  #
  # Boundary point matrices
  #
  Er0 = sparse([1], [1], [1], Nrp, Nrp)
  ErN = sparse([Nrp], [Nrp], [1], Nrp, Nrp)
  Es0 = sparse([1], [1], [1], Nsp, Nsp)
  EsN = sparse([Nsp], [Nsp], [1], Nsp, Nsp)

  er0 = sparse([1  ], [1], [1], Nrp, 1)
  erN = sparse([Nrp], [1], [1], Nrp, 1)
  es0 = sparse([1  ], [1], [1], Nsp, 1)
  esN = sparse([Nsp], [1], [1], Nsp, 1)

  er0T = sparse([1], [1  ], [1], 1, Nrp)
  erNT = sparse([1], [Nrp], [1], 1, Nrp)
  es0T = sparse([1], [1  ], [1], 1, Nsp)
  esNT = sparse([1], [Nsp], [1], 1, Nsp)

  #
  # Store coefficient matrices as matrices
  #
  crs0 = sparse(Diagonal(crs[1:Nrp]))
  crsN = sparse(Diagonal(crs[Nrp*Ns .+ (1:Nrp)]))
  csr0 = sparse(Diagonal(csr[1   .+ Nrp*(0:Ns)]))
  csrN = sparse(Diagonal(csr[Nrp .+ Nrp*(0:Ns)]))

  #
  # Surface mass matrices
  #
  H1 = Hs
  H1I = HsI

  H2 = Hs
  H2I = HsI

  H3 = Hr
  H3I = HrI

  H4 = Hr
  H4I = HrI

  #
  # Penalty terms
  #
  if p == 2
    l = 2
    β = 0.363636363
    α = 1 / 2
  elseif p == 4
    l = 4
    β = 0.2505765857
    α = 17 / 48
  elseif p == 6
    l = 7
    β = 0.1878687080
    α = 13649 / 43200
  else
    error("unknown order")
  end

  ψmin = reshape((crr + css - sqrt.((crr - css).^2 + 4crs.^2)) / 2, Nrp, Nsp)
  @assert minimum(ψmin) > 0

  hr = 2 / Nr
  hs = 2 / Ns

  ψ1 = ψmin[  1, :]
  ψ2 = ψmin[Nrp, :]
  ψ3 = ψmin[:,   1]
  ψ4 = ψmin[:, Nsp]
  for k = 2:l
    ψ1 = min.(ψ1, ψmin[k, :])
    ψ2 = min.(ψ2, ψmin[Nrp+1-k, :])
    ψ3 = min.(ψ3, ψmin[:, k])
    ψ4 = min.(ψ4, ψmin[:, Nsp+1-k])
  end
  τ1 = (2τscale / hr) * (crr[  1, :].^2 / β + crs[  1, :].^2 / α) ./ ψ1
  τ2 = (2τscale / hr) * (crr[Nrp, :].^2 / β + crs[Nrp, :].^2 / α) ./ ψ2
  τ3 = (2τscale / hs) * (css[:,   1].^2 / β + crs[:,   1].^2 / α) ./ ψ3
  τ4 = (2τscale / hs) * (css[:, Nsp].^2 / β + crs[:, Nsp].^2 / α) ./ ψ4

  τ1 = sparse(1:Nsp, 1:Nsp, τ1)
  τ2 = sparse(1:Nsp, 1:Nsp, τ2)
  τ3 = sparse(1:Nrp, 1:Nrp, τ3)
  τ4 = sparse(1:Nrp, 1:Nrp, τ4)

  C̃1 =  (Sr0 + Sr0T) + ((csr0 * Qs + QsT * csr0) ⊗ Er0) + ((τ1 * H1) ⊗ Er0)
  C̃2 = -(SrN + SrNT) - ((csrN * Qs + QsT * csrN) ⊗ ErN) + ((τ2 * H2) ⊗ ErN)
  C̃3 =  (Ss0 + Ss0T) + (Es0 ⊗ (crs0 * Qr + QrT * crs0)) + (Es0 ⊗ (τ3 * H3))
  C̃4 = -(SsN + SsNT) - (EsN ⊗ (crsN * Qr + QrT * crsN)) + (EsN ⊗ (τ4 * H4))

  # TODO: Fix minus sign (reverse of the paper)
  G1 = -(Is ⊗ er0T) * Sr0 - ((csr0 * Qs) ⊗ er0T)
  G2 = +(Is ⊗ erNT) * SrN + ((csrN * Qs) ⊗ erNT)
  G3 = -(es0T ⊗ Ir) * Ss0 - (es0T ⊗ (crs0 * Qr))
  G4 = +(esNT ⊗ Ir) * SsN + (esNT ⊗ (crsN * Qr))

  F1 = G1' - ((τ1 * H1) ⊗ er0)
  F2 = G2' - ((τ2 * H2) ⊗ erN)
  F3 = G3' - (es0 ⊗ (τ3 * H3))
  F4 = G4' - (esN ⊗ (τ4 * H4))

  HfI_F1T = H1I * G1 - (τ1 ⊗ er0')
  HfI_F2T = H2I * G2 - (τ2 ⊗ erN')
  HfI_F3T = H3I * G3 - (es0' ⊗ τ3)
  HfI_F4T = H4I * G4 - (esN' ⊗ τ4)

  HfI_G1 = H1I * G1
  HfI_G2 = H2I * G2
  HfI_G3 = H3I * G3
  HfI_G4 = H4I * G4

  M̃ = Ã + C̃1 + C̃2 + C̃3 + C̃4

  # Modify the operator to handle the boundary conditions
  bctype=(BC_LOCKED_INTERFACE, BC_LOCKED_INTERFACE, BC_LOCKED_INTERFACE, BC_LOCKED_INTERFACE)
  F = (F1, F2, F3, F4)
  τ = (τ1, τ2, τ3, τ4)
  HfI = (H1I, H2I, H3I, H4I)
  # Modify operators for the BC
  for lf = 1:4
    if LFToB[lf] == BC_NEUMANN
      M̃ -= F[lf] * (Diagonal(1 ./ (diag(τ[lf]))) * HfI[lf]) * F[lf]'
    elseif !(LFToB[lf] == BC_DIRICHLET ||
             LFToB[lf] == BC_LOCKED_INTERFACE ||
             LFToB[lf] >= BC_JUMP_INTERFACE)
      error("invalid bc")
    end
  end
  bctype=(LFToB[1], LFToB[2], LFToB[3], LFToB[4])

  # (E, V) = eigen(Matrix(M̃))
  # println((minimum(E), maximum(E)))
  JH = sparse(1:Np, 1:Np, view(J, :)) * (Hs ⊗ Hr)
  (M̃ = M̃,
   F = (F1, F2, F3, F4),
   HfI_FT = (HfI_F1T, HfI_F2T, HfI_F3T, HfI_F4T),
   HfI_G = (HfI_G1, HfI_G2, HfI_G3, HfI_G4),
   coord = metrics.coord,
   facecoord = metrics.facecoord,
   JH = JH,
   sJ = metrics.sJ,
   nx = metrics.nx,
   ny = metrics.ny,
   Hf = (H1, H2, H3, H4),
   HfI = (H1I, H2I, H3I, H4I),
   τ = (τ1, τ2, τ3, τ4),
   bctype=bctype)
end
#}}}

#{{{ gloλoperator: Build the trace operators
function gloλoperator(lop, vstarts, FToB, FToE, FToLF, EToO, EToS, Nr, Ns)
  nelems = length(lop)
  nfaces = length(FToB)
  Nλp = zeros(Int64, nfaces)
  FToλstarts = Array{Int64, 1}(undef, nfaces + 1)
  FToλstarts[1] = 1
  IT = Array{Int64,1}(undef,0)
  JT = Array{Int64,1}(undef,0)
  VT = Array{Float64,1}(undef,0)
  VD = Array{Float64,1}(undef,0)
  for f = 1:nfaces
    if FToB[f] == BC_DIRICHLET || FToB[f] == BC_NEUMANN
      FToλstarts[f+1] = FToλstarts[f]
      continue
    end
    (em, ep) = FToE[:, f]
    (fm, fp) = FToLF[:, f]
    Nλp[f] = (fm <= 2 ? Ns[em]+1 : Nr[em]+1)
    @assert Nλp[f] == (fp <= 2 ? Ns[ep]+1 : Nr[ep]+1)
    FToλstarts[f+1] = FToλstarts[f] + Nλp[f]

    @assert EToO[fm, em] && EToS[fm, em] == 1
    Fm = lop[em].F[fm]
    # swap I and J to get transpose
    (Je, Ie, Ve) = findnz(Fm)
    IT = [IT; Ie .+ (FToλstarts[f] - 1)]
    JT = [JT; Je .+ (vstarts[em] - 1)]
    VT = [VT; Ve]

    @assert EToS[fp, ep] == 2
    Fp = lop[ep].F[fp]
    # swap I and J to get transpose
    (Je, Ie, Ve) = findnz(Fp)
    # if element and face orientation do not match, then flip
    if EToO[fp, ep]
      IT = [IT; Ie .+ (FToλstarts[f] - 1)]
      τm = Vector(diag(lop[em].τ[fm]))
      τp = Vector(diag(lop[ep].τ[fp]))
    else
      IT = [IT; FToλstarts[f+1] .- Ie]
      τm = Vector(diag(lop[em].τ[fm]))
      τp = Vector(diag(rot180(lop[ep].τ[fp])))
    end
    JT = [JT; Je .+ (vstarts[ep] - 1)]
    VT = [VT; Ve]

    Hf = Vector(diag(lop[em].Hf[fm]))
    VD = [VD; Hf .* (τm + τp)]

  end
  λNp = FToλstarts[nfaces+1]-1
  VNp = vstarts[nelems+1]-1
  FbarT = sparse(IT, JT, VT, λNp, VNp)
  # Ttranspose = sparse(JT, IT, VT, VNp, λNp)
  (FToλstarts, FbarT, VD)
end
#}}}

#{{{ volbcarray()
function locbcarray_mod!(ge, lop, LFToB, bc_Dirichlet, bc_Neumann,
                     bcargs = ())
  F = lop.F
  (xf, yf) = lop.facecoord
  Hf = lop.Hf
  sJ = lop.sJ
  nx = lop.nx
  ny = lop.ny
  τ = lop.τ
  ge[:] .= 0
  for lf = 1:4
    if LFToB[lf] == BC_DIRICHLET
      vf = bc_Dirichlet(lf, xf[lf], yf[lf], bcargs...)
    elseif LFToB[lf] == BC_NEUMANN
      gN = bc_Neumann(lf, xf[lf], yf[lf], nx[lf], ny[lf], bcargs...)
      vf = sJ[lf] .* gN ./ diag(τ[lf])
    elseif LFToB[lf] == BC_LOCKED_INTERFACE
      continue # nothing to do here
    else
      error("invalid bc")
    end
    ge[:] -= F[lf] * vf
  end
end

#{{{ volbcarray()
function locbcarray!(ge, gδe, lop, LFToB, bc_Dirichlet, bc_Neumann, in_jump,
                     bcargs = ())
  F = lop.F
  (xf, yf) = lop.facecoord
  Hf = lop.Hf
  sJ = lop.sJ
  nx = lop.nx
  ny = lop.ny
  τ = lop.τ
  ge[:] .= 0
  for lf = 1:4
    if LFToB[lf] == BC_DIRICHLET
      vf = bc_Dirichlet(lf, xf[lf], yf[lf], bcargs...)
    elseif LFToB[lf] == BC_NEUMANN
      gN = bc_Neumann(lf, xf[lf], yf[lf], nx[lf], ny[lf], bcargs...)
      vf = sJ[lf] .* gN ./ diag(τ[lf])
    elseif LFToB[lf] == BC_LOCKED_INTERFACE
      continue # nothing to do here
    elseif LFToB[lf] >= BC_JUMP_INTERFACE
      # In this case we need to add in half the jump
      vf = in_jump(lf, xf[lf], yf[lf], bcargs...) / 2
      gδe[lf][:] -= Hf[lf] * τ[lf] * vf
    else
      error("invalid bc")
    end
    ge[:] -= F[lf] * vf
  end
end
#}}}

#{{{ computetraction
function computetraction_mod(lop, lf, u, δ)
  HfI_FT = lop.HfI_FT[lf]
  τf = lop.τ[lf]
  sJ = lop.sJ[lf]

  return (HfI_FT * u + τf * (δ .- δ / 2)) ./ sJ
end

#{{{ computetraction
function computetraction(lop, lf, u, λ, δ)
  HfI_FT = lop.HfI_FT[lf]
  τf = lop.τ[lf]
  sJ = lop.sJ[lf]

  return (HfI_FT * u + τf * (λ .- δ / 2)) ./ sJ
end
#}}}

#{{{ volsourcearray()
function locsourcearray!(ge, source, lop, volargs = ())

  (xloc, yloc) = lop.coord
  JHloc = lop.JH
  ge[:] += JHloc * source(xloc[:], yloc[:], volargs...)

end

#}}}

#{{{
struct SBPLocalOperator1{T<:Real, S<:Factorization}
  offset::Array{Int64,1}
  H::Array{T,1}
  X::Array{T,1}
  Y::Array{T,1}
  E::Array{Int64,1}
  F::Array{S,1}
  SBPLocalOperator1{T,S}(vstarts::Array{Int64,1}, H::Array{T,1}, X::Array{T,1},
                         Y::Array{T,1}, E::Array{Int64,1},
                         F::Array{S,1}) where {T<:Real, S<:Factorization} =
  new(vstarts, H, X, Y, E, F)
end

function SBPLocalOperator1(lop, Nr, Ns, factorization)
  nelems = length(lop)
  vstarts = Array{Int64, 1}(undef, nelems + 1)
  vstarts[1] = 1
  Np = Array{Int64, 1}(undef, nelems)
  VH = Array{Float64,1}(undef,0)
  X = Array{Float64,1}(undef,0)
  Y = Array{Float64,1}(undef,0)
  E = Array{Int64,1}(undef,0)
  FTYPE = typeof(factorization(sparse([1],[1],[1.0])))
  factors = Array{FTYPE, 1}(undef, nelems)
  for e = 1:nelems
    # Fill arrays to build global sparse matrix
    Np[e] = (Nr[e]+1)*(Ns[e]+1)
    vstarts[e+1] = vstarts[e] + Np[e]

    # Global "mass" matrix
    JH = lop[e].JH
    VH = [VH;Vector(diag(JH))]

    # global coordinates and element number array (needed for jump)
    (x,y) = lop[e].coord
    X = [X;x[:]]
    Y = [Y;y[:]]
    E = [E;e * ones(Int64, Np[e])]

    factors[e] = factorization(lop[e].M̃)
  end
  VNp = vstarts[nelems+1]-1 # total number of volume points

  SBPLocalOperator1{Float64, FTYPE}(vstarts, VH, X, Y, E, factors)
end
#}}}

function LocalGlobalOperators(lop, Nr, Ns, FToB, FToE, FToLF, EToO, EToS,
                              factorization)
  M = SBPLocalOperator1(lop, Nr, Ns, factorization)
  (FToλstarts, FbarT, D) = gloλoperator(lop, M.offset, FToB, FToE, FToLF, EToO,
                                        EToS, Nr, Ns)
  (M, FbarT, D, M.offset, FToλstarts)
end

function bcstarts(FToB, FToE, FToLF, bctype, Nr, Ns)
  nfaces = length(FToB)
  bcstarts = Array{Int64, 1}(undef, nfaces + 1)
  bcstarts[1] = 1
  for f = 1:nfaces
    if FToB[f] ∈ bctype
      e  = FToE[1,f]
      lf = FToLF[1,f]
      bcstarts[f+1] = bcstarts[f] + (lf ∈ (1,2) ? Ns[e] : Nr[e]) + 1
    else
      bcstarts[f+1] = bcstarts[f]
    end
  end
  bcstarts
end

function LocalToGLobalRHS!(b, g, gδ, u, M, FbarT, vstarts)
  u .= 0
  @inbounds for e = 1:length(M)
    if maximum(abs.(g[vstarts[e]:(vstarts[e+1]-1)])) > 0
      @views u[vstarts[e]:(vstarts[e+1]-1)] = (M[e] \
                                               g[vstarts[e]:(vstarts[e+1]-1)])
    end
  end
  mul!(b, FbarT, u)
  @. b = gδ - b
end

#{{{ assembleλmatrix: Schur complement system
function assembleλmatrix(FToλstarts, vstarts, EToF, FToB, F, D, FbarT)
  nfaces = length(FToλstarts)-1
  nelems = length(vstarts)-1
  λNp = FToλstarts[nfaces+1]-1
  sz = λNp

  for e = 1:nelems
    lλs = Array{Int64, 1}(undef, 4)
    for lf = 1:4
      f = EToF[lf,e]
      lλs[lf] = FToλstarts[f+1] - FToλstarts[f]
    end
    for lf = 1:4
      sz += lλs[lf]*sum(lλs)
    end
  end
  Ie = Array{Int64, 1}(undef, sz)
  Je = Array{Int64, 1}(undef, sz)
  Ve = Array{Float64, 1}(undef, sz)
  Ie[1:λNp] = 1:λNp
  Je[1:λNp] = 1:λNp
  Ve[1:λNp] = D
  offset = λNp
  Fbar = FbarT'
  for e = 1:nelems
    # println((e, nelems))
    vrng = vstarts[e]:(vstarts[e+1]-1)
    for lf = 1:4
      f = EToF[lf,e]
      if FToB[f] == BC_LOCKED_INTERFACE || FToB[f] >= BC_JUMP_INTERFACE
        λrng = FToλstarts[f]:(FToλstarts[f+1]-1)
        B = -(Matrix(F[e]' \ Fbar[vrng, λrng]))
        for lf2 = 1:4
          f2 = EToF[lf2,e]
          if FToB[f2] == BC_LOCKED_INTERFACE || FToB[f2] >= BC_JUMP_INTERFACE
            λrng2 = FToλstarts[f2]:(FToλstarts[f2+1]-1)
            C = -(FbarT[λrng2, vrng] * B)
            λblck = λrng*ones(Int64, 1, length(λrng2))
            λblck2 = ones(Int64, length(λrng), 1) * λrng2'
            last = length(λrng) * length(λrng2)
            Ie[offset.+(1:last)] = λblck[:]
            Je[offset.+(1:last)] = λblck2[:]
            Ve[offset.+(1:last)] = -C'[:]
            offset += last
          end
        end
      end
    end
  end
  @assert offset == sz
  B = sparse(Ie, Je, Ve, λNp, λNp)
  @assert B ≈ B'
  # println((λNp * λNp, nnz(B), nnz(B) / λNp^2))
  B
end

#}}}

# {{{ Constructor for inp files
function read_inp_2d(T, S, filename::String; bc_map=1:10000)
  # {{{ Read in the file
  f = try
    open(filename)
  catch
    error("InpRead cannot open \"$filename\" ")
  end
  lines = readlines(f)
  close(f)
  # }}}

  # {{{ Read in nodes
  str = "NSET=ALLNODES"
  linenum = SeekToSubstring(lines, str);
  linenum > 0 || error("did not find: $str")
  num_nodes = 0
  for l = linenum+1:length(lines)
    occursin(r"^\s*[0-9]*\s*,.*", lines[l]) ? num_nodes+=1 : break
  end
  Vx = fill(S(NaN), num_nodes)
  Vy = fill(S(NaN), num_nodes)
  Vz = fill(S(NaN), num_nodes)
  for l = linenum .+ (1:num_nodes)
    node_data = split(lines[l], r"\s|,", keepempty=false)
    (node_num, node_x, node_y, node_z) = try
      (parse(T, node_data[1]),
       parse(S, node_data[2]),
       parse(S, node_data[3]),
       parse(S, node_data[4]))
    catch
      error("cannot parse line $l: \"$(lines[l])\" ")
    end

    Vx[node_num] = node_x
    Vy[node_num] = node_y
    Vz[node_num] = node_z
  end
  # }}}

  # {{{ Read in Elements
  str = "ELEMENT"
  linenum = SeekToSubstring(lines, str);
  num_elm = 0
  while linenum > 0
    for l = linenum .+ (1:length(lines))
      occursin(r"^\s*[0-9]*\s*,.*", lines[l]) ? num_elm+=1 : break
    end
    linenum = SeekToSubstring(lines, str; first=linenum+1)
  end
  num_elm > 0 || error("did not find any element")

  EToV = fill(T(0), 4, num_elm)
  EToBlock = fill(T(0), num_elm)
  linenum = SeekToSubstring(lines, str);
  while linenum > 0
    foo = split(lines[linenum], r"[^0-9]", keepempty=false)
    B = parse(T, foo[end])
    for l = linenum .+ (1:num_elm)
      elm_data = split(lines[l], r"\s|,", keepempty=false)
      # read into z-order
      (elm_num, elm_v1, elm_v2, elm_v4, elm_v3) = try
        (parse(T, elm_data[1]),
         parse(T, elm_data[2]),
        parse(T, elm_data[3]),
        parse(T, elm_data[4]),
        parse(T, elm_data[5]))
      catch
        break
      end
      EToV[:, elm_num] = [elm_v1, elm_v2, elm_v3, elm_v4]
      EToBlock[elm_num] = B
    end
    linenum = SeekToSubstring(lines, str; first=linenum+1)
  end
  # }}}

  # {{{ Determine connectivity
  EToF = fill(T(0), 4, num_elm)

  VsToF = Dict{Tuple{Int64, Int64}, Int64}()
  numfaces = 0
  for e = 1:num_elm
    for lf = 1:4
      if lf == 1
        Vs = (EToV[1, e], EToV[3, e])
      elseif lf == 2
        Vs = (EToV[2, e], EToV[4, e])
      elseif lf == 3
        Vs = (EToV[1, e], EToV[2, e])
      elseif lf == 4
        Vs = (EToV[3, e], EToV[4, e])
      end
      if Vs[1] > Vs[2]
        Vs = (Vs[2], Vs[1])
      end
      if haskey(VsToF, Vs)
        EToF[lf, e] = VsToF[Vs]
      else
        numfaces = numfaces + 1
        EToF[lf, e] = VsToF[Vs] = numfaces
      end
    end
  end
  #}}}

  # {{{ Read in side set info
  FToB = Array{T, 1}(undef, numfaces)
  fill!(FToB, BC_LOCKED_INTERFACE)
  linenum = SeekToSubstring(lines, "\\*ELSET")
  inp_to_zorder = [3,  2, 4, 1]
  while linenum > 0
    foo = split(lines[linenum], r"[^0-9]", keepempty=false)
    (bc, face) = try
      (parse(T, foo[1]),
       parse(T, foo[2]))
    catch
      error("cannot parse line $linenum: \"$(lines[linenum])\" ")
    end
    bc = bc_map[bc]
    face = inp_to_zorder[face]
    for l = linenum+1:length(lines)
      if !occursin(r"^\s*[0-9]+", lines[l])
        break
      end
      elms = split(lines[l], r"\s|,", keepempty=false)
      for elm in elms
        elm = try
          parse(T, elm)
        catch
          error("cannot parse line $linenum: \"$(lines[l])\" ")
        end
        if bc == 3
          bc = BC_LOCKED_INTERFACE
        end
        FToB[EToF[face, elm]] = bc
        @assert (bc == BC_DIRICHLET || bc == BC_NEUMANN ||
                 bc == BC_LOCKED_INTERFACE || bc >= BC_JUMP_INTERFACE)
      end
    end
    linenum = SeekToSubstring(lines, "\\*ELSET"; first=linenum+1)
  end
  # }}}

  ([Vx Vy]', EToV, EToF, FToB, EToBlock)
end
read_inp_2d(filename;kw...) = read_inp_2d(Int64, Float64, filename;kw...)

function SeekToSubstring(lines, substring; first=1)
  for l = first:length(lines)
    if occursin(Regex(".*$(substring).*"), lines[l])
      return l
    end
  end
  return -1
end

# }}}

function plot_connectivity(verts, EToV)
  Lx = extrema(verts[1,:])
  Lx = (floor(Int, Lx[1]), ceil(Int, Lx[2]))
  Ly = extrema(verts[2,:])
  Ly = (floor(Int, Ly[1]), ceil(Int, Ly[2]))
  width = Lx[2] - Lx[1]
  height = Ly[2] - Ly[1]
  plt = Plot(BrailleCanvas(80, ceil(Int, 40 * height / width),
                           origin_x = Lx[1], origin_y = Ly[1],
                           width = width, height = height))

  annotate!(plt, :l, nrows(plt.graphics), string(Ly[1]), color = :light_black)
  annotate!(plt, :l, 1, string(Ly[2]), color = :light_black)
  annotate!(plt, :bl, string(Lx[1]), color = :light_black)
  annotate!(plt, :br, string(Lx[2]), color = :light_black)
  for e = 1:size(EToV, 2)
    (v1, v2, v3, v4) = EToV[1:4, e]
    x = verts[1, [v1 v2 v4 v3 v1]][:]
    y = verts[2, [v1 v2 v4 v3 v1]][:]
    lineplot!(plt, x, y)
  end
  title!(plt, "connectivity")
  display(plt)
end

function plot_blocks(lop)
  Lx = (floatmax(), -floatmax())
  Ly = (floatmax(), -floatmax())
  for e = 1:length(lop)
    (x, y) = lop[e].coord
    Lxe = extrema(x)
    Lye = extrema(y)
    Lx = (min(Lx[1], Lxe[1]), max(Lx[2], Lxe[2]))
    Ly = (min(Ly[1], Lye[1]), max(Ly[2], Lye[2]))
  end

  Lx = (floor(Int, Lx[1]), ceil(Int, Lx[2]))
  Ly = (floor(Int, Ly[1]), ceil(Int, Ly[2]))

  width = Lx[2] - Lx[1]
  height = Ly[2] - Ly[1]
  plt = Plot(BrailleCanvas(80, ceil(Int, 40 * height / width),
                           origin_x = Lx[1], origin_y = Ly[1],
                           width = width, height = height))

  annotate!(plt, :l, nrows(plt.graphics), string(Ly[1]), color = :light_black)
  annotate!(plt, :l, 1, string(Ly[2]), color = :light_black)
  annotate!(plt, :bl, string(Lx[1]), color = :light_black)
  annotate!(plt, :br, string(Lx[2]), color = :light_black)

  for e = 1:length(lop)
    (xf, yf) = lop[e].facecoord
    bctype = lop[e].bctype
    for lf = 1:length(xf)
      if bctype[lf] == BC_LOCKED_INTERFACE
        lineplot!(plt, xf[lf], yf[lf], color=:blue)
      elseif bctype[lf] == BC_DIRICHLET
        lineplot!(plt, xf[lf], yf[lf], color=:green)
      elseif bctype[lf] == BC_DIRICHLET
        lineplot!(plt, xf[lf], yf[lf], color=:yellow)
      else
        lineplot!(plt, xf[lf], yf[lf], color=:red)
      end
    end
  end
  title!(plt, "mesh")
  display(plt)
end

function rateandstate(V, psi, σn, ϕ, η, a, V0)
  Y = (1 ./ (2 .* V0)) .* exp.(psi ./ a)
  f = a .* asinh.(V .* Y)
  dfdV  = a .* (1 ./ sqrt.(1 + (V .* Y).^2)) .* Y

  g    = σn .* f    + η .* V - ϕ
  dgdV = σn .* dfdV + η
  (g, dgdV)
end

function newtbndv(func, xL, xR, x; ftol = 1e-6, maxiter = 500, minchange=0,
                  atolx = 1e-4, rtolx = 1e-4)
  (fL, _) = func(xL)
  (fR, _) = func(xR)
  if fL .* fR > 0
    return (typeof(x)(NaN), typeof(x)(NaN), -maxiter)
  end

  (f, df) = func(x)
  dxlr = xR - xL

  for iter = 1:maxiter
    dx = -f / df
    x  = x + dx

    if x < xL || x > xR || abs(dx) / dxlr < minchange
      x = (xR + xL) / 2
      dx = (xR - xL) / 2
    end

    (f, df) = func(x)

    if f * fL > 0
      (fL, xL) = (f, x)
    else
      (fR, xR) = (f, x)
    end
    dxlr = xR - xL

    if abs(f) < ftol && abs(dx) < atolx + rtolx * (abs(dx) + abs(x))
      return (x, f, iter)
    end
  end
  return (x, f, -maxiter)
end

max_stars_repo_name: brittany-erickson/HybridSBP

max_issues_repo_name: brittany-erickson/HybridSBP

max_forks_repo_name: brittany-erickson/HybridSBP


text:
@testset verbose = true "utility tests" begin
    @safetestset "Discard Burnin" begin
        using DifferentialEvolutionMCMC, Test, Random, Parameters, Distributions
        Random.seed!(29542)
        N = 10
        k = rand(Binomial(N, .5))
        data = (N = N,k = k)

        prior_loglike(θ) = logpdf(Beta(1, 1), θ)

        sample_prior() = rand(Beta(1, 1))

        bounds = ((0,1),)

        function loglike(data, θ)
            return logpdf(Binomial(data.N, θ), data.k)
        end

        names = (:θ,)

        model = DEModel(;
            sample_prior,
            prior_loglike,
            loglike,
            data,
            names
        )

        burnin = 1500
        n_iter = 3000

        de = DE(; sample_prior, Np=4, bounds, burnin, discard_burnin=false)

        chains = sample(model, de, n_iter)
        @test length(chains) == n_iter

        de = DE(; sample_prior, Np=4, bounds, burnin)
        chains = sample(model, de, n_iter)
        @test length(chains) == burnin
    end

    @safetestset "reset!" begin
        using DifferentialEvolutionMCMC, Test
        import DifferentialEvolutionMCMC: reset!

        p1 = Particle(Θ = [[.7,.5,.1],.4,.6])
        p2 = Particle(Θ = [[.9,.8,.5],.7,.8])
        idx = [[true,false,false],false,true]
        reset!(p1, p2, idx)

        @test p1.Θ[1][1] ≠ p2.Θ[1][1]
        @test p1.Θ[1][2] == p2.Θ[1][2]
        @test p1.Θ[1][3] == p2.Θ[1][3]
        @test p1.Θ[2] == p2.Θ[2]
        @test p1.Θ[3] ≠ p2.Θ[3]

        p1 = Particle(Θ = [[.7 .5;.1 .3],.4,.6])
        p2 = Particle(Θ = [[.9 .8;.5 .2],.7,.8])
        idx = [[true false; false true],false,true]
        reset!(p1, p2, idx)

        @test p1.Θ[1][1,1] ≠ p2.Θ[1][1,1]
        @test p1.Θ[1][1,2] == p2.Θ[1][1,2]
        @test p1.Θ[1][2,1] == p2.Θ[1][2,1]
        @test p1.Θ[1][2,2] ≠ p2.Θ[1][2,2]
        @test p1.Θ[1][3] == p2.Θ[1][3]
        @test p1.Θ[2] == p2.Θ[2]
        @test p1.Θ[3] ≠ p2.Θ[3]
    end

    @testset "projection" begin
        using Test, DifferentialEvolutionMCMC
        import DifferentialEvolutionMCMC: project

        # for example, see: https://www.youtube.com/watch?v=xSu-0xcRBo8&ab_channel=FireflyLectures
        proj(x1, x2) = (x1' * x2) / (x2' * x2) * x2
        x1 = [-1.0,4.0]
        x2 = [2.0,7.0]
        p1 = Particle(Θ = x1)
        p2 = Particle(Θ = x2)
        p3 = project(p1, p2)
        correct = proj(x1, x2)

        @test correct ≈ [52/53,182/53]
        @test p3.Θ ≈ correct

        x1 = [[-1.0,],4.0]
        x2 = [[2.0,],7.0]
        p1 = Particle(Θ = x1)
        p2 = Particle(Θ = x2)
        p3 = project(p1, p2)
        @test vcat(p3.Θ...) ≈ correct
    end

    @safetestset "Migration" begin
        using DifferentialEvolutionMCMC, Test, Random, Parameters, Distributions
        import DifferentialEvolutionMCMC: select_groups, select_particles, shift_particles!, sample_init

        Random.seed!(459) #Random.seed!(0451) #

        function equal(p1::Particle, p2::Particle)
            fields = fieldnames(Particle)
            for field in fields
                if getfield(p1, field) != getfield(p2, field)
                    println(field)
                    return false
                end
            end
            return true
        end

        N = 10
        k = rand(Binomial(N, .5))
        data = (N = N,k = k)

        prior_loglike(θ) = logpdf(Beta(1, 1), θ)

        sample_prior() = rand(Beta(1, 1))

        bounds = ((0,1),)

        function loglike(data, θ)
            return logpdf(Binomial(data.N, θ), data.k)
        end

        names = (:θ,)

        model = DEModel(;
            sample_prior,
            prior_loglike,
            loglike,
            data,
            names
        )

        burnin = 1500
        n_iter = 3000

        de = DE(; sample_prior, Np=4, bounds, burnin, discard_burnin=false)

        groups = sample_init(model, de, n_iter)
        sub_group = select_groups(de, groups)
        c_groups = deepcopy(groups)
        c_sub_group = deepcopy(sub_group)
        p_idx,particles = select_particles(sub_group)
        c_particles = deepcopy(particles)
        shift_particles!(sub_group, p_idx, particles)
        gidx = 1:length(sub_group)
        cidx = circshift(gidx, 1)
        cp_idx = circshift(p_idx, 1)
        for (i,c,p,cp) in zip(gidx, p_idx, cidx, cp_idx)
            @test sub_group[i][c].Θ == c_sub_group[p][cp].Θ
        end
        ridx = [4,3]
        for (i,c,r) in zip(1:2, p_idx[1:2], ridx)
            @test sub_group[i][c].Θ == groups[r][c].Θ
        end
    end

    @safetestset "particle operations" begin
        using DifferentialEvolutionMCMC, Test, Random, Distributions
        Random.seed!(29542)

        p1 = Particle(Θ = [1.,2.0])
        pr = p1 + 2
        @test pr.Θ ≈ [3,4]

        p1 = Particle(Θ = [1.,2.0])
        pr = p1 * 4
        @test pr.Θ ≈ [4,8]

        p1 = Particle(Θ = [1.,2.0])
        p2 = Particle(Θ = [1.,2.0])
        pr = p1 + p2
        @test pr.Θ ≈ [2,4]

        p1 = Particle(Θ = [1.,2.0])
        p2 = Particle(Θ = [1.,2.0])
        pr = 3 * (p1 + p2)
        @test pr.Θ ≈ [6,12]

        p1 = Particle(Θ = [1.,2.0])
        p2 = Particle(Θ = [-2.,3.0])
        pr = 3 * (p1 - p2)
        @test pr.Θ ≈ [9,-3]

        p1 = Particle(Θ = [1.,2.0])
        p2 = Particle(Θ = [-2.,3.0])
        p3 = Particle(Θ = [-2.,3.0])
        pr = 3 * (p1 - p2) + p3
        @test pr.Θ ≈ [7,0]

        p1 = Particle(Θ = [1.,2.0])
        b = Uniform(-.1, .1)
        pr = p1 + b
        @test p1.Θ ≈ pr.Θ atol = .2 # cummulative error
        @test p1.Θ ≠ pr.Θ
    end
end

ext: jl


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

mutable struct DnsNameAvailabilityResult <: SwaggerModel
    available::Any # spec type: Union{ Nothing, Bool } # spec name: available

    function DnsNameAvailabilityResult(;available=nothing)
        o = new()
        validate_property(DnsNameAvailabilityResult, Symbol("available"), available)
        setfield!(o, Symbol("available"), available)
        o
    end
end # type DnsNameAvailabilityResult

const _property_map_DnsNameAvailabilityResult = Dict{Symbol,Symbol}(Symbol("available")=>Symbol("available"))
const _property_types_DnsNameAvailabilityResult = Dict{Symbol,String}(Symbol("available")=>"Bool")
Base.propertynames(::Type{ DnsNameAvailabilityResult }) = collect(keys(_property_map_DnsNameAvailabilityResult))
Swagger.property_type(::Type{ DnsNameAvailabilityResult }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_DnsNameAvailabilityResult[name]))}
Swagger.field_name(::Type{ DnsNameAvailabilityResult }, property_name::Symbol) =  _property_map_DnsNameAvailabilityResult[property_name]

function check_required(o::DnsNameAvailabilityResult)
    true
end

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


text:
"""
    jacobian(fdm, f, x...)

Approximate the Jacobian of `f` at `x` using `fdm`. Results will be returned as a
`Matrix{<:Real}` of `size(length(y_vec), length(x_vec))` where `x_vec` is the flattened
version of `x`, and `y_vec` the flattened version of `f(x...)`. Flattening performed by
[`to_vec`](@ref).
"""
function jacobian(fdm, f, x::Vector{<:Real}; len=nothing)
    len !== nothing && Base.depwarn(
        "`len` keyword argument to `jacobian` is no longer required " *
        "and will not be permitted in the future.",
         :jacobian
    )
    ẏs = map(eachindex(x)) do n
        return fdm(zero(eltype(x))) do ε
            xn = x[n]
            x[n] = xn + ε
            ret = copy(first(to_vec(f(x))))  # copy required incase `f(x)` returns something that aliases `x`
            x[n] = xn  # Can't do `x[n] -= ϵ` as floating-point math is not associative
            return ret
        end
    end
    return (reduce(hcat, ẏs), )
end

function jacobian(fdm, f, x; len=nothing)
    x_vec, from_vec = to_vec(x)
    return jacobian(fdm, f ∘ from_vec, x_vec; len=len)
end

function jacobian(fdm, f, xs...; len=nothing)
    return ntuple(length(xs)) do k
        jacobian(fdm, x->f(replace_arg(x, xs, k)...), xs[k]; len=len)[1]
    end
end

replace_arg(x, xs::Tuple, k::Int) = ntuple(p -> p == k ? x : xs[p], length(xs))

"""
    _jvp(fdm, f, x::Vector{<:Real}, ẋ::AbstractVector{<:Real})

Convenience function to compute `jacobian(f, x) * ẋ`.
"""
function _jvp(fdm, f, x::Vector{<:Real}, ẋ::Vector{<:Real})
    return fdm(ε -> f(x .+ ε .* ẋ), zero(eltype(x)))
end

"""
    jvp(fdm, f, xẋs::Tuple{Any, Any}...)

Compute a Jacobian-vector product with any types of arguments for which [`to_vec`](@ref)
is defined. Each 2-`Tuple` in `xẋs` contains the value `x` and its tangent `ẋ`.
"""
function jvp(fdm, f, (x, ẋ)::Tuple{Any, Any})
    x_vec, vec_to_x = to_vec(x)
    _, vec_to_y = to_vec(f(x))
    return vec_to_y(_jvp(fdm, x_vec->to_vec(f(vec_to_x(x_vec)))[1], x_vec, to_vec(ẋ)[1]))
end
function jvp(fdm, f, xẋs::Tuple{Any, Any}...)
    x, ẋ = collect(zip(xẋs...))
    return jvp(fdm, xs->f(xs...), (x, ẋ))
end

"""
    j′vp(fdm, f, ȳ, x...)

Compute an adjoint with any types of arguments `x` for which [`to_vec`](@ref) is defined.
"""
function j′vp(fdm, f, ȳ, x)
    x_vec, vec_to_x = to_vec(x)
    ȳ_vec, _ = to_vec(ȳ)
    return (vec_to_x(_j′vp(fdm, x -> first(to_vec(f(vec_to_x(x)))), ȳ_vec, x_vec)), )
end

j′vp(fdm, f, ȳ, xs...) = j′vp(fdm, xs->f(xs...), ȳ, xs)[1]

function _j′vp(fdm, f, ȳ::Vector{<:Real}, x::Vector{<:Real})
    isempty(x) && return eltype(ȳ)[] # if x is empty, then so is the jacobian and x̄
    return transpose(first(jacobian(fdm, f, x))) * ȳ
end

"""
    grad(fdm, f, xs...)

Compute the gradient of `f` for any `xs` for which [`to_vec`](@ref) is defined.
"""
grad(fdm, f, xs...) = j′vp(fdm, f, 1, xs...)  # `j′vp` with seed of 1


text:
using BinaryBuilder

name = "xrootdgo"
version = v"0.31.1"

sources = [
    GitSource("https://github.com/go-hep/hep/",
              "c31af820c54f9d88ecdf40fe46a8bcf3a59f6f25"),
    DirectorySource("./bundled")
]

# Bash recipe for building across all platforms
script = raw"""
install_license ${WORKSPACE}/srcdir/LICENSE
cd $WORKSPACE/srcdir/

mkdir clib
cp $WORKSPACE/srcdir/main.go clib/main.go
mkdir -p ${libdir}
go get github.com/google/uuid
CGO_ENABLED=1 go build -buildmode=c-shared -o ${libdir}/xrootdgo.${dlext} clib/main.go
"""

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

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

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

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies;
               compilers=[:c, :go], julia_compat="1.6", preferred_gcc_version=v"6")


text:
"""
Cumulative Distribution Function of a Standard Gaussian Random Variable

		y=normcdf(x)

Where:\n
		x         = point of evaluation.

		y      = probability that a standard gaussian random variable is below x.

# Example
```julia-repl
julia> normcdf(0.0)
0.5
```
"""
function normcdf(x::number) where {number <: Number}
	return (1.0+erf(x/sqrt(2.0)))/2.0;
end

"""
Probability Distribution Function of a Standard Gaussian Random Variable

		y=normpdf(x)

Where:\n
		x         = point of evaluation.

		y      = value.

# Example
```julia-repl
julia> normpdf(0.0)
0.3989422804014327
```
"""
function normpdf(x::number) where {number <: Number}
	return exp(-0.5*x.^2)/sqrt(2*pi);
end

"""
Black & Scholes Price for European Options

		Price=blsprice(S0,K,r,T,σ,d=0.0,FlagIsCall=true)

Where:\n
		S0         = Value of the Underlying.
		K          = Strike Price of the Option.
		r          = Zero Rate.
		T          = Time to Maturity of the Option.
		σ          = Implied Volatility
		d          = Implied Dividend of the Underlying.
		FlagIsCall = true for Call Options, false for Put Options.

		Price      = price of the European Option.

# Example
```julia-repl
julia> blsprice(10.0,10.0,0.01,2.0,0.2,0.01)
1.1023600107733191
```
"""
function blsprice(S0::num1,K::num2,r::num3,T::num4,σ::num5,d::num6=0.0,FlagIsCall::Bool=true) where {num1 <: Number ,num2 <: Number ,num3 <: Number ,num4 <: Number ,num5 <: Number ,num6 <: Number}
	blscheck(S0,K,r,T,σ,d);
	d1=(log(S0/K)+(r-d+σ*σ*0.5)*T)/(σ*sqrt(T));
	d2=d1-σ*sqrt(T);
	Price=0.0;
	if FlagIsCall
	Price=S0*exp(-d*T)*normcdf(d1)-K*exp(-r*T)*normcdf(d2);
	else
	Price=-S0*exp(-d*T)*normcdf(-d1)+K*exp(-r*T)*normcdf(-d2);
	end
return Price;
end

"""
Black Price for European Options

		Price=blkprice(F0,K,r,T,σ,FlagIsCall=true)

Where:\n
		F0         = Value of the Forward.
		K          = Strike Price of the Option.
		r          = Zero Rate.
		T          = Time to Maturity of the Option.
		σ          = Implied Volatility.
		FlagIsCall = true for Call Options, false for Put Options.

		Price      = price of the European Option.

# Example
```julia-repl
julia> blkprice(10.0,10.0,0.01,2.0,0.2)
1.1023600107733191
```
"""
function blkprice(F0::num1,K::num2,r::num3,T::num4,σ::num5,FlagIsCall::Bool=true) where {num1 <: Number ,num2 <: Number ,num3 <: Number ,num4 <: Number ,num5 <: Number}
	blscheck(F0,K,r,T,σ);
	Price=blsprice(F0,K,r,T,σ,r,FlagIsCall);
return Price;
end

"""
Black & Scholes Delta for European Options

		Δ=blsdelta(S0,K,r,T,σ,d=0.0,FlagIsCall=true)

Where:\n
		S0         = Value of the Underlying.
		K          = Strike Price of the Option.
		r          = Zero Rate.
		T          = Time to Maturity of the Option.
		σ          = Implied Volatility
		d          = Implied Dividend of the Underlying.
		FlagIsCall = true for Call Options, false for Put Options.

		Δ          = delta of the European Option.

# Example
```julia-repl
julia> blsdelta(10.0,10.0,0.01,2.0,0.2,0.01)
0.5452173371920436
```
"""
function blsdelta(S0::num1,K::num2,r::num3,T::num4,σ::num5,d::num6=0.0,FlagIsCall::Bool=true) where {num1 <: Number ,num2 <: Number ,num3 <: Number ,num4 <: Number ,num5 <: Number ,num6 <: Number}
	blscheck(S0,K,r,T,σ,d);
	d1=(log(S0/K)+(r-d+σ*σ*0.5)*T)/(σ*sqrt(T));
	Δ=0.0;
	if FlagIsCall
	Δ=exp(-d*T)*normcdf(d1);
	else
	Δ=-exp(-d*T)*normcdf(-d1);
	end
return Δ;
end

"""
Black & Scholes Gamma for European Options

		Γ=blsgamma(S0,K,r,T,σ,d=0.0,FlagIsCall=true)

Where:\n
		S0         = Value of the Underlying.
		K          = Strike Price of the Option.
		r          = Zero Rate.
		T          = Time to Maturity of the Option.
		σ          = Implied Volatility
		d          = Implied Dividend of the Underlying.
		FlagIsCall = true for Call Options, false for Put Options.

		Γ          = gamma of the European Option.

# Example
```julia-repl
julia> blsgamma(10.0,10.0,0.01,2.0,0.2,0.01)
0.13687881535712826
```
"""
function blsgamma(S0::num1,K::num2,r::num3,T::num4,σ::num5,d::num6=0.0,FlagIsCall::Bool=true) where {num1 <: Number ,num2 <: Number ,num3 <: Number ,num4 <: Number ,num5 <: Number ,num6 <: Number}
	blscheck(S0,K,r,T,σ,d);
	d1=(log(S0/K)+(r-d+σ*σ*0.5)*T)/(σ*sqrt(T));
	Γ=exp(-d*T)*normpdf(d1)/(S0*σ*sqrt(T));
return Γ;
end

"""
Black & Scholes Vega for European Options

		ν=blsvega(S0,K,r,T,σ,d=0.0,FlagIsCall=true)

Where:\n
		S0         = Value of the Underlying.
		K          = Strike Price of the Option.
		r          = Zero Rate.
		T          = Time to Maturity of the Option.
		σ          = Implied Volatility
		d          = Implied Dividend of the Underlying.
		FlagIsCall = true for Call Options, false for Put Options.

		ν          = vega of the European Option.

# Example
```julia-repl
julia> blsvega(10.0,10.0,0.01,2.0,0.2,0.01)
5.475152614285131
```
"""
function blsvega(S0::num1,K::num2,r::num3,T::num4,σ::num5,d::num6=0.0,FlagIsCall::Bool=true) where {num1 <: Number ,num2 <: Number ,num3 <: Number ,num4 <: Number ,num5 <: Number ,num6 <: Number}
	blscheck(S0,K,r,T,σ,d);
	d1=(log(S0/K)+(r-d+σ*σ*0.5)*T)/(σ*sqrt(T));
	ν=S0*exp(-d*T)*normpdf(d1)*sqrt(T);
return ν;
end

"""
Black & Scholes Rho for European Options

		ρ=blsrho(S0,K,r,T,σ,d=0.0,FlagIsCall=true)

Where:\n
		S0         = Value of the Underlying.
		K          = Strike Price of the Option.
		r          = Zero Rate.
		T          = Time to Maturity of the Option.
		σ          = Implied Volatility
		d          = Implied Dividend of the Underlying.
		FlagIsCall = true for Call Options, false for Put Options.

		ρ          = rho of the European Option.

# Example
```julia-repl
julia> blsrho(10.0,10.0,0.01,2.0,0.2,0.01)
8.699626722294234
```
"""
function blsrho(S0::num1,K::num2,r::num3,T::num4,σ::num5,d::num6=0.0,FlagIsCall::Bool=true) where {num1 <: Number ,num2 <: Number ,num3 <: Number ,num4 <: Number ,num5 <: Number ,num6 <: Number}
	blscheck(S0,K,r,T,σ,d);
	d2=(log(S0/K)+(r-d-σ*σ*0.5)*T)/(σ*sqrt(T));
	if FlagIsCall
	ρ=K*exp(-r*T)*normcdf(d2)*T;
	else
		ρ=-K*exp(-r*T)*normcdf(-d2)*T;
	end
return ρ;
end

"""
Black & Scholes Theta for European Options

		Θ=blstheta(S0,K,r,T,σ,d=0.0,FlagIsCall=true)

Where:\n
		S0         = Value of the Underlying.
		K          = Strike Price of the Option.
		r          = Zero Rate.
		T          = Time to Maturity of the Option.
		σ          = Implied Volatility
		d          = Implied Dividend of the Underlying.
		FlagIsCall = true for Call Options, false for Put Options.

		Θ          = theta of the European Option.

# Example
```julia-repl
julia> blstheta(10.0,10.0,0.01,2.0,0.2,0.01)
-0.26273403060652334
```
"""
function blstheta(S0::num1,K::num2,r::num3,T::num4,σ::num5,d::num6=0.0,FlagIsCall::Bool=true) where {num1 <: Number ,num2 <: Number ,num3 <: Number ,num4 <: Number ,num5 <: Number ,num6 <: Number}
	blscheck(S0,K,r,T,σ,d);
	sqrtT       = sqrt(T);
	σ_sqrtT = σ .* sqrtT;

	d1 = (log(S0 ./ K) + (r - d + σ.^2 / 2) .* T)./ σ_sqrtT;

	shift = -exp(-d .* T) .* S0 .* normpdf(d1) .* σ / 2 ./ sqrtT;
	t1    = r .* K   .* exp(-r .* T);
	t2    = d .* S0 .* exp(-d .* T);
	Θ=0.0;
	if FlagIsCall
		Θ=shift - t1 .*      normcdf(d1 - σ_sqrtT)  + t2 .*  normcdf(d1)     ;
	else
		Θ=shift + t1 .* (1 - normcdf(d1 - σ_sqrtT)) + t2 .* (normcdf(d1) - 1);
	end
return Θ;
end

"""
Black & Scholes Lambda for European Options

		Λ=blslambda(S0,K,r,T,σ,d=0.0,FlagIsCall=true)

Where:\n
		S0         = Value of the Underlying.
		K          = Strike Price of the Option.
		r          = Zero Rate.
		T          = Time to Maturity of the Option.
		σ          = Implied Volatility
		d          = Implied Dividend of the Underlying.
		FlagIsCall = true for Call Options, false for Put Options.

		Λ          = lambda of the European Option.

# Example
```julia-repl
julia> blslambda(10.0,10.0,0.01,2.0,0.2,0.01)
4.945909973725978
```
"""
function blslambda(S0::num1,K::num2,r::num3,T::num4,σ::num5,d::num6=0.0,FlagIsCall::Bool=true) where {num1 <: Number ,num2 <: Number ,num3 <: Number ,num4 <: Number ,num5 <: Number ,num6 <: Number}
	blscheck(S0,K,r,T,σ,d);
	Price=blsprice(S0,K,r,T,σ,d,FlagIsCall);
	Δ=blsdelta(S0,K,r,T,σ,d,FlagIsCall);
	Λ=Δ*S0/Price;
return Λ;
end

"Check input for Black Scholes Formula"
function blscheck(S0::num1,K::num2,r::num3,T::num4,σ::num5,d::num6=0.0) where {num1 <: Number ,num2 <: Number ,num3 <: Number ,num4 <: Number ,num5 <: Number ,num6 <: Number}
lesseq(x::Complex,y::Complex) = real(x) <= real(y)
lesseq(x,y) = x <= y
if (lesseq(S0,zero(num1)))
	error("Spot Price Cannot Be Negative")
elseif (lesseq(K,zero(num2)))
	error("Strike Price Cannot Be Negative")
elseif (lesseq(T,zero(num4)))
	error("Time to Maturity Cannot Be Negative")
elseif (lesseq(σ,zero(num5)))
	error("Volatility Cannot Be Negative")
end
return;
end

"Brent Method: Scalar Equation Solver"
function brentMethod(f::Function, x0::Number, x1::Number,xtol::AbstractFloat=1e-14, ytol::AbstractFloat=1e-15)
		if xtol<0.0
		  throw(ErrorException("x tollerance cannot be negative"));
		end
	if ytol<0.0
		  throw(ErrorException("y tollerance cannot be negative"));
		end
		EPS = eps(x0)
	maxiter=80;
		y0 = f(x0)
		y1 = f(x1)
		if abs(y0) < abs(y1)
		    # Swap lower and upper bounds.
		    x0, x1 = x1, x0
		    y0, y1 = y1, y0
		end
		x2 = x0
		y2 = y0
		x3 = x2
		bisection = true
		for _ in 1:maxiter
		    # x-tolerance.
		    if abs(x1-x0) < xtol
		        return x1
		    end

		    # Use inverse quadratic interpolation if f(x0)!=f(x1)!=f(x2)
		    # and linear interpolation (secant method) otherwise.
		    if abs(y0-y2) > ytol && abs(y1-y2) > ytol
		        x = x0*y1*y2/((y0-y1)*(y0-y2)) +
		            x1*y0*y2/((y1-y0)*(y1-y2)) +
		            x2*y0*y1/((y2-y0)*(y2-y1))
		    else
		        x = x1 - y1 * (x1-x0)/(y1-y0)
		    end

		    # Use bisection method if satisfies the conditions.
		    delta = abs(2EPS*abs(x1))
		    min1 = abs(x-x1)
		    min2 = abs(x1-x2)
		    min3 = abs(x2-x3)
		    if (x < (3x0+x1)/4 && x > x1) ||
		       (bisection && min1 >= min2/2) ||
		       (!bisection && min1 >= min3/2) ||
		       (bisection && min2 < delta) ||
		       (!bisection && min3 < delta)
		        x = (x0+x1)/2
		        bisection = true
		    else
		        bisection = false
		    end
		    y = f(x)
		    # y-tolerance.
		    if abs(y) < ytol
		        return x
		    end
		    x3 = x2
		    x2 = x1
		    if sign(y0) != sign(y)
		        x1 = x
		        y1 = y
		    else
		        x0 = x
		        y0 = y
		    end
		    if abs(y0) < abs(y1)
		        # Swap lower and upper bounds.
		        x0, x1 = x1, x0
		        y0, y1 = y1, y0
		    end
		end
		throw(ErrorException("Max iteration exceeded, possible wrong result"));
end

export blsimpv
"""
Black & Scholes Implied Volatility for European Options

		σ=blsimpv(S0,K,r,T,Price,d=0.0,FlagIsCall=true,xtol=1e-14,ytol=1e-15)

Where:\n
		S0         = Value of the Underlying.
		K          = Strike Price of the Option.
		r          = Zero Rate.
		T          = Time to Maturity of the Option.
		Price      = Price of the Option.
		d          = Implied Dividend of the Underlying.
		FlagIsCall = true for Call Options, false for Put Options.

		σ          = implied volatility of the European Option.

# Example
```julia-repl
julia> blsimpv(10.0,10.0,0.01,2.0,2.0)
0.3433730534290586
```
"""
function blsimpv(S0::num1,K::num2,r::num3,T::num4,Price::num5,d::num6=0.0,FlagIsCall::Bool=true,xtol::Real=1e-14,ytol::Real=1e-15) where {num1 <: Real ,num2 <: Real ,num3 <: Real ,num4 <: Real ,num5 <: Real ,num6 <: Real}
if (Price< num5(0))
	throw(ErrorException("Option Price Cannot Be Negative"));
end
blscheck(S0,K,r,T,0.1,d);
f(x)=(blsprice(S0,K,r,T,x,d,FlagIsCall)-Price);
σ=brentMethod(f,0.001,1.2,xtol,ytol);

return σ;
end

export blkimpv
"""
Black Implied Volatility for European Options

		σ=blkimpv(F0,K,r,T,Price,FlagIsCall=true,xtol=1e-14,ytol=1e-15)

Where:\n
		F0         = Value of the Forward.
		K          = Strike Price of the Option.
		r          = Zero Rate.
		T          = Time to Maturity of the Option.
		Price      = Price of the Option.
		FlagIsCall = true for Call Options, false for Put Options.

		σ          = implied volatility of the European Option.

# Example
```julia-repl
julia> blkimpv(10.0,10.0,0.01,2.0,2.0)
0.36568658096623635
```
"""
function blkimpv(F0::num1,K::num2,r::num3,T::num4,Price::num5,FlagIsCall::Bool=true,xtol::Real=1e-14,ytol::Real=1e-15) where {num1 <: Real ,num2 <: Real ,num3 <: Real ,num4 <: Real ,num5 <: Real}
	blscheck(F0,K,r,T,0.1);
	σ=blsimpv(F0,K,r,T,Price,r,FlagIsCall,xtol,ytol)
	return σ;
end

##	ADDITIONAL Function

"""
Black & Scholes Psi for European Options

		Ψ=blspsi(S0,K,r,T,σ,d=0.0,FlagIsCall=true)

Where:\n
		S0         = Value of the Underlying.
		K          = Strike Price of the Option.
		r          = Zero Rate.
		T          = Time to Maturity of the Option.
		σ          = Implied Volatility
		d          = Implied Dividend of the Underlying.
		FlagIsCall = true for Call Options, false for Put Options.

		Ψ          = psi of the European Option.

# Example
```julia-repl
julia> blspsi(10.0,10.0,0.01,2.0,0.2,0.01)
-10.904346743840872
```
"""
function blspsi(S0::num1,K::num2,r::num3,T::num4,σ::num5,d::num6=0.0,FlagIsCall::Bool=true) where {num1 <: Number ,num2 <: Number ,num3 <: Number ,num4 <: Number ,num5 <: Number ,num6 <: Number}
	blscheck(S0,K,r,T,σ,d);
	d1=(log(S0/K)+(r-d+σ*σ*0.5)*T)/(σ*sqrt(T));
	if FlagIsCall
	Ψ=-S0*exp(-d*T)*normcdf(d1)*T;
	else
		Ψ=S0*exp(-d*T)*normcdf(-d1)*T;
	end
return Ψ;
end

"""
Black & Scholes Vanna for European Options

		Vanna=blsvanna(S0,K,r,T,σ,d=0.0,FlagIsCall=true)

Where:\n
		S0         = Value of the Underlying.
		K          = Strike Price of the Option.
		r          = Zero Rate.
		T          = Time to Maturity of the Option.
		σ          = Implied Volatility
		d          = Implied Dividend of the Underlying.
		FlagIsCall = true for Call Options, false for Put Options.

		Vanna        = vanna of the European Option.

# Example
```julia-repl
julia> blsvanna(10.0,10.0,0.01,2.0,0.2,0.01)
0.2737576307142566
```
"""
function blsvanna(S0::num1,K::num2,r::num3,T::num4,σ::num5,d::num6=0.0,FlagIsCall::Bool=true) where {num1 <: Number ,num2 <: Number ,num3 <: Number ,num4 <: Number ,num5 <: Number ,num6 <: Number}
	blscheck(S0,K,r,T,σ,d);
	d1=(log(S0/K)+(r-d+σ*σ*0.5)*T)/(σ*sqrt(T));
	d2=d1-σ*sqrt(T);
	Vanna=-exp(-d*T)*normpdf(d1)*d2/σ;
return Vanna;
end


text:
using Test
using LibPNG
using ColorTypes
using FixedPointNumbers
using ImageCore
using Logging
using Random

logger = ConsoleLogger(stdout, Logging.Debug)
global_logger(logger)

include("images.jl")

tmpdir = joinpath(tempdir(), "LibPNG")
try
    rm(tmpdir, recursive = true)
catch
    @error "Unable to remove temp directory at: $(tmpdir)"
end

max_stars_repo_name: FugroRoames/LibPNG

max_issues_repo_name: FugroRoames/LibPNG

max_forks_repo_name: FugroRoames/LibPNG


text:
using BenchmarkTools

using SIMD
using ChipSort

const SUITE = BenchmarkGroup()

function randa(T, K...)
    data = reshape(valloc(T, div(32, sizeof(T)), prod(K)), K...)
    data .= rand(T, K...)
    data
end

T=Int32
K=2^7
V=2^3
J=2^3
len = V*J*K
SUITE["Array8kInt32"] = BenchmarkGroup(["size-medium", "Int32"])
SUITE["Array8kInt32"]["JuliaStd"] = @benchmarkable sort!(data) setup=(data = randa($T, $len))
SUITE["Array8kInt32"]["ChipSort"] = @benchmarkable chipsort_medium!(data, Val(V), Val(J), Val(K)) setup=(data = randa($T, $len))
SUITE["Array8kInt32"]["CombSort"] = @benchmarkable combsort!(data) setup=(data = randa($T, $len))
SUITE["Array8kInt32"]["InsertionSort"] = @benchmarkable insertion_sort!(data) setup=(data = randa($T, $len))

Km=2^5
Vm=8
Jm=8
vj=Vm*Jm
multi_sort!(data) = for n in 1:Km
    sort!(@view data[:,n])
end
multi_chipsort_small!(data) = for n in 1:Km
    chipsort_small!((@view data[:,n]), Val(8), Val(8))
end
multi_combsort!(data) = for n in 1:Km
    combsort!(@view data[:,n])
end
multi_insertion_sort!(data) = for n in 1:Km
    insertion_sort!(@view data[:,n])
end
SUITE["Array$(Km)x$vj$T"] = BenchmarkGroup(["size-small", "Int32"])
SUITE["Array$(Km)x$vj$T"]["JuliaStd"] = @benchmarkable multi_sort!(data) setup=(data = randa($T, $vj, $Km))
SUITE["Array$(Km)x$vj$T"]["ChipSort"] = @benchmarkable multi_chipsort_small!(data) setup=(data = randa($T, $vj, $Km))
SUITE["Array$(Km)x$vj$T"]["CombSort"] = @benchmarkable multi_combsort!(data) setup=(data = randa($T, $vj, $Km))
SUITE["Array$(Km)x$vj$T"]["InsertionSort"] = @benchmarkable multi_insertion_sort!(data) setup=(data = randa($T, $vj, $Km))

T=Int32
Ka=2^5
Va=2^3
Ja=2^12
lena = Va*Ja*Ka
SUITE["Array1MInt32"] = BenchmarkGroup(["size-large", "Int32"])
SUITE["Array1MInt32"]["JuliaStd"] = @benchmarkable sort!(data) setup=(data = randa($T, $lena))
SUITE["Array1MInt32"]["ChipSort"] = @benchmarkable chipsort_large(data, Val(Va), Val(Ka)) setup=(data = randa($T, $lena))
# SUITE["Array1MInt32"]["CombSort"] = @benchmarkable combsort!(data) setup=(data = randa($T, $len))
# SUITE["Array1MInt32"]["InsertionSort"] = @benchmarkable insertion_sort!(data) setup=(data = randa($T, $len))


text:
using Test
using QuantumOptics
using LinearAlgebra

@testset "steadystate" begin

ωc = 1.2
ωa = 0.9
g = 1.0
γ = 0.5
κ = 1.1
η = 1.5

T = Float64[0.,1.]

fockbasis = FockBasis(10)
spinbasis = SpinBasis(1//2)
basis = tensor(spinbasis, fockbasis)

sx = sigmax(spinbasis)
sy = sigmay(spinbasis)
sz = sigmaz(spinbasis)
sp = sigmap(spinbasis)
sm = sigmam(spinbasis)

Ha = embed(basis, 1, 0.5*ωa*sz)
Hc = embed(basis, 2, ωc*number(fockbasis))
Hint = sm ⊗ create(fockbasis) + sp ⊗ destroy(fockbasis)
H = Ha + Hc + Hint
Hdense = dense(H)

Ja = embed(basis, 1, sqrt(γ)*sm)
Ja2 = embed(basis, 1, sqrt(0.5*γ)*sp)
Jc = embed(basis, 2, sqrt(κ)*destroy(fockbasis))
J = [Ja, Jc]
Jdense = map(dense, J)

Ψ₀ = spinup(spinbasis) ⊗ fockstate(fockbasis, 2)
ρ₀ = dm(Ψ₀)
ψ0_p = fockstate(fockbasis, 0)
ρ0_p = dm(ψ0_p)

tout, ρt = timeevolution.master([0,100], ρ₀, Hdense, Jdense; reltol=1e-7)

tss, ρss = steadystate.master(Hdense, Jdense; tol=1e-4)
@test tracedistance(ρss[end], ρt[end]) < 1e-3

ρss = steadystate.eigenvector(Hdense, Jdense)
@test tracedistance(ρss, ρt[end]) < 1e-6

ρss = steadystate.eigenvector(Hdense, Jdense, nev = 1)
@test tracedistance(ρss, ρt[end]) < 1e-6

ρss = steadystate.eigenvector(H, sqrt(2).*J; rates=0.5.*ones(length(J)))
@test tracedistance(ρss, ρt[end]) < 1e-3

ρss = steadystate.eigenvector(H, sqrt(2).*J; rates=0.5.*ones(length(J)), nev = 1)
@test tracedistance(ρss, ρt[end]) < 1e-3

ρss = steadystate.eigenvector(liouvillian(Hdense, Jdense))
@test tracedistance(ρss, ρt[end]) < 1e-6

@test_throws TypeError steadystate.eigenvector(H, J; ncv="a")

ev, ops = steadystate.liouvillianspectrum(Hdense, Jdense)
@test tracedistance(ρss, ops[1]) < 1e-12
@test ev[sortperm(abs.(ev))] == ev

ev, ops = steadystate.liouvillianspectrum(H, sqrt(2).*J; rates=0.5.*ones(length(J)), nev = 1)
@test tracedistance(ρss, ops[1]/tr(ops[1])) < 1e-12
@test ev[sortperm(abs.(ev))] == ev

# Compute steady-state photon number of a driven cavity (analytically: η^2/κ^2)
Hp = η*(destroy(fockbasis) + create(fockbasis))
Jp = [sqrt(2κ)*destroy(fockbasis)]
n_an = η^2/κ^2

tss,ρss = steadystate.master(Hp, Jp; rho0=ρ0_p, tol=1e-4)
nss = expect(create(fockbasis)*destroy(fockbasis), ρss[end])
@test n_an - real(nss) < 1e-3

ρss = steadystate.eigenvector(Hp, Jp)
nss = expect(create(fockbasis)*destroy(fockbasis), ρss)
@test n_an - real(nss) < 1e-3

ρss = steadystate.eigenvector(dense(Hp), map(dense, Jp))
nss = expect(create(fockbasis)*destroy(fockbasis), ρss)
@test n_an - real(nss) < 1e-3

# Test error messages
#@test_throws ErrorException steadystate.eigenvector(sx, [sm])
function fout_wrong(t, x)
  @assert x == t
end
@test_throws AssertionError steadystate.master(Hdense, Jdense; fout=fout_wrong)

end # testset


text:
"""
    LegendreDerivativeOperator{T<:Real, Grid}

A derivative operator on a nonperiodic Lobatto-Legendre grid with scalar type
`T` computing the first derivative using a Legendre expansion.
"""
struct LegendreDerivativeOperator{T<:Real} <: AbstractDerivativeOperator{T}
    jac::T
    Δx::T
    grid::Vector{T} # N nodes, including the left and the right boundary
    basis::LobattoLegendre{T}

    function LegendreDerivativeOperator(xmin::T, xmax::T, basis::LobattoLegendre{T}) where {T<:Real}
        grid = map_from_canonical.(basis.nodes, xmin, xmax, basis)
        jac = 2 / (xmax - xmin)
        Δx = inv(jac)

        new{T}(jac, Δx, grid, basis)
    end
end

"""
    LegendreDerivativeOperator(xmin::T, xmax::T, N::Int) where {T<:Real}

Construct the `LegendreDerivativeOperator` on a uniform grid between `xmin` and
`xmax` using `N` nodes and `N-1` Legendre modes.
"""
function LegendreDerivativeOperator(xmin::T, xmax::T, N::Int) where {T<:Real}
    @argcheck N >= 2

    basis = LobattoLegendre(N-1, T)

    LegendreDerivativeOperator(xmin, xmax, basis)
end

function legendre_derivative_operator(xmin::T, xmax::T, N::Int) where {T<:Real}
    LegendreDerivativeOperator(xmin, xmax, N)
end

derivative_order(D::LegendreDerivativeOperator) = 1
LinearAlgebra.issymmetric(D::LegendreDerivativeOperator) = false
integrate(func, u, D::LegendreDerivativeOperator) = D.Δx*integrate(func, u, D.basis)
function evaluate_coefficients(u, D::LegendreDerivativeOperator,
                               npoints=2*size(D,2)+1)
    evaluate_coefficients(u, D.basis, npoints)
end

function Base.show(io::IO, D::LegendreDerivativeOperator{T}) where {T}
    x = grid(D)
    print(io, "First derivative operator {T=", T, "} \n")
    print(io, "on the Lobatto Legendre nodes in [", first(x), ", ", last(x),
                "] using ", length(x), " nodes. \n")
end

function mul!(dest::AbstractVector{T}, D::LegendreDerivativeOperator, u::AbstractVector{T}) where {T}
    @unpack jac, basis = D
    N, _ = size(D)
    @boundscheck begin
        @argcheck N == length(u)
        @argcheck N == length(dest)
    end

    mul!(dest, basis.D, u)
    dest .*= jac

    nothing
end

function left_boundary_weight(D::LegendreDerivativeOperator)
    @inbounds retval = D.Δx * D.basis.weights[1]
    retval
end

function right_boundary_weight(D::LegendreDerivativeOperator)
    @inbounds retval = D.Δx * D.basis.weights[end]
    retval
end

function ConstantFilter(basis::LobattoLegendre{T}, filter, tmp::Array) where {T}
    Np1 = length(basis.nodes)
    coefficients = Array{T}(undef, Np1)
    set_filter_coefficients!(coefficients, filter)
    modal2nodal = legendre_vandermonde(basis)
    nodal2modal = FactorisationWrapper(factorize(modal2nodal))

    ConstantFilter(coefficients, nodal2modal, modal2nodal, tmp, filter)
end

function ConstantFilter(basis::LobattoLegendre{T}, filter, TmpEltype::DataType=T) where {T}
    Np1 = length(basis.nodes)
    tmp = Array{TmpEltype}(undef, Np1)

    ConstantFilter(basis, filter, tmp)
end

"""
    ConstantFilter(D::LegendreDerivativeOperator, filter, TmpEltype=T)

Create a modal filter with constant parameters adapted to the Legendre
derivative operator `D` with parameters given by the filter function `filter`.
"""
function ConstantFilter(D::LegendreDerivativeOperator{T}, filter, TmpEltype=T) where {T}
    ConstantFilter(D.basis, filter, TmpEltype)
end

# TODO: LegendreSuperSpectralViscosity


text:
## object-oriented Regex interface ##

include("pcre.jl")

const DEFAULT_OPTS = PCRE.JAVASCRIPT_COMPAT | PCRE.UTF8 | PCRE.NO_UTF8_CHECK

type Regex
    pattern::ByteString
    options::Uint32
    regex::Ptr{Void}
    extra::Ptr{Void}

    function Regex(pattern::String, options::Integer)
        pattern = bytestring(pattern)
        options = uint32(options)
        if (options & ~PCRE.OPTIONS_MASK) != 0
            error("invalid regex options: $options")
        end
        re = compile(new(pattern, options, C_NULL, C_NULL))
        finalizer(re,
            function(re::Regex)
                re.extra != C_NULL && PCRE.free_study(re.extra)
                re.regex != C_NULL && PCRE.free(re.regex)
            end)
        re
    end
end

function Regex(pattern::String, flags::String)
    options = DEFAULT_OPTS
    for f in flags
        options |= f=='i' ? PCRE.CASELESS  :
                   f=='m' ? PCRE.MULTILINE :
                   f=='s' ? PCRE.DOTALL    :
                   f=='x' ? PCRE.EXTENDED  :
                   error("unknown regex flag: $f")
    end
    Regex(pattern, options)
end
Regex(pattern::String) = Regex(pattern, DEFAULT_OPTS)

function compile(regex::Regex)
    if regex.regex == C_NULL
        regex.regex = PCRE.compile(regex.pattern, regex.options & PCRE.COMPILE_MASK)
        regex.extra = PCRE.study(regex.regex, PCRE.STUDY_JIT_COMPILE)
    end
    regex
end

macro r_str(pattern, flags...) Regex(pattern, flags...) end
macro r_mstr(pattern, flags...) Regex(pattern, flags...) end

copy(r::Regex) = r

function show(io::IO, re::Regex)
    imsx = PCRE.CASELESS|PCRE.MULTILINE|PCRE.DOTALL|PCRE.EXTENDED
    if (re.options & ~imsx) == DEFAULT_OPTS
        print(io, 'r')
        print_quoted_literal(io, re.pattern)
        if (re.options & PCRE.CASELESS ) != 0; print(io, 'i'); end
        if (re.options & PCRE.MULTILINE) != 0; print(io, 'm'); end
        if (re.options & PCRE.DOTALL   ) != 0; print(io, 's'); end
        if (re.options & PCRE.EXTENDED ) != 0; print(io, 'x'); end
    else
        print(io, "Regex(")
        show(io, re.pattern)
        print(io, ',')
        show(io, re.options)
        print(io, ')')
    end
end

# TODO: map offsets into non-ByteStrings back to original indices.
# or maybe it's better to just fail since that would be quite slow

immutable RegexMatch
    match::SubString{UTF8String}
    captures::Vector{Union(Nothing,SubString{UTF8String})}
    offset::Int
    offsets::Vector{Int}
end

function show(io::IO, m::RegexMatch)
    print(io, "RegexMatch(")
    show(io, m.match)
    if !isempty(m.captures)
        print(io, ", ")
        for i = 1:length(m.captures)
            print(io, i, "=")
            show(io, m.captures[i])
            if i < length(m.captures)
                print(io, ", ")
            end
        end
    end
    print(io, ")")
end

function ismatch(r::Regex, s::String, offset::Integer=0)
    compile(r)
    PCRE.exec(r.regex, r.extra, bytestring(s), offset, r.options & PCRE.EXECUTE_MASK, false)
end
function ismatch(r::Regex, s::SubString, offset::Integer=0)
    compile(r)
    PCRE.exec(r.regex, r.extra, s, offset, r.options & PCRE.EXECUTE_MASK, false)
end

function match(re::Regex, str::UTF8String, idx::Integer, add_opts::Uint32=uint32(0))
    opts = re.options & PCRE.EXECUTE_MASK | add_opts
    compile(re)
    m, n = PCRE.exec(re.regex, re.extra, str, idx-1, opts, true)
    if isempty(m); return nothing; end
    mat = SubString(str, m[1]+1, m[2])
    cap = Union(Nothing,SubString{UTF8String})[
            m[2i+1] < 0 ? nothing : SubString(str, m[2i+1]+1, m[2i+2]) for i=1:n ]
    off = Int[ m[2i+1]::Int32+1 for i=1:n ]
    RegexMatch(mat, cap, m[1]+1, off)
end

match(re::Regex, str::Union(ByteString,SubString), idx::Integer, add_opts::Uint32=uint32(0)) =
    match(re, utf8(str), idx, add_opts)

match(r::Regex, s::String) = match(r, s, start(s))
match(r::Regex, s::String, i::Integer) =
    error("regex matching is only available for bytestrings; use bytestring(s) to convert")

function matchall(re::Regex, str::UTF8String, overlap::Bool=false)
    regex = compile(re).regex
    extra = re.extra
    n = length(str.data)
    matches = SubString{UTF8String}[]
    offset = int32(0)
    opts = re.options & PCRE.EXECUTE_MASK
    opts_nonempty = opts | PCRE.ANCHORED | PCRE.NOTEMPTY_ATSTART
    prevempty = false
    ovec = Array(Int32, 3)
    while true
        result = ccall((:pcre_exec, :libpcre), Int32,
                       (Ptr{Void}, Ptr{Void}, Ptr{Uint8}, Int32,
                       Int32, Int32, Ptr{Int32}, Int32),
                       regex, extra, str, n,
                       offset, prevempty ? opts_nonempty : opts, ovec, 3)

        if result < 0
            if prevempty && offset < n
                offset = int32(nextind(str, offset + 1) - 1)
                prevempty = false
                continue
            else
                break
            end
        end

        push!(matches, SubString(str, ovec[1]+1, ovec[2]))
        prevempty = offset == ovec[2]
        if overlap
            if !prevempty
                offset = int32(nextind(str, offset + 1) - 1)
            end
        else
            offset = ovec[2]
        end
    end
    matches
end

matchall(re::Regex, str::Union(ByteString,SubString), overlap::Bool=false) =
    matchall(re, utf8(str), overlap)

function search(str::Union(ByteString,SubString), re::Regex, idx::Integer)
    if idx > nextind(str,endof(str))
        throw(BoundsError())
    end
    opts = re.options & PCRE.EXECUTE_MASK
    compile(re)
    m, n = PCRE.exec(re.regex, re.extra, str, idx-1, opts, true)
    isempty(m) ? (0:-1) : ((m[1]+1):prevind(str,m[2]+1))
end
search(s::String, r::Regex, idx::Integer) =
    error("regex search is only available for bytestrings; use bytestring(s) to convert")
search(s::String, r::Regex) = search(s,r,start(s))

immutable RegexMatchIterator
    regex::Regex
    string::UTF8String
    overlap::Bool

    function RegexMatchIterator(regex::Regex, string::String, ovr::Bool=false)
        new(regex, string, ovr)
    end
end
compile(itr::RegexMatchIterator) = (compile(itr.regex); itr)
eltype(itr::RegexMatchIterator) = RegexMatch
start(itr::RegexMatchIterator) = match(itr.regex, itr.string, 1, uint32(0))
done(itr::RegexMatchIterator, prev_match) = (prev_match == nothing)

# Assumes prev_match is not nothing
function next(itr::RegexMatchIterator, prev_match)
    prevempty = isempty(prev_match.match)

    if itr.overlap
        if !prevempty
            offset = nextind(itr.string, prev_match.offset)
        else
            offset = prev_match.offset
        end
    else
        offset = prev_match.offset + endof(prev_match.match)
    end

    opts_nonempty = uint32(PCRE.ANCHORED | PCRE.NOTEMPTY_ATSTART)
    while true
        mat = match(itr.regex, itr.string, offset,
                    prevempty ? opts_nonempty : uint32(0))

        if mat === nothing
            if prevempty && offset <= length(itr.string.data)
                offset = nextind(itr.string, offset)
                prevempty = false
                continue
            else
                break
            end
        else
            return (prev_match, mat)
        end
    end
    (prev_match, nothing)
end

function eachmatch(re::Regex, str::String, ovr::Bool=false)
    RegexMatchIterator(re,str,ovr)
end

eachmatch(re::Regex, str::String) = RegexMatchIterator(re,str)

# miscellaneous methods that depend on Regex being defined

filter!(r::Regex, v) = filter!(x->ismatch(r,x), v)
filter(r::Regex, v)  = filter(x->ismatch(r,x), v)

filter!(r::Regex, d::Dict) = filter!((k,v)->ismatch(r,k),d)
filter(r::Regex,  d::Dict) = filter!(r,copy(d))

# Don't serialize the pointers
function serialize(s, r::Regex)
    serialize_type(s, typeof(r))
    serialize(s, r.pattern)
    serialize(s, r.options)
end

function deserialize(s, t::Type{Regex})
    pattern = deserialize(s)
    options = deserialize(s)
    Regex(pattern, options)
end

max_stars_repo_name: porterjamesj/julia

max_issues_repo_name: porterjamesj/julia

max_forks_repo_name: porterjamesj/julia


text:
module FrostHelperDirectionalPuffer

using ..Ahorn, Maple

@mapdef Entity "FrostHelper/DirectionalPuffer" DirectionalPuffer(x::Integer, y::Integer, right::Bool=true, directory::String="objects/puffer/", explodeDirection::String="Right", dashRecovery::Integer=1, static::Bool=false, noRespawn::Bool=false, respawnTime::Number=2.5)

const placements = Ahorn.PlacementDict(
    "Directional Puffer (FrostHelper, Right)" => Ahorn.EntityPlacement(
        DirectionalPuffer,
        "point",
        Dict{String, Any}(
            "right" => true,
            "explodeDirection" => "Right"
        )
    ),
    "Directional Puffer (FrostHelper, Left)" => Ahorn.EntityPlacement(
        DirectionalPuffer,
        "point",
        Dict{String, Any}(
            "right" => false,
            "explodeDirection" => "Left"
        )
    ),
    "Directional Puffer (FrostHelper, None)" => Ahorn.EntityPlacement(
        DirectionalPuffer,
        "point",
        Dict{String, Any}(
            "explodeDirection" => "None"
        )
    )
)

const explodeDirections = String[
    "Left",
    "Right",
    "Both",
    "None",
]

Ahorn.editingOptions(entity::DirectionalPuffer) = Dict{String, Any}(
    "explodeDirection" => explodeDirections
)

function getSprite(entity::DirectionalPuffer)
    return string(get(entity, "directory", "objects/puffer/"), "idle00")
end

function Ahorn.selection(entity::DirectionalPuffer)
    x, y = Ahorn.position(entity)
    scaleX = get(entity, "right", false) ? 1 : -1

    return Ahorn.getSpriteRectangle("objects/puffer/idle00", x, y, sx=scaleX)
end

function Ahorn.render(ctx::Ahorn.Cairo.CairoContext, entity::DirectionalPuffer, room::Maple.Room)
    scaleX = get(entity, "right", false) ? 1 : -1

    Ahorn.drawSprite(ctx, getSprite(entity), 0, 0, sx=scaleX)
end

function Ahorn.flipped(entity::DirectionalPuffer, horizontal::Bool)
    if horizontal
        entity.right = !entity.right

        return entity
    end
end

end

max_stars_repo_name: JaThePlayer/FrostHelper

max_issues_repo_name: JaThePlayer/FrostHelper

max_forks_repo_name: JaThePlayer/FrostHelper


text:
####################
# Type Definitions #
####################
    abstract AbstractBasis{S<:AbstractStructure} <: AbstractQuantum{S}
    abstract AbstractFiniteBasis{S<:AbstractStructure} <: AbstractBasis{S}
    abstract AbstractInfiniteBasis{S<:AbstractStructure} <: AbstractBasis{S}

#############
# Functions #
#############
    # Note: All B<:AbstractBasis types should implement the following:
    #
    # checkcoeffs(coeffs, dim, basis::B) -> checks whether a coefficient
    #                                       array is valid for use with
    #                                       the given basis. This is used
    #                                       by QuArray to ensure that the
    #                                       coefficient array is not malformed
    #                                       with respect to the input bases.
    #                                       The second argument, `dim`, specifies
    #                                       the dimension of the coefficient
    #                                       array which corresponds to the provided
    #                                       basis.
    #
    # nfactors(basis::B) -> the number of factor bases for `basis`
    #
    # tensor(a::B, b::B) -> Take the tensor product of these two bases. This
    #                       function should optimally return a basis of same
    #                       type as the input bases. We can and should also
    #                       implement promote_rules/conversion methods between
    #                       bases.
    #
    # structure{S}(::Type{MyBasisType{S}}) -> returns S<:AbstractStructure for
    #                                         the provided basis type.

    checkcoeffs(coeffs, dim, basis::AbstractBasis) = error("checkcoeffs(coeffs, dim, ::$(typeof(basis))) must be defined!")

    for basis=(:AbstractBasis, :AbstractFiniteBasis, :AbstractInfiniteBasis)
        @eval begin
            structure{S}(::Type{($basis){S}}) = S
        end
    end

######################
# Include Statements #
######################
    include("finitebasis.jl")
    include("labelbasis.jl")

export AbstractBasis,
    AbstractFiniteBasis,
    AbstractInfiniteBasis,
    checkcoeffs,
    structure


text:
export jdqr

import LinearAlgebra.GeneralizedSchur

"""
Compute an approximate partial Schur decomposition of a square matrix A

    partial_schur, _ = jdqr(A, pairs = 5)

Accepts the following keyword arguments:

- `pairs`: number of Schur vectors and values requested
- `subspace_dimensions`: a range `min:max` that determines the size of the search subspace.
  In every outer iteration the search subspace is expanded to the `max` size and then
  trimmed to the `min` size. When the search subspace is trimmed, the best approximate Schur
  vectors are kept. A large search subspace speeds up convergence, but is computationally
  more expensive. It is advisable to set `min` a bit larger than `pairs`.
- `max_iter`: maximum number of steps. If the algorithm exceeds `max_iter` steps, it might
  not have converged for all `pairs` of Schur vectors.
- `target`: determines which eigenvalues to find. Options: Near(σ) where σ is a complex
  number in the complex plane near which we want to find eigenvalues; LargestMagnitude(σ),
  SmallestMagnitude(σ), LargestRealPart(σ), SmallestRealPart(σ), LargestImaginaryPart(σ),
  SmallestImaginaryPart(σ), where σ is a complex number / an educated guess.
- `tolerance`: the accuracy at which Schur vectors are considered converged, based on
  the residual norm.
- `solver`: An iterative correction equation solver used for the expansion of the search
  subspace.
- `verbosity=0` silent, `verbosity==1` prints a progress bar, `verbosity==2` prints more detailed debug output.
"""
function jdqr(
    A;
    solver::CorrectionSolver = BiCGStabl(size(A, 1), max_mv_products = 10, l = 2),
    pairs::Integer = 5,
    subspace_dimensions::AbstractUnitRange{<:Integer} = 10:20,
    max_iter::Integer = 200,
    target::Target = Near(0.0 + 0im),
    tolerance::Float64 = sqrt(eps(real(eltype(A)))),
    T::Type = ComplexF64,
    verbosity::Number = 0,
    verbose::Bool = false
)
    # `verbose = true` overrides verbosity only when verbosity is not set
    verbosity = verbosity == 0 && verbose ? 2 : verbosity

    solver_reltol = one(real(T))
    residuals::Vector{real(T)} = []

    n = size(A, 1)

    # V's vectors span the search space
    # AV will store (A - tI) * V, without any orthogonalization
    # W will hold AV, but with its columns orthogonal: AV = W * MA
    V = SubSpace(zeros(T, n, last(subspace_dimensions)))
    AV = SubSpace(zeros(T, n, last(subspace_dimensions)))
    W = SubSpace(zeros(T, n, last(subspace_dimensions)))

    # Temporaries (trying to reduces #allocations here)
    temporary = zeros(T, n, last(subspace_dimensions))

    # Projected matrices
    MA = ProjectedMatrix(zeros(T, last(subspace_dimensions), last(subspace_dimensions)))
    M = ProjectedMatrix(zeros(T, last(subspace_dimensions), last(subspace_dimensions)))

    # k is the number of converged eigenpairs
    k = 0

    # m is the current dimension of the search subspace
    m = 0

    pschur = PartialSchur(zeros(T, n, pairs), T)

    # Current eigenvalue
    λ = zero(T)

    # Current residual vector
    r = Vector{T}(undef, n)

    iter = 1

    harmonic_ritz_values::Vector{Vector{T}} = []
    converged_ritz_values::Vector{Vector{T}} = []

    local F::GeneralizedSchur
    local lastλ = T(NaN)
    local resnorm = real(T(Inf))

    progress = verbosity==1 ? Progress(max_iter) : nothing

    while k ≤ pairs && iter ≤ max_iter
        solver_reltol /= 2

        m += 1

        ### Expand
        resize!(V, m)
        resize!(AV, m)
        resize!(W, m)
        resize!(M, m)
        resize!(MA, m)

        if iter == 1
            rand!(V.curr)
        else
            solve_deflated_correction!(solver, A, V.curr, pschur.all, λ, r, solver_reltol)
        end

        # Search space is orthonormalized
        orthogonalize_and_normalize!(V.prev, V.curr, zeros(T, m - 1), DGKS)

        # AV is just the product (A - τI)V
        mul!(AV.curr, A, V.curr)
        axpy!(-target.τ, V.curr, AV.curr)

        # Expand W with (A - τI)V, and then orthogonalize
        copyto!(W.curr, AV.curr)

        # Orthogonalize w.r.t. the converged Schur vectors
        just_orthogonalize!(pschur.locked, W.curr, DGKS)

        # Orthonormalize W[:, m] w.r.t. previous columns of W
        MA.matrix[m, m] = orthogonalize_and_normalize!(W.prev, W.curr, view(MA.matrix, 1 : m - 1, m), DGKS)

        # Update the right-most column and last row of M = W' * V
        M.matrix[m, m] = dot(W.curr, V.curr)
        @inbounds for i = 1 : m - 1
            M.matrix[i, m] = dot(view(W.basis, :, i), V.curr)
            M.matrix[m, i] = dot(W.curr, view(V.basis, :, i))
        end

        # Assert orthogonality of V and W
        # Assert W * MA = (I - QQ') * (A - τI) * V
        # Assert that M = W' * V
        # @assert norm(W.all'W.all - I) < 1e-12
        # @assert norm(V.all'V.all - I) < 1e-12
        # @assert norm(M.curr - W.all'V.all) < 1e-12
        # @assert norm(W.all * MA.curr - AV.all + (pschur.locked * (pschur.locked'AV.all))) < pairs * tolerance

        search = true

        ### Extract
        while search

            F = schur(MA.curr, M.curr)
            λ = extract_harmonic!(F, V, AV, r, pschur, target.τ)
            resnorm = norm(r)

            # Convergence history of the harmonic Ritz values
            push!(harmonic_ritz_values, F.alpha ./ F.beta)
            push!(converged_ritz_values, copy(pschur.values))
            push!(residuals, resnorm)

            verbosity > 1 && println("Residual = ", resnorm)

            # A Ritz vector is converged
            if resnorm ≤ tolerance
                verbosity > 1 && println("Found an eigenvalue ", λ)
                lastλ = λ
                push!(converged_ritz_values[iter], λ)

                push!(pschur.values, λ)
                lock!(pschur)

                k += 1

                # Are we done yet?
                k == pairs && return pschur, harmonic_ritz_values, converged_ritz_values, residuals

                # Now remove this Schur vector from the search space.
                keep = 2 : m
                shrink!(temporary, V, view(F.right, :, keep), m - 1)
                shrink!(temporary, AV, view(F.right, :, keep), m - 1)
                shrink!(temporary, W, view(F.left, :, keep), m - 1)
                shrink!(M, view(F.T, keep, keep), m - 1)
                shrink!(MA, view(F.S, keep, keep), m - 1)

                m -= 1

                # TODO: Can the search space become empty? Probably, but not likely.
                m == 0 && throw(ErrorException("Search space became empty: TODO."))

                solver_reltol = one(real(T))
            else
                search = false
            end
        end

        ### Restart
        if m == last(subspace_dimensions)
            verbosity > 1 && println("Shrinking the search space.")

            # Move min_dimension of the smallest harmonic Ritz values up front
            keep = 1 : first(subspace_dimensions)
            schur_sort!(SmallestMagnitude(0.0+0im), F, keep)

            shrink!(temporary, V, view(F.right, :, keep), first(subspace_dimensions))
            shrink!(temporary, AV, view(F.right, :, keep), first(subspace_dimensions))
            shrink!(temporary, W, view(F.left, :, keep), first(subspace_dimensions))
            shrink!(M, view(F.T, keep, keep), first(subspace_dimensions))
            shrink!(MA, view(F.S, keep, keep), first(subspace_dimensions))
            m = first(subspace_dimensions)
        end

        iter += 1
        isnothing(progress) ||
            ProgressMeter.next!(progress,
                                showvalues=[(:Residual, resnorm),
                                            (:λ, lastλ),
                                            (:Pairs, "$(k)/$(pairs)")])
    end

    pschur, harmonic_ritz_values, converged_ritz_values, residuals
end

function extract_harmonic!(F::GeneralizedSchur, V::SubSpace, AV::SubSpace, r, pschur, τ)
    # Compute the Schur decomp to find the harmonic Ritz values

    # Move the smallest harmonic Ritz value up front
    schur_sort!(SmallestMagnitude(0.0 + 0im), F, 1)

    # Pre-ritz vector
    y = view(F.Z, :, 1)

    # Ritz vector
    mul!(pschur.active, V.all, y)

    # Rayleigh quotient λ = approx eigenvalue s.t. if r = (A-τI)u - λ * u, then r ⟂ u
    λ = conj(F.beta[1]) * F.alpha[1]

    # Residual r = (A - τI)u - λ * u = AV*y - λ * u
    mul!(r, AV.all, y)
    axpy!(-λ, pschur.active, r)

    # Orthogonalize w.r.t. Q
    just_orthogonalize!(pschur.locked, r, DGKS)

    # Assert that the residual is perpendicular to the Ritz vector
    # @assert abs(dot(r, pschur.active)) < 1e-5

    λ + τ
end

max_forks_repo_name: haampie/jacobi-davidson-julia


text:
# ---
# title: 1129. Shortest Path with Alternating Colors
# id: problem1129
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Breadth-first Search, Graph
# link: <https://leetcode.com/problems/shortest-path-with-alternating-colors/description/>
# hidden: true
# ---
#
# Consider a directed graph, with nodes labelled `0, 1, ..., n-1`.  In this
# graph, each edge is either red or blue, and there could be self-edges or
# parallel edges.
#
# Each `[i, j]` in `red_edges` denotes a red directed edge from node `i` to node
# `j`.  Similarly, each `[i, j]` in `blue_edges` denotes a blue directed edge
# from node `i` to node `j`.
#
# Return an array `answer` of length `n`, where each `answer[X]` is the length
# of the shortest path from node `0` to node `X` such that the edge colors
# alternate along the path (or `-1` if such a path doesn't exist).
#
#
#
# **Example 1:**
#
#
#
#     Input: n = 3, red_edges = [[0,1],[1,2]], blue_edges = []
#     Output: [0,1,-1]
#
#
# **Example 2:**
#
#
#
#     Input: n = 3, red_edges = [[0,1]], blue_edges = [[2,1]]
#     Output: [0,1,-1]
#
#
# **Example 3:**
#
#
#
#     Input: n = 3, red_edges = [[1,0]], blue_edges = [[2,1]]
#     Output: [0,-1,-1]
#
#
# **Example 4:**
#
#
#
#     Input: n = 3, red_edges = [[0,1]], blue_edges = [[1,2]]
#     Output: [0,1,2]
#
#
# **Example 5:**
#
#
#
#     Input: n = 3, red_edges = [[0,1],[0,2]], blue_edges = [[1,0]]
#     Output: [0,1,1]
#
#
#
#
# **Constraints:**
#
#   * `1 <= n <= 100`
#   * `red_edges.length <= 400`
#   * `blue_edges.length <= 400`
#   * `red_edges[i].length == blue_edges[i].length == 2`
#   * `0 <= red_edges[i][j], blue_edges[i][j] < n`
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
###
# Generates a histogram for SAA, Alpha LOO, and Alpha OR
# Used for illustrative picture in presentation
###
using Distributions, Random, DelimitedFiles
include("../src/JS_SAA_main.jl")

# Generates data that can be used to create a histogram of costs
# ps are half uniform, and half slightly concentrated at 1/d
# all data is synthetic
function runTest(numRuns, K, outPath; usePoisson=true)
	s = .95
	d = 10
	N = 20
	Random.seed!(8675309)
	p0 = ones(d)/d
	ps = rand(Dirichlet(ones(d)), floor(Int, K/2))
	ps = [ps rand(Dirichlet(3 * ones(d)), K - floor(Int, K/2)) ]
	@assert size(ps, 2) == K
	@assert size(ps, 1) == d
	alpha_grid = range(0, stop=50, length=75)
	lams = ones(K)

	cs, xs = JS.genNewsvendorsDiffSupp(repeat(1:d, inner=(1, K)), s, K)

	f = open("$(outPath)_$(K)_$(numRuns).csv", "w")

	#Header
	writedlm(f, ["Run" "Method" "TruePerf" "time" "alpha"], ',')

	#Compute the full-info value once for reference
	t =
	  @elapsed full_info = JS.zstar(xs, cs, ps, lams)
	println(full_info)
	writedlm(f, [1 "FullInfo" full_info t 0.], ',')

	Nhats = fill(N, K)
	for iRun = 1:numRuns
		if usePoisson
			Nhats = rand(Poisson(N), K)
		end
		mhats = JS.sim_path(ps, Nhats);

		#for data-driven shrinkage anchor
		phat_avg = vec(mean(mhats ./ Nhats', dims=2))

		#Gen the SAA const
		t =
		  @elapsed perf_SAA = JS.zbar(xs, cs, p0, 0., mhats, ps, lams)
		writedlm(f, [iRun "SAA" perf_SAA t 0.0], ',')

		#Gen the Oracle const
		t =
		  @elapsed alphaOR, min_indx, or_alpha_curve = JS.oracle_alpha(xs, cs, mhats, ps, lams, p0, alpha_grid)
		writedlm(f, [iRun "Oracle" or_alpha_curve[min_indx] t alphaOR], ',')

		#Gen the LOO cost with 1/d shrinkage
		t =
		  @elapsed alphaLOO, min_indx, looUnsc_curve = JS.loo_alpha(xs, cs, mhats, p0, alpha_grid)
		writedlm(f, [iRun "LOO_unif" or_alpha_curve[min_indx] t alphaLOO], ',')

		#Gen the LOO cost with the phatAvg shrinkage
		t =
		  @elapsed alphaLOO, min_indx, looUnsc_curve = JS.loo_alpha(xs, cs, mhats, phat_avg, alpha_grid)
		writedlm(f, [iRun "LOO_avg" or_alpha_curve[min_indx] t alphaLOO], ',')

		flush(f)
	end #endRun
	close(f)
end
runTest(10, 10, "../Results/tempFile", usePoisson=false)

#Run once with K = 1 for initial graph in presentation
runTest(1000, 1, "../Results/singleKUnifNewsvendor", usePoisson=false)

#Run again with K = 1000 for comparison graph in presentation
runTest(1000, 1000, "../Results/singleKUnifNewsvendor", usePoisson=false)

max_stars_repo_name: CausalML/DataPooling

max_issues_repo_name: CausalML/DataPooling

max_forks_repo_name: CausalML/DataPooling


text:
using Pkg

packages_to_test = [
    ("https://github.com/anthofflab/MimiRICE2010.jl.git", "mimi-next", "MimiRICE2010"),
    ("https://github.com/fund-model/MimiFUND.jl.git", "mimi-next", "MimiFUND")
]

for (pkg_url, pkg_rev, pkg_name) in packages_to_test
    mktempdir() do folder
        Pkg.activate(folder)

        Pkg.develop(PackageSpec(path=joinpath(@__DIR__, "..", "..")))

        Pkg.add(PackageSpec(url=pkg_url, rev=pkg_rev))

        Pkg.test(pkg_name)
    end
end


text:
using JSON
using JuMP
#using GLPKMathProgInterface
using Gurobi
using Chen1995
using Kurpel2018
using Check3DPackings

#solver = GLPKSolverMIP()
solver = GurobiSolver()

for mod in [Kurpel2018, Chen1995]
  @testset "$mod" begin
    inst_dir = "../../../instances/hbd_basic_tests/"
    #test_instances = filter(x -> endswith(x, ".json"), readdir(inst_dir))
    test_instances = ["test3.json"]
    @testset "$f" for f in test_instances
      json = JSON.parse(open(x -> read(x, String), inst_dir * f))
      pqr = map(l -> json[l], split("pqr", ""))
      LWH = map(l -> json[l], split("LWH", ""))
      model = mod.build(Model(solver = solver), pqr..., LWH...)
      #writeLP(model, "./$f.lp", genericnames = true)
      status = solve(model, suppress_warnings = true)
      @test string(status) == json["expected_status"]
      if status == :Optimal
        d, xyz, pqr′ = mod.extract_solution(model, pqr)
        @test !has_violations(d, xyz..., pqr′..., LWH...)
        @test all(is_orientation_of.(pqr..., pqr′...))
      end
    end
  end
end


text:
#function DemoMNIST()
    # Training a deep autoencoder on MNIST
    # The method uses Nesterov's accelerated gradient, with minibatches
    # (c) David Barber, University College London 2015

    #TODO: GPU version

    StartCode()
    bias=ADvariable()
h=ADvariable()
W=ADvariable()
        h2=rectlinAXplusBias(W,h,bias)
    loss=sumSquare(h2)
    net=EndCode()

    #instantiate root node values:
net.value[bias]=rand(2,1);
net.value[W]=rand(2,2);
net.value[h]=rand(2,3)

    net=compile(net) # compile the DAG and preallocate memory

@gpu CUDArt.init([0])
@gpu net=convert(net,"GPU")
gradcheck(net)
    println("Training:")
    error=Array(Float64,0)
    ParsToUpdate=Parameters(net)
    velo=NesterovInit(net)
    minibatchstart=1 # starting datapoint for the minibatch
    for iter=1:TrainingIts
        LearningRate=0.5/(1+iter/500)
        minibatchstart,minibatch=GetBatch(minibatchstart,BatchSize,Ntrain)
#        net.value[ytrain]=data[:,minibatch] # select batch of data
# @gpu       net.value[ytrain]=CudaArray(data[:,minibatch]) # select batch of data
        ADforward!(net;exclude=[ypred]) # can ignore ypred during training
        ADbackward!(net)
        push!(error,extract(net.value[net.FunctionNode])[1])
        println("iteration $iter: training loss = $(error[iter])")
        for par in ParsToUpdate
            #GradientDescentUpdate!(net.value[par],net.gradient[par],LearningRate)
            NesterovGradientDescentUpdate!(net.value[par],net.gradient[par],velo[par],LearningRate,iter)
        end
    end

    # Get the predictions:
    ADforward!(net;exclude=[])
    plot(error); title("training loss")
    figure()
    for i=1:20 # plot the reconstructions for a few datapoints
        p=imshow([reshape(net.value[h[1]][:,i],28,28)'  reshape(net.value[ypred.index][:,i],28,28)'],interpolation="none",cmap=ColorMap("gray"))
        display(p)
        println("press key to continue")
        readline(STDIN)
    end

max_stars_repo_name: davidbarber/AutoDiff

max_issues_repo_name: davidbarber/AutoDiff

max_forks_repo_name: davidbarber/AutoDiff


text:
include("gear.jl")
include("minlp.jl")
include("nlp1.jl")
include("nlp2.jl")
include("pool1.jl")
include("bad_expressions.jl")

include("moi_test.jl")


text:
using Parameters: @with_kw

calc_direction(_u, J, _L, Q) = det(vcat(J, rawtangentmat(Q)))
# rawtangentmat(Q) == Q[end:end, 1] --- tangent w/o det-fixing

find_simple_bifurcation!(cache, opts, sbint) =
    find_simple_bifurcation!(cache, opts, sbint.u0, sbint.u1,
                             sbint.direction)

find_simple_bifurcation!(cache, opts, args...) =
    find_zero!(cache, opts, calc_direction, args...)

@with_kw struct FindZeroInputError <: Exception
    f0
    f1
    u0
    u1
    f
end

Base.showerror(io::IO, e::FindZeroInputError) = print(io, """
Test function `f` evaluated at `u0` and `u1` have the same sign.
f0 = $(e.f0)
f1 = $(e.f1)
u0 = $(e.u0)
u1 = $(e.u1)
f = $(e.f)
""")

@with_kw struct ZeroNotFoundError <: Exception
    msg = "Zero not found"
    u::AbstractVector
    f::Real
    h::Real
end

function Base.showerror(io::IO, e::ZeroNotFoundError)
    println(io, e.msg)
    println(io, "h = ", e.h)
    println(io, "f = ", e.f)
    println(io, "u =")
    show(IOContext(io,
                   :displaysize => displaysize(io),
                   :limit => true),
         MIME("text/plain"),
         e.u')
    println(io)
end

function find_zero!(cache, opts, f, u0, u1, direction)
    prob_cache = cache.prob_cache
    H = cache.H
    J = cache.J
    Q = cache.Q
    rtol = opts.rtol
    atol = opts.atol

    H, J = residual_jacobian!(H, J, u1, prob_cache)
    A = vcat(J, _zeros(J, 1, size(J, 2)))  # TODO: improve
    L, Q = _lq!(A)
    f1 = f(u1, J, L, Q)

    H, J = residual_jacobian!(H, J, u0, prob_cache)
    A = vcat(J, _zeros(J, 1, size(J, 2)))  # TODO: improve
    L, Q = _lq!(A)
    tJ = tangent(L, Q)
    f0 = f(u0, J, L, Q)

    if f0 * f1 > 0
        throw(FindZeroInputError(f0, f1, u0, u1, f))
    end

    fu = f0
    u = u0
    h = norm(u0 .- u1) / 2
    for _ in 1:opts.max_adaptations
        # predictor
        v = u .+ direction * h .* tJ

        # corrector
        for _ in 1:opts.max_corrector_steps
            w, _, H, L, Q, J = corrector_step!(H, J, Q, v, prob_cache)
            if isalmostzero(H, atol)
                cache.corrector_success = true
                @goto corrector_success
            end
            v = w
        end
        throw(ZeroNotFoundError(msg="Corrector failed", u=v, f=fu, h=h))

        @label corrector_success
        tJv = tangent(L, Q)
        fv = f(v, J, L, Q)
        @assert all(isfinite.(v))

        # secant
        h = - fv / (fv - fu) * h
        u = v
        fu = fv

        if abs(h) < opts.h_zero
            return v, tJv, L, Q, J
        end

        if tJ ⋅ tJv < 0
            direction *= -1
        end
        tJ = tJv
    end
    throw(ZeroNotFoundError(u=u, f=fu, h=h))
end


text:
import .CUDA

CUDA.allowscalar(false)

function refresh(
    rng::Union{AbstractRNG, AbstractVector{<:AbstractRNG}},
    ::FullMomentumRefreshment,
    h::Hamiltonian,
    z::PhasePoint{TA}
) where {T<:AbstractFloat,TA<:CUDA.CuArray{<:T}}
    r = CUDA.CuArray{T, 2}(undef, size(h.metric)...)
    CUDA.CURAND.randn!(r)
    return phasepoint(h, z.θ, r)
end

# TODO: Ideally this should be merged with the CPU version. The function is
#       essentially the same but sampling requires a custom call to CUDA.
function mh_accept_ratio(
    rng::Union{AbstractRNG, AbstractVector{<:AbstractRNG}}, Horiginal::TA, Hproposal::TA,
) where {T<:AbstractFloat, TA<:CUDA.CuArray{<:T}}
    α = min.(one(T), exp.(Horiginal .- Hproposal))
    # NOTE: There is a chance that sharing the RNG over multiple
    #       chains for accepting / rejecting might couple
    #       the chains. We need to revisit this more rigirously
    #       in the future. See discussions at
    #       https://github.com/TuringLang/AdvancedHMC.jl/pull/166#pullrequestreview-367216534
    r = CUDA.CuArray{T, 1}(undef, length(Horiginal))
    CUDA.CURAND.rand!(r)
    accept = r .< α
    return accept, α
end


text:
module MatrixLMnet

using MatrixLM
import MatrixLM.calc_preds, MatrixLM.calc_preds!,
    MatrixLM.calc_resid, MatrixLM.calc_resid!

using DataFrames
using Distributed
using LinearAlgebra
using LinearAlgebra.BLAS
import LinearAlgebra.mul!, LinearAlgebra.norm, LinearAlgebra.dot,
     LinearAlgebra.diagm, LinearAlgebra.eigmin, LinearAlgebra.eigmax
using Random
using Statistics
using MLBase

export Response, Predictors, RawData, get_X, get_Z, get_Y, contr,
    add_intercept, remove_intercept, shuffle_rows, shuffle_cols,
    cd!, cd_active!, ista!, fista!, fista_bt!, admm!,
    mlmnet, Mlmnet, coef, coef_2d, predict, fitted, resid,
    mlmnet_perms,
    make_folds, make_folds_conds, mlmnet_cv, Mlmnet_cv, mlmnet_cv,
    avg_mse, lambda_min, avg_prop_zero, mlmnet_cv_summary

# Helper functions
include("utilities/std_helpers.jl")
include("mlmnet/mlmnet_helpers.jl")

# Optimizing algorithms
include("methods/cd.jl")
include("methods/ista.jl")
include("methods/fista.jl")
include("methods/fista_bt.jl")
include("methods/admm.jl")

# Top level functions that call Elastic Net algorithms using warm starts
include("mlmnet/mlmnet.jl")

# Predictions and residuals
include("utilities/predict.jl")

# Permutations
include("mlmnet/mlmnet_perms.jl")

# Cross-validation
include("crossvalidation/mlmnet_cv_helpers.jl")
include("crossvalidation/mlmnet_cv.jl")
include("crossvalidation/mlmnet_cv_summary.jl")

end