text:
using Revise, LinearAlgebra, StaticArrays, PyCall, FiniteDifferences, OptiMode, Plots, ChainRules, Zygote, Plots, HDF5, Dates  #, GLMakie, AbstractPlotting,
 #FFTW, BenchmarkTools, LinearMaps, IterativeSolvers, Roots, GeometryPrimitives
# pyplot()
# pygui()
SHM3 = SHermitianCompact{3,Float64,6}
mp = pyimport("meep")
mpb = pyimport("meep.mpb")

"""
################################################################################
#																			   #
#		         Define Utility Functions for Comparison of Julia 			   #
#		                        Solver with MPB 	                 		   #
#																			   #
################################################################################
"""

function write_sweep(sw_name;
                    data_dir="/home/dodd/data/OptiMode/grad_ng_p_rwg_dense/",
                    dt_fmt=dateformat"Y-m-d--H-M-S",
                    extension=".h5",
                    kwargs...)
    timestamp = Dates.format(now(),dt_fmt)
    fname = sw_name * "_" *  timestamp * extension
    @show fpath = data_dir * fname
    h5open(fpath, "cw") do file
        for (data_name,data) in kwargs
            write(file, string(data_name), data)
        end
    end
    return fpath
end

function read_sweep(sw_name;
                    data_dir="/home/dodd/data/OptiMode/grad_ng_p_rwg_dense/",
                    dt_fmt=dateformat"Y-m-d--H-M-S",
                    extension=".h5",
                    sw_keys=["ws","ts","p0","p̄_AD","p̄_FD","p̄_SJ"]
                    )
    # choose most recently timestamped file matching sw_name tag and extension
    fname = sort(  filter(x->(prod(split(x,"_")[begin:end-1])==prod(split(sw_name,"_"))),
                        readdir(data_dir));
                by=file->DateTime(split(file[begin:end-length(extension)],"_")[end],dt_fmt)
            )[end]
    @show fpath = data_dir * fname
    ds_data = h5open(fpath, "r") do file
        @show ds_keys = keys(file)
        ds_data = Dict([k=>read(file,k) for k in sw_keys]...)
    end
    return ds_data
end

# Generic data collection/loading fn use examples:
# fpath_test = write_sweep("wt";ws,ts,p0,p̄_AD,p̄_FD,p̄_SJ)
# ds_test = read_sweep("wt")

function get_Δs_mpb(ms)
    Δx, Δy, Δz = [(Δ == 0.0 ? 1.0 : Δ) for Δ in ms.geometry_lattice.size.__array__()]
end

function get_Ns_mpb(ms)
    Nx, Ny, Nz = [max(NN, 1) for NN in Int.(ms.resolution .* ms.geometry_lattice.size)]
end

function get_xyz_mpb(ms)
    Δx, Δy, Δz = get_Δs_mpb(ms) # [ (Δ==0. ? 1. : Δ) for Δ in ms_size.__array__() ]
    Nx, Ny, Nz = get_Ns_mpb(ms)
    x = ((Δx / Nx) .* (0:(Nx-1))) .- Δx / 2.0
    y = ((Δy / Ny) .* (0:(Ny-1))) .- Δy / 2.0
    z = ((Δz / Nz) .* (0:(Nz-1))) .- Δz / 2.0
    return x, y, z
end

function get_ε⁻¹_mpb(ms)
    x, y, z = get_xyz_mpb(ms)
    Nx = length(x)
    Ny = length(y)
    Nz = length(z)
    ε⁻¹ = Array{Float64,5}(undef, (3, 3, Nx, Ny, Nz))
    for i = 1:Nx, j = 1:Ny, k = 1:Nz
        ε⁻¹[:, :, i, j, k] .= real(ms.get_epsilon_inverse_tensor_point(mp.Vector3(
            x[i],
            y[j],
            z[k],
        )).__array__())
    end
    return ε⁻¹
end

function get_ε_mpb(ms)
    x, y, z = get_xyz_mpb(ms)
    Nx = length(x)
    Ny = length(y)
    Nz = length(z)
    ε = Array{Float64,5}(undef, (3, 3, Nx, Ny, Nz))
    for i = 1:Nx, j = 1:Ny, k = 1:Nz
        ε[:, :, i, j, k] .= real(ms.get_epsilon_inverse_tensor_point(mp.Vector3(
            x[i],
            y[j],
            z[k],
        )).inverse().__array__())
    end
    return ε
end

function ε_slices(
    ms;
    yind = 64,
    zind = 1,
    size = (800, 800),
    xlims = nothing,
    ylims = nothing,
    marker = :dot)
    x, y, z = get_xyz_mpb(ms)
    ε = get_ε_mpb(ms)
    p_diag = plot(
        x,
        ε[1, 1, :, yind, zind],
        label = "11",
        xlabel = "x [μm]",
        ylabel = "εᵢᵢ (diagonal elements)";
        xlims,
        ylims,
        marker,
    )
    plot!(p_diag, x, ε[2, 2, :, yind, zind], label = "22"; xlims, ylims, marker)
    plot!(p_diag, x, ε[3, 3, :, yind, zind], label = "33"; xlims, ylims, marker)
    p_offdiag = plot(
        x,
        ε[1, 2, :, yind, zind],
        label = "12",
        xlabel = "x [μm]",
        ylabel = "εᵢⱼ (off-diag. elements)";
        xlims,
        ylims,
        marker,
    )
    plot!(p_offdiag, x, ε[1, 3, :, yind, zind], label = "13"; xlims, ylims, marker)
    plot!(p_offdiag, x, ε[2, 3, :, yind, zind], label = "23"; xlims, ylims, marker)
    l = @layout [
        a
        b
    ]
    plot(p_diag, p_offdiag, layout = l, size = size)
end

"""
################################################################################
#																			   #
#		         Compare MPB and Julia solutions and gradients of 			   #
#		                Ridge Waveguide Dispersion Model 	           		   #
#																			   #
################################################################################
"""

"""
Default design parameters for ridge waveguide. Both MPB and OptiMode functions
should intake data in this format for convenient apples-to-apples comparison.
"""
p_def = [
    1.55,               #   wavelength              `λ`             [μm]
    1.7,                #   top ridge width         `w_top`         [μm]
    0.7,                #   ridge thickness         `t_core`        [μm]
    π / 14.0,           #   ridge sidewall angle    `θ`             [radian]
    2.4,                #   core index              `n_core`        [1]
    1.4,                #   substrate index         `n_subs`        [1]
]

pω_def = [
    1.45,               #   propagation constant    `kz`            [μm⁻¹]
    1.7,                #   top ridge width         `w_top`         [μm]
    0.7,                #   ridge thickness         `t_core`        [μm]
    π / 14.0,           #   ridge sidewall angle    `θ`             [radian]
    2.4,                #   core index              `n_core`        [1]
    1.4,                #   substrate index         `n_subs`        [1]
    0.5,                #   vacuum gap at boundaries `edge_gap`     [μm]
]

"""
MPB functions for ridge waveguide data
"""

function ms_rwg_mpb(
    p = p_def;
    Δx = 6.0,
    Δy = 4.0,
    Δz = 1.0,
    Nx = 128,
    Ny = 128,
    Nz = 1,
    edge_gap = 0.5)
    λ, w, t_core, θ, n_core, n_subs = p
    ω = 1 / λ
    nk = 10
    n_bands = 1
    res = mp.Vector3((Nx / Δx), (Ny / Δy), 1) #mp.Vector3(Int(Nx/Δx),Int(Ny/Δy),1) #16
    n_guess = 0.9 * n_core
    n_min = n_subs
    n_max = n_core
    t_subs = (Δy - t_core - edge_gap) / 2.0
    c_subs_y = -Δy / 2.0 + edge_gap / 2.0 + t_subs / 2.0
    # Set up MPB modesolver, use find-k to solve for one eigenmode `H` with prop. const. `k` at specified temporal freq. ω
    k_pts = mp.interpolate(nk, [mp.Vector3(0.05, 0, 0), mp.Vector3(0.05 * nk, 0, 0)])
    lat = mp.Lattice(size = mp.Vector3(Δx, Δy, 0))
    tanθ = tan(θ)
    tcore_tanθ = t_core * tanθ
    w_bottom = w + 2 * tcore_tanθ
    verts = [
        mp.Vector3(w / 2.0, t_core / 2.0, -5.0),
        mp.Vector3(-w / 2.0, t_core / 2.0, -5.0),
        mp.Vector3(-w_bottom / 2.0, -t_core / 2.0, -5.0),
        mp.Vector3(w_bottom / 2.0, -t_core / 2.0, -5.0),
    ]
    # verts = [ mp.Vector3(-w/2., -t_core/2., -5.),mp.Vector3(w, 2*t_core, -5.), mp.Vector3(w, -t_core/2., -5.)  ]
    core = mp.Prism(
        verts,
        10.0,
        axis = mp.Vector3(0.0, 0.0, 1.0),
        material = mp.Medium(index = n_core),
    )
    subs = mp.Block(
        size = mp.Vector3(Δx - edge_gap, t_subs, 10.0),
        center = mp.Vector3(0, c_subs_y, 0),
        material = mp.Medium(index = n_subs),
    )
    ms = mpb.ModeSolver(
        geometry_lattice = lat,
        geometry = [core, subs],
        k_points = k_pts,
        resolution = res,
        num_bands = n_bands,
        default_material = mp.vacuum,
        deterministic = true,
    )
    ms.init_params(mp.NO_PARITY, false)
    return ms
end

function nng_rwg_mpb(
    p = p_def;
    Δx = 6.0,
    Δy = 4.0,
    Δz = 1.0,
    Nx = 128,
    Ny = 128,
    Nz = 1,
    edge_gap = 0.5,
    band_idx = 1,
    tol = 1e-8)
    λ, w, t_core, θ, n_core, n_subs = p
    ω = 1 / λ
    nk = 10
    n_bands = 1
    res = mp.Vector3((Nx / Δx), (Ny / Δy), 1) #mp.Vector3(Int(Nx/Δx),Int(Ny/Δy),1) #16
    n_guess = 0.9 * n_core
    n_min = n_subs
    n_max = n_core
    ms = ms_rwg_mpb(p; Δx, Δy, Δz, Nx, Ny, Nz, edge_gap)
    kz = ms.find_k(
        mp.NO_PARITY,             # parity (meep parity object)
        ω,                    # ω at which to solve for k
        band_idx,                        # band_min (find k(ω) for bands
        band_idx,                        # band_max  band_min:band_max)
        mp.Vector3(0, 0, 1),      # k direction to search
        tol,                     # fractional k error tolerance
        n_guess * ω,              # kmag_guess, |k| estimate
        n_min * ω,                # kmag_min (find k in range
        n_max * ω,               # kmag_max  kmag_min:kmag_max)
    )[1]
    neff = kz / ω
    ng = 1 / ms.compute_one_group_velocity_component(mp.Vector3(0, 0, 1), band_idx)
    return [neff, ng]
end

function ε_rwg_mpb(
    p = p_def;
    Δx = 6.0,
    Δy = 4.0,
    Δz = 1.0,
    Nx = 128,
    Ny = 128,
    Nz = 1,
    edge_gap = 0.5,
    band_idx = 1,
    tol = 1e-8)
    λ, w, t_core, θ, n_core, n_subs = p
    n_guess = 0.9 * n_core
    ms = ms_rwg_mpb(p; Δx, Δy, Δz, Nx, Ny, Nz, edge_gap)
    ms.solve_kpoint(mp.Vector3(0, 0, 1) * n_guess / λ)
    return get_ε_mpb(ms)
end

function ε⁻¹_rwg_mpb(
    p = p_def;
    Δx = 6.0,
    Δy = 4.0,
    Δz = 1.0,
    Nx = 128,
    Ny = 128,
    Nz = 1,
    edge_gap = 0.5,
    band_idx = 1,
    tol = 1e-8)
    λ, w, t_core, θ, n_core, n_subs = p
    n_guess = 0.9 * n_core
    ms = ms_rwg_mpb(p; Δx, Δy, Δz, Nx, Ny, Nz, edge_gap)
    ms.solve_kpoint(mp.Vector3(0, 0, 1) * n_guess / λ)
    return get_ε⁻¹_mpb(ms)
end

∇nng_rwg_mpb_FD(
    p;
    Δx = 6.0,
    Δy = 4.0,
    Δz = 1.0,
    Nx = 128,
    Ny = 128,
    Nz = 1,
    edge_gap = 0.5,
    band_idx = 1,
    tol = 1e-8,
    nFD = 2) = FiniteDifferences.jacobian(
    central_fdm(nFD, 1),
    x -> nng_rwg_mpb(x; Δx, Δy, Δz, Nx, Ny, Nz, edge_gap, band_idx, tol),
    p,)[1]

wt_rwg_mpb(wt; NN = 128, bi = 1, tol = 1e-8) = nng_rwg_mpb(
    [p_def[1], wt[1], wt[2], p_def[4:6]...];
    Nx = NN,
    Ny = NN,
    band_idx = bi,
    tol,
)
∇wt_rwg_mpb_FD(wt; NN = 128, bi = 1, tol = 1e-8, nFD = 2) =
    FiniteDifferences.jacobian(central_fdm(nFD, 1), x -> wt_rwg_mpb(x; NN, bi, tol), wt)[1]

function wt_rwg_mpb_sweep(ws, ts, Ns; bi = 1, tol = 1e-8)
    nw = length(ws)
    nt = length(ts)
    nN = length(Ns)
    nng_mpb = zeros(Float64, (2, nw, nt, nN))
    ∇nng_mpb = zeros(Float64, (2, 2, nw, nt, nN))
    for wind = 1:nw, tind = 1:nt, Nind = 1:nN
        nng_mpb[:, wind, tind, Nind] =
            wt_rwg_mpb([ws[wind], ts[tind]]; NN = Ns[Nind], bi, tol)
        @views ∇nng_mpb[:, :, wind, tind, Nind] .=
            ∇wt_rwg_mpb_FD([ws[wind], ts[tind]]; NN = Ns[Nind], bi, tol)
    end
    return nng_mpb, ∇nng_mpb
end

