text:
using Revise

using DataFrames, Statistics, Distributions, StatsBase,
  LinearAlgebra, BenchmarkTools, Random, Optim, Zygote, UnPack

MFloat64 = Union{Missing, Float64}

Random.seed!(11)

# a simple AR(1) model
struct AR1{T<:Real, Tσ²<:Union{Nothing,Real}}
  ϕ::T
  b::T
  μ::T
  σ²::Tσ² #note- this is σ²e, not the overall σ² of the series
end

AR1(;ϕ, μ=missing, b=missing, σ²=nothing) = AR1(
  ϕ, coalesce(b, μ*(1-ϕ)), coalesce(μ, b/(1-ϕ)),σ²)

centerar1(ar1::AR1) = AR1(ar1.ϕ, 0.0, 0.0, ar1.σ²)

function estimatear1ols(x::Vector{T}) where T
  x̃ = @view x[2:end]
  Lx̃ = @view x[1:(end-1)]

  ϕ = cov(x̃, Lx̃)/(var(Lx̃))
  b = mean(x̃) - ϕ * mean(Lx̃)
  μ = b/(1-ϕ)

  ε = similar(x)
  ε[1] = x[1] - μ
  ε[2:end] = x̃ .- ϕ .* Lx̃
  σ² = var(ε)

  @assert b ≈ mean(x̃) - ϕ * mean(Lx̃)
  @assert ϕ ≈ cov(x̃ .- mean(x̃), Lx̃ .- mean(Lx̃))/var(
    Lx̃ .- mean(Lx̃))

  return AR1(;ϕ,b,σ²)
end

function estimatear1mle(x; tol=10^-10)

  #start by centering
  N=length(x)
  μ = mean(x)
  x̃ = @view(x[2:end]) .- μ
  Lx̃ = @view(x[1:(end-1)]) .- μ

  x₁ = x[1]
  #computes the next iteration of ϕ given a previous ϕ
  cσ²(ϕ) = (x₁^2*(1-ϕ^2) + sum((x̃ .- ϕ*Lx̃).^2))/N
  function objϕ(ϕ)
    σ² = cσ²(ϕ)
    return (sum((x̃ .- ϕ * Lx̃) .* Lx̃) - (σ²/(1-ϕ^2) - x₁^2)*ϕ)^2
  end

  ϕ = optimize(objϕ, -1.01, 1.01, Brent()) |> Optim.minimizer

  return AR1(;ϕ, μ, σ² = cσ²(ϕ))
end

function estimatear1mlecheck(x; tol=10^-10)

  #start by centering
  N=length(x)
  μ = mean(x)
  x̃ = @view(x[2:end]) .- μ
  Lx̃ = @view(x[1:(end-1)]) .- μ

  x₁ = x[1]

  l(ϕ, σ²) = (-N/2*log(2*π*σ²)
    + 1/2*log(1-ϕ^2)
    - x₁^2*(1-ϕ^2)/(2*σ²)
    - sum((x̃ .- ϕ .* Lx̃).^2)/(2*σ²)) * -1

  l(Φ) = l(Φ[1], Φ[2])

  ϕ₀ = cov(x̃, Lx̃)/var(Lx̃) #start by guessing the least squares solution
  ε₀= x̃ - ϕ₀*Lx̃
  σ²₀ = var(ε₀)

  ϕ, σ² = optimize(l,  [-1.0,0.0], [1.0,Inf], [ϕ₀, σ²₀], Fminbox(LBFGS()),
    autodiff = :forward) |> Optim.minimizer

  return AR1(;ϕ, μ, σ²)
end

AR1(x::AbstractVector; estimationmethod=estimatear1mle) = estimationmethod(x)

function checkmath2(x)
  if length(x) > 250
    x = x[1:250] |> deepcopy
    @warn "length(x) > 250 so truncating length to 250. This is only for the purposes
    of computing the full (not log) likelihood"
  end

  N=length(x)
  μ = mean(x)
  x̃ = @view(x[2:end]) .- μ
  Lx̃ = @view(x[1:(end-1)]) .- μ

  ϕ₀ = 0.75
  σ²₀ = 2.2
  x₁ = x̃[1]

  L(ϕ, σ²) =  (1-ϕ^2)^(1/2)/σ²^(1/2)*pdf(Normal(), x₁/(σ²/(1-ϕ^2))^(1/2))*prod(
    ((x̃ᵢ,Lx̃ᵢ)->1/σ²^(1/2)*pdf(Normal(), (x̃ᵢ-ϕ*Lx̃ᵢ)/σ²^(1/2))).(x̃,Lx̃))

  l(ϕ, σ²) = (-N/2*log(2*π*σ²)
    + 1/2*log(1-ϕ^2)
    - x₁^2*(1-ϕ^2)/(2*σ²)
    - sum((x̃ .- ϕ .* Lx̃).^2)/(2*σ²))

  @assert log(L(ϕ₀, σ²₀)) ≈ l(ϕ₀, σ²₀,) "log(L(ϕ, σ²)): $(log(L(ϕ₀, σ²₀)))
    l(ϕ, σ²): $(l(ϕ₀, σ²₀))"

  #now verify the derivitive and conditional wrt sigma^2
  lσ²(ϕ=ϕ₀, σ²=σ²₀,) = (-N/(2σ²)
    + x₁^2*(1-ϕ^2)/(2σ²^2)
    + sum((x̃ .- ϕ .* Lx̃).^2)/(2σ²^2))
  lσ²z(ϕ=ϕ₀, σ²=σ²₀,) = gradient(σ²->l(ϕ, σ²,), σ²)[1]
  lσ²n() =  (l(ϕ₀, σ²₀+1e-6,)-l(ϕ₀, σ²₀,))/1e-6
  @assert lσ²() ≈ lσ²z() "lσ²(...): $(lσ²())
    lσ²z(...): $(lσ²z())"

  cσ²(ϕ) = (x₁^2*(1-ϕ^2) + sum((x̃ .- ϕ*Lx̃).^2))/N
  @assert 1 + lσ²z(ϕ₀, cσ²(ϕ₀)) ≈ 1.0 "cσ²(ϕ₀): $(cσ²(ϕ₀))
    lσ²z(ϕ₀, cσ²(ϕ₀)): $(lσ²z(ϕ₀, cσ²(ϕ₀)))"

  #finally, verify the derivitive wrt ϕ
  lϕ(ϕ=ϕ₀, σ²=σ²₀) = (-ϕ/(1-ϕ^2) +
    x₁^2*ϕ/σ² +
    sum((x̃ .- ϕ * Lx̃) .* Lx̃)/σ²)
  lϕz(ϕ=ϕ₀, σ²=σ²₀,) = gradient(ϕ->l(ϕ, σ²), ϕ)[1]
  @assert lϕ() ≈ lϕz() "lϕ(...): $(lϕ())
    lϕz(...): $(lϕz())"

  ϕcerr(ϕ, σ²=σ²₀)= sum((x̃ .- ϕ * Lx̃) .* Lx̃) - (σ²/(1-ϕ^2) - x₁^2)*ϕ
  cϕ(σ²=σ²₀,) = optimize(ϕ->ϕcerr(ϕ)^2, -1.01, 1.01, Brent(), abs_tol=10^-14) |> Optim.minimizer
  @assert isapprox(lϕz(cϕ()), 0.0,atol=1e-5) "cϕ(ϕ₀): 2$(cϕ())
    lϕz(ϕcerr()): $(lϕz(ϕcerr()))"
  #lμ()
end

#simulates an AR1 series
function simulatear1!(ar1::AR1, x::AbstractVector)
  @unpack ϕ, b, σ², μ = ar1

  σ² === nothing && throw(" σ² must not be nothing")

  #first simulate the centered distribution
  #we will do this all in place
  σ=σ²^0.5
  rand!(Normal(0,σ), x)
  x[1] *= 1/(1-ϕ^2)^0.5 # the first entry is simulated from the unconditional distribution
  for t ∈ 2:length(x)
    x[t] += ϕ * x[t-1]
  end
  x .+= μ

  return x
end

#which approach is faster?
function simulatear1(ar1::AR1, T::Int)
  @unpack ϕ, σ², μ = ar1

  σ² === nothing && throw(" σ² must not be nothing")
  #T::Int = length(x)

  #first simulate the centered distribution
  #we will do this all in place
  σ=σ²^0.5
  x = rand(Normal(), T) .* σ
  x[1] *= 1/(1-ϕ^2)^0.5
  for t ∈ 2:T
    x[t] += ϕ * x[t-1]
  end
  x .+= μ

  return x
end

simulatear1!(ar1::AR1{T,Tσ²}, N::Int) where {T,Tσ²} = simulatear1!(ar1, Vector{T}(undef, N))

compareoncear1xy!(xar1::AR1, x::AbstractVector, yar1::AR1,  y::AbstractVector,
  comparefunc::Tcomparefunc
  ) where Tcomparefunc = comparefunc(simulatear1!(xar1,x), simulatear1!(yar1,y))

#this runs the simulations- basically a parametric bootstrap under the null
#statfunc is computed for each pass, while aggfunc is computed on the statistics
#at the end of the pass
function comparear1undernull(x::Tx,y::AbstractVector,
    statfunc::Tstatfunc, aggfunc::Taggfunc; Nsimulations::Int,
    checkstationarity::Bool = true, stationaritytol::Float64=10^-4
    ) where {Tx, Tstatfunc, Taggfunc}
  xar1 = AR1(x) |> centerar1
  yar1 = AR1(y) |> centerar1

  results = Dict(s => Vector{MFloat64}(undef, Nsimulations) for s ∈ keys(statfunc(x,y)))

  #this assumes stationary results are meaningless
  if checkstationarity && max(abs(xar1.ϕ), abs(yar1.ϕ)) ≥ 1.0-stationaritytol
    xar1 = AR1(NaN, NaN, xar1.μ, xar1.σ²)
    yar1 = AR1(NaN, NaN, xar1.μ, xar1.σ²)
  end

  xs = [similar(x) for t ∈ 1:Threads.nthreads()]
  ys = [similar(y) for t ∈ 1:Threads.nthreads()]

  #draw a sample and compute the aggregate function, sotring the reuslts in an array
  Threads.@threads for n ∈ 1:Nsimulations
    t::Int = Threads.threadid()
    res = compareoncear1xy!(xar1, xs[t], yar1, ys[t], statfunc)

    #this structure allows us to compute multiple statistics
    for s ∈ keys(res)
      @assert results[s][n] ≡ missing
      results[s][n] = res[s]
    end
  end

  return (; Dict(s=>aggfunc(skipmissing(results[s]), s) for s ∈ keys(results))...)
end

function comparearima11undernull(x::Tx,y::AbstractVector,
    statfunc::Tstatfunc, aggfunc::Taggfunc; Nsimulations::Int,
    checkstationarity::Bool = true, stationaritytol::Float64=10^-4) where {Tx, Tstatfunc, Taggfunc}

    #difference everything
    dx = x[2:end] .- x[1:(end-1)]
    dy = x[2:end] .- y[1:(end-1)]

    dxar1 = AR1(dx) |> centerar1
    dyar1 = AR1(dy) |> centerar1

    #make integration funcitons
    #ix(δx) = cumsum([x[1]; δx])
    #iy(δy) = cumsum([y[1]; δy])

    #dstatfunc(δx, δy) = statfunc(δx |> ix, δy |> iy)

    results = Dict(s => Vector{MFloat64}(undef, Nsimulations) for s ∈ keys(statfunc(dx,dy)))

    #this assumes stationary results are meaningless
    if checkstationarity && max(abs(dxar1.ϕ), abs(dyar1.ϕ)) ≥ 1.0-stationaritytol
      dxar1 = AR1(NaN, NaN, dxar1.μ, dxar1.σ²)
      dyar1 = AR1(NaN, NaN, dxar1.μ, dxar1.σ²)
    end

    dxs = [similar(dx) for t ∈ 1:Threads.nthreads()]
    dys = [similar(dy) for t ∈ 1:Threads.nthreads()]

    #draw a sample and compute the aggregate function, sotring the reuslts in an array
    Threads.@threads for n ∈ 1:Nsimulations
      t::Int = Threads.threadid()
      res = compareoncear1xy!(dxar1, dxs[t], dyar1, dys[t], statfunc)

      #this structure allows us to compute multiple statistics
      for s ∈ keys(res)
        @assert results[s][n] ≡ missing
        results[s][n] = res[s]
      end
    end

    return (; Dict(s=>aggfunc(skipmissing(results[s]), s) for s ∈ keys(results))...)
end

nonparametricnullp(v::AbstractVector{T}, x) where T<:Real = (
  x ≥ 0 ? sum(v .> x)/length(v) : sum(v .< x)/length(v))
nonparametricnullp(v, x) = nonparametricnullp(v |> skipmissing |> collect, x)

function testar1(;N)
  ϕ = 0.8
  ε = rand(Normal(), N) + rand(Normal(1,1), N)
  x = missings(Float64, N)
  x[1] = ε[1]

  global x = x
  for t ∈ 2:N
    x[t] = ϕ*x[t-1] + ε[t]
  end
  @info "ϕ=$ϕ, b=$(mean(ε)), μ=$(mean(x))"

  checkmath2(x)

  #=estolsar1 = @btime AR1($x, estimationmethod=estimatear1ols)
  @info "ols conditional: $estolsar1"

  estmlear1 = @btime AR1($x, estimationmethod=estimatear1mle)
  @info "mle: $estmlear1"

  estmle2ar1 = @btime AR1($x, estimationmethod=estimatear1mlecheck)
  @info "mle2: $estmle2ar1"=#

  actar1=AR1(;ϕ, b=0.5, σ²=2)
  xsim = simulatear1(actar1, 10^6)
  @assert isapprox(mean(xsim), 2.5, atol=0.01)  "mean(xsim): $(mean(xsim))"
  @assert isapprox(var(xsim), 2/(1-ϕ^2), atol=0.01) "var(xsim): $(var(xsim))"
  ar1sim = AR1(xsim)
  @info "ar1sim: $(ar1sim)"
  @assert isapprox(ar1sim.ϕ, ϕ, atol=0.01)
  @assert isapprox(ar1sim.b, 0.5, atol=0.01)
  @assert isapprox(ar1sim.μ, 2.5, atol=0.01)
  @assert isapprox(ar1sim.σ², 2, atol=0.01)

  xsim! = simulatear1!(actar1, Vector{Float64}(undef,10^6))
  @assert isapprox(mean(xsim), 2.5, atol=0.01)  "mean(xsim): $(mean(xsim))"
  @assert isapprox(var(xsim), 2/(1-ϕ^2), atol=0.01) "var(xsim): $(var(xsim))"
  ar1sim = AR1(xsim!)
  @info "ar1sim: $(ar1sim)"
  @assert isapprox(ar1sim.ϕ, ϕ, atol=0.01)
  @assert isapprox(ar1sim.b, 0.5, atol=0.01)
  @assert isapprox(ar1sim.μ, 2.5, atol=0.01)
  @assert isapprox(ar1sim.σ², 2, atol=0.01)

  Nsimulations=10^4
  dmeantest = 0.2
  #recall the compare func will pass the name of the relevant stat func when it is called
  aggfunc(v, s) = nonparametricnullp(v, Dict(:cor=>0.1, :dmean=>dmeantest)[s])
  statfunc(x,y) = Dict(:cor=>cor(x,y), :dmean=>mean(x)-mean(y))

  @info "std(x): $(std(x))"
  @info "expected std(x): $((2/(1-ϕ^2))^0.5)"
  Sdmean = (2*2/(1-ϕ^2)/length(x))^0.5
  Epdmean = cdf(Normal(), - dmeantest/Sdmean)
  res = comparear1undernull(x, deepcopy(x), statfunc, aggfunc; Nsimulations)
  @info "compareundernull result: $res (Epdmean: $Epdmean)"
  y = deepcopy(x)

  @btime comparear1undernull($x, $y, $statfunc, $aggfunc; Nsimulations=$Nsimulations)
  throw("simulation complete")

  ix = cumsum(x .+ 0.2)
  iy = cumsum(x .+ 0.2)
  res = comparear1undernull(ix, deepcopy(ix), statfunc, aggfunc; Nsimulations)
  @info "compareundernull dmean: $res"
  res = comparearima11undernull(ix, deepcopy(ix), statfunc, aggfunc; Nsimulations)
  @info "compareundernull dmean: $res"

