text:
#Assume all local Hamiltonians have the same nrange
function Base.:*(H::Hᶜ{MPO}, v::ITensor)
  ∑h = H.∑h
  Hᴸ = H.Hᴸ
  Hᴿ = H.Hᴿ
  ψ = H.ψ
  ψ′ = dag(ψ)'
  Nsites = nsites(ψ)
  range_∑h = nrange(∑h, 1)
  n = H.n
  l = linkinds(only, ψ.AL)
  l′ = linkinds(only, ψ′.AL)
  r = linkinds(only, ψ.AR)
  r′ = linkinds(only, ψ′.AR)
  s = siteinds(only, ψ)
  s′ = siteinds(only, ψ′)
  δˢ(n) = δ(dag(s[n]), prime(s[n]))
  δʳ(n) = δ(dag(r[n]), prime(r[n]))
  δˡ(n) = δ(l[n], l′[n])

  #Build the contribution of the left environment
  Hᶜᴸv = v * Hᴸ[n] * δʳ(n)
  #Build the contribution of the right environment
  Hᶜᴿv = v * δˡ(n) * Hᴿ[n]
  #We now start building terms where C overlap with the local Hamiltonian
  # We start with the tensor AL[n] - v - AR[n+1] ... AR[n + range_∑h - 1]
  Hᶜʰv = v * ψ.AL[n] * δˡ(n - 1) * ∑h[n][1] * ψ′.AL[n] #left extremity
  common_sites = findsites(ψ, ∑h[n])
  idx = 2 #list the sites Σh, we start at 2 because n is already taken into account
  for k in 1:(range_∑h - 2)
    if n + k == common_sites[idx]
      Hᶜʰv = Hᶜʰv * ψ.AR[n + k] * ∑h[n][idx] * ψ′.AR[n + k]
      idx += 1
    else
      Hᶜʰv = Hᶜʰv * ψ.AR[n + k] * δˢ(n + k) * ψ′.AR[n + k]
    end
  end
  Hᶜʰv =
    Hᶜʰv *
    ψ.AR[n + range_∑h - 1] *
    δʳ(n + range_∑h - 1) *
    ∑h[n][end] *
    ψ′.AR[n + range_∑h - 1] #right most extremity
  #Now we are building contributions of the form AL[n - j] ... AL[n] - v - AR[n + 1] ... AR[n + range_∑h - 1 - j]
  for j in 1:(range_∑h - 2)
    temp_Hᶜʰv = ψ.AL[n - j] * δˡ(n - 1 - j) * ∑h[n - j][1] * ψ′.AL[n - j]
    common_sites = findsites(ψ, ∑h[n - j])
    idx = 2
    for k in 1:j
      if n - j + k == common_sites[idx]
        temp_Hᶜʰv = temp_Hᶜʰv * ψ.AL[n - j + k] * ∑h[n - j][idx] * ψ′.AL[n - j + k]
        idx += 1
      else
        temp_Hᶜʰv = temp_Hᶜʰv * ψ.AL[n - j + k] * δˢ(n - j + k) * ψ′.AL[n - j + k]
      end
    end
    # Finished the AL part
    temp_Hᶜʰv = temp_Hᶜʰv * v
    for k in (j + 1):(range_∑h - 2)
      if n - j + k == common_sites[idx]
        temp_Hᶜʰv = temp_Hᶜʰv * ψ.AR[n - j + k] * ∑h[n - j][idx] * ψ′.AR[n - j + k]
        idx += 1
      else
        temp_Hᶜʰv = temp_Hᶜʰv * ψ.AR[n - j + k] * δˢ(n - j + k) * ψ′.AR[n - j + k]
      end
    end
    temp_Hᶜʰv =
      temp_Hᶜʰv *
      ψ.AR[n - j + range_∑h - 1] *
      δʳ(n - j + range_∑h - 1) *
      ∑h[n - j][end] *
      ψ′.AR[n - j + range_∑h - 1]
    Hᶜʰv = Hᶜʰv + temp_Hᶜʰv
  end
  Hᶜv = Hᶜᴸv + Hᶜʰv + Hᶜᴿv
  return Hᶜv * dag(δˡ(n)) * dag(δʳ(n))
end

function Base.:*(H::Hᴬᶜ{MPO}, v::ITensor)
  ∑h = H.∑h
  Hᴸ = H.Hᴸ
  Hᴿ = H.Hᴿ
  ψ = H.ψ
  ψ′ = dag(ψ)'
  Nsites = nsites(ψ)
  range_∑h = nrange(∑h, 1)
  n = H.n
  l = linkinds(only, ψ.AL)
  l′ = linkinds(only, ψ′.AL)
  r = linkinds(only, ψ.AR)
  r′ = linkinds(only, ψ′.AR)
  s = siteinds(only, ψ)
  s′ = siteinds(only, ψ′)

  δˢ(n) = δ(dag(s[n]), prime(s[n]))
  δʳ(n) = δ(dag(r[n]), prime(r[n]))
  δˡ(n) = δ(l[n], l′[n])

  #Contribution of the left environment
  Hᴬᶜᴸv = v * Hᴸ[n - 1] * δˢ(n) * δʳ(n)
  #Contribution of the right environment
  Hᴬᶜᴿv = v * δˡ(n - 1) * δˢ(n) * Hᴿ[n]
  #We now start building terms where AC overlap with the local Hamiltonian
  # We start with the tensor v - AR[n+1] ... AR[n + range_∑h - 1]
  Hᴬᶜʰv = v * δˡ(n - 1) * ∑h[n][1]
  common_sites = findsites(ψ, ∑h[n])
  idx = 2 #list the sites Σh, we start at 2 because n is already taken into account
  for k in 1:(range_∑h - 2)
    if n + k == common_sites[idx]
      Hᴬᶜʰv = Hᴬᶜʰv * ψ.AR[n + k] * ∑h[n][idx] * ψ′.AR[n + k]
      idx += 1
    else
      Hᴬᶜʰv = Hᴬᶜʰv * ψ.AR[n + k] * δˢ(n + k) * ψ′.AR[n + k]
    end
  end
  Hᴬᶜʰv =
    Hᴬᶜʰv *
    ψ.AR[n + range_∑h - 1] *
    δʳ(n + range_∑h - 1) *
    ∑h[n][end] *
    ψ′.AR[n + range_∑h - 1] #rightmost extremity
  #Now we are building contributions of the form AL[n - j] ... AL[n-1] - v - AR[n + 1] ... AR[n + range_∑h - 1 - j]
  for j in 1:(range_∑h - 1)
    temp_Hᴬᶜʰv = ψ.AL[n - j] * δˡ(n - j - 1) * ∑h[n - j][1] * ψ′.AL[n - j]
    common_sites = findsites(ψ, ∑h[n - j])
    idx = 2
    for k in 1:(j - 1)
      if n - j + k == common_sites[idx]
        temp_Hᴬᶜʰv = temp_Hᴬᶜʰv * ψ.AL[n - j + k] * ∑h[n - j][idx] * ψ′.AL[n - j + k]
        idx += 1
      else
        temp_Hᴬᶜʰv = temp_Hᴬᶜʰv * ψ.AL[n - j + k] * δˢ(n - j + k) * ψ′.AL[n - j + k]
      end
    end
    #Finished with AL, treating the center AC = v
    if j == range_∑h - 1
      temp_Hᴬᶜʰv = temp_Hᴬᶜʰv * v * ∑h[n - j][end] * δʳ(n - j + range_∑h - 1)
    else
      if n == common_sites[idx] #need to check whether we need to branch v
        temp_Hᴬᶜʰv = temp_Hᴬᶜʰv * v * ∑h[n - j][idx]
        idx += 1
      else
        temp_Hᴬᶜʰv = temp_Hᴬᶜʰv * v * δˢ(n)
      end
      for k in (j + 1):(range_∑h - 2)
        if n + k - j == common_sites[idx]
          temp_Hᴬᶜʰv = temp_Hᴬᶜʰv * ψ.AR[n + k - j] * ∑h[n - j][idx] * ψ′.AR[n + k - j]
          idx += 1
        else
          temp_Hᴬᶜʰv = temp_Hᴬᶜʰv * ψ.AR[n + k - j] * δˢ(n + k - j) * ψ′.AR[n + k - j]
        end
      end
      temp_Hᴬᶜʰv =
        temp_Hᴬᶜʰv *
        ψ.AR[n + range_∑h - 1 - j] *
        δʳ(n - j + range_∑h - 1) *
        ∑h[n - j][end] *
        ψ′.AR[n + range_∑h - 1 - j]
    end
    Hᴬᶜʰv = Hᴬᶜʰv + temp_Hᴬᶜʰv
  end
  Hᴬᶜv = Hᴬᶜᴸv + Hᴬᶜʰv + Hᴬᶜᴿv
  return Hᴬᶜv * dag(δˡ(n - 1)) * dag(δˢ(n)) * dag(δʳ(n))
end

function left_environment(∑h::InfiniteSum{MPO}, ψ::InfiniteCanonicalMPS; tol=1e-15)
  Nsites = nsites(ψ)
  range_∑h = nrange(∑h, 1)
  ψᴴ = dag(ψ)
  ψ′ = ψᴴ'
  ψ̃ = prime(linkinds, ψᴴ)

  l = linkinds(only, ψ.AL)
  l′ = linkinds(only, ψ′.AL)
  r = linkinds(only, ψ.AR)
  r′ = linkinds(only, ψ′.AR)
  s = siteinds(only, ψ)
  δʳ(n) = δ(dag(r[n]), prime(r[n]))
  δˡ(n) = δ(l[n], l′[n])
  δˢ(n) = δ(dag(s[n]), prime(s[n]))
  hᴸ = Vector{ITensor}(undef, Nsites)
  for k in 1:Nsites
    hᴸ[k] =
      δˡ(k - range_∑h) *
      ψ.AL[k - range_∑h + 1] *
      ∑h[(k - range_∑h + 1, k - range_∑h + 2)][1] *
      ψ′.AL[k - range_∑h + 1]
    common_sites = findsites(ψ, ∑h[(k - range_∑h + 1, k - range_∑h + 2)])
    idx = 2
    for j in 2:range_∑h
      if k - range_∑h + j == common_sites[idx]
        hᴸ[k] =
          hᴸ[k] *
          ψ.AL[k - range_∑h + j] *
          ∑h[(k - range_∑h + 1, k - range_∑h + 2)][idx] *
          ψ′.AL[k - range_∑h + j]
        idx += 1
      else
        hᴸ[k] =
          hᴸ[k] * ψ.AL[k - range_∑h + j] * δˢ(k - range_∑h + j) * ψ′.AL[k - range_∑h + j]
      end
    end
  end
  hᴸ = InfiniteMPS(hᴸ)
  eᴸ = [(hᴸ[k] * ψ.C[k] * δʳ(k) * ψ′.C[k])[] for k in 1:Nsites]
  for k in 1:Nsites
    # TODO: remove `denseblocks` once BlockSparse + DiagBlockSparse is supported
    hᴸ[k] -= eᴸ[k] * denseblocks(δ(inds(hᴸ[k])))
  end

  𝕙ᴸ = copy(hᴸ)
  # TODO restrict to the useful ones only?
  for n in 1:Nsites
    for k in 1:(Nsites - 1)
      temp = copy(hᴸ[n - k])
      for kp in reverse(0:(k - 1))
        temp = temp * ψ.AL[n - kp] * ψ̃.AL[n - kp]
      end
      𝕙ᴸ[n] = temp + 𝕙ᴸ[n]
    end
  end
  Hᴸ = left_environment(hᴸ, 𝕙ᴸ, ψ; tol=tol)

  return Hᴸ, eᴸ
end

function right_environment(∑h::InfiniteSum{MPO}, ψ::InfiniteCanonicalMPS; tol=1e-15)
  Nsites = nsites(ψ)
  range_∑h = nrange(∑h, 1)
  ψᴴ = dag(ψ)
  ψ′ = ψᴴ'
  ψ̃ = prime(linkinds, ψᴴ)

  l = linkinds(only, ψ.AL)
  l′ = linkinds(only, ψ′.AL)
  r = linkinds(only, ψ.AR)
  r′ = linkinds(only, ψ′.AR)
  s = siteinds(only, ψ)
  δʳ(n) = δ(dag(r[n]), prime(r[n]))
  δˡ(n) = δ(l[n], l′[n])
  δˢ(n) = δ(dag(s[n]), prime(s[n]))

  hᴿ = Vector{ITensor}(undef, Nsites)
  for k in 1:Nsites
    hᴿ[k] = ψ.AR[k + range_∑h] * δʳ(k + range_∑h) * ∑h[k + 1][end] * ψ′.AR[k + range_∑h]
    common_sites = findsites(ψ, ∑h[k + 1])
    idx = length(common_sites) - 1
    for j in (range_∑h - 1):-1:1
      if k + j == common_sites[idx]
        hᴿ[k] = hᴿ[k] * ψ.AR[k + j] * ∑h[k + 1][idx] * ψ′.AR[k + j]
        idx -= 1
      else
        hᴿ[k] = hᴿ[k] * ψ.AR[k + j] * δˢ(k + j) * ψ′.AR[k + j]
      end
    end
  end
  hᴿ = InfiniteMPS(hᴿ)
  eᴿ = [(hᴿ[k] * ψ.C[k] * δˡ(k) * ψ′.C[k])[] for k in 1:Nsites]
  for k in 1:Nsites
    hᴿ[k] -= eᴿ[k] * denseblocks(δ(inds(hᴿ[k])))
  end

  𝕙ᴿ = copy(hᴿ)
  #TODO restrict to the useful ones only
  for n in 1:Nsites
    for k in 1:(Nsites - 1)
      temp = copy(hᴿ[n + k])
      for kp in reverse(1:k)
        temp = temp * ψ.AR[n + kp] * ψ̃.AR[n + kp]
      end
      𝕙ᴿ[n] = temp + 𝕙ᴿ[n]
    end
  end
  Hᴿ = right_environment(hᴿ, 𝕙ᴿ, ψ; tol=tol)
  return Hᴿ, eᴿ
end

#In principle, could share even more code with vumps_mpo or with parallel
function tdvp_iteration_sequential(
  solver::Function,
  ∑h::InfiniteSum{MPO},
  ψ::InfiniteCanonicalMPS;
  (ϵᴸ!)=fill(1e-15, nsites(ψ)),
  (ϵᴿ!)=fill(1e-15, nsites(ψ)),
  time_step,
  solver_tol=(x -> x / 100),
)
  Nsites = nsites(ψ)
  range_∑h = nrange(∑h, 1)
  ϵᵖʳᵉˢ = max(maximum(ϵᴸ!), maximum(ϵᴿ!))
  _solver_tol = solver_tol(ϵᵖʳᵉˢ)

  ψ = copy(ψ)
  C̃ = InfiniteMPS(Vector{ITensor}(undef, Nsites))
  Ãᶜ = InfiniteMPS(Vector{ITensor}(undef, Nsites))
  Ãᴸ = InfiniteMPS(Vector{ITensor}(undef, Nsites))
  Ãᴿ = InfiniteMPS(Vector{ITensor}(undef, Nsites))
  eᴸ = Vector{Float64}(undef, Nsites)
  eᴿ = Vector{Float64}(undef, Nsites)

  for n in 1:Nsites
    ψᴴ = dag(ψ)
    ψ′ = ψᴴ'
    # XXX: make this prime the center sites
    ψ̃ = prime(linkinds, ψᴴ)

    l = linkinds(only, ψ.AL)
    l′ = linkinds(only, ψ′.AL)
    r = linkinds(only, ψ.AR)
    r′ = linkinds(only, ψ′.AR)
    s = siteinds(only, ψ)
    δʳ(n) = δ(dag(r[n]), prime(r[n]))
    δˡ(n) = δ(l[n], l′[n])
    δˢ(n) = δ(dag(s[n]), prime(s[n]))

    Hᴸ, eᴸ = left_environment(∑h, ψ; tol=_solver_tol)
    Hᴿ, eᴿ = right_environment(∑h, ψ; tol=_solver_tol)

    Cvalsₙ₋₁, Cvecsₙ₋₁, Cinfoₙ₋₁ = solver(
      Hᶜ(∑h, Hᴸ, Hᴿ, ψ, n - 1), time_step, ψ.C[n - 1], _solver_tol
    )
    Cvalsₙ, Cvecsₙ, Cinfoₙ = solver(Hᶜ(∑h, Hᴸ, Hᴿ, ψ, n), time_step, ψ.C[n], _solver_tol)
    Avalsₙ, Avecsₙ, Ainfoₙ = solver(
      Hᴬᶜ(∑h, Hᴸ, Hᴿ, ψ, n), time_step, ψ.AL[n] * ψ.C[n], _solver_tol
    )

    C̃[n - 1] = Cvecsₙ₋₁
    C̃[n] = Cvecsₙ
    Ãᶜ[n] = Avecsₙ

    Ãᴸ[n] = ortho_polar(Ãᶜ[n], C̃[n])
    Ãᴿ[n] = ortho_polar(Ãᶜ[n], C̃[n - 1])
    # Update state for next iteration
    #ψ = InfiniteCanonicalMPS(Ãᴸ, C̃, Ãᴿ)
    ψ.AL[n] = Ãᴸ[n]
    ψ.AR[n] = Ãᴿ[n]
    ψ.C[n - 1] = C̃[n - 1]
    ψ.C[n] = C̃[n]
  end

  for n in 1:Nsites
    ϵᴸ![n] = norm(Ãᶜ[n] - Ãᴸ[n] * C̃[n])
    ϵᴿ![n] = norm(Ãᶜ[n] - C̃[n - 1] * Ãᴿ[n])
  end
  return ψ, (eᴸ, eᴿ)
end

function tdvp_iteration_parallel(
  solver::Function,
  ∑h::InfiniteSum{MPO},
  ψ::InfiniteCanonicalMPS;
  (ϵᴸ!)=fill(1e-15, nsites(ψ)),
  (ϵᴿ!)=fill(1e-15, nsites(ψ)),
  time_step,
  solver_tol=(x -> x / 100),
)
  Nsites = nsites(ψ)
  range_∑h = nrange(∑h, 1)
  ϵᵖʳᵉˢ = max(maximum(ϵᴸ!), maximum(ϵᴿ!))
  _solver_tol = solver_tol(ϵᵖʳᵉˢ)
  ψᴴ = dag(ψ)
  ψ′ = ψᴴ'
  # XXX: make this prime the center sites
  ψ̃ = prime(linkinds, ψᴴ)

  l = linkinds(only, ψ.AL)
  l′ = linkinds(only, ψ′.AL)
  r = linkinds(only, ψ.AR)
  r′ = linkinds(only, ψ′.AR)
  s = siteinds(only, ψ)
  δʳ(n) = δ(dag(r[n]), prime(r[n]))
  δˡ(n) = δ(l[n], l′[n])
  δˢ(n) = δ(dag(s[n]), prime(s[n]))
  Hᴸ, eᴸ = left_environment(∑h, ψ; tol=_solver_tol)
  Hᴿ, eᴿ = right_environment(∑h, ψ; tol=_solver_tol)

  C̃ = InfiniteMPS(Vector{ITensor}(undef, Nsites))
  Ãᶜ = InfiniteMPS(Vector{ITensor}(undef, Nsites))
  for n in 1:Nsites
    Cvalsₙ, Cvecsₙ, Cinfoₙ = solver(Hᶜ(∑h, Hᴸ, Hᴿ, ψ, n), time_step, ψ.C[n], _solver_tol)
    Avalsₙ, Avecsₙ, Ainfoₙ = solver(
      Hᴬᶜ(∑h, Hᴸ, Hᴿ, ψ, n), time_step, ψ.AL[n] * ψ.C[n], _solver_tol
    )

    C̃[n] = Cvecsₙ
    Ãᶜ[n] = Avecsₙ
  end

  Ãᴸ = InfiniteMPS(Vector{ITensor}(undef, Nsites))
  Ãᴿ = InfiniteMPS(Vector{ITensor}(undef, Nsites))
  for n in 1:Nsites
    Ãᴸ[n] = ortho_polar(Ãᶜ[n], C̃[n])
    Ãᴿ[n] = ortho_polar(Ãᶜ[n], C̃[n - 1])
  end

  for n in 1:Nsites
    ϵᴸ![n] = norm(Ãᶜ[n] - Ãᴸ[n] * C̃[n])
    ϵᴿ![n] = norm(Ãᶜ[n] - C̃[n - 1] * Ãᴿ[n])
  end
  return InfiniteCanonicalMPS(Ãᴸ, C̃, Ãᴿ), (eᴸ, eᴿ)
end


text:
using Test
include("../src/rmstest.jl")

@testset "All Unit test for RMS" begin
for (root,dirs,files) in walkdir("../src/")
    for file in files
        if startswith(file,"Test")
            include(joinpath(root,file));
        end
    end
end
end;

ext: jl


text:
open("day06/day06.in") do io
    answers = Dict()
    answers_count = 0
    group_count = 0
    for line in readlines(io)
        line = line * " "
        if cmp(line, " ") == 0
            for (item, value) in answers
                if value == group_count
                    answers_count += 1
                end
            end
            answers = Dict()
            group_count = 0
        else
            for ch in line
                if ch != ' '
                    answers[ch] = get(answers, ch, 0) + 1
                end
            end
            group_count += 1
        end
    end
    for (item, value) in answers
        if value == group_count
            answers_count += 1
        end
    end
    println(answers_count)
end

max_stars_repo_name: detrin/advent-of-code

max_issues_repo_name: detrin/advent-of-code

max_forks_repo_name: detrin/advent-of-code


text:
module TestAvgFilter1DConstProp

using TapirBenchmarks
using Test

@testset for n in [3, 5, 100, 2^8]
    ys0, xs = avgfilter1d_constprop_setup(n)
    ys1 = avgfilter1d_constprop_seq!(copy(ys0), xs)
    @test avgfilter1d_constprop_threads!(copy(ys0), xs) == ys1
    @test avgfilter1d_constprop_tapir_dac!(copy(ys0), xs) == ys1
    @test avgfilter1d_constprop_tapir_seq!(copy(ys0), xs) == ys1
end

end


text:
using Documenter, ChipSort

makedocs(
    sitename="ChipSort.jl",
    pages = [
        "index.md",
        "theory.md",
        "algo-detail.md",
        "api.md",
        "benchmark.md"
    ],
    # format = Documenter.HTML(
    #     prettyurls = get(ENV, "CI", nothing) == "true"
    # )
)

deploydocs(
    repo = "github.com/nlw0/ChipSort.jl.git",
)


text:
function geometry(params)

       N = params.Nϕ
       L = params.Lϕ
    ϕmin = params.ϕmin

    D,y = cheb(N);
      ϕ = L / 2 * (y .+ 1);
      D = 2 / L * D;
     D2 = D*D;
      ϕ = ϕ .+ ϕmin

    return D, D2, ϕ
end


text:
module DIP3e

using Printf
using Pkg.Artifacts
const artifacts_toml = abspath(@__DIR__, "..", "Artifacts.toml")
const artifact_name = "dip3e_data"
const artifacts_url_prefix = "http://www.imageprocessingplace.com/downloads_V3/dip3e_downloads/dip3e_book_images/"
const artifact_lists = append!(map(idx->@sprintf("DIP3E_CH%02d_Original_Images.zip", idx), 1:12),
                               ["DIP3E_CH08_Original_Video_Clips.zip",
                                "DIP3E_Problem_Figures.zip"])

using ZipFile

function get_artifact_dirpath(ch_idx::Integer)
    artifact_name = @sprintf "DIP3E_CH%02d" ch_idx
    file_hash = artifact_hash(artifact_name, artifacts_toml)
    if isnothing(file_hash) || !artifact_exists(file_hash)
        new_hash = create_artifact() do artifact_dir
            zip_filename = artifact_name * "_Original_Images.zip"
            @info "downloading and unzipping dataset $(zip_filename)"
            src_url = joinpath(artifacts_url_prefix, zip_filename)
            download_path = joinpath(artifact_dir, zip_filename)
            download(src_url, download_path)
            isfile(download_path) || error("failed to download $(download_path), please try it again.")
            r = ZipFile.Reader(download_path)
            for f in r.files
                fname = basename(f.name)
                isempty(fname) && continue
                write(joinpath(artifact_dir, basename(f.name)), read(f))
            end
            close(r)
            rm(download_path)
        end
        if isnothing(file_hash)
            bind_artifact!(artifacts_toml, artifact_name, new_hash, lazy=true)
        end
        file_hash = new_hash
    end
    return joinpath(artifact_path(file_hash))
end

"""
    get_filepath(img_name)

Return the path to the image file with `img_name` in its name.

# Examples

```julia
using DIPTutorials
DIP3e.get_filepath("0101") # figure 1.1
```

See also: [`DIP3e.list_filenames`](@ref)
"""
function get_filepath(img_name::AbstractString)
    m = match(r"(\d+)", img_name)
    isnothing(m) && error("Unrecognized file name \"$(img_name)\".")
    ch_idx = parse(Int, m.captures[1][1:2])
    files = list_filenames(ch_idx; abspath=true)

    idx = findfirst(x->occursin(img_name, x), files)
    if isnothing(idx)
        error("image with \"$(img_name)\" in its name does not exist, please check possible filename with `DIP3e.list_filenames($(ch_idx))` ")
    end
    return files[idx]
end

"""
    list_filenames(idx::Integer; abspath=false)

Return a list of filenames available in chapter `idx`.

If `abspath` is `true`, then each item is a path that you can directly load.

# Examples

```julia
using DIPTutorials
DIP3e.list_filenames(1) # all figures in chapter 1
```

See also: [`DIP3e.get_filepath`](@ref)
"""
function list_filenames(ch_idx::Integer; abspath=false)
    root = get_artifact_dirpath(ch_idx)
    files = filter(x->endswith(lowercase(x), "tif"), readdir(root))
    return abspath ? joinpath.(root, files) : files
end

end


text:
"""
    Proof of concept for Caesar-ROS integration
    (check Caesar wiki for details/instructions)

    To do:
    - re-enable JSON replies
        s- periodic export of factor graph object

    To run:
    - source /opt/ros/noetic/setup.bash
    - cd ~/thecatkin_ws
        - source devel/setup.bash in all 3 terminals
    - run roscore in one terminal
    - Make sure the rosbag is in ~/data/Marine/lidar_radar.bag
    - and julia RExFeed.jl in another terminal.
"""

## Prepare python version
using Distributed
# addprocs(4)

using Pkg
Distributed.@everywhere using Pkg

Distributed.@everywhere begin
  ENV["PYTHON"] = "/usr/bin/python3"
  Pkg.build("PyCall")
end

using PyCall
Distributed.@everywhere using PyCall

## INIT
using RobotOS

# standard types
@rosimport sensor_msgs.msg: PointCloud2
@rosimport sensor_msgs.msg: NavSatFix
# @rosimport nmea_msgs.msg: Sentence
# seagrant type
# OMG, rosmsg list is awesome by the way...
# Also rosnode info
@rosimport seagrant_msgs.msg: radar

rostypegen()
# No using needed because we're specifying by full name.
# using .sensor_msgs.msg
# using .seagrant_msgs.msg

## Load Caesar mmisam stuff
using Caesar
using RoME
using DistributedFactorGraphs

using JSON2
using DistributedFactorGraphs.DocStringExtensions
using Dates

##

# /gps/fix              10255 msgs    : sensor_msgs/NavSatFix
# /gps/nmea_sentence    51275 msgs    : nmea_msgs/Sentence
# /radar_0               9104 msgs    : seagrant_msgs/radar
# /radar_pointcloud_0    9104 msgs    : sensor_msgs/PointCloud2
# /velodyne_points      20518 msgs    : sensor_msgs/PointCloud2

# function handleGPS(msg, fg)
# end

"""
    $TYPEDEF
Quick placeholder for the system state - we're going to use timestamps to align all the data.
"""
mutable struct SystemState
    curtimestamp::Float64
    cur_variable::Union{Nothing, DFGVariable}
    var_index::Int
    lidar_scan_index::Int
    max_lidar::Int
    SystemState() = new(-1000, nothing, 0, 0, 3)
end

"""
    $(SIGNATURES)
Update the system state variable if the timestamp has changed (increment variable)
"""
function updateVariableIfNeeded(fg::AbstractDFG, systemstate::SystemState, newtimestamp::Float64)
    # Make a new variable if so.
    if systemstate.curtimestamp == -1000 || systemstate.cur_variable === nothing || systemstate.curtimestamp < newtimestamp
        systemstate.curtimestamp = newtimestamp
        systemstate.cur_variable = addVariable!(fg, Symbol("x$(systemstate.var_index)"), Pose2, timestamp = unix2datetime(newtimestamp))
        systemstate.var_index += 1
        systemstate.lidar_scan_index = 0
    end
    return nothing
end

# TODO much room for improvement, and must be consolidated with RobotOS.jl
function _unpackROSMsgType(T::Type, msgdata)
    msgT = T()
    msgT.header = msgdata[2][:header]
    fnms = fieldnames(T)
    for nm in fnms
        # FIXME, still have to resolve cases where fieldname header is actually used
        nm == :header ? continue : nothing
        _tostr = JSON2.write(msgdata[2][nm])
        try
            setfield!(msgT, nm, JSON2.read(_tostr, typeof(getfield(msgT, nm))))
        catch
            @warn "not able to unpack ROS message field" T nm typeof(getfield(msgT, nm))
            @debug _tostr
        end
    end
    msgT
