text:
using LinearAlgebra
using KrylovKit

export AbstractLattice, SquareLattice
abstract type AbstractLattice end
struct SquareLattice <: AbstractLattice end

export VUMPSRuntime, SquareVUMPSRuntime

# NOTE: should be renamed to more explicit names
"""
    VUMPSRuntime{LT}

a struct to hold the tensors during the `vumps` algorithm, containing
- `d × d × d × d` `M` tensor
- `D × d × D` `AL` tensor
- `D × D`     `C` tensor
- `D × d × D` `AR` tensor
- `D × d × D` `FL` tensor
- `D × d × D` `FR` tensor
and `LT` is a AbstractLattice to define the lattice type.
"""
struct VUMPSRuntime{LT,T,N,AT<:AbstractArray{T,N},ET,CT}
    M::AT
    AL::ET
    C::CT
    AR::ET
    FL::ET
    FR::ET
    function VUMPSRuntime{LT}(M::AT, AL::AbstractArray{T}, C::AbstractArray{T}, AR::AbstractArray{T},
        FL::AbstractArray{T}, FR::AbstractArray{T}) where {LT<:AbstractLattice,T,N,AT<:AbstractArray{T,N}}
        new{LT,T,N,AT,typeof(AL),typeof(C)}(M,AL,C,AR,FL,FR)
    end
end

const SquareVUMPSRuntime{T,AT} = VUMPSRuntime{SquareLattice,T,4,AT}
SquareVUMPSRuntime(M::AT,AL,C,AR,FL,FR) where {T,AT<:AbstractArray{T, 4}} = VUMPSRuntime{SquareLattice}(M,AL,C,AR,FL,FR)

getD(rt::VUMPSRuntime) = size(rt.AL, 1)
getd(rt::VUMPSRuntime) = size(rt.M, 1)

@doc raw"
    SquareVUMPSRuntime(M::AbstractArray{T,4}, env::Val, χ::Int)

create a `SquareVUMPSRuntime` with M-tensor `M`. The AL,C,AR,FL,FR
tensors are initialized according to `env`. If `env = Val(:random)`,
the A is initialized as a random D×d×D tensor,and AL,C,AR are the corresponding
canonical form. FL,FR is the left and right environment:
```
┌── AL─       ┌──        ─ AR──┐         ──┐
│   │         │            │   │           │
FL─ M ─  = λL FL─        ─ M ──FR   = λR ──FR
│   │         │            │   │           │
┕── AL─       ┕──        ─ AR──┘         ──┘
```

# example

```jldoctest; setup = :(using ADVUMPS)
julia> rt = SquareVUMPSRuntime(randn(2,2,2,2), Val(:random), 4);

julia> size(rt.AL) == (4,2,4)
true

julia> size(rt.C) == (4,4)
true
```
"
function SquareVUMPSRuntime(M::AbstractArray{T,4}, env, D::Int; verbose = false) where T
    return SquareVUMPSRuntime(M, _initializect_square(M, env, D; verbose = verbose)...)
end

function _initializect_square(M::AbstractArray{T,4}, env::Val{:random}, D::Int; verbose = false) where T
    d = size(M,1)
    A = rand(T,D,d,D)
    AL, = leftorth(A)
    C, AR = rightorth(AL)
    _, FL = leftenv(AL, M)
    _, FR = rightenv(AR, M)
    verbose && print("random initial vumps environment-> ")
    AL,C,AR,FL,FR
end

function _initializect_square(M::AbstractArray{T,4}, chkp_file::String, D::Int; verbose = false) where T
    env = load(chkp_file)["env"]
    verbose && print("vumps environment load from $(chkp_file) -> ")
    env.AL,env.C,env.AR,env.FL,env.FR
end

function vumps(rt::VUMPSRuntime; tol::Real, maxiter::Int, verbose = false)
    # initialize
    olderror = Inf

    stopfun = StopFunction(olderror, -1, tol, maxiter)
    rt, err = fixedpoint(res->vumpstep(res...), (rt, olderror), stopfun)
    verbose && println("vumps done@step: $(stopfun.counter), error=$(err)")
    return rt
end

function vumpstep(rt::VUMPSRuntime,err)
    # global backratio = 1.0
    # Zygote.@ignore print(round(-log(10,backratio)),' ')
    M,AL,C,AR,FL,FR= rt.M,rt.AL,rt.C,rt.AR,rt.FL,rt.FR
    AC = Zygote.@ignore ein"asc,cb -> asb"(AL,C)
    _, AC = ACenv(AC, FL, M, FR)
    _, C = Cenv(C, FL, FR)
    AL, AR, _, _ = ACCtoALAR(AC, C)
    _, FL = leftenv(AL, M, FL)
    _, FR = rightenv(AR, M, FR)

    ##### avoid gradient explosion for too many iterations #####
    # M = backratio .* M + Zygote.@ignore (1-backratio) .* M
    # AL = backratio .* AL + Zygote.@ignore (1-backratio) .* AL
    # C = backratio .* C +  Zygote.@ignore (1-backratio) .* C
    # AR = backratio .* AR + Zygote.@ignore (1-backratio) .* AR
    # FL = backratio .* FL + Zygote.@ignore (1-backratio) .* FL
    # FR = backratio .* FR + Zygote.@ignore (1-backratio) .* FR

    err = error(AL,C,FL,M,FR)
    return SquareVUMPSRuntime(M, AL, C, AR, FL, FR), err
end

safesign(x::Number) = iszero(x) ? one(x) : sign(x)
"""
    qrpos(A)

Returns a QR decomposition, i.e. an isometric `Q` and upper triangular `R` matrix, where `R`
is guaranteed to have positive diagonal elements.
"""
qrpos(A) = qrpos!(copy(A))
function qrpos!(A)
    F = qr!(A)
    Q = Matrix(F.Q)
    R = F.R
    phases = safesign.(diag(R))
    rmul!(Q, Diagonal(phases))
    lmul!(Diagonal(conj!(phases)), R)
    return Q, R
end