end

@time testar1(N=1000)

max_stars_repo_name: clintonTE/CCA

max_issues_repo_name: clintonTE/CCA

max_forks_repo_name: clintonTE/CCA


text:
Base.inv(x::T) where T <: CasadiSymbolicObject = pycall(casadi.inv, T, x)
Base.sqrt(x::T) where T <: CasadiSymbolicObject = pycall(casadi.sqrt, T, x)
Base.sin(x::T) where T <: CasadiSymbolicObject = pycall(casadi.sin, T, x)
Base.cos(x::T) where T <: CasadiSymbolicObject = pycall(casadi.cos, T, x)
Base.vec(x::T) where T <: CasadiSymbolicObject = pycall(casadi.vec, T, x)
Base.transpose(x::T) where T <: CasadiSymbolicObject = pycall(casadi.transpose, T, x)

Base.size(x::CasadiSymbolicObject) = x.size()
Base.reshape(x::T, t::Tuple{Int, Int}) where T <: CasadiSymbolicObject = pycall(casadi.reshape, T, x, t)

ext: jl


text:
function _contains(convex_sets::COSMO.CompositeConvexSet, t::Type{<:COSMO.AbstractConvexCone})
  for set in convex_sets.sets
    if typeof(set) <: t
      return true
    end
  end
  return false
end

function chordal_decomposition!(ws::COSMO.Workspace)
  # do nothing if no psd cones present in the problem
  if !_contains(ws.p.C, DecomposableCones{Float64})
    ws.settings.decompose = false
    return nothing
  end
  ws.ci = ChordalInfo{Float64}(ws.p)

  find_sparsity_patterns!(ws)

  if ws.ci.num_decomposable > 0
    # find transformation matrix H and new composite convex set
    find_decomposition_matrix!(ws)

    # augment the original system
    augment_system!(ws)
    pre_allocate_variables!(ws)
  else
    ws.settings.decompose = false
  end
  nothing
end

# analyse PSD cone constraints for chordal sparsity pattern
function find_sparsity_patterns!(ws)
  row_ranges = get_set_indices(ws.p.C.sets)
  sp_ind = 1
  for (k, C) in enumerate(ws.p.C.sets)
    psd_row_range = row_ranges[k]
    csp = find_aggregate_sparsity(ws.p.A, ws.p.b, psd_row_range, C)
    sp_ind = analyse_sparsity_pattern!(ws.ci, csp, ws.p.C.sets, C, k, sp_ind)
  end
end

function analyse_sparsity_pattern!(ci, csp, sets, C::PsdCone{T}, k, sp_ind) where {T <: Real}
  if length(csp) < C.dim
    ordering = find_graph!(ci, csp, C.sqrt_dim, C)
    ci.sp_arr[sp_ind] = COSMO.SparsityPattern(ci.L, C.sqrt_dim, ordering)
    push!(ci.psd_cones_ind, k)
    ci.num_decomposable += 1
    return sp_ind + 1
  else
   sets[k] = COSMO.DensePsdCone{T}(C.dim)
   return sp_ind
  end
end

function analyse_sparsity_pattern!(ci, csp, sets, C::PsdConeTriangle{T}, k, sp_ind) where {T <: Real}
  if length(csp) < C.dim
    ordering = find_graph!(ci, csp, C.sqrt_dim, C)
    ci.sp_arr[sp_ind] = COSMO.SparsityPattern(ci.L, C.sqrt_dim, ordering)
    push!(ci.psd_cones_ind, k)
    ci.num_decomposable += 1
    return sp_ind + 1
  else
   sets[k] = COSMO.DensePsdConeTriangle{T}(C.dim)
    return sp_ind
  end
end

analyse_sparsity_pattern!(ci, csp, sets, C::AbstractConvexSet, k, sp_ind) = sp_ind

function nz_rows(a::SparseMatrixCSC, ind::UnitRange{Int64}, DROP_ZEROS_FLAG::Bool)
  DROP_ZEROS_FLAG && dropzeros!(a)
  active = falses(length(ind))
  for r in a.rowval
    if in(r, ind)
      active[r - ind.start + 1] = true
    end
  end
  return findall(active)
end

function number_of_overlaps_in_rows(A::SparseMatrixCSC)
  # sum the entries row-wise
  numOverlaps = sum(A, dims = 2)
  ri = findall(x -> x > 1, numOverlaps)
  return ri, numOverlaps[ri]
end

function find_aggregate_sparsity(A, b, ind, C::DecomposableCones{ <: Real})
  AInd = nz_rows(A, ind, false)
  # commonZeros = AInd[find(x->x==0,b[AInd])]
  bInd = findall(x -> x != 0, view(b, ind))
  commonNZeros = union(AInd, bInd)
  return commonNZeros
end
find_aggregate_sparsity(A, b, ind, C::AbstractConvexSet) = Int64[]

function vec_to_mat_ind(ind::Int64, n::Int64)
  ind > n^2 && error("Index ind out of range.")
  ind == 1 && (return 1, 1)

  r = ind % n

  if r == 0
    j = Int(ind / n)
    i = n
  else
    j = Int(floor(ind / n) + 1)
    i = r
  end
  return i, j
end

# Converts the matrix element (i, j) of A ∈ m x n into the corresponding linear index of v = vec(A)
function mat_to_vec_ind(i::Int64, j::Int64, m::Int64)
  (i > m || i <= 0 || j <= 0) && throw(BoundsError("Indices outside matrix bounds."))
  return (j - 1) * m + i
end

# Converts the matrix element (i, j) of A ∈ m x n into the corresponding linear index of v = svec(A, ::UpperTriangle)
function mat_to_svec_ind(i::Int64, j::Int64)
  if i <= j
    return div((j - 1) * j, 2) + i
  else
    return div((i - 1) * i, 2) + j
  end
end

# function svec_to_mat_ind(k::Int64)
#   j = isqrt(2 * k)
#   i = k - div((j - 1) * j, 2)
#   return i, j
# end

# function finds the transformation matrix H to decompose the vector s into its parts and stacks them into sbar
function find_decomposition_matrix!(ws)

  # allocate H and new decomposed cones
  n = COSMO.find_H_col_dimension(ws.p.C.sets, ws.ci.sp_arr)
  H_I = zeros(Int64, n)

  # find number of decomposed and total sets and allocate structure for new compositve convex set
  num_total, num_new_psd_cones = COSMO.num_cone_decomposition(ws)
  # decomposed_psd_cones = Array{COSMO.PsdCone}(undef, 0)
  C_new = Array{COSMO.AbstractConvexSet{Float64}}(undef, num_total)
  C_new[1] = COSMO.ZeroSet{Float64}(ws.ci.originalM)

  # loop over all convex sets and fill H and composite convex set accordingly
  row = 1
  entry = 1
  sp_ind = 1
  set_ind = 2
  for (kkk, C) in enumerate(ws.p.C.sets)
    set_ind, sp_ind, entry = COSMO.decompose!(H_I, C_new, set_ind, C, entry, row, ws.ci.sp_arr, sp_ind)
    row += C.dim
  end
  ws.p.C = COSMO.CompositeConvexSet(C_new)
  ws.ci.H = sparse(H_I, collect(1:n), ones(n))
end

function decompose!(H_I::Vector{Int64}, C_new, set_ind::Int64, C::COSMO.AbstractConvexSet, entry::Int64, row::Int64, sp_arr::Array{SparsityPattern}, sp_ind::Int64)
  #H[row:row + C.dim - 1, col:col + C.dim - 1] = sparse(1.0I, C.dim, C.dim)
  for i = 1:C.dim
    H_I[entry] = row + i - 1
    entry += 1
  end
  C_new[set_ind] = C

  return set_ind + 1, sp_ind, entry
end

get_blk_length(nBlk::Int64, C::COSMO.PsdCone{Float64}) = isqrt(nBlk)
function get_blk_length(nBlk::Int64, C::COSMO.PsdConeTriangle{Float64})
    d = isqrt(nBlk)
    return div((d + 1) * d, 2)
end
get_blk_rows(d::Int64, C::COSMO.PsdCone{Float64}) = d^2
get_blk_rows(d::Int64, C::COSMO.PsdConeTriangle{Float64}) = d

function decompose!(H_I::Vector{Int64}, C_new, set_ind::Int64,  C::DecomposableCones{ <: Real}, entry::Int64, row_start::Int64, sp_arr::Array{SparsityPattern}, sp_ind::Int64)
  sparsity_pattern = sp_arr[sp_ind]
  sntree = sparsity_pattern.sntree
  # original matrix size
  original_size = C.sqrt_dim

  for iii = 1:num_cliques(sntree)
    # new stacked size
    block_length = get_blk_length(sntree.nBlk[iii], C)

    c = COSMO.get_clique(sntree, iii)
    # the graph and tree algorithms determined the cliques vertices of an AMD-permuted matrix. Since the location of the data hasn't changed in reality, we have to map the clique vertices back
    c = map(v -> sparsity_pattern.ordering[v], c)
    sort!(c)
    entry = COSMO.add_subblock_map!(H_I, c, block_length, original_size, entry, row_start, C)

    # create and add new cone for subblock
    num_rows = get_blk_rows(block_length, C)
    C_new[set_ind] = typeof(C)(num_rows)
    set_ind += 1
  end
  return set_ind, sp_ind + 1, entry
end

# fills the corresponding entries of H for clique c
function add_subblock_map!(H_I::Vector{Int64}, clique_vertices::Array{Int64}, block_size::Int64, original_size::Int64, entry::Int64, row_start::Int64, ::PsdCone{<: Real})
  for vj in clique_vertices
    for vi in clique_vertices
      row = mat_to_vec_ind(vi, vj, original_size)
      H_I[entry] = row_start + row - 1
      entry += 1
    end
  end
  return entry::Int64
end

function add_subblock_map!(H_I::Vector{Int64}, clique_vertices::Array{Int64}, block_dim::Int64, original_size::Int64, entry::Int64,  row_start::Int64, ::PsdConeTriangle{<: Real})
  for vj in clique_vertices
    for vi in clique_vertices
      if vi <= vj
        row = mat_to_svec_ind(vi, vj)
        H_I[entry] = row_start + row - 1
        entry += 1
      end
    end
  end
  return entry::Int64
end

function find_H_col_dimension(sets, sp_arr)
  sum_cols = 0
  sp_arr_ind = 1
  for C in sets
    dim, sp_arr_ind = decomposed_dim(C, sp_arr, sp_arr_ind)
    sum_cols += dim
  end
  return sum_cols::Int64
end

decomposed_dim(C::AbstractConvexSet, sp_arr::Array{SparsityPattern}, sp_arr_ind::Int64) = (C.dim, sp_arr_ind)
function decomposed_dim(C::PsdCone{<: Real}, sp_arr::Array{SparsityPattern}, sp_arr_ind::Int64)
  dim = sum(sp_arr[sp_arr_ind].sntree.nBlk)
  return dim::Int64, (sp_arr_ind + 1)::Int64
end

function decomposed_dim(C::PsdConeTriangle{<: Real}, sp_arr::Array{SparsityPattern}, sp_arr_ind::Int64)
  dim = sum(map(d -> div((d + 1) * d, 2), map(x -> isqrt(x), sp_arr[sp_arr_ind].sntree.nBlk)))
  return dim::Int64, (sp_arr_ind + 1)::Int64
end

function num_cone_decomposition(ws)
  num_sets = length(ws.p.C.sets)
  num_old_psd_cones = length(ws.ci.psd_cones_ind)
  num_new_psd_cones = 0
  for iii = 1:num_old_psd_cones
    sp = ws.ci.sp_arr[iii]
    num_new_psd_cones += COSMO.num_cliques(sp.sntree)
  end
  ws.ci.num_decom_psd_cones = ws.ci.num_psd_cones - ws.ci.num_decomposable + num_new_psd_cones
  # the total number is the number of original non-psd cones + number of decomposed psd cones + 1 zeroset for the augmented system
  num_total = num_sets - num_old_psd_cones + num_new_psd_cones + 1
  return num_total, num_new_psd_cones
end

function augment_system!(ws)
  _augment!(ws.p, ws.ci.H)
  nothing
end

function _augment!(problem, H::SparseMatrixCSC)
  mH, nH = size(H)
  m, n = size(problem.A)
  problem.P = blockdiag(problem.P, spzeros(nH, nH))
  problem.q = vec([problem.q; zeros(nH)])
  problem.A = [problem.A H; spzeros(nH, n) -sparse(1.0I, nH, nH)]
  problem.b = vec([problem.b; zeros(nH)])
  problem.model_size[1] = size(problem.A, 1)
  problem.model_size[2] = size(problem.A, 2)
  nothing
end

function reverse_decomposition!(ws::COSMO.Workspace, settings::COSMO.Settings)
  mO = ws.ci.originalM
  nO = ws.ci.originalN

  H = ws.ci.H
  vars = Variables{Float64}(mO, nO, ws.ci.originalC)
  vars.x .= ws.vars.x[1:nO]
  vars.s  .= SplitVector{Float64}(H * ws.vars.s[mO + 1:end], ws.ci.originalC)

  # fill dual variables such that μ_k  = H_k μ for k=1,...,p
  fill_dual_variables!(ws, vars)

  # if user requests, perform positive semidefinite completion on entries of μ that were not in the decomposed blocks
  ws.vars = vars
  ws.p.C = ws.ci.originalC
  settings.complete_dual && psd_completion!(ws)

  return nothing
end

function fill_dual_variables!(ws::COSMO.Workspace, vars::COSMO.Variables)
  mO = ws.ci.originalM
  H = ws.ci.H

  # this performs the operation μ = sum H_k^T *  μ_k causing an addition of (identical valued) overlapping blocks
  vars.μ .= H * ws.vars.μ[mO + 1:end]

  # # to remove the overlaps we take the average of the values for each overlap by dividing by the number of blocks that overlap in a particular entry, i.e. number of 1s in each row of H
  rowInd, nnzs = number_of_overlaps_in_rows(H)

  for iii=1:length(rowInd)
    ri = rowInd[iii]
    vars.μ[ri] .= vars.μ[ri] / nnzs[iii]
  end
end

# The psd entries of μ that correspond to the zeros in s are not constrained by the problem
# however, in order to make the dual psd cone positive semidefinite we have to do a
# positive semidefinite completion routine to choose the values
function psd_completion!(ws::COSMO.Workspace)

  # loop over psd cones
  row_ranges = get_set_indices(ws.p.C.sets)
  sp_ind = 1
  for (kkk, C) in enumerate(ws.p.C.sets)
    sp_ind = complete!(ws.vars.μ, C, ws.ci.sp_arr, sp_ind, row_ranges[kkk])
  end

  return nothing
end

complete!(μ::AbstractVector, ::AbstractConvexSet, sp_arr::Array{SparsityPattern}, sp_ind::Int64, rows::UnitRange{Int64}) = sp_ind

function complete!(μ::AbstractVector, C::PsdCone{<: Real}, sp_arr::Array{SparsityPattern}, sp_ind::Int64, rows::UnitRange{Int64})
  sp = sp_arr[sp_ind]
  M = reshape(view(μ, rows), C.sqrt_dim, C.sqrt_dim)
  psd_complete!(M, C.sqrt_dim, sp.sntree, sp.ordering)
  return sp_ind + 1
end

function complete!(μ::AbstractVector, C::PsdConeTriangle{<: Real}, sp_arr::Array{SparsityPattern}, sp_ind::Int64, rows::UnitRange{Int64})
  sp = sp_arr[sp_ind]

  μ_view = view(μ, rows)
  populate_upper_triangle!(C.X, μ_view, 1. / sqrt(2))
  psd_complete!(C.X, C.sqrt_dim, sp.sntree, sp.ordering)
  extract_upper_triangle!(C.X, μ_view, sqrt(2))
  return sp_ind + 1
end