function plot_ε_rwg_xslices_N(
    Ns,
    p = p_def;
    Δx = 6.0,
    Δy = 4.0,
    Δz = 1.0,
    Nz = 1,
    edge_gap = 0.5,
    size = (800, 800),
    xlims = nothing,
    ylims = nothing,
    marker = :dot,
    cs = (:red, :blue, :green, :purple, :orange, :magenta))

    nN = length(Ns)
    ms = ms_rwg_mpb(p; Δx, Δy, Δz, Nx = Ns[1], Ny = Ns[1], Nz, edge_gap)
    x, y, z = get_xyz_mpb(ms)
    iy = Int(ceil(Ns[1] / 2))
    εs1 = Array{Float64,3}(undef, (3, 3, Ns[1]))
    for ix = 1:Ns[1]
        εs1[:, :, ix] .= real(ms.get_epsilon_inverse_tensor_point(mp.Vector3(
            x[ix],
            y[iy],
            z[1],
        )).inverse().__array__())
    end
    p_diag = plot(
        x,
        εs1[1, 1, :],
        c = cs[1],
        label = "11",
        xlabel = "x [μm]",
        ylabel = "εᵢᵢ (diagonal elements)",
        legend = :bottomright;
        xlims,
        ylims,
        marker,
    )
    plot!(p_diag, x, εs1[2, 2, :], c = cs[2], label = "22"; xlims, ylims, marker)
    plot!(p_diag, x, εs1[3, 3, :], c = cs[3], label = "33"; xlims, ylims, marker)
    p_offdiag = plot(
        x,
        εs1[1, 2, :],
        c = cs[4],
        label = "12",
        xlabel = "x [μm]",
        ylabel = "εᵢⱼ (off-diag. elements)",
        legend = :bottomright;
        xlims,
        ylims,
        marker,
    )
    plot!(p_offdiag, x, εs1[1, 3, :], c = cs[5], label = "13"; xlims, ylims, marker)
    plot!(p_offdiag, x, εs1[2, 3, :], c = cs[6], label = "23"; xlims, ylims, marker)
    for Nind = 2:nN
        NN = Ns[Nind]
        ms = ms_rwg_mpb(p; Δx, Δy, Δz, Nx = NN, Ny = NN, Nz, edge_gap)
        x, y, z = get_xyz_mpb(ms)
        iy = Int(ceil(NN[1] / 2))
        εs = Array{Float64,3}(undef, (3, 3, NN))
        for ix = 1:NN
            εs[:, :, ix] .= real(ms.get_epsilon_inverse_tensor_point(mp.Vector3(
                x[ix],
                y[iy],
                z[1],
            )).inverse().__array__())
        end
        plot!(p_diag, x, εs[1, 1, :], c = cs[1], label = nothing; xlims, ylims, marker)
        plot!(p_diag, x, εs[2, 2, :], c = cs[2], label = nothing; xlims, ylims, marker)
        plot!(p_diag, x, εs[3, 3, :], c = cs[3], label = nothing; xlims, ylims, marker)
        plot!(p_offdiag, x, εs[1, 2, :], c = cs[4], label = nothing; xlims, ylims, marker)
        plot!(p_offdiag, x, εs[1, 3, :], c = cs[5], label = nothing; xlims, ylims, marker)
        plot!(p_offdiag, x, εs[2, 3, :], c = cs[6], label = nothing; xlims, ylims, marker)
    end
    l = @layout [
        a
        b
    ]
    plot(p_diag, p_offdiag, layout = l, size = size)
end

function ε_rwg_xslices_N(
    Ns,
    p = p_def;
    Δx = 6.0,
    Δy = 4.0,
    Δz = 1.0,
    Nz = 1,
    edge_gap = 0.5)
    nN = length(Ns)
    ms = ms_rwg_mpb(p; Δx, Δy, Δz, Nx = Ns[end], Ny = Ns[end], Nz, edge_gap)
    x, y, z = get_xyz_mpb(ms)
    # iy = Int(ceil(Ns[end] / 2))
    εs = Array{Float64,4}(undef, (3, 3, Ns[end], nN))
    for ix = 1:Ns[end]
        εs[:, :, ix, nN] .= real(ms.get_epsilon_inverse_tensor_point(mp.Vector3(
            x[ix],
            0.0,
            0.0,
        )).inverse().__array__())
    end
    for Nind = 1:nN-1
        NN = Ns[Nind]
        ms = ms_rwg_mpb(p; Δx, Δy, Δz, Nx = NN, Ny = NN, Nz, edge_gap)
        # xtemp, y, z = get_xyz_mpb(ms)
        # iy = Int(ceil(NN[1] / 2))
        for ix = 1:Ns[end]
            εs[:, :, ix, Nind] .= real(ms.get_epsilon_inverse_tensor_point(mp.Vector3(
                x[ix],
                0.0,
                0.0,
            )).inverse().__array__())
        end
    end
    return x, εs
end

function ε_rwg_xslices_w(
    ws;
    p = [
        1.55,               #   wavelength              `λ`             [μm]
        1.7,                #   top ridge width         `w_top`         [μm]
        0.7,                #   ridge thickness         `t_core`        [μm]
        π / 10.0,                 #   ridge sidewall angle    `θ`             [radian]
        2.4,                #   core index              `n_core`        [1]
        1.4,                #   substrate index         `n_subs`        [1]
    ],
    Δx = 6.0,
    Δy = 4.0,
    Δz = 1.0,
    Nx = 256,
    Ny = 256,
    Nz = 1,
    edge_gap = 0.5)
    nw = length(ws)
    ms = ms_rwg_mpb([
        1.55,               #   wavelength              `λ`             [μm]
        ws[1],                #   top ridge width         `w_top`         [μm]
        0.7,                #   ridge thickness         `t_core`        [μm]
        π / 10.0,                 #   ridge sidewall angle    `θ`             [radian]
        2.4,                #   core index              `n_core`        [1]
        1.4,                #   substrate index         `n_subs`        [1]
    ]; Δx, Δy, Δz, Nx, Ny, Nz, edge_gap)
    x, y, z = get_xyz_mpb(ms)
    εs = Array{Float64,4}(undef, (3, 3, length(x), nw))
    for ix = 1:Nx
        εs[:, :, ix, 1] .= real(ms.get_epsilon_inverse_tensor_point(mp.Vector3(
            x[ix],
            0.0,
            0.0,
        )).inverse().__array__())
    end
    for wind = 2:length(ws)
        ww = ws[wind]
        ms = ms_rwg_mpb(
            [
                1.55,               #   wavelength              `λ`             [μm]
                ws[wind],                #   top ridge width         `w_top`         [μm]
                0.7,                #   ridge thickness         `t_core`        [μm]
                π / 10.0,                 #   ridge sidewall angle    `θ`             [radian]
                2.4,                #   core index              `n_core`        [1]
                1.4,                #   substrate index         `n_subs`        [1]
            ];
            Δx,
            Δy,
            Δz,
            Nx,
            Ny,
            Nz,
            edge_gap,
        )
        # xtemp, y, z = get_xyz_mpb(ms)
        # iy = Int(ceil(NN[1] / 2))
        for ix = 1:Nx
            εs[:, :, ix, wind] .= real(ms.get_epsilon_inverse_tensor_point(mp.Vector3(
                x[ix],
                0.0,
                0.0,
            )).inverse().__array__())
        end
    end
    return x, εs
end

"""
OptiMode functions for ridge waveguide data
"""
function nngω_rwg_OM(p::Vector{Float64} = pω_def;
                    Δx = 6.0,
                    Δy = 4.0,
                    Δz = 1.0,
                    Nx = 128, #16,
                    Ny = 128, #16,
                    Nz = 1,
                    band_idx = 1,
                    tol = 1e-8)
                    # kz, w, t_core, θ, n_core, n_subs, edge_gap = p
                    # nng_tuple = solve_nω(kz,ridge_wg(w,t_core,θ,edge_gap,n_core,n_subs,Δx,Δy),Δx,Δy,Δz,Nx,Ny,Nz;tol)
                    nng_tuple = solve_nω(p[1],ridge_wg(p[2],p[3],p[4],p[7],p[5],p[6],Δx,Δy),Δx,Δy,Δz,Nx,Ny,Nz;tol)
                    [nng_tuple[1],nng_tuple[2]]
end

function nng_rwg_OM(p = p_def;
                    Δx = 6.0,
                    Δy = 4.0,
                    Δz = 1.0,
                    Nx = 128,
                    Ny = 128,
                    Nz = 1,
                    edge_gap = 0.5,
                    band_idx = 1,
                    tol = 1e-8)
                    λ, w, t_core, θ, n_core, n_subs = p
                    ω = 1/λ
                    # solve_n(ω,ridge_wg(w,t_core,edge_gap,n_core,n_subs,dropgrad(Δx),dropgrad(Δy)),dropgrad(Δx),dropgrad(Δy),dropgrad(Δz),dropgrad(Nx),dropgrad(Ny),dropgrad(Nz))
                    nng_tuple = solve_n(ω,ridge_wg(w,t_core,θ,edge_gap,n_core,n_subs,Δx,Δy),Δx,Δy,Δz,Nx,Ny,Nz;tol)
                    [nng_tuple[1],nng_tuple[2]]
end

function ε_rwg_OM(p = p_def;
                    Δx = 6.0,
                    Δy = 4.0,
                    Δz = 1.0,
                    Nx = 128,
                    Ny = 128,
                    Nz = 1,
                    edge_gap = 0.5,
                    npix_sm = 2 )
                    λ, w, t_core, θ, n_core, n_subs = p
                    # εₛ(ridge_wg(w,t_core,edge_gap,n_core,n_subs,dropgrad(Δx),dropgrad(Δy)),dropgrad(Δx),dropgrad(Δy),dropgrad(Nx),dropgrad(Ny);npix_sm)
                    reshape(
                        permutedims(
                            εₛ(ridge_wg(w,t_core,θ,edge_gap,n_core,n_subs,Δx,Δy),Δx,Δy,Nx,Ny;npix_sm),
                            (3,4,1,2),
                        ),
                        (3,3,Nx,Ny,1),
                    )
end

function ε⁻¹_rwg_OM(p = p_def;
                    Δx = 6.0,
                    Δy = 4.0,
                    Δz = 1.0,
                    Nx = 128,
                    Ny = 128,
                    Nz = 1,
                    edge_gap = 0.5,
                    npix_sm = 2 )
                    λ, w, t_core, θ, n_core, n_subs = p
                    # εₛ⁻¹(ridge_wg(w,t_core,edge_gap,n_core,n_subs,dropgrad(Δx),dropgrad(Δy)),dropgrad(Δx),dropgrad(Δy),dropgrad(Nx),dropgrad(Ny);npix_sm)
                    reshape(
                        permutedims(
                            εₛ⁻¹(ridge_wg(w,t_core,θ,edge_gap,n_core,n_subs,Δx,Δy),Δx,Δy,Nx,Ny;npix_sm),
                            (3,4,1,2),
                        ),
                        (3,3,Nx,Ny,1),
                    )
end

∇nngω_rwg_OM_FD(
    p;
    Δx = 6.0,
    Δy = 4.0,
    Δz = 1.0,
    Nx = 128,
    Ny = 128,
    Nz = 1,
    edge_gap = 0.5,
    band_idx = 1,
    tol = 1e-8,
    nFD = 3) = FiniteDifferences.jacobian(
    central_fdm(nFD, 1),
    x -> nngω_rwg_OM(x; Δx, Δy, Δz, Nx, Ny, Nz, edge_gap, band_idx, tol),
    p,)[1]

∇nng_rwg_OM_FD(
    p;
    Δx = 6.0,
    Δy = 4.0,
    Δz = 1.0,
    Nx = 128,
    Ny = 128,
    Nz = 1,
    edge_gap = 0.5,
    band_idx = 1,
    tol = 1e-8,
    nFD = 2) = FiniteDifferences.jacobian(
    central_fdm(nFD, 1),
    x -> nng_rwg_OM(x; Δx, Δy, Δz, Nx, Ny, Nz, edge_gap, band_idx, tol),
    p,)[1]

wtω_rwg_OM(wt; NN = 128, bi = 1, tol = 1e-8) = nngω_rwg_OM(
    [pω_def[1], wt[1], wt[2], pω_def[4:7]...];
    Nx = NN,
    Ny = NN,
    band_idx = bi,
    tol,
)
∇wtω_rwg_OM_FD(wt; NN = 128, bi = 1, tol = 1e-8, nFD = 3) =
    FiniteDifferences.jacobian(central_fdm(nFD, 1), x -> wtω_rwg_OM(x; NN, bi, tol), wt)[1]