end

"""
    $SIGNATURES

Message callback for /radar_0.
"""
function handleRadar!(msg::seagrant_msgs.msg.radar, fg::AbstractDFG, systemstate::SystemState)
    @info "handleRadar" maxlog=10

    # if systemstate.cur_variable === nothing
    #     return nothing
    # end
    # Update the variable if needed
    timestamp = Float64(msg.header.stamp.secs) + Float64(msg.header.stamp.nsecs)/1.0e9
    updateVariableIfNeeded(fg, systemstate, timestamp)
    @info "[$timestamp] RADAR sample on $(systemstate.cur_variable.label)"

    # Make a data entry in the graph - use JSON2 to just write this (really really verbosely)
    ade,adb = addData!(fg, :radar, systemstate.cur_variable.label, :RADAR, Vector{UInt8}(JSON2.write(msg)))

end

"""
    $SIGNATURES

Message callback for Radar pings. Adds a variable to the factor graph and appends the scan as a bigdata element.
"""
function handleRadarPointcloud!(msg::sensor_msgs.msg.PointCloud2, fg::AbstractDFG, systemstate::SystemState)
    @info "handleRadarPointcloud" maxlog=10

    if systemstate.cur_variable === nothing
        # Keyed by the radar, skip if we don't have a variable yet.
        return nothing
    end
    timestamp = Float64(msg.header.stamp.secs) + Float64(msg.header.stamp.nsecs)/1.0e9
    @info "[$timestamp] RADAR pointcloud sample on $(systemstate.cur_variable.label)"

    # Make a data entry in the graph
    ade,adb = addData!(fg, :radar, systemstate.cur_variable.label, :RADARPC, Vector{UInt8}(JSON2.write(msg)),  mimeType="/radar_pointcloud_0;dataformat=Float32*[[X,Y,Z]]*12")
end

function _handleRadarPointcloud!(msg, args::Tuple)
    # sensor_msgs.msg.PointCloud2
    @info "_pointcloud"
end

"""
    $SIGNATURES

Message callback for LIDAR point clouds. Adds a variable to the factor graph and appends the scan as a bigdata element.
Note that we're just appending all the LIDAR scans to the variables because we are keying by RADAR.
"""
function handleLidar!(msg::sensor_msgs.msg.PointCloud2, fg::AbstractDFG, systemstate::SystemState)
    @info "handleLidar" maxlog=10
    # Compare systemstate and add the LIDAR scans if we want to.
    if systemstate.cur_variable === nothing
        return nothing
    end
    timestamp = Float64(msg.header.stamp.secs) + Float64(msg.header.stamp.nsecs)/1.0e9
    @info "[$timestamp] LIDAR pointcloud sample on $(systemstate.cur_variable.label) (sample $(systemstate.lidar_scan_index+1))"

    # Check if we have enough LIDAR's for this variable
    if systemstate.lidar_scan_index >= systemstate.max_lidar
        @warn "Ditching LIDAR sample for this variable, already have enough..."
        return nothing
    end

    # Make a data entry in the graph
    ade,adb = addData!(fg, :lidar, systemstate.cur_variable.label, Symbol("LIDAR$(systemstate.lidar_scan_index)"), Vector{UInt8}(JSON2.write(msg)), mimeType="/velodyne_points;dataformat=Float32*[[X,Y,Z]]*32")

    # NOTE: If JSON, then do this to get to Vector{UInt8} - # byteData = Vector{UInt8}(JSON2.write(xyzLidarF32))

    # Increment LIDAR scan count for this timestamp
    systemstate.lidar_scan_index += 1
end

"""
    $SIGNATURES

Message callback for Radar pings. Adds a variable to the factor graph and appends the scan as a bigdata element.
"""
function handleGPS!(msg::sensor_msgs.msg.NavSatFix, fg::AbstractDFG, systemstate::SystemState)
    @info "handleGPS" maxlog=10
    if systemstate.cur_variable === nothing
        # Keyed by the radar, skip if we don't have a variable yet.
        return nothing
    end
    timestamp = Float64(msg.header.stamp.secs) + Float64(msg.header.stamp.nsecs)/10^9
    # Update the variable if needed
    # updateVariableIfNeeded(fg, systemstate, timestamp)
    @info "[$timestamp] GPS sample on $(systemstate.cur_variable.label)"

    if :GPS in listDataEntries(fg, systemstate.cur_variable.label)
        @warn "GPS sample on $(systemstate.cur_variable.label) already exist, dropping"
        return nothing
    end

    io = IOBuffer()
    JSON2.write(io, msg)
    ade,adb = addData!(fg, :gps_fix, systemstate.cur_variable.label, :GPS, take!(io),  mimeType="application/json", description="JSON2.read(IOBuffer(datablob))")

end

function _handleRadar!(msgdata, args::Tuple)
    msgT = _unpackROSMsgType(seagrant_msgs.msg.radar, msgdata)
    handleRadar!(msgT, args...)
end

function _handleLidar!(msgdata, args::Tuple)
    msgT = _unpackROSMsgType(sensor_msgs.msg.PointCloud2, msgdata)
    handleLidar!(msgT, args...)
end

function _handleGPS!(msgdata, args)
    msgT = _unpackROSMsgType(sensor_msgs.msg.NavSatFix, msgdata)
    handleGPS!(msgT, args...)
end

function main(;iters::Integer=50)
    dfg_datafolder = "/tmp/caesar/rex"
    if isdir(dfg_datafolder)
        rm(dfg_datafolder; force=true, recursive=true)
    end
    mkdir(dfg_datafolder)

    @info "Hit CTRL+C to exit and save the graph..."

    init_node("rex_feed")
    # find the bagfile
    bagfile = joinpath(ENV["HOME"],"data","Marine","basin_radar.bag") # lidar_radar.bag
    bagSubscriber = RosbagSubscriber(bagfile)

    # Initialization
    fg = initfg()

    ds = FolderStore{Vector{UInt8}}(:radar, "$dfg_datafolder/data/radar")
    addBlobStore!(fg, ds)

    ds = FolderStore{Vector{UInt8}}(:gps_fix, "$dfg_datafolder/data/gps")
    addBlobStore!(fg, ds)

    # add if you want lidar also
    ds = FolderStore{Vector{UInt8}}(:lidar, "$dfg_datafolder/data/lidar")
    addBlobStore!(fg, ds)

    # System state
    systemstate = SystemState()

    # Enable and disable as needed.
    # Skipping LIDAR because those are huge...
    # radar_sub = Subscriber{seagrant_msgs.msg.radar}("/radar_0", handleRadar!, (fg, systemstate), queue_size = 10)
    radar_sub = bagSubscriber("/radar_0", _handleRadar!, (fg, systemstate))

    # radarpc_sub = Subscriber{sensor_msgs.msg.PointCloud2}("/radar_pointcloud0", handleRadarPointcloud!, (fg, systemstate), queue_size = 10)
    # lidar_sub = Subscriber{sensor_msgs.msg.PointCloud2}("/velodyne_points", handleLidar!, (fg,systemstate), queue_size = 10)
    # gps_sub = Subscriber{sensor_msgs.msg.NavSatFix}("/gps/fix", handleGPS!, (fg, systemstate), queue_size = 10)
    gps_sub = bagSubscriber("/gps/fix", _handleGPS!, (fg, systemstate))

    @info "subscribers have been set up; entering main loop"
    # loop_rate = Rate(20.0)
    while loop!(bagSubscriber)
        iters -= 1
        iters < 0 ? break : nothing
    end

    @info "Exiting"
    # After the graph is built, for now we'll save it to drive to share.
    # Save the DFG graph with the following:
    @info "Saving DFG to $dfg_datafolder/dfg"
    saveDFG(fg, "$dfg_datafolder/dfg")

end

##

main(iters=100)

# 275

## after the graph is saved it can be loaded and the datastores retrieved

dfg_datafolder = "/tmp/rex"

fg = loadDFG("$dfg_datafolder/dfg")

ds = FolderStore{Vector{UInt8}}(:radar, "$dfg_datafolder/data/radar")
addBlobStore!(fg, ds)

ds = FolderStore{Vector{UInt8}}(:gps_fix, "$dfg_datafolder/data/gps")
addBlobStore!(fg, ds)

# add if you want lidar also
ds = FolderStore{Vector{UInt8}}(:lidar, "$dfg_datafolder/data/lidar")
addBlobStore!(fg, ds)


text:
module ELPA

typealias  Graph{V,W} Dict{V, Dict{V, W}}

export addnode!, deletenode!, addedge!, deleteedge!, ne,
    lpa_addnode!, lpa_deletenode!, lpa_addedge!, lpa_deleteedge!, loadgml, savegml, modularity, Graph
# package code goes here
include("graph.jl")
include("evolving.jl")
include("graphio.jl")
include("utils.jl")

end # module


text:
using PyPlot
using PyCall
using Latexify
export plot_2D_phase_diagram_interactive

#################################################
###interactive plotting
#################################################

"Make 2-D coordinate arrays for vectorized evaluations of 2-D scalar/vector fields over 2-D grids"
meshgrid(x, y) = [i for i in x, j in 1:length(y)], [j for i in 1:length(x), j in y]

"label each location in parameter space by string stating which solutions are stable/unstable. Consider only physical solutions"
param_label(res::Result) = vec(join.(replace!.([[join.(string(el)) for el in bit_string[phys_string]]
                            for (bit_string,phys_string) in zip(res.classes["stable"],res.classes["physical"])], ("true"=>"S"),("false"=>"U"))));

"Argument corresponding to the closest value of an array to a given point"
_find_nearest(array,value) = argmin(abs.(array.-value))

_reshape_for_plot(arr,L,M,N) = reshape(reduce(hcat,reduce(hcat,arr)),L,M,N) #some useful transformation from matrices of vectors to higher dimensional tensors, for plotting

_get_var_name_labels(res::Result) = string.(res.problem.variables)

"Preprocess solution arrays to be plotted. Extract parameter information for labelling and setting axes of interactive plots"
function _get_interactive_plot_variables(res::Result,cut_type; string_f, marker_classification)
    x,y = collect(keys(res.swept_parameters))

    X,Y = collect(values(res.swept_parameters))

    gx,gy = collect(meshgrid(X, Y)) #artificial grid points
    var_names = _get_var_name_labels(res) #variable names for plot labels

    nvars     = length(res.solutions[1][1])
    nsolsmax  = length(res.solutions[1])

    sol_type, not_sol_type = _classify_plot_data(res, marker_classification)

    #prepare solution arrays for plotting
    if cut_type =="transform"
        solutions     = [transform_solutions(res, string) for string in string_f]
        physical_sols = [filter_solutions.(solution, res.classes["physical"]) for solution in solutions]

        Ys = [real.(filter_solutions.(physical_sol, res.classes[sol_type]))  for physical_sol in physical_sols]
        Yu = [real.(filter_solutions.(physical_sol, [.!el for el in res.classes[sol_type]])) for physical_sol in physical_sols]
    else
        physical_sols = filter_solutions.(res.solutions, res.classes["physical"])
        Ys = real.(filter_solutions.(physical_sols, res.classes[sol_type]))
        Yu = real.(filter_solutions.(physical_sols, [.!el for el in res.classes[sol_type]]))
    end

    return nvars,nsolsmax,Ys,Yu,x,y,X,Y,gx,gy,var_names,sol_type,not_sol_type
end

"Set up axes and plot an invisible grid of points to be used as reference for subsequent hovering labels"
function _get_interactive_plot_axes(x,y,gx,gy,var_names,cut_dim,cut_type,nvars,nsolsmax_physical,sol_type,not_sol_type; string_f)
    if  cut_type=="solutions"
        N_panels = nvars + 1
        lab = [sol_type,not_sol_type]
    elseif cut_type=="jacobian_eigenvalues"
        N_panels = nsolsmax_physical + 1
        lab = ["real","imag"]
    elseif cut_type == "transform"
        N_panels = length(string_f) + 1
        lab =   [sol_type,not_sol_type]
    end

    f,ax = plt.subplots(1,N_panels,figsize=(5*N_panels,5))

    sc = ax[1].scatter(gx',gy', s=1,alpha=0.) #create an invisible grid to attribute hovering labels, the transposes are introduce to match imshow ordering
    ax[1].set_xlabel(Latexify.latexify(x),fontsize=24)
    ax[1].set_ylabel(Latexify.latexify(y),fontsize=24)

    #hovering label preparation
    annot = ax[1].annotate("", xy=(0,0), xytext=(20,20),textcoords="offset points",
                    bbox=Dict(("boxstyle"=>"round"), ("fc"=>"w")),
                    arrowprops=Dict("arrowstyle"=>"->"))
    annot.set_visible(false)
    annot.get_bbox_patch().set_facecolor("w") #settings for hovering label box
    annot.get_bbox_patch().set_alpha(0.8)

    for l in 1:N_panels-1 #preparation of plot labels, no need to update these onclick
        if cut_dim=="1" #select direction along which solution cut_dim will be drawn
            ax[l+1].set_xlabel(Latexify.latexify(string(x)),fontsize=24);
        elseif cut_dim=="2"
            ax[l+1].set_xlabel(Latexify.latexify(string(y)),fontsize=24);
        end

        if cut_type=="solutions"
            ax[l+1].set_ylabel(Latexify.latexify(var_names[l]),fontsize=24);
        elseif cut_type=="jacobian_eigenvalues"
            ax[l+1].set_title(string("solution ", l),fontsize=12,loc="left");
        elseif cut_type=="transform"
            ax[l+1].set_title(Latexify.latexify(string_f[l]))
        end
    end
    return sc,ax,f,annot,lab
end

"Plotter of Jacobian eigenvalues along a cut of a 2D solutions. TODO: condense this with  plot_1D_jacobian_eigenvalues"
function _plot_2D_solutions_jacobian_cut(res::Result,filtered_sol,parameter_cut,Z::Vector{Float64},ax::PyCall.PyObject)

     # pre-evaluate the Jacobian with fixed parameters
    fixed_pair     = parameter_cut[2]
    fixed_params = Dict(k=>parse(ComplexF64,string(v))  for (k,v) in pairs(merge(res.fixed_parameters,Dict(fixed_pair))))
    fixed_params = Num_to_Variable(fixed_params)

    Jac = HomotopyContinuation.evaluate(HomotopyContinuation.jacobian(res.problem.system),collect(keys(fixed_params))=>collect(values(fixed_params))) #TODO make this work with HarmonicBalance.Jacobian

    swept_p      = Num_to_Variable(parameter_cut[1][1]) #swept parameter symbolic variable
    swept_values = parameter_cut[1][2] #swept parameter values

     # evaluate the numerical jacobian at the solutions. NaN elements in solutions need to be removed. These are introduced by filter_solutionss.
    Js = [HomotopyContinuation.evaluate(Jac,
            res.problem.system.variables => filtered_sol[idz,:],
            swept_p=> swept_values[idz]) for idz in 1:length(Z) if any(isnan.(filtered_sol[idz,:]).==false)]

    evals = [eigvals(J) for J in Js]
    Z_plot = [Z[idz] for idz in 1:length(Z) if any(isnan.(filtered_sol[idz,:]).==false)]
    ax.plot(Z_plot,real.(evals),"-")
    ax.plot(Z_plot,imag.(evals),"--")