# positive semidefinite completion (from Vandenberghe - Chordal Graphs..., p. 362)
# input: A - positive definite completable matrix
function psd_complete!(A::AbstractMatrix, N::Int64, sntree::SuperNodeTree, p::Array{Int64})

  ip = zeros(Int64, length(p))
  ip[p] = 1:length(p)

  # permutate matrix based on ordering p (p must be a vector type), W is in the order that the cliques are based on
  W = A[p, p]
  W = Matrix(W)
  num_cliques = COSMO.num_cliques(sntree)
  # go through supernode tree in descending order (given a post-ordering)
  snd_post_order = sntree.snd_post
  for j = (num_cliques - 1):-1:1
    snd_id = snd_post_order[j]

    # in order to obtain ν, α the vertex numbers of the supernode are mapped to the new position of the permuted matrix
    # index set of snd(i) sorted using the numerical ordering i,i+1,...i+ni
    ν = get_snd(sntree, snd_id)
    clique = get_clique(sntree, snd_id)
    # index set containing the elements of col(i) \ snd(i) sorted using numerical ordering σ(i)
    α = filter(x -> !in(x, ν), clique)

    # index set containing the row indizes of the lower-triangular zeros in column i (i: representative index) sorted by σ(i)
    i = ν[1]
    η = collect(i + 1:1:N)
    # filter out elements in lower triangular part of column i that are non-zero
    filter!(x -> !in(x, α), η)

    Waa = view(W, α, α)
    Waa = view(W, α, α)
    Wαν = view(W, α, ν)
    Wηα = view(W, η, α)

    Y = Waa \ Wαν

    W[η, ν] =  Wηα * Y
    # symmetry condition
    W[ν, η] = view(W, η, ν)'
  end

  # invert the permutation
  A[:, :] = W[ip, ip]
end


text:
module RulesDotProductIn2Test

using Test
using ReactiveMP
using Random
import ReactiveMP: @test_rules
import LinearAlgebra: dot

@testset "rules:typeof(dot):in2" begin

     @testset "Belief Propagation: (m_out::UnivariateNormalDistributionsFamily, m_in1::PointMass)" begin

        @test_rules [ with_float_conversions = true ] typeof(dot)(:in2, Marginalisation) [
            (input = (m_out = NormalMeanVariance(2.0, 2.0), m_in1 = PointMass(-1.0), meta = NoCorrection()), output = NormalWeightedMeanPrecision(-1.0, 0.5)),
            (input = (m_out = NormalMeanPrecision(1.0, 1.0), m_in1 = PointMass(-2.0), meta = NoCorrection()), output = NormalWeightedMeanPrecision(-2.0, 4.0)),
            (input = (m_out = NormalWeightedMeanPrecision(2.0, 1.0), m_in1 = PointMass(-1.0), meta = NoCorrection()), output = NormalWeightedMeanPrecision(-2.0, 1.0)),
        ]

        @test_rules [ with_float_conversions = true ] typeof(dot)(:in2, Marginalisation) [
            (input = (m_out = NormalMeanVariance(2.0, 2.0), m_in1 = PointMass(-1.0), meta = TinyCorrection()), output = NormalWeightedMeanPrecision(-1.0, 0.5)),
            (input = (m_out = NormalMeanPrecision(1.0, 1.0), m_in1 = PointMass(-2.0), meta = TinyCorrection()), output = NormalWeightedMeanPrecision(-2.0, 4.0)),
            (input = (m_out = NormalWeightedMeanPrecision(2.0, 1.0), m_in1 = PointMass(-1.0), meta = TinyCorrection()), output = NormalWeightedMeanPrecision(-2.0, 1.0)),
        ]

        @test_rules [ with_float_conversions = true ] typeof(dot)(:in2, Marginalisation) [
            (input = (m_out = NormalMeanVariance(2.0, 2.0), m_in1 = PointMass(0.0), meta = NoCorrection()), output = NormalWeightedMeanPrecision(0.0, 0.0)),
            (input = (m_out = NormalMeanPrecision(1.0, 1.0), m_in1 = PointMass(0.0), meta = NoCorrection()), output = NormalWeightedMeanPrecision(0.0, 0.0)),
            (input = (m_out = NormalWeightedMeanPrecision(2.0, 1.0), m_in1 = PointMass(0.0), meta = NoCorrection()), output = NormalWeightedMeanPrecision(0.0, 0.0)),
        ]

        @test_rules [ with_float_conversions = true ] typeof(dot)(:in2, Marginalisation) [
            (input = (m_out = NormalMeanVariance(2.0, 2.0), m_in1 = PointMass(0.0), meta = TinyCorrection()), output = NormalWeightedMeanPrecision(0.0, tiny)),
            (input = (m_out = NormalMeanPrecision(1.0, 1.0), m_in1 = PointMass(0.0), meta = TinyCorrection()), output = NormalWeightedMeanPrecision(0.0, tiny)),
            (input = (m_out = NormalWeightedMeanPrecision(2.0, 1.0), m_in1 = PointMass(0.0), meta = TinyCorrection()), output = NormalWeightedMeanPrecision(0.0, tiny)),
        ]

    end

    @testset "Belief Propagation: (m_out::UnivariateNormalDistributionsFamily, m_in1::PointMass)" begin

        @test_rules [ with_float_conversions = true ] typeof(dot)(:in2, Marginalisation) [
            (input = (m_out = NormalMeanVariance(2.0, 2.0), m_in1 = PointMass([-1.0, 1.0]), meta = NoCorrection()), output = MvNormalWeightedMeanPrecision([-1.0, 1.0], [0.5 -0.5; -0.5 0.5])),
            (input = (m_out = NormalMeanPrecision(1.0, inv(2.0)), m_in1 = PointMass([2.0, 1.0]), meta = NoCorrection()), output = MvNormalWeightedMeanPrecision([1.0, 0.5], [2.0 1.0; 1.0 0.5])),
            (input = (m_out = NormalWeightedMeanPrecision(1.0, 1.0), m_in1 = PointMass([-1.0, 3.0]), meta = NoCorrection()), output = MvNormalWeightedMeanPrecision([-1.0, 3.0], [1.0 -3.0; -3.0 9.0])),
        ]

        @test_rules [ with_float_conversions = false ] typeof(dot)(:in2, Marginalisation) [
            (input = (m_out = NormalMeanVariance(2.0, 2.0), m_in1 = PointMass([-1.0, 1.0]), meta = TinyCorrection()), output = MvNormalWeightedMeanPrecision([-1.0, 1.0], [0.5 + tiny -0.5; -0.5 0.5 + tiny])),
            (input = (m_out = NormalMeanPrecision(1.0, inv(2.0)), m_in1 = PointMass([2.0, 1.0]), meta = TinyCorrection()), output = MvNormalWeightedMeanPrecision([1.0, 0.5], [2.0 + tiny 1.0; 1.0 0.5 + tiny])),
            (input = (m_out = NormalWeightedMeanPrecision(1.0, 1.0), m_in1 = PointMass([-1.0, 3.0]), meta = TinyCorrection()), output = MvNormalWeightedMeanPrecision([-1.0, 3.0], [1.0 + tiny -3.0; -3.0 9.0 + tiny])),
        ]

    end

end

end


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

using Markdown
using InteractiveUtils

# ╔═╡ cc2cf6c0-8e92-11eb-1255-53074ed0e59c
using ShortCodes

# ╔═╡ 3fcaf0b4-8f47-11eb-2145-311f5d76757f
begin
	using .MyOtherModule: goo
	goo
end

# ╔═╡ f1aa59ae-8e94-11eb-12bc-3b1d9ccc314b
html"""

"""

# ╔═╡ 7400844c-8e8f-11eb-27d0-859ea6dba632
md"""
# The Julia Programming Language

![Julia Logo](https://julialang.org/assets/infra/logo.svg)
"""

# ╔═╡ de79e82c-8e8f-11eb-0513-c7df38cb3d30
md"""
# Why Julia?
"""

# ╔═╡ 03a8af1e-8e90-11eb-2d94-e52e8c2be959
md"""
# Interpreted Language: Pros

- Simple: no need to learn the concept of memory, type system, etc.
- Flexible: almost everything can be modified at runtime
- ...
"""

# ╔═╡ 61a2883a-8e90-11eb-26bb-49925c4ad4e9
md"""
# Interpreted Language: Cons

- Slow: hard to optimize, a compiler won't be able to gain enough information for optimization
- Less error proof: flexibility make it hard to check errors
- ...
"""

# ╔═╡ b039ac32-8e90-11eb-08ec-8710712591f7
md"""
# Static Compiled Language: Pros

- Fast: designed to be compiled, easier to optimize
- Error Proof: the static type system can help programmer detect more potential errors
- ...
"""

# ╔═╡ d9287bda-8e90-11eb-1c46-af0555ed5158
md"""
# Static Compiled Language: Cons

- Complex: usually a static language has more concepts and more complicated type system comparing to interpret language
- Not Flexible: not everything can be changed at runtime, makes it less convenient to prototype with
- ...
"""

# ╔═╡ c763fe62-90d7-11eb-2483-31a70c8487cc
md"""
## The two language problem

In traditional software development, one usually uses a static compiled language along with an interpreted language
to get both efficiency and flexibility, e.g most of the python projects uses C++ and Python.

However, using two languages creates a lot of problems:

1. the project will be more complicated: now you need to learn two languages to develop it
2. the project will be less generic than one language: how these two languages communicate is via C ABI, which has no generic type information (we will learn what is this later)
3. when someone write a prototype using the interpreted language, it can be a lot extra work to convert it to the static compiled language for efficiency.
"""

# ╔═╡ 232da2be-8e91-11eb-3289-0d9bf19a225e
md"""
## Just-in-time compilation: compile right before you run it!

One solution to accelerate interpreted language is via just-in-time compilation.

*some runtime values become compile time*

![](https://source.android.com/devices/tech/dalvik/images/jit-profile-comp.png)

However, a language designed for interpreter is usually hard to accelerate and compile in general.
"""

# ╔═╡ e17692c6-8e91-11eb-3573-c98d429a5403
md"""
## Combine Static Compiled Language with Interpreted Language

*a language designed for just-in-time compilation*

[Julia: A Fresh Approach to
Numerical Computing (DOI:10.1137/141000671)](https://julialang.org/assets/research/julia-fresh-approach-BEKS.pdf)
"""

# ╔═╡ 363accfa-8e92-11eb-1a5b-6fac3df0290c
md"""

## Julia: A Differentiable Programmable Language

*Julia is currently the only option for whole language differentibility with rich numerical ecosystem*

[The Software 2.0](https://karpathy.medium.com/software-2-0-a64152b37c35) - Andrej Karpathy
"""

# ╔═╡ c9959dd6-8e92-11eb-04ed-7b8ee3717c28
Twitter(1152288749243420672)

# ╔═╡ 9bfff004-8e8f-11eb-0a21-1db0b13d1b12
md"""
# Before we start

- how do I setup my Julia environment?
- where can I get help?
"""

# ╔═╡ 7028360e-8e93-11eb-0ab6-3f50c4f864a5
md"""
# Download Julia: julialang.org
"""

# ╔═╡ 923c58a6-8e93-11eb-15cd-f7cac665afad
html"""

"""

# ╔═╡ 74e8356e-8e97-11eb-1d7f-29b54ca49dc4
html"""

"""

# ╔═╡ 30267782-8e98-11eb-368a-15bac39281cf
md"""
## Download Editor

- Visual Studio Code (Recommended): [https://code.visualstudio.com/](https://code.visualstudio.com/)
- Atom: [https://atom.io/](https://atom.io/)
"""

# ╔═╡ 65c1388e-8e98-11eb-036e-6be0f9ceeae8
md"""
## Install Plugin (VS Code)
"""

# ╔═╡ 8f7710a2-8e98-11eb-02af-69a34e81c9cf
html"""

"""

# ╔═╡ c2b5415e-8e99-11eb-2cbf-c1392fb5075e
md"""
# Install your first Julia package: the Pluto Notebook

most of the meterials of this lecture will be written as Pluto notebooks, so you will need to install Pluto notebooks, this is also the first Julia package you will install.
"""

# ╔═╡ 100a2712-8e9a-11eb-057b-cf84805d9f82
html"""
please open Julia's interactive session (known as REPL) and press ] key in the REPL to use the package mode, then type the following command
"""

# ╔═╡ 1c9fec26-8e9a-11eb-1f50-3512e5f1a36f
md"""
```
(@v1.6) pkg> add Pluto
```
"""

# ╔═╡ 6fb8a292-8e9a-11eb-3614-2bac1c3c7ebc
html"""
then press  (backspace) button and run the following command
"""

# ╔═╡ 767c63dc-8e9a-11eb-299d-257d7d961df2
md"""
```julia
julia> using Pluto; Pluto.run()
```

you will see a webpage pops out in your browser, then you can create a new notebook to play with Julia
"""

# ╔═╡ bd5bc9d6-8e98-11eb-2bb4-4bc600b4e887
md"""
## Where do I get help?

- checkout the docstring!
- read the official documentation: [https://docs.julialang.org/en/v1/](https://docs.julialang.org/en/v1/)
- ask your TA/instructor

besides above, the Julia programming language has an active community:

1. write a post on the discourse: [https://discourse.julialang.org/](https://discourse.julialang.org/)
2. join Julia Zulip chat group: [https://julialang.zulipchat.com/register/](https://julialang.zulipchat.com/register/)
3. join Julia Slack chat group: [https://julialang.org/slack/](https://julialang.org/slack/)
"""

# ╔═╡ d3df0b42-8fff-11eb-0fab-2bd5070ca3ff
md"""
# Where do I find packages?

- checkout [juliahub.com](https://juliahub.com/ui/Home)
- Google your keywords along with JuliaLang
- search on Github

# Where are the package documentations?

- (again) checkout [juliahub.com](https://juliahub.com/ui/Home)
- on the package README
- on the package webpage
"""

# ╔═╡ 9d29f678-8e99-11eb-2fac-595bd3335f64
md"""
# The Julia Basics
"""

# ╔═╡ 51920c44-8e9b-11eb-2ffa-6bb90b32cb76
md"""
## Variables

you can create a variable like other scripting language by assigning a value to it
"""

# ╔═╡ 65b492d2-8e9b-11eb-1f7e-6ff3d12b2acc
x = 1

# ╔═╡ 83192d6a-8e9b-11eb-2b3a-e7d120398f48
md"""
every variable has a type, you can check it using `typeof`
"""

# ╔═╡ b3e93bda-8e9b-11eb-177f-c3e6006a3409
typeof(x)

# ╔═╡ b8dfa5c0-8e9b-11eb-3378-51fe3adafa47
md"""
By default Julia displays the output of the last operation. (You can suppress the output by adding ; (a semicolon) at the end.)
"""

# ╔═╡ c662d416-8e9b-11eb-3309-abf865336ab6
y = 1;

# ╔═╡ ccfa44fa-8e9b-11eb-2d50-79ad7c21f35e
md"""
### Functions
In Julia, you can also define short-form, one-line functions using = (assignment) similar to how you write things mathematically.
"""

# ╔═╡ dd27a734-8e9b-11eb-10cd-697750fc8d3c
f(x) = 2x

# ╔═╡ e8ae462e-8e9b-11eb-2910-b38ea3de5252
md"""
Typing the function's name gives information about the function. To call it we must use parentheses:
"""

# ╔═╡ faa11e56-8e9b-11eb-0a29-2b33770da084
md"""
For longer functions we use the following syntax with the function keyword and end:
"""

# ╔═╡ 018a5106-8e9c-11eb-0d13-a5542c3fcd4b
function g(x, y)
	z = x + y
	return z^2
end

# ╔═╡ 0a93166e-8e9c-11eb-2fd5-1d80e13be555
md"""
# Control Flows
In Julia, there are for, if and while, they look like the following
"""

# ╔═╡ 1645283c-8e9c-11eb-204c-d93911df6dfa
begin
	s = 0
	for i in 1:10
		s += 1
	end
end

# ╔═╡ 1f7ca574-8e9c-11eb-19f8-931e6de060ed
md"""
here `begin...end` is used to wrap a code block as Pluto required, you don't have to write it outside Pluto, we can now check the value of s by typing it again
"""