function wtω_rwg_OM_sweep(sw_name, ws, ts, Ns;
    bi = 1,
    tol = 1e-8,
    data_dir="/home/dodd/data/OptiMode/mpb_compare_rwg/",
    dt_fmt=dateformat"Y-m-d--H-M-S",
    extension=".h5",
    )
    timestamp = Dates.format(now(),dt_fmt)
    fname = sw_name * "_" *  timestamp * extension
    @show fpath = data_dir * fname
    nw = length(ws)
    nt = length(ts)
    nN = length(Ns)
    nng_OM = zeros(Float64, (2, nw, nt, nN))
    ∇nng_OM_FD = zeros(Float64, (2, 2, nw, nt, nN))
    ∇nng_OM_AD = zeros(Float64, (2, 2, nw, nt, nN))
    for wind = 1:nw, tind = 1:nt, Nind = 1:nN
        println("wind: $wind of $nw,  tind: $tind of $nt,  Nind: $Nind of $nN")
        wt = [ws[wind], ts[tind]]
        nng,nng_pb = Zygote.pullback(x -> wtω_rwg_OM(x; NN = Ns[Nind], bi, tol), wt)
        nng_OM[:, wind, tind, Nind] = nng
        @views ∇nng_OM_AD[:,:,wind,tind,Nind] .= [ real(nng_pb([1.,0.])[1])'    # = [   ∂n/∂w   ∂n/∂t
                                                   real(nng_pb([0.,1.])[1])' ]  #       ∂ng/∂w  ∂ng/∂t  ]
        @views ∇nng_OM_FD[:, :, wind, tind, Nind] .=
            ∇wtω_rwg_OM_FD([ws[wind], ts[tind]]; NN = Ns[Nind], bi, tol)
        h5open(fpath, "w") do file
            @write file nng_OM
            @write file ws
            @write file ts
            @write file Ns
            @write file tol
            write(file, "band_index", bi)
            write(file, "grad_nng_OM_FD", ∇nng_OM_FD)
            write(file, "grad_nng_OM_AD", ∇nng_OM_AD)
        end
    end
    return nng_OM, ∇nng_OM_AD, ∇nng_OM_FD
end

wt_rwg_OM(wt; NN = 128, bi = 1, tol = 1e-8) = nng_rwg_OM(
    [p_def[1], wt[1], wt[2], p_def[4:6]...];
    Nx = NN,
    Ny = NN,
    band_idx = bi,
    tol,
)
∇wt_rwg_OM_FD(wt; NN = 128, bi = 1, tol = 1e-8, nFD = 3) =
    FiniteDifferences.jacobian(central_fdm(nFD, 1), x -> wt_rwg_OM(x; NN, bi, tol), wt)[1]

function wt_rwg_OM_sweep(ws, ts, Ns; bi = 1, tol = 1e-8)
    nw = length(ws)
    nt = length(ts)
    nN = length(Ns)
    nng_OM = zeros(Float64, (2, nw, nt, nN))
    ∇nng_OM = zeros(Float64, (2, 2, nw, nt, nN))
    for wind = 1:nw, tind = 1:nt, Nind = 1:nN
        nng_OM[:, wind, tind, Nind] =
            wt_rwg_OM([ws[wind], ts[tind]]; NN = Ns[Nind], bi, tol)
        @views ∇nng_OM[:, :, wind, tind, Nind] .=
            ∇wt_rwg_OM_FD([ws[wind], ts[tind]]; NN = Ns[Nind], bi, tol)
    end
    return nng_OM, ∇nng_OM
end

##
# @assert size(ε_rwg_mpb(p_def)) == size(ε_rwg_OM(p_def))
# @assert size(ε⁻¹_rwg_mpb(p_def)) == size(ε⁻¹_rwg_OM(p_def))

# ∇wtω_rwg_OM_FD([1.5, 0.7]; NN = 64, bi = 1, tol = 1e-8, nFD = 3)
# nng_OM, ∇nng_OM_AD, ∇nng_OM_FD = wtω_rwg_OM_sweep(1.3:0.2:1.7, .5:0.2:.7, 2 .^(5:7); bi = 1, tol = 1e-8)
# ∇nng_err_OM = abs.(∇nng_OM_FD .- ∇nng_OM_AD) ./ abs.(∇nng_OM_FD)
##
function compare_ε11_rwg(p=p_def;
    Δx = 6.0,
    Δy = 4.0,
    Δz = 1.0,
    Nx = 128,
    Ny = 128,
    Nz = 1,
    edge_gap = 0.5,
    )
    eps_mpb = ε_rwg_mpb(p;Δx,Δy,Δz,Nx,Ny,Nz,edge_gap)
    eps_OM = ε_rwg_OM(p;Δx,Δy,Δz,Nx,Ny,Nz,edge_gap)
    x = ((Δx / Nx) .* (0:(Nx-1))) .- Δx / 2.0
    y = ((Δy / Ny) .* (0:(Ny-1))) .- Δy / 2.0
    hm_eps11_mpb = heatmap(x,
        y,
        eps_mpb[1,1,:,:,1]',
        xlabel = "x [μm]",
        ylabel = "y [μm]",
        title = "Smoothed ε₁₁, MPB",
        aspect_ratio=:equal,
        )
    hm_eps11_OM = heatmap(x,
        y,
        eps_OM[1,1,:,:,1]',
        xlabel = "x [μm]",
        ylabel = "y [μm]",
        title = "Smoothed ε₁₁, OptiMode",
        aspect_ratio=:equal,
        )
    l = @layout[  a  b  ]
    plot(hm_eps11_mpb,hm_eps11_OM,layout=l)
end
# compare_ε11_rwg(p_def)
##
# Ns = Int.(ceil.(2 .^ collect(6:0.5:9)))
# Ns = Int.(ceil.(2 .^ collect(4.4:0.02:8)))
# x_slices_N, ε_xslices_N = ε_rwg_xslices_N(Ns, [
#     1.55,               #   wavelength              `λ`             [μm]
#     1.7,                #   top ridge width         `w_top`         [μm]
#     0.7,                #   ridge thickness         `t_core`        [μm]
#     π / 10.0,              #   ridge sidewall angle    `θ`             [radian]
#     2.4,                #   core index              `n_core`        [1]
#     1.4,                #   substrate index         `n_subs`        [1]
# ])
# # Makie.surface(x_slices[75:105],(.01 .* Ns),ε_xslices[1,1,75:105,:])
# ixmin, ixmax = 80, 102
# sN11 = Plots.surface(Ns, x_slices_N[ixmin:ixmax], ε_xslices_N[1, 1, ixmin:ixmax, :])
# sN22 = Plots.surface(Ns, x_slices_N[ixmin:ixmax], ε_xslices_N[2, 2, ixmin:ixmax, :])
# sN12 = Plots.surface(Ns, x_slices_N[ixmin:ixmax], ε_xslices_N[1, 2, ixmin:ixmax, :])
# sN13 = Plots.surface(Ns, x_slices_N[ixmin:ixmax], ε_xslices_N[1, 3, ixmin:ixmax, :])
# sN23 = Plots.surface(Ns, x_slices_N[ixmin:ixmax], ε_xslices_N[2, 3, ixmin:ixmax, :])
#
# ws = collect(1.0:0.005:1.2)
# x_slices_w, ε_xslices_w = ε_rwg_xslices_w(ws)
# #ixmin, ixmax = 97, 113 #50,56
# ixmin, ixmax = 1, 256
# # sw11 = Plots.surface(Ns,x_slices_w,ε_xslices_w[1,1,:,:])
# sw11 = Plots.surface(ws, x_slices_w[ixmin:ixmax], ε_xslices_w[1, 1, ixmin:ixmax, :])
# sw22 = Plots.surface(ws, x_slices_w[ixmin:ixmax], ε_xslices_w[2, 2, ixmin:ixmax, :])
# sw33 = Plots.surface(ws, x_slices_w[ixmin:ixmax], ε_xslices_w[3, 3, ixmin:ixmax, :])
# sw12 = Plots.surface(ws, x_slices_w[ixmin:ixmax], ε_xslices_w[1, 2, ixmin:ixmax, :])
# sw13 = Plots.surface(ws, x_slices_w[ixmin:ixmax], ε_xslices_w[1, 3, ixmin:ixmax, :])
# sw23 = Plots.surface(ws, x_slices_w[ixmin:ixmax], ε_xslices_w[2, 3, ixmin:ixmax, :])

#
#
# using Plots
# pyplot()
# pygui()
# p_xslices = ε_rwg_xslices_N(Ns, p_def, xlims = (-1.1, -0.8))
#
# p_slices64 = ε_slices(ms_rwg_mpb(p_def; Nx = 64, Ny = 64), xlims = (-1.1, -0.8), yind = 32)
# p_slices128 =
#     ε_slices(ms_rwg_mpb(p_def; Nx = 128, Ny = 128), xlims = (-1.1, -0.8), yind = 64)
# p_slices256 =
#     ε_slices(ms_rwg_mpb(p_def; Nx = 256, Ny = 256), xlims = (-1.1, -0.8), yind = 128)
#
# p_slices256.subplots[1]
#
# nng_rwg_mpb(p_def)
# ∇nng_rwg_mpb_FD(p_def)
# wt_rwg_mpb([1.0, 0.7]; NN = 64, bi = 1, tol = 1e-5)
# grad_nng = ∇wt_rwg_mpb_FD([1.0, 0.7]; NN = 64, bi = 1, tol = 1e-5)

# function wt_rwg_sweep(ws,ts,Ns;bi=1,tol=1e-5,name="rwg_ng_bm_N",path="/home/dodd/github/OptiMode/test/")
#     bm = [ @benchmark(outs[$i] = ngom_rwg($ks,$w,$t_core,$edge_gap,$n_core,$n_subs,$Δx,$Δy,$Δz,$Ns[$i],$Ns[$i],1))  for i=1:length(Ns) ]
#     ns = [o[1] for o in outs]
#     ngs = [o[2] for o in outs]
#     # bm = [ @benchmark(ngom_rwg(1.6,$w,$t_core,$edge_gap,$n_core,$n_subs,$Δx,$Δy,$Δz,$nn,$nn,1))  for nn=Ns ]
#     df = DataFrame( N = Ns,
#                     Nsq = (Ns.^2),
#                     n = ns,
#                     ng = ngs,
#                     time = (time.(bm).*1e-9),
#                     gctime_frac = (gctime.(bm) ./ time.(bm)),
#                     mem = memory.(bm),
#                     allocs = allocs.(bm),
#                    )
#     CSV.write(*(path,name,".csv"), df)
#     # plot_bm_single(df)
#     df
# end

########################################
##### w-t-N sweep 1 (MPB only) #########
########################################
# ws = collect(0.8:0.1:1.7)
# ts = collect(0.5:0.1:1.3)
# Ns = Int.(2 .^(collect(5:8)))
# nng_mpb, ∇nng_mpb = wt_rwg_mpb_sweep(ws,ts,Ns;bi=1)
# ws1 = collect(0.8:0.1:1.7)
# ts1 = collect(0.5:0.1:1.3)
# Ns1 = Int.(2 .^(collect(5:8)))
# nng_mpb1 = copy(nng_mpb)
# ∇nng_mpb1 = copy(∇nng_mpb)
# n_mpb1 = nng_mpb1[1,:,:,:]
# ng_mpb1 = nng_mpb1[2,:,:,:]
# ∂n∂w_mpb1 = ∇nng_mpb1[1,1,:,:,:]
# ∂n∂t_mpb1 = ∇nng_mpb1[1,2,:,:,:]
# ∂ng∂w_mpb1 = ∇nng_mpb1[2,1,:,:,:]
# ∂ng∂t_mpb1 = ∇nng_mpb1[2,2,:,:,:]
# mpb_fname1="rwg_mpb_wt1.h5"
# # path="/home/dodd/github/OptiMode/test/"
# h5open(mpb_fname1, "cw") do file
#     write(file, "ws", ws1)  # alternatively, say "@write file A"
#     write(file, "ts", ts1)
#     write(file, "Ns", Ns1)
#     write(file, "nng", nng_mpb1)
#     write(file, "∇nng", ∇nng_mpb1)
#     write(file, "∂n∂w", ∂n∂w_mpb1)
#     write(file, "∂n∂t", ∂n∂t_mpb1)
#     write(file, "∂ng∂w", ∂ng∂w_mpb1)
#     write(file, "∂ng∂t", ∂ng∂t_mpb1)
# end
# ∇nng_mpb1_read = h5open(mpb_fname1, "r") do file
#     read(file, "∇nng")
# end

function plot_wt_sweep(ws, ts, Ns, nng, ∇nng; tind_w = 4, size = (800, 800))
    n = nng[1, :, :, :]
    ng = nng[2, :, :, :]
    ∂n∂w = ∇nng[1, 1, :, :, :]
    ∂n∂t = ∇nng[1, 2, :, :, :]
    ∂ng∂w = ∇nng[2, 1, :, :, :]
    ∂ng∂t = ∇nng[2, 2, :, :, :]
    # w slices
    p_n_w = plot(
        ws,
        n[:, tind_w, 1],
        label = "N=$(Ns[1])",
        xlabel = "top width [μm]",
        ylabel = "effective index n",
        legend = :bottomright,
    )
    p_ng_w = plot(
        ws,
        ng[:, tind_w, 1],
        label = "N=$(Ns[1])",
        xlabel = "top width [μm]",
        ylabel = "eff. group index ng",
        legend = false,
    )
    p_∂n∂w_w = plot(
        ws,
        ∂n∂w[:, tind_w, 1],
        label = "N=$(Ns[1])",
        xlabel = "top width [μm]",
        ylabel = "∂n/∂w [μm⁻¹]",
        legend = false,
    )
    p_∂ng∂w_w = plot(
        ws,
        ∂ng∂w[:, tind_w, 1],
        label = "N=$(Ns[1])",
        xlabel = "top width [μm]",
        ylabel = "∂ng/∂w [μm⁻¹]",
        legend = false,
    )
    p_∂n∂t_w = plot(
        ws,
        ∂n∂t[:, tind_w, 1],
        label = "N=$(Ns[1])",
        xlabel = "top width [μm]",
        ylabel = "∂n/∂t [μm⁻¹]",
        legend = false,
    )
    p_∂ng∂t_w = plot(
        ws,
        ∂ng∂t[:, tind_w, 1],
        label = "N=$(Ns[1])",
        xlabel = "top width [μm]",
        ylabel = "∂ng/∂t [μm⁻¹]",
        legend = false,
    )
    for Nind = 2:length(Ns)
        plot!(
            p_n_w,
            ws,
            n[:, tind_w, Nind],
            label = "N=$(Ns[Nind])",
            xlabel = "top width [μm]",
            ylabel = "effective index n",
        )
        plot!(
            p_ng_w,
            ws,
            ng[:, tind_w, Nind],
            label = "N=$(Ns[Nind])",
            xlabel = "top width [μm]",
            ylabel = "eff. group index ng",
        )
        plot!(
            p_∂n∂w_w,
            ws,
            ∂n∂w[:, tind_w, Nind],
            label = "N=$(Ns[Nind])",
            xlabel = "top width [μm]",
            ylabel = "∂n/∂w [μm⁻¹]",
        )
        plot!(
            p_∂ng∂w_w,
            ws,
            ∂ng∂w[:, tind_w, Nind],
            label = "N=$(Ns[Nind])",
            xlabel = "top width [μm]",
            ylabel = "∂ng/∂w [μm⁻¹]",
        )
        plot!(
            p_∂n∂t_w,
            ws,
            ∂n∂t[:, tind_w, Nind],
            label = "N=$(Ns[Nind])",
            xlabel = "top width [μm]",
            ylabel = "∂n/∂t [μm⁻¹]",
        )
        plot!(
            p_∂ng∂t_w,
            ws,
            ∂ng∂t[:, tind_w, Nind],
            label = "N=$(Ns[Nind])",
            xlabel = "top width [μm]",
            ylabel = "∂ng/∂t [μm⁻¹]",
        )
    end
    l_w = @layout [
        a b
        c d
        e f
    ]
    plot(p_n_w, p_ng_w, p_∂n∂w_w, p_∂ng∂w_w, p_∂n∂t_w, p_∂ng∂t_w, layout = l_w, size = size)
end

# plot_wt_sweep(ws,ts,Ns,nng_mpb,∇nng_mpb;tind_w=2)


text:
printstyled("  RandSeis utils\n", color=:light_green)

fs_range = exp10.(range(-6, stop=4, length=50))
fc_range = exp10.(range(-4, stop=2, length=20))

for fs in fs_range
  for fc in fc_range
    getbandcode(fs, fc=fc)
  end
end

# This should do it
D = Dict{String,Any}()
for i = 1:10
  RandSeis.pop_rand_dict!(D, 1000)
end

# Similarly for the yp2 codes
for i = 1:100
  i,c,u = RandSeis.getyp2codes('s', true)
  @test isa(i, Char)
  @test isa(c, Char)
  @test isa(u, String)
end

R = RandSeis.randResp(8)
@test length(R.z) == length(R.p) == 8

for i = 1:10000
  i,c,u = RandSeis.getyp2codes('s', false)
  @test isa(i, Char)
  @test isa(c, Char)
  @test isa(u, String)
end

printstyled("  randSeis*\n", color=:light_green)
for i = 1:10
  randSeisChannel()
  randSeisData()
  randSeisHdr()
  randSeisEvent()
end

printstyled("  namestrip, namestrip!\n", color=:light_green)
str = String(0x00:0xff)
S = randSeisData(3)
S.name[2] = str

for key in keys(bad_chars)
  test_str = namestrip(str, key)
  @test length(test_str) == 256 - (32 + length(bad_chars[key]))
end
redirect_stdout(out) do
  test_str = namestrip(str, "Nonexistent List")
end
namestrip!(S)
@test length(S.name[2]) == 210

ext: jl


text:
module Test_connections

using DescriptorSystems
using LinearAlgebra
using Test

println("Test_connections")
@testset "test_connections" begin

@testset "Standard state space" begin
# CONTINUOUS
C_111 = dss([1], [2], [3], [4])
C_211 = dss(eye(2), [1; 2], [1 0], [0])
C_212 = dss(eye(2), [1; 2], eye(2), [0; 0])
C_221 = dss(eye(2), [1 0; 0 2], [1 0], [0 0])
C_222 = dss(eye(2), [1 0; 0 2], eye(2), zeros(Int,2,2))
C_022 = dss(4*eye(2))

# DISCRETE
D_111 = dss([1], [2], [3], [4], Ts = 0.005)
D_211 = dss(eye(2), [1; 2], [1 0], [0], Ts = 0.005)
D_212 = dss(eye(2), [1;2], eye(2), [0; 0], Ts = 0.005)
D_221 = dss(eye(2), [1 0; 0 2], [1 0], [0 0], Ts = 0.005)
D_222 = dss(eye(2), [1 0; 0 2], eye(2), zeros(Int,2,2), Ts = 0.005)
D_022 = dss(4*eye(2), Ts = 0.005)

@test [C_111 C_221] == dss(eye(3), [2 0 0; 0 1 0; 0 0 2], [3 1 0], [4 0 0])
@test [C_111; C_212] == dss(eye(3), [2; 1; 2], [3 0 0; 0 1 0; 0 0 1], [4; 0; 0])
@test append(C_111, C_211) == dss(eye(3), [2 0; 0 1; 0 2], [3 0 0; 0 1 0], [4 0; 0 0])
@test [C_022 C_222] == dss(eye(2), [0 0 1 0; 0 0 0 2], [1 0; 0 1], [4 0 0 0; 0 4 0 0])
@test horzcat(C_022, C_222) == dss(eye(2), [0 0 1 0; 0 0 0 2], [1 0; 0 1], [4 0 0 0; 0 4 0 0])
@test [C_022; C_222] == dss(eye(2), [1 0; 0 2], [0 0; 0 0; 1 0; 0 1], [4 0; 0 4; 0 0; 0 0])
@test vertcat(C_022, C_222) == dss(eye(2), [1 0; 0 2], [0 0; 0 0; 1 0; 0 1], [4 0; 0 4; 0 0; 0 0])

@test [D_111 D_221] == dss(eye(3), [2 0 0; 0 1 0; 0 0 2], [3 1 0], [4 0 0], Ts = 0.005)
@test [D_111; D_212] == dss(eye(3), [2; 1; 2], [3 0 0; 0 1 0; 0 0 1], [4; 0; 0], Ts = 0.005)
@test append(D_111, D_211) == dss(eye(3), [2 0; 0 1; 0 2], [3 0 0; 0 1 0], [4 0; 0 0], Ts = 0.005)
@test [D_022 D_222] == dss(eye(2), [0 0 1 0; 0 0 0 2], [1 0; 0 1], [4 0 0 0; 0 4 0 0], Ts = 0.005)
@test horzcat(D_022, D_222) == dss(eye(2), [0 0 1 0; 0 0 0 2], [1 0; 0 1], [4 0 0 0; 0 4 0 0], Ts = 0.005)
@test [D_022; D_222] == dss(eye(2), [1 0; 0 2], [0 0; 0 0; 1 0; 0 1], [4 0; 0 4; 0 0; 0 0], Ts = 0.005)
@test vertcat(D_022, D_222) == dss(eye(2), [1 0; 0 2], [0 0; 0 0; 1 0; 0 1], [4 0; 0 4; 0 0; 0 0], Ts = 0.005)

@test series(C_111, C_212) == C_212*C_111
@test parallel(C_111, C_211) == C_111 + C_211

sys = rss(3,4,5); K = rand(2,3); ind = 3:4; out = 2:4;
@test iszero(feedback(feedback(sys,K,ind,out),K,ind,out; negative = false)-sys,atol=1.e-7)

# Errors
@test_throws ErrorException [C_111 D_111]                 # Sampling time mismatch
@test_throws ErrorException [C_111; D_111]                # Sampling time mismatch
@test_throws ErrorException append(C_111, D_111)          # Sampling time mismatch
@test_throws ErrorException [C_111 C_212]                 # Dimension mismatch
@test_throws ErrorException [C_111; C_221]                # Dimension mismatch

# hcat and vcat for StateSpace and Matrix
A = [-1.1 -1.2; -1.3 -1.4]
B = [1 2; 3 4]
C = [5 6; 7 8]
D = [1 0; 0 1]
P = dss(A, B, C, D)
@test [P fill(2.5, 2, 1)] == dss(A, [B fill(0, 2, 1)], C, [D fill(2.5, 2, 1)])
@test [fill(2.5, 2, 1) P] == dss(A, [fill(0, 2, 1) B], C, [fill(2.5, 2, 1) D])
@test [P; fill(2.5, 1, 2)] == dss(A, B, [C; fill(0, 1, 2)], [D; fill(2.5, 1, 2)])
@test [fill(2.5, 1, 2); P] == dss(A, B, [fill(0, 1, 2); C], [fill(2.5, 1, 2); D])

# hcat and vcat for StateSpace and Number
P = dss(-1.0, 2.0, 3.0, 4.0)
@test [P 2.5] == dss(-1.0, [2.0 0.0], 3.0, [4.0 2.5])
@test [2.5 P] == dss(-1.0, [0.0 2.0], 3.0, [2.5 4.0])
@test [P; 2.5] == dss(-1.0, 2.0, [3.0; 0.0], [4.0; 2.5])
@test [2.5; P] == dss(-1.0, 2.0, [0.0; 3.0], [2.5; 4.0])

#
@test [2.5I P 3.5] == dss(-1.0, [0.0 2.0 0.0], 3.0, [2.5 4.0 3.5])
@test [2.5; P; 3.5I] == dss(-1.0, 2.0, [0.0; 3.0; 0.0], [2.5; 4.0; 3.5])

# Concatenation of discrete system with constant
@test [D_111 1.] == dss([1.0], [2.0 0.0], [3.0], [4.0 1.0], Ts = 0.005)
@test [1. D_111] == dss([1.0], [0.0 2.0], [3.0], [1.0 4.0], Ts = 0.005)
@test [D_111 I] == dss([1.0], [2.0 0.0], [3.0], [4.0 1.0], Ts = 0.005)
@test [2I D_111] == dss([1.0], [0.0 2.0], [3.0], [2.0 4.0], Ts = 0.005)
# Type and sample time
@test [D_111 1.] isa DescriptorStateSpace{Float64}
@test [D_111 I].Ts == 0.005
# Continuous version
@test [C_111 1.] == dss([1.0], [2.0 0.0], [3.0], [4.0 1.0])
@test [1. C_111] == dss([1.0], [0.0 2.0], [3.0], [1.0 4.0])
@test [C_111 1.] isa DescriptorStateSpace{Float64}
@test [C_111 5.0I].Ts == 0.0

# Concatenation of discrete system with matrix
@test [D_222 fill(1.5, 2, 2)] == [D_222 dss(fill(1.5, 2, 2),Ts = 0.005)]
@test [C_222 fill(1.5, 2, 2)] == [C_222 dss(fill(1.5, 2, 2))]

# hvcat numbers
@test [C_111 1.5; 2 3I] ==
    [C_111 dss(1.5); dss(2.0) dss(3.0)]
@test [D_111 1.5; 2 3] ==
    [D_111 dss(1.5,Ts = 0.005); dss(2.0,Ts = 0.005) dss(3.0,Ts = 0.005)]
# hvcat matrices
@test [C_222 fill(1.5, 2, 2); fill(2, 2, 2) fill(3, 2, 2)] ==
    [C_222 dss(fill(1.5, 2, 2)); dss(fill(2, 2, 2)) dss(fill(3, 2, 2))]
@test [D_222 fill(1.5, 2, 2); fill(2, 2, 2) fill(3, 2, 2)] ==
    [D_222 dss(fill(1.5, 2, 2),Ts = 0.005); dss(fill(2, 2, 2),Ts = 0.005) dss(fill(3, 2, 2),Ts = 0.005)]

end # standard state space

@testset "Descriptor state space" begin
# CONTINUOUS
C_111 = dss([1], [0], [2], [3], [4])
C_211 = dss(eye(2), [1 0; 0 0], [1; 2], [1 0], [0])
C_212 = dss(eye(2), I, [1; 2], eye(2), [0; 0])
C_221 = dss(eye(2), zeros(2,2), [1 0; 0 2], [1 0], [0 0])
C_222 = dss(eye(2), 2*eye(2), [1 0; 0 2], eye(2), zeros(Int,2,2))
C_022 = dss(4*eye(2))
CS_211 = dss(eye(2), [1; 2], [1 0], [0])

# DISCRETE
D_111 = dss([1], [0], [2], [3], [4], Ts = 0.005)
D_211 = dss(eye(2), [1 0; 0 0], [1; 2], [1 0], [0], Ts = 0.005)
D_212 = dss(eye(2), I, [1;2], eye(2), [0; 0], Ts = 0.005)
D_221 = dss(eye(2), zeros(2,2), [1 0; 0 2], [1 0], [0 0], Ts = 0.005)
D_222 = dss(eye(2), 2*eye(2), [1 0; 0 2], eye(2), zeros(Int,2,2), Ts = 0.005)
D_022 = dss(4*eye(2), Ts = 0.005)
DS_211 = dss(eye(2), [1; 2], [1 0], [0], Ts = 0.005)

@test [C_111 C_221] == dss(eye(3), zeros(3,3), [2 0 0; 0 1 0; 0 0 2], [3 1 0], [4 0 0])
@test [C_111; C_212] == dss(eye(3), [ 0 0 0 ; 0 1 0; 0 0 1], [2; 1; 2], [3 0 0; 0 1 0; 0 0 1], [4; 0; 0])
@test append(C_111, C_211) == dss(eye(3), [ 0 0 0 ; 0 1 0; 0 0 0], [2 0; 0 1; 0 2], [3 0 0; 0 1 0], [4 0; 0 0])
@test append(C_111, CS_211) == dss(eye(3), [ 0 0 0 ; 0 1 0; 0 0 1], [2 0; 0 1; 0 2], [3 0 0; 0 1 0], [4 0; 0 0])
@test append(I, C_211) == dss(eye(2), [1 0; 0 0], [0 1; 0 2], [0 0; 1 0], [1 0; 0 0])
@test [C_022 C_222] == dss(eye(2), 2*eye(2), [0 0 1 0; 0 0 0 2], [1 0; 0 1], [4 0 0 0; 0 4 0 0])
@test [C_211 CS_211] == dss(eye(4), [1 0 0 0;0 0 0 0; 0 0 1 0; 0 0 0 1], [1 0;2 0;0 1;0 2], [1 0 1 0], [0 0])
@test [C_211 I] == dss(eye(2), [1 0 ;0 0], [1 0;2 0], [1 0], [0 1])
@test [I C_211] == dss(eye(2), [1 0 ;0 0], [0 1;0 2], [1 0], [1 0])

@test horzcat(C_022, C_222) == dss(eye(2), 2*eye(2), [0 0 1 0; 0 0 0 2], [1 0; 0 1], [4 0 0 0; 0 4 0 0])
@test [C_022; C_222] == dss(eye(2), 2*eye(2), [1 0; 0 2], [0 0; 0 0; 1 0; 0 1], [4 0; 0 4; 0 0; 0 0])
@test vertcat(C_022, C_222) == dss(eye(2), 2*eye(2), [1 0; 0 2], [0 0; 0 0; 1 0; 0 1], [4 0; 0 4; 0 0; 0 0])

@test [D_111 D_221] == dss(eye(3), zeros(3,3), [2 0 0; 0 1 0; 0 0 2], [3 1 0], [4 0 0], Ts = 0.005)
@test [D_111; D_212] == dss(eye(3), [ 0 0 0 ; 0 1 0; 0 0 1], [2; 1; 2], [3 0 0; 0 1 0; 0 0 1], [4; 0; 0], Ts = 0.005)
@test append(D_111, D_211) == dss(eye(3), [ 0 0 0 ; 0 1 0; 0 0 0], [2 0; 0 1; 0 2], [3 0 0; 0 1 0], [4 0; 0 0], Ts = 0.005)
@test append(D_111, DS_211) == dss(eye(3), [ 0 0 0 ; 0 1 0; 0 0 1], [2 0; 0 1; 0 2], [3 0 0; 0 1 0], [4 0; 0 0], Ts = 0.005)
@test [D_022 D_222] == dss(eye(2), 2*eye(2), [0 0 1 0; 0 0 0 2], [1 0; 0 1], [4 0 0 0; 0 4 0 0], Ts = 0.005)
@test horzcat(D_022, D_222) == dss(eye(2), 2*eye(2), [0 0 1 0; 0 0 0 2], [1 0; 0 1], [4 0 0 0; 0 4 0 0], Ts = 0.005)
@test [D_022; D_222] == dss(eye(2),2*eye(2),  [1 0; 0 2], [0 0; 0 0; 1 0; 0 1], [4 0; 0 4; 0 0; 0 0], Ts = 0.005)
@test vertcat(D_022, D_222) == dss(eye(2), 2*eye(2), [1 0; 0 2], [0 0; 0 0; 1 0; 0 1], [4 0; 0 4; 0 0; 0 0], Ts = 0.005)
@test [C_211; CS_211] == dss(eye(4), [1 0 0 0;0 0 0 0; 0 0 1 0; 0 0 0 1], [1;2;1;2], [1 0 0 0; 0 0 1 0], [0; 0])
@test [C_211; I] == dss(eye(2), [1 0 ;0 0], [1;2], [1 0; 0 0], [0; 1])
@test [I; C_211] == dss(eye(2), [1 0 ;0 0], [1;2], [0 0; 1 0], [1; 0])

@test series(C_111, C_212) == C_212*C_111
@test series(C_222, [1 1]) == dss([1 1])*C_222
@test series([1, 1],C_222) == C_222*dss([1, 1])
@test parallel(C_111, C_211) == C_111 + C_211
@test parallel([1,2], C_212) == dss([1,2]) + C_212
@test parallel(C_212,[1,2]) == C_212 + dss([1,2])

# Errors
@test_throws ErrorException [C_111 D_111]                 # Sampling time mismatch
@test_throws ErrorException [C_111; D_111]                # Sampling time mismatch
@test_throws ErrorException append(C_111, D_111)          # Sampling time mismatch
@test_throws ErrorException [C_111 C_212]                 # Dimension mismatch
@test_throws ErrorException [C_111; C_221]                # Dimension mismatch

# hcat and vcat for Descriptor StateSpace and Matrix
A = [-1.1 -1.2; -1.3 -1.4]
E = [-1. -2.; 0. -4]
B = [1 2; 3 4]
B = [1 2; 3 4]
C = [5 6; 7 8]
D = [1 0; 0 1]
P = dss(A, E, B, C, D)
@test [P fill(2.5, 2, 1)] == dss(A, E, [B fill(0, 2, 1)], C, [D fill(2.5, 2, 1)])
@test [fill(2.5, 2, 1) P] == dss(A, E, [fill(0, 2, 1) B], C, [fill(2.5, 2, 1) D])
@test [P; fill(2.5, 1, 2)] == dss(A, E, B, [C; fill(0, 1, 2)], [D; fill(2.5, 1, 2)])
@test [fill(2.5, 1, 2); P] == dss(A, E, B, [fill(0, 1, 2); C], [fill(2.5, 1, 2); D])

# hcat and vcat for Descriptor StateSpace and Number
P = dss(-1.0, 0.0, 2.0, 3.0, 4.0)
@test [P 2.5] == dss(-1.0, 0.0, [2.0 0.0], 3.0, [4.0 2.5])
@test [2.5 P] == dss(-1.0, 0.0, [0.0 2.0], 3.0, [2.5 4.0])
@test [P; 2.5] == dss(-1.0, 0.0, 2.0, [3.0; 0.0], [4.0; 2.5])
@test [2.5; P] == dss(-1.0, 0.0, 2.0, [0.0; 3.0], [2.5; 4.0])

#
@test [2.5I P 3.5] == dss(-1.0, 0.0, [0.0 2.0 0.0], 3.0, [2.5 4.0 3.5])
@test [2.5; P; 3.5I] == dss(-1.0, 0.0, 2.0, [0.0; 3.0; 0.0], [2.5; 4.0; 3.5])

# Concatenation of discrete system with constant
@test [D_111 1.] == dss([1.0], [0], [2.0 0.0], [3.0], [4.0 1.0], Ts = 0.005)
@test [1. D_111] == dss([1.0], [0], [0.0 2.0], [3.0], [1.0 4.0], Ts = 0.005)
@test [D_111 I] == dss([1.0], [0], [2.0 0.0], [3.0], [4.0 1.0], Ts = 0.005)
@test [2I D_111] == dss([1.0], [0], [0.0 2.0], [3.0], [2.0 4.0], Ts = 0.005)
# Type and sample time
@test [D_111 1.] isa DescriptorStateSpace{Float64}
@test [D_111 I].Ts == 0.005
# Continuous version
@test [C_111 1.] == dss([1.0], [0], [2.0 0.0], [3.0], [4.0 1.0])
@test [1. C_111] == dss([1.0], [0], [0.0 2.0], [3.0], [1.0 4.0])
@test [C_111 1.] isa DescriptorStateSpace{Float64}
@test [C_111 5.0I].Ts == 0.0

# Concatenation of discrete system with matrix
@test [D_222 fill(1.5, 2, 2)] == [D_222 dss(fill(1.5, 2, 2),Ts = 0.005)]
@test [C_222 fill(1.5, 2, 2)] == [C_222 dss(fill(1.5, 2, 2))]

# hvcat numbers
@test [C_111 1.5; 2 3I] ==
    [C_111 dss(1.5); dss(2.0) dss(3.0)]
@test [D_111 1.5; 2 3] ==
    [D_111 dss(1.5,Ts = 0.005); dss(2.0,Ts = 0.005) dss(3.0,Ts = 0.005)]
# hvcat UniformScaling
@test [C_111 I; 3I] == [C_111 dss(eye(Float64,1)); dss(3*eye(Float64,2,2))]
@test [C_111 I; I 3I] == [C_111 dss(1); dss(1) dss(3*eye(1))]

# hvcat matrices
@test [C_222 fill(1.5, 2, 2); fill(2, 2, 2) fill(3, 2, 2)] ==
    [C_222 dss(fill(1.5, 2, 2)); dss(fill(2, 2, 2)) dss(fill(3, 2, 2))]
@test [D_222 fill(1.5, 2, 2); fill(2, 2, 2) fill(3, 2, 2)] ==
    [D_222 dss(fill(1.5, 2, 2),Ts = 0.005); dss(fill(2, 2, 2),Ts = 0.005) dss(fill(3, 2, 2),Ts = 0.005)]

# type piracy tests
hcat(rand(1,2), [1], I) #OK
hcat(rand(1,2), 1, I)   #NOK
hcat(1, 1, I)           #NOK
vcat(I,[1, 's'])        #OK
vcat([1],[1, 's'])      #OK
vcat(1,[1, 's'])        #NOK
[1 1; I]                #NOK
[[1 1] [1 's']]         #OK
[[1 1] [1 's'];I]       #OK
[1 1 [1 's'];I]         #NOK

end # descriptor state space

end
end # module


text:
module TestMPNModel
  using Test
  using EcologicalNetworks

  @test typeof(mpnmodel(10, 0.3, 0.2)) <: UnipartiteNetwork

  @test richness(mpnmodel(10, 0.3, 0.2)) == 10
  @test richness(mpnmodel(10, 0.4, 0.5)) == 10
  @test richness(mpnmodel(20, 0.35, 0.9)) == 20

  # Need to design test, and nh model doesnt have any argument error yet.
  # @test_throws ArgumentError mnpmodel(10, 99)
  # @test_throws ArgumentError mnpmodel(10, 51)
  # @test_throws ArgumentError mnpmodel(10, 50)
  # @test_throws ArgumentError mnpmodel(10, 200)
  # @test_throws ArgumentError mnpmodel(10, 0)

  #empirical_foodweb = EcologicalNetworks.nz_stream_foodweb()[1]
  #@test typeof(nichemodel(empirical_foodweb)) <: UnipartiteNetwork

end


text:
using SafeTestsets
@safetestset "My f tests" begin include("my_f_test.jl")
end
# 2x+3y


text:
using CompileBot

botconfig = BotConfig(
    "Plots",
    version = ["1.6", "1.7", "1.8", "nightly"],  # <<< keep these versions in sync with .github/workflows/SnoopCompile.yml
    # else_version = "nightly",
)


text:
"""
    Zoom <: Augmentor.ImageOperation

Description
--------------

Scales the image height and image width by the specified
`factors`, but crops the image such that the original size is
preserved.

The provided `factors` can either be numbers or vectors of
numbers.

- If numbers are provided, then the operation is deterministic
  and will always scale the input image with the same factors.

- In the case vectors are provided, then each time the operation
  is applied a valid index is sampled and the elements
  corresponding to that index are used as scaling factors.

In contrast to [`Scale`](@ref) the size of the output image is
the same as the size of the input image, while the content is
scaled the same way. The same effect could be achieved by
following a [`Scale`](@ref) with a [`CropSize`](@ref), with the
caveat that one would need to know the exact size of the input
image before-hand.

Usage
--------------

    Zoom(factors)

    Zoom(factors...)

Arguments
--------------

- **`factors`** : `NTuple` or `Vararg` of `Real` or
    `AbstractVector` that denote the scale factor(s) for each
    array dimension. If only one variable is specified it is
    assumed that height and width should be scaled by the same
    factor(s).

See also
--------------

[`Scale`](@ref), [`Resize`](@ref), [`augment`](@ref)

Examples
--------------

```julia
using Augmentor
img = testpattern()

# half the image size
augment(img, Zoom(0.5))

# uniformly scale by a random factor from 1.2, 1.3, or 1.4
augment(img, Zoom([1.2, 1.3, 1.4]))

# scale by either 0.5x0.7 or by 0.6x0.8
augment(img, Zoom([0.5, 0.6], [0.7, 0.8]))
```
"""
struct Zoom{N,T<:AbstractVector} <: ImageOperation
    factors::NTuple{N,T}

    function Zoom{N}(factors::NTuple{N,T}) where {N,T<:AbstractVector}
        eltype(T) <: Real || throw(ArgumentError("The specified factors must be vectors of Real. Actual: $T"))
        n = length(factors[1])
        n > 0 || throw(ArgumentError("The specified factors must all have a length greater than 0"))
        all(f->length(f)==n, factors) || throw(ArgumentError("The specified factors must all be of the same length for each dimension. Actual: $(factors)"))
        new{N,T}(factors)
    end
end
Zoom() = throw(MethodError(Zoom, ()))
Zoom(::Tuple{}) = throw(MethodError(Zoom, ((),)))
Zoom(factors...) = Zoom(factors)
Zoom(factor::Union{AbstractVector,Real}) = Zoom((factor, factor))
Zoom(factors::NTuple{N,Any}) where {N} = Zoom(map(vectorize, factors))
Zoom(factors::NTuple{N,Range}) where {N} = Zoom{N}(promote(factors...))
function Zoom(factors::NTuple{N,AbstractVector}) where N
    Zoom{N}(map(Vector{Float64}, factors))
end
function (::Type{Zoom{N}})(factors::NTuple{N,Any}) where N
    Zoom(map(vectorize, factors))
end

@inline supports_affineview(::Type{<:Zoom}) = true
@inline supports_eager(::Type{<:Zoom}) = false

randparam(op::Zoom, imgs::Tuple) = randparam(op, imgs[1])

function randparam(op::Zoom, img::AbstractArray{T,N}) where {T,N}
    i = safe_rand(1:length(op.factors[1]))
    ntuple(j -> Float64(op.factors[j][i]), Val{N})
end

function toaffinemap(op::Zoom{2}, img::AbstractMatrix, idx)
    @inbounds tfm = recenter(@SMatrix([Float64(idx[1]) 0.; 0. Float64(idx[2])]), center(img))
    tfm
end

function applylazy(op::Zoom, img::AbstractArray, idx)
    applyaffineview(op, prepareaffine(img), idx)
end

function applyaffineview(op::Zoom{N}, img::AbstractArray{T,N}, idx) where {T,N}
    wv = invwarpedview(img, toaffinemap(op, img, idx), indices(img))
    direct_view(wv, indices(img))
end

function applyaffineview(op::Zoom{N}, v::SubArray{T,N,<:InvWarpedView}, idx) where {T,N}
    tinv = toaffinemap(op, v, idx)
    img = parent(v)
    nidx = ImageTransformations.autorange(img, tinv)
    wv = InvWarpedView(img, tinv, map(unionrange, nidx, indices(img)))
    view(wv, v.indexes...)
end

function showconstruction(io::IO, op::Zoom)
    fct = length(op.factors[1]) == 1 ? map(first,op.factors) : op.factors
    print(io, typeof(op).name.name, '(', join(map(string, fct),", "), ')')
end

function Base.show(io::IO, op::Zoom{N}) where N
    if get(io, :compact, false)
        str = join(map(t->join(round_if_float(t,2),"×"), collect(zip(op.factors...))), ", ")
        if length(op.factors[1]) == 1
            print(io, "Zoom by $(str)")
        else
            print(io, "Zoom by I ∈ {$(str)}")
        end
    else
        fct = length(op.factors[1]) == 1 ? map(first,op.factors) : op.factors
        print(io, typeof(op).name, "{$N}($(fct))")
    end
end


text:
# Exercise 17-1
struct MyTime
    hour :: Int64
    minute :: Int64
    second :: Int64
end

function timetoint(time::MyTime)
    minutes = time.hour * 60 + time.minute
    seconds = minutes * 60 + time.second
end

function inttotime(seconds::Int64)
    (minutes, second) = divrem(seconds, 60)
    hour, minute = divrem(minutes, 60)
    MyTime(hour, minute, second)
end

# Exercise 17-2
mutable struct Point
    x :: Number
    y :: Number
end

Point(a) = Point(a, 0)
Point() = Point(0, 0)

# println(Point())
# println(Point(7.2))
# println(Point(2.5, 3))
# println(Point(2.5, "abc"))

# Exercise 17-3
import Base.+

function +(p1::Point, p2::Point)
    return Point(p1.x + p2.x, p1.y + p2.y)
end

function +(p1::Tuple, p2::Point)
    return Point(p1[1] + p2.x, p1[2] + p2.y)
end

function +(p1::Point, p2::Tuple)
    return Point(p1.x + p2[1], p1.y + p2[2])
end

a = Point(1, 2.0)
b = Point(-1.5, 3.7)
c = (6.5, -4.5)
# println(a + b)
# println(a + c)
# println(c + b)

# Exercise 17-4
using Printf
import Base.+

struct MyTimeMidnight
    seconds :: Int64
end

function Base.show(io::IO, time::MyTimeMidnight)
    (minutes, second) = divrem(time.seconds, 60)
    hour, minute = divrem(minutes, 60)
    @printf(io, "%02d:%02d:%02d", hour, minute, second)
end

function increment(time::MyTimeMidnight, seconds::Int64)
    return MyTimeMidnight(time.seconds + seconds)
    # return inttotime(timetoint(time) + seconds)
end

function +(t1::MyTimeMidnight, t2::MyTimeMidnight)
    return MyTimeMidnight(t1.seconds + t2.seconds)
end

function +(time::MyTimeMidnight, seconds::Int64)
    increment(time, seconds)
end

function +(seconds::Int64, time::MyTimeMidnight)
    time + seconds
end

function ex17_4()
    t1 = MyTimeMidnight(54321)
    t2 = MyTimeMidnight(10000)
    println(t1)
    println(t2)
    println(t1 + t2)
    println(t1 + 1000)
    println(2000 + t1)
end

# ex17_4()

# Exercise 17-5
struct Kangaroo
    pouchcontents :: Array
    function Kangaroo()
        new([])
    end
end

function putinpouch(k::Kangaroo, obj)
    push!(k.pouchcontents, obj)
end

function Base.show(io::IO, k::Kangaroo)
    println(io, "$(k.pouchcontents)")
end

function ex17_5()
    kanga = Kangaroo()
    roo = Kangaroo()
    putinpouch(roo, "abc")
    putinpouch(roo, 123)
    putinpouch(kanga, roo)
    putinpouch(kanga, (1.2, 2.3, "123"))
    println(kanga)
end

# ex17_5()


text:
"""
    space(::SiteType"Electron";
          conserve_qns = false,
          conserve_sz = conserve_qns,
          conserve_nf = conserve_qns,
          conserve_nfparity = conserve_qns,
          qnname_sz = "Sz",
          qnname_nf = "Nf",
          qnname_nfparity = "NfParity")

Create the Hilbert space for a site of type "Electron".

Optionally specify the conserved symmetries and their quantum number labels.
"""
function ITensors.space(
  ::SiteType"Electron";
  conserve_qns=false,
  conserve_sz=conserve_qns,
  conserve_nf=conserve_qns,
  conserve_nfparity=conserve_qns,
  qnname_sz="Sz",
  qnname_nf="Nf",
  qnname_nfparity="NfParity",
  # Deprecated
  conserve_parity=nothing,
)
  if !isnothing(conserve_parity)
    conserve_nfparity = conserve_parity
  end
  if conserve_sz && conserve_nf
    return [
      QN((qnname_nf, 0, -1), (qnname_sz, 0)) => 1
      QN((qnname_nf, 1, -1), (qnname_sz, +1)) => 1
      QN((qnname_nf, 1, -1), (qnname_sz, -1)) => 1
      QN((qnname_nf, 2, -1), (qnname_sz, 0)) => 1
    ]
  elseif conserve_nf
    return [
      QN(qnname_nf, 0, -1) => 1
      QN(qnname_nf, 1, -1) => 2
      QN(qnname_nf, 2, -1) => 1
    ]
  elseif conserve_sz
    return [
      QN((qnname_sz, 0), (qnname_nfparity, 0, -2)) => 1
      QN((qnname_sz, +1), (qnname_nfparity, 1, -2)) => 1
      QN((qnname_sz, -1), (qnname_nfparity, 1, -2)) => 1
      QN((qnname_sz, 0), (qnname_nfparity, 0, -2)) => 1
    ]
  elseif conserve_nfparity
    return [
      QN(qnname_nfparity, 0, -2) => 1
      QN(qnname_nfparity, 1, -2) => 2
      QN(qnname_nfparity, 0, -2) => 1
    ]
  end
  return 4
end

ITensors.val(::ValName"Emp", ::SiteType"Electron") = 1
ITensors.val(::ValName"Up", ::SiteType"Electron") = 2
ITensors.val(::ValName"Dn", ::SiteType"Electron") = 3
ITensors.val(::ValName"UpDn", ::SiteType"Electron") = 4
ITensors.val(::ValName"0", st::SiteType"Electron") = val(ValName("Emp"), st)
ITensors.val(::ValName"↑", st::SiteType"Electron") = val(ValName("Up"), st)
ITensors.val(::ValName"↓", st::SiteType"Electron") = val(ValName("Dn"), st)
ITensors.val(::ValName"↑↓", st::SiteType"Electron") = val(ValName("UpDn"), st)

ITensors.state(::StateName"Emp", ::SiteType"Electron") = [1.0, 0, 0, 0]
ITensors.state(::StateName"Up", ::SiteType"Electron") = [0.0, 1, 0, 0]
ITensors.state(::StateName"Dn", ::SiteType"Electron") = [0.0, 0, 1, 0]
ITensors.state(::StateName"UpDn", ::SiteType"Electron") = [0.0, 0, 0, 1]
ITensors.state(::StateName"0", st::SiteType"Electron") = state(StateName("Emp"), st)
ITensors.state(::StateName"↑", st::SiteType"Electron") = state(StateName("Up"), st)
ITensors.state(::StateName"↓", st::SiteType"Electron") = state(StateName("Dn"), st)
ITensors.state(::StateName"↑↓", st::SiteType"Electron") = state(StateName("UpDn"), st)

function ITensors.op!(Op::ITensor, ::OpName"Nup", ::SiteType"Electron", s::Index)
  Op[s' => 2, s => 2] = 1.0
  return Op[s' => 4, s => 4] = 1.0
end

function ITensors.op!(Op::ITensor, ::OpName"Ndn", ::SiteType"Electron", s::Index)
  Op[s' => 3, s => 3] = 1.0
  return Op[s' => 4, s => 4] = 1.0
end

function ITensors.op!(Op::ITensor, ::OpName"Nupdn", ::SiteType"Electron", s::Index)
  return Op[s' => 4, s => 4] = 1.0
end

function ITensors.op!(Op::ITensor, ::OpName"Ntot", ::SiteType"Electron", s::Index)
  Op[s' => 2, s => 2] = 1.0
  Op[s' => 3, s => 3] = 1.0
  return Op[s' => 4, s => 4] = 2.0
end

function ITensors.op!(Op::ITensor, ::OpName"Cup", ::SiteType"Electron", s::Index)
  Op[s' => 1, s => 2] = 1.0
  return Op[s' => 3, s => 4] = 1.0
end

function ITensors.op!(Op::ITensor, ::OpName"Cdagup", ::SiteType"Electron", s::Index)
  Op[s' => 2, s => 1] = 1.0
  return Op[s' => 4, s => 3] = 1.0
end

function ITensors.op!(Op::ITensor, ::OpName"Cdn", ::SiteType"Electron", s::Index)
  Op[s' => 1, s => 3] = 1.0
  return Op[s' => 2, s => 4] = -1.0
end

function ITensors.op!(Op::ITensor, ::OpName"Cdagdn", ::SiteType"Electron", s::Index)
  Op[s' => 3, s => 1] = 1.0
  return Op[s' => 4, s => 2] = -1.0
end

function ITensors.op!(Op::ITensor, ::OpName"Aup", ::SiteType"Electron", s::Index)
  Op[s' => 1, s => 2] = 1.0
  return Op[s' => 3, s => 4] = 1.0
end

function ITensors.op!(Op::ITensor, ::OpName"Adagup", ::SiteType"Electron", s::Index)
  Op[s' => 2, s => 1] = 1.0
  return Op[s' => 4, s => 3] = 1.0
end

function ITensors.op!(Op::ITensor, ::OpName"Adn", ::SiteType"Electron", s::Index)
  Op[s' => 1, s => 3] = 1.0
  return Op[s' => 2, s => 4] = 1.0
end

function ITensors.op!(Op::ITensor, ::OpName"Adagdn", ::SiteType"Electron", s::Index)
  Op[s' => 3, s => 1] = 1.0
  return Op[s' => 4, s => 2] = 1.0
end

function ITensors.op!(Op::ITensor, ::OpName"F", ::SiteType"Electron", s::Index)
  Op[s' => 1, s => 1] = +1.0
  Op[s' => 2, s => 2] = -1.0
  Op[s' => 3, s => 3] = -1.0
  return Op[s' => 4, s => 4] = +1.0
end

function ITensors.op!(Op::ITensor, ::OpName"Fup", ::SiteType"Electron", s::Index)
  Op[s' => 1, s => 1] = +1.0
  Op[s' => 2, s => 2] = -1.0
  Op[s' => 3, s => 3] = +1.0
  return Op[s' => 4, s => 4] = -1.0
end

function ITensors.op!(Op::ITensor, ::OpName"Fdn", ::SiteType"Electron", s::Index)
  Op[s' => 1, s => 1] = +1.0
  Op[s' => 2, s => 2] = +1.0
  Op[s' => 3, s => 3] = -1.0
  return Op[s' => 4, s => 4] = -1.0
end

function ITensors.op!(Op::ITensor, ::OpName"Sz", ::SiteType"Electron", s::Index)
  Op[s' => 2, s => 2] = +0.5
  return Op[s' => 3, s => 3] = -0.5
end

function ITensors.op!(Op::ITensor, ::OpName"Sᶻ", st::SiteType"Electron", s::Index)
  return op!(Op, OpName("Sz"), st, s)
end

function ITensors.op!(Op::ITensor, ::OpName"Sx", ::SiteType"Electron", s::Index)
  Op[s' => 2, s => 3] = 0.5
  return Op[s' => 3, s => 2] = 0.5
end

function ITensors.op!(Op::ITensor, ::OpName"Sˣ", st::SiteType"Electron", s::Index)
  return op!(Op, OpName("Sx"), st, s)
end

function ITensors.op!(Op::ITensor, ::OpName"S+", ::SiteType"Electron", s::Index)
  return Op[s' => 2, s => 3] = 1.0
end

function ITensors.op!(Op::ITensor, ::OpName"S⁺", st::SiteType"Electron", s::Index)
  return op!(Op, OpName("S+"), st, s)
end
function ITensors.op!(Op::ITensor, ::OpName"Sp", st::SiteType"Electron", s::Index)
  return op!(Op, OpName("S+"), st, s)
end
function ITensors.op!(Op::ITensor, ::OpName"Splus", st::SiteType"Electron", s::Index)
  return op!(Op, OpName("S+"), st, s)
end

function ITensors.op!(Op::ITensor, ::OpName"S-", ::SiteType"Electron", s::Index)
  return Op[s' => 3, s => 2] = 1.0
end

function ITensors.op!(Op::ITensor, ::OpName"S⁻", st::SiteType"Electron", s::Index)
  return op!(Op, OpName("S-"), st, s)
end
function ITensors.op!(Op::ITensor, ::OpName"Sm", st::SiteType"Electron", s::Index)
  return op!(Op, OpName("S-"), st, s)
end
function ITensors.op!(Op::ITensor, ::OpName"Sminus", st::SiteType"Electron", s::Index)
  return op!(Op, OpName("S-"), st, s)
end

ITensors.has_fermion_string(::OpName"Cup", ::SiteType"Electron") = true
ITensors.has_fermion_string(::OpName"Cdagup", ::SiteType"Electron") = true
ITensors.has_fermion_string(::OpName"Cdn", ::SiteType"Electron") = true
ITensors.has_fermion_string(::OpName"Cdagdn", ::SiteType"Electron") = true


text:
export project_cardinality!

function project_cardinality!(
                             x::Union{Vector{TF},Vector{Complex{TF}}},
                             k::TI,
                             ) where {TF<:Real,TI<:Integer}
  """
  Project the vector x onto the set of vectors with cardinality (l0 'norm') less then or equal to k.
  project m onto {m | card(m) <= k} : x = argmin_x 1/2 ||x-m||2^2 s.t. card(x)<=k
  """

#sort_ind = sortperm( a, by=abs, rev=true)
#val_max=x[sort_ind[1:k]]
#x=zeros(length(x))
#x[sort_ind[1:k]]=val_max

#alternative
sort_ind = sortperm( x, by=abs, rev=true)
x[sort_ind[k+1:end]] .= 0.0
  return x
end

function project_cardinality!(
                             x    ::Array{TF,2},  #matrix, project each row or colum
                             k    ::TI,               #maximum cardinality
                             mode ::Tuple{String,String},    #
                             return_vec=true ::Bool
                             ) where {TF<:Real,TI<:Integer}
  """
  Project each of the columns/rows of the model in matrix form onto the set of vectors with cardinality (l0 'norm') less then or equal to k.
  project m onto {m | card(m)<= k} : x = argmin_x 1/2 ||x-m||2^2 s.t. card(x)<=k
  """

#sort_ind = sortperm( a, by=abs, rev=true)
#val_max=x[sort_ind[1:k]]
#x=zeros(length(x))
#x[sort_ind[1:k]]=val_max

#alternative
if mode[1] == "fiber"
  if mode[2] == "x"
  Threads.@threads for i=1:size(x,2)
      sort_ind = sortperm( view(x,:,i), by=abs, rev=true)
  @inbounds x[sort_ind[k+1:end],i] .= 0.0
    end
  elseif mode[2] == "z"
  Threads.@threads for i=1:size(x,1)
      sort_ind = sortperm( view(x,i,:), by=abs, rev=true)
  @inbounds x[i,sort_ind[k+1:end]] .= 0.0
    end
  end
else
  error("for 2D models, the mode of application for project_cardinality! needs to be (fiber,x) or (fiber,y). Or, provide the model as a vector")
end

if return_vec==true
  return vec(x)
else
  return x
end

end

function project_cardinality!(
                             x    ::Array{TF,3},
                             k    ::TI,
                             mode ::Tuple{String,String},
                             return_vec=true ::Bool
                             ) where {TF<:Real,TI<:Integer}
  """
  Project each of the x or y or z fibers of the model in tensor form onto the set of vectors with cardinality (l0 'norm') less then or equal to k.
  project m onto {m | card(m)<= k} : x = argmin_x 1/2 ||x-m||2^2 s.t. card(x)<=k
  """

(n1,n2,n3) = deepcopy(size(x))
#sort_ind = sortperm( a, by=abs, rev=true)
#val_max=x[sort_ind[1:k]]
#x=zeros(length(x))
#x[sort_ind[1:k]]=val_max

#alternative

#Fiber based projection for 3D tensor
if mode[1] == "fiber"
  if mode[2] == "x"
    for i=1:n2
    Threads.@threads for j=1:n3
          #sort_ind = sortperm( x[:,i], by=abs, rev=true)
          sort_ind = sortperm( view(x,:,i,j), by=abs, rev=true)
  @inbounds x[sort_ind[k+1:end],i,j] .= 0.0
        end
    end
  elseif mode[2] == "z"
    Threads.@threads for i=1:n1
    for j=1:n2
        #sort_ind = sortperm( x[:,i], by=abs, rev=true)
        sort_ind = sortperm( view(x,i,j,:), by=abs, rev=true)
  @inbounds x[i,j,sort_ind[k+1:end]] .= 0.0
      end
    end
  elseif mode[2] == "y"
    for i=1:n1
    Threads.@threads for j=1:n3
        #sort_ind = sortperm( x[:,i], by=abs, rev=true)
        sort_ind = sortperm( view(x,i,:,j), by=abs, rev=true)
  @inbounds x[i,sort_ind[k+1:end],j] .= 0.0
      end
    end
  end

elseif mode[1] == "slice" #Slice based projection for 3D tensor
  #code currently does not mutate the input for slice projections for 3D tensors...
  #permute and reshape
  if mode[2] == "x"
    x = permutedims(x,[2;3;1])
    x = reshape(x,n2*n3,n1)
  elseif mode[2] == "y"
    x = permutedims(x,[1;3;2])
    x = reshape(x,n1*n3,n2)
  elseif mode[2] == "z"
    x = reshape(x,n1*n2,n3) #there are n3 slices of dimension n1*n2
  end

  #project, same for all modes because we permuted and reshaped already
  Threads.@threads for i=1:size(x,2)
    sort_ind = sortperm( view(x,:,i), by=abs, rev=true)
    @inbounds x[sort_ind[k+1:end],i] .= 0.0
  end

  #reverse reshape and permute back
  if mode[2] == "x"
    x = reshape(x,n2,n3,n1)
    x = permutedims(x,[3;1;2]);
  elseif mode[2] == "y"
    x = reshape(x,n1,n3,n2)
    x = permutedims(x,[1;3;2]);
  elseif mode[2] == "z"
    x = reshape(x,n1,n2,n3) #there are n3 slices of dimension n1*n2
  end
end #if slice/fiber mode

if return_vec==true
  x=vec(x)
end
return x
end


text:
using MakieLayout
using Makie

begin
    scene = Scene(resolution = (1000, 1000));
    screen = display(scene)
    campixel!(scene);

    maingl = GridLayout(scene, 1, 3;
        addedcolgaps = Fixed(0),
        alignmode = Outside(30)
    )

    maingl[1, 2] = LAxis(scene)
    maingl[1, 3] = LAxis(scene)

    guigl = maingl[1, 1] = GridLayout(1, 3; height=Auto(false))

    # guigl[1, 1] = LButton(scene, 100, 30, "button")
    guigl[1, 1] = LText(scene, text="HelloWorld", halign=:left)
    guigl[1, 2] = LText(scene, text="Blablo", halign=:left)
    guigl[1, 3] = LButton(scene, width = 120, height = 30, label = "Bliblu")

    guigl[2, 1] = LText(scene, text="Mamama", halign=:left)
    guigl[2, 2] = LText(scene, text="Mimimi", halign=:left)
    guigl[2, 3] = LButton(scene, width = 120, height = 30, label = "Momomo")

    guigl[3, 1] = LText(scene, text="VeerrryLoong", halign=:left)
    guigl[3, 2] = LText(scene, text="Blablo", halign=:left)
    guigl[3, 3] = LButton(scene, width = 120, height = 30, label = "Short")

end


text:
#=
Copyright (c) 2015, Intel Corporation

Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions are met:

    * Redistributions of source code must retain the above copyright notice,
      this list of conditions and the following disclaimer.
    * Redistributions in binary form must reproduce the above copyright
      notice, this list of conditions and the following disclaimer in the
      documentation and/or other materials provided with the distribution.
    * Neither the name of Intel Corporation nor the names of its contributors
      may be used to endorse or promote products derived from this software
      without specific prior written permission.

THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE
FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
=#

export MemoryDataLayer

type MemoryDataEnsemble{N,M} <: DataEnsemble
    name         :: Symbol
    neurons      :: Array{DataNeuron,N}
    value        :: Array{Float32, M}  # M != N because of batch dimension
    connections  :: Vector{Connection}
    phase        :: Phase
    net_subgroup :: Cint
end

function MemoryDataEnsemble{N, M}(name::Symbol, neurons::Array{DataNeuron,N}, value::Array{Float32, M}, phase::Phase)
    MemoryDataEnsemble{N,M}(name, neurons, value, Connection[], phase, convert(Cint, 1))
end

function forward{N}(ens::MemoryDataEnsemble, data::Array{Float32,N}, net::Net, phase::Phase)
    if net.time_steps > 1
        data[:] = ens.value[[Colon() for _ in ndims(ens)]..., :, net.curr_time_step]
    else
        data[:] = ens.value[:]
    end
end

function backward{N}(ens::MemoryDataEnsemble, data::Array{Float32,N}, net::Net, phase::Phase)
end

function MemoryDataLayer(net::Net, name::Symbol, shape::Tuple; phase=TrainTest)
    data_neurons = Array(DataNeuron, shape...)
    for i in 1:length(data_neurons)
        data_neurons[i] = DataNeuron(0.0)
    end
    shape = [shape...]
    push!(shape, batch_size(net))
    if net.time_steps > 1
        push!(shape, net.time_steps)
    end
    value = Array(Float32, shape...)
    ens = MemoryDataEnsemble(name, data_neurons, value, phase)
    add_ensemble(net, ens)
    ens, value
end


text:
# -----------------------------------------------------------
#
# Evaluating the 1D Galerkin Basis functions
#
# -----------------------------------------------------------

const K_max = 10

# Efficiency criticality: HIGH
# Central to the efficiency of the package

# Performing fastmath on an array of coefficients
# of x^k and |x|^k, giving the value at a specified point
# supported on [-1, 1]
@fastmath function array2poly(v::AbstractArray{T, 1}, x::Real) where T <: Real
    if abs(x)>1
        return zero(T)
    end
    n = length(v)
    k = div(n,2)
    s = zero(T)
    @inbounds for i in k:-1:1
        # Using Horner's method
        s *= x
        s += v[i] + flipsign(v[i+k], x)
    end
    return s
end

function array2poly(v::AbstractArray{T, 1}) where T <: Real
    return (x-> array2poly(v,x))
end

#precompute Legendre
leg_coeffs = legendre(K_max)

# Legendre polynomial supported on [-1, 1]
function LegendreP(k, x)
    k<=K_max || throw(DomainError())
    return array2poly(leg_coeffs[k+1], x)
end

function LegendreP(k)
    k<=K_max || throw(DomainError())
    return array2poly(leg_coeffs[k+1])
end

# precomputing the DG functions
# (Thanks to Alex Arslan for this.)
# Meaning: [baisis size k][coefficient, mode]
# Lengths: [K_max][2*k, k]
const dg_coeffs = Vector{Matrix{Float64}}(undef, K_max)

for k in 1:K_max
    Q = dg_basis(k)
    @assert length(Q) == k
    @assert all(length(Q[mode]) == 2*k for mode in 1:k)
    dg_coeffs[k] = Matrix{Float64}(undef, 2*k, k)
    for mode in 1:k
        dg_coeffs[k][:,mode] = Q[mode][:]
    end
end

# This is the dg basis function corresponding a given mode
# supported on [-1, 1]
function h(k, mode, x)
    mode<=k || throw(DomainError())
    return array2poly((@view dg_coeffs[k][:,mode]), x)
end

function h(k, mode)
    mode<=k || throw(DomainError())
    return array2poly((@view dg_coeffs[k][:,mode]))
end

# ----------------------------------------------
#
# Methods relating to the 1D Position Basis:
#
# ----------------------------------------------

# Legendre polynomial on [0, 1]
function leg(mode::Int, x::T) where T <: Real
    return sqrt(2.0)*LegendreP(mode-1, 2*x-1)
end

# level >= 1
# cell in 1:1<<level
# mode in 1:k
# Legendre polynomial on [(cell-1)*h, cell*h] with h = 1/(1<<level)
function basis(level::Int, cell::Int, mode::Int, x::T) where T <: Real
    return leg(mode, (1<<level)*x - (cell-1)) * (2.0)^(level/2)
end

function basis(level::Int, cell::Int, mode::Int)
    return x->basis(level, cell, mode, x)
end

# Perform numerical integration to get position basis coeffs
#
# Never used except if user demands to evolve wave equation
# using only the position basis, never going through heir
function pos_vcoeffs_DG(k::Int, level::Int, f::Function;
                        rtol = REL_TOL, atol = ABS_TOL, maxevals = MAX_EVALS)
    vcoeffs = Array{Float64}(undef, (1<<level)*(k))
    i = 1
    for cell in 1:(1<<level)
        for mode in 1:k
            fcn            = x->(basis(level, cell, mode,x)*f(x))
            left_bound  = (cell-1)/(1<<level)
            right_bound = (cell)/(1<<level)
            vcoeffs[i]  = hquadrature(fcn, left_bound, right_bound; atol=atol)[1]
            i += 1
        end
    end
    return vcoeffs
end

# ----------------------------------------------
# Building the hier2pos matrix:
# ----------------------------------------------

# Given an array of the type above, of coefficients
# for both both x^k and |x|^k, then given a choice of
# side (right vs. left) this becomes just a pure
# polynomial array of half the length
function convert_polyarray(v::AbstractArray{T, 1}, side = "left") where T <: Real
    n = Int(length(v)//2)
    if side == "left"
        return [v[i]-v[n+i] for i in 1:n]
    elseif side == "right"
        [v[i]+v[n+i] for i in 1:n]
    else throw(ArgumentError())
    end
end

# Gives the polynomial array corresponding to p(x/c)
function scale_polyarray(v::AbstractArray{T, 1}, c::Real) where T <: Real
    n = length(v)
    return [v[i] * (1//c^(i-1)) for i in 1:n]
end

# Gives the polynomial array corresponding to p(x-a)
function shift_polyarray(v::AbstractArray{T, 1}, a::Real) where T <: Real
    n = length(v)
    v_new = zeros(v)
    for i::Int in 1:n
        for j::Int in 1:i
            v_new[j] += v[i] * (-a)^(i-j) * binomial(i-1, j-1)
        end
    end
    return v_new
end

# Performs symbolic integration of the polynomials
# represented by v and w on the interval [a, b]
function integrate_polyarray(v::AbstractArray{T, 1}, w::AbstractArray{T, 1}; a::Real = 0, b::Real = 1) where T <: Real
    ans = zero(T)
    for i::Int in 1:length(v)
        for j::Int in 1:length(w)
            ans += v[i] * w[j] * ((b^(i+j-1)-a^(i+j-1))//(i+j-1))
        end
    end
    return ans
end

# With a bit of case work, we can do an explicit
# symbolic integration of the elements in the
# hierarchical basis against the legendre polynomials
# defining the position basis:
function hier2pos(k::Int, max_level::Int, level::Int, cell::Int, mode::Int)

    ans = Real[]

    # Compute the support and discontinuity
    # of the hierarchical basis element
    left_point  = (cell-1)//(1<<max(0, level-1))
    mid_point   = (cell- 1//2)//(1<<max(0, level-1))
    right_point = cell//(1<<max(0, level-1))

    # Level 0 of the hierarchical basis is a continous
    # Legendre polynomial
    if level == 0
        vl = vr = sqrt(2)*shift_polyarray(
                                scale_polyarray(leg_coeffs[mode][1:k],
                                                1//2),
                                            1//2)

    # The higher DG levels are given by two polynomials
    # defined on neighboring (right & left) intervals of size h/2,
    # which would need to be integrated against individually
    else
        h = 1//(1<<max(0, level-1))
        vl = shift_polyarray(
                scale_polyarray(
                    convert_polyarray((@view dg_coeffs[k][:,mode]), "left"),
                                1//2),
                            1//2)
        vl = shift_polyarray(scale_polyarray(vl, h), (cell-1)*h)*sqrt(1<<(level))

        vr = shift_polyarray(
                scale_polyarray(
                    convert_polyarray((@view dg_coeffs[k][:,mode]), "right"),
                                 1//2),
                            1//2)
        vr = shift_polyarray(scale_polyarray(vr, h), (cell-1)*h)*sqrt(1<<(level))
    end

    h = 1//(1<<max_level)
    for i in 1:1<<max_level
        for p in 1:k
            # Check if our hier function is supported on this interval
            if (i-1)*h < left_point || i*h > right_point
                push!(ans, 0)
                continue
            end

            # Obtain the polynomial array for the position basis function
            unit_leg = shift_polyarray(
                            scale_polyarray(leg_coeffs[p][1:k], 1//2), 1//2)
            pos_element = shift_polyarray(
                                scale_polyarray(unit_leg, h), (i-1)*h)
            pos_element *= sqrt(1<<(max_level+1)) # (2.0)^((max_level+1)/2)

            # Perform symbolic integration depending on which
            # interval the position basis function overlaps with the
            # DG basis element
            if (i-1)*h < mid_point
                val = integrate_polyarray(vl, pos_element; a=(i-1)*h, b=i*h)
            else
                val = integrate_polyarray(vr, pos_element; a=(i-1)*h, b=i*h)
            end
            push!(ans, val)
        end
    end
    return ans
end

# Construct the change of basis matrix from hierarchical
# to position basis
function hier2pos(k::Int, max_level::Int; atol=ABS_TOL)
    j = 1
    I = Int[]
    J = Int[]
    V = Float64[]
    for level in 0:max_level
        for cell in 1:(1<<max(0, level-1))
            for mode in 1:k
                ans = pos_vcoeffs_DG(k, max_level, v(k, level, cell, mode))
                for i in 1:length(ans)
                    if abs(ans[i]) > atol
                        push!(I, i)
                        push!(J, j)
                        push!(V, ans[i])
                    end
                end
                j += 1
            end
        end
    end
    return sparse(I, J, V, k * (1<<max_level), k * (1<<max_level), +)
end

function pos2hier(k::Int, max_level::Int; atol=ABS_TOL)
    return hier2pos(k, max_level; atol=atol)'
end


text:
@testset "indexing_Primitive_buffer" begin
    for i ∈ 1:N_OUTER
        len, dtype, lpad, b, v = rand_primitive_buffer()
        p = Primitive{dtype}(b, lpad+1, len)
        # integer indices
        for j ∈ 1:N_IDX_CHECK
            k = rand(1:len)
            @test p[k] == v[k]
        end
        # AbstractVector{<:Integer} indices
        for j ∈ 1:N_IDX_CHECK
            idx = rand(1:len, rand(1:MAX_IDX_LEN))
            @test p[idx] == v[idx]
        end
        @test p[:] == v
    end
end

@testset "indexing_Primitive_construct" begin
    for i ∈ 1:N_OUTER
        dtype = rand(PRIMITIVE_ELTYPES)
        len = rand(1:MAX_VECTOR_LENGTH)
        v = rand(dtype, len)
        p = Primitive(v)
        # integer indices
        for j ∈ 1:N_IDX_CHECK
            k = rand(1:len)
            @test p[k] == v[k]
        end
        # AbstractVector{<:Integer} indices
        for j ∈ 1:N_IDX_CHECK
            idx = rand(1:len, rand(1:MAX_IDX_LEN))
            @test p[idx] == v[idx]
        end
        @test p[:] == v
    end
end

@testset "indexing_Primitive_setindex" begin
    for i ∈ 1:N_OUTER
        dtype = rand(PRIMITIVE_ELTYPES)
        len = rand(1:MAX_VECTOR_LENGTH)
        v = rand(dtype, len)
        p = Primitive(v)
        for j ∈ 1:N_IDX_CHECK
            k = rand(1:len)
            x = rand(dtype)
            p[k] = x
            @test p[k] == x
        end
        for j ∈ 1:N_IDX_CHECK
            idx = unique(rand(1:len, rand(1:MAX_IDX_LEN)))
            x = rand(dtype, length(idx))
            p[idx] = x
            @test p[idx] == x
        end
        # contiguous
        for j ∈ 1:N_IDX_CHECK
            a, b = extrema(rand(1:len, 2))
            idx = a:b
            x = rand(dtype, length(idx))
            p[idx] = x
            @test p[idx] == x
        end
    end
end

@testset "indexing_Primitive_empty" begin
    for i ∈ 1:N_OUTER
        dtype = rand(PRIMITIVE_ELTYPES)
        p = Primitive(dtype[])
        @test isempty(p)
        @test_throws BoundsError p[1]
    end
end

@testset "indexing_NullablePrimitive_buffer" begin
    for i ∈ 1:N_OUTER
        len, dtype, bmask, lpad, b, v = rand_nullableprimitive_buffer()
        p = NullablePrimitive{dtype}(b, 1+lpad, 1+bmask+lpad, len)
        # integer indices
        for j ∈ 1:N_IDX_CHECK
            k = rand(1:len)
            @test (ismissing(p[k]) && ismissing(v[k])) || (p[k] == v[k])
        end
        # AbstractVector{<:Integer} indices
        for j ∈ 1:N_IDX_CHECK
            idx = rand(1:len, rand(1:MAX_IDX_LEN))
            sp = p[idx]
            sv = v[idx]
            @test sp ≅ sv
        end
        @test p[:] ≅ v
    end
end

@testset "indexing_NullablePrimitive_construct" begin
    for i ∈ 1:N_OUTER
        dtype = rand(PRIMITIVE_ELTYPES)
        len = rand(1:MAX_VECTOR_LENGTH)
        v = randmissings(dtype, len)
        p = NullablePrimitive(v)
        # integer indices
        for j ∈ 1:N_IDX_CHECK
            k = rand(1:len)
            @test p[k] ≅ v[k]
        end
        for j ∈ 1:N_IDX_CHECK
            idx = rand(1:len, rand(1:MAX_IDX_LEN))
            @test p[idx] ≅ v[idx]
        end
        @test p[:] ≅ v
    end
end

@testset "indexing_NullablePrimitive_setindex" begin
    for i ∈ 1:N_OUTER
        dtype = rand(PRIMITIVE_ELTYPES)
        len = rand(1:MAX_VECTOR_LENGTH)
        v = randmissings(dtype, len)
        p = NullablePrimitive(v)
        for j ∈ 1:N_IDX_CHECK
            k = rand(1:len)
            x = randmissing(dtype)
            p[k] = x
            @test p[k] ≅ x
        end
        for j ∈ 1:N_IDX_CHECK
            idx = unique(rand(1:len, rand(1:MAX_IDX_LEN)))
            x = randmissings(dtype, length(idx))
            p[idx] = x
            @test p[idx] ≅ x
        end
        # contiguous
        for j ∈ 1:N_IDX_CHECK
            a, b = extrema(rand(1:len, 2))
            idx = a:b
            x = rand(dtype, length(idx))
            p[idx] = x
            @test p[idx] ≅ x
        end
    end
end

@testset "indexing_NullablePrimitive_empty" begin
    for i ∈ 1:N_OUTER
        dtype = rand(PRIMITIVE_ELTYPES)
        p = NullablePrimitive(dtype[])
        @test isempty(p)
        @test_throws BoundsError p[1]
    end
end

@testset "indexing_List_buffer" begin
    len  = 5
    offstype = rand(OFFSET_ELTYPES)
    offs = convert(Vector{offstype}, [0,7,15,18,24,30])
    vals = convert(Vector{UInt8}, codeunits("fireα walk∀ with🐱 me🐟"))
    valspad = zeros(UInt8, Arrow.paddinglength(length(vals)))
    lpad = randpad()
    rpad = randpad()
    b = vcat(lpad, convert(Vector{UInt8}, reinterpret(UInt8, offs)), vals, valspad, rpad)
    l = List{String,offstype}(b, 1+length(lpad), 1+length(lpad)+sizeof(offstype)*length(offs),
                              len, UInt8, length(vals))
    @test offsets(l)[:] == offs
    @test values(l)[:] == vals
    @test l[1] == "fireα "
    @test l[2] == "walk∀ "
    @test l[3] == "wit"
    @test l[4] == "h🐱 "
    @test l[5] == "me🐟"
    @test l[[1,3,5]] == ["fireα ", "wit", "me🐟"]
    @test l[[false,true,false,true,false]] == ["walk∀ ", "h🐱 "]
    @test l[:] == ["fireα ", "walk∀ ", "wit", "h🐱 ", "me🐟"]
end

@testset "indexing_List_construct" begin
    for i ∈ 1:N_OUTER
        len = rand(1:MAX_VECTOR_LENGTH)
        v = String[randstring(rand(0:MAX_STRING_LENGTH)) for i ∈ 1:len]
        l = List(v)
        for j ∈ 1:N_IDX_CHECK
            k = rand(1:len)
            @test l[k] == v[k]
        end
        for j ∈ 1:N_IDX_CHECK
            idx = rand(1:len, rand(1:MAX_IDX_LEN))
            @test l[idx] == v[idx]
        end
        @test l[:] == v
    end
end

@testset "indexing_List_empty" begin
    for i ∈ 1:N_OUTER
        l = List(fill("", i))
        @test l[rand(1:i)] == ""
    end
end

@testset "indexing_NullableList_buffer" begin
    len = 7
    offstype = rand(OFFSET_ELTYPES)
    offs = convert(Vector{offstype}, [0,4,9,9,14,14,17,21])
    offs = convert(Vector{offstype}, [0,9,17,17,28,28,31,35])
    vals = convert(Vector{UInt8}, codeunits("kirk 🚀η spockbones, 💀ncc1701"))
    valspad = zeros(UInt8, Arrow.paddinglength(length(vals)))
    pres = Bool[true,true,false,true,false,true,true]
    mask = Arrow.bitpackpadded(pres)
    lpad = randpad()
    rpad = randpad()
    b = vcat(lpad, mask, convert(Vector{UInt8}, reinterpret(UInt8, offs)), vals, valspad, rpad)
    l = NullableList{String,offstype}(b, 1+length(lpad), 1+length(lpad)+length(mask),
                                      1+length(lpad)+length(mask)+sizeof(offstype)*length(offs),
                                      len, UInt8, length(vals))
    @test offsets(l)[:] == offs
    @test values(l)[:] == vals
    @test l[1] == "kirk 🚀"
    @test l[2] == "η spock"
    @test ismissing(l[3])
    @test l[4] == "bones, 💀"
    @test ismissing(l[5])
    @test l[6] == "ncc"
    @test l[7] == "1701"
    @test l[[2,5,4]] ≅ ["η spock", missing, "bones, 💀"]
    @test l[[true,false,true,false,false,false,false]] ≅ ["kirk 🚀", missing]
    @test l[:] ≅ ["kirk 🚀", "η spock", missing, "bones, 💀", missing, "ncc", "1701"]
end

@testset "indexing_NullableList_construct" begin
    for i ∈ 1:N_OUTER
        len = rand(1:MAX_VECTOR_LENGTH)
        v = Union{String,Missing}[rand(Bool) ? missing : randstring(rand(0:MAX_STRING_LENGTH))
                                  for i ∈ 1:(len-1)]
        # won't work if all strings are empty or missing
        push!(v, "a")
        l = NullableList(v)
        for j ∈ 1:N_IDX_CHECK
            k = rand(1:len)
            @test l[k] ≅ v[k]
        end
        for j ∈ 1:N_IDX_CHECK
            idx = rand(1:len, rand(1:MAX_IDX_LEN))
            @test l[idx] ≅ v[idx]
        end
        @test l[:] ≅ v
    end
end

@testset "indexing_NullableList_empty" begin
    for i ∈ 1:N_OUTER
        l = NullableList(fill("", i))
        @test l[rand(1:i)] == ""
    end
end

@testset "indexing_BitPrimitive_buffer" begin
    len = 10
    bits = vcat(UInt8[0xfd,0x02], zeros(UInt8, 6))
    lpad = randpad()
    rpad = randpad()
    b = vcat(lpad, bits, rpad)
    v = Bool[true,false,true,true,true,true,true,true,false,true]
    p = BitPrimitive(b, 1+length(lpad), 10)
    for j ∈ 1:N_IDX_CHECK
        k = rand(1:len)
        @test p[k] == v[k]
    end
    for j ∈ 1:N_IDX_CHECK
        idx = rand(1:len, rand(1:MAX_IDX_LEN))
        @test p[idx] == v[idx]
    end
    @test p[:] == v
end

@testset "indexing_BitPrimitive_construct" begin
    for i ∈ 1:N_OUTER
        len = rand(1:MAX_VECTOR_LENGTH)
        v = rand(Bool, len)
        p = BitPrimitive(v)
        for j ∈ 1:N_IDX_CHECK
            k = rand(1:len)
            @test p[k] == v[k]
        end
        for j ∈ 1:N_IDX_CHECK
            idx = rand(1:len, rand(1:MAX_IDX_LEN))
            @test p[idx] == v[idx]
        end
        @test p[:] == v
    end
end

@testset "indexing_NullableBitPrimitive_buffer" begin
    len = 7
    bits = vcat(UInt8[0x4f], zeros(UInt8, 7))
    bmask = vcat(UInt8[0x79], zeros(UInt8, 7))
    lpad = randpad()
    rpad = randpad()
    b = vcat(lpad, bmask, bits, rpad)
    v = [true,missing,missing,true,false,false,true]
    p = NullableBitPrimitive(b, 1+length(lpad), 1+length(lpad)+length(bmask), length(v))
    for j ∈ 1:N_IDX_CHECK
        k = rand(1:len)
        @test p[k] ≅ v[k]
    end
    for j ∈ 1:N_IDX_CHECK
        idx = rand(1:len, rand(1:MAX_IDX_LEN))
        @test p[idx] ≅ v[idx]
    end
    @test p[:] ≅ v
end

@testset "indexing_NullableBitPrimitive_construct" begin
    for i ∈ 1:N_OUTER
        len = rand(1:MAX_VECTOR_LENGTH)
        pres = rand(Bool, len)
        vals = rand(Bool, len)
        v = Union{Bool,Missing}[pres[i] ? missing : vals[i] for i ∈ 1:len]
        p = NullableBitPrimitive(v)
        for j ∈ 1:N_IDX_CHECK
            k = rand(1:len)
            @test p[k] ≅ v[k]
        end
        for j ∈ 1:N_IDX_CHECK
            idx = rand(1:len, rand(1:MAX_IDX_LEN))
            @test p[idx] ≅ v[idx]
        end
        @test p[:] ≅ v
    end
end

@testset "indexing_DictEncoding_buffer" begin
    len = 7
    refs = Primitive(Int32[0,1,2,1,0,3,2])
    data = List(["fire", "walk", "with", "me"])
    d = DictEncoding(refs, data)
    @test references(d) == refs
    @test levels(d) == data
    @test d[1] == "fire"
    @test d[2] == "walk"
    @test d[3] == "with"
    @test d[4] == "walk"
    @test d[5] == "fire"
    @test d[6] == "me"
    @test d[7] == "with"
    @test d[[1,4,3,6]] == ["fire", "walk", "with", "me"]
    @test d[[true,false,false,false,false,false,true]] == ["fire", "with"]
    @test d[:] == ["fire", "walk", "with", "walk", "fire", "me", "with"]
end

@testset "indexing_DictEncoding_construct" begin
    v = [-999, missing, 55, -999, 42]
    d = DictEncoding(categorical(v))
    pool = CategoricalPool{Int,Int32}([-999, 55, 42])
    ref = CategoricalArray{Union{Int,Missing},1}(Int32[1,0,2,1,3], pool)
    @test typeof(d.refs) == NullablePrimitive{Int32}
    @test typeof(d.pool) == Primitive{Int64}
    @test d[1] == -999
    @test ismissing(d[2])
    @test d[3] == 55
    @test d[4] == -999
    @test d[5] == 42
    @test d[[1,3,5]] ≅ v[[1,3,5]]
    @test d[[false,true,false,true,false]] ≅ v[[false,true,false,true,false]]
    @test d[1:end] ≅ v
    @test categorical(d) ≅ ref
end


text:
using Printf

function readNums(filePath)
    lines = filter(l -> length(l) > 1, split(open(f->read(f, String), filePath), "\n"))
    nums = map(ln -> [parse(Int8, char) for char in split(ln, "")], lines)
    return hcat(nums...)'
end

function solution1(nums)
    mask = trues(size(nums))
    # left to right
    for column = 1:(size(nums, 2)-1)
        mask[:,column] = mask[:,column] .& (nums[:,column] .< nums[:,column+1])
    end
    # right to left
    for column = 2:size(nums, 2)
        mask[:,column] = mask[:,column] .& (nums[:,column] .< nums[:,column-1])
    end
    # top to bottom
    for row = 1:(size(nums, 1)-1)
        mask[row,:] = mask[row,:] .& (nums[row,:] .< nums[row+1,:])
    end
    # bottom to top
    for row = 2:size(nums, 1)
        mask[row,:] = mask[row,:] .& (nums[row,:] .< nums[row-1,:])
    end
    return sum(nums[mask] .+ 1)
end

function solution2(nums)
    borders = nums .== 9
    function bfs(p_row, p_col)
        # out of bounds
        if (p_row < 1 || p_row > size(borders, 1) || p_col < 1 || p_col > size(borders, 2))
            return 0
        end
        # already visited or border
        if borders[p_row, p_col]
            return 0
        end
        borders[p_row, p_col] = 1
        return (1
                + bfs(p_row + 1, p_col    )
                + bfs(p_row - 1, p_col    )
                + bfs(p_row    , p_col + 1)
                + bfs(p_row    , p_col - 1))
    end
    basin_counts = zeros(Int64, 0)
    for col = 1:(size(borders, 2))
        for row = 1:(size(borders, 1))
            res = bfs(row, col)
            if res != 0
                append!(basin_counts, res)
            end
        end
    end
    largest3 = sort(basin_counts)[end-2:end]
    return largest3[1] * largest3[2] * largest3[3]
end

if length(ARGS) < 1
	println("Provide path of a file you'd like to solve")
	exit(1)
end

nums = readNums(ARGS[1])
println("Day 9:");
Printf.@printf "Solution 1: %d\n" solution1(nums)
Printf.@printf "Solution 2: %d\n" solution2(nums)


text:
abstract type StorageMode end
export MtStorageMode, Shared, Managed, Private
export ReadUsage, WriteUsage, ReadWriteUsage

# Metal Has 4 storage types
# Shared  -> Buffer in Host memory, accessed by the GPU. Requires no sync
# Managed -> Mirrored memory buffers in host and GPU. Requires syncing
# Private -> Memory in Device, not accessible by Host.
# Memoryless -> iOS stuff. ignore it
module AS
import ..StorageMode

struct Shared      <: StorageMode end
struct Managed     <: StorageMode end
struct Private     <: StorageMode end
struct Memoryless  <: StorageMode end
end

const CPUStorage = Union{AS.Shared,AS.Managed}
Base.convert(::Type{MtStorageMode}, ::Type{AS.Shared})     = MtStorageModeShared
Base.convert(::Type{MtStorageMode}, ::Type{AS.Managed})    = MtStorageModeManaged
Base.convert(::Type{MtStorageMode}, ::Type{AS.Private})    = MtStorageModePrivate
Base.convert(::Type{MtStorageMode}, ::Type{AS.Memoryless}) = MtStorageModeMemoryless

Base.convert(::Type{MtResourceOptions}, ::Type{AS.Shared})     = MtResourceStorageModeShared
Base.convert(::Type{MtResourceOptions}, ::Type{AS.Managed})    = MtResourceStorageModeManaged
Base.convert(::Type{MtResourceOptions}, ::Type{AS.Private})    = MtResourceStorageModePrivate
Base.convert(::Type{MtResourceOptions}, ::Type{AS.Memoryless}) = MtResourceStorageModeMemoryless

# Broken because multiple zeros
#@enum_without_prefix MtResourceOptions MtResource
const MtlResourceOptions    = MtResourceOptions

const Shared                = MtResourceStorageModeShared
const Managed               = MtResourceStorageModeManaged
const Private               = MtResourceStorageModePrivate
const Memoryless            = MtResourceStorageModeMemoryless

const DefaultCPUCache       = MtResourceCPUCacheModeDefaultCache
const CombinedWriteCPUCache = MtResourceCPUCacheModeWriteCombined

const DefaultTracking       = MtResourceHazardTrackingModeDefault
const Untracked             = MtResourceHazardTrackingModeUntracked
const Tracked               = MtResourceHazardTrackingModeTracked

const Default               = DefaultCPUCache

Base.:(==)(a::MtlResourceOptions, b::MtStorageMode) =
    (UInt(a) >> 4) == UInt(b)

Base.:(==)(a::MtStorageMode, b::MtlResourceOptions) =
    b == a

##################
Base.convert(::Type{MtResourceUsage}, val::Integer)    = MtResourceUsage(val)

const ReadUsage = MtResourceUsageRead
const WriteUsage = MtResourceUsageWrite
const ReadWriteUsage = convert(MtResourceUsage, MtResourceUsageRead | MtResourceUsageWrite)


text:
using Distributions

xObs     = [18 132; 45 178]
rowSums  = [sum(xObs[i,:]) for i in 1:2]
colSums  = [sum(xObs[:,i]) for i in 1:2]
n        = sum(xObs)

rowProps = rowSums/n
colProps = colSums/n

xExpect  = [colProps[c]*rowProps[r]*n for r in 1:2, c in 1:2]
testStat = sum([(xObs[r,c]-xExpect[r,c])^2 / xExpect[r,c] for r in 1:2,c in 1:2])
pVal = ccdf(Chisq(1),testStat)

println("Chi-squared value: ", testStat)
println("P-value: ", pVal)

max_stars_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook

max_issues_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook

max_forks_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook