text:
using JLD2, SparseArrays

#=
Assoc Serialized for saving using JLD2
Note that saving would convert row and col types to number.
=#

# Delimiter for saving- using new line is safer than comma!
del = "\n"
maxchar = typemax(UInt16)-5

SparseOrTranspose = Union{AbstractSparseMatrix,Adjoint{<:Any,<:SparseMatrixCSC},Transpose{<:Any,<:SparseMatrixCSC}}

struct AssocSerial
    rowstr::Array{AbstractString,1}
    colstr::Array{AbstractString,1}
    valstr::Array{AbstractString,1}

    rownum::Array{AbstractString,1}
    colnum::Array{AbstractString,1}
    valnum::Array{AbstractString,1}

    A::SparseOrTranspose
end

function writeas(data::Assoc)
    #Get parts from the assoc
    row = getrow(data)
    col = getcol(data)
    val = getval(data)
    A   = dropzeros!(getadj(data))

    #Split the mapping array into string and number arrays for storage
    rowindex = searchsortedfirst(row,AbstractString,lt=isa)
    colindex = searchsortedfirst(col,AbstractString,lt=isa)
    valindex = searchsortedfirst(val,AbstractString,lt=isa)
    #test if there are string elements in the array
    if (rowindex == 1)
        rowstr = []
    else
        rowstr = chunkstring(join(row[1:(rowindex-1)],del)*del,maxchar)
    end

    if (colindex == 1)
        colstr = []
    else
        colstr = chunkstring(join(col[1:(colindex-1)],del)*del,maxchar)
    end

    if (valindex == 1)
        valstr = []
    else
        valstr = chunkstring(join(val[1:(valindex-1)],del)*del,maxchar)
    end

    #test if there are number elements in the array.  Note that by serialization all numbers will be converted to float.

    if (rowindex == size(row,1)+1)
        rownum = []
    else
        rownum = chunkstring(join(row[rowindex:end],del)*del,maxchar)
    end

    if (colindex == size(col,1)+1)
        colnum = []
    else
        colnum = chunkstring(join(col[colindex:end],del)*del,maxchar)
    end

    if (valindex == size(val,1)+1)
        valnum = []
    else
        valnum = chunkstring(join(val[valindex:end],del)*del,maxchar)
    end

    #Reconstitute converted parts into Serialized

    return AssocSerial(rowstr,colstr,valstr,rownum,colnum,valnum,A)
end

function readas(serData::AssocSerial)
    row = []
    # Last character is delimiter- in case saved with different delimiter
    if !isempty(serData.rowstr)
       row = vcat(row, split(join(serData.rowstr,"")[1:end-1],serData.rowstr[end][end]))
    end
    if !isempty(serData.rownum)
        row = vcat(row, map(float,split(join(serData.rownum,"")[1:end-1],serData.rownum[end][end])))
    end
    row = Array{Union{AbstractString,Number}}(row)

    col = []
    if !isempty(serData.colstr)
       col = vcat(col, split(join(serData.colstr,"")[1:end-1],serData.colstr[end][end]))
    end
    if !isempty(serData.colnum)
        col = vcat(col, map(float,split(join(serData.colnum,"")[1:end-1],serData.colnum[end][end])))
    end
    col = Array{Union{AbstractString,Number}}(col)

    val = []
    if !isempty(serData.valstr)
       val = vcat(val, split(join(serData.valstr,"")[1:end-1],serData.valstr[end][end]))
    end
    if !isempty(serData.valnum)
        val = vcat(val, map(float,split(join(serData.valnum,"")[1:end-1],serData.valnum[end][end])))
    end
    val = Array{Union{AbstractString,Number}}(val)

    return Assoc(row,col,val,serData.A)
end

function chunkstring(str,chunksize)
    res = []

    idx = 1:chunksize:length(str)

    for i in idx
        res = [res; str[i:min(length(str),i+chunksize-1)]]
    end

    return res
end

function saveassoc(savedir,A::Assoc)
    Awrite = writeas(A)
    @save savedir Awrite
end

function loadassoc(savedir)
    @load savedir Awrite
    return readas(Awrite)
end


text:
#VERSION >= v"0.4.0-dev+6521" && __precompile__()

"""
Julia wrappers around Let's be rational black76 calculator from Peter Jäckel
"""
module Black76

export normcdf, normpdf, px, impliedvol, ntv

const libpath = joinpath(dirname(@__FILE__), "..", "deps", "src", "liblbr")
"Cumulative standardized normal density"
normcdf(z::Real) = ccall((:norm_cdf, libpath), Float64, (Float64,), z)
"standard normal density"
normpdf(z::Real) = ccall((:norm_pdf, libpath), Float64, (Float64,), z)

"normalised black76 option time value"
function ntv(x::Real, w::Real)
    isnan(x) && return x
    isfinite(x) && return ccall((:normalised_black_call, libpath), Float64, (Float64, Float64), -abs(x), w)
    return 0
end

"Black76 European option value"
function px(fwd::Real, strike::Real, yf::Real, q::Real, iv::Real)
    ccall((:black, libpath), Float64, (Float64, Float64, Float64, Float64, Float64),
          fwd, strike, iv, yf, q)
end

"Black76 European option implied volatility"
function impliedvol(fwd::Real, strike::Real, yf::Real, q::Real, px::Real)
    ccall((:implied_volatility_from_a_transformed_rational_guess, libpath), Float64,
          (Float64, Float64, Float64, Float64, Float64), px, fwd, strike, yf, q)
end

end # module Black76

ext: jl


text:
"""
    nll(L::LPDO{MPS}, data::Matrix{Pair{String,Int}})

Compute the negative log-likelihood using an MPS ansatz
over a dataset `data`:

`nll ∝ -∑ᵢlog P(σᵢ)`
"""
function nll(L::LPDO{MPS}, data::Matrix{Pair{String,Int}})
  data = convertdatapoints(copy(data); state=true)
  ψ = L.X
  N = length(ψ)
  @assert N == size(data)[2]
  loss = 0.0
  s = siteinds(ψ)

  for n in 1:size(data)[1]
    x = data[n, :]
    ψx = dag(ψ[1]) * state(x[1], s[1])
    for j in 2:N
      ψ_r = dag(ψ[j]) * state(x[j], s[j])
      ψx = ψx * ψ_r
    end
    prob = abs2(ψx[])
    loss -= log(prob) / size(data)[1]
  end
  return loss
end

nll(ψ::MPS, data::Matrix{Pair{String,Int}}) = nll(LPDO(ψ), data)

"""
    nll(L::LPDO{MPO}, data::Matrix{Pair{String,Int}})

Compute the negative log-likelihood using an LPDO ansatz
over a dataset `data`:

`nll ∝ -∑ᵢlog P(σᵢ)`
"""
function nll(L::LPDO{MPO}, data::Matrix{Pair{String,Int}})
  data = convertdatapoints(copy(data); state=true)
  lpdo = L.X
  N = length(lpdo)
  loss = 0.0
  s = firstsiteinds(lpdo)
  for n in 1:size(data)[1]
    x = data[n, :]

    # Project LPDO into the measurement eigenstates
    Φdag = dag(copy(lpdo))
    for j in 1:N
      Φdag[j] = Φdag[j] = Φdag[j] * state(x[j], s[j])
    end

    # Compute overlap
    prob = inner(Φdag, Φdag)
    loss -= log(real(prob)) / size(data)[1]
  end
  return loss
end

"""
    PastaQ.gradlogZ(L::LPDO; sqrt_localnorms = nothing)
    PastaQ.gradlogZ(ψ::MPS; localnorms = nothing)

Compute the gradients of the log-normalization with respect
to each LPDO tensor component:

- `∇ᵢ = ∂ᵢlog⟨ψ|ψ⟩` for `ψ = M = MPS`
- `∇ᵢ = ∂ᵢlogTr(ρ)` for `ρ = M M†` , `ρ = LPDO`
"""
function gradlogZ(lpdo::LPDO; sqrt_localnorms=nothing)
  M = lpdo.X
  N = length(M)
  L = Vector{ITensor}(undef, N - 1)
  R = Vector{ITensor}(undef, N)

  if isnothing(sqrt_localnorms)
    sqrt_localnorms = ones(N)
  end

  # Sweep right to get L
  L[1] = dag(M[1]) * prime(M[1], "Link")
  for j in 2:(N - 1)
    L[j] = L[j - 1] * dag(M[j])
    L[j] = L[j] * prime(M[j], "Link")
  end
  Z = L[N - 1] * dag(M[N])
  Z = real((Z * prime(M[N], "Link"))[])
  # Sweep left to get R
  R[N] = dag(M[N]) * prime(M[N], "Link")
  for j in reverse(2:(N - 1))
    R[j] = R[j + 1] * dag(M[j])
    R[j] = R[j] * prime(M[j], "Link")
  end
  # Get the gradients of the normalization
  gradients = Vector{ITensor}(undef, N)
  gradients[1] = prime(M[1], "Link") * R[2] / (sqrt_localnorms[1] * Z)
  for j in 2:(N - 1)
    gradients[j] = (L[j - 1] * prime(M[j], "Link") * R[j + 1]) / (sqrt_localnorms[j] * Z)
  end
  gradients[N] = (L[N - 1] * prime(M[N], "Link")) / (sqrt_localnorms[N] * Z)

  return 2 * gradients, log(Z)
end

gradlogZ(ψ::MPS; localnorms=nothing) = gradlogZ(LPDO(ψ); sqrt_localnorms=localnorms)

"""
    PastaQ.gradnll(L::LPDO{MPS}, data::Array; sqrt_localnorms = nothing)
    PastaQ.gradnll(ψ::MPS, data::Array; localnorms = nothing)

Compute the gradients of the cross-entropy between the MPS probability
distribution of the empirical data distribution for a set of projective
measurements in different local bases. The probability of a single
data-point `σ = (σ₁,σ₂,…)` is :

`P(σ) = |⟨σ|Û|ψ⟩|²`

where `Û` is the depth-1 local circuit implementing the basis rotation.
The cross entropy function is

`nll ∝ -∑ᵢlog P(σᵢ)`

where `∑ᵢ` runs over the measurement data. Returns the gradients:

`∇ᵢ = - ∂ᵢ⟨log P(σ))⟩_data`
"""
function gradnll(L::LPDO{MPS}, data::Matrix{Pair{String,Int}}; sqrt_localnorms=nothing)
  data = convertdatapoints(copy(data); state=true)

  ψ = L.X
  N = length(ψ)

  s = siteinds(ψ)

  links = [linkind(ψ, n) for n in 1:(N - 1)]

  ElT = eltype(ψ[1])

  nthreads = Threads.nthreads()

  L = [Vector{ITensor}(undef, N) for _ in 1:nthreads]
  Lpsi = [Vector{ITensor}(undef, N) for _ in 1:nthreads]

  R = [Vector{ITensor}(undef, N) for _ in 1:nthreads]
  Rpsi = [Vector{ITensor}(undef, N) for _ in 1:nthreads]

  for nthread in 1:nthreads
    for n in 1:(N - 1)
      L[nthread][n] = ITensor(ElT, undef, links[n])
      Lpsi[nthread][n] = ITensor(ElT, undef, s[n], links[n])
    end
    Lpsi[nthread][N] = ITensor(ElT, undef, s[N])

    for n in N:-1:2
      R[nthread][n] = ITensor(ElT, undef, links[n - 1])
      Rpsi[nthread][n] = ITensor(ElT, undef, links[n - 1], s[n])
    end
    Rpsi[nthread][1] = ITensor(ElT, undef, s[1])
  end

  if isnothing(sqrt_localnorms)
    sqrt_localnorms = ones(N)
  end

  ψdag = dag(ψ)

  gradients = [[ITensor(ElT, inds(ψ[j])) for j in 1:N] for _ in 1:nthreads]

  grads = [[ITensor(ElT, undef, inds(ψ[j])) for j in 1:N] for _ in 1:nthreads]

  loss = zeros(nthreads)

  Threads.@threads for n in 1:size(data)[1]
    nthread = Threads.threadid()

    x = data[n, :]

    """ LEFT ENVIRONMENTS """
    L[nthread][1] .= ψdag[1] .* state(x[1], s[1])
    for j in 2:(N - 1)
      Lpsi[nthread][j] .= L[nthread][j - 1] .* ψdag[j]
      L[nthread][j] .= Lpsi[nthread][j] .* state(x[j], s[j])
    end
    Lpsi[nthread][N] .= L[nthread][N - 1] .* ψdag[N]
    ψx = (Lpsi[nthread][N] * state(x[N], s[N]))[]
    prob = abs2(ψx)
    loss[nthread] -= log(prob) / size(data)[1]

    """ RIGHT ENVIRONMENTS """
    R[nthread][N] .= ψdag[N] .* state(x[N], s[N])
    for j in reverse(2:(N - 1))
      Rpsi[nthread][j] .= ψdag[j] .* R[nthread][j + 1]
      R[nthread][j] .= Rpsi[nthread][j] .* state(x[j], s[j])
    end

    """ GRADIENTS """
    # TODO: fuse into one call to mul!
    grads[nthread][1] .= state(x[1], s[1]) .* R[nthread][2]
    gradients[nthread][1] .+= (1 / (sqrt_localnorms[1] * ψx)) .* grads[nthread][1]
    for j in 2:(N - 1)
      Rpsi[nthread][j] .= L[nthread][j - 1] .* state(x[j], s[j])
      # TODO: fuse into one call to mul!
      grads[nthread][j] .= Rpsi[nthread][j] .* R[nthread][j + 1]
      gradients[nthread][j] .+= (1 / (sqrt_localnorms[j] * ψx)) .* grads[nthread][j]
    end
    grads[nthread][N] .= L[nthread][N - 1] .* state(x[N], s[N])
    gradients[nthread][N] .+= (1 / (sqrt_localnorms[N] * ψx)) .* grads[nthread][N]
  end

  for nthread in 1:nthreads
    for g in gradients[nthread]
      g .= (-2 / size(data)[1]) .* g
    end
  end

  gradients_tot = [ITensor(ElT, inds(ψ[j])) for j in 1:N]
  loss_tot = 0.0
  for nthread in 1:nthreads
    gradients_tot .+= gradients[nthread]
    loss_tot += loss[nthread]
  end

  return gradients_tot, loss_tot
end

function gradnll(ψ::MPS, data::Matrix{Pair{String,Int}}; localnorms=nothing)
  return gradnll(LPDO(ψ), data; sqrt_localnorms=localnorms)
end