# ╔═╡ 2db39ddc-8e9c-11eb-302a-890b52609e00
s

# ╔═╡ 33cddbd0-8e9c-11eb-02cc-5fbfbc0f4037
md"""
Here, 1:10 is a **range** representing the numbers from 1 to 10:
"""

# ╔═╡ 3c08b980-8e9c-11eb-3b35-3d8f254d4e4f
typeof(1:10)

# ╔═╡ 44cfe19e-8e9c-11eb-311b-3119586254e7
md"""
the if else statement looks like the following
"""

# ╔═╡ 4b2e285a-8e9c-11eb-3fb2-23dc7a5652d8
if s < 10
	# do something
elseif 10 < s < 13
	# do something
else
	# do something
end

# ╔═╡ f5c8c320-8f45-11eb-37c9-9b6e20cef2fb
md"""
## Modules

Modules is a way to define namespace. You can declare a Julia module using the `module` keyword, this is useful when you want to distinguish different project definitions etc.
"""

# ╔═╡ 0f61adba-8f46-11eb-16ff-6b4a9207f2d8
module MyModule

foo(x) = 1
goo(x) = 2

end

# ╔═╡ 19c77032-8f46-11eb-0c6d-d70bf3de9e87
module MyOtherModule

foo(x) = 1
goo(x) = 2

end

# ╔═╡ 2519ae78-8f46-11eb-214f-8324e9ad80c7
MyModule.foo === MyOtherModule.foo

# ╔═╡ 2e7b2654-8f46-11eb-3dfd-cb9a705f1b3c
md"here `===` is used to check if `foo` is identical. We can overload a function in a specific module using the dot operator"

# ╔═╡ 53902570-8f46-11eb-23bc-8f97b854a025
MyModule.foo(x::Int) = 3

# ╔═╡ 6a3bd774-8f46-11eb-2159-211c4128056a
MyModule.foo

# ╔═╡ 6d5a71a4-8f46-11eb-2277-15a035dcca14
MyOtherModule.foo

# ╔═╡ 216e5eda-8f47-11eb-05ea-915fa17e1f80
md"""
You can import the names from a module using the `using` keyword
"""

# ╔═╡ d46fa57a-8f47-11eb-2310-23a5f3ea7a8d
md"here the dot means to import a module from local scope. In fact, Julia packages are also Julia modules, but in global scope, that's why you don't need the dot for packages"

# ╔═╡ 53388a2c-8e9c-11eb-3ceb-7100da65a4a0
md"""
## Matrix and Array
Julia carries its own Array type, if you use Python, it is similar to `numpy.array` in Python except:

- index starts from `1`

- the multi-dimensional index is column-wise

You can also have list comprehension:
"""

# ╔═╡ 65cf20ae-8e9c-11eb-20f1-9dd4f08a8103
[i for i in 1:10]

# ╔═╡ 6aa77704-8e9c-11eb-22ac-bdbd3d23583b
md"""
it works for multi-dimensional case too:
"""

# ╔═╡ 70be63f2-8e9c-11eb-35b6-3b8f3bffa0f4
[(i, j) for i in 1:10, j in 1:5]

# ╔═╡ 7722da64-8e9c-11eb-3c4d-d98f275cb706
md"""
most builtin functions follow the same convention as numpy or MATLAB, e.g you can create a random matrix using:
"""

# ╔═╡ 8166d48c-8e9c-11eb-2032-cb423df0bd61
rand(5, 5)

# ╔═╡ 86253606-8e9c-11eb-271c-9950a3ffcd90
md"""
if you have question about using a function, you can always type question mark `?` in your REPL following the function name

```julia
julia> ?rand
```
"""

# ╔═╡ 7c354c6c-90e0-11eb-1125-d13c2252863c
md"""
## Broadcast

Sometimes, we want to apply a function element-wise to a matrix, in Julia you can do it on ANY function you defined using the broadcast syntax
"""

# ╔═╡ 96065b90-90e0-11eb-1041-e7654ab0b592
B = rand(5, 5)

# ╔═╡ 9c66f4f4-90e0-11eb-3512-e9b6e9bd838f
sin.(B)

# ╔═╡ a4433278-90e0-11eb-3cb7-9b8fcbac8b36
cos.(B)

# ╔═╡ b7c10a02-90e0-11eb-2098-9f26461cd856
md"works for multi-arg function too"

# ╔═╡ ade6e702-90e0-11eb-1b49-e7deedaafe70
B .* B

# ╔═╡ 775b7928-8f46-11eb-1fd1-613158bb44c5
md"""
# Some useful standard libararies

Since Julia is a general purposed programming language, certain functionality needs to be imported via the corresponding standard library. Here are some frequently used stdlib in numerics:

- `LinearAlgebra`: provides functions like `eigen`, `qr`, `svd` etc.
- `Statistics`: provides functions like `mean`, `var` etc.
- `SparseArrays`: provides the sparse matrix support

to find out other standard libraries, checkout the official Julia documentation
"""

# ╔═╡ bc6678a0-8e9d-11eb-006b-4766c2f47db9
md"""
## Multiple Dispatch

Julia functions is overloadable, we can put type annotations to distinguish different kinds of Julia functions
"""

# ╔═╡ dec6a71c-8e9d-11eb-29a4-cdfd1c27c01b
f(x::Int) = x + 1

# ╔═╡ e43bfe90-8e9d-11eb-3f8d-d5fd9b5d96cc
f(x::Float64) = x + 2

# ╔═╡ f518b756-8e9d-11eb-2153-3f8b7bad9e26
md"""
this works on arbitrary number of arguments
"""

# ╔═╡ fe2103dc-8e9d-11eb-2b6e-e73204feef2b
f(x, y::Int) = x + y

# ╔═╡ eecc7206-8e9b-11eb-0b13-b7f613eb974d
f

# ╔═╡ f2ba5414-8e9b-11eb-2941-277569b5be73
f(2)

# ╔═╡ ead8e42a-8e9d-11eb-3205-c1f58cac9b34
f(1)

# ╔═╡ ed7a9ade-8e9d-11eb-0793-379a695cfbdd
f(1.0)

# ╔═╡ 034b9626-8e9e-11eb-19e9-69e8979b480f
f(1, 2)

# ╔═╡ 0c1aced0-8e9e-11eb-0915-c9b5887e6beb
md"""
## Composite Types

Composite types are composition of other type of data, we can define our own composite types to create an abstraction
"""

# ╔═╡ 27df65ec-8e9e-11eb-3942-efa2aa12ea47
struct Person
	name::String
	age::Int
end

# ╔═╡ 31b00b96-8e9e-11eb-2f97-35ba8a796dfe
md"""
here `struct` creates an **immutable** type called `Person` that contains two fields `name` and `age`, by **immutable** we mean the instance of such type cannot mutate is fields.
"""

# ╔═╡ 48ed15ea-8e9e-11eb-1c6e-c5dc33a5b64f
me = Person("Roger", 25)

# ╔═╡ ab86757a-8e9e-11eb-2cdb-b382a39980c8
md"""
If your TA Ejaaz wants to change the name to his, it will error:
"""

# ╔═╡ 70b2eeec-8e9e-11eb-2d1b-bd58953e23fa
me.name = "Ejaaz"

# ╔═╡ 7fe242f0-8e9e-11eb-112f-89a3e3f9deff
md"""
This is useful to prevent programmers from mutating certain fields by accident, it is always a good practice to declare things to be **immutable** if the object won't be mutated

However, if you need to mutate the object afterwards, you can declare it as
"""

# ╔═╡ d4f2e470-8e9e-11eb-2968-d9df58079088
mutable struct Salary
	amount::Real
end

# ╔═╡ fb3f5564-8e9e-11eb-1e2d-65eeb023c79a
mysalary = Salary(120.0)

# ╔═╡ 0628776c-8e9f-11eb-05b5-e5708ebb2f81
mysalary.amount = 300

# ╔═╡ 0d339d8e-8e9f-11eb-3c3d-a36ce456314f
mysalary

# ╔═╡ 1ac4a3c8-8e9f-11eb-0d59-4beb79165dce
md"""
# Abstract Types

Sometimes you may have a few types that have similar properties, then you could define an abstract type for them, and let your concrete types to subtype from it
"""

# ╔═╡ 219876ea-8ea0-11eb-0be5-877e04a796b2
abstract type Animal end

# ╔═╡ 73a7dda4-8ea0-11eb-2a21-319c04b4e76d
md"""
and we can have a `Cat` type that is a subtype of `Animal`
"""

# ╔═╡ 80d86da2-8ea0-11eb-21e0-13174e6914ea
struct Cat <: Animal
	name::String
end

# ╔═╡ 84cd7b02-8ea0-11eb-3877-6bf6ed736336
md"""
here the operator `<:` means "subtype of" and similarly a `Dog` type
"""

# ╔═╡ 9eb515a4-8ea0-11eb-3907-19a21e6a3284
struct Dog <: Animal
	name::String
	belongsto::Person
end

# ╔═╡ a9ff885c-8ea0-11eb-1e04-e9863bdc2ad5
md"""
and we know all the `Animal`s speak something, so we can declare a function `speak` without a method implemented first
"""

# ╔═╡ bfcc792e-8ea0-11eb-314e-3db7e02e4fea
function speak end

# ╔═╡ d229fd44-8ea0-11eb-2e6b-bdd201d78923
md"""
we know how `Cat` speaks "meow"
"""

# ╔═╡ dcc95a10-8ea0-11eb-0313-151c0569da8d
speak(::Cat) = "meow"

# ╔═╡ ea530fa0-8ea0-11eb-0848-754126df85a6
md"""
and dog barks
"""

# ╔═╡ ff3a9398-8ea0-11eb-0e8e-e77a7d5ca2c0
speak(::Dog) = "bark"

# ╔═╡ 0d436aa0-8ea1-11eb-289a-191346b662c8
md"""
and we assume all the subtype of `Animal` has a field `name`, if a specific type does not have the field `name`, we can later overload the `name` function for it.
"""

# ╔═╡ 20c63966-8ea1-11eb-17dc-c7099ed1b239
name(x::Animal) = x.name

# ╔═╡ 39d6f102-8ea1-11eb-26b8-3bfd382b9347
md"""
now since we know the name of an animal and what it speaks, without knowing the specific object, we can define a method as following
"""

# ╔═╡ 8ef9d29e-8ea1-11eb-3997-695badb02d48
who_speaks_what(x::Animal) = "$(name(x)) speaks $(speak(x))"

# ╔═╡ f94dca4e-8ea1-11eb-3b04-5b2c79fc37bb
md"""
here, the `$` operator means string interpolation, we use it to insert the return value of `name(x)` and `speak(x)` into a string
"""

# ╔═╡ a0b81c16-8ea1-11eb-2f29-b14299f72021
who_speaks_what(Cat("Kitty"))

# ╔═╡ ab89224a-8ea1-11eb-3a3f-332bb716a224
md"""
this is called **generic programming**, we will later introduce its power in a concrete physics problem
"""

# ╔═╡ 39adb89c-8ea2-11eb-1246-7d550ddd737f
md"""
## Parameterized Types

we now are able to implement a function without knowing a concrete type by using an abstract type, however, what if we want to define an object without knowing its concrete type? e.g
the complex number `Complex` is the composition of two real numbers, the real number can be an integer, a floating point or even a symbolic value.

With this motivation, we introduce the parameterized type for defining **generic types**.

A Julia type can be parameterized using braces `{}`, and we can define a type variable `T` that is upbounded by `Real`, which means type `T` must be subtype of the abstract type `Real`
"""

# ╔═╡ a4d56714-8ea2-11eb-3d6d-4d0726879e74
struct MyNaiveComplex{T <: Real} <: Number
	real::T
	imag::T
end

# ╔═╡ 3571c524-8ea3-11eb-28e4-db4dddbaf8db
md"""
we can create a concrete object by filling the corresponding values
"""

# ╔═╡ 4849d286-8ea3-11eb-0d80-3f30d7d4ac50
MyNaiveComplex(1.0, 2.0)

# ╔═╡ 4e495a26-8ea3-11eb-0b17-d9d1256b53ca
md"the Julia compiler will automatically infer that `T` is actually `Float64`"

# ╔═╡ b273822e-8ea3-11eb-2f18-27307a2bd565
md"we can define some functions on this type without knowing what is `T`"

# ╔═╡ c0342cb0-8ea3-11eb-1f18-072f68fa57f4
myreal(x::MyNaiveComplex) = x.real

# ╔═╡ d0dafa26-8ea3-11eb-270e-33b9a1457e04
myimag(x::MyNaiveComplex) = x.imag

# ╔═╡ 5986463c-8ea4-11eb-334d-9d09ad119cdd
md"""
## Pretty Printing

as you have noticed, the printing of our custom types looks ugly, we can define our own pretty printing to make things look better
"""

# ╔═╡ 79e1f7d2-8ea4-11eb-2328-f186849eca74
Base.show(io::IO, x::MyNaiveComplex) = print(io, x.real, "+", x.imag, "im")

# ╔═╡ a5f35578-8ea4-11eb-2f5a-256eb27a49b0
md"here we overloads a builtin Julia function `show` in module `Base` to customize the printing for our own composite type"

# ╔═╡ 90bd7f9e-8ea4-11eb-2c52-e91178c599c3
MyNaiveComplex(1.0, 2.0)

# ╔═╡ 98898fa6-8ea4-11eb-334a-71a88a70f74d
md"in Julia, we use `im` to represent the imaginary unit"

# ╔═╡ f6f429e8-8ea4-11eb-2ae8-33005a4e501a
md"""
## Simple Symbolic Computing

now we know what is pretty printing, we can implement a very simple symbolic computing program

first we define an abstract type `SymExpr` for symbolic expression
"""

# ╔═╡ d9f4f198-8ea8-11eb-1edd-a1f051464aae
abstract type SymExpr end

# ╔═╡ ec32967c-8ea9-11eb-3eac-753f07be2e97
md"and we define a `Variable` type to represent a symbolic variable, it contains a field called `value`"

# ╔═╡ a3301890-8ea8-11eb-288c-b7d78bfdbf4b
struct Variable <: SymExpr
	value
end

# ╔═╡ 0e37132c-8eaa-11eb-22fc-9fb1db8f20be
md"we define a symbolic term type to represent composite symbolic terms, in general an symbolic term is in fact a function call, so we store the corresponding function and its arguements to represent a symbolic term"

# ╔═╡ ad7279ec-8ea8-11eb-080a-6fb902f16957
struct Term{F} <: SymExpr
	f::F
	args::Vector{Any}
end

# ╔═╡ 37424cbe-8eaa-11eb-0fe1-a5a33c53dddb
md"here we use a type parameter `F` to dispatch on the function later since we may need to handle the expression differently for different functions, e.g we can define a `+` function as following"

# ╔═╡ c010c8ce-8ea8-11eb-1d63-a3f0d2a406e8
Base.:(+)(x::SymExpr, y::SymExpr) = Term(+, Any[x, y])

# ╔═╡ 5b3136e4-8eaa-11eb-2773-55bd7d580e77
md"now we can try to add up some symbolic variables"

# ╔═╡ f253efd0-8ea8-11eb-0f0f-7792d68a57a2
Variable(:x) + Variable(:y) + Variable(:z)

# ╔═╡ 663303a6-8eaa-11eb-07f2-b9c74b621ef8
md"however this doesn't work if we have some constants"

# ╔═╡ 7d3f01da-8eaa-11eb-3812-fb6378a04d3f
Variable(:x) + Variable(:y) + 1

# ╔═╡ 8972ad58-8eaa-11eb-3c6b-e786c2764313
md"thus we need to overload for constant numbers, for simplicity, we only handle real numbers here。"

# ╔═╡ d1fb0f04-8ea8-11eb-1d9e-2baecef5a3d0
Base.:(+)(x::SymExpr, y::Real) = Term(+, Any[x, y])

# ╔═╡ b778ee60-8eaa-11eb-0984-47151aa5d7d3
md"because `+` operation on real numbers is communitive, we need to define it for the other order"