end

"Take cut of filtered solutions once the user clicks on a 2D plot"
function _prepare_solution_cuts(ax::PyCall.PyObject,res::Result,Ys,Yu,cut_dim,cut_type,idx,idy,nvars,nsolsmax,X,Y)
    p1,p2 = collect(keys(res.swept_parameters))
    if cut_dim=="1"  #select direction along which solution cut_dim will be drawn
        ax.axhline(Y[idy],ls="--",c="w")
        if cut_type=="transform"
            solution_cut_s, solution_cut_u = reduce(hcat,Ys[:,idy]), reduce(hcat,Yu[:,idy])
        else
            solution_cut_s   = _reshape_for_plot(Ys[:,idy],nvars,nsolsmax,length(X))
            solution_cut_u   = _reshape_for_plot(Yu[:,idy],nvars,nsolsmax,length(X))
        end
        parameter_val  = [(p1=>res.swept_parameters[p1]), (p2=>res.swept_parameters[idy][2])]
        Z = X #cut parameter
    elseif cut_dim=="2"
        ax.axvline(X[idx],ls="--",c="w")
        if cut_type=="transform"
            solution_cut_s, solution_cut_u   = reduce(hcat,Ys[idx,:]), reduce(hcat,Yu[idx,:])
        else
            solution_cut_s  =  _reshape_for_plot(Ys[idx,:],nvars,nsolsmax,length(Y))
            solution_cut_u  =  _reshape_for_plot(Yu[idx,:],nvars,nsolsmax,length(Y))
        end
        parameter_val  = [(p2=>res.swept_parameters[p2]), (p1=>res.swept_parameters[idx][1])]
        Z = Y #cut parameter
    end
    return Z,solution_cut_s, solution_cut_u,parameter_val
end