"""
    PastaQ.gradnll(lpdo::LPDO{MPO}, data::Array; sqrt_localnorms = nothing)

Compute the gradients of the cross-entropy between the LPDO probability
distribution of the empirical data distribution for a set of projective
measurements in different local bases. The probability of a single
data-point `σ = (σ₁,σ₂,…)` is :

`P(σ) = ⟨σ|Û ρ Û†|σ⟩ = |⟨σ|Û M M† Û†|σ⟩ = |⟨σ|Û M`

where `Û` is the depth-1 local circuit implementing the basis rotation.
The cross entropy function is

`nll ∝ -∑ᵢlog P(σᵢ)`

where `∑ᵢ` runs over the measurement data. Returns the gradients:

`∇ᵢ = - ∂ᵢ⟨log P(σ))⟩_data`
"""
function gradnll(L::LPDO{MPO}, data::Matrix{Pair{String,Int}}; sqrt_localnorms=nothing)
  data = convertdatapoints(copy(data); state=true)
  lpdo = L.X
  N = length(lpdo)

  s = firstsiteinds(lpdo)

  links = [linkind(lpdo, n) for n in 1:(N - 1)]

  kraus = Index[]
  for j in 1:N
    push!(kraus, firstind(lpdo[j], "Purifier"))
  end

  ElT = eltype(lpdo[1])

  nthreads = Threads.nthreads()

  L = [Vector{ITensor}(undef, N) for _ in 1:nthreads]
  Llpdo = [Vector{ITensor}(undef, N) for _ in 1:nthreads]
  Lgrad = [Vector{ITensor}(undef, N) for _ in 1:nthreads]

  R = [Vector{ITensor}(undef, N) for _ in 1:nthreads]
  Rlpdo = [Vector{ITensor}(undef, N) for _ in 1:nthreads]

  Agrad = [Vector{ITensor}(undef, N) for _ in 1:nthreads]

  T = [Vector{ITensor}(undef, N) for _ in 1:nthreads]
  Tp = [Vector{ITensor}(undef, N) for _ in 1:nthreads]

  grads = [Vector{ITensor}(undef, N) for _ in 1:nthreads]
  gradients = [Vector{ITensor}(undef, N) for _ in 1:nthreads]

  for nthread in 1:nthreads
    for n in 1:(N - 1)
      L[nthread][n] = ITensor(ElT, undef, links[n]', links[n])
    end
    for n in 2:(N - 1)
      Llpdo[nthread][n] = ITensor(ElT, undef, kraus[n], links[n]', links[n - 1])
    end
    for n in 1:(N - 2)
      Lgrad[nthread][n] = ITensor(ElT, undef, links[n], kraus[n + 1], links[n + 1]')
    end
    Lgrad[nthread][N - 1] = ITensor(ElT, undef, links[N - 1], kraus[N])

    for n in N:-1:2
      R[nthread][n] = ITensor(ElT, undef, links[n - 1]', links[n - 1])
    end
    for n in (N - 1):-1:2
      Rlpdo[nthread][n] = ITensor(ElT, undef, links[n - 1]', kraus[n], links[n])
    end

    Agrad[nthread][1] = ITensor(ElT, undef, kraus[1], links[1]', s[1])
    for n in 2:(N - 1)
      Agrad[nthread][n] = ITensor(ElT, undef, links[n - 1], kraus[n], links[n]', s[n])
    end

    T[nthread][1] = ITensor(ElT, undef, kraus[1], links[1])
    Tp[nthread][1] = prime(T[nthread][1], "Link")
    for n in 2:(N - 1)
      T[nthread][n] = ITensor(ElT, undef, kraus[n], links[n], links[n - 1])
      Tp[nthread][n] = prime(T[nthread][n], "Link")
    end
    T[nthread][N] = ITensor(ElT, undef, kraus[N], links[N - 1])
    Tp[nthread][N] = prime(T[nthread][N], "Link")

    grads[nthread][1] = ITensor(ElT, undef, links[1], kraus[1], s[1])
    gradients[nthread][1] = ITensor(ElT, links[1], kraus[1], s[1])
    for n in 2:(N - 1)
      grads[nthread][n] = ITensor(ElT, undef, links[n], links[n - 1], kraus[n], s[n])
      gradients[nthread][n] = ITensor(ElT, links[n], links[n - 1], kraus[n], s[n])
    end
    grads[nthread][N] = ITensor(ElT, undef, links[N - 1], kraus[N], s[N])
    gradients[nthread][N] = ITensor(ElT, links[N - 1], kraus[N], s[N])
  end

  if isnothing(sqrt_localnorms)
    sqrt_localnorms = ones(N)
  end

  loss = zeros(nthreads)

  Threads.@threads for n in 1:size(data)[1]
    nthread = Threads.threadid()

    x = data[n, :]

    """ LEFT ENVIRONMENTS """
    T[nthread][1] .= lpdo[1] .* dag(state(x[1], s[1]))
    L[nthread][1] .= prime(T[nthread][1], "Link") .* dag(T[nthread][1])
    for j in 2:(N - 1)
      T[nthread][j] .= lpdo[j] .* dag(state(x[j], s[j]))
      Llpdo[nthread][j] .= prime(T[nthread][j], "Link") .* L[nthread][j - 1]
      L[nthread][j] .= Llpdo[nthread][j] .* dag(T[nthread][j])
    end
    T[nthread][N] .= lpdo[N] .* dag(state(x[N], s[N]))
    prob = L[nthread][N - 1] * prime(T[nthread][N], "Link")
    prob = prob * dag(T[nthread][N])
    prob = real(prob[])
    loss[nthread] -= log(prob) / size(data)[1]

    """ RIGHT ENVIRONMENTS """
    R[nthread][N] .= prime(T[nthread][N], "Link") .* dag(T[nthread][N])
    for j in reverse(2:(N - 1))
      Rlpdo[nthread][j] .= prime(T[nthread][j], "Link") .* R[nthread][j + 1]
      R[nthread][j] .= Rlpdo[nthread][j] .* dag(T[nthread][j])
    end

    """ GRADIENTS """
    Tp[nthread][1] .= prime(lpdo[1], "Link") .* dag(state(x[1], s[1]))
    Agrad[nthread][1] .= Tp[nthread][1] .* state(x[1], s[1])
    grads[nthread][1] .= R[nthread][2] .* Agrad[nthread][1]
    gradients[nthread][1] .+= (1 / (sqrt_localnorms[1] * prob)) .* grads[nthread][1]
    for j in 2:(N - 1)
      Tp[nthread][j] .= prime(lpdo[j], "Link") .* dag(state(x[j], s[j]))
      Lgrad[nthread][j - 1] .= L[nthread][j - 1] .* Tp[nthread][j]
      Agrad[nthread][j] .= Lgrad[nthread][j - 1] .* state(x[j], s[j])
      grads[nthread][j] .= R[nthread][j + 1] .* Agrad[nthread][j]
      gradients[nthread][j] .+= (1 / (sqrt_localnorms[j] * prob)) .* grads[nthread][j]
    end
    Tp[nthread][N] .= prime(lpdo[N], "Link") .* dag(state(x[N], s[N]))
    Lgrad[nthread][N - 1] .= L[nthread][N - 1] .* Tp[nthread][N]
    grads[nthread][N] .= Lgrad[nthread][N - 1] .* state(x[N], s[N])
    gradients[nthread][N] .+= (1 / (sqrt_localnorms[N] * prob)) .* grads[nthread][N]
  end

  for nthread in 1:nthreads
    for g in gradients[nthread]
      g .= (-2 / size(data)[1]) .* g
    end
  end

  gradients_tot = Vector{ITensor}(undef, N)
  gradients_tot[1] = ITensor(ElT, links[1], kraus[1], s[1])
  for n in 2:(N - 1)
    gradients_tot[n] = ITensor(ElT, links[n], links[n - 1], kraus[n], s[n])
  end
  gradients_tot[N] = ITensor(ElT, links[N - 1], kraus[N], s[N])

  loss_tot = 0.0
  for nthread in 1:nthreads
    gradients_tot .+= gradients[nthread]
    loss_tot += loss[nthread]
  end

  return gradients_tot, loss_tot
end

"""
    PastaQ.gradients(L::LPDO, data::Array; sqrt_localnorms = nothing)
    PastaQ.gradients(ψ::MPS, data::Array; localnorms = nothing)

Compute the gradients of the cost function:
`C = log(Z) - ⟨log P(σ)⟩_data`
"""
function gradients(L::LPDO, data::Matrix{Pair{String,Int}}; sqrt_localnorms=nothing, kwargs...)
  g_logZ, logZ = gradlogZ(L; sqrt_localnorms=sqrt_localnorms)
  g_nll, nll = gradnll(L, data; sqrt_localnorms=sqrt_localnorms)

  grads = g_logZ + g_nll
  loss = logZ + nll

  # TODO: check if this slows down the training
  # permute the gradients
  for j in 1:length(L)
    grads[j] = permute(grads[j], inds(L.X[j])...)
  end

  return grads, loss
end

gradients(ψ::MPS, data::Matrix{Pair{String,Int}}; localnorms=nothing) =
  gradients(LPDO(ψ), data; sqrt_localnorms=localnorms)


text:
# An order of magnitude faster than parse
function buf_to_uint(buf::Array{UInt8,1}, L::T) where T<:Integer
  o = one(UInt64)
  n = zero(UInt64)
  p = UInt64(L-o)
  i = zero(T)
  t = UInt64(10)
  while i < L
    i += one(T)
    n += UInt64(getindex(buf, i)-0x30)*t^p
    p -= o
  end
  return n
end

function buf_to_int(buf::Array{UInt8,1}, L::Int64)
  c = 10^(L-1)
  n = zero(Int64)
  for i = 1:L
    n += c*(buf[i]-0x30)
    c = div(c, 10)
  end
  return n
end

function buf_to_int(buf::Array{UInt8,1}, i::Int64, j::Int64)
  c = Int16(10)^Int16(i-j)
  n = zero(Int16)
  for p = j:i
    n += c * Int16(buf[p]-0x30)
    c = div(c, Int16(10))
  end
  return n
end

# Factor of two faster than parse(Float64, str) for a string
# Order of magnitude faster than parse(Float64, String[uu]) for a u8 array
function buf_to_double(buf::Array{UInt8,1}, L::Int64)
  read_state = 0x00
  mant_sgn = 1.0e0
  exp_sgn = 1.0e0

  j = one(Int64)
  zz = zero(Int64)
  z = zero(Int8)
  o = one(Int8)
  i = z
  imax = Int8(16)
  p0 = Int64(10)^Int64(imax)
  p = p0
  t = Int64(10)
  d = zz
  m = zz
  q = zz
  v = zz
  x = zz
  while j ≤ L
    c = buf[j]
    if read_state == 0x00
      if is_u8_digit(c) && (i < imax)
        p = div(p,t)
        i += o
        v += Int64(c-0x30)*p
        m = p
      elseif c == 0x2d # '-'
        mant_sgn = -1.0e0
      elseif c == 0x2b # '+'
        mant_sgn = 1.0e0
      elseif c == 0x2e # '.'
        # transition to state 1
        read_state = 0x01
        i = z
        p = p0
      elseif c in (0x45, 0x46, 0x65, 0x66) # 'E', 'F', 'e', 'f'
        # transition to state 2 (no decimal)
        read_state = 0x02
        i = z
        p = p0
      end
    elseif read_state == 0x01
      if is_u8_digit(c) && i < imax
        p = div(p,t)
        i += o
        d += Int64(c-0x30)*p
      elseif c in (0x45, 0x46, 0x65, 0x66)
        # transition to state 2
        read_state = 0x02
        i = z
        p = p0
      end
    else
      if is_u8_digit(c) && (i < imax)
        p = div(p,t)
        i += o
        x += Int64(c-0x30)*p
        q = p
      elseif c == 0x2d # '-'
        exp_sgn = -1.0e0
      elseif c == 0x2b # '+'
        exp_sgn = 1.0e0
      else
        break
      end
    end
    j += 1
  end
  v1 = Float64(div(v, m))
  v2 = Float64(d) / Float64(p0)
  if q > zz
    v3 = 10.0e0^(exp_sgn * Float64(div(x,q)))
  else
    v3 = one(Float64)
  end
  return mant_sgn * (v1 + v2) * v3
end


text:
# Description:
# Dimensions: 2

# The plot on the left shows the three-hump Camel function on its recommended input domain, and the plot on the right shows only a portion of this domain, to allow for easier viewing of the function's key characteristics. The function has three local minima.

# Input Domain:
# The function is usually evaluated on the square xi ∈ [-5, 5], for all i = 1, 2.

# Global Minimum:
# https://www.sfu.ca/~ssurjano/camel33.png

function threeHumpCamel(xx)

    x1 = xx[1];
    x2 = xx[2];

    term1 = 2 * x1^2;
    term2 = -1.05 * x1^4;
    term3 = x1^6 / 6;
    term4 = x1 * x2;
    term5 = x2^2;

    y = term1 + term2 + term3 + term4 + term5;

    return [y]
end


text:
# This file is a part of BAT.jl, licensed under the MIT License (MIT).

struct GradFunc{
    OF<:Function,
    UF<:Function,
    VS<:AbstractValueShape,
    GS<:AbstractValueShape,
    Alg<:DifferentiationAlgorithm
} <: GradientFunction
    _orig_f::OF
    _unshaped_f::UF
    _input_shape::VS
    _grad_shape::GS
    _diffalg::Alg
end

function GradFunc(f::Function, diffalg::DifferentiationAlgorithm)
    input_shape = varshape(f)
    unshaped_f = unshaped(f)
    grad_shape = gradient_shape(input_shape)
    GradFunc(f, unshaped_f, input_shape, grad_shape, diffalg)
end

function Base.show(io::IO, gf::GradFunc)
    print(io, Base.typename(typeof(gf)).name, "(")
    show(io, gf._orig_f)
    print(io, ", ")
    show(io, gf._diffalg)
    print(io, ")")
end

function Base.show(io::IO, M::MIME"text/plain", gf::GradFunc)
    print(io, Base.typename(typeof(gf)).name, "(")
    show(io, M, gf._orig_f)
    print(io, ", ")
    show(io, M, gf._diffalg)
    print(io, ")")
end

_fixup_var(shape::AbstractValueShape, v::Any) = v

_fixup_var(shape::NamedTupleShape, v::NamedTuple) = shape(unshaped(v, shape))

function (gf::GradFunc)(v::Any)
    input_shape = gf._input_shape
    v_shaped = _fixup_var(input_shape, v)
    v_unshaped = unshaped(v_shaped)

    grad_f_unshaped = similar(v_unshaped)

    value = unshaped_gradient!(grad_f_unshaped, gf._unshaped_f, v_unshaped, gf._diffalg)

    (value, gf._grad_shape(grad_f_unshaped))
end

function (gf::GradFunc)(::typeof(!), grad_f::Any, v::Any)
    input_shape = gf._input_shape
    v_shaped = _fixup_var(input_shape, v)
    v_unshaped = unshaped(v_shaped)

    if isnothing(grad_f)
        R(gf._unshaped_f(v_unshaped))
    else
        grad_f_unshaped = unshaped(grad_f, gf._grad_shape)
        unshaped_gradient!(grad_f_unshaped, gf._unshaped_f, v_unshaped, gf._diffalg)
    end
end

function valgradof(f::Function, algorithm::DifferentiationAlgorithm = vjp_algorithm(f))
    GradFunc(f, algorithm)
end


text:
const Operation = @NamedTuple{name::Symbol, type::DataType, broadcasted::Bool}

_nameof(x) = nameof(typeof(x))
_nameof(x::Function) = nameof(x)

"""
    Module

A data structure to store information to generate hardware for a circuit.
This structure can be manually modified if needed but typically [`@circuit`](@ref)
is used to auto-populate it.

Hardware generation traverses `dfg` and uses `handlers` to generate Verilog strings.

# Fields:
- `name::Symbol`: the name of the module
- `parameters::Dict{Symbol, Number}`: a map from the name of each parameter to its default value
- `submodules::Vector{Type}`: a list of submodule types
- `dfg::MetaDiGraph{Int, Float64}`: a data flow graph representing the circuit to be generated
- `handlers::Dict{Operation, AbstractHandler}`: a map from operation type to a hardware generation handler.

See also: [HW.generate](@ref)
"""
@kwdef mutable struct Module{T}
    fn::T
    name::Symbol = _nameof(fn)
    bitwidth::@NamedTuple{integral::Int, fractional::Int} = (integral = 1, fractional = 0)
    parameters::Dict{String, Any} = Dict{String, Any}()
    submodules::Vector{Type} = Type[]
    dfg::MetaDiGraph{Int, Float64} = MetaDiGraph()
    handlers::Dict{Any, Tuple{Any, Any}} = Dict{Any, Tuple{Any, Any}}()
end

Base.show(io::IO, m::Module) =
    print(io, "Module $(m.name) with $(length(m.parameters)) parameters and $(length(m.submodules)) submodules.")
Base.show(io::IO, ::MIME"text/plain", m::Module) = print(io, """
    Module $(m.name):
        Parameters:
            $(m.parameters)
        Submodules:
            $(m.submodules)
        Number of operations: $(nv(m.dfg))
        Number of inputs: $(length(map(x -> getinputs(m.dfg, x), getroots(m.dfg))))
        Number of outputs: $(length(map(x -> getoutputs(m.dfg, x), getbuds(m.dfg))))
    """)

getinputs(g::MetaDiGraph, v) = get_prop(g, v, :inputs)
getoutputs(g::MetaDiGraph, v) = get_prop(g, v, :outputs)
getoperator(g::MetaDiGraph, v) = get_prop(g, v, :operator)

findnode(g::MetaDiGraph, inputs, outputs, operator) =
    collect(filter_vertices(g, (g, v) -> all(getname.(getinputs(g, v)) .== inputs) &&
                                         all(getname.(getoutputs(g, v)) .== outputs) &&
                                         get_prop(g, v, :operator) == operator))
getroots(g::MetaDiGraph) = collect(filter_vertices(g, (g, v) -> isempty(inneighbors(g, v))))
getbuds(g::MetaDiGraph) = collect(filter_vertices(g, (g, v) -> isempty(outneighbors(g, v))))
getparents(g::MetaDiGraph, x::Net) =
    filter_vertices(g, (g, v) -> x ∈ getoutputs(g, v))
getchildren(g::MetaDiGraph, x::Net) =
    filter_vertices(g, (g, v) -> x ∈ getinputs(g, v))
function traverse(g::MetaDiGraph, vs::Vector{T}, visited = T[]) where T
    levelup = unique(reduce(vcat, map(v -> outneighbors(g, v), vs)))
    parents = filter(v -> all(x -> x ∈ visited, inneighbors(g, v)), levelup)

    return parents, union(parents, visited)
end

function addnode!(m::Module, inputs::Vector{Net}, outputs::Vector{Net}, op::Operation)
    add_vertex!(m.dfg)
    node = nv(m.dfg)
    set_prop!(m.dfg, node, :inputs, inputs)
    set_prop!(m.dfg, node, :outputs, outputs)
    set_prop!(m.dfg, node, :operator, op)

    for input in inputs
        for parent in getparents(m.dfg, input)
            add_edge!(m.dfg, parent, node)
        end
    end

    for output in outputs
        for child in getchildren(m.dfg, output)
            add_edge!(m.dfg, node, child)
        end
    end

    return m
end

function getnetlist(m::Module)
    netlist = Net[]
    for v in vertices(m.dfg)
        inputs = getinputs(m.dfg, v)
        outputs = getoutputs(m.dfg, v)
        append!(netlist, inputs)
        append!(netlist, outputs)
    end

    return netlist |> unique
end

function printdfg(m::Module)
    nodes = getroots(m.dfg)
    visited = nodes
    padding = ""

    while !isempty(nodes)
        for node in nodes
            inputs = getinputs(m.dfg, node)
            outputs = getoutputs(m.dfg, node)
            op = getoperator(m.dfg, node)
            println("$(padding)inputs: $inputs")
            println("$(padding)outputs: $outputs")
            println("$(padding)op: $op")
            println("")
        end

        nodes, visited = traverse(m.dfg, nodes, visited)
        padding *= "   "
    end
end

function _printnet(net::Net)
    reg = isreg(net) ? "reg" : ""
    bitlength = bitwidth(net) * prod(netsize(net)) - 1
    bitstr = (bitlength == 0) ? "" : "[$bitlength:0]"
    names = issigned(net) ? "$(name(net))_p, $(name(net))_m" : name(net)

    return join([reg, bitstr, names], " ")
end

_encodeparameter(buffer, name, value) =
    write(buffer, "parameter $name = $value;")
function _encodeparameter(buffer, name, value::AbstractArray)
    for (i, sz) in enumerate(size(value))
        write(buffer, "parameter $(name)_sz_$i = $sz;\n")
    end
    write(buffer, "parameter $name = {")
    write(buffer, join(value, ", "))
    write(buffer, "};")
end
function _encodeparameter(buffer, name, value::Tuple)
    write(buffer, "parameter $(name)_sz_1 = $(length(value));\n")
    write(buffer, "parameter $name = {")
    write(buffer, join(value, ", "))
    write(buffer, "};")
end

function _generatetopmatter(buffer, m::Module, netlist::Netlist)
    netlist = unique(netlist)
    inputs = filter(isinput, netlist)
    outputs = filter(isoutput, netlist)
    # parameters = filter(isparameter, netlist)
    internals = filter(isinternal, netlist)
    wires = filter(iswire, internals)
    regs = filter(isreg, internals)

    write(buffer, "module $(m.name)(CLK, nRST, ")

    write(buffer, join(map(inputs) do input
        issigned(input) ? "$(name(input))_p, $(name(input))_m" : name(input)
    end, ", "))
    write(buffer, ", ")
    write(buffer, join(map(outputs) do output
        issigned(output) ? "$(name(output))_p, $(name(output))_m" : name(output)
    end, ", "))
    write(buffer, ");\n\n")

    for (name, value) in m.parameters
        _encodeparameter(buffer, name, value)
        write(buffer, "\n")
    end
    write(buffer, "\n")

    write(buffer, "input CLK, nRST;\n")

    write(buffer, join(map(inputs) do input
        "input $(_printnet(input));"
    end, "\n"))
    write(buffer, "\n")
    write(buffer, join(map(outputs) do output
        "output $(_printnet(output));"
    end, "\n"))
    write(buffer, "\n")
    write(buffer, join(map(_printnet, regs), "\n"))
    write(buffer, "\n")
    write(buffer, join(map(wires) do wire
        "wire $(_printnet(wire));"
    end, "\n"))
    write(buffer, "\n\n")

    return buffer
end

function _sync_nodes!(nets, netlist)
    for (i, net) in enumerate(nets)
        j = find(netlist, net)
        if !isempty(j)
            nets[i] = only(netlist[j])
        end
    end

    return nets
end

"""
    HW.generate(m::Module, netlist::Netlist)
    HW.generate(m::Module, f)
    HW.generate(c::Tuple{Module, Function}, dut, args...)

Generate the Verilog implementation of the module.
Users will most likely call the last method form above.

# Fields:
- `m::Module`: the module to generate
- `netlist::Netlist`: the netlist for the circuit being generated
- `f`: a closure with one argument (a netlist) that calls a runtime information extraction function
- `c::Tuple{Module, Function}`: the tuple returned by [`@circuit`](@ref)
- `dut`: an instance of the circuit struct
- `args`: example arguments to circuit
"""
function generateverilog(io::IO, m::Module)
    outstr = ""
    netlist = getnetlist(m)

    # start at inputs
    nodes = getroots(m.dfg)
    visited = nodes

    buffer = (io isa IOBuffer) ? IOBuffer() : mktemp()[2]
    while !isempty(nodes)
        # for each node, invoke the appropriate handler
        for node in nodes
            inputs = getinputs(m.dfg, node)
            outputs = getoutputs(m.dfg, node)
            operator = getoperator(m.dfg, node)

            handler = gethandler(operator.broadcasted, operator.type, jltypeof.(inputs)...)
            handler, state = get!(m.handlers, typeof(handler), (handler, init_state(handler)))

            set_prop!(m.dfg, node, :inputs, _sync_nodes!(inputs, netlist))
            _, state = handler(buffer, netlist, state, inputs, outputs)
            m.handlers[typeof(handler)] = (handler, state)
            set_prop!(m.dfg, node, :outputs, _sync_nodes!(outputs, netlist))
        end

        # move one level up in the DFG
        nodes, visited = traverse(m.dfg, nodes, visited)
    end

    # set outputs
    for v in getbuds(m.dfg)
        for output in getoutputs(m.dfg, v)
            setclass!(netlist, output, :output)
        end
    end

    # print top matter
    _generatetopmatter(io, m, netlist)
    # print main matter
    seekstart(buffer)
    for line in eachline(buffer)
        write(io, line)
        write(io, "\n")
    end
    write(io, "\nendmodule\n")

    return io
end
generateverilog(m::Module) = generateverilog(IOBuffer(), m)


text:
struct Prior{F <: AbstractFloat}
  K::Int
  p::Beta
  gamma::Beta
  eta::Dirichlet
  mu::Normal
  a_omega::F
  nu::LogNormal
  psi::Normal
  tau::Gamma
end

function compute_prior_mu(data::Data)
  yfinite = [data.yC_finite; data.yT_finite]
  return Normal(mean(yfinite), std(yfinite))
end

function Prior(K; p=Beta(100, 100), gamma=Beta(1,1), eta=nothing,
               mu=Normal(0, 3), a_omega=2.5, tau=Gamma(0.5, 1),
               nu=LogNormal(0.5, 0.1), # nu=LogNormal(1.6, 0.4),
               psi=Normal(-1, 1), data=nothing)
  eta == nothing && (eta = Dirichlet(K, 1/K))
  data == nothing || (mu = compute_prior_mu(data))
  return Prior(K, p, gamma, eta, mu, a_omega, nu, psi, tau)
end

#= TEST:
prior = Prior(3)
for fn in fieldnames(Prior)
  println("$(fn) => $(getfield(prior, fn))")
end
=#

max_stars_repo_name: luiarthur/CytofDensityEstimation

max_issues_repo_name: luiarthur/CytofDensityEstimation

max_forks_repo_name: luiarthur/CytofDensityEstimation


text:
##############################################################################
##
## isnested
##
##############################################################################

function isnested(fe::FixedEffect, prefs)
    entries = zeros(eltype(prefs), fe.n)
    for (feref, pref) in zip(fe.refs, prefs)
        if entries[feref] == 0
            # it's a new level, create entry
            entries[feref] = pref
        elseif entries[feref] != pref
            # not nested: for the same level in a, two different levels in b
            return false
        end
    end
    return true
end

##############################################################################
##
## Number of distinct values
##
##############################################################################

function nunique(fe::FixedEffect)
    out = zeros(Int, fe.n)
    for ref in fe.refs
        out[ref] += 1
    end
    sum(x -> x > 0, out)
end

##############################################################################
##
## Drop Singletons
##
##############################################################################

function drop_singletons!(esample, fe::FixedEffect)
    cache = zeros(Int, fe.n)
    for i in eachindex(esample)
        if esample[i]
            cache[fe.refs[i]] += 1
        end
    end
    for i in eachindex(esample)
        if esample[i] && cache[fe.refs[i]] <= 1
            esample[i] = false
        end
    end
end

##############################################################################
##
## Subet and Make sure Interaction if Vector{Float64} (instead of missing)
##
##############################################################################

# index and convert interaction Vector{Float64,, Missing} to Vector{Missing}
function _subset(fe::FixedEffect, esample)
    interaction = convert(AbstractVector{Float64}, fe.interaction[esample])
    FixedEffect{typeof(fe.refs), typeof(interaction)}(fe.refs[esample], interaction, fe.n)
end


text:
function pcd!(rbm::ConvRBM, data::AbstractArray;
    batchsize::Int = 1,
    epochs::Int = 1,
    optim = Flux.ADAM(),
    history::MVHistory = MVHistory(),
    wts = nothing,
    steps::Int = 1,
    vm::AbstractArray = initial_fantasy_v(rbm, data, batchsize),
    callback = Returns(nothing)
)
    @assert length(vsizes(rbm, data).batch_size) == 1 # checks sizes and require single batch dim
    @assert isnothing(wts) || _nobs(data) == _nobs(wts)
    for epoch in 1:epochs
        batches = RBMs.minibatches(data, wts; batchsize)
        Δt = @elapsed for (batch_idx, (vd, wd)) in enumerate(batches)
            vm .= RBMs.sample_v_from_v(rbm, vm; steps)
            gs = Zygote.gradient(rbm) do rbm
                contrastive_divergence(rbm, vd, vm; wd)
            end
            # removes keys with no gradient
            ∂ = Base.structdiff(only(gs), NamedTuple{(:stride, :pad, :dilation, :groups, :pool)})
            Δ = RBMs.update!(∂, rbm, optim) # update step
            callback(; rbm, epoch, batch_idx, vm, vd, wd, history, ∂, Δ)
            RBMs.update!(rbm, Δ)
            push!(history, :∂, RBMs.gradnorms(∂))
            push!(history, :Δ, RBMs.gradnorms(∂))
        end
        push!(history, :epoch, epoch)
        push!(history, :Δt, Δt)
        push!(history, :vm, copy(vm))
        @debug "epoch $epoch/$epochs ($(round(Δt, digits=2))s)"
    end
    return history
end

function RBMs.update!(∂::NamedTuple, x::ConvRBM, optim)
    for (k, g) in pairs(∂)
        if hasproperty(x, k)
            RBMs.update!(g, getproperty(x, k), optim)
        else
            g .= 0
        end
    end
    return ∂
end

function RBMs.update!(x::ConvRBM, ∂::NamedTuple)
    for (k, Δ) in pairs(∂)
        hasproperty(x, k) && RBMs.update!(getproperty(x, k), Δ)
    end
    return x
end

function contrastive_divergence(
    rbm::ConvRBM, vd::AbstractArray, vm::AbstractArray; wd = nothing, wm = nothing
)
    Fd = RBMs.mean_free_energy(rbm, vd; wts = wd)
    Fm = RBMs.mean_free_energy(rbm, vm; wts = wm)
    return Fd - Fm
end

function RBMs.mean_free_energy(rbm::ConvRBM, v::AbstractArray; wts = nothing)::Number
    F = RBMs.free_energy(rbm, v)
    if isempty(vsizes(rbm, v).batch_size)
        wts::Nothing
        return F
    else
        return RBMs.wmean(F; wts)
    end
end

function initial_fantasy_v(rbm::ConvRBM, data::AbstractArray, batchsize::Int)
    inputs = falses(size(data)[1:(end - 1)]..., batchsize)
    vm = RBMs.transfer_sample(visible(rbm), inputs)
    return oftype(data, vm)
end


text:
# include("../../src/SFFM.jl")
using LinearAlgebra, Plots, SFFM

## define the model(s)
include("exampleModelDef.jl")

## section 4.3: the marginal stationary distribution of X
## mesh
Δtemp = 0.4
Nodes = collect(approxBounds[1, 1]:Δtemp:approxBounds[1, 2])
nBases = 2
Basis = "lagrange"
mesh = SFFM.DGMesh(
    approxModel,
    Nodes,
    nBases,
    Basis = Basis,
)
let q = SFFM.plot(approxModel)
    ## analytic version for comparison
    # construction
    Ψₓ = SFFM.PsiFunX(approxModel)
    ξₓ = SFFM.MakeXiX(approxModel, Ψₓ)
    pₓ, πₓ, Πₓ, Kₓ = SFFM.StationaryDistributionX( approxModel, Ψₓ, ξₓ)

    # evaluate the distribution
    analyticX = SFFM.SFFMDensity(
        [pₓ[:];0;0],
        πₓ(SFFM.CellNodes(mesh)),
        SFFM.CellNodes(mesh),
    )

    # plot it
    q = SFFM.plot!(q,
        approxModel,
        mesh,
        analyticX,
        color = :red,
        label = "Analytic",
        marker = :x,
        seriestype = :line,
        jitter = 0.5,
    )

    # now DG it
    c = 0
    colours = [:green;:blue]
    for nBases in 1:2
        c = c+1
        mesh = SFFM.DGMesh(
            approxModel,
            Nodes,
            nBases,
            Basis=Basis,
        )

        # compute the marginal via DG
        All = SFFM.MakeAll( approxModel, mesh, approxType = "projection")
        Ψ = SFFM.PsiFun(All.D)

        # the distribution of X when Y first returns to 0
        ξ = SFFM.MakeXi(All.B, Ψ)

        marginalX, p, K = SFFM.MakeLimitDistMatrices(
            All.B,
            All.D,
            All.R.RDict,
            Ψ,
            ξ,
            mesh,
            approxModel,
        )
        # convert marginalX to a distribution for plotting
        Dist = SFFM.Coeffs2Dist(
            approxModel,
            mesh,
            marginalX,
            SFFM.SFFMDensity,
        )
        # plot it
        q = SFFM.plot!(
            q,
            approxModel,
            mesh,
            Dist,
            color = colours[c],
            label = "DG: N_k = "*string(nBases),
            seriestype = :line,
            jitter = 0.5,
        )

    end
    q = plot!(windowsize=(600,600))
    titles = ["Phase 11" "Phase 10" "Phase 01" "Phase 00"]
    for sp in 1:4
        q = plot!(
            subplot = sp,
            xlims = (-0.5,5),
            title = titles[sp],
            xlabel = "x",
            ylabel = "Density / Probability",
        )
    end
    display(q)
    # savefig(pwd()*"/examples/paperNumerics/dump/marginalStationaryDistX.png")
end

max_stars_repo_name: angus-lewis/SFFM

max_issues_repo_name: angus-lewis/SFFM

max_forks_repo_name: angus-lewis/SFFM


text:
################
# AbstNS
################

abstract type AbstNS end

################
# NSdict0
################

const _NSdict0 = Dict{Symbol, Function}()

_NSdict0[:_fixed] = ns -> ns.__fix_lck[1]
_NSdict0[:_lcked] = ns -> ns.__fix_lck[2]
_NSdict0[:_frzed] = ns -> all(ns.__fix_lck)
_NSdict0[:_keys] = ns -> Tuple(Base.keys(ns.__dict))
_NSdict0[:_vals] = ns -> Tuple(x.obj for x ∈ Base.values(ns.__dict))
_NSdict0[:_keyvals] = ns -> (; zip(ns._keys, ns._vals)...)
_NSdict0[:_printkeyvals] =
    ns -> ([println(k, ": ", v) for (k,v) ∈ pairs(ns._keyvals)]; return)
_NSdict0[:_printkeytypes] =
    ns -> ([println(k, ": ", typeof(v)) for (k,v) ∈ pairs(ns._keyvals)]; return)
_NSdict0[:_fix] = ns -> (ns.__fix_lck[1] = true; ns)
_NSdict0[:_unfix] = ns -> (ns.__fix_lck[1] = false; ns)
_NSdict0[:_lck] = ns -> (ns.__fix_lck[2] = true;  ns)
_NSdict0[:_unlck] = ns -> (ns.__fix_lck[2] = false; ns)
_NSdict0[:_frz] = ns -> (ns.__fix_lck[1] = ns.__fix_lck[2] = true; ns)
_NSdict0[:_unfrz] = ns -> (ns.__fix_lck[1] = ns.__fix_lck[2] = false; ns)
_NSdict0[:_clr] = ns -> (ns.del(); ns)
_NSdict0[:_copy] = ns -> deepcopy(ns)
_NSdict0[:_cst_keys] =
    ns -> (d = ns.__dict; [k for k ∈ Base.keys(d) if isa(d[k], NScst_item)])
_NSdict0[:_noncst_keys] =
    ns -> (d = ns.__dict; [k for k ∈ Base.keys(d) if isa(d[k], NSnoncst_item)])
_NSdict0[:import] = ns ->
    ((g::AbstNS, a::Vararg{Symbol}; exclude=[], deep=false) ->
     begin
         d = ns.__dict
         gd = g.__dict
         if deep
             if length(a) > 0
                 for k ∈ a
                     k ∉ exclude && (d[k] = deepcopy(gd[k]))
                 end
             else
                 for (k, v) ∈ pairs(gd)
                     k ∉ exclude && (d[k] = deepcopy(v))
                 end
             end
         else
             if length(a) > 0
                 for k ∈ a
                     k ∉ exclude && (d[k] = gd[k])
                 end
             else
                 for (k, v) ∈ pairs(gd)
                     k ∉ exclude && (d[k] = v)
                 end
             end
         end
         ns
     end)
_NSdict0[:copyout] = ns ->
    ((a::Vararg{Symbol}; exclude=[], deep=false) ->
     begin
         g = typeof(ns)()
         d = ns.__dict
         gd = g.__dict
         if deep
             if length(a) > 0
                 for k ∈ a
                     k ∉ exclude &&
                         (gd[k] = (ns.iscst(k)
                                   ? NScst_item
                                   : NSnoncst_item)(deepcopy(d[k].obj)))
                 end
             else
                 for (k, w) ∈ pairs(d)
                     k ∉ exclude &&
                         (gd[k] = (ns.iscst(k)
                                   ? NScst_item
                                   : NSnoncst_item)(deepcopy(w.obj)))
                 end
             end
         else
             if length(a) > 0
                 for k ∈ a
                     k ∉ exclude &&
                         (gd[k] = (ns.iscst(k)
                                   ? NScst_item
                                   : NSnoncst_item)(d[k].obj))
                 end
             else
                 for (k, w) ∈ pairs(d)
                     k ∉ exclude &&
                         (gd[k] = (ns.iscst(k)
                                   ? NScst_item
                                   : NSnoncst_item)(w.obj))
                 end
             end
         end
         g
     end)
_NSdict0[:copyfrom] = ns ->
    ((g::AbstNS, a::Vararg{Symbol}; exclude=[], deep=false) ->
     # ns.import(g.copyout(a...; exclude=exclude, deep=deep))
     begin
         d = ns.__dict
         gd = g.__dict
         if deep
             if length(a) > 0
                 for k ∈ a
                    k ∉ exclude &&
                         (d[k] = (g.iscst(k)
                                  ? NScst_item
                                  : NSnoncst_item)(deepcopy(gd[k].obj)))
                 end
             else
                 for (k, w) ∈ pairs(gd)
                     k ∉ exclude &&
                         (d[k] = (g.iscst(k)
                                  ? NScst_item
                                  : NSnoncst_item)(deepcopy(w.obj)))
                 end
             end
         else
             if length(a) > 0
                 for k ∈ a
                     k ∉ exclude &&
                         (d[k] = (g.iscst(k)
                                  ? NScst_item
                                  : NSnoncst_item)(gd[k].obj))
                 end
             else
                 for (k, w) ∈ pairs(gd)
                     k ∉ exclude &&
                         (d[k] = (g.iscst(k)
                                  ? NScst_item
                                  : NSnoncst_item)(w.obj))
                 end
             end
         end
         ns
     end)
# g.export()
#     : export all properties from g to new ns
# g.export(:a, :b, :c)
#     : export properties :a, :b, :c from g to new ns
_NSdict0[:export] = ns ->
    ((a::Vararg{Symbol}; exclude=[], deep=false) ->
     begin
         g = typeof(ns)()
         d = ns.__dict
         gd = g.__dict
         if deep
             if length(a) > 0
                 for k ∈ a
                     k ∉ exclude && (gd[k] = deepcopy(d[k]))
                 end
             else
                 for (k, v) ∈ pairs(d)
                     k ∉ exclude && (gd[k] = deepcopy(v))
                 end
             end
         else
             if length(a) > 0
                 for k ∈ a
                     k ∉ exclude && (gd[k] = d[k])
                 end
             else
                 for (k, v) ∈ pairs(ns.__dict)
                     k ∉ exclude && (gd[k] = v)
                 end
             end
         end
         g
     end)
_NSdict0[:deepimport] = ns ->
    ((g::AbstNS, a::Vararg{Symbol}; exclude=[]) ->
     ns.import(g, a...; exclude=exclude, deep=true))
_NSdict0[:deepexport] = ns ->
    ((a::Vararg{Symbol}; exclude=[]) -> ns.export(a...; exclude=exclude, deep=true))

# g.tontpl()
#     : export all properties from g to a NamedTuple
# g.tontpl(:a, :b, :c)
#     : export properties :a, :b, :c from g to a NamedTuple
_NSdict0[:tontpl] = ns ->
    ((a::Vararg{Symbol}; exclude=[], deep=false) ->
    begin
        gd = Dict()
         if deep
             if length(a) > 0
                 for k ∈ a
                     k ∉ exclude && (gd[k] = deepcopy(Base.getproperty(ns,k)))
                 end
             else
                 for k ∈ ns._keys
                     k ∉ exclude && (gd[k] = deepcopy(Base.getproperty(ns,k)))
                 end
             end
         else
             if length(a) > 0
                 for k ∈ a
                     k ∉ exclude && (gd[k] = Base.getproperty(ns,k))
                 end
             else
                 for k ∈ ns._keys
                     k ∉ exclude && (gd[k] = Base.getproperty(ns,k))
                 end
             end
         end
        (;gd...)
        end);

# g.load("x.ns")
#     : load "x.ns" and import all properties from it
# g.load("x.ns", :a, :b, :c)
#     : load "x.ns" and import properties :a, :b, :c from it
_NSdict0[:load] = ns ->
    ((filename::AbstractString,
      atr::Vararg{Symbol};
      exclude=[],
      forcename=false) ->
     begin
         if !forcename &&
            (length(filename) < length("a.ns") ||
             filename[end-length(".ns")+1:end] != ".ns")
             filename = filename * ".ns"
         end
         ns.import(Serialization.deserialize(filename),
                   atr...;
                   exclude=exclude)
     end)
_NSdict0[:save] = ns ->
    ((filename::AbstractString,
      atr::Vararg{Symbol};
      exclude=[],
      forcename=false) ->
     begin
         if !forcename &&
            (length(filename) < length("a.ns") ||
             filename[end-length(".ns")+1:end] != ".ns")
             filename = filename * ".ns"
         end
         length(atr) == 0 && (atr = ns._keys)
         atr = [k for k ∈ ns._keys if k ∉ exclude]
         g = ns.copyout(atr...; exclude=exclude)
         Serialization.serialize(filename, g)
         g
     end)
_NSdict0[:haskey] = ns -> NShaskey(ns)
_NSdict0[:iscst] = ns ->
    (key -> (!haskey(ns.__dict, key) &&
             error("this NS does not have a key named '$(atr)'." );
             isa(ns.__dict[key], NScst_item)))
_NSdict0[:del] = ns -> NSdel(ns)
_NSdict0[:cstize] = ns -> NScstize(ns)
_NSdict0[:decstize] = ns -> NSdecstize(ns)

# tags
_NSdict0[:cst] = ns -> NScst(ns)
_NSdict0[:prp] = ns -> NSTag{:prp, false}(ns)
_NSdict0[:fnc] = ns -> NSTag{:fnc, false}(ns)
_NSdict0[:mth] = ns -> NSTag{:mth, false}(ns)
_NSdict0[:dfn] = ns -> NSTag{:dfn, false}(ns)
_NSdict0[:req] = ns -> NSTag{:req, false}(ns)


text:
function GeneralFluxLimiter(func_Ψ,u,grid::FVMPDEGrid,dim)
    ϵ = 1.0e-8
    indices = grid.indices
    Δx = grid.dX
    Δx_C = grid.dX_C

    func_r(Uₓ₁,Uₓ₂) = (Uₓ₁ .+ ϵ) ./ (Uₓ₂ .+ ϵ)
    func_U(Δxᵢ,r,uₓ) = (Δxᵢ ./ 2) .* func_Ψ(r) .* (uₓ .+ ϵ)

    d = length(size(u)) - 1
    #nx,ny,nz,nvars = size(u)

    minIndex = fill(Tuple(ones(Int8,d)),size(indices))
    maxIndex = fill(Tuple(size(u)[1:end-1]),size(indices))

    adderArray = zeros(Int8,d)
    adderArray[dim] = 1
    adder = CartesianIndex(Tuple(adderArray))

    ViewUᵢ₋₂ = CartesianIndex.(max.(Tuple.(indices .- adder .- adder),minIndex))
    ViewUᵢ₋₁ = CartesianIndex.(max.(Tuple.(indices .- adder         ),minIndex))
    ViewUᵢ₊₁ = CartesianIndex.(min.(Tuple.(indices .+ adder         ),maxIndex))
    ViewUᵢ₊₂ = CartesianIndex.(min.(Tuple.(indices .+ adder .+ adder),maxIndex))

    #println(size( (u[ViewUᵢ₊₂,:] .- u[ViewUᵢ₊₁,:])))
    #println(size(Δx))
    #println(size(Δx_C))

    Δx_extended = [Δx[1];Δx[1];Δx;Δx[end];Δx[end]]

    Uₓi₊¾ = (u[ViewUᵢ₊₂,:] .- u[ViewUᵢ₊₁,:]) ./ Δx_extended[4:end-0]    # ¾ : 3/2
    Uₓi₊½ = (u[ViewUᵢ₊₁,:] .- u) ./ Δx_extended[3:end-1]
    Uₓi₋½ = (u .- u[ViewUᵢ₋₁,:]) ./ Δx_extended[2:end-2]
    Uₓi₋¾ = (u[ViewUᵢ₋₁,:] .- u[ViewUᵢ₋₂,:]) ./ Δx_extended[1:end-3]

    Δx_C_extended = [Δx_C[1];Δx_C;Δx_C[end]]
    Uᴸi₊½ = u[ViewUᵢ₊₁,:] .- func_U(Δx_C_extended[3:end-0],func_r(Uₓi₊½,Uₓi₊¾),Uₓi₊¾)
    Uᴸi₋½ = u             .- func_U(Δx_C_extended[2:end-1],func_r(Uₓi₋½,Uₓi₊½),Uₓi₊½)
    Uᴿi₊½ = u             .+ func_U(Δx_C_extended[2:end-1],func_r(Uₓi₊½,Uₓi₋½),Uₓi₋½)
    Uᴿi₋½ = u[ViewUᵢ₋₁,:] .+ func_U(Δx_C_extended[1:end-2],func_r(Uₓi₋½,Uₓi₋¾),Uₓi₋¾)

    Uₓᴸ = (Uᴸi₊½ .- Uᴸi₋½) ./ Δx_C
    Uₓᴿ = (Uᴿi₊½ .- Uᴿi₋½) ./ Δx_C
    #return Uᴿi₊½, Uᴿi₋½, Uᴸi₊½ , Uᴸi₋½
    #return Uₓi₊¾, Uₓi₊½, Uₓi₋½, Uₓi₋¾
    return Uₓᴸ, Uₓᴿ
end

################## LIMITERS ####################

function κ_Scheme(r,κ)
    α₁ = (1+κ)/2
    α₂ = (1-κ)/2
    return (α₁ .*  r .+ α₂)
end

function uw1(r)
    return zeros(size(r))
end

function uw2(r)
    return ones(size(r))
end

function uw3(r)
    return κ_Scheme(r,1/3)
end

function uw4(r)
    return κ_Scheme(r,1/2)
end

function scd(r)
    return r
end

function fr(r)
    return κ_Scheme(r,0)
end

function kn(r)
    return @. max(0,min(2*r,min(1/3 + 2 /3 * r,2)))
end

function sb(r)
    return @. max(0,min(2 * r,1),min(r,2))
end

function mm(r)
    return @. max(0,min(r,1))
end

function mu(r)
    return @. max(0,min(2 * r,(r+1)/2,2))
end

function ha(r)
    return @. (r + abs(r))/(r+1)
end

function va1(r)
    return @. r*(r + 1)/(r^2+1)
end

function va2(r)
    return @. r*2/(r^2+1)
end

function vl(r)
    return @. (r + abs(r))/(1 + abs(r))
end

function op(r)
    return @. 3*r*(r + 1)/(r^2+r+1)/2
end

function hc(r)
    return @. 1.5*(r + abs(r))/(r+2)
end

function hq(r)
    return @. 2*(r + abs(r))/(r+3)
end

function cm(r)
    return @. max(0,r*(3*r+1)/(r+1)^2)
end

function mc(r)
    return @. max(0,min(2*r,0.5*(r+1),2))
end

function sm(r)
    return @. max(0,min(2*r,(0.75*r+0.25),4))
end

function um(r)
    return @. max(0,min(2*r,(0.25*r+0.75),2))
end


text:
# Factor Graph OS type utilities
#  IIF methods should direclty detect extended types from user import
# of convert in their namespace

import DistributedFactorGraphs: AbstractPointParametricEst

manikde!(pts::AbstractArray{Float64,2}, vartype::InferenceVariable) = manikde!(pts, getManifolds(vartype))
manikde!(pts::AbstractArray{Float64,2}, vartype::Type{<:InferenceVariable}) = manikde!(pts, getManifolds(vartype))

"""
    $SIGNATURES

Return N=100 measurement samples for a factor in `<:AbstractDFG`.
"""
function getMeasurements(dfg::AbstractDFG, fsym::Symbol, N::Int=100)
  fnc = getFactorFunction(dfg, fsym)
  getSample(fnc, N)
end

"""
    $SIGNATURES

Get graph node (variable or factor) dimension.
"""
getDimension(var::DFGVariable) = getSofttype(var).dims
getDimension(fct::DFGFactor) = solverData(fct).fnc.zDim

"""
    $SIGNATURES

Get the ParametricPointEstimates---based on full marginal belief estimates---of a variable in the distributed factor graph.

DevNotes
- TODO update for manifold subgroups.
"""
function calcVariablePPE(var::DFGVariable,
                         softt::InferenceVariable;
                         solveKey::Symbol=:default,
                         method::Type{MeanMaxPPE}=MeanMaxPPE  )::MeanMaxPPE
  #
  P = getKDE(var)
  manis = getManifolds(softt) # getManifolds(vnd)
  ops = buildHybridManifoldCallbacks(manis)
  Pme = getKDEMean(P) #, addop=ops[1], diffop=ops[2]
  Pma = getKDEMax(P, addop=ops[1], diffop=ops[2])
  suggested = zeros(getDimension(var))
  # TODO standardize after AMP3D
  @assert length(manis) == getDimension(var)
  for i in 1:length(manis)
    mani = manis[i]
    if mani == :Euclid
      suggested[i] = Pme[i]
    elseif mani == :Circular
      suggested[i] = Pma[i]
    else
      error("Unknown manifold to find PPE, $softt, $mani")
    end
  end
  MeanMaxPPE(solveKey, suggested, Pma, Pme, now())
end
# function calcVariablePPE!(retval::Vector{Float64},
#                           var::DFGVariable,
#                           softt::InferenceVariable;
#                           method::Type{MeanMaxPPE}=MeanMaxPPE )::Nothing
#   #
#   P = getKDE(var)
#   manis = getManifolds(softt) # getManifolds(vnd)
#   ops = buildHybridManifoldCallbacks(manis)
#   Pme = getKDEMean(P, addop=ops[1], diffop=ops[2])
#   Pma = getKDEMax(P, addop=ops[1], diffop=ops[2])
#   for i in 1:length(manis)
#     mani = manis[i]
#     if mani == :Euclid
#       retval[i] = Pme[i]
#     elseif mani == :Circular
#       retval[i] = Pma[i]
#     else
#       error("Unknown manifold to find PPE, $softt, $mani")
#     end
#   end
#   nothing
# end
# """
#     $SIGNATURES
#
# Get the ParametricPointEstimates---based on full marginal belief estimates---of a variable in the distributed factor graph.
# """
# function calcVariablePPE(var::DFGVariable,
#                          softt::InferenceVariable;
#                          method::Type{<:AbstractPointParametricEst}=MeanMaxPPE  )::Vector{Float64}
#   #
#   # vect = zeros(softt.dims)
#   mmppe = calcVariablePPE(MeanMaxPPE, var, softt, method=method)
#   return mmppe.suggested
# end

# calcVariablePPE!(retvec::Vector{Float64}, var::DFGVariable; method::Type{<:AbstractPointParametricEst}=MeanMaxPPE) = calcVariablePPE!(retvec, var, getSofttype(var), method=method)
calcVariablePPE(var::DFGVariable; method::Type{<:AbstractPointParametricEst}=MeanMaxPPE, solveKey::Symbol=:default) = calcVariablePPE(var, getSofttype(var), method=method, solveKey=solveKey)
function calcVariablePPE(dfg::AbstractDFG, sym::Symbol; method::Type{<:AbstractPointParametricEst}=MeanMaxPPE, solveKey::Symbol=:default )
  var = getVariable(dfg, sym)
  calcVariablePPE(var, getSofttype(var), method=method, solveKey=solveKey)
end

"""
    $(SIGNATURES)

Return array of all variable nodes connected to the last `n` many poses (`:x*`).

Example:

```julia
# Shallow copy the tail end of poses from a factor graph `fg1`
vars = lsRear(fg1, 5)
fg1_r5 = subgraphFromVerts(fg1, vars)
```
"""
function lsRear(fgl::FactorGraph, n::Int=1)
  @warn "lsRear in current form is not gauranteed to work right, use sort(ls(fg, r\"x\")[end]) instead.  Also see sortVarNested"
  lasts = ls(fgl)[1][(end-n):end]
  syms = ls(fgl, lasts)
  union(lsf.(fgl, syms)[:]...)
end

"""
    $SIGNATURES

Return `Vector{Symbol}` of landmarks attached to vertex vsym in `fgl::FactorGraph`.
"""
function landmarks(fgl::FactorGraph, vsym::Symbol)
  fsyms = ls(fgl, vsym)
  lms = Symbol[]
  for fs in fsyms
    for varv = lsf(fgl, fs)
      if string(varv)[1] == 'l'
        push!(lms, varv)
      end
    end
  end
  lms
end

function evalLikelihood(fg::FactorGraph, sym::Symbol, point::Vector{Float64})
  p = getVertKDE(fg, sym)
  Ndim(p) == length(point) ? nothing : error("point (dim=$(length(point))) must have same dimension as belief (dim=$(Ndim(p)))")
  evaluateDualTree(p, reshape(point,:,1))[1]
end

# Evaluate the likelihood of an Array{2} of points on the marginal belief of some variable
# note the dimensions must match
function evalLikelihood(fg::FactorGraph, sym::Symbol, points::Array{Float64,2})
  p = getVertKDE(fg, sym)
  Ndim(p) == size(points,1) ? nothing : error("points (dim=$(size(points,1))) must have same dimension as belief (dim=$(Ndim(p)))")
  evaluateDualTree(p, (points))
end

function setThreadModel!(fgl::FactorGraph;model=IncrementalInference.SingleThreaded)
  for (key, id) in fgl.fIDs
    solverData(getFactor(fgl, key)).fnc.threadmodel = model
  end
  nothing
end

# not sure if and where this is still being used
function _evalType(pt::String)::Type
    try
        getfield(Main, Symbol(pt))
    catch ex
        io = IOBuffer()
        showerror(io, ex, catch_backtrace())
        err = String(take!(io))
        error("_evalType: Unable to locate factor/distribution type '$pt' in main context (e.g. do a using on all your factor libraries). Please check that this factor type is loaded into main. Stack trace = $err")
    end
end

"""
    $(SIGNATURES)

Print the maximum point values form all variables approximate marginals in the factor graph.
The full marginal can be recovered for example `X0 = getVertKDE(fg, :x0)`.
"""
function printgraphmax(fgl::FactorGraph)
    verts = union(ls(fgl)...)
    map(v -> println("$v : $(getKDEMax(getVertKDE(fgl, v)))"), verts);
end

"""
    $SIGNATURES

Return interger index of desired variable element.

Example
-------
```julia
pp = RoME.Point2()
getIdx(pp, :posY) # = 2
```

Internal Notes
--------------
- uses number i < 100 for index number, and
- uses +100 offsets to track the minibatch number of the requested dimension
"""
function getIdx(pp::T, sym::Symbol, i::Int=0)::Tuple{Int, Int} where {T <: Tuple}
  # i > 99 ? error("stop") : nothing
  i-=100
  for p in pp
    i,j = getIdx(p, sym, i)
    if i > 0
      return i, j
    end
  end
  return i,-1
end
getIdx(pp::Symbol, sym::Symbol, i::Int=0)::Tuple{Int, Int} = pp==sym ? (abs(i)%100+1, div(abs(i)-100,100)) : (i-1, div(abs(i)-100,100))
function getIdx(pp::V, sym::Symbol, i::Int=0)::Tuple{Int, Int} where {V <: InferenceVariable}
  return getIdx(pp.dimtype, sym)
end

"""
   $SIGNATURES

Display the content of `VariableNodeData` to console for a given factor graph and variable tag`::Symbol`.

Dev Notes
- TODO split as two show macros between AMP and DFG
"""
function showVariable(fgl::G, vsym::Symbol) where G <: AbstractDFG
  vert = DFG.getVariable(fg, vsym)
  vnd = solverData(vert)
  println("label: $(vert.label), exVertexId: $(vert.index)")
  println("tags: $( haskey(vert.attributes, string(:tags)) ? vert.attributes[string(:tags)] : string(:none))")
  println("size marginal samples $(size(getVal(vnd)))")
  println("kde bandwidths: $(getBW(vnd)[:,1])")
  println("kde mean: $(round.(getKDEMean(getKDE(vnd)),digits=4))")
  println("kde max: $(round.(getKDEMax(getKDE(vnd)),digits=4))")
  println()
  vnd
end

"""
    $SIGNATURES

Return `::Bool` on whether this variable has been marginalized.
"""
isMarginalized(vert::DFGVariable) = solverData(vert).ismargin
isMarginalized(dfg::AbstractDFG, sym::Symbol) = isMarginalized(DFG.getVariable(dfg, sym))

"""
    $SIGNATURES

Return bool on whether a certain factor has user defined multihypothesis.

Related

getMultihypoDistribution
"""
isMultihypo(fct::DFGFactor) = isa(solverData(fct).fnc.hypotheses, Distribution)

"""
    $SIGNATURES

Return the categorical distributed used for multihypothesis selection in a factor.

Related

isMultihypo
"""
getMultihypoDistribution(fct::DFGFactor) = solverData(fct).fnc.hypotheses

"""
    $SIGNATURES

Free all variables from marginalization.
"""
function dontMarginalizeVariablesAll!(fgl::G) where G <: AbstractDFG
  fgl.solverParams.isfixedlag = false
  fgl.solverParams.qfl = 9999999999
  for sym in ls(fgl)
    solverData(getVariable(fgl, sym)).ismargin = false
  end
  nothing
end

"""
    $SIGNATURES

Free all variables from marginalization.

Related

dontMarginalizeVariablesAll!
"""
function unfreezeVariablesAll!(fgl::G) where G <: AbstractDFG
  dontMarginalizeVariablesAll!(fgl)
end

"""
    $SIGNATURES

Reset initialization flag on all variables in `::FactorGraphs`.

Notes
- Numerical values remain, but inference will overwrite since init flags are now `false`.
"""
function resetVariableAllInitializations!(fgl::FactorGraph)
  vsyms = ls(fgl)
  for sym in vsyms
    setVariableInitialized!(getVariable(fgl, sym), :false)
  end
  nothing
end

function convert(::Type{Tuple{BallTreeDensity,Float64}},
                 p::EasyMessage )
  (AMP.manikde!(p.pts, p.bws, p.manifolds), p.inferdim)
end

function convert(::Type{EasyMessage},
                 bel::Tuple{BallTreeDensity,Float64},
                 manifolds::T) where {T <: Tuple}
  EasyMessage(getPoints(bel[1]), getBW(bel[1])[:,1], manifolds, bel[2])
end

#


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

using Markdown
using InteractiveUtils

# ╔═╡ 0b868b30-6420-11eb-18f9-1592f7542fe6
begin
    using PlutoUI
    using Plots
    using Printf
	using SpecialFunctions
end

# ╔═╡ 4c4fff40-64c8-11eb-2c45-57f4793ee997
md"_Parabolic Partial Differential Equation, version 1_"

# ╔═╡ 56d1f810-64c8-11eb-3e0d-81c0bbafdc62
md"""
## **Parabolic Partial Differential Equation**
#### **Description:**
This notebook computes the Parabolic Partial Differential Equation with 4 different discretization methods (2 explicit - 2 implicit)

		FTCS(Δt,Δy,ν)
		DuFortFrankel(Δt,Δy,ν)
		Laasonen(Δt,Δy,ν)
		CrankNicolson(Δt,Δy,ν)

If the arguments are missing, it will use the default values.

		FTCS(Δt=0.002,Δy=0.001,ν=0.000217)
		DuFortFrankel(Δt=0.002,Δy=0.001,ν=0.000217)
		Laasonen(Δt=0.002,Δy=0.001,ν=0.000217)
		CrankNicolson(Δt=0.002,Δy=0.001,ν=0.000217)

#### **Parabolic Partial Differential Equation**
In this problem, velocity profile between 2 infinite plates will be investigated. The bottom plates is moving with U₀ and the upper plate is stationary. The equation for the system can be represented by:

$$\frac{∂u}{∂t}=α\frac{∂²u}{∂x²}$$

Boundary conditions for the system:

$$u=U₀=1 \hspace{12pt} y=0$$
$$u=0    \hspace{36pt} y=h$$

where h is the distance between two plates. In this problem, it is taken as 0.04m. The partial system of equation is solved using following discritizations:

#### **Forward Time Central Space (FTCS) Method**

It is a first-order in time and second-order in space. The discritization for this partial equation:

$$u_{i}^{n+1}=u_{i}^{n}+d(u_{i+1}^{n}-2u_{i}^{n}+u_{i-1}^{n})$$

where

$$d=\frac{α(Δt)}{Δx^{2}}$$

For the stability condition:

$$\frac{α(Δt)}{Δx^{2}}=d\leqslant0.5$$

#### **DuFort-Frankel Method**

It is a second-order in time and second-order in space. The discritization for this partial equation:

$$(1+2d)u_{i}^{n+1}=(1-2d)u_{i}^{n-1}+2d(u_{i+1}^{n}+u_{i-1}^{n})$$

where

$$d=\frac{α(Δt)}{Δx^{2}}$$

For the stability condition:

$$\frac{α(Δt)}{Δx^{2}}\leqslant0.5$$

This method requires 2 sets of initial data so in this notebook, FTCS is used to obtain data set at t=Δt.

#### **Laasonen Method**

It is a first-order in time and second-order in space. The discritization for this partial equation:

$$du_{i-1}^{n+1}-(1+2d)u_{i}^{n+1}+du_{i+1}^{n+1}=-u_{i}^{n}$$

where

$$d=\frac{α(Δt)}{Δx^{2}}$$

This method is unconditionally stable.

#### **Crank Nicholson Method**

It is a second-order in time and second-order in space. The discritization for this partial equation:

$$-\frac{1}{2}du_{i-1}^{n+1}-(1+d)u_{i}^{n+1}-\frac{1}{2}du_{i+1}^{n+1}=\frac{1}{2}du_{i-1}^{n}+(1-d)u_{i}^{n}+\frac{1}{2}du_{i+1}^{n}$$

where

$$d=\frac{α(Δt)}{Δx^{2}}$$

This method is unconditionally stable.

Details of schemes:
Computational Fluid Dynamics, Volume I, Hoffmann - Chiang

Feel free to ask questions!

*Furkan Oz*
foz@okstate.edu
"""

# ╔═╡ 29683fb0-6419-11eb-05d0-2d7392bc8afa
"""
Computes the Forward Time Central Space for the Parabolic Partial Differential Equation
This notebook computes the Parabolic Partial Differential Equation
		FTCS(Δt,Δy,ν)
If the arguments are missing, it will use the default values.

		FTCS(Δt=0.002,Δy=0.001,ν=0.000217)
Furkan Oz,
foz@okstate.edu,

"""
function FTCS(Δt=0.002,Δy=0.001,ν=0.000217)
	#ν=0.000217 # viscosity [m²/s]
	h = 0.04 # Distance between 2 plate [m]
	T = 1.25 # Total time [s]
	#Δt = 0.002 # Time increment
	Nt = T/Δt+1 # Number of time step
	Nt = round(Int32,Nt)

	#Δy = 0.001 # Spacing between node points
	y = 0:Δy:h # y-coordinates
	Nx = h/Δy+1 # Number of Node points
	Nx = round(Int32,Nx)

	U₀ = 1 # Plate speed
	d = ν*Δt/(Δy)^2
	#Initializing
	u = zeros(Nx,Nt)

	# Boundary Conditions
	u[1,:] = [U₀ for i=1:Nt]
	u[Nx,:] = [0 for i=1:Nt]

	# Initial Conditions
	u[:,1] = [0 for i=1:Nx]
	u[1,1] = U₀
	for ii=1:Nt-1
		for i=2:Nx-1
			u[i,ii+1] = u[i,ii]+d*(u[i+1,ii]-2*u[i,ii]+u[i-1,ii])
		end
	end

	return u,y

end

# ╔═╡ 36be1c60-6442-11eb-204a-b37677c1acfd
"""
Computes the DuFort Frankel method for the Parabolic Partial Differential Equation
This notebook computes the Parabolic Partial Differential Equation
		DuFortFrankel(Δt,Δy,ν)
If the arguments are missing, it will use the default values.

		DuFortFrankel(Δt=0.002,Δy=0.001,ν=0.000217)
Furkan Oz,
foz@okstate.edu,

"""
function DuFortFrankel(Δt=0.002,Δy=0.001,ν=0.000217)
	#ν=0.000217 # viscosity [m²/s]
	h = 0.04 # Distance between 2 plate [m]
	T = 1.25 # Total time [s]
	#Δt = 0.002 # Time increment
	Nt = T/Δt+1 # Number of time step
	Nt = round(Int32,Nt)

	#Δy = 0.001 # Spacing between node points
	y = 0:Δy:h # y-coordinates
	Nx = h/Δy+1 # Number of Node points
	Nx = round(Int32,Nx)

	U₀ = 1 # Plate speed
	d = ν*Δt/(Δy)^2
	#Initializing
	u = zeros(Nx,Nt)

	# Boundary Conditions
	u[1,:] = [U₀ for i=1:Nt]
	u[Nx,:] = [0 for i=1:Nt]

	# Initial Conditions
	u[:,1] = [0 for i=1:Nx]
	u[1,1] = U₀
	for ii=1:1
		for i=2:Nx-1
			u[i,ii+1] = u[i,ii]+d*(u[i+1,ii]-2*u[i,ii]+u[i-1,ii])
		end
	end
	for ii=2:Nt-1
		for i=2:Nx-1
			u[i,ii+1] = ((1-2*d)*u[i,ii-1]+2*d*(u[i+1,ii]+u[i-1,ii]))/(1+2*d)
		end
	end
	return u,y

end

# ╔═╡ 899d6a70-6448-11eb-2870-692c6e93c421
"""
Computes the Laasonen method for the Parabolic Partial Differential Equation
This notebook computes the Parabolic Partial Differential Equation
		Laasonen(Δt,Δy,ν)
If the arguments are missing, it will use the default values.

		Laasonen(Δt=0.002,Δy=0.001,ν=0.000217)
Furkan Oz,
foz@okstate.edu,

"""
function Laasonen(Δt=0.002,Δy=0.001,ν=0.000217)
	#ν=0.000217 # viscosity [m²/s]
	h = 0.04 # Distance between 2 plate [m]
	T = 1.25 # Total time [s]
	#Δt = 0.002 # Time increment
	Nt = T/Δt+1 # Number of time step
	Nt = round(Int32,Nt)

	#Δy = 0.001 # Spacing between node points
	y = 0:Δy:h # y-coordinates
	Nx = h/Δy+1 # Number of Node points
	Nx = round(Int32,Nx)

	U₀ = 1 # Plate speed
	d = ν*Δt/(Δy)^2

	#Initializing
	A = zeros(Nx,Nx)
	b = zeros(Nx)
	u = zeros(Nx,Nt)

	# Boundary Conditions
	A[1,1] = 1
	b[1] = U₀
	A[Nx,Nx] = 1
	b[Nx] = 0

	# Initial Conditions
	u[:,1] = [0 for i=1:Nx]
	u[1,1] = U₀

	for ii=2:Nt-1
		# Boundary Conditions
		A[1,1] = 1
		b[1] = U₀
		A[Nx,Nx] = 1
		b[Nx] = 0
		for i=2:Nx-1
			A[i,i-1] = d
			A[i,i+0] =-1-2*d
			A[i,i+1] = d
			b[i] =-u[i,ii]
		end
		ulin = A\b
		u[1,ii+1] = b[1]
		u[Nx,ii+1] = b[Nx]
		u[2:Nx-1,ii+1] = [ulin[i] for i=2:Nx-1]
	end
	return u,y
end

# ╔═╡ 9c410110-644c-11eb-3bd9-93935242ac98
"""
Computes the Crank Nicolson method for the Parabolic Partial Differential Equation
This notebook computes the Parabolic Partial Differential Equation
		CrankNicolson(Δt,Δy,ν)
If the arguments are missing, it will use the default values.

		CrankNicolson(Δt=0.002,Δy=0.001,ν=0.000217)
Furkan Oz,
foz@okstate.edu,

"""
function CrankNicolson(Δt=0.002,Δy=0.001,ν=0.000217)
	#ν=0.000217 # viscosity [m²/s]
	h = 0.04 # Distance between 2 plate [m]
	T = 1.25 # Total time [s]
	#Δt = 0.002 # Time increment
	Nt = T/Δt+1 # Number of time step
	Nt = round(Int32,Nt)

	#Δy = 0.001 # Spacing between node points
	y = 0:Δy:h # y-coordinates
	Nx = h/Δy+1 # Number of Node points
	Nx = round(Int32,Nx)

	U₀ = 1 # Plate speed
	d = ν*Δt/(Δy)^2

	#Initializing
	A = zeros(Nx,Nx)
	b = zeros(Nx)
	u = zeros(Nx,Nt)

	# Boundary Conditions
	A[1,1] = 1
	b[1] = U₀
	A[Nx,Nx] = 1
	b[Nx] = 0

	# Initial Conditions
	u[:,1] = [0 for i=1:Nx]
	u[1,1] = U₀

	for ii=3:Nt-1
		# Boundary Conditions
		A[1,1] = 1
		b[1] = U₀
		A[Nx,Nx] = 1
		b[Nx] = 0
		for i=2:Nx-1
			A[i,i-1] =-0.5*d
			A[i,i+0] = 1+d
			A[i,i+1] =-0.5*d
			b[i] = 0.5*d*u[i+1,ii]+(1-d)*u[i,ii]+0.5*d*u[i-1,ii]
		end
		ulin = A\b
		u[1,ii+1] = b[1]
		u[Nx,ii+1] = b[Nx]
		u[2:Nx-1,ii+1] = [ulin[i] for i=2:Nx-1]
	end
	return u,y
end

# ╔═╡ 8939c23e-642a-11eb-2de6-ade5ae2f4b4d
function Analytical(t,Δy=0.001,ν=0.000217)
	#ν=0.000217 # viscosity [m²/s]
	h = 0.04 # Distance between 2 plate [m]
	y = 0:Δy:h
	Nx = h/Δy+1 # Number of Node points
	Nx = round(Int32,Nx)
	U₀ = 1 # Plate speed
	η = [y[i]/(2*sqrt(ν*t)) for i=1:Nx]
	η₁= h/(2*sqrt(ν*t))
	u = [U₀*(erfc(η[i])-erfc(2*η₁-η[i])+erfc(2*η₁+η[i])-erfc(4*η₁-η[i])+erfc(4*η₁+η[i])-erfc(8*η₁-η[i])+erfc(8*η₁+η[i])-erfc(16*η₁-η[i])) for i=1:Nx]
	return u,y
end

# ╔═╡ 70d8c3ee-641f-11eb-1d57-a943a361c397
begin
	uFTCS,yFTCS = FTCS()
	plot([uFTCS[:,1],uFTCS[:,90],uFTCS[:,180],uFTCS[:,270],uFTCS[:,360],uFTCS[:,450],uFTCS[:,540]],yFTCS,
			title = "FTCS - \\Deltat = 0.002 - \\Deltay = 0.001",
			label = ["t=0s" "t=0.18s" "t=0.36s" "t=0.54s" "t=0.72s" "t=0.90s" "t=1.08s"],
			legend = :topright,
			ylims = (0,0.04),
			xlabel = "U",
			ylabel = "y",
			linewidth = 2,
			linecolor = :black,
			markershape = :auto,
			markercolor = :red,
			)
end

# ╔═╡ 95697270-6427-11eb-02be-bf3b75fbbd33
begin
	u2,y2 = FTCS(0.005,0.001)

	plot([u2[:,1],u2[:,50],u2[:,100],u2[:,150],u2[:,200]],y2,
			title = "FTCS - \\Deltat = 0.005 - \\Deltay = 0.001",
			label = ["t=0s" "t=0.25s" "t=0.50s" "t=0.75s" "t=1.00s"],
			legend = :topright,
			ylims = (0,0.04),
			xlabel = "U",
			ylabel = "y",
			linewidth = 2,
			linecolor = :black,
			markershape = :auto,
			markercolor = :red,
			)
end

# ╔═╡ b1e5b270-643f-11eb-1a37-e3227c598177
begin
	ua1,ya1 = Analytical(0.18,0.001)
	uDF,yDF = DuFortFrankel()
	uLaasonen,yLaasonen = Laasonen()
	uCrankNicolson,yCrankNicolson = CrankNicolson()
	plot([ua1-uFTCS[:,90],ua1-uDF[:,90],ua1-uLaasonen[:,90],ua1-uCrankNicolson[:,90]],[ya1],
				title = "Error - t = 0.18 - \\Deltay = 0.001",
				label = ["FTCS" "D/F" "Laasonen" "CrankNicolson"],
				legend = :topright,
				xlims = (-0.005,0.01),
				ylims = (0,0.04),
				xlabel = "U",
				ylabel = "y",
				linewidth = 2,
				linecolor = :black,
				markershape = :auto,
				markercolor = :red,
				)
end

# ╔═╡ 8c11e680-6440-11eb-076a-0b1e169ef37a
begin
	ua2,ya2 = Analytical(1.08,0.001)
	plot([ua2-uFTCS[:,540],ua2-uDF[:,540],ua2-uLaasonen[:,540],ua2-uCrankNicolson[:,540]],[ya2],
				title = "Error - t = 1.08 -\\Deltay = 0.001",
				label = ["FTCS" "D/F" "Laasonen" "CrankNicolson"],
				legend = :topright,
				xlims = (-0.0005,0.0015),
				ylims = (0,0.04),
				xlabel = "U",
				ylabel = "y",
				linewidth = 2,
				linecolor = :black,
				markershape = :auto,
				markercolor = :red,
				)
end

# ╔═╡ 33a79720-644e-11eb-0eca-e794b9b87ccf
begin
	ua3,ya3 = Analytical(1.00,0.001)
	uLaasonen005,yLaasonen005 = Laasonen(0.005)
	uLaasonen01,yLaasonen01 = Laasonen(0.01)
	uLaasonen1,yLaasonen1 = Laasonen(0.1)
	uLaasonen2,yLaasonen2 = Laasonen(0.2)
	plot([ua3-uLaasonen005[:,200],ua3-uLaasonen01[:,100],ua3-uLaasonen1[:,10],ua3-uLaasonen2[:,5]],[ya2],
					title = "Error - Laasonen - t = 1.00 -\\Deltay = 0.001",
					label = ["\\Deltat=0.005" "\\Deltat=0.01" "\\Deltat=0.1" "\\Deltat=0.2"],
					legend = :topright,
					xlims = (0,0.2),
					ylims = (0,0.04),
					xlabel = "U",
					ylabel = "y",
					linewidth = 2,
					linecolor = :black,
					markershape = :auto,
					markercolor = :red,
					)
end

# ╔═╡ Cell order:
# ╟─0b868b30-6420-11eb-18f9-1592f7542fe6
# ╟─4c4fff40-64c8-11eb-2c45-57f4793ee997
# ╟─56d1f810-64c8-11eb-3e0d-81c0bbafdc62
# ╟─29683fb0-6419-11eb-05d0-2d7392bc8afa
# ╟─36be1c60-6442-11eb-204a-b37677c1acfd
# ╟─899d6a70-6448-11eb-2870-692c6e93c421
# ╟─9c410110-644c-11eb-3bd9-93935242ac98
# ╟─8939c23e-642a-11eb-2de6-ade5ae2f4b4d
# ╟─70d8c3ee-641f-11eb-1d57-a943a361c397
# ╟─95697270-6427-11eb-02be-bf3b75fbbd33
# ╟─b1e5b270-643f-11eb-1a37-e3227c598177
# ╟─8c11e680-6440-11eb-076a-0b1e169ef37a
# ╟─33a79720-644e-11eb-0eca-e794b9b87ccf


text:
"""

```
centerofcoordinates(coor::AbstractVector{T}) where T
```

Computes the center of coordinates of a vector.

"""
function centerofcoordinates(coor::AbstractVector{T}) where T
  cm = zeros(T)
  n = length(coor)
  for i in 1:n
    cm = cm + coor[i]
  end
  cm = cm / n
  return cm
end


text:
# Forked from https://github.com/mohamed82008/InplaceQR.jl.git
module InplaceQR
include("QR/qr_fact.jl")
include("QR/qr_api.jl")
end # module


text:
using LinearForms

export call, @mm

import Base.call
function call(u::HilbertVector, f, params...)
    HilbertVector{field(u)}(u.idx, u.space, [(f, params...); u.opstack])
end

macro mm(xp)
    transposecalls!(xp)
end


text:
using ADatasets
using Iterators
using StatsBase
using Iterators
using Lazy
using DataFrames
using IPMeasures
using SemiSupervised
using FluxExtensions
using Flux
using kNN
using MLDataPattern
import ADatasets: runtest
import FluxExtensions: layerbuilder, freeze

idir = filter(isdir,["/Users/tpevny/Work/Data/datasets/numerical","/mnt/output/data/datasets/numerical","/opt/output/data/datasets/numerical"])[1];
odir = filter(isdir,["/Users/tpevny/Work/Julia/results/datasets","/mnt/output/results/datasets","/opt/output/results/datasets"])[1];

function fit(data,layers,hidden,zdim,β)
  idim = size(data.data[1], dims = 1)
  σ2 = SemiSupervised.bandwidthofkde(randn(zdim,1000),20)
	m = SemiSupervised.VAE(FluxExtensions.layerbuilder(idim,hidden,2*zdim,3,"relu","linear","Dense"),
		FluxExtensions.layerbuilder(zdim,hidden,idim,3,"relu","linear","Dense"), β)
  m = to32(m);
  mf = freeze(m)
  opt = Flux.Optimise.ADAM(params(m));

  FluxExtensions.learn((x) -> SemiSupervised.loss(m,getobs(x)),opt,data,() -> (),10000)

  # let's estimate the distance between the prior and posterior using MMD (zmmd)
  # and calibration distance on prior and prior
  x = IPMeasures.samplecolumns(data.data[1],1000)
  n = size(x, dims = 2)
  z = SemiSupervised.gaussiansample(SemiSupervised.hsplitsoftp(mf.q(x))...)
  kσ2 = SemiSupervised.bandwidthofkde(z,20)
	kde = x -> SemiSupervised.kde(x,z,kσ2)
	df = DataFrame(
		layers = layers,
		hidden = hidden,
		zdim = zdim,
		β = β,
		zmmd = IPMeasures.mmdg(z, randn(zdim,n), 0.5/σ2),
		ummd = IPMeasures.mmdg(randn(zdim,n), randn(zdim,n), 0.5/σ2),
		recerror =  Flux.mse(x,mf.g(z)) / n,
		likelihood =  - mean(SemiSupervised.log_normal(x, mf.g(z), β)),
		train_px = mean(SemiSupervised.px(mf,x, 100,mf.β)),
		train_pxis = mean(SemiSupervised.pxis(mf, x , 100, mf.β)),
		train_kde_pxis = mean(SemiSupervised.pxis(mf, x, kde, 100, mf.β)),
		train_pxvita = mean(SemiSupervised.pxvita(mf,x,mf.β))
		)
	(m, kde), df
end

predicts = [(m, x) -> - SemiSupervised.px(m[1],x,100,m[1].β)[:],
	(m, x) -> - SemiSupervised.pxis(m[1],x,100, m[1].β)[:],
	(m, x) -> - SemiSupervised.pxis(m[1], x, m[2], 100, m[1].β)[:],
	(m, x) -> - SemiSupervised.pxvita(m[1],x,m[1].β)[:]]
prnames = Symbol.(["px","pxis","kde_pxis","vita"])

layers, hidden, zdim, β = 3, 32, 32, 2.0

train, test, clusterdness = ADatasets.makeset(ADatasets.loaddataset("breast-cancer-wisconsin","easy",idir)..., 0.75,"low")

subspace(x,ft) = (x[1][ft,:],x[2])
d = size(train[1], dims = 1)
df = mapreduce(vcat, [[i,j] for i in 2:d for j in 1:i-1]) do ft
	trn, tst = subspace(train,ft), subspace(test,ft)
	data = RandomBatches((trn[1],),100,10000)
	m,info = fit(data, layers, hidden, zdim, β)
	mf = FluxExtensions.freeze(m)
	df = DataFrame(ADatasets.evaluate(predicts, mf, tst)')
	names!(df,prnames)

	model = kNN.KNNAnomaly(trn[1],:kappa)
	ks = [1, 3, 5, 9, 13, 17, 21]
	aucs = map(k -> ADatasets.evaluate(x -> kNN.predict(model,x,k),tst), ks)
	i = argmax(aucs)
	df[:knn] = aucs[i]
	df[:k] = ks[i]
	df[:i] = ft[1]
	df[:j] = ft[2]
	display(df)
	df
end

function rankvec(x)
	y = zeros(Int,size(x))
	y[sortperm(x, rev = true)] .= 1:length(y)
	y ./ length(y)
end

ft = [1,2]
trn, tst = subspace(train,ft), subspace(test,ft)
data = RandomBatches((trn[1],),100,10000)
m,info = fit(data, layers, hidden, zdim, β)
mf = FluxExtensions.freeze(m)
z1 = SemiSupervised.pxis(mf[1], tst[1], mf[2], 100, m[1].β)[:]
model = kNN.KNNAnomaly(trn[1],:kappa)
z2 = kNN.predict(model,tst[1],21)
scatter(tst[1][1,:],tst[1][2,:],zcolor = rankvec(z1), marker = 3*tst[2])
scatter(tst[1][1,:],tst[1][2,:],zcolor = rankvec(z2), marker = 3*tst[2])


text:
using LightGraphs

import LightGraphs: src, dst

struct ElementaryPathInstance{T}
  graph::AbstractGraph{T}
  costs::Dict{Edge{T}, Float64}
  src::T
  dst::T
end

graph(i::ElementaryPathInstance{T}) where T = i.graph
src(i::ElementaryPathInstance{T}) where T = i.src
dst(i::ElementaryPathInstance{T}) where T = i.dst

# dimension(i::ElementaryPathInstance{T}) where T = ne(graph(i))
cost(i::ElementaryPathInstance{T}, u::T, v::T) where T = i.costs[Edge(u, v)]

struct ElementaryPathSolution{T}
  instance::ElementaryPathInstance{T}
  path::Vector{Edge{T}}
  states::Dict{T, Float64}
  solutions::Dict{T, Vector{Edge{T}}}
end

function lp_dp(i::ElementaryPathInstance{T}) where T # I.e. Bellman-Ford algorithm. Assumption: no positive-cost cycle in the graph.
  V = Dict{T, Float64}()
  S = Dict{T, Vector{Edge{T}}}()

  # Initialise.
  for v in vertices(graph(i))
    V[v] = -Inf
    S[v] = Edge{T}[]
  end
  V[src(i)] = 0.0

  # Dynamic part.
  for _ in 1:ne(graph(i))
    changes = false # Stop the algorithm as soon as it no more makes progress.

    for e in edges(graph(i))
      u, v = src(e), dst(e)
      w = cost(i, u, v)

      # If using the solution to the currently explored subproblem would
      # lead to a cycle, skip it.
      if any(src(e) == v for e in S[u])
        continue
      end

      # Compute the maximum: is passing through w advantageous?
      if V[u] + w > V[v]
        V[v] = V[u] + w
        S[v] = vcat(S[u], Edge(u, v))

        changes = true
      end
    end

    if ! changes
      break
    end
  end

  # Checking existence of negative-cost cycle.
  for e in edges(graph(i))
    u, v = src(e), dst(e)
    w = cost(i, u, v)

    if V[u] + w > V[v]
      @warn("The graph contains a positive-cost cycle around edge $(u) -> $(v).")
    end
  end

  return ElementaryPathSolution(i, S[dst(i)], V, S)
end


text:
module ACSetInterface
export ACSet, acset_schema,
  nparts, parts, has_part, has_subpart, subpart, incident,
  add_part!, add_parts!, set_subpart!, set_subparts!, rem_part!, rem_parts!,
  copy_parts!, copy_parts_only!, disjoint_union, tables, pretty_tables, @acset

using StaticArrays: StaticArray

using ..Syntax: GATExpr, args
using ..Theories: dom, codom
using ...Meta: strip_lines

using PrettyTables: pretty_table
using Tables

abstract type ACSet end

"""
Get the schema of an acset at runtime.
"""
function acset_schema end

""" Number of parts of given type in an acset.
"""
function nparts end

""" Parts of given type in an acset.
"""
@inline parts(acs, type) = 1:nparts(acs, type)

""" Whether an acset has a part with the given name.
"""
function has_part end

@inline has_part(acs, type::Symbol, part::Int) = 1 <= part <= nparts(acs, type)
@inline has_part(acs, type::Symbol, part::AbstractVector{Int}) =
  let n=nparts(acs, type); [ 1 <= x <= n for x in part ] end

function has_subpart end

""" Get subpart of part in acset.

Both single and vectorized access are supported, with a view of the underlying
data being returned in the latter case. Chaining, or composition, of subparts is
also supported. For example, given a vertex-attributed graph `g`,

```
subpart(g, e, [:src, :vattr])
```

returns the vertex attribute of the source vertex of the edge `e`. As a
shorthand, subparts can also be accessed by indexing:

```
g[e, :src] == subpart(g, e, :src)
```

Be warned that indexing with lists of subparts works just like `subpart`:
`g[e,[:src,:vattr]]` is equivalent to `subpart(g, e, [:src,:vattr])`. This
convention differs from DataFrames but note that the alternative interpretation
of `[:src,:vattr]` as two independent columns does not even make sense, since
they have different domains (belong to different tables).
"""
function  subpart end
@inline Base.@propagate_inbounds subpart(acs, part, name) = view_or_slice(subpart(acs, name), part)

function view_or_slice end
@inline view_or_slice(x::AbstractVector, i::Union{Integer,StaticArray}) = x[i]
@inline view_or_slice(x::AbstractVector, ::Colon) = x
@inline Base.@propagate_inbounds view_or_slice(x::AbstractVector, i) = @view x[i]

@inline subpart(acs, expr::GATExpr{:generator}) = subpart(acs, first(expr))
@inline subpart(acs, expr::GATExpr{:id}) = parts(acs, first(dom(expr)))

function subpart(acs, part, names::AbstractVector{Symbol})
  foldl(names, init=part) do part, name
    subpart(acs, part, name)
  end
end
subpart(acs, part, expr::GATExpr{:compose}) =
  subpart(acs, part, subpart_names(expr))

subpart(acs, names::AbstractVector{Symbol}) =
  subpart(acs, subpart(acs, names[1]), names[2:end])
subpart(acs, expr::GATExpr{:compose}) = subpart(acs, subpart_names(expr))

subpart_names(expr::GATExpr{:generator}) = Symbol[first(expr)]
subpart_names(expr::GATExpr{:id}) = Symbol[]
subpart_names(expr::GATExpr{:compose}) =
  mapreduce(subpart_names, vcat, args(expr))

""" Get superparts incident to part in acset.

If the subpart is indexed, this takes constant time; otherwise, it takes linear
time. As with [`subpart`](@ref), both single and vectorized access, as well as
chained access, are supported. Note that sequences of morphisms are supplied in
the usual left-to-right order, so that

```
incident(g, x, [:src, :vattr])
```

returns the list of all edges whose source vertex has vertex attribute `x`.

Note that when the subpart is indexed, this function returns a view of the
underlying index, which should not be mutated. To ensure that a fresh copy is
returned, regardless of whether indexing is enabled, set the keyword argument
`copy=true`.
"""
function incident end

function incident(acs, part, names::AbstractVector{Symbol};
                  copy::Bool=false)
  # Don't need to pass `copy` because copy will be made regardless.
  foldr(names, init=part) do name, part
    reduce(vcat, incident(acs, part, name), init=Int[])
  end
end

incident(acs, part, expr::GATExpr; kw...) =
  incident(acs, part, subpart_names(expr); kw...)

@inline add_part!(acs, type; kw...) = add_part!(acs, type, (;kw...))

""" Add part of given type to acset, optionally setting its subparts.

Returns the ID of the added part.

See also: [`add_parts!`](@ref).
"""
@inline function add_part!(acs, type::Symbol, kw)
  part = only(add_parts!(acs,type,1))
  try
    set_subparts!(acs, part, kw)
  catch e
    rem_part!(acs, type, part)
    rethrow(e)
  end
  part
end

""" Add parts of given type to acset, optionally setting their subparts.

Returns the range of IDs for the added parts.

See also: [`add_part!`](@ref).
"""
function add_parts! end

@inline add_parts!(acs, type::Symbol, n::Int; kw...) = add_parts!(acs, type, n, (;kw...))

@inline function add_parts!(acs, type::Symbol, n::Int, kw)
  parts = add_parts!(acs, type, n)
  try
    set_subparts!(acs, parts, kw)
  catch e
    rem_parts!(acs, type, parts)
    rethrow(e)
  end
  parts
end

""" Mutate subpart of a part in a C-set.

Both single and vectorized assignment are supported.

See also: [`set_subparts!`](@ref).
"""
function set_subpart! end
@inline set_subpart!(acs, name, vals) = set_subpart!(acs, :, name, vals)

@inline set_subpart!(acs, ::Colon, name, vals) =
  set_subpart!(acs, 1:length(subpart(acs,name)), name, vals)

# Inlined for the same reason as `subpart`.

@inline function set_subpart!(acs, part::AbstractVector{Int}, name, vals)
  broadcast(part, vals) do part, vals
    set_subpart!(acs, part, name, vals)
  end
end

""" Mutate subparts of a part in a C-set.

Both single and vectorized assignment are supported.

See also: [`set_subpart!`](@ref).
"""
@inline @generated function set_subparts!(acs, part, kw::NamedTuple{keys}) where {keys}
  Expr(:block,[:(set_subpart!(acs, part, $(Expr(:quote, name)), kw.$name)) for name in keys]...)
end

@inline set_subparts!(acs, part; kw...) = set_subparts!(acs, part, (;kw...))

""" Remove part from a C-set.

The part is removed using the "pop and swap" strategy familiar from
[Graphs.jl](https://github.com/JuliaGraphs/Graphs.jl), where the "removed" part
is actually replaced by the last part, which is then deleted. This strategy has
important performance benefits since only the last part must be assigned a new
ID, as opposed to assigning new IDs to *every* part following the removed part.

The removal operation is *not* recursive. When a part is deleted, any superparts
incident to it are retained, but their subparts become undefined (equal to the
integer zero). For example, in a graph, if you call `rem_part!` on a vertex, the
edges incident the `src` and/or `tgt` vertices of the edge become undefined but
the edge itself is not deleted.

Indexing has both positive and negative impacts on performance. On the one hand,
indexing reduces the cost of finding affected superparts from linear time to
constant time. On the other hand, the indices of subparts must be updated when
the part is removed. For example, in a graph, indexing `src` and `tgt` makes
removing vertices faster but removing edges (slightly) slower.

See also: [`rem_parts!`](@ref).
"""
function rem_part! end

""" Remove parts from a C-set.

The parts must be supplied in sorted order, without duplicates.

See also: [`rem_part!`](@ref).
"""
@inline function rem_parts!(acs, type, parts)
  issorted(parts) || error("Parts to be removed must be in sorted order")
  for part in Iterators.reverse(parts)
    rem_part!(acs, type, part)
  end
end

""" Copy parts from a C-set to a C′-set.

The selected parts must belong to both schemas. All subparts common to the
selected parts, including data attributes, are preserved. Thus, if the selected
parts form a sub-C-set, then the whole sub-C-set is preserved. On the other
hand, if the selected parts do *not* form a sub-C-set, then some copied parts
will have undefined subparts.

TODO: handle colons
"""
function copy_parts! end

copy_parts!(to, from, obs::Tuple) =
  copy_parts!(to, from, NamedTuple{obs}((:) for ob in obs))

copy_parts!(to, from; kw...) = copy_parts!(to, from, (;kw...))

""" Copy parts from a C-set to a C′-set, ignoring all non-data subparts.

The selected parts must belong to both schemas. Attributes common to both
schemas are also copied, but no other subparts are copied.

See also: [`copy_parts!`](@ref).
"""
function copy_parts_only! end

function disjoint_union(acs1, acs2)
  acs = copy(acs1)
  copy_parts!(acs, acs2)
  acs
end

"""
Get a named tuple of Tables.jl-compatible tables from an acset
"""
function tables end

""" Display an ACSet prettily

This works for any acset that implements tables
"""

# Fancy Displaying

function pretty_tables(io::IO, acs::ACSet; kw...)
  options = merge((nosubheader=true, show_row_number=true), (; kw...))
  for (ob, table) in pairs(tables(acs))
    # Note: PrettyTables will not print tables with zero rows.
    if !(isempty(Tables.columnnames(table)) || Tables.rowcount(table) == 0)
      pretty_table(io, table, row_number_column_title=string(ob); options...)
    end
  end
end

pretty_tables(acs::ACSet; kw...) = pretty_tables(stdout, acs; kw...)

function Base.show(io::IO, ::MIME"text/plain", acs::T) where T <: ACSet
  print(io, T)
  print(io, " with elements ")
  join(io, ["$ob = $(parts(acs,ob))" for ob in keys(tables(acs))], ", ")
  println(io)
  pretty_tables(io, acs)
end

function Base.show(io::IO, ::MIME"text/html", acs::T) where T <: ACSet
  println(io, "
")
  print(io, "")
  print(io, T)
  print(io, " with elements ")
  join(io, ["$ob = $(parts(acs,ob))" for ob in keys(tables(acs))], ", ")
  println(io, "")
  pretty_tables(io, acs, backend=Val(:html), standalone=false)
  println(io, "
")
end

macro acset(head, body)
  @assert body.head == :block
  vals = Expr(:call, :(Dict{Symbol,Any}))
  for l in strip_lines(body).args
    @assert l.head == :(=)
    push!(vals.args, :($(Expr(:quote, l.args[1])) => $(l.args[2])))
  end
  :(init_acset($(esc(head)), $(esc(vals))))
end

function init_acset(T::Type{<:ACSet}, initvals::Dict{Symbol,Any})
  acs = T()
  s = acset_schema(acs)
  ob_specs = filter((kv) -> kv[1] ∈ s.obs, pairs(initvals))
  hom_specs = filter((kv) -> kv[1] ∈ s.homs, pairs(initvals))
  attr_specs = filter((kv) -> kv[1] ∈ s.attrs, pairs(initvals))
  for (k,v) in ob_specs
    add_parts!(acs, k, Int(v))
  end
  for (k,v) in Iterators.flatten((hom_specs, attr_specs))
    set_subpart!(acs, :, k, collect_nonvector(v))
  end
  acs
end
collect_nonvector(v::AbstractVector) = v
collect_nonvector(v) = collect(v)

end

max_forks_repo_name: epatters/CompCat


text:
abstract type Gate end

struct QulacsGate <: Gate
    pyobj::PyObject
end

function RX(target::Int, angle::Float64)
    QulacsGate(qulacs.gate.RX(target, angle))
end

function RY(target::Int, angle::Float64)
    QulacsGate(qulacs.gate.RY(target, angle))
end

function RZ(target::Int, angle::Float64)
    QulacsGate(qulacs.gate.RZ(target, angle))
end

function get_matrix(gate::QulacsGate)
    gate.pyobj.get_matrix()
end

max_stars_repo_name: Sakurairihito/QCMaterialNew

max_issues_repo_name: sakurairihito/QCMaterialNew

max_forks_repo_name: sakurairihito/QCMaterialNew


text:
using Chain
using Utils
using OffsetArrays

dirs = Dict(
    'D' => (1, 0),
    'U' => (-1, 0),
    'L' => (0, -1),
    'R' => (0, 1),
)

const Num = Int32

function parse_dir(item)
    d = item[1]
    parsed = parse(Num, item[2:end])
    CartesianIndex(dirs[d] .* parsed)
end

Input = @NamedTuple{wires::Array{Array{CartesianIndex,1},1}, min::CartesianIndex{2}, max::CartesianIndex{2}}

inputfile = "$(replace(@__FILE__, r".jl$" => "")).txt"
function load(text)::Input
    mi = CartesianIndex(0, 0)
    ma = CartesianIndex(0, 0)
    re = Array{CartesianIndex,1}[]
    for line in per_line(text, false)
        pos = CartesianIndex(0, 0)
        inner = CartesianIndex[]
        for item in per_split(line, ',')
            parsed = parse_dir(item)
            pos += parsed
            mi = min(mi, pos)
            ma = max(ma, pos)
            push!(inner, parsed)
        end
        push!(re, inner)
    end
    Input((wires=re, min=mi, max=ma))
end
fload() = inputfile |> slurp |> load
input = fload()

manhatten(input) = mapreduce(j -> abs(input[j]), +, (1,2))

function drawlines(input::Input = fload())
    mi::CartesianIndex = input.min
    rows = input.max[1] - mi[1] + 1
    cols = input.max[2] - mi[2] + 1

    map(input.wires) do wire
        steps = OffsetArray(zeros(Num, rows, cols), mi[1] - 1, mi[2] - 1)
        start = CartesianIndex(0, 0)
        increment = 0

        for move in wire
            stop = start + move
            if start[1] == stop[1]
                same = 1
                diff =  2
            else
                same = 2
                diff =  1
            end
            step = start[diff] < stop[diff] ? 1 : -1
            for i in start[diff]:step:stop[diff]-step
                cursor = CartesianIndex(same == 1 ? start[1] : i, same == 2 ? start[2] : i)
                getindex(steps, cursor) > 0 || setindex!(steps, getindex(steps, cursor) + increment, cursor)
                increment += 1
            end
            start = stop
        end

        crossings = steps .> 0
        (; crossings, steps)
    end
end

function part1(input::Input = fload())
    @chain input begin
        drawlines
        mapreduce(i -> i.crossings, (a,b) -> a .& b, _)
        findall
        minimum(manhatten, _)
    end
end

function part2(input::Input = fload())
    lines = drawlines(input)
    crossings = mapreduce(i -> i.crossings, (a,b) -> a .& b, lines)
    mapreduce(i -> i.steps[crossings], +, lines) |> minimum
end

println("Closest intersection to the origin: ", part1())
println("Earliest intersection after the origin: ", part2())

max_stars_repo_name: CmdQ/AoC


text:
using ModelingToolkit
using Test

# Derivatives
@parameters t σ ρ β
@variables x(t) y(t) z(t)
@derivatives D'~t D2''~t

@test isequal(expand_derivatives(D(t)), 1)
@test isequal(expand_derivatives(D(D(t))), 0)

dsin = D(sin(t))
@test isequal(expand_derivatives(dsin), cos(t))

dcsch = D(csch(t))
@test isequal(expand_derivatives(dcsch), simplify_constants(coth(t) * csch(t) * -1))

@test isequal(expand_derivatives(D(-7)), 0)
@test isequal(expand_derivatives(D(sin(2t))), simplify_constants(cos(2t) * 2))
@test isequal(expand_derivatives(D2(sin(t))), simplify_constants(-sin(t)))
@test isequal(expand_derivatives(D2(sin(2t))), simplify_constants(sin(2t) * -4))
@test isequal(expand_derivatives(D2(t)), 0)
@test isequal(expand_derivatives(D2(5)), 0)

# Chain rule
dsinsin = D(sin(sin(t)))
@test isequal(expand_derivatives(dsinsin), cos(sin(t))*cos(t))

d1 = D(sin(t)*t)
d2 = D(sin(t)*cos(t))
@test isequal(expand_derivatives(d1), t*cos(t)+sin(t))
@test isequal(expand_derivatives(d2), simplify_constants(cos(t)*cos(t)+sin(t)*(-1*sin(t))))

eqs = [0 ~ σ*(y-x),
       0 ~ x*(ρ-z)-y,
       0 ~ x*y - β*z]
sys = NonlinearSystem(eqs,[x,y,z],[σ,ρ,β])
jac = calculate_jacobian(sys)
@test isequal(jac[1,1], σ*-1)
@test isequal(jac[1,2], σ)
@test isequal(jac[1,3], 0)
@test isequal(jac[2,1], ρ-z)
@test isequal(jac[2,2], -1)
@test isequal(jac[2,3], x*-1)
@test isequal(jac[3,1], y)
@test isequal(jac[3,2], x)
@test isequal(jac[3,3], -1*β)

# Variable dependence checking in differentiation
@variables a(t) b(a)
@test !isequal(D(b), 0)

@test isequal(expand_derivatives(D(x * y)), simplify_constants(y*D(x) + x*D(y)))
@test_broken isequal(expand_derivatives(D(x * y)), simplify_constants(D(x)*y + x*D(y)))

@test isequal(expand_derivatives(D(2t)), 2)
@test isequal(expand_derivatives(D(2x)), 2D(x))


text:
using RDatasets
data = dataset("datasets", "iris")
first(data, 4)

typeof(data)

array_data = convert(Array, data)
typeof(array_data)

matrix_data = convert(Matrix, data[!, 1:4])
typeof(matrix_data)

data[!, 1:4] # Choosing columns
data[1:3, 1] # Choosing rows

names(data)

# Accessing colums several ways
data.Species
data[!, :Species]
data[!, [:Species, :SepalLength]]
unique(data[!, :Species])
levels(data[!, :Species]) # For categorical

max_stars_repo_name: coinslab/ComputationalCognitiveModeling

max_issues_repo_name: coinslab/ComputationalCognitiveModeling

max_forks_repo_name: coinslab/ComputationalCognitiveModeling


text:
module RunSettings
    #
    # This module contains things needed to control one run e.g. the output destination, number of households to use andd so on.
    #
    using Parameters

    using ScottishTaxBenefitModel
    using .Definitions
    using UUIDs

    export
        Settings,
        MT_Routing,
        uc_full,
        lmt_full,
        modelled_phase_in,

        IneqIncomeMeasure,
        TargetBCIncomes,

        bhc_net_income,
        eq_bhc_net_income,
        ahc_net_income,
        eq_ahc_net_income,
        DataIncomeSource,
        ds_hbai,
        ds_frs,

        ahc_hh,
        bhc_hh,
        total_bens,
        total_taxes

    @enum TargetBCIncomes ahc_hh bhc_hh total_bens total_taxes

    @enum MT_Routing uc_full lmt_full modelled_phase_in
    @enum IneqIncomeMeasure bhc_net_income eq_bhc_net_income ahc_net_income eq_ahc_net_income
    @enum DataIncomeSource ds_hbai ds_frs
    @with_kw mutable struct Settings
        uuid :: UUID = UUID("c2ae9c83-d24a-431c-b04f-74662d2ba07e")
        uid :: Int = 1 # placeholder for maybe a user somewhere
        run_name :: String = "default_run"
        start_year :: Int = 2015
        end_year :: Int = 2018
        scotland_full :: Bool = true
        weighted :: Bool = false
        data_dir :: String = MODEL_DATA_DIR
        household_name = "model_households_scotland"
        people_name    = "model_people_scotland"
        dump_frames :: Bool = false
        num_households :: Int = 0
        num_people :: Int = 0
        prices_file = "indexes.tab"
        to_y :: Int = 2022
        to_q :: Int = 2
        output_dir :: String = joinpath(tempdir(),"output")
        means_tested_routing :: MT_Routing = uc_full
        poverty_line :: Real = -1.0
        ineq_income_measure  :: IneqIncomeMeasure = eq_bhc_net_income
        growth :: Real = 0.02 # for time to exit poverty
        income_data_source :: DataIncomeSource = ds_hbai
        do_marginal_rates  :: Bool = false
        do_replacement_rates :: Bool = false
        replacement_rate_hours :: Int = 30
        # We jam on 68 here since we don't want changes to pension age
        # in the parameters to affect the numbers of people in
        # mr/rr calculations.
        mr_rr_upper_age :: Int = 68
        target_bc_income :: TargetBCIncomes = ahc_hh
        target_mr_rr_income :: TargetBCIncomes = ahc_hh
        mr_incr = 0.001
        requested_threads = 1
        impute_employer_pension = true
    end

end

max_stars_repo_name: grahamstark/ScottishTaxBenefitModel

max_issues_repo_name: grahamstark/ScottishTaxBenefitModel

max_forks_repo_name: grahamstark/ScottishTaxBenefitModel


text:
###############
# Prepare FMU #
###############

t_start = 0.0
t_stop = 1.0

myFMU = fmiLoad("SpringPendulum1D", ENV["EXPORTINGTOOL"], ENV["EXPORTINGVERSION"])
@test fmiIsCoSimulation(myFMU)
@test fmiIsModelExchange(myFMU)
comp = fmiInstantiate!(myFMU; loggingOn=false)
@test comp != 0
# choose FMU or FMUComponent
fmuStruct = nothing
envFMUSTRUCT = ENV["FMUSTRUCT"]
if envFMUSTRUCT == "FMU"
    fmuStruct = myFMU
elseif envFMUSTRUCT == "FMUCOMPONENT"
    fmuStruct = comp
end
sol = fmiSimulateCS(fmuStruct, t_start, t_stop)
@test sol.success
sol = fmiSimulateME(fmuStruct, t_start, t_stop)
@test sol.success
fmiUnload(myFMU)

myFMU = fmiLoad("SpringPendulum1D", ENV["EXPORTINGTOOL"], ENV["EXPORTINGVERSION"]; type=:ME)
@test myFMU.type == FMI.fmi2TypeModelExchange
comp = fmiInstantiate!(myFMU; loggingOn=false)
fmuStruct = nothing
envFMUSTRUCT = ENV["FMUSTRUCT"]
if envFMUSTRUCT == "FMU"
    fmuStruct = myFMU
elseif envFMUSTRUCT == "FMUCOMPONENT"
    fmuStruct = comp
end
sol = fmiSimulate(fmuStruct, t_start, t_stop)
@test sol.success
fmiUnload(myFMU)

myFMU = fmiLoad("SpringPendulum1D", ENV["EXPORTINGTOOL"], ENV["EXPORTINGVERSION"]; type=:CS)
@test myFMU.type == FMI.fmi2TypeCoSimulation
comp = fmiInstantiate!(myFMU; loggingOn=false)
fmuStruct = nothing
envFMUSTRUCT = ENV["FMUSTRUCT"]
if envFMUSTRUCT == "FMU"
    fmuStruct = myFMU
elseif envFMUSTRUCT == "FMUCOMPONENT"
    fmuStruct = comp
end
sol = fmiSimulate(fmuStruct, t_start, t_stop)
@test sol.success
fmiUnload(myFMU)