# ╔═╡ 0be349e8-8ea9-11eb-09fb-0b5d8c1f3642
Base.:(+)(x::Real, y::SymExpr) = Term(+, Any[x, y])

# ╔═╡ d81f8142-8eaa-11eb-2d8f-59f6393a3aa2
md"OK now this works!"

# ╔═╡ dd7dde40-8eaa-11eb-20fb-1541507c6bea
Variable(:x) + Variable(:y) + 1

# ╔═╡ a435b978-8eaa-11eb-08af-dfe9a3cecb46
md"now let's define a pretty printing to make it looks better"

# ╔═╡ 3a092504-8ea9-11eb-3117-033359cf9b73
Base.show(io::IO, x::Variable) = print(io, x.value)

# ╔═╡ 49b9d1a6-8ea9-11eb-2231-73d83734b41f
function Base.show(io::IO, x::Term{typeof(+)})
	for each in x.args[1:end-1]
		print(io, each, "+")
	end
	print(io, x.args[end])
end

# ╔═╡ 7281760c-8ea9-11eb-05b5-7d65a60c601d
Variable(:x) + Variable(:y) + Variable(:z) + 1

# ╔═╡ e5eb3738-8eaa-11eb-18f3-314a142a7429
md"""
the above program defines a very simple symbolic program, it does *no actual calculation* so far, but just record what we do using a tree data structure,
in some places, you may heard people call this tree structure as abstract/concrete syntax tree or algebra data type. We will later explore automatic differentiation
using this.
"""

# ╔═╡ 994d5a18-90df-11eb-0cbc-6df0aebb7e0a
md"""

## Generic Programming in Physics

![](https://github.com/Roger-luo/PIML/blob/master/notebooks/assets/tropical.png?raw=true)
"""

# ╔═╡ bc905b2c-90df-11eb-20cd-1395db871a99
md"""
## Implementing Tropical Numbers in Julia

you can define a tropical number in a few lines in Julia and use it with other packages (that makes sense mathematically)
"""

# ╔═╡ b771d9ea-90df-11eb-3e55-a9ecbd119cc7
begin
	struct Tropical{T} <: Number
	    n::T
	    Tropical{T}(x) where T = new{T}(T(x))
	    function Tropical(x::T) where T
	        new{T}(x)
	    end
	    function Tropical{T}(x::Tropical{T}) where T
	        x
	    end
	    function Tropical{T1}(x::Tropical{T2}) where {T1,T2}
	        new{T1}(T2(x.n))
	    end
	end

	Base.show(io::IO, t::Tropical) = Base.print(io, "$(t.n)ₜ")

	Base.:*(a::Tropical, b::Tropical) = Tropical(a.n + b.n)
	Base.:+(a::Tropical, b::Tropical) = Tropical(max(a.n, b.n))
end

# ╔═╡ e8187180-90df-11eb-043e-5b7847cc69fe
md"""
for a more complete definition of tropical numbers in Julia that handles corner cases, please refer to: [https://github.com/TensorBFS/TropicalNumbers.jl](https://github.com/TensorBFS/TropicalNumbers.jl)
"""

# ╔═╡ 30e07b24-90e0-11eb-249c-63f6d375cb86
A = Tropical.(rand(2, 2))

# ╔═╡ 406f2162-90e0-11eb-1d20-4936fcecad7a
A * A

# ╔═╡ 47d48e06-90e0-11eb-141c-b112c58c18db
md"""

### Question: How to do it in Python?

If you can't, why?
"""

# ╔═╡ 7aa0e478-90d7-11eb-1fc7-ddfc1c7d9284
md"""
## Write efficient Julia code

When you start caring about the efficiency, despite of algorithmic improvement and hardware improvement. There is a key concept in Julia
that you may frequently see.

In traditional software development, there is a separation between static compiled part and dynamic part by using two languages.
In julia the boundary of static compiled and dynamic interpreted is the **type stability**

In static compiled language, such as C++ (version later than C++11), one often needs to ensure the types can be correctly inferred.
This is true in Julia when you wish to check this, but the program will still work if the compiler fail to infer the type.
"""

# ╔═╡ 9d2f7258-8e9c-11eb-23ca-cf11c2a331d8
md"""
## Package management

Julia carries its own package manager, you can use it as a normal package:

```julia
julia> using Pkg
```

to install a pacakge, you can use

```julia
julia> Pkg.add("Yao")
```

to remove a pacakge, you can use

```julia
julia> Pkg.rm("Yao")
```
"""

# ╔═╡ e7064866-8e9c-11eb-0af6-85bdfd62fd8b
md"""
If you are using a REPL, there is a convenient Pkg mode to let you install packages directly. You can enter the Pkg mode by typing `]` key like how we install the Pluto package previously
"""

# ╔═╡ ec8c1dc4-8e9c-11eb-1881-7342f68b79a2
md"""
## Environments

All Julia program **runs inside an environment**, it is the global environment by default. It is usually recommended to run your notebook in a local environment, so we won't hit any version conflicts between different packages.
we will introduce this part more when we write our own automatic differentiation engine.
"""

# ╔═╡ ef361372-8f44-11eb-1fc7-21656813dfe9
md"""
## Jupyter Notebooks

Although, this course will be teached using Pluto notebooks and VS Code editor, you can also use Jupyter notebooks if you don't like Pluto notebooks.

*Fun fact: the name Jupyter actually means Julia, Python and R*

## Install Jupyter notebooks

you can install the Julia kernel for Jupyter notebook by install the `IJulia` package in your Julia REPL

```julia
pkg> add IJulia
```

and then you can either start your jupyter notebook using

```julia
julia> using IJulia

julia> notebook()
```

or from command line as you do for python, and choose the kernel to be Julia.
"""

# ╔═╡ 660c5dfa-8f45-11eb-0c2b-6de98f143ebb
md"""

# Learn More about the Julia programming language

- watch JuliaCon videos on YouTube: [https://www.youtube.com/channel/UC9IuUwwE2xdjQUT_LMLONoA](https://www.youtube.com/channel/UC9IuUwwE2xdjQUT_LMLONoA) or Bilibili: [https://space.bilibili.com/356692611?from=search&seid=3189679744055766688](https://space.bilibili.com/356692611?from=search&seid=3189679744055766688)
- join JuliaCon 2021! [https://juliacon.org/2021/](https://juliacon.org/2021/)
- check out the MIT computational thinking course: [https://computationalthinking.mit.edu/](https://computationalthinking.mit.edu/)
"""

# ╔═╡ Cell order:
# ╠═f1aa59ae-8e94-11eb-12bc-3b1d9ccc314b
# ╟─cc2cf6c0-8e92-11eb-1255-53074ed0e59c
# ╟─7400844c-8e8f-11eb-27d0-859ea6dba632
# ╟─de79e82c-8e8f-11eb-0513-c7df38cb3d30
# ╟─03a8af1e-8e90-11eb-2d94-e52e8c2be959
# ╟─61a2883a-8e90-11eb-26bb-49925c4ad4e9
# ╟─b039ac32-8e90-11eb-08ec-8710712591f7
# ╟─d9287bda-8e90-11eb-1c46-af0555ed5158
# ╟─c763fe62-90d7-11eb-2483-31a70c8487cc
# ╟─232da2be-8e91-11eb-3289-0d9bf19a225e
# ╟─e17692c6-8e91-11eb-3573-c98d429a5403
# ╟─363accfa-8e92-11eb-1a5b-6fac3df0290c
# ╟─c9959dd6-8e92-11eb-04ed-7b8ee3717c28
# ╟─9bfff004-8e8f-11eb-0a21-1db0b13d1b12
# ╟─7028360e-8e93-11eb-0ab6-3f50c4f864a5
# ╟─923c58a6-8e93-11eb-15cd-f7cac665afad
# ╟─74e8356e-8e97-11eb-1d7f-29b54ca49dc4
# ╟─30267782-8e98-11eb-368a-15bac39281cf
# ╟─65c1388e-8e98-11eb-036e-6be0f9ceeae8
# ╟─8f7710a2-8e98-11eb-02af-69a34e81c9cf
# ╟─c2b5415e-8e99-11eb-2cbf-c1392fb5075e
# ╟─100a2712-8e9a-11eb-057b-cf84805d9f82
# ╟─1c9fec26-8e9a-11eb-1f50-3512e5f1a36f
# ╟─6fb8a292-8e9a-11eb-3614-2bac1c3c7ebc
# ╟─767c63dc-8e9a-11eb-299d-257d7d961df2
# ╟─bd5bc9d6-8e98-11eb-2bb4-4bc600b4e887
# ╟─d3df0b42-8fff-11eb-0fab-2bd5070ca3ff
# ╟─9d29f678-8e99-11eb-2fac-595bd3335f64
# ╟─51920c44-8e9b-11eb-2ffa-6bb90b32cb76
# ╠═65b492d2-8e9b-11eb-1f7e-6ff3d12b2acc
# ╟─83192d6a-8e9b-11eb-2b3a-e7d120398f48
# ╠═b3e93bda-8e9b-11eb-177f-c3e6006a3409
# ╟─b8dfa5c0-8e9b-11eb-3378-51fe3adafa47
# ╠═c662d416-8e9b-11eb-3309-abf865336ab6
# ╟─ccfa44fa-8e9b-11eb-2d50-79ad7c21f35e
# ╠═dd27a734-8e9b-11eb-10cd-697750fc8d3c
# ╟─e8ae462e-8e9b-11eb-2910-b38ea3de5252
# ╠═eecc7206-8e9b-11eb-0b13-b7f613eb974d
# ╠═f2ba5414-8e9b-11eb-2941-277569b5be73
# ╟─faa11e56-8e9b-11eb-0a29-2b33770da084
# ╠═018a5106-8e9c-11eb-0d13-a5542c3fcd4b
# ╟─0a93166e-8e9c-11eb-2fd5-1d80e13be555
# ╠═1645283c-8e9c-11eb-204c-d93911df6dfa
# ╟─1f7ca574-8e9c-11eb-19f8-931e6de060ed
# ╠═2db39ddc-8e9c-11eb-302a-890b52609e00
# ╟─33cddbd0-8e9c-11eb-02cc-5fbfbc0f4037
# ╠═3c08b980-8e9c-11eb-3b35-3d8f254d4e4f
# ╟─44cfe19e-8e9c-11eb-311b-3119586254e7
# ╠═4b2e285a-8e9c-11eb-3fb2-23dc7a5652d8
# ╟─f5c8c320-8f45-11eb-37c9-9b6e20cef2fb
# ╠═0f61adba-8f46-11eb-16ff-6b4a9207f2d8
# ╠═19c77032-8f46-11eb-0c6d-d70bf3de9e87
# ╠═2519ae78-8f46-11eb-214f-8324e9ad80c7
# ╟─2e7b2654-8f46-11eb-3dfd-cb9a705f1b3c
# ╠═53902570-8f46-11eb-23bc-8f97b854a025
# ╠═6a3bd774-8f46-11eb-2159-211c4128056a
# ╠═6d5a71a4-8f46-11eb-2277-15a035dcca14
# ╟─216e5eda-8f47-11eb-05ea-915fa17e1f80
# ╠═3fcaf0b4-8f47-11eb-2145-311f5d76757f
# ╟─d46fa57a-8f47-11eb-2310-23a5f3ea7a8d
# ╟─53388a2c-8e9c-11eb-3ceb-7100da65a4a0
# ╠═65cf20ae-8e9c-11eb-20f1-9dd4f08a8103
# ╟─6aa77704-8e9c-11eb-22ac-bdbd3d23583b
# ╠═70be63f2-8e9c-11eb-35b6-3b8f3bffa0f4
# ╟─7722da64-8e9c-11eb-3c4d-d98f275cb706
# ╠═8166d48c-8e9c-11eb-2032-cb423df0bd61
# ╟─86253606-8e9c-11eb-271c-9950a3ffcd90
# ╟─7c354c6c-90e0-11eb-1125-d13c2252863c
# ╠═96065b90-90e0-11eb-1041-e7654ab0b592
# ╠═9c66f4f4-90e0-11eb-3512-e9b6e9bd838f
# ╠═a4433278-90e0-11eb-3cb7-9b8fcbac8b36
# ╟─b7c10a02-90e0-11eb-2098-9f26461cd856
# ╠═ade6e702-90e0-11eb-1b49-e7deedaafe70
# ╟─775b7928-8f46-11eb-1fd1-613158bb44c5
# ╟─bc6678a0-8e9d-11eb-006b-4766c2f47db9
# ╠═dec6a71c-8e9d-11eb-29a4-cdfd1c27c01b
# ╠═e43bfe90-8e9d-11eb-3f8d-d5fd9b5d96cc
# ╠═ead8e42a-8e9d-11eb-3205-c1f58cac9b34
# ╠═ed7a9ade-8e9d-11eb-0793-379a695cfbdd
# ╟─f518b756-8e9d-11eb-2153-3f8b7bad9e26
# ╠═fe2103dc-8e9d-11eb-2b6e-e73204feef2b
# ╠═034b9626-8e9e-11eb-19e9-69e8979b480f
# ╟─0c1aced0-8e9e-11eb-0915-c9b5887e6beb
# ╠═27df65ec-8e9e-11eb-3942-efa2aa12ea47
# ╟─31b00b96-8e9e-11eb-2f97-35ba8a796dfe
# ╠═48ed15ea-8e9e-11eb-1c6e-c5dc33a5b64f
# ╟─ab86757a-8e9e-11eb-2cdb-b382a39980c8
# ╠═70b2eeec-8e9e-11eb-2d1b-bd58953e23fa
# ╟─7fe242f0-8e9e-11eb-112f-89a3e3f9deff
# ╠═d4f2e470-8e9e-11eb-2968-d9df58079088
# ╠═fb3f5564-8e9e-11eb-1e2d-65eeb023c79a
# ╠═0628776c-8e9f-11eb-05b5-e5708ebb2f81
# ╠═0d339d8e-8e9f-11eb-3c3d-a36ce456314f
# ╟─1ac4a3c8-8e9f-11eb-0d59-4beb79165dce
# ╠═219876ea-8ea0-11eb-0be5-877e04a796b2
# ╟─73a7dda4-8ea0-11eb-2a21-319c04b4e76d
# ╠═80d86da2-8ea0-11eb-21e0-13174e6914ea
# ╟─84cd7b02-8ea0-11eb-3877-6bf6ed736336
# ╠═9eb515a4-8ea0-11eb-3907-19a21e6a3284
# ╟─a9ff885c-8ea0-11eb-1e04-e9863bdc2ad5
# ╠═bfcc792e-8ea0-11eb-314e-3db7e02e4fea
# ╟─d229fd44-8ea0-11eb-2e6b-bdd201d78923
# ╠═dcc95a10-8ea0-11eb-0313-151c0569da8d
# ╟─ea530fa0-8ea0-11eb-0848-754126df85a6
# ╠═ff3a9398-8ea0-11eb-0e8e-e77a7d5ca2c0
# ╟─0d436aa0-8ea1-11eb-289a-191346b662c8
# ╠═20c63966-8ea1-11eb-17dc-c7099ed1b239
# ╟─39d6f102-8ea1-11eb-26b8-3bfd382b9347
# ╠═8ef9d29e-8ea1-11eb-3997-695badb02d48
# ╟─f94dca4e-8ea1-11eb-3b04-5b2c79fc37bb
# ╠═a0b81c16-8ea1-11eb-2f29-b14299f72021
# ╟─ab89224a-8ea1-11eb-3a3f-332bb716a224
# ╟─39adb89c-8ea2-11eb-1246-7d550ddd737f
# ╠═a4d56714-8ea2-11eb-3d6d-4d0726879e74
# ╟─3571c524-8ea3-11eb-28e4-db4dddbaf8db
# ╠═4849d286-8ea3-11eb-0d80-3f30d7d4ac50
# ╟─4e495a26-8ea3-11eb-0b17-d9d1256b53ca
# ╟─b273822e-8ea3-11eb-2f18-27307a2bd565
# ╠═c0342cb0-8ea3-11eb-1f18-072f68fa57f4
# ╠═d0dafa26-8ea3-11eb-270e-33b9a1457e04
# ╟─5986463c-8ea4-11eb-334d-9d09ad119cdd
# ╠═79e1f7d2-8ea4-11eb-2328-f186849eca74
# ╟─a5f35578-8ea4-11eb-2f5a-256eb27a49b0
# ╠═90bd7f9e-8ea4-11eb-2c52-e91178c599c3
# ╟─98898fa6-8ea4-11eb-334a-71a88a70f74d
# ╟─f6f429e8-8ea4-11eb-2ae8-33005a4e501a
# ╠═d9f4f198-8ea8-11eb-1edd-a1f051464aae
# ╟─ec32967c-8ea9-11eb-3eac-753f07be2e97
# ╠═a3301890-8ea8-11eb-288c-b7d78bfdbf4b
# ╟─0e37132c-8eaa-11eb-22fc-9fb1db8f20be
# ╠═ad7279ec-8ea8-11eb-080a-6fb902f16957
# ╟─37424cbe-8eaa-11eb-0fe1-a5a33c53dddb
# ╠═c010c8ce-8ea8-11eb-1d63-a3f0d2a406e8
# ╟─5b3136e4-8eaa-11eb-2773-55bd7d580e77
# ╠═f253efd0-8ea8-11eb-0f0f-7792d68a57a2
# ╟─663303a6-8eaa-11eb-07f2-b9c74b621ef8
# ╠═7d3f01da-8eaa-11eb-3812-fb6378a04d3f
# ╟─8972ad58-8eaa-11eb-3c6b-e786c2764313
# ╠═d1fb0f04-8ea8-11eb-1d9e-2baecef5a3d0
# ╟─b778ee60-8eaa-11eb-0984-47151aa5d7d3
# ╠═0be349e8-8ea9-11eb-09fb-0b5d8c1f3642
# ╟─d81f8142-8eaa-11eb-2d8f-59f6393a3aa2
# ╠═dd7dde40-8eaa-11eb-20fb-1541507c6bea
# ╟─a435b978-8eaa-11eb-08af-dfe9a3cecb46
# ╠═3a092504-8ea9-11eb-3117-033359cf9b73
# ╠═49b9d1a6-8ea9-11eb-2231-73d83734b41f
# ╠═7281760c-8ea9-11eb-05b5-7d65a60c601d
# ╟─e5eb3738-8eaa-11eb-18f3-314a142a7429
# ╟─994d5a18-90df-11eb-0cbc-6df0aebb7e0a
# ╟─bc905b2c-90df-11eb-20cd-1395db871a99
# ╠═b771d9ea-90df-11eb-3e55-a9ecbd119cc7
# ╟─e8187180-90df-11eb-043e-5b7847cc69fe
# ╠═30e07b24-90e0-11eb-249c-63f6d375cb86
# ╠═406f2162-90e0-11eb-1d20-4936fcecad7a
# ╟─47d48e06-90e0-11eb-141c-b112c58c18db
# ╟─7aa0e478-90d7-11eb-1fc7-ddfc1c7d9284
# ╟─9d2f7258-8e9c-11eb-23ca-cf11c2a331d8
# ╟─e7064866-8e9c-11eb-0af6-85bdfd62fd8b
# ╟─ec8c1dc4-8e9c-11eb-1881-7342f68b79a2
# ╟─ef361372-8f44-11eb-1fc7-21656813dfe9
# ╟─660c5dfa-8f45-11eb-0c2b-6de98f143ebb