"""
    lqpos(A)

Returns a LQ decomposition, i.e. a lower triangular `L` and isometric `Q` matrix, where `L`
is guaranteed to have positive diagonal elements.
"""
lqpos(A) = lqpos!(copy(A))
function lqpos!(A)
    F = qr!(Matrix(A'))
    Q = Matrix(Matrix(F.Q)')
    L = Matrix(F.R')
    phases = safesign.(diag(L))
    lmul!(Diagonal(phases), Q)
    rmul!(L, Diagonal(conj!(phases)))
    return L, Q
end

"""
    leftorth(A, [C]; kwargs...)

Given an MPS tensor `A`, return a left-canonical MPS tensor `AL`, a gauge transform `C` and
a scalar factor `λ` such that ``λ AL^s C = C A^s``, where an initial guess for `C` can be
provided.
```
    ┌─AL─      ┌──
    │ │     =  │
    ┕─AL─      ┕──
```

"""
function leftorth(A, C = Matrix{eltype(A)}(I, size(A,1), size(A,1)); tol = 1e-12, maxiter = 100, kwargs...)
    λ2s, ρs, info = eigsolve(C'*C, 1, :LM; ishermitian = false, tol = tol, maxiter = 1, kwargs...) do ρ
        ρE = ein"cd,dsb,csa -> ab"(ρ, A, conj(A))
        return ρE
    end
    ρ = ρs[1] + ρs[1]'
    ρ ./= tr(ρ)
    # C = cholesky!(ρ).U
    # If ρ is not exactly positive definite, cholesky will fail
    F = svd!(ρ)
    C = lmul!(Diagonal(sqrt.(F.S)), F.Vt)
    _, C = qrpos!(C)

    D, d, = size(A)
    Q, R = qrpos!(reshape(C*reshape(A, D, d*D), D*d, D))
    AL = reshape(Q, D, d, D)
    λ = norm(R)
    rmul!(R, 1/λ)
    numiter = 1
    while norm(C-R) > tol && numiter < maxiter
        # C = R
        λs, Cs, info = eigsolve(R, 1, :LM; ishermitian = false, tol = tol, maxiter = 1, kwargs...) do X
            Y = ein"cd,dsb,csa -> ab"(X,A,conj(AL))
            return Y
        end
        _, C = qrpos!(Cs[1])
        # The previous lines can speed up the process when C is still very far from the correct
        # gauge transform, it finds an improved value of C by finding the fixed point of a
        # 'mixed' transfer matrix composed of `A` and `AL`, even though `AL` is also still not
        # entirely correct. Therefore, we restrict the number of iterations to be 1 and don't
        # check for convergence
        Q, R = qrpos!(reshape(C*reshape(A, D, d*D), D*d, D))
        AL = reshape(Q, D, d, D)
        λ = norm(R)
        rmul!(R, 1/λ)
        numiter += 1
    end
    C = R
    return real(AL), real(C), λ
end

"""
    rightorth(A, [C]; kwargs...)

Given an MPS tensor `A`, return a gauge transform C, a right-canonical MPS tensor `AR`, and
a scalar factor `λ` such that `λ C AR^s = A^s C`, where an initial guess for `C` can be
provided.
````
    ─ AR─┐     ──┐
      │  │  =    │
    ─ AR─┘     ──┘
````
"""
function rightorth(A, C = Matrix{eltype(A)}(I, size(A,1), size(A,1)); tol = 1e-12, maxiter = 100, kwargs...)
    AL, C, λ = leftorth(permutedims(A,(3,2,1)), permutedims(C,(2,1)); tol = tol,maxiter = maxiter, kwargs...)
    return permutedims(C,(2,1)), permutedims(AL,(3,2,1)), λ
end

"""
    leftenv(A, M, FL; kwargs)

Compute the left environment tensor for MPS `AL` and MPO `M`, by finding the left fixed point
of `AL - M - conj(AL)` contracted along the physical dimension.
```
┌── AL─       ┌──
│   │         │
FL─ M ─  = λL FL─
│   │         │
┕── AL─       ┕──
```
"""
function leftenv(AL, M, FL = rand(eltype(AL), size(AL,1), size(M,1), size(AL,1)); kwargs...)
    λs, FLs, info = eigsolve(FL -> ein"γcη,ηpβ,csap,γsα -> αaβ"(FL,AL,M,conj(AL)),FL, 1, :LM; ishermitian = false, kwargs...)
    # if length(λs) > 1 && norm(real(λs[1]) - real(λs[2])) < 1e-12
    #     @show λs
    #     if real(λs[1]) > 0
    #         return real(λs[1]), real(FLs[1])
    #     else
    #         return real(λs[2]), real(FLs[2])
    #     end
    # end
    return real(λs[1]), real(FLs[1])
end

"""
    rightenv(A, M, FR; kwargs...)

Compute the right environment tensor for MPS `AR` and MPO `M`, by finding the right fixed point
of `AR - M - conj(AR)`` contracted along the physical dimension.
```
 ─ AR──┐         ──┐
   │   │           │
 ─ M ──FR   = λR ──FR
   │   │           │
 ─ AR──┘         ──┘
```
"""
function rightenv(AR, M, FR = randn(eltype(AR), size(AR,1), size(M,3), size(AR,1)); kwargs...)
    λs, FRs, info = eigsolve(FR -> ein"αpγ,γcη,ascp,βsη -> αaβ"(AR,FR,M,conj(AR)), FR, 1, :LM; ishermitian = false, kwargs...)
    # if length(λs) > 1 && norm(real(λs[1]) - real(λs[2])) < 1e-12
    #     @show λs
    #     if real(λs[1]) > 0
    #         return real(λs[1]), real(FRs[1])
    #     else
    #         return real(λs[2]), real(FRs[2])
    #     end
    # end
    return real(λs[1]), real(FRs[1])
end

"""
    λ, FL4 = bigleftenv(AL, M, FL4 = rand(eltype(AL), size(AL,1), size(M,1), size(M,1), size(AL,1)); kwargs...)

Compute the left environment tensor for MPS `AL` and MPO `M`, by finding the left fixed point
of `AL - M - M - conj(AL)` contracted along the physical dimension.
```
┌── AL─       ┌──
│   │         │
│── M ─       │──
FL4 │    = λL FL4
│── M ─       │──
│   │         │
┕── AL─       ┕──
```
"""
function bigleftenv(AL, M, FL4 = rand(eltype(AL), size(AL,1), size(M,1), size(M,1), size(AL,1)); kwargs...)
    λFL4s, FL4s, info = eigsolve(FL4 -> ein"dcba,def,ckge,bjhk,aji -> fghi"(FL4,AL,M,M,conj(AL)), FL4, 1, :LM; ishermitian = false, kwargs...)
    # @show λFL4s
    return real(λFL4s[1]), real(FL4s[1])
end

"""
    λ, FR4 = bigrightenv(AR, M, FR4 = randn(eltype(AR), size(AR,1), size(M,3), size(M,3), size(AR,1)); kwargs...)

Compute the right environment tensor for MPS `AR` and MPO `M`, by finding the right fixed point
of `AR - M - conj(AR)`` contracted along the physical dimension.
```
 ─ AR──┐         ──┐
   │   │           │
 ─ M ──│         ──│
   │  FR4   = λR  FR4
 ─ M ──│         ──│
   │   │           │
 ─ AR──┘         ──┘
```
"""
function bigrightenv(AR, M, FR4 = randn(eltype(AR), size(AR,1), size(M,3), size(M,3), size(AR,1)); kwargs...)
    λFR4s, FR4s, info = eigsolve(FR4 -> ein"fghi,def,ckge,bjhk,aji -> dcba"(FR4,AR,M,M,conj(AR)), FR4, 1, :LM; ishermitian = false, kwargs...)
    # @show λFR4s
    return real(λFR4s[1]), real(FR4s[1])
end

"""
Compute the up environment tensor for MPS `FL`,`FR` and MPO `M`, by finding the up fixed point
    of `FL - M - FR` contracted along the physical dimension.
````
┌── AC──┐
│   │   │           ┌── AC──┐
FL─ M ──FR  =  λAC  │   │   │
│   │   │

````
"""
function ACenv(AC, FL, M, FR;kwargs...)
    λs, ACs, _ = eigsolve(AC -> ein"αaγ,γpη,asbp,ηbβ -> αsβ"(FL,AC,M,FR), AC, 1, :LM; ishermitian = false, kwargs...)
    # if length(λs) > 1 && norm(real(λs[1]) - real(λs[2])) < 1e-12
    #     @show λs
    #     if real(λs[1]) > 0
    #         return real(λs[1]), real(ACs[1])
    #     else
    #         return real(λs[2]), real(ACs[2])
    #     end
    # end
    return real(λs[1]), real(ACs[1])
end

"""
Compute the up environment tensor for MPS `FL` and `FR`, by finding the up fixed point
    of `FL - FR` contracted along the physical dimension.
````
┌──C──┐
│     │          ┌──C──┐
FL─── FR  =  λC  │     │
│     │

````
"""
function Cenv(C, FL, FR;kwargs...)
    λs, Cs, _ = eigsolve(C -> ein"αaγ,γη,ηaβ -> αβ"(FL,C,FR), C, 1, :LM; ishermitian = false, kwargs...)
    # if length(λs) > 1 && norm(real(λs[1]) - real(λs[2])) < 1e-12
    #     @show λs
    #     if real(λs[1]) > 0
    #         return real(λs[1]), real(Cs[1])
    #     else
    #         return real(λs[2]), real(Cs[2])
    #     end
    # end
    return real(λs[1]), real(Cs[1])
end

"""
    AL, AR, errL, errR = ACCtoALAR(AC, C)

QR factorization to get `AL` and `AR` from `AC` and `C`

````
──AL──C──  =  ──AC──  = ──C──AR──
  │             │            │
````
"""
function ACCtoALAR(AC, C)
    D, d, = size(AC)

    QAC, RAC = qrpos(reshape(AC,(D*d, D)))
    QC, RC = qrpos(C)
    AL = reshape(QAC*QC', (D, d, D))
    errL = norm(RAC-RC)
    # @show errL

    LAC, QAC = lqpos(reshape(AC,(D, d*D)))
    LC, QC = lqpos(C)
    AR = reshape(QC'*QAC, (D, d, D))
    errR = norm(LAC-LC)
    # @show errR
    return AL, AR, errL, errR
end

"""
    err = error(AL,C,FL,M,FR)

Compute the error through all environment `AL,C,FL,M,FR`

````
        ┌── AC──┐
        │   │   │           ┌── AC──┐
MAC1 =  FL─ M ──FR  =  λAC  │   │   │
        │   │   │

        ┌── AC──┐
        │   │   │           ┌──C──┐
MAC2 =  FL─ M ──FR  =  λAC  │     │
        │   │   │
        ┕── AL─

── MAC1 ──    ≈    ── AL ── MAC2 ──
    │                 │
````
"""
function error(AL,C,FL,M,FR)
    AC = ein"asc,cb -> asb"(AL,C)
    MAC = ein"αaγ,γpη,asbp,ηbβ -> αsβ"(FL,AC,M,FR)
    MAC -= ein"asd,cpd,cpb -> asb"(AL,conj(AL),MAC)
    err = norm(MAC)
end


text:
# This is an example of two synaptically-coupled single-compartment neurons with Hodgkin-Huxley dynamics.
# It's written for illustrative purposes and uses the high-level interface of ModelingToolkit.jl as it already exists today.

using DifferentialEquations, ModelingToolkit, OrdinaryDiffEq, Unitful

import Unitful: Voltage

# Channel gating functions
function αm(V::Voltage)
    V = ustrip(Float64, u"mV", V)
    return (0.01*(V + 55))/(1 - exp(-(V + 55)/10))
end

function βn(V::Voltage)
    V = ustrip(Float64, u"mV", V)
    return .125 * exp(-(V + 65)/80)
end

function αm(V::Voltage)
    V = ustrip(Float64, u"mV", V)
    (0.1*(V + 40))/(1 - exp(-(V + 40)/10))
end

function βm(V::Voltage)
    V = ustrip(Float64, u"mV", V)
    4*exp(-(V + 65)/18)
end

function αh(V)
    V = ustrip(Float64, u"mV", V)
    0.07*exp(-(V+65)/20)
end

function βh(V)
    V = ustrip(Float64, u"mV", V)
    1/(1 + exp(-(V + 35)/10))
end

# Steady-state channel gating
m∞(V) = αm(V)/(αm(V) + βm(V))
h∞(V) = αh(V)/(αh(V) + βh(V))
n∞(V) = αn(V)/(αn(V) + βn(V))

@parameters t
params = @parameters Ie gl gna gk gs El Ena Ek Cm
vars = @variables V(t) m(t) h(t) n(t)
D = Differential(t)
# HH equations
eqs = [D(V) ~ (Ie - gl*(V - El) - gna*m^3*h*(V - Ena) - gk*n^4*(V - Ek))/Cm,
       D(m) ~ αm(V) * (1 - m) - βm(V)*m,
       D(h) ~ αh(V) * (1 - h) - βh(V)*h,
       D(n) ~ αn(V) * (1 - n) - βn(V)*n]

sys = ODESystem(eqs, t, [vars...], [])

Vinit = -65.0u"mV"
radius = 0.0025u"cm" # given in cm
area = 4*pi*(radius)^2 # area in cm²

# Set applied current to zero initially
u0 = [V => Vinit
      m => m∞(Vinit)
      h => h∞(Vinit)
      n => n∞(Vinit)]

p = [Ie => 0.8e-3u"mA" # tonic excitability to force first neuron to fire repetitively
     gl  => 0.3u"mS/cm^2" * area
     gna => 120.0u"mS/cm^2" * area
     gk  => 36.0u"mS/cm^2" * area
     gs  => 0.025u"mS/cm^2" * area
     El  => -54.4u"mV"
     Ena => 50.0u"mV"
     Ek  => -77.0u"mV"
     Cm  => 1.0u"µF/cm^2" * area]

tspan = (0.0u"ms", 250.0u"ms")

prob = ODEProblem(sys, u0, tspan, p)
sol = solve(prob, KenCarp47())
plot(sol, vars=(0,1))
plot!(sol, vars=(0,2))

ext: jl

max_stars_repo_name: wsphillips/FutureNeuralSimulation

max_issues_repo_name: wsphillips/FutureNeuralSimulation

max_forks_repo_name: wsphillips/FutureNeuralSimulation


text:
using ProblemParser
using Utils

using Chain
using Lazy
using Underscores

const file = find_input(@__FILE__)
const input = parse(FirstRest(Blocks(),
        Split(",", Convert(Int8)),
        Blocks(Rectangular(Lines(Split(Convert(Int8)))))
    ), file |> slurp)

function cross(board, num)
    for i in eachindex(board)
        val = board[i]
        if val == num
            board[i] = -1
            return
        end
    end
end

function winner(board)
    along(rowcol) = any(curry(all, ==(-1)), rowcol(board))
    along(eachcol) || along(eachrow)
end

score(board, num) = sum(filter(>=(0), board)) * num

function part1(input)
    boards = deepcopy(input[2])
    for draw in input[1], board in boards
        cross(board, draw)
        winner(board) && return score(board, draw)
    end
end
assertequal(part1(input), 8136)

function part2(input)
    boards::Vector{Matrix{Int64}} = deepcopy(input[2])
    for draw in input[1]
        next::typeof(boards) = []
        for board in boards
            cross(board, draw)
            if !winner(board)
                push!(next, board)
            elseif length(boards) == 1
                return score(board, draw)
            end
        end
        boards = next
    end
end
assertequal(part2(input), 12738)

max_stars_repo_name: CmdQ/AoC


text:
x = ccall((:add, "libccall"), Int32, (Int32, Int32), 1, 2)
println(1, " + ", 2, " = ", x)

y = [1.0 2.0 3.0; 4.0 5.0 6.0; 7.0 8.0 9.0]
println("Before increment: ", y)
ccall((:incr_array, "libccall"), Cvoid, (Array{Float64},), y)
println("After increment: ", y)

max_stars_repo_name: gobuki/jlrs

max_issues_repo_name: gobuki/jlrs

max_forks_repo_name: gobuki/jlrs


text:
# Practical Machine learning
# Decision Trees example
# Chapter 5

using Base.Test
using DecisionTree

n,m = 10^3, 5 ;
features = rand(n,m);
weights = rand(-1:1,m);
labels = _int(features * weights);

println("\n##### nfoldCV Classification Tree #####")
accuracy = nfoldCV_tree(labels, features, 0.9, 3)
@test mean(accuracy) > 0.7

println("\n##### nfoldCV Classification Forest #####")
accuracy = nfoldCV_forest(labels, features, 2, 10, 3)
@test mean(accuracy) > 0.7

println("\n##### nfoldCV Adaboosted Stumps #####")
nfoldCV_stumps(labels, features, 7, 3)


text:
import Agtor.@def

"""Represents generic farm infrastructure."""
abstract type Infrastructure <: AgComponent end

@def infrastructure_fields begin
    name::String

    # capital per ha or total implementation cost.
    # up to each implementing component to correctly calculate
    # maintenance, etc.
    capital_cost::Union{AgParameter, Float64}

    # num years maintenance occurs, and
    # assumed proportion of capital cost
    minor_maintenance_schedule::Union{Int64, Float64, AgParameter}
    major_maintenance_schedule::Union{Int64, Float64, AgParameter}
    minor_maintenance_rate::Union{Float64, AgParameter}
    major_maintenance_rate::Union{Float64, AgParameter}
end

"""Calculate maintenance costs.

Warning: This can be on a per hectare basis or given as a total.
"""
function maintenance_cost(infra::Infrastructure, year_step::Int64)::Float64
    minor=infra.minor_maintenance_schedule
    major=infra.major_maintenance_schedule
    maintenance_cost::Float64 = 0.0
    if year_step % major == 0
        maintenance_cost = infra.capital_cost * infra.major_maintenance_rate
    elseif year_step % minor == 0
        maintenance_cost = infra.capital_cost * infra.minor_maintenance_rate
    end

    return maintenance_cost
end

max_stars_repo_name: ConnectedSystems/agtor-jl

max_issues_repo_name: ConnectedSystems/agtor-jl

max_forks_repo_name: ConnectedSystems/agtor-jl


text:
using Clang
using Clang.LibClang.Clang_jll

const LBT_ROOT = abspath(dirname(dirname(@__DIR__)))

# Steal default include path from `cc`
function get_default_include_path()
    io = IOBuffer()
    run(pipeline(`cc -E -Wp,-v -`, stdout=devnull, stderr=io, stdin=devnull))
    lines = split(String(take!(io)), "\n")
    start_idx = findfirst(l -> startswith(l, "#include <") && endswith(l, "search starts here:"), lines)
    stop_idx = findfirst(l -> startswith(l, "End of search list"), lines[start_idx+1:end])
    return strip.(lines[start_idx+1:start_idx+stop_idx-1])
end

function extract_symbol_names(include_dir::String = joinpath(LBT_ROOT, "include", "LP64", "x86_64-linux-gnu"))
    # Collect all headers in that directory (nonrecursive because I'm lazy)
    headers = [joinpath(include_dir, header) for header in readdir(include_dir) if endswith(header, ".h")]

    # Build `-I` arguments to all the places we'll need
    clang_args = String[
        # Include standard locations
        ["-I$(inc)" for inc in get_default_include_path()]...,

        # Include our own include directory
        "-I$(include_dir)",
    ]

    # Parse the headers into a `Clang.jl` context
    @info("Extracting symbols from $(length(headers)) headers in $(include_dir)")
    ctx = DefaultContext()
    parse_headers!(ctx, headers, args=clang_args, includes=headers)

    # Collect all symbol names defined within these translation units
    symbol_names = String[]
    for trans_unit in ctx.trans_units
        root_cursor = getcursor(trans_unit)
        push!(ctx.cursor_stack, root_cursor)
        header = spelling(root_cursor)
        # loop over all of the child cursors and wrap them, if appropriate.
        ctx.children = children(root_cursor)
        for (i, child) in enumerate(ctx.children)
            # Skip everything except function declarations
            if !isa(child, CLFunctionDecl)
                continue
            end
            # Skip functions that aren't defined in this header
            if filename(child) != header
                continue
            end

            push!(symbol_names, spelling(child))
        end
        empty!(ctx.api_buffer)  # clean up api_buffer for the next header
    end

    # Drop duplicate `__`-prefixed functions
    @info("Extracted $(length(symbol_names)) symbols")
    return sort(symbol_names)
end

# By default, we extract symbols corresponding to an LP64 (32-bit) x86_64-linux-gnu build of OpenBLAS
symbol_names = extract_symbol_names()

open(joinpath(LBT_ROOT, "src", "exported_funcs.inc"), "w") do io
    println(io, "// This file generated by `ext/clang/generate_func_list.jl`")
    println(io, "#ifndef EXPORTED_FUNCS")
    println(io, "#define EXPORTED_FUNCS(XX) \\")
    for name in symbol_names
        println(io, "    XX($(name)) \\")
    end
    println(io)
    println(io, "#endif")
end

max_stars_repo_name: Gnimuc/libblastrampoline

max_issues_repo_name: Gnimuc/libblastrampoline

max_forks_repo_name: Gnimuc/libblastrampoline


text:
using Cooccurence
using Test

@testset "Cooccurence.jl" begin
    # Write your tests here.
end


text:
using Test
using ThreeBodyTB
using Suppressor

#basic loading

#name of test directory
TESTDIR=ThreeBodyTB.TESTDIR

@testset "basic crystal dftout tests" begin
    types=["Li"]
    pos=zeros((1,3))
    A=[ [10.0 0 0]; [0 10.0 0]; [ 0 0 10.0]]

    c=makecrys(A, pos, types, units="Bohr")

    @test c.A == A
    @test c.coords == pos
    @test c.types == types

    c2=makecrys("$TESTDIR/POSCAR_testing")
    @test c2.types == ["Cl", "N","N","Sr", "Sr"]

    c3=makecrys("$TESTDIR/qe_testing.in")
    @test c3.types == ["Si", "C"]

    forces=ones((1,3))
    energy=0.0
    energy_smear = 0.0
    d=makedftout(c, energy, energy_smear, forces, zeros(3,3))
    @test energy == d.energy
    @test forces == d.forces

    d2=makedftout(A,pos,types, energy, energy_smear, forces, zeros(3,3))
    @test d.forces == d2.forces

    #supposed to throw an error
    @suppress begin
        forces2=zeros((2,3))
        @test_throws ErrorException("Error forces") makedftout(c, energy,energy_smear, forces2, zeros(3,3))
    end

    dft_out = ThreeBodyTB.QE.loadXML("$TESTDIR/dimer.save/")

    convert_ryd_ha = 2.0

    @test convert_ryd_ha*-1.428445472333963e1 ≈ dft_out.energy
    @test convert_ryd_ha*-6.425031283667987e-2 ≈ dft_out.bandstruct.efermi
    @test convert_ryd_ha*-1.832485909355843e0 ≈ dft_out.bandstruct.eigs[1,1]

end

@testset "advanced crystal tests" begin
    @suppress begin
        types=["Li"]
        pos=zeros((1,3))
        A=[ [10.0 0 0]; [0 10.0 0]; [ 0 0 10.0]]

        c=makecrys(A, pos, types, units="Bohr")

#        println(c)
        c2 = c * [2,2,2]

        ThreeBodyTB.CrystalMod.generate_random_distortion(c2,0.01,0.01)
        ThreeBodyTB.CrystalMod.write_poscar(c2,"$TESTDIR/POSCAR_tmp")
        rm("$TESTDIR/POSCAR_tmp")

        ThreeBodyTB.CrystalMod.write_efs(c, 0.0, zeros(1,3), zeros(3,3), "$TESTDIR/qe.out")
        rm("$TESTDIR/qe.out")

        @test 1 == 1

    end
end

if false

    #run QE . You actually need quantum espresso working to run this test, not appropriate for CI.
    @testset "run QE" begin
        types=["Li"]
        pos=zeros((1,3))
        A=[ [6.0 0 0]; [0 6.0 0]; [ 0 0 6.0]]

        c=makecrys(A, pos, types, units="Bohr")

        d = ThreeBodyTB.DFT.runSCF(c, nprocs=4)

        @test -14.35143737 ≈ d.energy

    end
end


text:
import GitLab
using Base.Test
myauth = GitLab.authenticate(ENV["GITLAB_AUTH"]) # don't hardcode your access tokens!
println("Authentication successful")
options = Dict("private_token" => myauth.token)

myrepo = GitLab.repo_by_name("TestProject1"; headers=options)
## @show myrepo

r = GitLab.star(myrepo; params=options)
## @show r
## r.status of 201 means a new star and a response of 304 means that it was already present
@test  r.status == 201 || r.status == 304

r = GitLab.unstar(myrepo; params=options)
## @show r
@test r.status == 200

println("Activity Tests Done !!!")


text:
using Documenter, MRphy

makedocs(
    modules = [MRphy],
    sitename="MRphy.jl",
    authors = "Tianrui Luo",
)

deploydocs(
    repo = "github.com/tianrluo/MRphy.jl.git"
)


text:
#!/usr/bin/env julia

##########################################################################
# run_convergence.jl
#
# Run all convergence simulations.
##########################################################################
include("run_convergence-ibc.jl")
include("run_convergence-indoors.jl")
include("run_convergence-randomlargescale.jl")
include("run_convergence-triangularmacro.jl")
include("run_convergence-triangularhetnet.jl")


text:
module ConicBenchmarkUtilities

using GZip

using SparseArrays
using LinearAlgebra

# this is required because findnz does not support arrays by default in julia v1
function SparseArrays.findnz(A::AbstractMatrix)
    I = findall(!iszero, A)
    return (getindex.(I, 1), getindex.(I, 2), A[I])
end

export readcbfdata, cbftompb, mpbtocbf, writecbfdata
export remove_zero_varcones, socrotated_to_soc, remove_ints_in_nonlinear_cones, dualize

include("cbf_input.jl")
include("cbf_output.jl")
include("mpb.jl")
include("preprocess_mpb.jl")
include("convex_to_cbf.jl")
include("jump_to_cbf.jl")

end # module

max_stars_repo_name: mlubin/ConicBenchmarkUtilities

max_issues_repo_name: mlubin/ConicBenchmarkUtilities

max_forks_repo_name: mlubin/ConicBenchmarkUtilities


text:
# Parity operators based on MPB code:
# https://github.com/NanoComp/mpb/blob/master/src/maxwell/maxwell_constraints.c
using LoopVectorization
function _𝓟x!(Ho::AbstractArray{Complex{T},3},Hi::AbstractArray{Complex{T},3},grid::Grid{2}) where T<:Real
    Nx,Ny = size(grid)
	myzero = zero(eltype(Ho))
	temp1 = zero(eltype(Ho))
	temp2 = zero(eltype(Ho))
    @avx for iy ∈ 1:Ny
		# Ho[1,1,iy] = -Hi[1,1,iy]
		# Ho[2,1,iy] = myzero #zero(Complex{T})
		# Ho[1,1,iy] = Hi[1,1,iy] #myzero
		# Ho[2,1,iy] = Hi[2,1,iy] #zero(Complex{T})
		for ix ∈ 1:(Nx÷2)-1, l in 0:0 #1:((Nx÷2)-1), l in 0:0
			temp1 		= 0.5*( Hi[1+l,ix+1,iy] - Hi[1+l,Nx-ix+1,iy] )
			temp2 		= 0.5*( Hi[2+l,ix+1,iy] + Hi[2+l,Nx-ix+1,iy] )
			Ho[1+l,ix+1,iy] 	= temp1
			Ho[2+l,ix+1,iy] 	= temp2
			Ho[1+l,Nx-ix+1,iy] = -temp1
			Ho[2+l,Nx-ix+1,iy] = temp2
		end
	end
	return Ho
end

function _𝓟x̄!(Ho::AbstractArray{Complex{T},3},Hi::AbstractArray{Complex{T},3},grid::Grid{2}) where T<:Real
    Nx,Ny = size(grid)
	temp1 = zero(eltype(Ho))
	temp2 = zero(eltype(Ho))
	myzero = zero(eltype(Ho))
    @avx for iy ∈ 1:Ny
		# Ho[1,1,iy] = Hi[1,1,iy]
		# Ho[2,1,iy] = Hi[2,1,iy] #myzero #zero(Complex{T})
		# Ho[1,1,iy] = myzero
		# Ho[2,1,iy] = -Hi[2,1,iy] #zero(Complex{T})
		for ix ∈ 1:(Nx÷2)-1, l in 0:0 # 1:((Nx÷2)-1), l in 0:0
			temp1 		= 0.5*( Hi[1+l,ix+1,iy] + Hi[1+l,Nx-ix+1,iy] )
			temp2 		= 0.5*( Hi[2+l,ix+1,iy] - Hi[2+l,Nx-ix+1,iy] )
			Ho[1+l,ix+1,iy] 	= temp1
			Ho[2+l,ix+1,iy] 	= temp2
			Ho[1+l,Nx-ix+1,iy] = temp1
			Ho[2+l,Nx-ix+1,iy] = -temp2
		end
	end
	return Ho
end
𝓟x!(grid::Grid{2}) = LinearMap{ComplexF64}((Ho,Hi) -> vec(_𝓟x!(reshape(Ho,(2,grid.Nx,grid.Ny)),reshape(Hi,(2,grid.Nx,grid.Ny)),grid)),*(2,grid.Nx,grid.Ny),ishermitian=false,ismutating=true)
𝓟x!(grid::Grid{3}) = LinearMap{ComplexF64}((Ho,Hi) -> vec(_𝓟x!(reshape(Ho,(2,grid.Nx,grid.Ny,grid.Nz)),reshape(Hi,(2,grid.Nx,grid.Ny,grid.Nz)),grid)),*(2,grid.Nx,grid.Ny,grid.Nz),ishermitian=false,ismutating=true)
𝓟x̄!(grid::Grid{2}) = LinearMap{ComplexF64}((Ho,Hi) -> vec(_𝓟x̄!(reshape(Ho,(2,grid.Nx,grid.Ny)),reshape(Hi,(2,grid.Nx,grid.Ny)),grid)),*(2,grid.Nx,grid.Ny),ishermitian=false,ismutating=true)
𝓟x̄!(grid::Grid{3}) = LinearMap{ComplexF64}((Ho,Hi) -> vec(_𝓟x̄!(reshape(Ho,(2,grid.Nx,grid.Ny,grid.Nz)),reshape(Hi,(2,grid.Nx,grid.Ny,grid.Nz)),grid)),*(2,grid.Nx,grid.Ny,grid.Nz),ishermitian=false,ismutating=true)

##

eigind=1
geom = rwg(p)
Hr = reshape(HH[:,eigind],(2,Ns...));
E0 = E⃗(k,Hr,ω,geom,grid; svecs=false, normalized=false);
H0 = H⃗(k,Hr,ω,geom,grid; svecs=false, normalized=false);

HrPp = copy(Hr)
_𝓟x!(HrPp,copy(Hr),ms.grid);
# HrPp = reshape(HPp,(2,Ns...)); similar(Hr);
HrPn = copy(Hr)
_𝓟x̄!(HrPn,copy(Hr),ms.grid);
# HrPn = reshape(HPp,(2,Ns...)); similar(Hr);

H⃗1[] = H⃗(k1,Hr,ω,geom,grid; svecs=false, normalized=false)
Hx1[] = real(H⃗1[][1,:,:])
Hy1[] = real(H⃗1[][2,:,:])

H⃗Pp[] = H⃗(k1,HrPp,ω,geom,grid; svecs=false, normalized=false)
HxPp[] = real(H⃗Pp[][1,:,:])
HyPp[] = real(H⃗Pp[][2,:,:])

H⃗Pn[] = H⃗(k1,HrPn,ω,geom,grid; svecs=false, normalized=false)
HxPn[] = real(H⃗Pn[][1,:,:])
HyPn[] = real(H⃗Pn[][2,:,:])

E⃗1[] = E⃗(k1,Hr,ω,geom,grid; svecs=false, normalized=false)
Ex1[] = real(E⃗1[][1,:,:])
Ey1[] = real(E⃗1[][2,:,:])

E⃗Pp[] = E⃗(k1,HrPp,ω,geom,grid; svecs=false, normalized=false)
ExPp[] = real(E⃗Pp[][1,:,:])
EyPp[] = real(E⃗Pp[][2,:,:])

E⃗Pn[] = E⃗(k1,HrPn,ω,geom,grid; svecs=false, normalized=false)
ExPn[] = real(E⃗Pn[][1,:,:])
EyPn[] = real(E⃗Pn[][2,:,:])

Hs = reshape([  Hx1,  Hy1, HxPp ,HyPp ,  HxPn ,HyPn ],(2,3))
Es = reshape([  Ex1,  Ey1, ExPp ,EyPp ,  ExPn ,EyPn ],(2,3))
Hmagmax = [lift(A->maximum(abs,A),X) for X in Hs]
Emagmax = [lift(A->maximum(abs,A),X) for X in Es]

fig = Figure()
xs = x(ms.grid)
ys = y(ms.grid)

pos_H = [fig[i,j] for i=1:2,j=1:3]
ax_H = [Axis(pos[1,1]) for pos in pos_H]
pos_E = [fig[i,j] for i=3:4,j=1:3]
ax_E = [Axis(pos[1,1]) for pos in pos_E]

# nx = sqrt.(getindex.(inv.(ms.M̂.ε⁻¹),1,1))
cmaps_H = [:diverging_bwr_40_95_c42_n256, :diverging_protanopic_deuteranopic_bwy_60_95_c32_n256]
cmaps_E = [:diverging_bkr_55_10_c35_n256, :diverging_bky_60_10_c30_n256]
label_base = reshape( ["x","y","xPp","yPp","xPn","yPn",], (2,3))
labels_H = "H".*label_base
labels_E = "E".*label_base
heatmaps_H = [heatmap!(ax_H[i,j], xs, ys, Hs[i,j],colormap=cmaps_H[i],label=labels_H[i,j],colorrange=(-to_value(Hmagmax[i,j]),to_value(Hmagmax[i,j]))) for i=1:2,j=1:3]
heatmaps_E = [heatmap!(ax_E[i,j], xs, ys, Es[i,j],colormap=cmaps_E[i],label=labels_E[i,j],colorrange=(-to_value(Emagmax[i,j]),to_value(Emagmax[i,j]))) for i=1:2,j=1:3]
cbars_H = [Colorbar(pos_H[i,j][1, 2], heatmaps_H[i,j],  width=20 ) for i=1:2,j=1:3]
cbars_E = [Colorbar(pos_E[i,j][1, 2], heatmaps_E[i,j],  width=20 ) for i=1:2,j=1:3]

map( (axx,ll)->text!(axx,ll,position=(-1.4,1.1),textsize=0.7,color=:black), ax_H, labels_H )
map( (axx,ll)->text!(axx,ll,position=(-1.4,1.1),textsize=0.7,color=:white), ax_E, labels_E )
ax_all = vcat(ax_H,ax_E)
hidexdecorations!.(ax_all[1:end-1,:])
hideydecorations!.(ax_all[:,2:end])
[axx.xlabel= "x [μm]" for axx in ax_all[1:end-1,:]]
[axx.ylabel= "y [μm]" for axx in ax_all[:,1]]
[ axx.aspect=DataAspect() for axx in ax_all ]
linkaxes!(ax_all...)

txt= fig[5,2] = indicator(fig)
# for axx in ax_all
# 	on(mouseposition) do mpos
#
# end

fig

##

function _𝓟x!(Ho::AbstractArray{Complex{T},4},Hi::AbstractArray{Complex{T},4},grid::Grid{3}) where T<:Real
    Nx,Ny,Nz = size(grid)
	temp1 = zero(eltype(Ho))
	temp2 = zero(eltype(Ho))
    @avx for iz ∈ 1:Nz, iy ∈ 1:Ny, ix ∈ 1:((Nx÷2)-1), l in 0:0
		temp1 		= 0.5*( Hi[1+l,ix+1,iy,iz] - Hi[1+l,Nx-ix+1,iy,iz] )
		temp2 		= 0.5*( Hi[2+l,ix+1,iy,iz] + Hi[2+l,Nx-ix+1,iy,iz] )
		Ho[1+l,ix+1,iy,iz] 	= temp1
		Ho[2+l,ix+1,iy,iz] 	= temp2
		Ho[1+l,Nx-ix+1,iy,iz] = -temp1
		Ho[2+l,Nx-ix+1,iy,iz] = temp2
	end
	return Ho
end

function _𝓟x̄!(Ho::AbstractArray{Complex{T},4},Hi::AbstractArray{Complex{T},4},grid::Grid{3}) where T<:Real
    Nx,Ny,Nz = size(grid)
	temp1 = zero(eltype(Ho))
	temp2 = zero(eltype(Ho))
    @avx for iz ∈ 1:Nz, iy ∈ 1:Ny, ix ∈ 1:((Nx÷2)-1), l in 0:0
		temp1 		= 0.5*( Hi[1+l,ix+1,iy,iz] + Hi[1+l,Nx-ix+1,iy,iz] )
		temp2 		= 0.5*( Hi[2+l,ix+1,iy,iz] - Hi[2+l,Nx-ix+1,iy,iz] )
		Ho[1+l,ix+1,iy,iz] 	= temp1
		Ho[2+l,ix+1,iy,iz] 	= temp2
		Ho[1+l,Nx-ix+1,iy,iz] = temp1
		Ho[2+l,Nx-ix+1,iy,iz] = -temp2
	end
	return Ho
end

##
Nx = 128
Ny = 128

for ix ∈ 0:(Nx÷2) #1:((Nx÷2)-1), l in 0:0
	ij = ix #(iy-1) * Ny + ix-1
    ij2 = (ix > 0 ? Nx - ix : 0) # (iy-1) * Ny + ((ix-1) > 0 ? Nx - (ix-1) : 0)
	println("")
	println("ij, ij2: $ij , $ij2")
	println("")
end

for ix ∈ 1:(Nx÷2)+1 #1:((Nx÷2)-1), l in 0:0
	# ix1 = ix + 1
	# ix2 = (ix > 0 ? Nx - ix : 0) + 1
	ix1 = ix
	ix2 = (ix > 1 ? Nx - (ix-2) : 1)
	println("")
	println("ix1, ix2: $ix1 , $ix2")
	println("")
end

##
ms = ModeSolver(1.45, rwg(p), grid, nev=4)
ω = 0.75
k,H⃗ = solve_k(ω,rwg(p),grid)
ms.H⃗ |> size
##
H1 = copy(ms.H⃗[:,1])
H2 = copy(ms.H⃗[:,2])
H3 = copy(ms.H⃗[:,3])
H4 = copy(ms.H⃗[:,4])
Px! = 𝓟x!(ms.grid)
Px̄! = 𝓟x̄!(ms.grid)

H1c = copy(H1)
Px!*H1c

H1c = copy(H1)
mul!(H1c,Px!,H1c)

H1c = copy(H1)
maximum(abs2.(mul!(H1c,Px!,H1c) .- H1))

H1c = copy(H1)
maximum(abs2.(mul!(H1c,Px̄!,H1c) .- H1))

H1c = copy(H1)
maximum(abs2.(Px!*H1c .- H1))

H1c = copy(H1)
maximum(abs2.(Px̄!*H1c .- H1))

H1c = copy(H1)
maximum(abs2.(Px!2(Hrc,gr,-1) .- Hr))
using LinearAlgebra
mul!(H1c,Px!,H1c)

vec(_𝓟x!(reshape(H1c,(2,grid.Nx,grid.Ny)),grid))
33
33
33

##
# function _𝓟y!(H::AbstractArray{Complex{T},3},grid::Grid{2},parity=1) where T<:Real
#     Nx,Ny = size(grid)
# 	temp1 = zero(eltype(H))
# 	temp2 = zero(eltype(H))
#     @avx for iy ∈ 1:((Ny÷2)-1), ix ∈ 1:Nx, l in 0:0
# 		temp1 		= 0.5*( H[1+l,ix,iy+1] - parity * H[1+l,ix,Ny-iy+1] )
# 		temp2 		= 0.5*( H[2+l,ix,iy+1] + parity * H[2+l,ix,Ny-iy+1] )
# 		H[1+l,ix,iy+1] 	= temp1
# 		H[2+l,ix,iy+1] 	= temp2
# 		H[1+l,ix,Ny-iy+1] = -parity 	* 	temp1
# 		H[2+l,ix,Ny-iy+1] = parity 		* 	temp2
# 	end
# 	return H
# end
#
# function _𝓟y!(H::AbstractArray{Complex{T},4},grid::Grid{3},parity=1) where T<:Real
#     Nx,Ny,Nz = size(grid)
# 	temp1 = zero(eltype(H))
# 	temp2 = zero(eltype(H))
#     @avx for iz ∈ 1:Nz, iy ∈ 1:((Ny÷2)-1), ix ∈ 1:Nx, l in 0:0
# 		temp1 		= 0.5*( H[1+l,ix,iy+1,iz] - parity * H[1+l,ix,Ny-iy+1,iz] )
# 		temp2 		= 0.5*( H[2+l,ix,iy+1,iz] + parity * H[2+l,ix,Ny-iy+1,iz] )
# 		H[1+l,ix,iy+1,iz] 	= temp1
# 		H[2+l,ix,iy+1,iz] 	= temp2
# 		H[1+l,ix,Ny-iy+1,iz] = -parity 	* 	temp1
# 		H[2+l,ix,Ny-iy+1,iz] = parity 		* 	temp2
# 	end
# 	return H
# end
#
# function _𝓟ₓ!(H::StructArray{Complex{T},3},grid::Grid{2},parity=1) where T<:Real
# 	𝓟ₓ!(H.re,grid,parity)
# 	𝓟ₓ!(H.im,grid,parity)
# 	return nothing
# end

## test code
#
# Hrc = copy(Hr)
# maximum(abs2.(Px!1(Hrc,gr,-1) .- Hr))
# Hrc = copy(Hr)
# maximum(abs2.(Px!2(Hrc,gr,-1) .- Hr))
#
# Hr2c = copy(Hr2)
# maximum(abs2.(Px!1(Hr2c,gr,1) .- Hr2))
# Hr2c = copy(Hr2)
# maximum(abs2.(𝓟ₓ!3(Hr2c,gr,1) .- Hr2))
# # maximum(abs2.(Px!2(Hr2c,gr,1) .- Hr2))
#
# Hrstc = copy(Hrst)
# maximum(abs2.(Px!1(Hrstc,gr,-1) .- Hrst))
# Hrstc = copy(Hrst)
# maximum(abs2.(𝓟ₓ!(Hrstc,gr,-1) .- Hrst))
#
# @btime $𝓟ₓ!3($Hr2c,$gr,1)
# @btime 𝓟ₓ!($Hrstc,$gr,-1)
#
# 3
# 3
#
#  # test
# ##
##
using GLMakie
using GLMakie: lines!
H⃗ = OptiMode.H⃗
rwg(x) = ridge_wg_partial_etch(x[1],x[2],x[3],x[4],0.5,LNx,SiO₂,Δx,Δy) # partially etched ridge waveguide with dispersive materials, x[3] is partial etch fraction of top layer, x[3]*x[2] is etch depth, remaining top layer thickness = x[2]*(1-x[3]).
p = [
       1.7,                #   top ridge width         `w_top`         [μm]
       0.7,                #   ridge thickness         `t_core`        [μm]
       0.5,                #   ridge thickness         `t_core`        [μm]
       π / 14.0,           #   ridge sidewall angle    `θ`             [radian]
               ];
geom = rwg(p)
# ms = ModeSolver(1.45, geom, grid, nev=4,constraint=𝓟x!(grid))
ms = ModeSolver(1.45, geom, grid, nev=4,) # constraint=𝓟x!(grid))
ω = 0.75
k,Hv = solve_k(ms,ω,rwg(p),nev=4,eigind=1)
k_old = 1.5167250170625406
k_new = k
k1 = copy(k)
abs(k_old - k_new) / k_new
HrN = Node(copy(reshape(ms.H⃗[:,1],(2,Ns...))))
##
using ArrayInterface
using RecursiveArrayTools
gr = ms.grid
shape = size(grid)
m⃗ = ms.M̂.m⃗ |> vec |> VectorOfArray
n⃗ = ms.M̂.n⃗ |> vec |> VectorOfArray
gv = g⃗(gr) |> vec |> VectorOfArray
xv = x⃗(gr) |> vec |> VectorOfArray
(xx,yy,zz),(gx,gy,gz),(nx,ny,nz),(mx,my,mz) = vecarr_to_vectors.((xv,gv,n⃗,m⃗))
xxr,yyr,zzr,gxr,gyr,gzr,nxr,nyr,nzr,mxr,myr,mzr = reshape.((xx,yy,zz,gx,gy,gz,nx,ny,nz,mx,my,mz),(shape,))

# length of eigenvector (H-field in recip space/transverse pol. basis is always N(grid)*2 = (N(grid)=Nx×Ny×Nz recip. lattice grid points) × (2 transverse polarizations for each plane wave)
N(gr)*2 == size(HH,1)
# x-mirror symmetry of flattened (by `vec`) x-value array (at each xyz point)
# @assert xx[end:-1:begin] == -xx
# x-mirror symmetry of flattened (by `vec`) y-value array (at each xyz point)
# @assert yy[end:-1:begin] == -yy
# x-mirror symmetry of flattened (by `vec`) gx-value array (at each xyz point)
# @assert gx[end:-1:begin] ≈ (-gx .- (1.0/gr.Δx))
# lines!(axes_g[1],-gx[end:-1:begin],color=:green)

# @assert mx[1:Nx÷2][2:end] == -mx[Nx:-1:Nx÷2+1][1:end-1]
# @assert mxs[1:Nx÷2][2:end] == -mxs[Nx:-1:Nx÷2+1][1:end-1]
#
# @assert gxr[1:Nx÷2,:][2:end,:] == -gxr[Nx:-1:Nx÷2+1,:][1:end-1,:]
# @assert mxr[1:Nx÷2,:][2:end,:] == -mxr[Nx:-1:Nx÷2+1,:][1:end-1,:]
#
# @assert gyr[:,1:Ny÷2][:,2:end] == -gyr[:,Ny:-1:Ny÷2+1][:,1:end-1]
# @assert myr[:,1:Ny÷2][:,2:end] == -myr[:,Ny:-1:Ny÷2+1][:,1:end-1]

# odd x-parity example
maximum(abs2.(Hr))
maximum(abs2.(Hr[1,1:Nx÷2,:][2:end,:] .- Hr[1,Nx:-1:Nx÷2+1,:][1:end-1,:]))	# <--- these are zero everywhere for odd parity
maximum(abs2.(Hr[1,1:Nx÷2,:][2:end,:] .+ Hr[1,Nx:-1:Nx÷2+1,:][1:end-1,:]))	# <--- these are nonzero everywhere for odd parity
maximum(abs2.(Hr[2,1:Nx÷2,:][2:end,:] .- Hr[2,Nx:-1:Nx÷2+1,:][1:end-1,:]))	# <--- these are nonzero everywhere for odd parity
maximum(abs2.(Hr[2,1:Nx÷2,:][2:end,:] .+ Hr[2,Nx:-1:Nx÷2+1,:][1:end-1,:]))	# <--- these are zero everywhere for odd parity

# direct indexing
Hr[1,2:Nx÷2,:] == Hr[1,1:Nx÷2,:][2:end,:]
Hr[1,Nx:-1:Nx÷2+2,:] == Hr[1,Nx:-1:Nx÷2+1,:][1:end-1,:]

# odd parity with direct indexing
parity = -1	# -1 => odd, 1 => even
Hrsym = similar(Hr)
Hrsym[1,2:Nx÷2,:] .= 0.5*( Hr[1,2:Nx÷2,:] - parity * Hr[1,Nx:-1:Nx÷2+2,:] )
Hrsym[2,2:Nx÷2,:] .= 0.5*( Hr[2,2:Nx÷2,:] + parity * Hr[2,Nx:-1:Nx÷2+2,:] )
Hrsym[1,Nx:-1:Nx÷2+2,:] .= -parity * Hrsym[1,2:Nx÷2,:]
Hrsym[2,Nx:-1:Nx÷2+2,:] .=  parity * Hrsym[2,2:Nx÷2,:]

maximum(abs2.(Hrsym))
maximum(abs2.(Hrsym[1,1:Nx÷2,:][2:end,:] .- Hrsym[1,Nx:-1:Nx÷2+1,:][1:end-1,:]))	# <--- these are zero everywhere for odd parity
maximum(abs2.(Hrsym[1,1:Nx÷2,:][2:end,:] .+ Hrsym[1,Nx:-1:Nx÷2+1,:][1:end-1,:]))	# <--- these are nonzero everywhere for odd parity
maximum(abs2.(Hrsym[2,1:Nx÷2,:][2:end,:] .- Hrsym[2,Nx:-1:Nx÷2+1,:][1:end-1,:]))	# <--- these are nonzero everywhere for odd parity
maximum(abs2.(Hrsym[2,1:Nx÷2,:][2:end,:] .+ Hrsym[2,Nx:-1:Nx÷2+1,:][1:end-1,:]))	# <--- these are zero everywhere for odd parity
# these are zero everywhere for when the parity operation left Hrsym ≈ Hr:
maximum(abs2.(Hrsym[1,2:Nx÷2,:] .- Hr[1,2:Nx÷2,:]))
maximum(abs2.(Hrsym[1,Nx:-1:Nx÷2+2,:] .- Hr[1,Nx:-1:Nx÷2+2,:]))
maximum(abs2.(Hrsym[2,2:Nx÷2,:] .- Hr[2,2:Nx÷2,:]))
maximum(abs2.(Hrsym[2,Nx:-1:Nx÷2+2,:] .- Hr[2,Nx:-1:Nx÷2+2,:]))

# @assert Hr[1,1:Nx÷2,:][2:end,:] == -Hr[1,Nx:-1:Nx÷2+1,:][1:end-1,:]
# @assert Hr[2,:,1:Ny÷2][:,2:end] == Hr[2,:,Ny:-1:Ny÷2+1][:,1:end-1]

# even x-parity example
Hr2 = reshape(HH[:,1],(2,Nx,Ny))
maximum(abs2.(Hr2))
maximum(abs2.(Hr2[1,1:Nx÷2,:][2:end,:] .- Hr2[1,Nx:-1:Nx÷2+1,:][1:end-1,:]))	# <--- these are nonzero everywhere for even parity
maximum(abs2.(Hr2[1,1:Nx÷2,:][2:end,:] .+ Hr2[1,Nx:-1:Nx÷2+1,:][1:end-1,:]))	# <--- these are zero everywhere for even parity
maximum(abs2.(Hr2[2,1:Nx÷2,:][2:end,:] .- Hr2[2,Nx:-1:Nx÷2+1,:][1:end-1,:]))	# <--- these are zero everywhere for even parity
maximum(abs2.(Hr2[2,1:Nx÷2,:][2:end,:] .+ Hr2[2,Nx:-1:Nx÷2+1,:][1:end-1,:]))	# <--- these are nonzero everywhere for even parity

# even parity with direct indexing
parity = 1	# -1 => odd, 1 => even
Hr2sym = similar(Hr2)
Hr2sym[1,2:Nx÷2,:] .= 0.5*( Hr2[1,2:Nx÷2,:] - parity * Hr2[1,Nx:-1:Nx÷2+2,:] )
Hr2sym[2,2:Nx÷2,:] .= 0.5*( Hr2[2,2:Nx÷2,:] + parity * Hr2[2,Nx:-1:Nx÷2+2,:] )
Hr2sym[1,Nx:-1:Nx÷2+2,:] .= -parity * Hr2sym[1,2:Nx÷2,:]
Hr2sym[2,Nx:-1:Nx÷2+2,:] .=  parity * Hr2sym[2,2:Nx÷2,:]

maximum(abs2.(Hr2sym))
maximum(abs2.(Hr2sym[1,1:Nx÷2,:][2:end,:] .- Hr2sym[1,Nx:-1:Nx÷2+1,:][1:end-1,:]))	# <--- these are zero everywhere for odd parity
maximum(abs2.(Hr2sym[1,1:Nx÷2,:][2:end,:] .+ Hr2sym[1,Nx:-1:Nx÷2+1,:][1:end-1,:]))	# <--- these are nonzero everywhere for odd parity
maximum(abs2.(Hr2sym[2,1:Nx÷2,:][2:end,:] .- Hr2sym[2,Nx:-1:Nx÷2+1,:][1:end-1,:]))	# <--- these are nonzero everywhere for odd parity
maximum(abs2.(Hr2sym[2,1:Nx÷2,:][2:end,:] .+ Hr2sym[2,Nx:-1:Nx÷2+1,:][1:end-1,:]))	# <--- these are zero everywhere for odd parity
# these are zero everywhere for when the parity operation left Hr2sym ≈ Hr2:
maximum(abs2.(Hr2sym[1,2:Nx÷2,:] .- Hr2[1,2:Nx÷2,:]))
maximum(abs2.(Hr2sym[1,Nx:-1:Nx÷2+2,:] .- Hr2[1,Nx:-1:Nx÷2+2,:]))
maximum(abs2.(Hr2sym[2,2:Nx÷2,:] .- Hr2[2,2:Nx÷2,:]))
maximum(abs2.(Hr2sym[2,Nx:-1:Nx÷2+2,:] .- Hr2[2,Nx:-1:Nx÷2+2,:]))

##
gxs = fftshift(gxr)|>vec
mxs = fftshift(mxr)|>vec
nxs = fftshift(nxr)|>vec
gys = fftshift(gyr)|>vec
mys = fftshift(myr)|>vec
nys = fftshift(nyr)|>vec
lines!(axes_g[1],gxs,color=:black)
lines!(axes_g[2],gys,color=:black)
lines!(axes_m[1],mxs,color=:black)
lines!(axes_m[2],mys,color=:black)
lines!(axes_n[1],nxs,color=:black)
lines!(axes_n[2],nys,color=:black)
##
# x-mirror symmetry of flattened (by `vec`) gy-value array (at each xyz point)
@assert gy[end:-1:begin] ≈ (-gy .- (1.0/gr.Δy))

##
fig = Figure()
axes_xy = fig[1:2,1] = [Axis(fig) for i=1:2]
axes_g = fig[1:2,2] = [Axis(fig) for i=1:2]
axes_m = fig[1,3:5] = [Axis(fig) for i=1:3]
axes_n = fig[2,3:5] = [Axis(fig) for i=1:3]
lxys = map((axx,xy,cc,ll)->(lines!(axx,xy,color=cc,label=ll);axislegend(axx)), axes_xy, [xx,yy], [:purple,:black],["x","y"])
lgs = map((axx,gg,cc,ll)->(lines!(axx,gg,color=cc,label=ll);axislegend(axx)), axes_g, [gx,gy], [:magenta,:cyan],["gx","gy"])
lns = map((axx,nn,cc,ll)->(lines!(axx,nn,color=cc,label=ll);axislegend(axx)), axes_n, [nx,ny,nz], [:red,:blue,:green],["nx","ny","nz"])
lms = map((axx,mm,cc,ll)->(lines!(axx,mm,color=cc,label=ll);axislegend(axx)), axes_m, [mx,my,mz], [:red,:blue,:green],["mx","my","mz"])
linkxaxes!(vcat(axes_xy,axes_g,axes_m,axes_n)...)
fig
##

##
HH = copy(ms.H⃗) # recip. space H-field soln. copy
Ns = size(ms.grid)
eigind=1
Hr = reshape(HH[:,eigind],(2,Ns...))
geom = rwg(p)
HrN[] = Hr; #reshape(HH[:,eigind],(2,Ns...));
E0 = E⃗(k,to_value(HrN),ω,geom,grid; svecs=false, normalized=false);
H0 = H⃗(k,to_value(HrN),ω,geom,grid; svecs=false, normalized=false);

E⃗1	= Node(E0)
Ex1	= Node(E0[1,:,:])
Ey1	= Node(E0[1,:,:])
E⃗Pp	= Node(E0)
ExPp	= Node(E0[1,:,:])
EyPp	= Node(E0[1,:,:])
E⃗Pn	= Node(E0)
ExPn	= Node(E0[1,:,:])
EyPn	= Node(E0[1,:,:])

H⃗1	= Node(H0)
Hx1	= Node(H0[1,:,:])
Hy1	= Node(H0[1,:,:])
H⃗Pp	= Node(H0)
HxPp	= Node(H0[1,:,:])
HyPp	= Node(H0[1,:,:])
H⃗Pn	= Node(H0)
HxPn	= Node(H0[1,:,:])
HyPn	= Node(H0[1,:,:])

##
HPp = vec(copy(Hr))
mul!(HPp,𝓟x!(ms.grid),vec(copy(Hr)));
# mul!(HPp,Pxs,vec(copy(Hr)));
HrPp = reshape(HPp,(2,Ns...));

HPn = vec(copy(Hr))
mul!(HPn,𝓟x̄!(ms.grid),vec(copy(Hr)));
# mul!(HPn,Pxbs,vec(copy(Hr)));
HrPn = reshape(HPn,(2,Ns...));

HrN = Node(Hr)

H⃗1[] = H⃗(k1,Hr,ω,geom,grid; svecs=false, normalized=false)
Hx1[] = real(H⃗1[][1,:,:])
Hy1[] = real(H⃗1[][2,:,:])

H⃗Pp[] = H⃗(k1,HrPp,ω,geom,grid; svecs=false, normalized=false)
HxPp[] = real(H⃗Pp[][1,:,:])
HyPp[] = real(H⃗Pp[][2,:,:])

H⃗Pn[] = H⃗(k1,HrPn,ω,geom,grid; svecs=false, normalized=false)
HxPn[] = real(H⃗Pn[][1,:,:])
HyPn[] = real(H⃗Pn[][2,:,:])

E⃗1[] = E⃗(k1,Hr,ω,geom,grid; svecs=false, normalized=false)
Ex1[] = real(E⃗1[][1,:,:])
Ey1[] = real(E⃗1[][2,:,:])

E⃗Pp[] = E⃗(k1,HrPp,ω,geom,grid; svecs=false, normalized=false)
ExPp[] = real(E⃗Pp[][1,:,:])
EyPp[] = real(E⃗Pp[][2,:,:])

E⃗Pn[] = E⃗(k1,HrPn,ω,geom,grid; svecs=false, normalized=false)
ExPn[] = real(E⃗Pn[][1,:,:])
EyPn[] = real(E⃗Pn[][2,:,:])

# on(HrN) do HrN
#
# 	HrPp = copy(to_value(HrN))
# 	_𝓟x!(HrPp,copy(Hr),ms.grid);
# 	# HrPp = reshape(HPp,(2,Ns...)); similar(Hr);
# 	HrPn = copy(to_value(HrN))
# 	_𝓟x̄!(HrPn,copy(Hr),ms.grid);
# 	# HrPn = reshape(HPp,(2,Ns...)); similar(Hr);
#
# 	H⃗1[] = H⃗(k1,to_value(HrN),ω,geom,grid; svecs=false, normalized=false)
# 	Hx1[] = real(H⃗1[][1,:,:])
# 	Hy1[] = real(H⃗1[][2,:,:])
#
#
# 	H⃗Pp[] = H⃗(k1,HrPp,ω,geom,grid; svecs=false, normalized=false)
# 	HxPp[] = real(H⃗Pp[][1,:,:])
# 	HyPp[] = real(H⃗Pp[][2,:,:])
#
# 	H⃗Pn[] = H⃗(k1,HrPn,ω,geom,grid; svecs=false, normalized=false)
# 	HxPn[] = real(H⃗Pn[][1,:,:])
# 	HyPn[] = real(H⃗Pn[][2,:,:])
#
# 	E⃗1[] = E⃗(k1,Hr,ω,geom,grid; svecs=false, normalized=false)
# 	Ex1[] = real(E⃗1[][1,:,:])
# 	Ey1[] = real(E⃗1[][2,:,:])
#
# 	E⃗Pp[] = E⃗(k1,HrPp,ω,geom,grid; svecs=false, normalized=false)
# 	ExPp[] = real(E⃗Pp[][1,:,:])
# 	EyPp[] = real(E⃗Pp[][2,:,:])
#
# 	E⃗Pn[] = E⃗(k1,HrPn,ω,geom,grid; svecs=false, normalized=false)
# 	ExPn[] = real(E⃗Pn[][1,:,:])
# 	EyPn[] = real(E⃗Pn[][2,:,:])
#
# 	Hs .= reshape([  Hx1,  Hy1, HxPp ,HyPp ,  HxPn ,HyPn ],(2,3))
# 	Es .= reshape([  Ex1,  Ey1, ExPp ,EyPp ,  ExPn ,EyPn ],(2,3))
# 	Hmagmax .= [lift(A->maximum(abs,A),X) for X in Hs]
# 	Emagmax .= [lift(A->maximum(abs,A),X) for X in Es]
# end
##
# using LoopVectorization
# Pxr = zeros(eltype(Hr),size(Hr))
# Px̄r = zeros(eltype(Hr),size(Hr))
# @avx for iy ∈ 1:Ny
# 	for ix ∈ 2:(Nx÷2), l in 0:0 #1:((Nx÷2)-1), l in 0:0
# 		ix1 = ix
# 		ix2 = (ix > 1 ? Nx - (ix-2) : 1)
# 		Pxr[1+l,ix1,iy] = 1.0
# 		Pxr[2+l,ix1,iy] = 1.0
# 		Pxr[1+l,ix2,iy] = -1.0
# 		Pxr[2+l,ix2,iy] = 1.0
# 	end
# end
#
# @avx for iy ∈ 1:Ny
# 	for ix ∈ 2:(Nx÷2), l in 0:0 #1:((Nx÷2)-1), l in 0:0
# 		ix1 = ix
# 		ix2 = (ix > 1 ? Nx - (ix-2) : 1)
# 		Px̄r[1+l,ix1,iy] = 1.0
# 		Px̄r[2+l,ix1,iy] = 1.0
# 		Px̄r[1+l,ix2,iy] = 1.0
# 		Px̄r[2+l,ix2,iy] = -1.0
# 	end
# end
#
# Pxv = vec(Pxr) #.* inv(N(grid))
# Px̄v = vec(Px̄r) #.* inv(N(grid))

##
#
# H⃗Px = H⃗(k1,Pxr,ω,geom,grid; svecs=false, normalized=false)
# HPxxr =  sign.(real(H⃗Px[1,:,:]) )
# HPxyr =  sign.(real(H⃗Px[2,:,:]) )
# HPxzr =  sign.(real(H⃗Px[3,:,:]) )
# HPxxi = sign.( imag(H⃗Px[1,:,:]) )
# HPxyi =  sign.(imag(H⃗Px[2,:,:]) )
# HPxzi =  sign.(imag(H⃗Px[3,:,:]) )
#
# fig = Figure()
#
# Hs = reshape([  HPxxr, HPxxi, HPxyr, HPxyi, HPxzr, HPxzi,    ],(2,3))
# Hmagmax = [maximum(abs,A) for A in Hs]
# cmaps_H = [:diverging_bwr_40_95_c42_n256, :diverging_protanopic_deuteranopic_bwy_60_95_c32_n256]
# pos_H = [fig[i,j] for i=1:2,j=1:3]
# ax_H = [Axis(pos[1,1]) for pos in pos_H]
# label_base = reshape( ["xr","xi", "yr", "yi", "zr", "zi",], (2,3))
# labels_H = "H".*label_base
# heatmaps_H = [heatmap!(ax_H[i,j], xs, ys, Hs[i,j],colormap=cmaps_H[i],label=labels_H[i,j],colorrange=(-Hmagmax[i,j],Hmagmax[i,j])) for i=1:2,j=1:3]
# cbars_H = [Colorbar(pos_H[i,j][1, 2], heatmaps_H[i,j],  width=20 ) for i=1:2,j=1:3]
# map( (axx,ll)->text!(axx,ll,position=(-1.4,1.1),textsize=0.7,color=:black), ax_H, labels_H )
# ax_all = vcat(ax_H)
# hidexdecorations!.(ax_all[1:end-1,:])
# hideydecorations!.(ax_all[:,2:end])
# [axx.xlabel= "x [μm]" for axx in ax_all[1:end-1,:]]
# [axx.ylabel= "y [μm]" for axx in ax_all[:,1]]
# [ axx.aspect=DataAspect() for axx in ax_all ]
# linkaxes!(ax_all...)
#
# txt= fig[5,2] = indicator(fig)
# # for axx in ax_all
# # 	on(mouseposition) do mpos
# #
# # end
#
# fig

## interactions

function indicator(ax::Axis,ob)
    register_interaction!(ax, :indicator) do event::MouseEvent, axis
    if event.type === MouseEventTypes.over
        ob[] = event.data
    end
    end
end
function indicator(grid::GridLayout,ob)
    foreach(Axis,grid;recursive=true) do ax
    indicator(ax,ob)
    end
end
function indicator(grid::GridLayout)
    ob = Observable(Point2f0(0.,0.))
    indicator(grid,ob)
    ob
end
function indicator(fig,args...; tellwidth=false, kwargs...)
    Label(
        fig,
        lift(ind->"x: $(ind[1])  y: $(ind[2])",indicator(fig.layout)),
        args...; tellwidth=tellwidth, kwargs...
    )
end

## H and E fields
using AbstractPlotting: mouseposition
using AbstractPlotting

Hs = reshape([  Hx1,  Hy1, HxPp ,HyPp ,  HxPn ,HyPn ],(2,3))
Es = reshape([  Ex1,  Ey1, ExPp ,EyPp ,  ExPn ,EyPn ],(2,3))
Hmagmax = [lift(A->maximum(abs,A),X) for X in Hs]
Emagmax = [lift(A->maximum(abs,A),X) for X in Es]

fig = Figure()
xs = x(ms.grid)
ys = y(ms.grid)

pos_H = [fig[i,j] for i=1:2,j=1:3]
ax_H = [Axis(pos[1,1]) for pos in pos_H]
pos_E = [fig[i,j] for i=3:4,j=1:3]
ax_E = [Axis(pos[1,1]) for pos in pos_E]

# nx = sqrt.(getindex.(inv.(ms.M̂.ε⁻¹),1,1))
cmaps_H = [:diverging_bwr_40_95_c42_n256, :diverging_protanopic_deuteranopic_bwy_60_95_c32_n256]
cmaps_E = [:diverging_bkr_55_10_c35_n256, :diverging_bky_60_10_c30_n256]
label_base = reshape( ["x","y","xPp","yPp","xPn","yPn",], (2,3))
labels_H = "H".*label_base
labels_E = "E".*label_base
heatmaps_H = [heatmap!(ax_H[i,j], xs, ys, Hs[i,j],colormap=cmaps_H[i],label=labels_H[i,j],colorrange=(-to_value(Hmagmax[i,j]),to_value(Hmagmax[i,j]))) for i=1:2,j=1:3]
heatmaps_E = [heatmap!(ax_E[i,j], xs, ys, Es[i,j],colormap=cmaps_E[i],label=labels_E[i,j],colorrange=(-to_value(Emagmax[i,j]),to_value(Emagmax[i,j]))) for i=1:2,j=1:3]
cbars_H = [Colorbar(pos_H[i,j][1, 2], heatmaps_H[i,j],  width=20 ) for i=1:2,j=1:3]
cbars_E = [Colorbar(pos_E[i,j][1, 2], heatmaps_E[i,j],  width=20 ) for i=1:2,j=1:3]

map( (axx,ll)->text!(axx,ll,position=(-1.4,1.1),textsize=0.7,color=:black), ax_H, labels_H )
map( (axx,ll)->text!(axx,ll,position=(-1.4,1.1),textsize=0.7,color=:white), ax_E, labels_E )
ax_all = vcat(ax_H,ax_E)
hidexdecorations!.(ax_all[1:end-1,:])
hideydecorations!.(ax_all[:,2:end])
[axx.xlabel= "x [μm]" for axx in ax_all[1:end-1,:]]
[axx.ylabel= "y [μm]" for axx in ax_all[:,1]]
[ axx.aspect=DataAspect() for axx in ax_all ]
linkaxes!(ax_all...)

txt= fig[5,2] = indicator(fig)
# for axx in ax_all
# 	on(mouseposition) do mpos
#
# end

fig

## E fields only

# Es = reshape([  Ex1,  Ey1, ExPp ,EyPp ,  ExPn ,EyPn ],(2,3))
# fig = Figure()
#
#
# xs = x(ms.grid)
# ys = y(ms.grid)
#
# ax_nx = fig[1,2] = Axis(fig)
# ax_E = fig[2:3,1:3] = reshape( [ Axis(fig) for i=1:6 ], (2,3)) #, title = t) for t in [ "|Eₓ|² @ ω", "|Eₓ|² @ 2ω" ] ]
#
# nx = sqrt.(getindex.(inv.(ms.M̂.ε⁻¹),1,1))
#
# # cmaps_E = [:linear_ternary_red_0_50_c52_n256, :linear_ternary_blue_0_44_c57_n256]
# cmaps_E = [:berlin, :lisbon]
# labels_E = ["rel. |Eₓ|² @ ω","rel. |Ey|² @ ω"]
# heatmaps_E = [heatmap!(ax_E[i,j], xs, ys, Es[i,j],colormap=cmaps_E[i],label=labels_E[i]) for i=1:2,j=1:3]
# hm_nx = heatmap!(ax_nx,xs,ys,nx;colormap=:viridis)
# text!(ax_nx,"nₓ",position=(1.4,1.1),textsize=0.7,color=:white)
# text!(ax_E[1,1],"real Eₓ (ω)",position=(-1.4,1.1),textsize=0.7,color=:white)
# text!(ax_E[2,1],"real Ey (ω)",position=(-1.7,1.1),textsize=0.7,color=:white)
# text!(ax_E[1,2],"real EₓPp (ω)",position=(-1.4,1.1),textsize=0.7,color=:white)
# text!(ax_E[2,2],"real EyPp (ω)",position=(-1.7,1.1),textsize=0.7,color=:white)
# text!(ax_E[1,3],"real EₓPn (ω)",position=(-1.4,1.1),textsize=0.7,color=:white)
# text!(ax_E[2,3],"real EyPn (ω)",position=(-1.7,1.1),textsize=0.7,color=:white)
# cbar_nx = fig[1,4] = Colorbar(fig,hm_nx ) #,label="nₓ")
# cbar_Ex = fig[2,4] = Colorbar(fig,heatmaps_E[1]) #,label="rel. |Eₓ|² @ ω")
# cbar_Ey = fig[3,4] = Colorbar(fig,heatmaps_E[2]) #,label="rel. |Eₓ|² @ 2ω")
# for cb in [cbar_Ex, cbar_Ey, cbar_nx]
#     cb.width=30
#     cb.height=Relative(2/3)
# end
# # label, format
# # hidexdecorations!(ax_E[1])
# hidexdecorations!(ax_nx)
#
# ax_E[1].ylabel = "y [μm]"
# ax_E[2].ylabel = "y [μm]"
# ax_E[2].xlabel = "x [μm]"
#
# ax_nx.aspect=DataAspect()
# all_axes = vcat(ax_nx,ax_E)
# [ aa.aspect=DataAspect() for aa in all_axes ]
# linkxaxes!(all_axes...)
# fig


text:
module Tracking
    using DocStringExtensions, GNSSSignals, StaticArrays, TrackingLoopFilters
    using Unitful: upreferred, Hz, dBHz, ms
    import Base.zero, Base.length

    export
        get_early,
        get_prompt,
        get_late,
        get_correlator,
        get_carrier_doppler,
        get_carrier_phase,
        get_code_doppler,
        get_code_phase,
        get_early_late_sample_shift,
        get_secondary_code_or_bit_found,
        get_state,
        get_cn0,
        track,
        TrackingState,
        NumAnts,
        MomentsCN0Estimator,
        AbstractCN0Estimator,
        get_bits,
        get_num_bits,
        EarlyPromptLateCorrelator,
        VeryEarlyPromptLateCorrelator,
        SecondaryCodeOrBitDetector

    struct NumAnts{x}
    end

    NumAnts(x) = NumAnts{x}()

    include("cn0_estimation.jl")
    include("discriminators.jl")
    include("bit_buffer.jl")
    include("correlator.jl")
    include("secondary_code_or_bit_detector.jl")
    include("tracking_state.jl")
    include("tracking_results.jl")
    include("tracking_loop.jl")
    include("gpsl1.jl")
    include("gpsl5.jl")
    include("galileo_e1b.jl")
end


text:
# These `deserialize` methods are used to be able to deserialize intervals using the
# structure that was used before Intervals 1.3.

function Serialization.deserialize(s::AbstractSerializer, ::Type{Interval{T}}) where T
    left = deserialize(s)
    right = deserialize(s)
    inclusivity = deserialize(s)

    L = bound_type(first(inclusivity))
    R = bound_type(last(inclusivity))

    return Interval{T,L,R}(left, right)
end

function Serialization.deserialize(s::AbstractSerializer, ::Type{AnchoredInterval{P,T}}) where {P,T}
    anchor = deserialize(s)
    inclusivity = deserialize(s)

    L = bound_type(first(inclusivity))
    R = bound_type(last(inclusivity))

    return AnchoredInterval{P,T,L,R}(anchor)
end


text:
function test_domain(out::Float64)
    @test (out >= -1) && (out <= 1)
end

function test_domain(out::Array{Float64})
    @test all((out .>= -1) .& (out .<= 1))
end

function scale(ins::Float64)
    2.0 * ins - 1.0
end

function scale(ins::Array{Float64})
    (2.0 .* ins) .- 1.0
end

function test_controller(cfg::Dict, c::Controller)
    cell_state() = scale(rand(cfg["n_cell_state"]))
    cell_params() = scale(rand(cfg["n_cell_params"]))
    conn_state() = scale(rand(cfg["n_conn_state"]))
    conn_params() = scale(rand(cfg["n_conn_params"]))

    @testset "Cell division" begin
        outs = c.cell_division(cell_params())
        @test typeof(outs) == Bool
        @test length(outs) == 1
    end

    @testset "New cell params" begin
        outs = c.new_cell_params(cell_params())
        @test typeof(outs) == Array{Float64,1}
        @test length(outs) == cfg["n_cell_params"]
        test_domain(outs)
    end

    @testset "Cell death" begin
        outs = c.cell_death(cell_params())
        @test typeof(outs) == Bool
        @test length(outs) == 1
    end

    @testset "Cell state update" begin
        outs = c.cell_state_update(vcat(cell_params(), cell_state(), rand()))
        @test typeof(outs) == Array{Float64,1}
        @test length(outs) == cfg["n_cell_state"] + 1
        test_domain(outs)
    end

    @testset "Cell param update" begin
        outs = c.cell_param_update(vcat(cell_params(), cell_state()))
        @test typeof(outs) == Array{Float64,1}
        @test length(outs) == cfg["n_cell_params"]
        test_domain(outs)
    end

    @testset "Connect" begin
        outs = c.connect(vcat(cell_params(), cell_params(), rand()))
        @test typeof(outs) == Bool
        @test length(outs) == 1
    end

    @testset "New conn params" begin
        outs = c.new_conn_params(vcat(cell_params(), cell_params(), rand()))
        @test typeof(outs) == Array{Float64,1}
        @test length(outs) == cfg["n_conn_params"]
        test_domain(outs)
    end

    @testset "Disconnect" begin
        outs = c.disconnect(
            vcat(cell_params(), cell_params(), conn_params(), conn_state()))
        @test typeof(outs) == Bool
        @test length(outs) == 1
    end

    @testset "Conn state update" begin
        outs = c.conn_state_update(vcat(conn_params(), conn_state(), rand(2)))
        @test typeof(outs) == Array{Float64,1}
        @test length(outs) == cfg["n_conn_state"] + 1
        test_domain(outs)
    end

    @testset "Conn param update" begin
        outs = c.conn_param_update(
            vcat(cell_params(), cell_params(), conn_params(), conn_state()))
        @test typeof(outs) == Array{Float64,1}
        @test length(outs) == cfg["n_conn_params"]
        test_domain(outs)
    end
end

@testset "Default controller" begin
    cfg = Config()
    c = Controller(cfg)
    test_controller(cfg, c)
end

@testset "Random controller" begin
    cfg = Config()
    c = rand_controller(cfg)
    test_controller(cfg, c)
end

@testset "Const controller" begin
    cfg = Config()
    c = const_controller(cfg)
    test_controller(cfg, c)
end

@testset "Static controller" begin
    cfg = Config(Config())
    c = static_controller(cfg)
    test_controller(cfg, c)
end

# @testset "SNN controller" begin
#     cfg = Config(Config(), "cfg/snn.yaml")
#     c = snn_controller(cfg)
#     test_controller(cfg, c)
# end

@testset "CGP controller" begin
    include("../evo/cgp.jl")
    cfg = Config()
    c = cgp_controller(cfg)
    test_controller(cfg, c)
end

nothing


text:
include("../src/SpinMonteCarlo.jl")

using SpinMonteCarlo

const S = 0.5
const J = 1.0
const Gamma = 0.0
const Ls = [8]
const Ts = 0.5:0.5:10.01
const MCS = 8192
const Therm = MCS >> 3

params = Dict{String, Any}[]
for L in Ls
    for T in Ts
        push!(params,
              Parameter("Model"=>QuantumXXZ,
                        "Lattice"=>chain_lattice,
                        "L"=>L, "T"=>T, "J"=>J, "S"=>S, "Gamma"=>Gamma,
                        "Update Method"=>loop_update!,
                        "MCS"=>MCS, "Therm"=>Therm,
                        "Verbose"=>true,
                       ))
    end
end

obs = runMC(params)

const pnames = ["S", "J", "Gamma", "L", "T"]
const onames = ["Magnetization",
                "|Magnetization|",
                "Magnetization^2",
                "Magnetization^4",
                "Binder Ratio",
                "Susceptibility",
                "Connected Susceptibility",
                "Energy",
                "MCS per Second",
                "Time per MCS",
               ]

const io = open("res-QuantumXXZ.dat", "w")
i=1
for pname in pnames
    println(io, "# \$$i : $pname")
    i+=1
end
for oname in onames
    println(io, "# \$$i, $(i+1): $oname")
    i+=2
end

for (p,o) in zip(params, obs)
    for pname in pnames
        print(io, p[pname], " ")
    end
    for oname in onames
        @printf(io, "%.15f %.15f ", mean(o[oname]), stderror(o[oname]))
    end
    println(io)
end


text:
#### naive sampling

immutable CategoricalDirectSampler <: Sampleable{Univariate,Discrete}
    prob::Vector{Float64}

    function CategoricalDirectSampler(p::Vector{Float64})
        isempty(p) && error("p is empty.")
        new(p)
    end
end
ncategories(s::CategoricalDirectSampler) = length(s.prob)

function rand(s::CategoricalDirectSampler)
    p = s.prob
    n = length(p)
    i = 1
    c = p[1]
    u = rand()
    while c < u && i < n
        c += p[i += 1]
    end
    return i
end

##### Alias Table #####

immutable AliasTable <: Sampleable{Univariate,Discrete}
    accept::Vector{Float64}
    alias::Vector{Int}
    isampler::RandIntSampler
end
ncategories(s::AliasTable) = length(s.accept)

function AliasTable{T<:Real}(probs::AbstractVector{T})
    n = length(probs)
    n > 0 || error("The input probability vector is empty.")
    accp = Array(Float64, n)
    alias = Array(Int, n)
    StatsBase.make_alias_table!(probs, 1.0, accp, alias)
    AliasTable(accp, alias, RandIntSampler(n))
end

rand(s::AliasTable) =
    (i = rand(s.isampler); u = rand(); @inbounds r = u < s.accept[i] ? i : s.alias[i]; r)

show(io::IO, s::AliasTable) = @printf(io, "AliasTable with %d entries", ncategories(s))


text:
using BirdClef
data_path= get(ENV,"BIRDCLEF_DATA_DIR",joinpath(@__DIR__,"data"))
main(data_path)

max_stars_repo_name: BirdClef/birdclef

max_issues_repo_name: BirdClef/birdclef

max_forks_repo_name: BirdClef/birdclef


text:
# ---
# title: 989. Add to Array-Form of Integer
# id: problem989
# author: Qling
# date: 2021-05-10
# difficulty: Easy
# categories: Array
# link: <https://leetcode.com/problems/add-to-array-form-of-integer/description/>
# hidden: true
# ---
#
# For a non-negative integer `X`, the  _array-form of`X`_ is an array of its
# digits in left to right order.  For example, if `X = 1231`, then the array
# form is `[1,2,3,1]`.
#
# Given the array-form `A` of a non-negative integer `X`, return the array-form
# of the integer `X+K`.
#
#
#
# **Example 1:**
#
#
#
#     Input: A = [1,2,0,0], K = 34
#     Output: [1,2,3,4]
#     Explanation: 1200 + 34 = 1234
#
#
# **Example 2:**
#
#
#
#     Input: A = [2,7,4], K = 181
#     Output: [4,5,5]
#     Explanation: 274 + 181 = 455
#
#
# **Example 3:**
#
#
#
#     Input: A = [2,1,5], K = 806
#     Output: [1,0,2,1]
#     Explanation: 215 + 806 = 1021
#
#
# **Example 4:**
#
#
#
#     Input: A = [9,9,9,9,9,9,9,9,9,9], K = 1
#     Output: [1,0,0,0,0,0,0,0,0,0,0]
#     Explanation: 9999999999 + 1 = 10000000000
#
#
#
#
# **Note：**
#
#   1. `1 <= A.length <= 10000`
#   2. `0 <= A[i] <= 9`
#   3. `0 <= K <= 10000`
#   4. If `A.length > 1`, then `A[0] != 0`
#
#
## @lc code=start
using LeetCode

function add_to_array_form(num::Vector{Int}, k::Int)::Vector{Int}
    n = length(num)
    i, val, carry = n, 0, 0
    res = Int[]

    while i > 0 || k != 0
        x = i > 0 ? num[i] : 0
        y = k != 0 ? k % 10 : 0

        val = x + y + carry
        append!(res, val % 10)
        carry = val ÷ 10
        i, k = i - 1, k ÷ 10
    end

    if carry > 0
        append!(res, carry)
    end

    return reverse(res)
end
## @lc code=end


text:
function OMPTimeDomain(V::Array{Complex128,2},Precision::Int64,OrthoTol::Float64,gmethod,lmethod,ltol,gWallClockTime)
	set_bigfloat_precision(Precision)
	#Add a zero at the beginning of the trace to put make it
	#compatible with a strictly proper transfer function
	#which is what we are using to fit the data
	V = cat(1,zeros(Complex128,(1,size(V,2))),V)
	L = size(V,1)
	N = size(V,1)
	zarb = genZArb(L)
	zdbl = convert(Vector{Complex128},zarb)
	#pole vectors one at arbitrary precsion, one at dbl
	parb = zeros(ComplexPair{BigFloat},N-1)
	Larb = zeros(ComplexPair{BigFloat},N-1)
	Ldbl = zeros(ComplexPair{Float64},N-1)
	pdbl = zeros(Complex128,N-1)
	AICc = zeros(Float64,N-1)
    RSS = zeros(Float64,N-1)
	B = ones(ComplexPair{BigFloat},(L,N))/sqrt(L)
	Bdbl = ones(Complex128,(L,N))/sqrt(L)
	#Set up some fitting parameters for NLOPT
	opt = Opt(gmethod,2)
    lopt = Opt(lmethod,2)
    ftol_rel!(lopt, ltol) #local stopping criteria
    maxtime!(opt,gWallClockTime) #global stopping critera
    betalb = [0.01; 0]
    betaub = [0.99; (2*pi-eps(Float64))]
    lower_bounds!(opt,betalb)
    upper_bounds!(opt,betaub)
	terminate_fit = 0;
	k = 0;
	Resid = copy(V)
	while terminate_fit<1
		k += 1
		if k==1
			anm1 = zeros(Complex128,0)
		else
			anm1 = [pdbl[k-1]]
		end
		Bnm1 = convert(Vector{Complex128},B[:,k])
		max_objective!(opt,(x,g)->ompObj(x,g,zdbl,anm1,Bnm1,Resid))
		beta0 = genStart(betalb,betaub)
		(oVal,betaOpt,exitcode) = optimize(opt,beta0)
		println(oVal)
		parb[k] = convert(ComplexPair{BigFloat},betaOpt[1]*exp(im*betaOpt[2]))
		pdbl[k] = convert(Complex128,parb[k])
		(Bnext,Ldbl[k]) = phiNext(zdbl,[pdbl[1:k]],Bnm1)
		(B[:,k+1],Larb[k]) = phiNext(zarb,[parb[1:k]],B[:,k])
		Bdbl[:,k+1] = convert(Vector{Complex128},B[:,k+1])
		OrthogonalityMetric = maximum(abs(abs(Bnext'*Bdbl[:,1:k+1]) - cat(1,zeros(Float64,k),one(Float64))'))
		println(OrthogonalityMetric)
		if OrthogonalityMetric>OrthoTol
			#terminate the fit if we lose orthogonality when operating at double precision
			terminate_fit = 2;
		end
		Tdbl = ifft(Bdbl[:,1:k+1],1)*sqrt(L)
		Resid = V-Tdbl*Tdbl'*V
		RSS[k] = norm(Resid)
		#We have L-1 actual data points per column x size(V,2) columns, 2 numbers per element
		#And we have 2 parameters per pole, +1 for the variance estimate, plus 2 per linear coefficient
		AICc[k] = akaikeInformationCriterion(RSS[k],2*k+1+2*k*size(V,2),2*size(V,2)*(size(V,1)-1))
		if k==(N-1)
			terminate_fit = 1
		end
	end
    Tdbl = ifft(Bdbl[:,1:k+1],1)*sqrt(L)
    A = Tdbl[:,1:(k+1)]'*V
	return (zarb,parb,Larb,A,AICc,terminate_fit)
end

function OMPTimeDomainIndPoles(V::Array{Complex128,2},Precision::Int64,OrthoTol::Float64,gmethod,lmethod,ltol,gWallClockTime)
	for j=1:size(V,2)
		OMPTimeDomain(V,Precision,OrthoTol,gmethod,lmethod,ltol,gWallClockTime)
	end
end

function akaikeInformationCriterion(RSS,numParams,dataLen)
    #This gives the sample size corrected Akaike Information Criteria
    #commonly known as AICc in the literature
    K = numParams+1
    return dataLen*log(RSS/dataLen)+2*K+2*K*((K+1)/(dataLen-K-1))
end

function genStart(lb,ub)
	bstart = zero(lb)
	for k=1:length(lb)
		bstart[k] = (ub[k]-lb[k])*rand(1)[1]+lb[k]
	end
	return bstart
end

function ompObj(beta::Vector,grad::Vector,z::Vector{Complex128},anm1::Vector{Complex128},Bnm1::Vector{Complex128},D::Array{Complex128,2})
    R = beta[1];
    T = beta[2];
    p = [anm1; R*exp(im*T)]
    (Bn,L) = phiNext(z,p,Bnm1)
    Tn = ifft(Bn,1)*sqrt(length(z))
    Proj = Tn'*D;
    return sum(abs(Proj),2)[1]
end

function createFreqAxis(N::Int,dT)
	#For some time data with index values 0:(N-1) this function creates the corresponding frequency axis
	#in units of 1/time spacing
    fs = 1/dT
    df = fs/N
    f = [0:df:(fs-df)] - (fs-mod(N,2)*df)/2
    return f
end


text:
# TODO: compactify Bools to flags
# TODO: composite isbits conversion
# TODO: aggregates isbits conversion (Maybe just NTuple{N, UInt64})

import ConstructionBase
using Base.Meta: isexpr

"""
    subtypes_fun(::Val{T}) where T

return the conceptual subtypes in `Symbol`s.
"""
function subtypes_fun end

"""
    isa_type_fun(::Val{T}, Val{S}, x)::Bool where {T, S}

It checks if `x` is a `S <: T`, where `S` is a symbol and `T` is a type.
"""
function isa_type_fun end

"""
    rt_constructor(x)::Function where {T, S}

Return the runtime constructor of `x`
"""
function rt_constructor end

"""
    @compactify [show_methods=true] exprs
"""
macro compactify(exs...)
    kws = []
    arg = exs[end]
    for i in 1:length(exs)-1
        x = exs[i]
        if x isa Expr && x.head === :(=) # Keyword given of the form "foo=bar"
            if length(x.args) != 2
                return Expr(:call, :error, "Invalid keyword argument: $x")
            end
            push!(kws, x.args[1] => x.args[2])
        else
            return Expr(:call, :error, "@$fcn expects only one non-keyword argument")
        end
    end
    _compactify(__module__, arg; debug=false, kws...)
end

function getname(T)
    if T isa Symbol
        return T
    elseif T isa Expr
        return T.head == :curly ? T.args[1] : error("$T as type name is not supported")
    end
end
params(T) = isexpr(T, :curly) ? T.args[2:end] : ()

function _compactify(mod, block; debug=false, show_methods=true)
    isexpr(block, :block) || error("@compatify takes a block!")
    stmts = block.args
    hasabstract = false
    abstract2concrete = Dict() # abstract type name => ismutable, struct body, [(name, concete fields), ...]
    names = []
    # notation: T means abstract name, S means concrete name
    for ex in stmts; ex isa LineNumberNode && continue
        if isexpr(ex, :macrocall) && ex.args[1] === Symbol("@abstract")
            hasabstract = true
            struct_body = ex.args[3]
            @assert isexpr(struct_body, :struct)
            ismutable, T, fields = struct_body.args
            if isexpr(T, :<:) # if there's a super type
                T = T.args[1]
            end
            T = getname(T)
            T in names && error("$T struct is already defined")

            field2val = parse_default_value!(struct_body)
            push!(names, T)
            abstract2concrete[T] = ismutable, struct_body, field2val, []
        elseif isexpr(ex, :struct)
            struct_body = ex
            field2val = parse_default_value!(struct_body)
            ismutable, S, fields = struct_body.args
            isexpr(S, :(<:)) || error("$S must be a subtype of some @abstract type!")
            S_with_params, T_with_params = S.args
            S = getname(S_with_params)
            T = getname(T_with_params)
            S in names && error("$S struct is already defined")

            T in keys(abstract2concrete) || error("$T >: $S is not a @abstract type.")
            ismut, = abstract2concrete[T]
            ismut == ismutable || error("$S and $T should have the same mutability!")
            push!(names, S)
            fields = filter(x->!(x isa LineNumberNode), fields.args)
            # destructs fields of the form `a::T` to field name and type.
            fields = [isexpr(f, :(::)) ? (f.args[1], expr_to_type(mod, f.args[2])) : (f, Any) for f in fields]
            push!(abstract2concrete[T][end], (S, fields, field2val, S_with_params, T_with_params))
        else
            error("What is this? $ex")
        end
    end

    hasabstract || error("There must be at least one `@abstract`!")

    expr = Expr(:block)

    gensymidx = Ref(0)
    gensym = let gensymidx = gensymidx
        x -> Symbol("###$x###", (gensymidx[] += 1))
    end

    for (T, (ismutable, struct_body, T_field2val, Ss)) in pairs(abstract2concrete)
        EnumType = gensym(T)
        EnumNumType = Int32
        # S1=0 S2=1 ...
        # We need some very weird names so there won't be a name collision.
        # @enum doesn't work with var"##xy".
        enum_expr = :(@enum $EnumType::$EnumNumType)
        for i in 1:length(Ss)
            push!(enum_expr.args, :($(Symbol(:₋₃₋₁₂₉, T, :₋__, Ss[i][1], :₋₃₋₁₉₉₂₋₋)) = $(EnumNumType(i-1))))
        end
        push!(expr.args, enum_expr)

        push!(expr.args, struct_body)
        if debug
            @info "Parsed:"
            Base.print_array(stdout, Ss); println()
        end

        S2enum_num = Dict{Symbol, EnumNumType}(S => EnumNumType(i-1) for (i, (S, )) in enumerate(Ss))

        # S: struct name | f: field name | t: field type
        #
        # For non-isbits types, we just make them `Any` and type assert it in
        # the `getproperty` function.
        #
        # First: we initially have something like
        #   sz\S   S1       S2       S3       S4
        #  sizeof  8        8        1        1
        #  sizeof  8        Any      1        1
        #  sizeof  4        1        4        1
        #  sizeof  1        2        8        4
        #  sizeof  2        Any               Any
        #
        # Here, we sort the isbits' field size (the first dim) and separate
        # non-isbits fields. We then have:
        #   sz\S   S1       S2       S3       S4
        #  sizeof  8        8        8        4
        #  sizeof  8        2        4        1
        #  sizeof  4        1        1        1
        #  sizeof  2                 1        1
        #  sizeof  1
        #  Anys             x                 x
        #  Anys             x
        isbits_S_ft = []
        nonisbits_S_ft = []
        max_num_isbits = max_num_nonisbits = 0
        for S in Ss
            S, fts, S_field2val = S
            isbits_ft = []
            nonisbits_ft = []
            push!(isbits_S_ft, S => isbits_ft)
            push!(nonisbits_S_ft, S => nonisbits_ft)
            for (f, t) in fts
                if isbitstype(t)
                    push!(isbits_ft, (f, t))
                else
                    push!(nonisbits_ft, (f, t))
                end
            end
            max_num_isbits = max(max_num_isbits, length(isbits_ft))
            max_num_nonisbits = max(max_num_nonisbits, length(nonisbits_ft))
            sort!(isbits_ft, by=sizeof ∘ last, rev=true)
        end

        # common fields in the abstract type T
        common_fields = [
                         isexpr(f, :(::)) ? f.args[1] : f
                         for f in filter(x->!(x isa LineNumberNode), struct_body.args[3].args)
                        ]

        # This structure is needed to generate `getproperty` to simulate the
        # `Ss` types.
        S2fields = Dict{
                        Symbol, # S: concrete type's name. We need this to later build getproperty
                                          # :a => (:b, Any, Complex)
                        Dict{Symbol,Any}, # oldname => (newname, newtype, oldtype)
                       }()
        # This structure is needed to generate constructors.
        back_edges = Dict{Symbol,Any}() # newname => [(S, oldname, oldtype), ...]
        compact_fields = Symbol[]
        compact_types = []

        # We then just need to read-off the sorted fields by taking the max
        # horizonally.
        #
        #   sz\S   S1       S2       S3       S4    [compactified]
        #  sizeof  8        8        8        4           8
        #  sizeof  8        2        4        1           8
        #  sizeof  4        1        1        1           4
        #  sizeof  2                 1        1           2
        #  sizeof  1                                      1
        #  Anys             x                 x           x
        #  Anys             x                             x
        for i in 1:max_num_isbits
            siz = idx = 0
            for (j, (S, fts)) in enumerate(isbits_S_ft); i > length(fts) && continue
                f, t = fts[i]
                current_siz = sizeof(t)
                if siz < current_siz
                    siz = current_siz
                    idx = j
                end
            end
            S, fts = isbits_S_ft[idx]
            f, newtype = fts[i]
            newname = gensym(f)
            edgs = back_edges[newname] = []
            for (S, fts) in isbits_S_ft; i > length(fts) && continue
                f, t = fts[i]
                namemap = get!(() -> Dict{Symbol,Any}(), S2fields, S)
                namemap[f] = (newname, newtype, t)
                push!(edgs, (S, f, t))
            end
            push!(compact_fields, newname)
            push!(compact_types, newtype)
            push!(struct_body.args[end].args, :($newname::$newtype))
        end
        for i in 1:max_num_nonisbits
            newtype = Any
            newname = gensym("Any")
            edgs = back_edges[newname] = []
            for (S, fts) in nonisbits_S_ft; i > length(fts) && continue
                f, t = fts[i]
                namemap = get!(() -> Dict{Symbol,Any}(), S2fields, S)
                namemap[f] = (newname, newtype, t)
                push!(edgs, (S, f, t))
            end
            push!(compact_fields, newname)
            push!(compact_types, newtype)
            push!(struct_body.args[end].args, :($newname::$newtype))
        end

        debug && @info "" isbits_S_ft nonisbits_S_ft

        tagname = gensym("tag")
        tagname_q = Meta.quot(tagname)
        push!(struct_body.args[end].args, :($tagname::$EnumType))

        # build getproperties
        propnames = :(function (::$(typeof(ConstructionBase.getproperties)))(x::$T) end)
        push!(expr.args, propnames)
        body = propnames.args[end].args
        push!(body, Expr(:meta, :inline))
        ifold = expr
        for (S, namemap) in pairs(S2fields)
            # if we are simulating for type `S`.
            error_message = :($unreachable())
            condition = :($reinterpret($EnumNumType, $getfield(x, $tagname_q)) === $(S2enum_num[S]))
            uninitialized = expr === ifold
            names = vcat(common_fields, map(first, collect(pairs(namemap))))
            tup = Expr(:tuple)
            for n in names
                push!(tup.args, :($getproperty(x, $(QuoteNode(n)))))
            end
            nt = :($NamedTuple{$((names...,))}($tup))
            ifnew = Expr(ifelse(uninitialized, :if, :elseif), condition, nt)
            uninitialized ? push!(body, ifnew) : push!(ifold.args, ifnew)
            ifold = ifnew
        end

        # build rt_constructor
        rt_fun = :(function (::$(typeof(rt_constructor)))(x::$T) end)
        push!(expr.args, rt_fun)
        body = rt_fun.args[end].args
        push!(body, Expr(:meta, :inline))
        ifold = expr
        for (S, namemap) in pairs(S2fields)
            # if we are simulating for type `S`.
            error_message = :($unreachable())
            condition = :($reinterpret($EnumNumType, $getfield(x, $tagname_q)) === $(S2enum_num[S]))
            uninitialized = expr === ifold
            ifnew = Expr(ifelse(uninitialized, :if, :elseif), condition, S)
            uninitialized ? push!(body, ifnew) : push!(ifold.args, ifnew)
            ifold = ifnew
        end

        # build getproperty
        getprop = :(function (::$(typeof(getproperty)))(x::$T, s::$Symbol) end)
        push!(expr.args, getprop)
        body = getprop.args[end].args
        push!(body, Expr(:meta, :inline))
        ifold = expr
        ifold_original = ifold
        for cs in common_fields
            uninitialized = expr === ifold
            qncs = QuoteNode(cs)
            behavior = Expr(:call, getfield, :x, qncs)
            ifnew = Expr(ifelse(uninitialized, :if, :elseif), :(s === $qncs), behavior)
            uninitialized ? push!(body, ifnew) : push!(ifold.args, ifnew)
            uninitialized && (ifold_original = ifnew)
            ifold = ifnew
        end
        for (S, namemap) in pairs(S2fields)
            # if we are simulating for type `S`.
            behavior = expr
            behavior_og = behavior
            for (oldname, (newname, newtype, oldtype)) in pairs(namemap)
                uninitialized = expr === behavior
                newf = QuoteNode(newname)
                oldf = QuoteNode(oldname)
                condition = :(s === $oldf)
                behavior′ = :($getfield(x, $newf))
                if newtype === Any
                    # type assert doesn't really work here, because
                    # (1+1im)::Complex{Real} errors. Alternatively, we could add
                    # `convert(Complex{Real}, 1+1im)` but this is even worse as
                    # it allocates. So for the time being, we will not check the
                    # type and just return.
                    behavior′ = :($behavior′)
                else
                    @assert isbitstype(oldtype) && isbitstype(newtype)
                    behavior′ = :($reconstruct($oldtype, $behavior′)::$oldtype)
                end
                if uninitialized
                    behavior_og = behavior = Expr(:if, condition, behavior′)
                else
                    ifnew = Expr(:elseif, condition, behavior′)
                    push!(behavior.args, ifnew)
                    behavior = ifnew
                end
            end
            error_message = :($throw_no_field($(Val(S)), s))
            condition = :($reinterpret($EnumNumType, $getfield(x, $tagname_q)) === $(S2enum_num[S]))
            uninitialized = expr === ifold
            if behavior === expr
                behavior_og = error_message
            else
                push!(behavior.args, error_message)
            end
            ifnew = Expr(ifelse(uninitialized, :if, :elseif), condition, behavior_og)
            uninitialized ? push!(body, ifnew) : push!(ifold.args, ifnew)
            ifold = ifnew
        end
        @assert expr !== ifold "no getproperty matches?"

        # Let's now make the constructor
        construct_args = []
        append!(construct_args, common_fields)
        append!(construct_args, compact_fields)

        for (S, fts, S_field2val, S_with_params, T_with_params) in Ss
            parameters = Expr(:parameters)
            for f in common_fields
                push!(parameters.args, Expr(:kw, f, T_field2val[f]))
            end
            for (f, _) in fts
                push!(parameters.args, Expr(:kw, f, S_field2val[f]))
            end

            constructor = quote
                struct $S_with_params 1+1 end
                function (::Type{$S_with_params})($parameters) where {$(params(S_with_params)...)} end
            end
            constructor_body = constructor.args[end].args[end].args

            # We check if the compact field is native in the struct S. If it is
            # native, then we only need to translate oldname to the newname. If
            # it is not native, then we need to set it to some default value.
            for (newname, newtype) in zip(compact_fields, compact_types)
                isany = newtype === Any
                edgs = back_edges[newname]
                is_native = false
                for (S2, oldname, oldtype) in edgs
                    if S2 === S
                        is_native = true
                        if isany # don't call simulate_type for Anys
                            push!(constructor_body, :($newname = $oldname))
                        else
                            push!(constructor_body, :($newname = $simulate_type($newtype, $oldname)))
                        end
                        break
                    end
                end
                is_native && continue
                S2, oldname, oldtype = first(edgs)
                idx = findfirst(x->x[1] === S2, Ss)
                defval = Ss[idx][3][oldname]
                if isany # don't call simulate_type for Anys
                    push!(constructor_body, :($newname = $defval))
                else
                    push!(constructor_body, :($newname = $simulate_type($newtype, $defval)))
                end
            end

            # call the real constructor.
            construct_expr = Expr(:call, T_with_params)
            append!(construct_expr.args, construct_args)
            # the type tag is the last arg
            push!(construct_expr.args, Expr(:call, reinterpret, EnumType, S2enum_num[S]))
            push!(constructor_body, construct_expr)
            push!(expr.args, constructor)
        end

        if show_methods
            # Let's do pretty print
            pretty_print = :(function (::$(typeof(Base.show)))(io::$IO, ::$(MIME"text/plain"), obj::$T) end)
            body = pretty_print.args[end].args
            ifold = expr
            for (S, fts, S_field2val) in Ss
                uninitialized = ifold === expr
                enum = Expr(:call, reinterpret, EnumType, S2enum_num[S])
                condition = :($enum === $getfield(obj, $tagname_q))
                behavior = Expr(:call, print, :io, Meta.quot(S), "(")
                n = length(fts)
                for (i, (f, _)) in enumerate(fts)
                    f = Meta.quot(f)
                    push!(behavior.args, f)
                    push!(behavior.args, " = ")
                    push!(behavior.args, :($getproperty(obj, $f)))
                    i == n || push!(behavior.args, ", ")
                end
                push!(behavior.args, ")::")
                push!(behavior.args, T)
                ifnew = Expr(ifelse(uninitialized, :if, :elseif), condition, behavior)
                uninitialized ? push!(body, ifnew) : push!(ifold.args, ifnew)
                ifold = ifnew
            end
            push!(expr.args, pretty_print)
        end

        # Let's generate `subtypes`-like function.
        subtypes_fun_expr = :((::$(typeof(subtypes_fun)))(::$Val{T}) where {T<:$T} = $([x[1] for x in Ss]...,))
        push!(expr.args, subtypes_fun_expr)
        # Let's generate `isa S` checking functions
        for (S, fts, S_field2val) in Ss
            isa_fun = :(
                        function (::$(typeof(isa_type_fun)))(::$Val{T}, ::$Val{$(Meta.quot(S))}, x) where {T<:$T}
                            $reinterpret($EnumNumType, $getfield(x, $tagname_q)) === $(S2enum_num[S])
                        end
                       )
            push!(expr.args, isa_fun)
        end
    end
    expr = esc(expr)
    debug && print(expr)
    expr
end

@generated function reconstruct(::Type{T}, x::S) where {T,S}
    @assert isbitstype(T)
    @assert sizeof(T) ≤ sizeof(S)
    if sizeof(T) == 0
        return T.instance
    elseif sizeof(T) == sizeof(S)
        return :(reinterpret($T, x))
    else
        IS = Symbol(:UInt, 8*sizeof(S))
        IT = Symbol(:UInt, 8*sizeof(T))
        return :(reinterpret($T, reinterpret($IS, x) % $IT))
    end
end

@generated function simulate_type(::Type{T}, x::S) where {T,S}
    @assert isbitstype(T)
    @assert sizeof(T) ≥ sizeof(S)
    if sizeof(T) == 0
        return T.instance
    elseif sizeof(T) == sizeof(S)
        return :(reinterpret($T, x))
    else
        IS = Symbol(:UInt, 8*sizeof(S))
        IT = Symbol(:UInt, 8*sizeof(T))
        return :(reinterpret($T, $IT(reinterpret($IS, x))))
    end
end

@noinline unreachable() = error("unreachable reached")
@noinline throw_no_field(::Val{S}, s) where {S} = error("type $S has no field $s.")
@noinline throw_type_error(::Val{S}, s) where {S} = error("Expected $S, got $s::$(typeof(s)).")

function parse_default_value!(expr::Expr)
    @assert isexpr(expr, :struct)
    field2val = Dict{Symbol, Any}()
    body = expr.args[end].args
    for (i, ex) in enumerate(body); ex isa LineNumberNode && continue
        isexpr(ex, :(=)) || error("$ex doesn't have a default value!")
        ft, val = ex.args
        body[i] = ft
        f = isexpr(ft, :(::)) ? ft.args[1] : ft
        field2val[f] = val
    end
    field2val
end


text:
export get_filter_params, FilterBank1d, get_FilterBanks, FilterBank1dBlock

"""
    ξmax(Q)

Calculate the maximum central frequency of a Morlet wavelet.

# Arguments
 - `Q::Int64`: Number of wavelets per octave.
"""
function ξmax(Q::Int64)
    return max(1/(1+2^(3/Q)), 0.35)
end

"""
    σmin(σ0, J)

Calculate the bandwidth of the lowpass filter ϕ.

# Arguments
 - `σ0::Float64`: Bandwidth of lowpass filter ϕ at scale 0.
 - `J::Int64`: Maximum wavelet scale.
"""
function σmin(σ0::Float64, J::Int64)
    return σ0/2^J
end

"""
    σmax(ξmax, Q, rψ)

Calculate the bandwidth of a Morlet wavelet in a family of Q per octave.

# Arguments
 - `ξmax::Float64`: Maximum central frequency.
 - `Q::Int64`: Number of wavelets per octave.
 - `rψ::Float64`: Parameter controlling frequency width of filters.
"""
function σmax(ξmax::Float64, Q::Int64, rψ::Float64)
    return ξmax * ((1-2^(-1/Q))/(1+2^(-1/Q))) / (2*log(1/rψ))^(1/2)
end

"""
    max_dyadic_subsampling(ξ, σ, α)

Compute the maximal subsampling possible for a Gabor wavelet without aliasing.

# Arguments
 - `ξ::Float64`: Central frequency of the wavelet
 - `σ::Float64`: Bandwidth of the wavelet
 - `α::Float64`: Parameter controlling allowable error, with error∝1/α
"""
function max_dyadic_subsampling(ξ::Float64, σ::Float64, α::Float64)
    upper_bound = min(ξ+α*σ, 0.5)
    j = convert(Int, floor(-log2(upper_bound))) - 1
    return j
end

"""
    n_dyadic_steps(Q, σ_min, σ_max)

Get the number of dyadic steps between σ_min and σ_max given Q wavelets per octave.

# Arguments
 - `Q::Int64`: Number of wavelets per octave.
 - `σ_min::Float64`: Bandwidth of the lowpass filter
 - `σ_max::Float64`: Maximum filter bandwidth
"""
function n_dyadic_steps(Q::Int64, σ_min::Float64, σ_max::Float64)
    return convert(Int, ceil(Q*log2(σ_max/σ_min)))
end

"""
    get_filter_params(Q, J, σ0, rψ=rψ, α=α)

Get the central frequencies, bandwidths, and maximum subsampling for a family of Morlet wavelets.

# Arguments
 - `Q::Int64`: Number of wavelets per octave.
 - `J::Int64`: Maximum scale of the filters.
 - `σ_min::Float64`: Bandwidth of the low-pass filter.
 - `rψ`::Float64: Parameter controlling frequency width of filters.
 - `α::Float64`: Parameter controlling allowable error, with error∝1/α
"""
function get_filter_params(Q::Int64, σ_min::Float64; rψ::Float64=sqrt(0.5), α::Float64=5.0)
    @assert Q>=1 "At least one wavelet per octave is required."

    # maximum central frequency
    ξ_max = ξmax(Q)

    # maximum filter bandwidth
    σ_max = σmax(ξ_max, Q, rψ)

    # number of dyadic steps between σ_min and σ_max
    n = n_dyadic_steps(Q, σ_min, σ_max)

    # initialize vectors
    σ = collect(0.0:n+Q-2)
    ξ = collect(0.0:n+Q-2)
    j = Vector{Int64}(undef, n+Q-1)

    # σ_i = σ_max/2^(i/Q)
    σ[1:n] .= σ_max./2.0.^(σ[1:n]./Q)
    σ[n+1:end] .= σ_min
    # ξ_i = ξ_max/2^(i/Q)
    ξ[1:n] .= ξ_max./2.0.^(ξ[1:n]./Q)
    ξ[n+1:end] .= ξ[n].*((Q.-collect(1:Q-1))./Q)

    # js
    j .= max_dyadic_subsampling.(ξ, σ, α)

    return σ, ξ, j
end

"""
    get_filter_params(Q, J, σ0, rψ=rψ, α=α)

Get the central frequencies, bandwidths, and maximum subsampling for families of
Morlet wavelets for each layer of the scattering transform.

# Arguments
- `Q::Array{Int64}`: Number of wavelets per octave for each layer.
- `J::Int64`: Maximum scale of the filters.
- `σ0::Float64`: Bandwidth of the low-pass filter.
- `rψ`::Float64: Parameter controlling frequency width of filters.
- `α::Float64`: Parameter controlling allowable error, with error∝1/α
"""
function get_filter_params(Q::Vector{Int64}, J::Int64, σ0::Float64; rψ::Float64=sqrt(0.5), α::Float64=5.0)
    # minimum filter bandwidth
    σ_min = σ0/2^J

    # initialize arrays for the parameters for each layer
    ξs = Vector{Vector{Float64}}(undef, length(Q))
    σs = Vector{Vector{Float64}}(undef, length(Q))
    js = Vector{Vector{Float64}}(undef, length(Q))

    # get filter parameters for each layer
    for i=1:length(Q)
        σs[i], ξs[i], js[i] = get_filter_params(Q[i], σ_min, rψ=rψ, α=α)
    end

    return σs, ξs, js
end

abstract type AbstractFilterBank end

"""
    FilterBank1d

A structure for collections of 1D Morlet Wavelets and a Gaussian lowpass filter

## Fields: FilterBank1d
 | **Field** | **Description** |
 |:----------|:----------------|
 | :N        | log2(signal length) |
 | :J        | log2(max filter scale) |
 | :Q        | Number of wavelets per octave |
 | :ω        | Frequencies |
 | :ϕ        | Fourier transformed low-pass Gaussian filter |
 | :Λ        | Fourier transformed Morlet wavelets |
"""
struct FilterBank1d <: AbstractFilterBank
    N::Int # log2 length of the signal in the time domain
    J::Int # log2 maximum scale of filters
    Q::Int # Number of wavelets per octave

    ω::Vector{Float64} # frequencies
    ϕ::GaussianFilter # low-pass filter
    Λ::Vector{MorletWavelet} # Fourier-transformed Morlet wavelets stored in a vector

    function FilterBank1d(N::Int64, Q::Int64, J::Int64, σ0::Float64; rψ::Float64=sqrt(0.5), α::Float64=5.0)
        σs, ξs, js = get_filter_params(Q, σ0/2^J, rψ=rψ, α=α)
        ω = FFTW.fftfreq(2^N)

        Λ = Vector{MorletWavelet}(undef, length(σs))
        ϕ = GaussianFilter(2^N, σ0/2^J)
        for i=1:length(σs)
            γ = GaborWavelet(2^N, σs[i], ξs[i], js[i])
            Λ[i] = MorletWavelet(γ, ϕ)
        end

        return new(N, J, Q, ω, ϕ, Λ)
    end
end

"""
    FilterBank1dBlock

A structure for collections of 1D Morlet Wavelets and a Gaussian lowpass filter,
both stored as arrays.

## Fields: FilterBank1d
 | **Field** | **Description** |
 |:----------|:----------------|
 | :N        | log2(signal length) |
 | :J        | log2(max filter scale) |
 | :Q        | Number of wavelets per octave |
 | :σs       | Bandwidth of filters |
 | :ξs       | Central frequencies of filters |
 | :js       | Maximum acceptable subsampling |
 | :σg       | Bandwidth of gaussian lowpass filter |
 | :ω        | Frequencies |
 | :ϕ        | Fourier transformed low-pass Gaussian filter |
 | :Λ        | Fourier transformed Morlet wavelets |
"""
struct FilterBank1dBlock <: AbstractFilterBank
    N::Int64 # log2 length of filter in time domain
    J::Int # log2 maximum scale of filters
    Q::Int # Number of wavelets per octave

    # these would have been in the MorletWavelet objects
    σs::Vector{Float64}
    ξs::Vector{Float64}
    js::Vector{Int64}
    # this would have been in the GaussianFilter object
    σg::Float64

    ω::Vector{Float64} # frequencies
    ϕ::Vector{Float64} # low-pass filter
    Λ::Array{Float64,2} # Fourier-transformed Morlet wavelets stored in a block

    function FilterBank1dBlock(N::Int64, Q::Int64, J::Int64, σ0::Float64; rψ::Float64=sqrt(0.5), α::Float64=5.0)
        σs, ξs, js = get_filter_params(Q, σ0/2^J, rψ=rψ, α=α)
        ω = FFTW.fftfreq(2^N)

        Λ = Array{Float64,2}(undef, 2^N, length(σs))
        ϕ = GaussianFilter(2^N, σ0/2^J)
        for i=1:length(σs)
            γ = GaborWavelet(2^N, σs[i], ξs[i], js[i])
            Λ[:, i] = MorletWavelet(γ, ϕ).ψ
        end

        return new(N, J, Q, σs, ξs, js, ϕ.σ, ω, ϕ.ϕ, Λ)
    end

    function FilterBank1dBlock(F::FilterBank1d)
        Λ = Array{Float64,2}(undef, 2^F.N, length(F.Λ))
        σs = Vector{Float64}(undef, length(F.Λ))
        ξs = Vector{Float64}(undef, length(F.Λ))
        js = Vector{Int64}(undef, length(F.Λ))
        for i=1:length(F.Λ)
            Λ[:,i].=F.Λ[i].ψ
            σs[i]=F.Λ[i].σ
            ξs[i]=F.Λ[i].ξ
            js[i]=F.Λ[i].j
        end
        return new(F.N, F.J, F.Q, σs, ξs, js, F.ω, F.ϕ.ϕ, Λ)
    end
end

"""
    get_FilterBanks(N, Q, J, σ0, rψ=rψ, α=α)

Get filter banks for every layer.

# Arguments
- `N::Int64`: log2(signal length)
- `Q::Vector{Int64}`: Number of wavelets per octave for each layer.
- `J::Vector{Int64}`: Maximum scale of the filters for each layer
- `σ0::Vector{Float64}`: Bandwidth of the low-pass filter for each layer
- `rψ`::Float64: Parameter controlling frequency width of filters.
- `α::Float64`: Parameter controlling allowable error, with error∝1/α
- `typ::Any`: FilterBank1d or FilterBank1dBlock
"""
function get_FilterBanks(N::Int64, Q::Vector{Int64}, J::Vector{Int64}, σ0::Vector{Float64}; rψ::Float64=sqrt(0.5), α::Float64=5.0, typ::Type{T}=FilterBank1d) where T<:AbstractFilterBank
    FilterBanks=Vector{typ}(undef, length(Q))
    for i=1:length(Q)
        if typ==FilterBank1d
            FilterBanks[i] = FilterBank1d(N, Q[i], J[i], σ0[i], rψ=rψ, α=α)
        elseif typ==FilterBank1dBlock
            FilterBanks[i] = FilterBank1dBlock(N, Q[i], J[i], σ0[i], rψ=rψ, α=α)
        end
    end
    return FilterBanks
end
"""
    get_FilterBanks(N, Q, J, σ0, rψ=rψ, α=α)

Convenience method fo getting filter banks for each layer with common σ0.
"""
function get_FilterBanks(N::Int64, Q::Vector{Int64}, J::Vector{Int64}, σ0::Float64; rψ::Float64=sqrt(0.5), α::Float64=5.0, typ::Type{T}=FilterBank1d) where T<:AbstractFilterBank
    σ0s = repeat([σ0], length(Q))
    get_FilterBanks(N, Q, J, σ0s; rψ=rψ, α=α, typ=typ)
end
"""
    get_FilterBanks(N, Q, J, σ0, rψ=rψ, α=α)

Convenience method fo getting filter banks for each layer with common J.
"""
function get_FilterBanks(N::Int64, Q::Vector{Int64}, J::Int64, σ0::Vector{Float64}; rψ::Float64=sqrt(0.5), α::Float64=5.0, typ::Type{T}=FilterBank1d) where T<:AbstractFilterBank
    Js = repeat([J], length(Q))
    get_FilterBanks(N, Q, Js, σ0; rψ=rψ, α=α, typ=typ)
end
"""
    get_FilterBanks(N, Q, J, σ0, rψ=rψ, α=α)

Convenience method fo getting filter banks for each layer with common σ0 and J.
"""
function get_FilterBanks(N::Int64, Q::Vector{Int64}, J::Int64, σ0::Float64; rψ::Float64=sqrt(0.5), α::Float64=5.0, typ::Type{T}=FilterBank1d) where T<:AbstractFilterBank
    σ0s = repeat([σ0], length(Q))
    Js = repeat([J], length(Q))
    get_FilterBanks(N, Q, Js, σ0s; rψ=rψ, α=α, typ=typ)
end