"""
    plot_2D_phase_diagram_interactive(res::Result;
     observable="nsols",
      stable=false,nrows=2,
      ncols=2,cut_dim="1",
      cut_type="solutions",
      string_f=nothing,
      marker_classification="stable")

Interactive phase diagram of 2D solutions stored in `Result`.
This includes a clickable version of a `plot_2D_phase_diagram` for a given `observable` and extra panels containing 1D cuts of solutions, functions of solutions, or Jacobian eigenvalues.

Keyword arguments

- `observable`: reference observable to represent dynamical phases in the problem. If `observable="nsols"`, number of solutions for each point is shown.
   If instead `observable="binary"``, the result of classification of bistrings `[is_stable(solution_1),is_stable(solution_2),...]` is presented (see `classify_binaries!(Result)` function).
- `ncols`, `nrows`: number of rows and columns of the plot window.
- `cut_dim`: dimension along which 1D quantities will be calculated. `cut_dim="1"` (`cut_dim="2"`) takes a cut along the horizontal (vertical) parameter dimension of the 2D plot.
- `cut_type`: quantity to be represented along the 1D cut. If `cut_type=solutions`, steady state variables are shown with a panel per `Problem` variable.
    Else if `cut_type=jacobian eigenvalues`, Re and Im parts of complex Jacobian eigenvalues for each solution are shown  with a panel per solution.
    If instead `cut_type=transform`, functions of the solution variables passed to `string_f` (see below) are displayed.
- `string_f`: list of strings for transformed observables to be plotted in 1D when `cut_type=transform`, e.g. `string_f=["sqrt(u1^2 + v1^2)","sqrt(u2^2 + v2^2)"]` for `Problem` variables `u1,u2,v1,v2`.
- `marker_classification`: A class of the solutions (created by `classify_solutions!`) which is distinguished with different markers. Entering an inequality creates a new class "custom_class".
- `ylim`: vertical limits for 1D cuts
"""
function plot_2D_phase_diagram_interactive(res::Result; observable::String="nsols", stable::Bool=false,nrows::Int64=2,ncols::Int64=2,cut_dim::String="1",cut_type::String="solutions",string_f=nothing,marker_classification="stable",ylim=nothing)
    pygui(true) #opens a separate window

    _set_plotting_settings()
    rcParams = PyPlot.PyDict(PyPlot.matplotlib."rcParams")
    rcParams["text.usetex"] = false #cannot be set globally as it seems to conflict with pygui(true)

    names= param_label(res) #bitstring stable/unstable (S/U) labelling each physical solution.
    cut_types = ["solutions", "jacobian_eigenvalues","transform"]
    cut_type ∈ cut_types || error("Only the following types of 1D cuts are allowed:  ", cut_types)

    nvars,nsolsmax,Ys,Yu,x,y,X,Y,gx,gy,var_names,sol_type,not_sol_type = _get_interactive_plot_variables(res,cut_type,string_f=string_f,marker_classification=marker_classification)

    nsolsmax_physical = sum(any.(classify_branch(res, "physical"))) #number of physical solutions
    pretty_ylabels    = [_prettify_label(res,v)  for v in var_names]

    if !isnothing(string_f) #prettify ylabels, just for plotting, if transformed solutions are required
        string_f   = [_prettify_label(res,st) for st in string_f]
    end

    sc,ax,f,annot,lab = _get_interactive_plot_axes(x,y,gx,gy,pretty_ylabels,cut_dim,cut_type,nvars,nsolsmax_physical,sol_type,not_sol_type,string_f=string_f)

    length(vec(ax)) <= nrows*ncols || error("insufficient # of panels requested, please increase nrows or ncols") #sanity check before any plot is made

    _,im,Nmax = plot_2D_phase_diagram(res; stable=stable,observable=observable,ax=ax[1]) #skip saved data dictionary

    "Update annotations when cursor moves"
    function update_annot(ind)
        annot.xy = sc.get_offsets()[ind["ind"][1],:]
        text = names[ind["ind"][1]]
        annot.set_text(text)
    end

    "Action when mouse hovers over the figure: show string indicating stability/instability of solution"
    function hover(event)
        vis = annot.get_visible()
        if event[:inaxes] == ax[1]
            cont, ind = sc.contains(event)
            if cont
                update_annot(ind)
                annot.set_visible(true)
                f[:canvas].draw_idle()
            else
                if vis
                    annot.set_visible(false)
                    f[:canvas].draw_idle()
                end
            end
        end
    end

    f[:canvas][:mpl_connect]("motion_notify_event",hover)

    "Simple mouse click function to store coordinates and plot the corresponding cuts"
    function onclick(event)
        ix, iy = event[:xdata], event[:ydata]
        idx, idy = _find_nearest(X, ix), _find_nearest(Y, iy)  #find closest point to the mouse click.

        if length(ax[1][:lines])>0 #clear plot in current axis in case there is any
            ax[1][:lines] = []
        end

        #Z is the parameter along which cut is taken
        if cut_type!="transform"
            Z,solution_cut_s, solution_cut_u,parameter_val  = _prepare_solution_cuts(ax[1],res,Ys,Yu,cut_dim,cut_type,idx,idy,nvars,nsolsmax,X,Y) #ax[1] is the axis where obs_2D is displayed
        else
            solution_cut_s, solution_cut_u   = [], []
            for l=1:length(string_f)
                Z,sol_cut, sol_cut_u,parameter_val  = _prepare_solution_cuts(ax[1],res,Ys[l],Yu[l],cut_dim,cut_type,idx,idy,nvars,nsolsmax,X,Y)
                append!(solution_cut_s,[sol_cut])
                append!(solution_cut_u,[sol_cut_u])
            end
        end

        for l in 1:length(ax)-1
            if length(ax[l+1][:lines])>0 #clear plot in current axis in case there is any
                ax[l+1][:lines] = []
            end
        end

        if cut_type=="solutions"
            for l in 1:nvars
                ax[l+1].plot(Z, _squeeze!(solution_cut_s[l,:,:]'),lw=3)
                ax[l+1].plot(Z, _squeeze!(solution_cut_u[l,:,:]'),ls="--",lw=3)
            end
        elseif cut_type=="jacobian_eigenvalues"
            for l in 1:nsolsmax
                ax[l+1].axhline(0,ls=":",lw=4) #reference value for unstability test
                for sols in _squeeze!.([solution_cut_s[:,l,:]',solution_cut_u[:,l,:]'])
                    _plot_2D_solutions_jacobian_cut(res,sols,parameter_val,Z,ax[l+1])
                end
            end
        elseif cut_type=="transform"
            for l in 1:length(string_f)
                ax[l+1].plot(Z,_squeeze!(solution_cut_s[l]'),ls="-",lw=3)
                ax[l+1].plot(Z,_squeeze!(solution_cut_u[l]'),ls="--",lw=3)
            end
        end

        if !isnothing(ylim)
            for (l,ax_) in enumerate(ax)
                if l>1
                    ax_.set_ylim(ylim)
                end
            end
        end

        legend_elements = [plt.Line2D([1], [1], linestyle="-", color="k", label=lab[1], markerfacecolor="k", markersize=5),
                           plt.Line2D([1], [1], linestyle="--", color="k", label=lab[2],markerfacecolor="k", markersize=5)]
        ax[end].legend(handles=legend_elements,loc="best")
    end

    f[:canvas][:mpl_connect]("button_press_event", onclick)

    on_close(event) = pygui(false) #When closing the interactive plotting window, subsequent figures are static by default

    f[:canvas][:mpl_connect]("close_event",on_close)

    ax = resize_axes!(f,ax,nrows,ncols)
    if nrows>1
        _prepare_colorbar(f,ax,im,Nmax)
    else
        _prepare_colorbar(f,ax[end],im,Nmax)
    end
end


text:
"""
```julia
arguments(call_expr::Expr)
```
For internal use.

Retrieves the argument expressions from a call expression.
"""
function arguments(call_expr::Expr)
  if call_expr.head != :call
    throw(ErrorException("Malformed function declaration"))
  end
  call_expr = copy(call_expr)
  call_expr.head = :tuple
  shift!(call_expr.args)
  return call_expr
end

"""
```julia
GtkBuilderAidData(builder::GtkBuilder, handlers::FunctionTable)

```

"""
mutable struct GtkBuilderAidData
  builder::Union{GtkBuilder, Void}
  handlers::FunctionTable
  userdata
  wpipe::IO
end

"""
```julia
aid(builder::Union{GtkBuilder, Void) = aid.builder, userdata = aid.userdata; wpipe::IO = aid.wpipe)
```

The callable form of GtkBuilderAidData

"""
function (x::GtkBuilderAidData)(builder::Union{GtkBuilder, Void} = x.builder, userdata = x.userdata; wpipe::IO = x.wpipe)
  if isa(builder, Void)
    throw(MethodError("GtkBuilderAidData object must have builder to be callable without builder object"))
  else
    connect_signals(builder, x.handlers, userdata; wpipe = wpipe)
  end
  builder
end

function (x::GtkBuilderAidData)(builderFile::String, userdata = x.userdata; wpipe::IO = x.wpipe)
  builder = GtkBuilder(filename = builderFile)
  x(builder, userdata; wpipe = wpipe)
  builder
end

"""
```julia
ApplicationUserdata(app::GtkApplication, builder::GtkBuilder, userdata)
```

"""
mutable struct ApplicationUserdata
  builder::GtkBuilder
  handlers::FunctionTable
  app::GtkApplication
  mainWindowId::String
  userdata
  wpipe::IO
end

# A helper callback function for application application start functionality
@guarded function qsActivateApp(widget, qsdata)
  connect_signals(qsdata.builder, qsdata.handlers, qsdata; wpipe = qsdata.wpipe)
  win = Gtk.GAccessor.object(qsdata.builder, qsdata.mainWindowId)
  push!(qsdata.app, win)
  showall(win)
  nothing
end

start_application(handlers::FunctionTable, appName::String, mainWindowId::String, builder::String, userdata; wpipe = Base.STDERR) =
  start_application(handlers, appName, mainWindowId, GtkBuilder(filename = builder), userdata; wpipe = wpipe)

function start_application(handlers::FunctionTable, appName::String, mainWindowId::String, builder::GtkBuilder, userdata; wpipe = Base.STDERR)
  app = GtkApplication(appName, 0)
  qsdata = ApplicationUserdata(builder, handlers, app, mainWindowId, userdata, wpipe)
  signal_connect(qsActivateApp, app, :activate, Void, (), false, qsdata)
  Gtk.register(app)
  # I forgot about this wonkiness, printing to stdout here is actually a necessary step
  println(join(("Starting Application:", appName), " "))
  run(app)
  println("Application Completed!")
end

"""
@GtkFunctionTable begin
  # ... callback functions ...
end

This macro allows using glade with Julia mostly as easy as working with Glade
in C is when coupled with `connect_signals` and the GtkBuilder object. The
macro creates a dictionary of all the callbacks that serves as a replacement
for the C function table that Gtk pulls callbacks from during usage with C.
"""
macro GtkFunctionTable(args...)
  length(args) < 1 && throw(ArgumentError("GtkBuilderAid macro requires at least one argument"))

  user_block = args[end]::Expr
  user_block.head != :block && throw(ArgumentError("The last argument to this macro must be a block"))

  # Analogous to function declarations of a C header file
  callbacks = Set{Symbol}();

  for entry in user_block.args

    if isa(entry, Expr)
      if entry.head == :macrocall
        # expand the macro
        expanded = macroexpand(entry)
        if expanded.head == :function
          entry.head = expanded.head
          entry.args = expanded.args
        end
      end
      if entry.head == :function && length(entry.args[1].args) >= 3
        fcall = entry.args[1]
        fname = fcall.args[1]
        push!(callbacks, fname)
      end

      if entry.head == :(=)
        left = entry.args[1]
        if isa(left, Expr) && left.head == :call && length(left.args) >= 3
          fname = left.args[1]
          push!(callbacks, fname)
        end
      end
    end
  end

  funcdata = Expr(:vect)
  for fname in callbacks
    push!(funcdata.args, Expr(:tuple, fname, string(fname)))
  end

  # Escape the modified user block
  quote
    $(esc(user_block))

    handlers = FunctionTable()
    for func in $(esc(funcdata))
      handlers[string(func[2])] = func[1]
    end

    handlers
  end

end

"""
```julia
@GtkBuilderAid directives()... begin
  # ... implementation ...
end
```

This macro is meant to make using glade with Julia as easy as working with
Glade in C is. From a staticly compiled language the function names are just
pulled from a compiled file but that option isn't available in Julia, at least
not cleanly.
"""
macro GtkBuilderAid(args...)
  length(args) < 1 && throw(ArgumentError("GtkBuilderAid macro requires at least one argument"))

  user_block = args[end]::Expr
  user_block.head != :block && throw(ArgumentError("The last argument to this macro must be a block"))

  # The whole directives thing was a mistake in retrospect
  lastDirective = length(args) - 1
  builder = nothing
  if length(args) >= 2 && isa(args[end - 1], AbstractString)
    lastDirective -= 1
    # Enables the pre-bound version
    filename = args[end - 1]
    if !isfile(filename)
      throw(ErrorException("Provided UI file, $filename, does not exist"))
    end
    builder = :(GtkBuilder(filename = $filename))
  end

  userdata = ()
  final_function_name = :genned_function
  for directive in args[1:lastDirective]
    # Only support expression-style directives
    if typeof(directive) <: Expr && directive.head == :call
      # A function call style directive
      if directive.args[1] == :function_name
        final_function_name = esc(directive.args[2])
      end

      if directive.args[1] == :userdata
        userdata = arguments(directive).args[1]
      end

      if directive.args[1] == :userdata_tuple
        # Creates a tuple from the arguments
        # and uses that as the userinfo argument
        userdata = arguments(directive)
      end

    else
      throw(ErrorException("Directives must be in the format of a function call"))
    end
  end

  block = quote
    handlers = @GtkFunctionTable $(esc(user_block))
    GtkBuilderAidData($builder, handlers, $(esc(userdata)), Base.STDERR)
  end

  quote
    $final_function_name = $block
  end
end


text:
module SentinelMissings
struct SentinelMissing{T,SV} <: Number
    x::T
end
Base.promote_rule(::Type{SentinelMissing{SM,SV}}, ::Union{T,Missing}) where {SM,SV,T} =
    Union{promote_type(SM, T),Missing}
Base.promote_rule(::Type{SentinelMissing{SM,SV}}, ::Type{T}) where {SM,SV,T} =
    Union{promote_type(SM, T),Missing}
Base.promote_rule(
    ::Type{SentinelMissing{SM,SV}},
    ::Type{SentinelMissing{SM,SV}},
) where {SM,SV} = Union{SM,Missing}
Base.convert(::Type{Union{T,Missing}}, sm::SentinelMissing) where {T} =
    ismissing(sm) ? missing : convert(T, sm[])
Base.convert(::Type{T}, sm::SentinelMissing) where {T<:Number} = convert(T, sm[])
Base.getindex(x::SentinelMissing) = ismissing(x) ? missing : x.x
Base.convert(::Type{Any}, sm::SentinelMissing) = sm[]
Base.ismissing(x::SentinelMissing) = isequal(x.x, senval(x))
Base.nonmissingtype(::Type{SentinelMissing{T,SV}}) where {T,SV} = T
Base.isless(x::SentinelMissing, y::SentinelMissing) = isless(x[], y[])
Base.isless(x::SentinelMissing, y::Number) = isless(x[], y)
Base.isless(y::Number, x::SentinelMissing) = isless(y, x[])
function Base.convert(::Type{<:T}, x::Number) where {T<:SentinelMissing}
    sv = senval(T)
    et = eltype(T)
    SentinelMissing{et,sv}(convert(et, x))
end
Base.convert(::Type{<:T}, x::Missing) where {T<:SentinelMissing} =
    SentinelMissing{eltype(T),senval(T)}(senval(T))
for op in (:(+), :(-), :(*), :(/), :(^))
    eval(
        :(
            Base.$(op)(x1::T, x2::T) where {T<:SentinelMissing{SM}} where {SM} =
                $(op)(convert(Union{SM,Missing}, x1), convert(Union{SM,Missing}, x2))
        ),
    )
end
Base.zero(T::Type{<:SentinelMissing{SM}}) where {SM} = T(zero(SM))
Base.one(T::Type{<:SentinelMissing{SM}}) where {SM} = T(one(SM))
Base.convert(::Type{<:T}, x::T) where {T<:SentinelMissing} = x
senval(::SentinelMissing{<:Any,SV}) where {SV} = SV
senval(::Type{<:SentinelMissing{<:Any,SV}}) where {SV} = SV
Base.eltype(::Type{<:SentinelMissing{T}}) where {T} = T
Base.show(io::IO, x::SentinelMissing) = print(io, x[])
Base.similar(x::AbstractArray{<:SentinelMissing{T,SV}}) where {T,SV} =
    zeros(Union{T,Missing}, size(x))
Base.similar(
    x::AbstractArray{<:SentinelMissing{T,SV}},
    dims::Tuple{Vararg{Int64,N}},
) where {T,SV,N} = zeros(Union{T,Missing}, dims)

"""
    as_sentinel(x, v)

Reinterprets a Number Array or a Number `x` so that values in x that equal v will be treated as missing.
This is done by reinterpreting the array as a `SentinelMissing` without copying the data.
"""
as_sentinel(x::AbstractArray{T}, v) where {T<:Number} =
    reinterpret(SentinelMissing{T,convert(T, v)}, x)
as_sentinel(x::T, v) where {T<:Number} = SentinelMissing{T,convert(T, v)}(x)
export as_sentinel
end


text:
include("Utils.jl")

struct Expanded
    a
    b
    c
end

function tofactored(F::Expanded)
    _validitycheck

    delta = calcdelta(F.a, F.b, F.c)

    x1 = (-F.b - sqrt(delta)) / 2F.a
    x2 = (-F.b + sqrt(delta)) / 2F.a

    return Factored(F.a, x1, x2)
end

function tovertex(F::Expanded)
    _validitycheck

    delta = calcdelta(F.a, F.b, F.c)

    p = -F.b / 2F.a
    q = -delta / 4F.a

    return Vertex(F.a, p, q)
end

function calcvalue(F::Expanded, x)
    return F.a * x^2 + F.b * x + F.c
end


text:
module navigation_system

using ReachabilityAnalysis, ModelingToolkit

const var = @variables x, y, v_x, v_y

function mode0()
    A = [0 0   1    0 ;
         0 0   0    1 ;
         0 0 -1.2  0.1;
         0 0  0.1 -1.2]
    B = [0, 0, -0.1, 1.2]
    invariant = UnionSetArray([HalfSpace(x <= 1),
                               HalfSpace(v_x <= 0),
                               HalfSpace(y <= 1),
                               HalfSpace(v_y <= 0)])
    return @system(x' = Ax + B, x ∈ invariant)
end

function mode1()
    A = [0 0   1    0 ;
         0 0   0    1 ;
         0 0 -1.2  0.1;
         0 0  0.1 -1.2]
    B = [0, 0, -4.8, 0.4]
    invariant = UnionSetArray([HalfSpace(x >= 1),
                               HalfSpace(v_x <= 0),
                               HalfSpace(y <= 1),
                               HalfSpace(v_y <= 0)])
    return @system(x' = Ax + B, x ∈ invariant)
end

function mode2()
    A = [0 0   1    0 ;
         0 0   0    1 ;
         0 0 -1.2  0.1;
         0 0  0.1 -1.2]
    B = [0, 0, 2.4, -0.2]
    invariant = UnionSetArray([HalfSpace(x <= 1),
                               HalfSpace(v_x <= 0),
                               HalfSpace(y >= 1),
                               HalfSpace(v_y >= 0)])
    @system(x' = Ax + B, x ∈ invariant)
end

function mode3()
    A = [0 0   1    0 ;
         0 0   0    1 ;
         0 0 -1.2  0.1;
         0 0  0.1 -1.2]
    B = [0, 0, 3.9, -3.9]
    invariant = UnionSetArray([HalfSpace(x >= 1),
                               HalfSpace(v_x >= 0),
                               HalfSpace(y >= 1),
                               HalfSpace(v_y >= 0)])
    @system(x' = Ax + B, x ∈ invariant)
end

function navigation_system_hybrid()
    automaton = LightAutomaton(4)
    add_transition!(automaton, 1, 2, 1)
    add_transition!(automaton, 2, 1, 2)
    add_transition!(automaton, 2, 4, 3)
    add_transition!(automaton, 4, 2, 4)
    add_transition!(automaton, 1, 3, 5)
    add_transition!(automaton, 3, 1, 6)
    add_transition!(automaton, 3, 4, 7)
    add_transition!(automaton, 4, 3, 8)

    modes = [mode0(), mode1(), mode2(), mode3()]

    ## transition 1 -> 2
    guard = HPolyhedron([HalfSpace(x >= 1), HalfSpace(v_x >= 0)])
    trans1 = ConstrainedIdentityMap(4, guard)
    ## transition 2 -> 1
    guard = HPolyhedron([HalfSpace(x <= 1), HalfSpace(v_x <= 0)])
    trans2 = ConstrainedIdentityMap(4, guard)
    ## transition 2 -> 4
    guard = HPolyhedron([HalfSpace(y >= 1), HalfSpace(v_y >= 0)])
    trans3 = ConstrainedIdentityMap(4, guard)
    ## transition 4 -> 2
    guard = HPolyhedron([HalfSpace(y <= 1), HalfSpace(v_y <= 0)])
    trans4 = ConstrainedIdentityMap(4, guard)
    ## transition 1 -> 3
    guard = HPolyhedron([HalfSpace(y >= 1), HalfSpace(v_y >= 0)])
    trans5 = ConstrainedIdentityMap(4, guard)
    ## transition 3 -> 1
    guard = HPolyhedron([HalfSpace(y <= 1), HalfSpace(v_y <= 0)])
    trans6 = ConstrainedIdentityMap(4, guard)
    ## transition 3 -> 4
    guard = HPolyhedron([HalfSpace(x >= 1), HalfSpace(v_x >= 0)])
    trans7 = ConstrainedIdentityMap(4, guard)
    ## transition 4 -> 3
    guard = HPolyhedron([HalfSpace(x <= 1), HalfSpace(v_x <= 0)])
    trans8 = ConstrainedIdentityMap(4, guard)
    resetmaps = [trans1, trans2, trans3, trans4, trans5, trans6, trans7, trans8]

    return HybridSystem(automaton, modes, resetmaps, [AutonomousSwitching()])
end

function model(X0)
    H = navigation_system_hybrid()
    return IVP(H, [(1, X0)])
end

end  # module

max_stars_repo_name: mforets/ReachabilityLibrary

max_issues_repo_name: mforets/ReachabilityLibrary

max_forks_repo_name: mforets/ReachabilityLibrary


text:
using EMISApproximateEquilibrium
const SR = EMISApproximateEquilibrium
using Combinatorics
using PowerSystems
const PSY = PowerSystems

build = Dict{Type{<:PSY.Component}, Array{AbstractString}}()
build[PSY.ThermalStandard] =  repeat( ["118_CC_1", "101_CT_1", "313_CC_1", "207_CT_1", "101_STEAM_4"], 3);
build[PSY.HydroEnergyReservoir] = ["122_HYDRO_5"];
build[PSY.RenewableDispatch] = repeat(["309_WIND_1", "101_PV_4", "313_PV_1", "324_PV_3"], 3);
build[PSY.GenericBattery] = repeat(["313_STORAGE_1"], 1)
future_config_dict = Dict();
buildable_devices = Dict()
buildable_devices[PSY.ThermalStandard] =
    SR.all_combinations(build, PSY.ThermalStandard);
buildable_devices[PSY.RenewableDispatch] =
    SR.all_combinations(build, PSY.RenewableDispatch);
buildable_devices[PSY.HydroEnergyReservoir] =
    SR.all_combinations(build, PSY.HydroEnergyReservoir);
buildable_devices[PSY.GenericBattery] =
    SR.all_combinations(build, PSY.GenericBattery);

future = SR.combvec_samples(10000,
    buildable_devices[PSY.ThermalStandard],
    buildable_devices[PSY.HydroEnergyReservoir],
    buildable_devices[PSY.RenewableDispatch],
    buildable_devices[PSY.GenericBattery],
)

build_no = 1
exp_name = "10k"
future_config_dict, sample_count = SR.update_config(future_config_dict, future, 1);
rm("./sample_configuration_$(exp_name)_$(build_no)/", force = true, recursive = true)
if !isdir(joinpath(pwd(), "sample_configuration_$(exp_name)_$(build_no)"))
    sample_configuration_dir = mkdir(joinpath(pwd(), "sample_configuration_$(exp_name)_$(build_no)"))
else
    sample_configuration_dir = (joinpath(pwd(), "sample_configuration_$(exp_name)_$(build_no)"))
end
SR.generate_configuration(sample_configuration_dir, future_config_dict)

include("./generate_jobs.jl")
generate_jobs(pwd(), "/home/sdalvi/work/RTS-GMLC", sample_count, build_no, exp_name)

#=
Jade commands to run
jade auto-config generic_command jobs.txt -c config.json
jade submit-jobs config.json \
    --output=output \
    --max-nodes=64 \
    --per-node-batch-size=1 \
    --hpc-config=hpc_config.toml
=#


text:
using PowerSpectra
using Healpix
using CSV
using Test
using LinearAlgebra
using DataFrames
using DelimitedFiles
using NPZ
using IdentityRanges

##
@testset "Mode Coupling Matrix TT" begin
    nside = 256
    mask = readMapFromFITS("data/example_mask_1.fits", 1, Float64)
    flat_beam = SpectralVector(ones(3*nside))
    flat_map = HealpixMap{Float64, RingOrder}(ones(nside2npix(nside)) )
    M = mcm(:TT, mask, mask; lmin=2)
    reference = readdlm("data/mcm_TT_diag.txt")
    @test all(reference .≈ diag(parent(M))[1:end-1])
    map1 = readMapFromFITS("data/example_map.fits", 1, Float64)
    pCl = SpectralVector(alm2cl(map2alm(map1 * mask)))
    Cl_hat = M \ pCl[IdentityRange(2:end)]
    reference_spectrum = readdlm("data/example_TT_spectrum.txt")
    @test all(reference_spectrum .≈ Cl_hat[2:end])
    @test all(reference_spectrum .≈ Cl_hat[2:end])
end

##
@testset "Mode Coupling Matrix Diag EE" begin
    nside = 256
    mask = readMapFromFITS("data/example_mask_1.fits", 1, Float64)
    flat_mask = HealpixMap{Float64, RingOrder}(ones(nside2npix(nside)) )
    M = mcm(:M⁺⁺, mask, mask)
    # factorized_mcm12 = lu(parent(M))
    reference = readdlm("data/mcm_EE_diag.txt")
    @test all(reference .≈ diag(parent(M))[3:767])
end

##
@testset "Mode Coupling Matrix Diag TE/ET" begin
    nside = 256
    mask = readMapFromFITS("data/example_mask_1.fits", 1, Float64)
    M = mcm(:TE, mask, mask)
    reference = readdlm("data/mcm_TE_diag.txt")
    @test all(reference .≈ diag(parent(M))[3:767])

    M = mcm(:ET, mask, mask)
    reference = readdlm("data/mcm_TE_diag.txt")
    @test all(reference .≈ diag(parent(M))[3:767])
end

##
@testset "Full Non-Trivial MCM" begin
    nside = 256
    mask1_T = readMapFromFITS("data/mask1_T.fits", 1, Float64)
    mask2_T = readMapFromFITS("data/mask2_T.fits", 1, Float64)
    mask1_P = readMapFromFITS("data/mask1_P.fits", 1, Float64)
    mask2_P = readMapFromFITS("data/mask2_P.fits", 1, Float64)
    unit_map = HealpixMap{Float64, RingOrder}(ones(nside2npix(nside)) )
    unit_beam = spectralones(0:(3nside-1))

    M = mcm(:TT, mask1_T, mask2_T)
    M_ref = npzread("data/mcmTT.npy")
    @test all(isapprox(parent(M)[3:end, 3:end], M_ref[3:end, 3:end], atol=1e-11))

    M = mcm(:TE, mask1_T, mask2_P)
    M_ref = npzread("data/mcmTE.npy")
    for k in 0:3nside
        @test all(isapprox(diag(parent(M), k)[3:end], diag(M_ref, k)[3:end]))
    end

    M = mcm(:ET, mask1_P, mask2_T)
    M_ref = npzread("data/mcmET.npy")

    for k in 0:3nside
        @test all(isapprox(diag(parent(M), k)[3:end], diag(M_ref, k)[3:end]))
    end

    M = mcm(:M⁺⁺, mask1_P, mask2_P)
    M_ref = npzread("data/mcmEE.npy")
    for k in 0:3nside
        @test all(isapprox(diag(parent(M), k)[3:end], diag(M_ref, k)[3:end]))
    end
end

max_stars_repo_name: xzackli/PlanckCov

max_issues_repo_name: xzackli/PlanckCov


text:
module TSA

using Statistics
using Distributions
using LinearAlgebra
using Random
using Flux

function DFT(x)
    N = length(x)
    X̂ = zeros(Complex,N)
    μ = collect(Int64,0:1:N-1)
    for i = 1:N
        X̂ .+= x[i]*(exp(-2π*im*(i-1)/N)).^μ
    end
    return X̂
end

function IDFT(x)
    N = length(x)
    X = zeros(Complex,N)
    μ = collect(Int64,0:1:N-1)
    for i = 1:N
        X .+= x[i]/N*(exp(2π*im*(i-1)/N)).^μ
    end
    return real(X)
end

function hamming(x)
    N = length(x)
    n = collect(Int64,0:1:N-1)
    wind = map(x->x+0.54,-0.46*cos.(2π*n/(N-1)))
    return wind
end

function hanning(x)
    N = length(x)
    n = collect(Int64,0:1:N-1)
    wind = map(x->x+0.5,-0.5*cos.(2π*n/(N-1)))
    return wind
end

function embedding(x,m,τ)
    N = length(x)
    X = zeros(N-(m-1)τ,m)
    for i = 1:m
        X[1:N-(m-1)τ,i] .= x[1+(i-1)τ:N-(m-1)τ+(i-1)τ]
    end
    return X
end

function FNN(x,τ,m)
    E = []
    for i = 1:m+1
        X = embedding(x,i,τ)
        a = zeros(size(X)[1])
        mindist = zeros(size(X)[1])
        for j = 1:size(X)[1]
            if mindist[j] == 0
                d = zeros(size(X)[1])
                for k = 1:size(X)[1]
                    d[k] = norm(X[j,:]-X[k,:])
                end
                d[j] = Inf
                k = argmin(d)
                mindist[j] = d[k]
                mindist[k] = d[k]
            end
        end

        X2 = embedding(x,i+1,τ)
        mindist2 = zeros(size(X2)[1])
        for j = 1:size(X2)[1]
            if mindist2[j] == 0
                d = zeros(size(X2)[1])
                for k = 1:size(X2)[1]
                    d[k] = norm(X2[j,:]-X2[k,:])
                end
                d[j] = Inf
                k = argmin(d)
                mindist2[j] = d[k]
                mindist2[k] = d[k]
            end
        end

        a = mindist2./mindist[1:size(X2)[1]]
        E = push!(E,sum(a)/size(X2)[1])
    end
    E1 = []
    E1 = E[2:length(E)]./E[1:length(E)-1]
    return E1
end

function TE(x,τ,m,k,q,M)
    Etrans = []
    for i = 1:m
        X = embedding(x,i,τ)
        d =zeros(size(X)[1]-τ,k+1)
        for j = 1:size(X)[1]-τ
            dist = zeros(size(X)[1]-τ)
            for l = 1:size(X)[1]-τ
                dist[l] = norm(X[j,:]-X[l,:])
            end

            d[j,1:k+1] .= sortperm(dist)[1:k+1]
        end
        Et = []
        d = Int.(d)
        for j = 1:q
            rn = map(Int,rand(1:size(X)[1]-τ,M))
            E = []
            for n in rn
                v = X[d[n,:].+τ,:] .- X[d[n,:],:]
                v̄ = sum(v,dims=1)/(k+1)
                E = push!(E,sum(abs2,(v.-v̄))/((k+1)*sum(abs2,v̄)))
            end
            Et = push!(Et,median(E))
        end
        Etrans = push!(Etrans,mean(Et))
    end
    return Etrans
end

function MaxLyapunov(x,m,τ,itr,k,meth=0,epoch=100)
    σ = std(x)
    X = embedding(x,m,τ)
    N = size(X)[1]-1
    ϵ1 = 0
    if meth == 0
        X̂ = copy(X[1:N,:])
        for j = 1:itr
            d = zeros(BigInt,N,k)
            for i = 1:N
                dist = zeros(N)
                for l = 1:N
                    dist[l] = norm(X̂[i,:]-X[l,:])
                end
                d[i,1:k] = sortperm(dist)[2:k+1]
                X̂[i,:] .= sum(X[d[i,:].+1,:])/k
            end
            if j==1
                ϵ1 = sqrt(sum((X[1+j:end,1].-X̂[1:end,1]).^2)/N)/σ
            end
        end
        ϵ = sqrt(sum((X[1+itr:end,1].-X̂[1:end-itr+1,1]).^2)/(size(X)[1]-itr))/σ
    elseif meth==1
        X̂ = copy(X[1:end-1,:])
        m2 = 2*m
        d1 = Dense(m,m2,sigmoid) |>f64
        d2 = Dense(m2,m) |>f64
        model = Chain(d1,d2)
        loss(x,x̂) = Flux.Losses.mse(model(x),x̂)
        ps = params(model)
        opt = Descent()
        xt = []
        yt = []
        for i = 1:N
            xt = push!(xt,X̂[i,:])
            yt = push!(yt,X[i+1,:])
        end
        for i = 1:epoch
            data = zip(xt,yt)
            Flux.train!(loss,ps,data,opt)
        end
        for j = 1:itr
            for i = 1:N-j
                X̂[i,:] = model(X̂[i,:])
            end
            if j==1
                ϵ1 = sqrt(sum((X[1+j:end,1].-X̂[1:N-j+1,1]).^2)/N)/σ
            end
        end
        ϵ = sqrt(sum((X[1+itr:end,1].-X̂[1:end-itr+1,1]).^2)/(size(X)[1]-itr))/σ
    end
    λ = log(ϵ/ϵ1)/(itr-1)
    return λ
end

function surrogate(x,m,τ,itr,c)
    N = length(x)
    E = []
    E0 = MaxLyapunov(x,m,τ,itr,1,1)
    for i = 1:c
        sX = copy(x)
        fsX = DFT(sX)
        fsX = fsX.*(exp(2*π*im)).^rand(Float64,N)
        sX = IDFT(fsX)
        sX = real(sX)
        E1 = MaxLyapunov(sX,m,τ,itr,1,1)
        E = push!(E,E1)
        end
    t = (E0-mean(E))/(std(E)/(c^0.5))
    return t
end
end


text:
"""
    Magic

`DataType` used for storing magic numbers, rook or bishop move masks, and respective shift and offset values for the computation of sliding moves.
"""
struct Magic
    mask::Bitboard
    num::Bitboard
    shift::UInt8
    #offset::Int32
end
Magic(mask::UInt64, num::UInt64, shift::Integer, offset::Integer) = Magic(Bitboard(mask), Bitboard(num), shift)#, offset)

const ROOK_MAGICS = @SVector [Magic(0x000101010101017e, 0x008000400020801a, 52, 0), Magic(0x000202020202027c, 0x0840004020001003, 53, 4096),
    Magic(0x000404040404047a, 0x8880200010018108, 53, 6144), Magic(0x0008080808080876, 0x0480040800801001, 53, 8192),
    Magic(0x001010101010106e, 0x0900080002050010, 53, 10240), Magic(0x002020202020205e, 0x0200100804010200, 53, 12288),
    Magic(0x004040404040403e, 0x0200010084080200, 53, 14336), Magic(0x008080808080807e, 0x0200024100821224, 52, 16384),
    Magic(0x0001010101017e00, 0x08a0802080004000, 53, 20480), Magic(0x0002020202027c00, 0x0006401004402000, 54, 22528),
    Magic(0x0004040404047a00, 0x0800802000100080, 54, 23552), Magic(0x0008080808087600, 0x4004800803100081, 54, 24576),
    Magic(0x0010101010106e00, 0x0106000810200600, 54, 25600), Magic(0x0020202020205e00, 0x2981000802040100, 54, 26624),
    Magic(0x0040404040403e00, 0x0004000804018210, 54, 27648), Magic(0x0080808080807e00, 0x0201000192026300, 53, 28672),
    Magic(0x00010101017e0100, 0x1038218000804000, 53, 30720), Magic(0x00020202027c0200, 0x1010004000200040, 54, 32768),
    Magic(0x00040404047a0400, 0x0108420012002881, 54, 33792), Magic(0x0008080808760800, 0x2098008008100080, 54, 34816),
    Magic(0x00101010106e1000, 0x4088010004090010, 54, 35840), Magic(0x00202020205e2000, 0x2000808002000400, 54, 36864),
    Magic(0x00404040403e4000, 0x028a0c0008021025, 54, 37888), Magic(0x00808080807e8000, 0x0000220000810064, 53, 38912),
    Magic(0x000101017e010100, 0x0040400480208000, 53, 40960), Magic(0x000202027c020200, 0x8040200040100042, 54, 43008),
    Magic(0x000404047a040400, 0x2000100080802000, 54, 44032), Magic(0x0008080876080800, 0x3090004040080401, 54, 45056),
    Magic(0x001010106e101000, 0x0100080080800400, 54, 46080), Magic(0x002020205e202000, 0x8230020080800400, 54, 47104),
    Magic(0x004040403e404000, 0x0002000200080104, 54, 48128), Magic(0x008080807e808000, 0x0c00010200004084, 53, 49152),
    Magic(0x0001017e01010100, 0x0000400082800020, 53, 51200), Magic(0x0002027c02020200, 0x0400200080804007, 54, 53248),
    Magic(0x0004047a04040400, 0x8400801000802001, 54, 54272), Magic(0x0008087608080800, 0x0018100080800800, 54, 55296),
    Magic(0x0010106e10101000, 0x0810310005004800, 54, 56320), Magic(0x0020205e20202000, 0x8000020080800400, 54, 57344),
    Magic(0x0040403e40404000, 0x2000302134000208, 54, 58368), Magic(0x0080807e80808000, 0x0000298402000645, 53, 59392),
    Magic(0x00017e0101010100, 0x0000800140018022, 53, 61440), Magic(0x00027c0202020200, 0x00184020100c4000, 54, 63488),
    Magic(0x00047a0404040400, 0x6210002804002000, 54, 64512), Magic(0x0008760808080800, 0x0001000c10010020, 54, 65536),
    Magic(0x00106e1010101000, 0x4201000800110004, 54, 66560), Magic(0x00205e2020202000, 0x2002000810020004, 54, 67584),
    Magic(0x00403e4040404000, 0xa40c21480a040090, 54, 68608), Magic(0x00807e8080808000, 0x0000004081020004, 53, 69632),
    Magic(0x007e010101010100, 0x3302410080002300, 53, 71680), Magic(0x007c020202020200, 0x9210002000401040, 54, 73728),
    Magic(0x007a040404040400, 0x0092110041200300, 54, 74752), Magic(0x0076080808080800, 0x0100100109002300, 54, 75776),
    Magic(0x006e101010101000, 0x000a800800240280, 54, 76800), Magic(0x005e202020202000, 0x0100040080020080, 54, 77824),
    Magic(0x003e404040404000, 0x0000506201080400, 54, 78848), Magic(0x007e808080808000, 0x1002209110440600, 53, 79872),
    Magic(0x7e01010101010100, 0x0140800300182241, 52, 81920), Magic(0x7c02020202020200, 0x8280184280220102, 53, 86016),
    Magic(0x7a04040404040400, 0x0005118042000a22, 53, 88064), Magic(0x7608080808080800, 0x008a000890204006, 53, 90112),
    Magic(0x6e10101010101000, 0x000200304844204a, 53, 92160), Magic(0x5e20202020202000, 0x108200080930041a, 53, 94208),
    Magic(0x3e40404040404000, 0x0482000100840842, 53, 96256), Magic(0x7e80808080808000, 0x0012192404430182, 52, 98304)]

const BISHOP_MAGICS = @SVector [Magic(0x0040201008040200, 0x4404700420508600, 58, 0), Magic(0x0000402010080400, 0x9120018401104008, 59, 64),
    Magic(0x0000004020100a00, 0x4004010401100084, 59, 96), Magic(0x0000000040221400, 0x0051040080062000, 59, 128),
    Magic(0x0000000002442800, 0x1124042010010092, 59, 160), Magic(0x0000000204085000, 0x1100882108001000, 59, 192),
    Magic(0x0000020408102000, 0x180a280248040000, 59, 224), Magic(0x0002040810204000, 0x0038440210900400, 58, 256),
    Magic(0x0020100804020000, 0x80074210e1010301, 59, 320), Magic(0x0040201008040000, 0x0020101012004e42, 59, 352),
    Magic(0x00004020100a0000, 0x8210100110411001, 59, 384), Magic(0x0000004022140000, 0x00008404108010a0, 59, 416),
    Magic(0x0000000244280000, 0x2060071040002804, 59, 448), Magic(0x0000020408500000, 0x6400020202211011, 59, 480),
    Magic(0x0002040810200000, 0x000c0d0111202000, 59, 512), Magic(0x0004081020400000, 0x8004050118020221, 59, 544),
    Magic(0x0010080402000200, 0x8005804088080910, 59, 576), Magic(0x0020100804000400, 0x6004a82001020200, 59, 608),
    Magic(0x004020100a000a00, 0x0101001004002041, 57, 640), Magic(0x0000402214001400, 0x8a1800a082004282, 57, 768),
    Magic(0x0000024428002800, 0x8806001012100098, 57, 896), Magic(0x0002040850005000, 0x0022810040504000, 57, 1024),
    Magic(0x0004081020002000, 0x4040800048480801, 59, 1152), Magic(0x0008102040004000, 0x0a42041304824120, 59, 1184),
    Magic(0x0008040200020400, 0x0204106085200814, 59, 1216), Magic(0x0010080400040800, 0x0890110682040111, 59, 1248),
    Magic(0x0020100a000a1000, 0x4000500001040080, 57, 1280), Magic(0x0040221400142200, 0x0000808018020102, 55, 1408),
    Magic(0x0002442800284400, 0x2401010010104008, 55, 1920), Magic(0x0004085000500800, 0x1052008004100080, 57, 2432),
    Magic(0x0008102000201000, 0x0000950104090801, 59, 2560), Magic(0x0010204000402000, 0x1200808000220864, 59, 2592),
    Magic(0x0004020002040800, 0x8188421201082050, 59, 2624), Magic(0x0008040004081000, 0x0035011001a05c00, 59, 2656),
    Magic(0x00100a000a102000, 0x0051480201104400, 57, 2688), Magic(0x0022140014224000, 0x4110040400180210, 55, 2816),
    Magic(0x0044280028440200, 0x0414140400001100, 55, 3328), Magic(0x0008500050080400, 0x4e60008100688040, 57, 3840),
    Magic(0x0010200020100800, 0x000121022ac40200, 59, 3968), Magic(0x0020400040201000, 0x00040042904a0094, 59, 4000),
    Magic(0x0002000204081000, 0x2948010421001040, 59, 4032), Magic(0x0004000408102000, 0x4060821011020200, 59, 4064),
    Magic(0x000a000a10204000, 0x0001004030008202, 57, 4096), Magic(0x0014001422400000, 0x0018020122080401, 57, 4224),
    Magic(0x0028002844020000, 0x0042ca41a2000400, 57, 4352), Magic(0x0050005008040200, 0x01c0212040810100, 57, 4480),
    Magic(0x0020002010080400, 0x1060480f41020041, 59, 4608), Magic(0x0040004020100800, 0x020102020a040040, 59, 4640),
    Magic(0x0000020408102000, 0x1082008220100402, 59, 4672), Magic(0x0000040810204000, 0x0102020202828102, 59, 4704),
    Magic(0x00000a1020400000, 0x4022108848080004, 59, 4736), Magic(0x0000142240000000, 0x0020400508480000, 59, 4768),
    Magic(0x0000284402000000, 0x1010108590440010, 59, 4800), Magic(0x0000500804020000, 0xd812040810a10140, 59, 4832),
    Magic(0x0000201008040200, 0x0191043000a20038, 59, 4864), Magic(0x0000402010080400, 0x8010020200620000, 59, 4896),
    Magic(0x0002040810204000, 0x0001008801084280, 58, 4928), Magic(0x0004081020400000, 0x0000060882015050, 59, 4992),
    Magic(0x000a102040000000, 0x8120000100809001, 59, 5024), Magic(0x0014224000000000, 0x0200100002104420, 59, 5056),
    Magic(0x0028440200000000, 0x10d100200821010a, 59, 5088), Magic(0x0050080402000000, 0x03002c0448100110, 59, 5120),
    Magic(0x0020100804020000, 0x14004005045c0042, 59, 5152), Magic(0x0040201008040200, 0x4002901008831040, 58, 5184)]

function subindex(occupied::Bitboard, magics::Magic)
    (((occupied & magics.mask).val * magics.num.val) >> magics.shift) + 1
end

# obtain the table index using the shift, occupancy mask, and magic number. Adding the offset.
function tableindex(occupied::Bitboard, magics::Magic)
    subindex(occupied, magics) + magics.offset
end

# generate all possible sliding moves in a given move direction, for a given occupancy.
function slidermoves(sqr::Bitboard, occupied::Bitboard, directions::SArray{Tuple{4},Function,1,4})
    mask = EMPTY
    for direction in directions
    	newsqr = sqr
		while !isempty(newsqr)
			newsqr = direction(newsqr)
    		mask ⊻= newsqr
    		!isempty(newsqr & occupied) && break
    	end
    end
    return mask
end

# initiate the attacks table for rook and bishop pieces
function initSlidingTable(table::Vector{Vector{Bitboard}}, magics::SArray{Tuple{64}, Magic, 1, 64}, directions::SArray{Tuple{4},Function,1,4})
	@inbounds for sqr in 1:64
		occupied = EMPTY
		table[sqr] = Vector{Bitboard}(undef, 2^count(magics[sqr].mask))
		for i in 1:(1 << count(magics[sqr].mask))
			idx = subindex(occupied, magics[sqr])
			table[sqr][idx] = slidermoves(Bitboard(sqr), occupied, directions)
			occupied = Bitboard(occupied.val - magics[sqr].mask.val) & magics[sqr].mask
		end
	end
	return table
end

max_stars_repo_name: J-Revell/Astellarn

max_issues_repo_name: J-Revell/Astellarn

max_forks_repo_name: J-Revell/Astellarn


text:
using Flux, DiffEqFlux, DelayDiffEq, Test

## Setup DDE to optimize
function delay_lotka_volterra(du,u,h,p,t)
  x, y = u
  α, β, δ, γ = p
  du[1] = dx = (α - β*y)*h(p,t-0.1)[1]
  du[2] = dy = (δ*x - γ)*y
end
h(p,t) = ones(eltype(p),2)
prob = DDEProblem(delay_lotka_volterra,[1.0,1.0],h,(0.0,10.0),constant_lags=[0.1])
p = param([2.2, 1.0, 2.0, 0.4])
function predict_fd_dde()
  diffeq_fd(p,sol->sol[1,:],101,prob,MethodOfSteps(Tsit5()),saveat=0.0:0.1:10.0)
end
loss_fd_dde() = sum(abs2,x-1 for x in predict_fd_dde())
@test_broken loss_fd_dde()
@test_broken Flux.back!(loss_fd_dde())

function predict_rd_dde()
  diffeq_rd(p,prob,MethodOfSteps(Tsit5()),saveat=0.1)[1,:]
end
loss_rd_dde() = sum(abs2,x-1 for x in predict_rd_dde())
loss_rd_dde()
Flux.back!(loss_rd_dde())


text:
# setup internal paths -
_PATH_TO_SRC = dirname(pathof(@__MODULE__))

# load external packages
using CSV
using DataFrames
using Dates
using Distributions
using Plots
using LinearAlgebra
using BSON
using BSON: @save

# load my codes
include(joinpath(_PATH_TO_SRC, "Types.jl"))
include(joinpath(_PATH_TO_SRC, "Expiration.jl"))
include(joinpath(_PATH_TO_SRC, "Factory.jl"))
include(joinpath(_PATH_TO_SRC, "Lattice.jl"))
include(joinpath(_PATH_TO_SRC, "Greeks.jl"))
include(joinpath(_PATH_TO_SRC, "Files.jl"))
include(joinpath(_PATH_TO_SRC, "Decisions.jl"))
include(joinpath(_PATH_TO_SRC, "Search.jl"))


text:
using NKLandscapes
using Base.Test


text:
module NewFESpacesTests
##
using Gridap.DOFBases: LagrangianDOFBasis

using Gridap, Test

using Gridap.Helpers

using Base.Cartesian

import Gridap: ∇

n = 2
model = CartesianDiscreteModel(partition=(n,))
diritag = [1]
order = 2
fespace = H1ConformingFESpace(Float64,model,order,diritag)

l = 1.0
ufun(x) = x[1]*(x[1]-2.0*l)
ufun_grad(x) = VectorValue(2*x[1]-2.0,)
bfun(x) = -2.0
∇(::typeof(ufun)) = ufun_grad

V = TestFESpace(fespace)
U = TrialFESpace(fespace,ufun)

# Define integration mesh and quadrature
trian = Triangulation(model)
quad = CellQuadrature(trian,degree=order*2)

# Define forms
a(v,u) = inner(∇(v), ∇(u))
b(v) = inner(v,bfun)

# Define Assembler
assem = SparseMatrixAssembler(V,U)

# Define the FEOperator
op = LinearFEOperator(a,b,V,U,assem,trian,quad)

# Define the FESolver
ls = BackslashSolver()
solver = LinearFESolver(ls)

# Solve!
uh = solve(solver,op)

# Define exact solution and error
u = CellField(trian,ufun)
e = u - uh

# Define norms to measure the error
l2(u) = inner(u,u)
h1(u) = a(u,u) + l2(u)

# Compute errors
el2 = sqrt(sum( integrate(l2(e),trian,quad) ))
eh1 = sqrt(sum( integrate(h1(e),trian,quad) ))

@test el2 < 1.e-8
@test eh1 < 1.e-8

# writevtk(trian,"onedres",cellfields=["uh"=>uh,"e"=>e])
##
end # module


text:
# ========================================================
# Methods to solve the minimum time Hamiltonian problem
# ========================================================

abstract type HamiltonianSolver end

include("brute_force.jl")
include("ford_fulkerson.jl")
include("local.jl")


text:
# Version 5.2
@show 1 & 2   # AND
@show 1 | 2   # OR
@show 1 ^ 2   # XOR -- for Julia 6.0 the operator is `⊻`
@show ~1      # NOT
@show 1 >>> 2 # SHIFT RIGHT (LOGICAL)
@show 1 >> 2  # SHIFT RIGHT (ARITMETIC)
@show 1 << 2  # SHIFT LEFT (ARITMETIC/LOGICAL)

A = BitArray([true, true, false, false, true])
@show A ror(A,1) ror(A,2) ror(A,5) # ROTATION RIGHT
@show rol(A,1) rol(A,2) rol(A,5) # ROTATION LEFT

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
module TestSetfield

import PerformanceTestTools
import Setfield
using Documenter: doctest

include("test_setindex.jl")
include("test_examples.jl")
include("test_setmacro.jl")
include("test_core.jl")
include("test_functionlenses.jl")
include("test_staticarrays.jl")
include("test_quicktypes.jl")
PerformanceTestTools.@include("perf.jl")

doctest(Setfield)

end  # module


text:
__precompile__()

"""
Module for dependent bootstrap procedures, by Colin T Bowers

Implemented bootstrap methods: \n
	- IID
	- Stationary
	- Moving Block
	- Circular Block
	- NoOverlapBlock

Implemented block length selection procedures: \n
	- Patton, Politis, and White (2009) Correction to Automatic Block Length Selection For The Dependent Bootstrap

Accepted input dataset types: \n
	- Vector{<:Number}
	- Matrix{<:Number} (where rows are observations and columns are variables)
	- Vector{Vector{<:Number}} (where elements of inner vectors are observations and outer vectors are variables)
	- DataFrame
	- TimeSeries.TimeArray{T,N} (only for N = 1 and N = 2)

Additional input dataset types are easily added. Please open an issue at https://github.com/colintbowers/DependentBootstrap.jl

The module has only a single exported type: \n
	- BootInput  <-- Core input type accepted by all exported functions. Typically constructed via keyword method. See ?BootInput for more detail.

All exported functions exhibit the following keyword signatures: \n
	- exported_func(data, bootinput::BootInput)
	- exported_func(data ; kwargs...)

Most users will be content to use the keyword argument method. In practice, this method wraps a
keyword argument BootInput constructor, which is then input to the exported function BootInput
method. For more detail on accepted keywords, see ?BootInput. All exported functions then use
the input dataset and bootstrap methodology described in BootInput in order to return the
appropriate statistics. A list of exported functions follows: \n
	- optblocklength <-- Estimate the optimal block length for the input dataset
	- dbootinds      <-- Get a vector of bootstrap resampling indices
	- dbootdata      <-- Get a vector of resampled datasets
	- dbootlevel1    <-- Get a vector of level 1 resampled bootstrap statistics
	- dboot          <-- Get the level 2 bootstrap statistic
	- dbootlevel2    <-- Identical to dboot. Included for completeness
	- dbootvar       <-- Wrapper on dboot that sets the level 2 statistic as the variance
	- dbootconf      <-- Wrapper on dboot that sets the level 2 statistic as a confidence interval

I use the phrases level 1 and level 2 statistics in this package in the same manner discussed in
Chapter 1 of Lahiri's textbook Resampling Methods for Dependent Data.

This package has an MIT license. Please see associated LICENSE.md file for more detail.
"""
module DependentBootstrap

using 	StatsBase, Distributions, DataFrames, TimeSeries

import 	Base: 	show

export 	BootInput,
		optblocklength,
		dbootinds,
		dbootdata,
		dbootlevel1,
		dbootlevel2,
		dboot,
		dbootvar,
		dbootconf

const NUM_RESAMPLE = 1000::Int #Default number of bootstrap resamples

#-----------------------------------------------------------------------------------------
# OLD CODE THAT USED Requires.jl TO LOAD DataFrames.jl and TimeSeries.jl
#-----------------------------------------------------------------------------------------
#In order to accommodate different types of datasets, we use the lazy loading
#features offered by the Requires package here, so that the modules in which these
#datasets are defined are not loaded unless they are actually needed. All functions
#that reference these modules are placed here inside __init__() (which is a
#requirement when using Requires package on Julia v0.7+). If you want to add a new
#dataset type to DependentBootstrap, it should be as simple as adding appropriate
#methods here. This module is structured so that nothing else should need to be done.
#The methods that need to be added are:
#    num_obs(data::T)::Int <- Number of observations in dataset
#	 num_var(data::T)::Int <- Number of variables in dataset
#    local_get_var(data::T, i::Int)::Vector <- Get the data associated with the ith variable and output as a Vector
#	 local_get_index(data::T, inds::Vector{Int})::T <- Resample data using resampling indices in inds
#For local_get_index, note that the output type should always match the input type of the dataset.
# function __init__()
# 	@require TimeSeries="9e3dc215-6440-5c97-bce1-76c03772f85e" begin
# 		(num_obs(data::TimeSeries.TimeArray{T,1})::Int) where {T} = size(data, 1)
# 		(num_obs(data::TimeSeries.TimeArray{T,2})::Int) where {T} = size(data, 1)
# 		(num_var(data::TimeSeries.TimeArray{T,1})::Int) where {T} = 1
# 		(num_var(data::TimeSeries.TimeArray{T,2})::Int) where {T} = size(data, 2)
# 		(local_get_var(data::TimeSeries.TimeArray{T,1}, i::Int)::Vector{T}) where {T} = i == 1 ? data.values[:] : error("Invalid index $(i) given data $(typeof(data)) with number of columns $(num_var(data))")
# 		(local_get_var(data::TimeSeries.TimeArray{T,2}, i::Int)::Vector{T}) where {T} = (1 <= i <= num_var(data)) ? data.values[:, i] : error("Invalid index $(i) given data $(typeof(data)) with number of columns $(num_var(data))")
# 		(local_get_index(data::TimeSeries.TimeArray{T,1}, inds::Vector{Int})::TimeSeries.TimeArray{T,1}) where {T} = TimeSeries.TimeArray(TimeSeries.timestamp(data), data.values[inds], TimeSeries.colnames(data) ; unchecked=true)
# 		(local_get_index(data::TimeSeries.TimeArray{T,2}, inds::Vector{Int})::TimeSeries.TimeArray{T,2}) where {T} = TimeSeries.TimeArray(TimeSeries.timestamp(data), data.values[inds, :], TimeSeries.colnames(data) ; unchecked=true)
# 	end
# 	@require DataFrames="a93c6f00-e57d-5684-b7b6-d8193f3e46c0" begin
# 		(num_obs(data::DataFrames.DataFrame)::Int) = size(data, 1)
# 		(num_var(data::DataFrames.DataFrame)::Int) = size(data, 2)
# 		(local_get_var(data::DataFrames.DataFrame, i::Int)) = (1 <= i <= num_var(data)) ? data[:, i] : error("Invalid index $(i) given data $(typeof(data)) with number of columns $(num_var(data))")
# 		(local_get_index(data::DataFrames.DataFrame, inds::Vector{Int})::DataFrames.DataFrame) = data[inds, :]
# 	end
# end
#-----------------------------------------------------------------------------------------

include("types.jl")
include("blocklength.jl")
include("bootinds.jl")
#include("tapered.jl")
include("core.jl")

end # module


text:
"""
Finds roots of 4th order polynomials of the form
```math
ax⁴ + bx³ + cx² + dx + e = 0
```
based on the Fast Quartic and Cubic Solver
https://github.com/NKrvavica/fqs
by Nino Krvavica.
"""
module QuarticSolver

function solvequadratic(a0, b0, c0)
    inva0 = 1/a0
    a, b = b0*inva0, c0*inva0

    a₀ = -0.5*a
    Δ = a₀^2 - b
    sqrtΔ = sqrt(Δ)

    r1 = a₀ - sqrtΔ
    r2 = a₀ + sqrtΔ

    return r1, r2
end

function solvecubic(a0, b0, c0, d0)
    inva0 = 1/a0
    a, b, c = b0*inva0, c0*inva0, d0*inva0

    third = 1/3
    thirda = a*third
    thirda² = thirda^2
    sqrt3 = sqrt(3)

    f = third*b - thirda²
    g = thirda*(2*thirda² - b) + c
    h = 0.25*g^2 + f^3  # discriminant (`Δ` or `d` in some papers)

    if f == g == h == 0  # this check seems a little strong
        r1 = -cbrt(c)
        return r1, r1, r1
    elseif h isa Complex || h <= 0   # casus irreducibilis
        j = sqrt(-f)
        k = acos(-0.5*g/j^3)
        m = cos(third*k)
        n = sqrt3*sin(third*k)
        r1 = 2*j*m - thirda
        r2 = -j*(m + n) - thirda
        r3 = -j*(m - n) - thirda
        return r1, r2, r3
    else
        sqrth = sqrt(h)
        S = cbrt(-0.5*g + sqrth)
        U = cbrt(-0.5*g - sqrth)
        SplusU = S + U
        SminusU = S - U
        tmp = SminusU*sqrt3*0.5im
        r1 = SplusU - thirda
        r2 = -0.5*SplusU - thirda + tmp
        r3 = -0.5*SplusU - thirda - tmp
        return r1, r2, r3
    end
end

function solvequartic(a0, b0, c0, d0, e0)
    inva0 = 1/a0
    a, b, c, d = b0*inva0, c0*inva0, d0*inva0, e0*inva0

    a₀ = 0.25*a
    a₀² = a₀^2

    # Subsidiary cubic equation
    p = 3*a₀² - 0.5*b
    q = a*a₀² - b*a₀ + 0.5*c
    r = 3*a₀²^2 - b*a₀² + c*a₀ - d

    # One root of cubic
    z0, _, _ = solvecubic(1, p, r, p*r-0.5*q^2)
    z0 = complex(z0)

    s = sqrt(2*p + 2*z0)
    s == 0 ? t = z0^2 + r : t = -q/s

    r0, r1 = solvequadratic(1, s, z0 + t)
    r2, r3 = solvequadratic(1, -s, z0 - t)

    return r0-a₀, r1-a₀, r2-a₀, r3-a₀
end

end # module

max_stars_repo_name: EP-Guy/QuarticSolver

max_issues_repo_name: EP-Guy/QuarticSolver

max_forks_repo_name: EP-Guy/QuarticSolver


text:
#Adapted from IndirectArray

struct Mat{T <: dtypes} <: AbstractArray{T,3}
    mat
    data_raw
    data

    @inline function Mat{T}(mat, data_raw::AbstractArray{T,3}) where {T <: dtypes}
        data = reinterpret(T, data_raw)
        new{T}(mat, data_raw, data)
    end

    @inline function Mat(data_raw::AbstractArray{T, 3}) where {T <: dtypes}
        data = reinterpret(T, data_raw)
        mat = nothing
        new{T}(mat, data_raw, data)
    end
end

function Base.deepcopy_internal(x::Mat{T}, y::IdDict) where {T}
    if haskey(y, x)
        return y[x]
    end
    ret = Base.copy(x)
    y[x] = ret
    return ret
end

Base.size(A::Mat) = size(A.data)
Base.axes(A::Mat) = axes(A.data)
Base.IndexStyle(::Type{Mat{T}}) where {T} = IndexCartesian()

Base.strides(A::Mat{T}) where {T} = strides(A.data)
Base.copy(A::Mat{T}) where {T} = Mat(copy(A.data_raw))
Base.pointer(A::Mat) = Base.pointer(A.data)

Base.unsafe_convert(::Type{Ptr{T}}, A::Mat{S}) where {T, S} = Base.unsafe_convert(Ptr{T}, A.data)

@inline function Base.getindex(A::Mat{T}, I::Vararg{Int,3}) where {T}
    @boundscheck checkbounds(A.data, I...)
    @inbounds ret = A.data[I...]
    ret
end

@inline function Base.setindex!(A::Mat, x, I::Vararg{Int,3})
    @boundscheck checkbounds(A.data, I...)
    A.data[I...] = x
    return A
end


text:
# This file was generated, do not modify it. # hide
coerce!(X, Count => Continuous)

xgb  = XGBR()
xgbm = machine(xgb, X, y)
fit!(xgbm, rows=train)

rms(y[test], MLJ.predict(xgbm, rows=test))

max_issues_repo_name: alan-turing-institute/MLJTutorials

max_forks_repo_name: alan-turing-institute/MLJTutorials


text:
function get_topk(probs, k = 5)
  T = eltype(probs)
  prob = Array{T, 1}()
  idx = Array{Int, 1}()
  while(k!=0)
    push!(idx, indmax(probs))
    push!(prob, probs[idx[end]])
    probs[idx[end]] = 0.0
    k -= 1
  end
  (prob, idx)
end

function one_hot_encode(preds, idx)
  one_hot = zeros(eltype(preds.data), size(preds)[1], 1)
  one_hot[idx ,1] = 255.0
  one_hot |> gpu
end

function save_gradient_images(gradient, file_name)
  gradient = gradient - minimum(gradient)
  gradient /= maximum(gradient)
  gradient = permutedims(squeeze(gradient, 4), (3, 2, 1))
  img = colorview(RGB{eltype(gradient)}, gradient)
  display(img)
  save(file_name, img)
end

function save_grayscale_gradient(gradient, file_name)
  gradient = gradient - minimum(gradient)
  gradient /= maximum(gradient)
  gradient = permutedims(squeeze(gradient, 4), (3, 2, 1))
  img = colorview(RGB{eltype(gradient)}, gradient)
  img = Gray.(img)
  display(img)
  save(file_name, img)
end

function positive_negative_saliency(gradient)
  pos_saliency = max.(zero(gradient), gradient) / maximum(gradient)
  neg_saliency = max.(zero(gradient), -gradient) / maximum(-gradient)
  (pos_saliency, neg_saliency)
end

im2arr_rgb(img) = permutedims(float.(channelview(imresize(img, (224, 224)))), (3, 2, 1))

im_mean = reshape([0.485, 0.456, 0.406], 1, 1, 3)

# The below code is a directly copied from https://github.com/avik-pal/DeepDream.jl
# There might be some code redundancy

function image_to_arr(img; preprocess = true)
  local x = img
  x = Float32.(channelview(img))
  x = permutedims(x, [3,2,1])
  if(preprocess)
    x = x .- im_mean
  end
  x = reshape(x, size(x,1), size(x,2), size(x,3), 1) * 255 |> gpu
end

function load_image(path, resize = false; size_save = true)
  img = load(path)
  if(size_save)
    global original_size = size(img)
  end
  if(resize)
    image_to_arr(imresize(img, (224, 224)))
  else
    image_to_arr(img)
  end
end

function generate_image(x, resize_original = false)
  x = reshape(x, size(x)[1:3]...)/255 |> cpu
  x = x .+ im_mean
  x = permutedims(x, [3,2,1])
  x -= minimum(x)
  x /= maximum(x)
  if resize_original
    imresize(colorview(RGB, x), original_size)
  else
    colorview(RGB, x)
  end
end

save_image(path, x) = save(path, generate_image(x, true))


text:
include("custom_log.jl")

include("symbol_func.jl")

include("norm_int.jl")

include("is_positive.jl")

include("find_roots.jl")

include("reorder_lists.jl")

include("romanize.jl")

include("pluck.jl")

include("is_inside.jl")

include("is_true.jl")

include("progress.jl")

include("filter_approx.jl")

include("safe_symbol.jl")

include("safe_get.jl")

include("is_nan.jl")

include("is_nothing.jl")


text:
module Day14

import ..data_dir # from parent module
input = read(joinpath(data_dir, "day14"), String)

export part1, part2

struct Insertion
    pattern::AbstractString
    insert::AbstractString
end

function Insertion(rule::AbstractString)
    Insertion(split(rule, " -> ")...)
end

function template_and_insertions(input::AbstractString)
    lines = split(input, "\n")
    template = lines[1]
    insertions = Insertion.(lines[3:end])
    return template, insertions
end

function insert_after(s::AbstractString, to_insert::AbstractString, index::Int64)
    s[1:index] * to_insert * s[index+1:end]
end

function substitute(template::AbstractString, insertions::Vector{Insertion})
    characters_to_insert = Vector{Tuple{Int64,AbstractString}}()
    for i = 1:length(template)-1
        current_pair = template[i] * template[i+1]
        for insertion in insertions
            if current_pair == insertion.pattern
                push!(characters_to_insert, (i, insertion.insert))
                break
            end
        end
    end

    for i = 1:length(characters_to_insert)
        index, to_insert = characters_to_insert[i]
        offset = i - 1
        index_with_offset = index + offset
        template = insert_after(template, to_insert, index_with_offset)
    end

    return template
end

function most_common_minus_least_common(template::AbstractString)
    letters = split(template, "")
    counts = [count(x -> x == letter, letters) for letter in unique(letters)]
    maximum(counts) - minimum(counts)
end

function part1(input = input)
    template, insertions = template_and_insertions(input)
    for i = 1:10
        template = substitute(template, insertions)
    end
    most_common_minus_least_common(template)
end

mutable struct CompressedTemplate
    pairs::Dict{Tuple{Char,Char},Int64}
    letter_counts::Dict{Char,Int64}
end

function CompressedTemplate(template::AbstractString)
    pairs = Dict{Tuple{Char,Char},Int64}()
    for i = 1:length(template)-1
        current_pair = (template[i], template[i+1])
        pairs[current_pair] = get(pairs, current_pair, 0) + 1
    end

    letters = split(template, "")
    letter_counts = Dict(only(letter) => count(x -> x == letter, letters) for letter in unique(letters))

    CompressedTemplate(pairs, letter_counts)
end

function pair(insertion::Insertion)
    left, right = split(insertion.pattern, "")
    (only(left), only(right))
end

function new_pairs(insertion::Insertion)
    left, right = split(insertion.pattern, "")
    # only converts to Char
    left, right, insert = only(left), only(right), only(insertion.insert)
    return (left, insert), (insert, right)
end

function substitute(template::CompressedTemplate, insertions::Vector{Insertion})
    new_template = deepcopy(template)
    for insertion in insertions
        current_pattern = pair(insertion)
        occurrences = get(template.pairs, current_pattern, 0)

        pair_left, pair_right = new_pairs(insertion)
        new_template.pairs[pair_left] = get(new_template.pairs, pair_left, 0) + occurrences
        new_template.pairs[pair_right] = get(new_template.pairs, pair_right, 0) + occurrences
        new_template.pairs[current_pattern] = get(new_template.pairs, current_pattern, 0) - occurrences

        new_letter = only(insertion.insert)
        new_template.letter_counts[new_letter] = get(new_template.letter_counts, new_letter, 0) + occurrences
    end
    new_template
end

function most_common_minus_least_common(template::CompressedTemplate)
    counts = values(template.letter_counts)
    maximum(counts) - minimum(counts)
end

function part2(input = input)
    template, insertions = template_and_insertions(input)
    template = CompressedTemplate(template)
    for i = 1:40
        template = substitute(template, insertions)
    end
    most_common_minus_least_common(template)
end

end #module


text:
module PositionVelocityTime

    using DocStringExtensions, Geodesy, GNSSDecoder, GNSSSignals, LinearAlgebra, Parameters
    using Unitful: s, Hz

    export  calc_PVT,
            is_sat_healthy_and_decodable,
            PVTSolution,
            sat_position_ECI_2_ECEF,
            sat_position_ECEF,
            user_position,
            SatelliteState

    """
    Struct of decoder, code- and carrierphase of satellite vehicle
    """
    @with_kw struct SatelliteState{CP <: Real}
        decoder_state::GNSSDecoderState
        code_phase::CP
        carrier_phase::CP = 0
    end

    """
    Calculates ECEF position of user

    $SIGNATURES
    ´sat_state´: satellite state, combining decoded data, code- and carrierphase

    This function calculates the position of the user in ECEF coordinates
    The implementation follows IS-GPS-200K Table 20-IV.
    """
    function calc_PVT(
        satellite_states::AbstractVector{SatelliteState{Float64}}
        )

        usable_satellite_states = filter(x -> is_sat_healthy_and_decodable(x.decoder_state), satellite_states)

        length(usable_satellite_states) >= 4 || throw(SmallData("Not enough usable SV Data"))
        sv_positions = map(x -> sat_position_ECEF(x), usable_satellite_states)
        pseudoranges = pseudo_ranges(usable_satellite_states)

        userpos = user_position(sv_positions, pseudoranges)

    end

    include("pseudo_range.jl")
    include("user_position.jl")
    include("sat_position.jl")
    include("sv_time.jl")
    include("errors.jl")
end #end module