max_stars_repo_name: Roger-luo/PIML

max_issues_repo_name: Roger-luo/PIML

max_forks_repo_name: Roger-luo/PIML


text:
function process_default_favicon(request::HTTP.Request, state::HandlerState)
    ico_contents = read(
        joinpath(ROOT_PATH, "src", "favicon.ico")
    )
    return HTTP.Response(
        200,
        ["Content-Type" => "image/x-icon"],
        body = ico_contents
        )
end


text:
Base.@kwdef struct Gradient{T}
    directions::Matrix{T}
    sequences::Vector{TimeProfile}
    values::Vector{T}
    values_type::String
end


text:
#    Copyright (c) 2021 J.A. Duffek
#    Copyright (c) 2000 D.M. Spink
#
#    This program is free software; you can redistribute it and/or modify
#    it under the terms of the GNU General Public License as published by
#    the Free Software Foundation; either version 2 of the License, or
#    (at your option) any later version.
#
#    This program is distributed in the hope that it will be useful,
#    but WITHOUT ANY WARRANTY; without even the implied warranty of
#    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
#    GNU General Public License for more details.
#
#    You should have received a copy of the GNU General Public License
#    along with this program; if not, see <http://www.gnu.org/licenses/>.

using Plots;

"""
    demotorus()

A second demonstration of surface construction.

# Examples:
```julia
julia> demotorus()
```
"""
function demotorus()
sphere = nrbrevolve(nrbcirc(1.0),[1.0;0.0;0.0]);
nrbplot(sphere,[20;20],
        title = "Ball and torus - surface construction by revolution.",
        c=:copper,
        cbar=false);
torus = nrbrevolve(nrbcirc(0.2,[0.9;1.0]),[1.0;0.0;0.0]);
nrbplot!(torus,[20;10],c=:copper,cbar=false);
nrbplot!(nrbtform(torus,vectrans([-1.8])),[20;10],c=:copper,cbar=false);
# TODO: come back later as pyplot backend does not has light/shading yet
#       furthermore grid lines should be off and interpolated shading should be
#       on
end # demotorus


text:
function prefilterforolscontributions(df::DataFrame;
  minpointspluskforβ::Int = MIN_POINTS_PLUS_K_FOR_BETA,
  returnfield::Symbol = :lreturn,
  essentialfields::Vector{Symbol}=[:Llreturnwin, :pcontributions])::DataFrame

  local subdf::SubDataFrame
  local ssubdf::SubDataFrame

  #make sure we ahve enough unique points for each organization
  df[:useforregressions] = trues(size(df,1))

  checkfields::Vector{Symbol} = [returnfield; essentialfields; :useforregressions]

  #loop over all unique eins
  for subdf ∈ groupby(df, :ein)

    ssubdf = view(subdf, completecases(subdf[checkfields]), checkfields)
    subN = minimum((f->size(unique(ssubdf[f]),1)).(essentialfields))

    #make sure we have enough points for the regression
    if subN < minpointspluskforβ + 2 #add two for the existing dof
      ssubdf[:useforregressions] .= false
    end
  end

  return df[df[:useforregressions],:]
end

function makeolscontributiontables(data::NCCSData; rerunprep::Bool = true)

  #NOTE: This is also called in the compute function
  #if nothing is changed it does not need to be run

  #can disable the below

  local displayedvars::Vector{Symbol}

  #this gets repeated a lot so we can save time doing it now
  regdf::DataFrame = prefilterforolscontributions(data.df)

  #Univariate Regressions
  displayedvars = [:LagReturn, :LagReturn3yr]
  olspanel(regdf, suffix = "contribv2",
    displayedvars = displayedvars,
    contentRowNames = (string).(displayedvars),
    regressionfunction = olsretcontributionL13,
    titleCaption = "Contributions vs Lagged Returns")

  displayedvars = [:LagContrib, :LagExp]
  olspanel(regdf, suffix = "retcon",
    displayedvars = displayedvars,
    contentRowNames = (string).(displayedvars),
    regressionfunction = olscontributionregretcon,
    titleCaption = "Returns vs Lagged Contributions")

  displayedvars = [:LagReturn, :LagReturn3yr, :LagExp]
  olspanel(regdf, suffix = "conret",
    displayedvars = displayedvars,
    contentRowNames = (string).(displayedvars),
    regressionfunction = olscontributionregconret,
    titleCaption = "Returns vs Lagged Contributions")

  #expense regressions
  #=displayedvars = [:LagReturn, :Lag2xReturn, :LagReturn3yr, :Lag2xReturn3yr]
  olspanel(regdf, suffix = "contrib2xlagv2",
    displayedvars = displayedvars,
    contentRowNames = (string).(displayedvars),
    regressionfunction = olsretcontributionLL13,
    titleCaption = "Contributions vs Twice-lagged Returns")=#

end

function olscontributiontables(data::NCCSData)
  makeolspaneltables(data)
end


text:
using Distributed
addprocs(5; exeflags=["--project"])
@everywhere using RenyiEntropyOfMagic
using Serialization

mc_M2 = pmap(0.1:0.1:2.5) do h
  @info "running" h
  tfim_renyi_entropy(20, h)
end
serialize("20.dat", mc_M2)

# using Yao
# L = 8
# tfim_h = -sum(kron(L, i=>Z, mod1(i+1, L)=>Z) for i in 1:L) - sum(put(L, i=>X) for i in 1:L)
# mat(Z)
# H = Matrix(mat(tfim_h))
# using LinearAlgebra

# eigmin(H)
# eigvals(H)
# psi = eigvecs(H)[:, 1]

# Os = enumerate_paulistring(L)

# -log(sum((psi' * kron(O...) * psi)^4 for O in Os)/2^L)/L

# exact_tfim_renyi_entropy(8, 1.0)

# As = get_matrices(8, 1.0)
# exact_renyi_entropy(As)

# renyi_entropy(As; nsamples=10^6)

# mc_M2 = [tfim_renyi_entropy(20, h) for h in 0.1:0.1:2.5]

# exact_M2 = [exact_tfim_renyi_entropy(290, h) for h in 0.1:0.1:2.5]

# tfim_renyi_entropy(20, 1.0)
# exact_tfim_renyi_entropy(8, 0.1)

# using CairoMakie

# lines(exact_M2)

# i = Index(2, "S=1/2")
# j = Index(2, "S=1/2")
# println(ITensors.op("X", i))


text:
## Minimum diagonal adjustment to a matrix to make it
## positive-definite, preserving the off-diagonal values.
#
# The problem is to factor Q into C*C', where each column of C
# has at most two non-zero entries. Consequently we have the
# equation
#
#     C_{ik}*C_{jk} = Q_{ij}
#
# where k is the column of C corresponding to (i,j).
# In addition to these paired columns, we also allow single-entry
# columns with the value c_i (note lowercase). Then we can write
#
#    (C*C')_{ii} = c_i^2 + \sum_k C_{ik}^2
#
# and we seek to minimize
#
#  P = \sum_i (C*C')_{ii}/|Q_{ii}|
#
# over the i with Q_{ii} != 0, subject to the constraint
#
#  (C*C')_{ii} >= |Q_{ii}|.
#
# Unlike the straight numerator, the scaled version used in P is
# invariant with respect to physical units.
#
# This minimization can be done analytically by introducing a ratio
#
#     r_k = C_{ik}/C_{jk}
#
# and computing the gradient of P wrt r.

"""
    Qp = possemidef(Q)

Return a variant of Q preserving the off-diagonal values and adjusting
the diagonals, as needed, to make the matrix positive semi-definite. Any
NaN entries of Q are treated as 0. Q is assumed to be symmetric and
only the lower triangle is used.
"""
function possemidef(Q::AbstractMatrix)
    ind1, ind2 = axes(Q)
    ind1 == ind2 || error("Matrix must be square")
    Qp = copy(Q)
    for i in ind1
        Qp[i,i] = 0
    end
    for j in ind2
        csum = zero(eltype(Q))
        qjj = abs(Q[j,j])
        for i = j+1:last(ind1)
            q = Q[i,j]
            ((q == 0) | isnan(q)) && continue
            qii = abs(Q[i,i])
            r = sign(q) * sqrt(qjj/qii)
            if r == 0 || !isfinite(r)
                r = oftype(r, sign(q))
            end
            Qp[j,j] += q * r
            Qp[i,i] += q / r
        end
    end
    for i in ind1
        # this is equivalent to adding c_i^2
        Qp[i,i] = max(Qp[i,i], abs(Q[i,i]))
    end
    return Qp
end


text:
function transition(time::Int64, system::Array{Int64,2}, parameters::Dict{String,Any})::Array{Float64,2}

    # initialize -
    T = zeros(2,2)
    A = system[time,1]
    B = system[time,2]

    # get parameters -
    K1 = parameters["K1"]
    n1 = parameters["n1"]
    K2 = parameters["K2"]
    n2 = parameters["n2"]

    # compute -
    p1 = A^n1/(K1^n1+A^n1)
    p2 = B^n2/(K2^n2+B^n2)

    # package -
    T[1,1] = 1 - p1
    T[1,2] = p1
    T[2,1] = p2
    T[2,2] = 1 - p2

    # return -
    return T
end

function machine(time::Int64, state::Int64, system::Array{Int64,2}, parameters::Dict{String,Any})::Int64

    # compute the probabilities -
    T = transition(time, system, parameters)

    # generate a random number -
    r = rand() # uniform from 0 -> 1

    # what are my choice?
    choices = 1.0 .- T[state,:]

    # what is the index of the max?
    idx_max = argmax(choices)
    p = choices[idx_max]

    if (r >= p)
        return argmax(choices)
    else
        return argmin(choices)
    end
end


text:
# This is the value of standardize used for the vast majority of tests:
const _standard = Ref(true)

if !_standard[]
    @info "Suppressing 2x2 block checks pending standardization fixes."
end

# Check that 2x2 diagonal blocks have standard form.
# This requires that tiny values have been cleaned out.
function checkblocks(A::StridedMatrix; debug=false)
    n = size(A,1)
    isok = true
    for j=1:n-1
        if A[j+1,j] != 0
            isok &= (A[j,j] == A[j+1,j+1])
            isok &= (A[j,j+1] != 0)
            isok &= (sign(A[j,j+1]) * sign(A[j+1,j]) < 0) # signs must differ (overflow?)
        end
        if debug && !isok
            print("problem block: ",A[j:j+1,j:j+1])
            (A[j,j] == A[j+1,j+1]) || print(" diag")
            (A[j,j+1] != 0) || print(" 0_super")
            (sign(A[j,j+1]) * sign(A[j+1,j]) < 0) || print(" trail_sign")
            println()
            return false
        end
    end
    isok
end

# This only applies to standard forms w/ co-sorted eigenvals
function checkeigvals(A::StridedMatrix{T}, w::StridedVector, tol; debug=false
                      ) where T <: Real
    n = size(A,1)
    ulp = eps(T)
    tiny = GenericSchur.safemin(T)
    isok = true
    for j=1:n
        isok &= (A[j,j] == real(w[j]))
    end
    if n>1
        if A[2,1] == 0
            isok &= (imag(w[1]) == 0)
        end
        if A[n,n-1] == 0
            isok &= (imag(w[n]) == 0)
        end
    end
    for j=1:n-1
        if A[j+1,j] != 0
            t = sqrt(abs(A[j+1,j]))*sqrt(abs(A[j,j+1]))
            cmp = max(ulp*t,tiny)
            isok &= (abs(imag(w[j]) - t) / cmp < tol)
            isok &= (abs(imag(w[j+1]) + t) / cmp < tol)
        elseif (j>1) && (A[j+1,j] == 0) && (A[j,j-1] == 0)
            isok &= (imag(w[j]) == 0)
        end
    end
    isok
end

function schurtest(A::Matrix{T}, tol; normal=false, standard=_standard[],
                   baddec=false,
                   ) where {T<:Real}
    n = size(A,1)
    ulp = eps(T)
    if T <: BlasFloat || (standard != GenericSchur._STANDARDIZE_DEFAULT)
        S = GenericSchur.gschur(A, standardize=standard)
    else
        # FIXME: no keywords allowed here; thanks, LinearAlgebra
        S = schur(A)
    end
    # test 1: S.T is upper quasi-triangular
    @test all(tril(S.T,-2) .== 0)
    # check complex 2x2 blocks
    if standard
        @test checkblocks(S.T)
    end
    # test 2: norm(A - S.Z * S.T * S.Z') / (n * norm(A) * ulp) < tol
    if baddec
        @test_broken norm(A - S.Z * S.T * S.Z') / (n * norm(A) * ulp) < tol
    else
        @test norm(A - S.Z * S.T * S.Z') / (n * norm(A) * ulp) < tol
    end
    # test 3: S.Z is orthogonal: norm(I - S.Z * S.Z') / (n * ulp) < tol
    @test norm(I - S.Z * S.Z') / (n * ulp) < tol
    # test 4: S.values are e.v. of T
    if standard
        @test checkeigvals(S.T, S.values, tol)
    end

    # It is tempting to check eigenvalues against LAPACK eigvals(A)
    # for suitable types, but the comparison is misleading without
    # detailed analysis because many of the test matrices are severly
    # non-normal.

    if normal
        # verify that Schur "diagonalizes" normal matrices
        # necessary but not really sufficient (complex blocks)
        @test norm(triu(S.T,2)) / (n * norm(A) * ulp) < tol
        # TODO: in this case comparison to LAPACK eigenvalues is legit
    end

    if standard
        Sc = triangularize(S)

        # test 1: S.T is upper triangular
        @test all(tril(Sc.T,-1) .== 0)
        # test 2: norm(A - S.Z * S.T * S.Z') / (n * norm(A) * ulp) < tol
        @test norm(A - Sc.Z * Sc.T * Sc.Z') / (n * norm(A) * ulp) < tol
        # test 3: S.Z is unitary: norm(I - S.Z * S.Z') / (n * ulp) < tol
        @test norm(I - Sc.Z * Sc.Z') / (n * ulp) < tol
        # test 4: S.values are e.v. of T
        @test all(csort(Sc.values) .== csort(diag(Sc.T)))
    end

end

# random orthogonal matrix
function rando(::Type{T},n) where {T<:Real}
    if T ∈ [Float16,Float32,Float64]
        A = randn(T,n,n)
        F = qr(A)
        return Matrix(F.Q)
    else
        # don't have normal deviates for other types, but need appropriate QR
        A = randn(Float64,n,n)
        F = qr(convert.(T,A))
        return Matrix(F.Q)
    end
end

Random.seed!(1234)

for T in [BigFloat,Float16]
    @testset "group $T" begin
        unfl = floatmin(T)
        ovfl = one(T) / unfl
        ulp = eps(T)
        ulpinv = one(T) / ulp
        rtulp = sqrt(ulp)
        rtulpi = one(T) / rtulp

        ens = [4,32]
        @testset "random general" begin
            for n in ens
                A = rand(T,n,n)
                schurtest(A,10)
            end
        end
        @testset "random orthogonal" begin
            for n in ens
                A = rando(T,n)
                schurtest(A,10,normal=true)
            end
        end

    end # group testset
end # type loop (non-BlasFloat)

@testset "Godunov" begin
    setprecision(BigFloat, 80) do
        A,v = godunov(BigFloat)
        S = schur(A)
        @test isapprox(csort(S.values),v,atol=2e-5)
        vnew = eigvals(A)
        @test isapprox(csort(vnew),v,atol=2e-5)
    end
end

for T in [Float64, Float32]
    @testset "group $T" begin

        unfl = floatmin(T)
        ovfl = one(T) / unfl
        ulp = eps(T)
        ulpinv = one(T) / ulp
        rtulp = sqrt(ulp)
        rtulpi = one(T) / rtulp

        # sizes of matrices to test
        ens = [4,32]

        itypemax = 10
        # On author's dev system, tests pass with tol=5 for most and 6 for all,
        # but generalization is risky.
        # (LAPACK uses 20, FWIW.)
        tols = fill(20,itypemax)

        simrconds = [1.0, rtulpi, 0.0]
        magns = [1.0, ovfl*ulp, unfl*ulpinv]

        @testset "Jordan" begin
            itype=3
            kmagn = [1]
            nj = length(kmagn)
            for n in ens
                for j in 1:nj
                    anorm = magns[kmagn[j]]
                    A = diagm(1 => ones(T,n-1)) + anorm * I
                    schurtest(A,tols[itype])
                end
            end
        end

        @testset "general, ev specified" begin
            itype=6
            kmagn =  [1,1,1,1,1,1,1,1,2,3]
            kmode =  [4,3,1,5,4,3,1,5,5,5]
            kconds = [1,1,1,1,2,2,2,2,2,2]
            nj = length(kmode)
            for n in ens
                for j in 1:nj
#               CALL DLATME( N, 'S', ISEED, WORK, IMODE, COND, CONE,
#     $                      ' ', 'T', 'T', 'T', RWORK, 4, CONDS, N, N, ANORM,
#     $                      A, LDA, WORK( 2*N+1 ), IINFO )
                    A = zeros(T,n,n)
                    simrcond = simrconds[kconds[j]]
                    imode = kmode[j]
                    anorm = magns[kmagn[j]]
                    rcond = ulpinv
                    verbosity[] > 1 &&
                        println("latme: n=$n $anorm, $imode, $rcond, $simrcond")
                    latme!(A,anorm,imode,rcond,simrcond)
                    schurtest(A,tols[itype])
                end
            end
        end

        @testset "diagonal, random ev" begin
            itype=7
            kmode = [6]
            kmagn = [1]
            nj = length(kmagn)
            for n in ens
                for j in 1:nj
     #           CALL ZLATMR( N, N, 'D', ISEED, 'N', WORK, 6, ONE, CONE,
     # $                      'T', 'N', WORK( N+1 ), 1, ONE,
     # $                      WORK( 2*N+1 ), 1, ONE, 'N', IDUMMA, 0, 0,
     # $                      ZERO, ANORM, 'NO', A, LDA, IWORK, IINFO )
                    A = zeros(T,n,n)
                    imode = kmode[j]
                    anorm = magns[kmagn[j]]
                    rcond = 1.0
                    latmr!(A,anorm,imode,rcond,kl=0,ku=0)
                    schurtest(A,tols[itype])
                end
            end
        end

        @testset "symmetric, random ev" begin
            itype=8
            kmagn = [1]
            kmode = [6]
            nj = length(kmagn)
            for n in ens
                for j in 1:nj
     # note sym='H'
     #           CALL ZLATMR( N, N, 'D', ISEED, 'H', WORK, 6, ONE, CONE,
     # $                      'T', 'N', WORK( N+1 ), 1, ONE,
     # $                      WORK( 2*N+1 ), 1, ONE, 'N', IDUMMA, N, N,
     # $                      ZERO, ANORM, 'NO', A, LDA, IWORK, IINFO )
                    A = zeros(T,n,n)
                    imode = kmode[j]
                    anorm = magns[kmagn[j]]
                    rcond = 1.0
                    latmr!(A,anorm,imode,rcond,sym='H',kl=n,ku=n)
                    schurtest(A,tols[itype])
                end
            end
        end

        @testset "general, random ev" begin
            itype=9
            kmagn = [1,2,3]
            # kmode = [4,3,1] # data statement is a lie
            kmode = [6,6,6]
            nj = length(kmagn)
            for n in ens
                for j in 1:nj
     #           CALL ZLATMR( N, N, 'D', ISEED, 'N', WORK, 6, ONE, CONE,
     # $                      'T', 'N', WORK( N+1 ), 1, ONE,
     # $                      WORK( 2*N+1 ), 1, ONE, 'N', IDUMMA, N, N,
     # $                      ZERO, ANORM, 'NO', A, LDA, IWORK, IINFO )
            # zdrves zeros first 2 rows/cols & last row w/o explanation
            # presumably to test balancing, which we don't do, so skip it.
                    A = zeros(T,n,n)
                    imode = kmode[j]
                    anorm = magns[kmagn[j]]
                    rcond = 1.0
                    (verbosity[] > 1) &&
                        println("latmr: n=$n $anorm, $imode, $rcond")
                    latmr!(A,anorm,imode,rcond)
                    schurtest(A,tols[itype])
                    # throw this in for coverage
                    schurtest(A,tols[itype],standard=!_standard[])
                end
            end
        end

        @testset "triangular, random ev" begin
            # trivial, not even run in zdrves
            itype=10
            kconds = [1]
            kmagn = [1]
            kmode = [1]
            nj = length(kmagn)
            for n in ens
                for j in 1:nj
     # note: kl=0
     #           CALL ZLATMR( N, N, 'D', ISEED, 'N', WORK, 6, ONE, CONE,
     # $                      'T', 'N', WORK( N+1 ), 1, ONE,
     # $                      WORK( 2*N+1 ), 1, ONE, 'N', IDUMMA, N, 0,
     # $                      ZERO, ANORM, 'NO', A, LDA, IWORK, IINFO )
#
                    A = zeros(T,n,n)
                    imode = 6
                    anorm = magns[kmagn[j]]
                    rcond = 1.0
                    latmr!(A,anorm,imode,rcond,kl=0)
                    schurtest(A,tols[itype])
                end
            end
        end

    end # group testset
end # type loop (BlasFloat)


text:
using EGR
using Redis

println("TestUniversal")

createOracleOutputLevel = 1
numEquivalentPasses = 0.9
algOutputLevel = 2
maxOutputNum=20

myREfunction(problem, opts, sd, expIndices,n) = false
myWriteFunction(problem, sd, opts, k, gnum, origWant, fromOutputsFunction) = false

for (gradientOracle, numVars, numTP, csDataType, LossFunctionString, myOutputter, L2reg, thisDataName, thisProblem) in Oracles

	println(" $thisDataName $LossFunctionString L2reg = $L2reg")

	maxG  = Int(round(numEquivalentPasses*numTP))

	expIndices=unique(round(Int, logspace(0,log10(maxG+1),maxOutputNum)))-1

	myOutputOpts =  OutputOpts(myOutputter,expIndices)

	myOptss(stepSizePower) = Opts(zeros(numVars,1); stepSizePower=stepSizePower, maxG=maxG, outputLevel=algOutputLevel, stepOutputLevel=0)

	chunkSize  = 2
	u = (k,I)-> k
	alg(thisProblem(Task(() -> getSequentialFinite(numTP, gradientOracle))),  myOptss(0), onlyAddBeta1(u,"uTest",numVars, chunkSize,false), myOutputOpts, myWriteFunction)
	s = (k,I)-> Int(floor(k==0 ? 0 : k))
	alg(thisProblem(Task(() -> getSequentialFinite(numTP, gradientOracle))),  myOptss(0), onlyUpdateBeta1(s, "sTest", numVars, int(floor(numTP/chunkSize)), chunkSize,false), myOutputOpts, myWriteFunction)
	alg(thisProblem(Task(() -> getSequentialFinite(numTP, gradientOracle))),  myOptss(0), uLinBeta1(1.0,numVars, chunkSize,false), myOutputOpts, myWriteFunction)
	alg(thisProblem(Task(() -> getSequentialFinite(numTP, gradientOracle))),  myOptss(0), uQuadBeta1(1.0,numVars, chunkSize,false), myOutputOpts, myWriteFunction)
	alg(thisProblem(Task(() -> getSequentialFinite(numTP, gradientOracle))),  myOptss(0), uExpBeta1(2.0,numVars, chunkSize,false), myOutputOpts, myWriteFunction)

end

println("TestUniversal successful!")
println()


text:
export Selector, Class, Object

# Selectors

selname(s::Ptr{Nothing}) =
  ccall(:sel_getName, Ptr{Cchar}, (Ptr{Nothing},),
        s) |> unsafe_string

struct Selector
  ptr::Ptr{Nothing}
  Selector(ptr::Ptr{Nothing}) = new(ptr)
end

unsafe_convert(::Type{Ptr{Nothing}}, sel::Selector) = sel.ptr

function Selector(name)
  Selector(ccall(:sel_registerName, Ptr{Nothing}, (Ptr{Cchar},),
                 pointer(string(name))))
end

macro sel_str(name)
  Selector(name)
end

name(sel::Selector) = selname(sel.ptr)

function show(io::IO, sel::Selector)
  print(io, "sel")
  show(io, string(name(sel)))
end

# Classes

struct Class
  ptr::Ptr{Nothing}
  Class(ptr::Ptr{Nothing}) = new(ptr)
end

unsafe_convert(::Type{Ptr{Nothing}}, class::Class) = class.ptr

classptr(name) = ccall(:objc_getClass, Ptr{Nothing}, (Ptr{Cchar},),
                       pointer(string(name)))

function Class(name)
  ptr = classptr(name)
  ptr == C_NULL && error("Couldn't find class $name")
  return Class(ptr)
end

classexists(name) = classptr(name) ≠ C_NULL

name(class::Class) =
  ccall(:class_getName, Ptr{Cchar}, (Ptr{Nothing},),
            class) |> unsafe_string |> Symbol

ismeta(class::Class) =
  ccall(:class_isMetaClass, Cint, (Ptr{Nothing},),
        class) |> int2bool

function super(class::Class)
  ptr = ccall(:class_getSuperclass, Ptr{Nothing}, (Ptr{Nothing},),
              class.ptr)
  ptr == C_NULL && return nothing
  Class(ptr)
end

function show(io::IO, class::Class)
  ismeta(class) && print(io, "^")
  print(io, name(class))
end

function methods(class::Class)
  count = Cuint[0]
  meths = ccall(:class_copyMethodList, Ptr{Ptr{Nothing}}, (Ptr{Nothing}, Ptr{Cuint}),
                class, count)
  meths′ = [unsafe_load(meths, i) for i = 1:count[1]]
  c_free(meths)
  meths = [ccall(:method_getName, Ptr{Nothing}, (Ptr{Nothing},), meth) for meth in meths′]
  return map(meth->selname(meth), meths)
end

# Objects

mutable struct Object
  ptr::Ptr{Nothing}
end

unsafe_convert(::Type{Ptr{Nothing}}, obj::Object) = obj.ptr

class(obj) =
  ccall(:object_getClass, Ptr{Nothing}, (Ptr{Nothing},),
        obj) |> Class

methods(obj::Object) = methods(class(obj))

show(io::IO, obj::Object) = print(io, class(obj), " Object")


text:
## 3D cases
#==========================================================================================#

"""
    print_info(wing :: Union{Wing, HalfWing})

Print the relevant geometric characteristics of a `HalfWing` or `Wing`.
"""
function print_info(wing :: Union{Wing, HalfWing}, head = ""; browser = false)
    labels = [ "Span (m)", "Area (m²)", "MAC (m)", "Aspect Ratio"]
    wing_info = [ info(wing)... ]
    data = [ labels wing_info ]
    header = [ head, "Value" ]
    h1 = Highlighter( (data,i,j) -> (j == 1), foreground = :blue, bold = true)
    if browser
        pretty_table(String, data, header, alignment = [:c, :c], tf = tf_html_minimalist, backend = :html, highlighters = HTMLHighlighter( (data,i,j) -> (j == 1), HTMLDecoration(font_weight = "bold")), formatters = ft_round(8))
    else
        pretty_table(data, header, alignment = [:c, :c], tf = tf_compact, highlighters = h1, vlines = :none, formatters = ft_round(8))
    end
end

VortexLattice.VLMState(fs :: Freestream; rho_ref = 1.225, r_ref = zeros(3), area_ref = 1, chord_ref = 1, span_ref = 1, name = "Aircraft") = VLMState(fs.V, fs.alpha, fs.beta, fs.omega, r_ref, rho_ref, area_ref, chord_ref, span_ref, name)

"""
    solve_case(horseshoe_panels :: Matrix{Panel3D}, normals, freestream :: Freestream, r_ref, ρ = 1.225; symmetry = false)

Evaluate a vortex lattice case given an array of `Panel3D`s with associated normal vectors (not necessarily the same as the panels' normals), a `Freestream`, reference density ``\\rho`` and reference point ``r_\\text{ref}`` for moments, with an option for symmetry.
"""
function solve_case(horseshoe_panels :: Array{<: Panel3D}, normals, U, α, β, Ω, rho_ref, area_ref, chord_ref, span_ref, r_ref)
    # Make horseshoes and collocation points
    horseshoes = horseshoe_line.(horseshoe_panels)

    # Solve system
    Γs = reshape(solve_system(horseshoes[:], normals[:], U, Ω), size(horseshoe_panels))

    # Compute forces and moments
    surface_forces, surface_moments, trefftz_force = case_dynamics(Γs, horseshoes, U, α, β, Ω, rho_ref, r_ref)

    # Compute aerodynamic coefficients
    nearfield_coeffs, farfield_coeffs, CFs, CMs = evaluate_coefficients(surface_forces, surface_moments, trefftz_force, U, α, β, Ω, rho_ref, area_ref, chord_ref, span_ref)

    nearfield_coeffs, farfield_coeffs, CFs, CMs, horseshoes, Γs
end

"""
    solve_case(wing :: Union{Wing, HalfWing}, freestream :: Freestream, ρ :: Real, r_ref; span_num :: Integer = 5, chord_num :: Integer = 10)

Evaluate a vortex lattice case given a `Wing` or `HalfWing` with a `Freestream`, reference density ``\\rho`` and reference point ``r_\\text{ref}`` for moments, ``n_s`` span-wise panels and ``n_c`` chord-wise panels.
"""
function solve_case(wing :: Union{Wing, HalfWing}, freestream :: Freestream; rho_ref = 1.225, r_ref = [0.25, 0, 0], area_ref = projected_area(wing), chord_ref = mean_aerodynamic_chord(wing), span_ref = span(wing), mu_ref = 1.5e-5, span_num :: Union{Integer, Vector{<: Integer}}, chord_num :: Integer, viscous = false, a_ref = 330., x_tr = 0.3, spacing = spanwise_spacing(wing))
    # Unpack Freestream
    U, α, β, Ω = aircraft_velocity(freestream), freestream.alpha, freestream.beta, freestream.omega

    # Determine spanwise panel distribution and spacing
    space     = ifelse(typeof(spacing) <: String, [spacing], spacing)
    span_nums = number_of_spanwise_panels(wing, span_num)

    # Compute panels and normals
    horseshoe_panels, camber_panels = vlmesh_wing(wing, span_nums, chord_num, space)
    normals = panel_normal.(camber_panels)

    # Compute forces and moments
    nearfield_coeffs, farfield_coeffs, CFs, CMs, horseshoes, Γs = solve_case(horseshoe_panels, normals, U, α, β, Ω, rho_ref, area_ref, chord_ref, span_ref, r_ref)

    # Viscous drag evaluation
    if viscous
        # Compute profile drag via wetted-area equivalent skin-friction method
        CDv = profile_drag_coefficient(wing, x_tr, freestream.V, rho_ref, a_ref, area_ref, mu_ref)

        # Add profile and total drag coefficients
        nf_coeffs = [ nearfield_coeffs[1] + CDv; CDv; nearfield_coeffs ]
        ff_coeffs = [  farfield_coeffs[1] + CDv; CDv; farfield_coeffs  ]
    else
        nf_coeffs = nearfield_coeffs
        ff_coeffs = farfield_coeffs
    end

    nf_coeffs, ff_coeffs, CFs, CMs, horseshoe_panels, normals, horseshoes, Γs
end

function solve_case(components :: Dict{String, Tuple{Matrix{Panel3D{T}}, Matrix{SVector{3,T}}}}, freestream :: Freestream; rho_ref = 1.225, r_ref = zeros(3), area_ref = 1, chord_ref = 1, span_ref = 1, name = "Aircraft", print = false, print_components = false) where T <: Real
    # Get panels
    meshes = values(components)

    # Flattening for VLM
    horseshoe_panels = first.(meshes)
    normals          = last.(meshes)
    horsies          = reduce(vcat, vec.(horseshoe_panels))
    normies          = reduce(vcat, vec.(normals))

    # Get required vortex lattice variables, i.e. horseshoes, collocation points and normals
    horseshoes = horseshoe_line.(horsies)
    horseshoes_arr = [ horseshoe_line.(horses) for horses in horseshoe_panels ]

    # Unpack Freestream
    U, α, β, Ω = aircraft_velocity(freestream), freestream.alpha, freestream.beta, freestream.omega

    # Solve system
    Γs = solve_system(horseshoes, normies, U, Ω)

    # Reshaping
    panel_sizes = size.(horseshoe_panels)
    panel_inds 	= [ 0; cumsum(prod.(panel_sizes)) ]
    Γs_arr 		= reshape_array(Γs, panel_inds, panel_sizes)

    # Compute forces and moments
    results = case_dynamics.(Γs_arr, horseshoes_arr, Ref(Γs), Ref(horseshoes), Ref(U), α, β, Ref(Ω), rho_ref, Ref(r_ref))

    # Components' non-dimensional forces and moments
    data = [ evaluate_coefficients(dyn..., U, α, β, Ω, rho_ref, area_ref, chord_ref, span_ref) for dyn in results ]

    nf_comp_coeffs = getindex.(data, 1)
    ff_comp_coeffs = getindex.(data, 2)
    CFs            = getindex.(data, 3)
    CMs            = getindex.(data, 4)

    # Aircraft's non-dimensional forces and moments
    nf_coeffs = reduce((x, y) -> x .+ y, nf_comp_coeffs) # Sum nearfield coefficients
    ff_coeffs = reduce((x, y) -> x .+ y, ff_comp_coeffs) # Sum farfield coefficients
    name_CFs  = reduce(vcat, vec.(CFs))                  # Collect surface force coefficients
    name_CMs  = reduce(vcat, vec.(CMs))                  # Collect surface moment coefficients

    # Printing
    if print_components; print_coefficients.(nf_comp_coeffs, ff_comp_coeffs, keys(components)) end
    if print; 			 print_coefficients(nf_coeffs, ff_coeffs, name) 					   end

    # Dictionary assembly
    name_data = (nf_coeffs, ff_coeffs, name_CFs, name_CMs, horsies, normies, horseshoes, Γs)
    comp_data = tuple.(nf_comp_coeffs, ff_comp_coeffs, CFs, CMs, horseshoe_panels, normals, horseshoes_arr, Γs_arr)

    names 	= [ name						 ; # Aircraft name
                (collect ∘ keys)(components) ] # Component names
    data  	= [ name_data ;	# Aircraft data
                comp_data ]	# Component data

    Dict(names .=> data)
end

streamlines(freestream :: Freestream, points, horseshoes, Γs, length, num_steps :: Integer) = VortexLattice.streamlines.(points, Ref(velocity(freestream)), Ref(freestream.omega), Ref(horseshoes), Ref(Γs), Ref(length), Ref(num_steps))

max_forks_repo_name: GodotMisogi/AeroMDAO


text:
#
# This file is part of the DiscreteEvents.jl Julia package, MIT license
#
# Paul Bayer, 2020
#
# This is a Julia package for discrete event simulation
#
# -------------------------------------------
#
# This is the first introductory pet example.
# We simulate the easy life of a pet in the morning with an event-based simulation.
#

using DiscreteEvents, Printf, Random

# First we define some data structures for pets:
abstract type PetState end
struct Sleeping  <: PetState end
struct Scuffing  <: PetState end
struct Running   <: PetState end

abstract type PetEvent end
struct GetWeary  <: PetEvent end
struct GetHungry <: PetEvent end
struct Scuff     <: PetEvent end
struct LeapUp    <: PetEvent end
struct Sleep     <: PetEvent end

mutable struct Pet
    clk::Clock
    name::String
    state::PetState
    speak::String
end

# Then we need an API for Pet since we do not want to access its state directly.
state(p) = p.state
speak(p, n) = @printf("%5.2f %s: %s\n", tau(p.clk), p.name, p.speak^n)

function setstate!(p::Pet, q::PetState)
    p.state = q
    @printf("%5.2f %s: %s\n", tau(p.clk), p.name, repr(p.state))
end

# We describe the behaviour of the pet in time with some state-transition functions:
function doit!(p::Pet, ::Sleeping, ::LeapUp)   # leap up after sleeping
    setstate!(p, Running())
    event!(p.clk, fun(doit!, p, fun(state, p), GetHungry()), after, 5*rand())
end

function doit!(p::Pet, ::Scuffing, ::LeapUp)   # leap up while scuffing
    setstate!(p, Running())
    event!(p.clk, fun(doit!, p, fun(state, p), GetWeary()), after, 2*rand())
end

function doit!(p::Pet, ::Running, ::GetHungry) # get hungry while running
    speak(p, 5)
    event!(p.clk, fun(doit!, p, fun(state, p), Scuff()), after, rand())
end

function doit!(p::Pet, ::Running, ::GetWeary)  # get weary while running
    speak(p, 2)
    event!(p.clk, fun(doit!, p, fun(state, p), Sleep()), after, 2*rand())
end

function doit!(p::Pet, ::Running, ::Scuff)     # scuff after running
    setstate!(p, Scuffing())
    event!(p.clk, fun(doit!, p, fun(state, p), LeapUp()), after, 2*rand())
end

function doit!(p::Pet, ::Running, ::Sleep)     # sleep after running
    setstate!(p, Sleeping())
    event!(p.clk, fun(doit!, p, fun(state, p), LeapUp()), after, 10*rand())
end

# Then we setup a clock and a pet, schedule the first event on it and run the clock:
Random.seed!(123)     # we seed the random number generator for reprodicibility
clk = Clock()
snoopy = Pet(clk, "Snoopy", Sleeping(), "huff")
event!(clk, fun(doit!, snoopy, snoopy.state, LeapUp()), after, 5)
run!(clk, 25)


text:
module WaveToy

using DifferentialEquations
using Fun1d
# using Random

struct Wave{S,T} <: AbstractArray{T,2}
    ϕ::GridFun{S,T}
    ψ::GridFun{S,T}
end

# Iteration
Base.IteratorSize(::Type{<:Wave}) = Iterators.HasShape{2}()
Base.eltype(::Type{Wave{S,T}}) where {S,T} = T
Base.isempty(x::Wave) = isempty(x.ϕ)
function Base.iterate(x::Wave, state...)
    return iterate(Iterators.flatten((x.ϕ, x.ψ)), state...)
end
# Base.length(x::Wave) = 2 * length(x.ϕ)
Base.size(x::Wave) = (length(x.ϕ), 2)
Base.size(x::Wave, d) = size(x)[d]

# Indexing
function lin2cart(x::Wave, i::Number)
    n = length(x.ϕ)
    return (i - 1) % n + 1, (i - 1) ÷ n + 1
end
Base.firstindex(x::Wave) = error("not implemented")
Base.getindex(x::Wave, i) = getindex(x, i.I...)
Base.getindex(x::Wave, i::Number) = getindex(x, lin2cart(x, i)...)
Base.getindex(x::Wave, i, j) = getindex((x.ϕ, x.ψ)[j], i)
Base.lastindex(x::Wave) = error("not implemented")
Base.setindex!(x::Wave, v, i) = setindex!(x, v, i.I...)
Base.setindex!(x::Wave, v, i::Number) = setindex!(x, v, lin2cart(x, i))
Base.setindex!(x::Wave, v, i, j) = setindex!((x.ϕ, x.ψ)[j], v, i)

# Abstract Array
Base.IndexStyle(::Wave) = IndexCartesian()
Base.similar(x::Wave) = Wave(similar(x.ϕ), similar(x.ψ))
function Base.similar(x::Wave, ::Type{T}) where {T}
    return Wave(similar(x.ϕ, T), similar(x.ψ, T))
end
Base.similar(x::Wave, ::Dims) = similar(x)
Base.similar(x::Wave, ::Dims, ::Type{T}) where {T} = similar(x, T)

# Broadcasting
Base.BroadcastStyle(::Type{<:Wave}) = Broadcast.ArrayStyle{Wave}()
function Base.similar(bc::Broadcast.Broadcasted{Broadcast.ArrayStyle{Wave}},
                      ::Type{T}) where {T}
    x = find_wave(bc)
    return similar(x, T)
end
find_wave(bc::Base.Broadcast.Broadcasted) = find_wave(bc.args)
find_wave(args::Tuple) = find_wave(find_wave(args[1]), Base.tail(args))
find_wave(x) = x
find_wave(::Tuple{}) = nothing
find_wave(a::Wave, rest) = a
find_wave(::Any, rest) = find_wave(rest)

# Others
function Base.map(fun, x::Wave, ys::Wave...)
    return Wave(map(fun, x.ϕ, (y.ϕ for y in ys)...),
                map(fun, x.ψ, (y.ψ for y in ys)...))
end

# function Base.rand(rng::AbstractRNG, ::Random.SamplerType{Wave{T}}) where {T}
#     return Wave{T}(rand(rng, T), rand(rng, T))
# end

Base.zero(::Type{<:Wave}) = error("not implemented")
Base.zero(x::Wave) = Wave(zero(x.ϕ), zero(x.ψ))

Base.:+(x::Wave) = map(+, x)
Base.:-(x::Wave) = map(-, x)

Base.:+(x::Wave, y::Wave) = map(+, x, y)
Base.:-(x::Wave, y::Wave) = map(-, x, y)

Base.:*(x::Wave, a::Number) = map(b -> b * a, x)
Base.:*(a::Number, x::Wave) = map(b -> a * b, x)
Base.:/(x::Wave, a::Number) = map(b -> b / a, x)
Base.:\(a::Number, x::Wave) = map(b -> a \ b, x)

################################################################################

function setup(::Type{S}) where {S}
    domain = Domain{S}(0, 1)
    grid = Grid(domain, 20)
    return grid
end

function init(::Type{T}, grid::Grid) where {T}
    # println("WaveToy.init")
    # ∂ₜₜu = ∂ₓₓu
    # u(t,x) = f(t+x) + g(t-x)
    # u(t,x) = sin(π t) cos(π x)
    ϕ = project(T, grid, x -> 0)
    ψ = project(T, grid, x -> π * cospi(x))
    return Wave(ϕ, ψ)
end

function rhs(state::Wave, param, t)
    # println("WaveToy.rhs t=$t")
    ϕdot = state.ψ
    ψdot = deriv2(state.ϕ)
    staterhs = Wave(ϕdot, ψdot)
    return staterhs::Wave
end

function main()
    T = Float64
    grid = setup(T)
    state = init(T, grid)::Wave
    tspan = T[0, 1]
    atol = eps(T)^(T(3) / 4)
    prob = ODEProblem(rhs, state, tspan)
    sol = solve(prob, Tsit5(); abstol=atol)
    return sol
end

# include("examples/wavetoy.jl")
# sol = WaveToy.main();
# using Plots
# plot([sol(t).ϕ for t in LinRange(0, 1, 11)])

################################################################################

function setup_sphere(::Type{S},points) where {S}
    domain = Domain{S}(0, 5)
    grid = Grid(domain, points)
    return grid
end

function init_sphere(::Type{T}, grid::Grid, t::T) where {T}
    # println("WaveToy.init_sphere")
    # ∂ₜₜu = ∂ᵣᵣu + 2/r ∂ᵣu
    # u(t,r) = 1/r f(t+r) + 1/r g(t-r)
    # u(t,r) = 1/r (exp(-(t-r)²) - exp(-(t+r)²))
    rmin = sqrt(eps(T))
    function u(t, r)
        if abs(r) <= rmin
            return (4 + 4 / T(3) * (-3 + 2 * t^2) * r^2) * t * exp(-t^2)
        else
            return 1 / r * (exp(-(t - r)^2) - exp(-(t + r)^2))
        end
    end
    function u̇(t, r)
        if abs(r) <= rmin
            return (4 * (1 - 2 * t^2) +
                    4 / T(3) * (-3 + 12 * t^2 - 4 * t^4) * r^2) * exp(-t^2)
        else
            return -2 / r *
                   ((t - r) * exp(-(t - r)^2) - (t + r) * exp(-(t + r)^2))
        end
    end
    ϕ = project(T, grid, r -> u(t, r))
    ψ = project(T, grid, r -> u̇(t, r))
    return Wave(ϕ, ψ)
end

function rhs_sphere(state::Wave, param, t)
    # println("WaveToy.rhs_sphere t=$t")
    ϕdot = state.ψ
    ψdot = deriv2_sphere(state.ϕ)
    staterhs = Wave(ϕdot, ψdot)
    return staterhs::Wave
end

function main_sphere(points)
    T = Float64
    grid = setup_sphere(T,points)
    t0 = T(0)
    state = init_sphere(T, grid, t0)::Wave
    t1 = T(3)
    prob = ODEProblem(rhs_sphere, state, (t0, t1))
    atol = eps(T)^(T(3) / 4)
    sol = solve(prob, Tsit5(); abstol=atol)
    return sol
end

# include("examples/wavetoy.jl")
# sol = WaveToy.main_sphere();
# using Plots
# plot([sol(t).ϕ for t in LinRange(0, 3, 11)])

end