text:
## Filter hash construct core

function fink_filter_hash(c, L; nx=256, wd=2, t=1, shift=false, Omega=false, safety_on=true, wd_cutoff=1)
    # -------- compute the filter bank
    filt, hash = fink_filter_bank(c, L; nx=nx, wd=wd, t=t, shift=shift, Omega=Omega, safety_on=safety_on, wd_cutoff=wd_cutoff)

    # -------- list of non-zero pixels
    flist = fink_filter_list(filt)

    # -------- pack everything you need into the info structure
    hash["filt_index"] = flist[1]
    hash["filt_value"] = flist[2]

    # -------- compute matrix that projects iso coeffs, add to hash
    S1_iso_mat = S1_iso_matrix(hash)
    hash["S1_iso_mat"] = S1_iso_mat
    S2_iso_mat = S2_iso_matrix(hash)
    hash["S2_iso_mat"] = S2_iso_mat
    hash["num_iso_coeff"] = size(S1_iso_mat)[1] + size(S2_iso_mat)[1] + 2
    hash["num_coeff"] = size(S1_iso_mat)[2] + size(S2_iso_mat)[2] + 2

    return hash
end

function fink_filter_hash_p(c, L; nx=256, wd=2, t=1, shift=false, Omega=false, safety_on=true, wd_cutoff=1, p=2)
    # -------- compute the filter bank
    filt, hash = fink_filter_bank_p(c, L; nx=nx, wd=wd, t=t, shift=shift, Omega=Omega, safety_on=safety_on, wd_cutoff=wd_cutoff, p=p)

    # -------- list of non-zero pixels
    flist = fink_filter_list(filt)

    # -------- pack everything you need into the info structure
    hash["filt_index"] = flist[1]
    hash["filt_value"] = flist[2]

    # -------- compute matrix that projects iso coeffs, add to hash
    S1_iso_mat = S1_iso_matrix(hash)
    hash["S1_iso_mat"] = S1_iso_mat
    S2_iso_mat = S2_iso_matrix(hash)
    hash["S2_iso_mat"] = S2_iso_mat
    hash["num_iso_coeff"] = size(S1_iso_mat)[1] + size(S2_iso_mat)[1] + 2
    hash["num_coeff"] = size(S1_iso_mat)[2] + size(S2_iso_mat)[2] + 2

    return hash
end

function fink_filter_hash_gpu(c, L; nx=256, wd=2, t=1, shift=false, Omega=false, safety_on=true,threeD=false,cz=1,nz=256)

    # -------- compute the filter bank
    filt, hash = fink_filter_bank(c, L; nx=nx, wd=wd, t=t, shift=shift, Omega=Omega, safety_on=safety_on)

    flist = fink_filter_list(filt)
    # -------- pack everything you need into the info structure
    hash["filt_index"] = flist[1]
    hash["filt_value"] = flist[2]

    if threeD
        hash=fink_filter_bank_3dizer(hash,cz,nz=nz)
    end

    hash=list_to_box(hash,dim=threeD ? 3 : 2)

    # send to gpu
    CUDA.reclaim()
    hash["filt_vals"] = [CUDA.CuArray(val) for val in hash["filt_vals"]]
    CUDA.reclaim()

    #Not implemented
    # -------- compute matrix that projects iso coeffs, add to hash
    #S1_iso_mat = S1_iso_matrix(hash)
    hash["S1_iso_mat"] = nothing
    #S2_iso_mat = S2_iso_matrix(hash)
    hash["S2_iso_mat"] = nothing

    return hash
end

function fink_filter_hash_gpu_sparse(c, L; nx=256, wd=2, t=1, shift=false, Omega=false, safety_on=true, wd_cutoff=1, prec=Float32)
    filt, hash = eqws.fink_filter_bank(c, L; nx=nx, wd=wd, t=t, shift=shift, Omega=Omega, safety_on=safety_on, wd_cutoff=wd_cutoff)

    Nf = size(filt)[3]
    hash["Nf"] = Nf

    gpu_filt = CUDA.zeros(Float64,size(filt))
    #send to GPU by convert to sparse CPU, transfer to sparse GPU, convert to dense GPU
    for i=1:Nf
        gpu_filt[:,:,i] .= CUDA.CuArray(CuSparseMatrixCSC(sparse(filt[:,:,i])))
    end
    hash["gpu_filt"] = gpu_filt

    S1_iso_mat = CuSparseMatrixCSC(S1_iso_matrix_gpu(hash, prec))
    hash["S1_iso_mat"] = S1_iso_mat
    S2_iso_mat = CuSparseMatrixCSC(S2_iso_matrix_gpu(hash, prec))
    hash["S2_iso_mat"] = S2_iso_mat
    hash["num_iso_coeff"] = size(S1_iso_mat)[1] + size(S2_iso_mat)[1] + 2
    hash["num_coeff"] = size(S1_iso_mat)[2] + size(S2_iso_mat)[2] + 2

    return hash
end

function fink_filter_hash_gpu_fake(c, L; nx=256, wd=2, t=1, shift=false, Omega=false, safety_on=true,threeD=false,cz=1,nz=256)

    # -------- compute the filter bank
    filt, hash = fink_filter_bank(c, L; nx=nx, wd=wd, t=t, shift=shift, Omega=Omega, safety_on=safety_on)

    flist = fink_filter_list(filt)
    # -------- pack everything you need into the info structure
    hash["filt_index"] = flist[1]
    hash["filt_value"] = flist[2]

    if threeD
        hash=fink_filter_bank_3dizer(hash,cz,nz=nz)
    end

    hash=list_to_box(hash,dim=threeD ? 3 : 2)

    #Not implemented
    # -------- compute matrix that projects iso coeffs, add to hash
    #S1_iso_mat = S1_iso_matrix(hash)
    hash["S1_iso_mat"] = nothing
    #S2_iso_mat = S2_iso_matrix(hash)
    hash["S2_iso_mat"] = nothing

    return hash
end

function fink_filter_bank_3dizer(hash, cz; nz=256, tz=1, Omega3d=false)
    @assert hash["t"]==2 #can remove later when we add support for tz=2
    @assert tz==1

    # -------- set parameters
    nx = convert(Int64,hash["npix"])
    n2d = size(hash["filt_value"])[1]
    dz = nz/2-1
    J = size(hash["j_value"])[1]
    L = size(hash["theta_value"])[1]
    psi_ind_L2d = hash["psi_ind_L"]

    im_scale = convert(Int8,log2(nz))
    # -------- number of bins in radial direction (size scales)
    K = (im_scale-3)*cz + 1
    normj = 1/sqrt(cz)

    Omega2d = haskey(hash, "Omega_index")

    if Omega2d
        tot_num_filt = n2d*K + 3
    else
        tot_num_filt = n2d*K + 1
    end

    # -------- allocate output array of zeros
    J_L_K     = zeros(Int32, tot_num_filt,3)
    psi_ind_L = zeros(Int32, tot_num_filt)
    psi_index = zeros(Int32, J+1, L, K+1)
    k_value   = zeros(Float64, K)
    logr      = zeros(Float64,nz)
    filt_temp = zeros(Float64,nx,nx,nz)
    F_z       = zeros(Float64,nx,nx,nz)
    #F_p       = zeros(Float64,nx,nx,nz)

    filtind = fill(CartesianIndex{3}[], tot_num_filt)
    filtval = fill(Float64[], tot_num_filt)

    for z = 1:convert(Int64,dz+1)
        sz = mod(z+dz,nz)-dz-1
        z2 = sz^2
        logr[z] = 0.5*log2(max(1,z2))
    end

    @inbounds for k_ind = 1:K
        k = (k_ind-1)/cz
        k_value[k_ind] = k  # store for later
        krad  = im_scale-k-2
        Δk    = abs.(logr.-krad)
        kmask = (Δk .<= 1)
        k_count = count(kmask)
        k_vals = findall(kmask)

        # -------- radial part
        #I have somehow dropped a factor of sqrt(2) which has been reinserted as a 0.5... fix my brain
        @views F_z[:,:,kmask].= reshape(cos.(Δk[kmask] .* (π/2)),1,1,k_count)
        @inbounds for index = 1:n2d
            p_ind = hash["filt_index"][index]
            p_filt = hash["filt_value"][index]
            #F_p[p_ind,:] .= p_filt
            f_ind    = k_ind + (index-1)*K
            @views filt_tmp = p_filt.*F_z[p_ind,kmask]
            #temp_ind = findall((F_p .> 1E-13) .& (F_z .> 1E-13))
            #@views filt_tmp = F_p[temp_ind].*F_z[temp_ind]
            ind = findall(filt_tmp .> 1E-13)
            @views filtind[f_ind] = map(x->CartesianIndex(p_ind[x[1]][1],p_ind[x[1]][2],k_vals[x[2]]),ind)
            #filtind[f_ind] = temp_ind[ind]
            @views filtval[f_ind] = filt_tmp[ind]
            @views J_L_K[f_ind,:] = [hash["J_L"][index,1],hash["J_L"][index,2],k_ind-1]
            psi_index[hash["J_L"][index,1]+1,hash["J_L"][index,2]+1,k_ind] = f_ind
            psi_ind_L[f_ind] = psi_ind_L2d[index]
        end
    end

    filter_power = zeros(nx,nx,nz)
    center_power = zeros(nx,nx,nz)
    for j = 1:J
        for l = 1:L
            for k = 1:K
                index = psi_index[j,l,k]
                filter_power[filtind[index]] .+= filtval[index].^2
            end
        end
    end

    for k = 1:K
        index = psi_index[J+1,1,k]
        filter_power[filtind[index]] .+= filtval[index].^2
    end

    #### ------- This is phi0 containing the plane
    # phi contains power near k=0 not yet accounted for
    # for plane half-covered (tz=1), add other half-plane

    if tz == 1
        filter_power .+= circshift(filter_power[:,:,end:-1:1],(0,0,1))
    end

    center_power = fftshift(filter_power)
    phi_cen = zeros(nx, nx, nz)
    phi_cen[:,:,nz÷2:nz÷2+2] = center_power[:,:,nz÷2+3].*ones(1,1,3)
    phi_cen_shift = fftshift(phi_cen)
    ind = findall(phi_cen_shift .> 1E-13)
    val = phi_cen_shift[ind]

    # -------- before adding ϕ to filter bank, renormalize ψ if t=1
    if tz==1 filtval[1:n2d*K] .*= sqrt(2.0) end  # double power for half coverage

    # -------- add result to filter array
    filtind[n2d*K + 1] = ind
    filtval[n2d*K + 1] = sqrt.(val)
    J_L_K[n2d*K + 1,:] = [J+1,0,K+1]
    psi_index[J+1,1,K+1] = n2d*K + 1
    psi_ind_L[n2d*K + 1] = 0

    if Omega3d
        #### ------- This is omega0 containing the plane

        omega0 = zeros(nx, nx, nz)
        omega0[:,:,nz÷2:nz÷2+2] .= ones(nx,nx,3)
        omega0 .-= phi_cen

        omega0_shift = fftshift(omega0)
        ind = findall(omega0_shift .> 1E-13)
        val = omega0_shift[ind]

        # -------- add result to filter array
        filtind[n2d*K + 2] = ind
        filtval[n2d*K + 2] = sqrt.(val)
        J_L_K[n2d*K + 2,:] = [J+1,1,K+1]
        psi_index[J+1,2,K+1] = n2d*K + 2
        psi_ind_L[n2d*K + 2] = 0

        #### ------- This is omega3 around the edges

        filter_power = zeros(nx,nx,nz)
        center_power = zeros(nx,nx,nz)
        for index = 1:n2d*K + 2
            filter_power[filtind[index]] .+= filtval[index].^2
        end

        if tz == 1
            filter_power .+= circshift(filter_power[:,:,end:-1:1],(0,0,1))
        end

        center_power = ones(nx,nx,nz) .- filter_power
        ind = findall(center_power .> 1E-13)
        val = center_power[ind]

        # -------- add result to filter array
        filtind[n2d*K + 3] = ind
        filtval[n2d*K + 3] = sqrt.(val)
        J_L_K[n2d*K + 3,:] = [J+1,2,K+1]
        psi_index[J+1,3,K+1] = n2d*K + 3
        psi_ind_L[n2d*K + 3] = 0
    end

    # -------- metadata dictionary
    info=Dict()
    info["npix"]         = hash["npix"]
    info["j_value"]      = hash["j_value"]
    info["theta_value"]  = hash["theta_value"]
    info["2d_psi_index"]    = hash["psi_index"]
    info["2d_phi_index"]    = hash["phi_index"]
    info["2d_J_L"]          = hash["J_L"]
    info["2d_t"]           = hash["t"]
    info["2d_wd"]           = hash["wd"]
    info["2d_fs_center_r"]  = hash["fs_center_r"]
    info["2d_filt_index"]   = hash["filt_index"]
    info["2d_filt_value"]   = hash["filt_value"]
    if Omega2d
        info["Omega_index"]   = hash["Omega_index"]
    end

    info["nz"]              = nz
    info["cz"]              = cz
    info["J_L_K"]           = psi_index
    info["psi_index"]       = psi_index
    info["k_value"]         = k_value
    info["filt_index"]      = filtind
    info["filt_value"]      = filtval
    info["tz"]             = tz
    info["Omega3d"]        = Omega3d
    info["psi_ind_L"]       = psi_ind_L

    # -------- compute matrix that projects iso coeffs, add to hash
    S1_iso_mat = S1_iso_matrix3d(info)
    info["S1_iso_mat"] = S1_iso_mat
    S2_iso_mat = S2_iso_matrix3d(info)
    info["S2_iso_mat"] = S2_iso_mat
    #

    return info
end

## Filter bank utilities

function fink_filter_bank(c, L; nx=256, wd=2, t=1, shift=false, Omega=false, safety_on=true, wd_cutoff=1)
    #c     - sets the scale sampling rate (1 is dyadic, 2 is half dyadic)
    #L     - number of angular bins (usually 8*t or 16*t)
    #wd    - width of the wavelets (default 1, wd=2 for a double covering)
    #t    - plane coverage (default 1, full 2pi 2)
    #shift - shift in θ by 1/2 of the θ spacing
    #Omega - true= append Omega filter (all power beyond Nyquist) so the sum of filters is 1.0

    # -------- assertion errors to make sure arguments are reasonable
    #@test wd <= L/2

    # -------- set parameters
    dθ   = t*π/L
    θ_sh = shift ? dθ/2 : 0.0
    dx   = nx/2-1

    im_scale = convert(Int8,log2(nx))
    # -------- number of bins in radial direction (size scales)
    J = (im_scale-3)*c + 1
    normj = 1/sqrt(c)

    # -------- allocate output array of zeros
    filt      = zeros(nx, nx, J*L+(Omega ? 2 : 1))
    psi_index = zeros(Int32, J, L)
    psi_ind_in= zeros(Int32, J*L+(Omega ? 2 : 1), 2)
    psi_ind_L = zeros(Int32, J*L+(Omega ? 2 : 1))
    theta     = zeros(Float64, L)
    j_value   = zeros(Float64, J)
    info=Dict{String,Any}()

    # -------- compute the required wd
    j_rad_exp = zeros(J)
    for j_ind = 1:J
        j = (j_ind-1)/c
        jrad  = im_scale-j-2
        j_rad_exp[j_ind] = 2^(jrad)
    end

    wd_j = max.(ceil.(wd_cutoff.*L./(t.*π.*j_rad_exp)),wd)

    if !safety_on
        wd_j.=wd
    end

    # loop over wd from small to large
    ## there is some uneeded redundancy added by doing this esp in l loop
    for wd in sort(unique(wd_j))
        # -------- allocate theta and logr arrays
        θ    = zeros(nx, nx)
        logr = zeros(nx, nx)

        wdθ  = wd*dθ
        norm = 1.0/(sqrt(wd))
        # -------- loop over l
        for l = 0:L-1
            θ_l        = dθ*l+θ_sh
            theta[l+1] = θ_l

        # -------- allocate anggood BitArray
            anggood = falses(nx, nx)

        # -------- loop over pixels
            for x = 1:nx
                sx = mod(x+dx,nx)-dx -1    # define sx,sy so that no fftshift() needed
                for y = 1:nx
                    sy = mod(y+dx,nx)-dx -1
                    θ_pix  = mod(atan(sy, sx)+π -θ_l, 2*π)
                    θ_good = abs(θ_pix-π) <= wdθ

                    # If this is a pixel we might use, calculate log2(r)
                    if θ_good
                        anggood[y, x] = θ_good
                        θ[y, x]       = θ_pix
                        r2            = sx^2 + sy^2
                        logr[y, x]    = 0.5*log2(max(1,r2))
                    end
                end
            end
            angmask = findall(anggood)
        # -------- compute the wavelet in the Fourier domain
        # -------- the angular factor is the same for all j
            F_angular = norm .* cos.((θ[angmask].-π).*(L/(2*wd*t)))

        # -------- loop over j for the radial part
        #    for (j_ind, j) in enumerate(1/c:1/c:im_scale-2)
            j_ind_w_wd = findall(wd_j.==wd)
            for j_ind in j_ind_w_wd
                j = (j_ind-1)/c
                j_value[j_ind] = 1+j  # store for later
                jrad  = im_scale-j-2
                Δj    = abs.(logr[angmask].-jrad)
                rmask = (Δj .<= 1) #deprecating the 1/c to 1, constant width

        # -------- radial part
                F_radial = normj .* cos.(Δj[rmask] .* (π/2)) #deprecating c*π/2 to π/2
                ind      = angmask[rmask]
        #      Let's have these be (J,L) if you reshape...
        #        f_ind    = (j_ind-1)*L+l+1
                f_ind    = j_ind + l*J
                filt[ind, f_ind] = F_radial .* F_angular[rmask]
                psi_index[j_ind,l+1] = f_ind
                psi_ind_in[f_ind,:] = [j_ind-1,l]
                psi_ind_L[f_ind] = 1
            end
        end
    end

    # -------- phi contains power near k=0 not yet accounted for
    filter_power = (sum(filt.*filt, dims=3))[:,:,1]

    # -------- for plane half-covered (t=1), add other half-plane
    if t == 1
        filter_power .+= circshift(filter_power[end:-1:1,end:-1:1],(1,1))
    end

    # -------- compute power required to sum to 1.0
    i0 = round(Int16,nx/2-2)
    i1 = round(Int16,nx/2+4)
    center_power = 1.0 .- fftshift(filter_power)[i0:i1,i0:i1]
    zind = findall(center_power .< 1E-15)
    center_power[zind] .= 0.0  # set small numbers to zero
    phi_cen = zeros(nx, nx)
    phi_cen[i0:i1,i0:i1] = sqrt.(center_power)

    # -------- before adding ϕ to filter bank, renormalize ψ if t=1
    if t==1 filt .*= sqrt(2.0) end  # double power for half coverage

    # -------- add result to filter array
    phi_index  = J*L+1
    filt[:,:,phi_index] .= fftshift(phi_cen)
    psi_ind_in[phi_index,:] = [J,0]
    psi_ind_L[phi_index] = 0

    if Omega     # append a filter containing the rest (outside Nyquist)
        filter_power += filt[:,:,phi_index].^2
        edge_power    = 1.0 .- filter_power
        zind          = findall(edge_power .< 1E-15)
        edge_power[zind]     .= 0.0  # set small numbers to zero
        Omega_index           = J*L+2
        info["Omega_index"]   = Omega_index
        filt[:,:,Omega_index] = sqrt.(edge_power)
        psi_ind_in[Omega_index,:] = [J,1]
        psi_ind_L[Omega_index] = 0
    end

    # -------- metadata dictionary
    info["npix"]         = nx
    info["j_value"]      = j_value
    info["theta_value"]  = theta
    info["psi_index"]    = psi_index
    info["phi_index"]    = phi_index
    info["J_L"]          = psi_ind_in
    info["t"]            = t
    info["wd"]           = wd_j
    info["wd_cutoff"]    = wd_cutoff
    info["fs_center_r"]  = j_rad_exp
    info["psi_ind_L"]    = psi_ind_L
    info["p"]            = 2

    return filt, info
end

function fink_filter_bank_p(c, L; nx=256, wd=2, t=1, shift=false, Omega=false, safety_on=true, wd_cutoff=1, p=2)
    #c     - sets the scale sampling rate (1 is dyadic, 2 is half dyadic)
    #L     - number of angular bins (usually 8*t or 16*t)
    #wd    - width of the wavelets (default 1, wd=2 for a double covering)
    #t    - plane coverage (default 1, full 2pi 2)
    #shift - shift in θ by 1/2 of the θ spacing
    #Omega - true= append Omega filter (all power beyond Nyquist) so the sum of filters is 1.0

    # -------- assertion errors to make sure arguments are reasonable
    #@test wd <= L/2

    # -------- set parameters
    dθ   = t*π/L
    θ_sh = shift ? dθ/2 : 0.0
    dx   = nx/2-1

    im_scale = convert(Int8,log2(nx))
    # -------- number of bins in radial direction (size scales)
    J = (im_scale-3)*c + 1
    normj = 1/sqrt(c)

    # -------- allocate output array of zeros
    filt      = zeros(nx, nx, J*L+(Omega ? 2 : 1))
    psi_index = zeros(Int32, J, L)
    psi_ind_in= zeros(Int32, J*L+(Omega ? 2 : 1), 2)
    psi_ind_L = zeros(Int32, J*L+(Omega ? 2 : 1))
    theta     = zeros(Float64, L)
    j_value   = zeros(Float64, J)
    info=Dict{String,Any}()

    # -------- compute the required wd
    j_rad_exp = zeros(J)
    for j_ind = 1:J
        j = (j_ind-1)/c
        jrad  = im_scale-j-2
        j_rad_exp[j_ind] = 2^(jrad)
    end

    wd_j = max.(ceil.(wd_cutoff.*L./(t.*π.*j_rad_exp)),wd)

    if !safety_on
        wd_j.=wd
    end

    # loop over wd from small to large
    ## there is some uneeded redundancy added by doing this esp in l loop
    for wd in sort(unique(wd_j))
        # -------- allocate theta and logr arrays
        θ    = zeros(nx, nx)
        logr = zeros(nx, nx)

        wdθ  = wd*dθ
        norm = 1.0/((wd).^(1/p))
        # -------- loop over l
        for l = 0:L-1
            θ_l        = dθ*l+θ_sh
            theta[l+1] = θ_l

        # -------- allocate anggood BitArray
            anggood = falses(nx, nx)

        # -------- loop over pixels
            for x = 1:nx
                sx = mod(x+dx,nx)-dx -1    # define sx,sy so that no fftshift() needed
                for y = 1:nx
                    sy = mod(y+dx,nx)-dx -1
                    θ_pix  = mod(atan(sy, sx)+π -θ_l, 2*π)
                    θ_good = abs(θ_pix-π) <= wdθ

                    # If this is a pixel we might use, calculate log2(r)
                    if θ_good
                        anggood[y, x] = θ_good
                        θ[y, x]       = θ_pix
                        r2            = sx^2 + sy^2
                        logr[y, x]    = 0.5*log2(max(1,r2))
                    end
                end
            end
            angmask = findall(anggood)
        # -------- compute the wavelet in the Fourier domain
        # -------- the angular factor is the same for all j
            F_angular = norm .* cos.((θ[angmask].-π).*(L/(2*wd*t))).^(2/p)

        # -------- loop over j for the radial part
        #    for (j_ind, j) in enumerate(1/c:1/c:im_scale-2)
            j_ind_w_wd = findall(wd_j.==wd)
            for j_ind in j_ind_w_wd
                j = (j_ind-1)/c
                j_value[j_ind] = 1+j  # store for later
                jrad  = im_scale-j-2
                Δj    = abs.(logr[angmask].-jrad)
                rmask = (Δj .<= 1) #deprecating the 1/c to 1, constant width

        # -------- radial part
                F_radial = normj .* cos.(Δj[rmask] .* (π/2)).^(2/p) #deprecating c*π/2 to π/2
                ind      = angmask[rmask]
        #      Let's have these be (J,L) if you reshape...
        #        f_ind    = (j_ind-1)*L+l+1
                f_ind    = j_ind + l*J
                filt[ind, f_ind] = F_radial .* F_angular[rmask]
                psi_index[j_ind,l+1] = f_ind
                psi_ind_in[f_ind,:] = [j_ind-1,l]
                psi_ind_L[f_ind] = 1
            end
        end
    end

    # -------- phi contains power near k=0 not yet accounted for
    filter_power = (sum(filt.^p, dims=3))[:,:,1]

    # -------- for plane half-covered (t=1), add other half-plane
    if t == 1
        filter_power .+= circshift(filter_power[end:-1:1,end:-1:1],(1,1))
    end

    # -------- compute power required to sum to 1.0
    i0 = round(Int16,nx/2-2)
    i1 = round(Int16,nx/2+4)
    center_power = 1.0 .- fftshift(filter_power)[i0:i1,i0:i1]
    zind = findall(center_power .< 1E-15)
    center_power[zind] .= 0.0  # set small numbers to zero
    phi_cen = zeros(nx, nx)
    phi_cen[i0:i1,i0:i1] = (center_power).^(1/p)

    # -------- before adding ϕ to filter bank, renormalize ψ if t=1
    if t==1 filt .*= 2.0^(1/p) end  # double power for half coverage

    # -------- add result to filter array
    phi_index  = J*L+1
    filt[:,:,phi_index] .= fftshift(phi_cen)
    psi_ind_in[phi_index,:] = [J,0]
    psi_ind_L[phi_index] = 0

    if Omega     # append a filter containing the rest (outside Nyquist)
        filter_power += filt[:,:,phi_index].^p
        edge_power    = 1.0 .- filter_power
        zind          = findall(edge_power .< 1E-15)
        edge_power[zind]     .= 0.0  # set small numbers to zero
        Omega_index           = J*L+2
        info["Omega_index"]   = Omega_index
        filt[:,:,Omega_index] = (edge_power).^(1/p)
        psi_ind_in[Omega_index,:] = [J,1]
        psi_ind_L[Omega_index] = 0
    end

    # -------- metadata dictionary
    info["npix"]         = nx
    info["j_value"]      = j_value
    info["theta_value"]  = theta
    info["psi_index"]    = psi_index
    info["phi_index"]    = phi_index
    info["J_L"]          = psi_ind_in
    info["t"]            = t
    info["wd"]           = wd_j
    info["wd_cutoff"]    = wd_cutoff
    info["fs_center_r"]  = j_rad_exp
    info["psi_ind_L"]    = psi_ind_L
    info["p"]            = p

    return filt, info
end

## Make a list of non-zero pixels for the (Fourier plane) filters

function fink_filter_list(filt)
    (ny,nx,Nf) = size(filt)

    # Allocate output arrays
    filtind = fill(CartesianIndex{2}[], Nf)
    filtval = fill(Float64[], Nf)

    # Loop over filters and record non-zero values
    for l=1:Nf
        f = @view filt[:,:,l]
        ind = findall(f .> 1E-13)
        val = f[ind]
        filtind[l] = ind
        filtval[l] = val
    end
    return [filtind, filtval]
end

function list_to_box(hash;dim=2)#this modifies the hash
    (Nf,) = size(hash["filt_index"])
    nx=hash["npix"]
    if dim==3 nz=hash["nz"] end

    # Allocate output arrays
    filtmms = Array{Int64}(undef,dim,2,Nf)
    filtvals = Array{Any}(undef,Nf)#is this a terrible way to allocate memory?
    if dim==2
        filt=zeros(Float64,nx,nx)
    else
        filt=zeros(Float64,nx,nx,nz)
    end
    # Loop over filters and record non-zero values
    for l=1:Nf
        filt.=0.
        f_i = hash["filt_index"][l]
        v_i=hash["filt_value"][l]
        for i = 1:length(f_i)
            filt[f_i[i]]= v_i[i]
        end
        filt = FFTW.fftshift(filt)

        ind = findall(filt .> 1E-13)
        if dim==2
            ind=getindex.(ind,[1 2])
        else
            ind=getindex.(ind,[1 2 3])
        end
        mins=minimum(ind,dims=1)
        maxs=maximum(ind,dims=1)

        filtmms[:,1,l] = mins
        filtmms[:,2,l] = maxs
        if dim==2
            filtvals[l]=filt[mins[1]:maxs[1],mins[2]:maxs[2]]
        else
            filtvals[l]=filt[mins[1]:maxs[1],mins[2]:maxs[2],mins[3]:maxs[3]]
        end
    end

    hash["filt_index"] = nothing
    hash["filt_value"] = nothing

    hash["filt_mms"] = filtmms
    hash["filt_vals"]=filtvals
    return hash
end

max_stars_repo_name: andrew-saydjari/DHC

max_issues_repo_name: andrew-saydjari/DHC

max_forks_repo_name: andrew-saydjari/DHC


text:
export ishom, ishet, heterozygosity, het

#TODO how to treat haploids?
"""
```
ishom(locus::T) where T <: GenoArray
ishom(locus::Genotype)
ishom(locus::Missing)
```
A series of methods to test if a locus or loci are homozygous and return `true` if
it is, `false` if it isn't, and `missing` if it's `missing`. The vector version
simply maps the function over the elements.
"""
@inline function ishom(locus::Genotype)
    # if the first equals all others, return true
    return all(@inbounds first(locus) .== locus)
end

ishom(locus::Missing) = missing

@inline function ishom(locus::T) where T<:GenoArray
    return @inbounds map(ishom, locus)
end

@inline function ishom(locus::T) where T<:Base.SkipMissing
    return @inbounds map(ishom, locus)
end

"""
    ishom(locus::Genotype, allele::Signed)
    ishom(loci::GenoArray, allele::Signed)
Returns `true` if the `locus`/`loci` is/are homozygous for the specified `allele`.
"""
function ishom(geno::T, allele::U) where T<:Genotype where U<:Integer
    ∈(allele, geno) & ishom(geno) ? true : false
end

ishom(geno::T, allele::U) where T<:GenoArray where U<:Integer = map(i -> ishom(i, allele), geno)

ishom(geno::Missing, allele::U) where U<:Integer = missing

"""
```
ishet(locus::T) where T <: GenoArray
ishet(locus::Genotype)
ishet(locus::Missing)
```
A series of methods to test if a locus or loci are heterozygous and return `true` if
it is, `false` if it isn't. The vector version simply broadcasts the function over the
elements. Under the hood, this function is simply `!ishom`.
"""
@inline function ishet(locus::Genotype)
    return !ishom(locus)
end

ishet(locus::Missing) = missing

@inline function ishet(locus::T) where T<:GenoArray
    return @inbounds map(ishet, locus)
end

@inline function ishet(locus::T) where T<:Base.SkipMissing
    return @inbounds map(ishet, locus)
end

"""
    ishet(locus::Genotype, allele::Signed)
    ishet(loci::GenoArray, allele::Signed)
Returns `true` if the `locus`/`loci` is/are heterozygous for the specified `allele`.
"""
function ishet(geno::T, allele::U) where T<:Genotype where U<:Integer
    ∈(allele, geno) & !ishom(geno) ? true : false
end

ishet(geno::T, allele::U) where T<:GenoArray where U<:Integer = map(i -> ishet(i, allele), geno)

ishet(geno::Missing, allele::U) where U<:Integer = missing

"""
    counthet(geno::T, allele::Int) where T<:GenoArray
    counthet(geno::T, allele::AbstractVector{U}) where T<:GenoArray where U<:Integer
Given a `GenoArray`, count the number of times `allele` appears in the
heterozygous state.
"""
function counthet(geno::T, allele::U) where T<:GenoArray where U<:Integer
    mapreduce(i -> ishet(i, allele), +, skipmissing(geno))
end

function counthet(geno::T, allele::AbstractVector{U}) where T<:GenoArray where U<:Integer
    tmp = skipmissing(geno)
    isempty(tmp) && return fill(0, length(allele))
    map(_allele -> mapreduce(i -> ishet(i, _allele), +, tmp), allele)
end

"""
    counthom(geno::T, allele::Int) where T<:GenoArray
Given a `GenoArray`, count the number of times `allele` appears in the
homozygous state.
"""
function counthom(geno::T, allele::U) where T<:GenoArray where U <: Signed
    mapreduce(i -> ishom(i, allele), +, skipmissing(geno))
end

function counthom(geno::T, allele::AbstractVector{U}) where T<:GenoArray where U<:Integer
    tmp = skipmissing(geno)
    isempty(tmp) && return fill(0, length(allele))
    map(_allele -> mapreduce(i -> ishom(i, _allele), +, tmp), allele)
end

"""
    gene_diversity_nei87(het_exp::Union{Missing,AbstractFloat}, het_obs::Union{Missing,AbstractFloat}, n::Union{Integer, Float64}, corr::Bool = true)
Calculate overall gene diversity with the adjustment/correction given by Nei:

Hₜ = 1 −sum(pbar²ᵢ + Hₛ/(ñ * np) − Het_obs/(2ñ*np))
- _ñ_ is the number of genotypes for a locus for a population
- _np_ is the number of genotypes of a locus across all populations
    - i.e. sum(_ñ_)
- _pbar²_ is the observed homozygosity of a locus for that population
- _Hₛ_ is the within population gene diversity given by:
    - Hₛ =  ñ/(ñ-1) * (1 - sum(pbar²ᵢ - Het_observed / 2ñ))

(Nei M. (1987) Molecular Evolutionary Genetics. Columbia University Press).
use `corr = false` to ignore sample-size correction `* n/(n-1)`.
"""
@inline function gene_diversity_nei87(het_exp::T, het_obs::T, n::Union{Integer,T}; corr::Bool = true) where T<: AbstractFloat
    corr_val = corr ? n/(n-1.0) : 1.0
    return @fastmath (het_exp - (het_obs/n/2.0)) * corr_val
end

@inline function gene_diversity_nei87(het_exp::AbstractFloat, het_obs::Missing, n::Union{Integer,AbstractFloat}; corr::Bool = true)
    return missing
end

@inline function gene_diversity_nei87(het_exp::Missing, het_obs::AbstractFloat, n::Union{Integer,AbstractFloat}; corr::Bool = true)
    return missing
end

@inline function gene_diversity_nei87(het_exp::Missing, het_obs::Missing, n::Union{Integer,AbstractFloat}; corr::Bool = true)
    return missing
end

"""
    hetero_o(data::T) where T <: GenoArray
Returns observed heterozygosity as a mean of the number of heterozygous genotypes, defined
as genotypes returning `true` for `ishet()`. This is numerically feasible because
`true` values are mathematically represented as `1`, whereas `false` are represented
as `0`.
"""
@inline function hetero_o(data::T) where T <: GenoArray
    adjusted_vector = ishet(data) |> skipmissing
    isempty(adjusted_vector) ? missing : mean(adjusted_vector)
end

"""
    hetero_e(allele_freqs::Vector{T}) where T <: GenoArray
Returns the expected heterozygosity of an array of genotypes,
calculated as 1 - sum of the squared allele frequencies.
"""
@inline function hetero_e(data::T) where T <: GenoArray
    all(ismissing.(data)) == true ? missing : 1.0 - mapreduce(i -> i^2, + , allele_freq_vec(data))
end

"""
    heterozygosity(data::PopData, by::String = "locus")
Calculate observed and expected heterozygosity in a `PopData` object. For loci,
heterozygosity is calculated in the Nei fashion, such that heterozygosity is
calculated as the average over heterozygosity per locus per population.
### Modes
- `"locus"` or `"loci"` : heterozygosity per locus (default)
- `"sample"` or `"ind"` or `"individual"` : heterozygosity per individual/sample
- `"population"` or `"pop"` : heterozygosity per population
## Example
heterozygosity(@nancycats, "population" )
"""
function heterozygosity(data::PopData, by::String = "locus")
    if by ∈ ["locus", "loci"]
        tmp = DataFrames.combine(
                groupby(data.loci, [:locus, :population]),
                :genotype => nonmissing => :n_tmp,
                :genotype => hetero_o => :het_pop_obs,
                :genotype => hetero_e => :het_pop_exp
            )
        return DataFrames.combine(
                groupby(tmp, :locus),
                :n_tmp => sum => :n,
                :het_pop_obs => (h_o -> mean(skipmissing(h_o))) => :het_obs,
                :het_pop_exp => (h_e -> mean(skipmissing(h_e))) => :het_exp
            )

    elseif lowercase(by) ∈  ["sample", "ind", "individual"]
        return DataFrames.combine(
                groupby(data.loci, :name),
                :genotype => nonmissing => :n,
                :genotype => hetero_o => :het_obs
            )

    elseif lowercase(by) ∈  ["pop", "population"]
        return DataFrames.combine(
                groupby(data.loci, :population),
                :genotype => nonmissing => :n,
                :genotype => hetero_o => :het_obs,
                :genotype => hetero_e => :het_exp
            )
    else
        error("please specify by = \"locus\", \"sample\", or \"population\"")
    end
end

const het = heterozygosity

#NOTE this is not intended to be performant. It's a convenience function.
"""
    het_sample(data::PopData, individual::String)
Calculate the observed heterozygosity for an individual in a `PopData` object.
"""
@inline function het_sample(data::PopData, individual::String)
    data.loci[data.loci.name .== individual, :genotype] |> hetero_o
end

ext: jl


text:
## Nucleotide IUPAC Regex ##

macro ntr_str(pattern, flags...)
  npat = UInt8[]
  flag_open = false
  for i in 1:length(pattern)
    elem = pattern[i]
    if elem=='[' && !flag_open && i!=1 && pattern[i-1]!='\\'
      flag_open = true
      push!(npat,elem)
    elseif elem=='[' && !flag_open && i==1
      flag_open = true
      push!(npat,elem)
    elseif flag_open && elem==']' && i!=1 && pattern[i-1]!='\\'
      flag_open = false
      push!(npat,elem)
    elseif haskey(_AMBIGUOUS_NUCLEIC_IUPAC,nucleotide(elem))
      if !flag_open push!(npat,'[') end
      append!(npat,_AMBIGUOUS_NUCLEIC_IUPAC[nucleotide(elem)])
      if !flag_open push!(npat,']') end
    else
      push!(npat,elem)
    end
  end
  Regex(ASCIIString(npat),flags...)
end

## Amino Acid IUPAC Regex ##

macro aar_str(pattern, flags...)
  npat = UInt8[]
  flag_open = false
  for i in 1:length(pattern)
    elem = pattern[i]
    if elem=='[' && !flag_open && i!=1 && pattern[i-1]!='\\'
      flag_open = true
      push!(npat,elem)
    elseif elem=='[' && !flag_open && i==1
      flag_open = true
      push!(npat,elem)
    elseif flag_open && elem==']' && i!=1 && pattern[i-1]!='\\'
      flag_open = false
      push!(npat,elem)
    elseif haskey(_AMBIGUOUS_AMINO_IUPAC,aminoacid(elem))
      if !flag_open push!(npat,'[') end
      append!(npat,_AMBIGUOUS_AMINO_IUPAC[aminoacid(elem)])
      if !flag_open push!(npat,']') end
    else
      push!(npat,elem)
    end
  end
  Regex(ASCIIString(npat),flags...)
end

## DNA to RNA ##

function dna2rna!(s::AbstractArray{Nucleotide})
  for i in 1:length(s)
    if Int(s[i]) | 32 == 116
      s[i] += 1
    end
  end
  s
end

dna2rna(s::AbstractArray{Nucleotide}) = dna2rna!(copy(s))

## RNA to DNA ##

function rna2dna!(s::AbstractArray{Nucleotide})
  for i in 1:length(s)
    if Int(s[i]) | 32 == 117
      s[i] -= 1
    end
  end
  s
end

rna2dna(s::AbstractArray{Nucleotide}) = rna2dna!(copy(s))

## Complement ##

complement!(s::AbstractArray{Nucleotide},d) = swap!(s, d)
complement(s::AbstractArray{Nucleotide},d) = swap(s,d)

## Reverse Complement ##

function reversecomplement!{T<:BioUnit}(s::AbstractArray{T},y::Array{T})
  i=1
  j=length(s)
  while i<=j
    I=s[i]
    J=s[j]
    s[i] = y[J]
    s[j] = y[I]
    i += 1
    j -= 1
  end
  s
end

reversecomplement{T<:BioUnit}(s::AbstractArray{T},y::Array{T}) = reversecomplement!(copy(s),y)

## percent of GC ##

function percentGC(seq::ASCIIString)
  len = length(seq)
  sum = 0
  for elem in seq
    if elem | 32 == 'c' || elem | 32 == 'g' || elem | 32 == 's'
      sum += 1
    end
  end
  sum/len
end

function percentGC(seq::Vector{Nucleotide})
  len = length(seq)
  sum = 0
  for elem in seq
    if elem | 32 == 'c' || elem | 32 == 'g' || elem | 32 == 's'
      sum += 1
    end
  end
  sum/len
end


text:
module GHWT_tf_2d

include("utils.jl")

using ..GraphSignal, ..GraphPartition, ..BasisSpecification, ..GHWT,  ..GHWT_2d, LinearAlgebra

include("common.jl")

include("partition_fiedler.jl")

export ghwt_tf_init_2d, ghwt_tf_bestbasis_2d, ghwt_tf_synthesis_2d, ghwt_tf_init_2d_Lindberg, eghwt_init_2d, eghwt_bestbasis_2d, eghwt_synthesis_2d, eghwt_init_2d_Lindberg, BS2loc

"""
    tf2d_init(tag::Matrix{Int},tag_r::Matrix{Int})

Storing the relation between the coefficient location in the 2D coefficients matrix of all levels with the (level, tag, region) information

### Input Arguments
* `tag`: The matrix of the tag information
* `tag_r`: The matrix of the region information

### Output Arguments
* `tag2ind`: tag2ind[(j,k,l)] is the location of the column of coefficients on level j, region k and tag l
* `ind2tag`: key and item of dicionary `tag2ind` interchanged

Copyright 2018 The Regents of the University of California

Implemented by Yiqun Shao (Adviser: Dr. Naoki Saito)
"""
function tf2d_init(tag::Matrix{Int},tag_r::Matrix{Int})
    (m,n) = size(tag)
    tag2ind = Dict{Tuple{Int,Int,Int},Int}()
    ind2tag = Dict{Int,Tuple{Int,Int,Int}}()
    for j = 1:n
        for i = 1:m
            tag2ind[(j,tag_r[i,j],tag[i,j])] = i + (j-1)*m
            ind2tag[i+(j-1)*m] = (j,tag_r[i,j],tag[i,j])
        end
    end
    return tag2ind, ind2tag
end

"""
    dmatrix_new,tag2ind_new,ind2tag_new,tag_tf = tf_core_2d_col(dmatrix::Matrix{Float64},tag2ind::Dict{Tuple{Int,Int,Int},Int},ind2tag::Dict{Int,Tuple{Int,Int,Int}},jmax::Integer)

One forward iteration of 2-d time-frequency adapted ghwt on the column direction

### Input Arguments
* `dmatrix`: The cost-functional values of coeffcients
* `tag2ind`: tag2ind[(j,k,l)] is the location of the column of coefficients on level j, region k and tag l
* `ind2tag`: key and item of dicionary `tag2ind` interchanged

### Output Arguments
* `dmatrix_new`: The cost-functional values of coeffcients
* `tag2ind_new`: tag2ind_new[(j,k,l)] is the location of the column of coefficients on level j, region k and tag l
* `ind2tag_new`: key and item of dicionary `tag2ind` interchanged
* `tag_tf`: recording the time-or-frequency information

Copyright 2018 The Regents of the University of California

Implemented by Yiqun Shao (Adviser: Dr. Naoki Saito)
"""
function tf_core_2d_col(dmatrix::Matrix{Float64},tag2ind::Dict{Tuple{Int,Int,Int},Int},ind2tag::Dict{Int,Tuple{Int,Int,Int}},jmax::Integer)
    # power of the cost-functional
    costp = 1
    p,q = size(dmatrix)
    dmatrix_new = zeros(p,q)

    # Initialization
    tag_tf = zeros(Bool,(p,q))
    tag2ind_new = Dict{Tuple{Int,Int,Int},Int}()
    ind2tag_new = Dict{Int,Tuple{Int,Int,Int}}()

    s = 1
    # Iterate through columns
    for i = 1:q
        # Find out the corresponding level,region,tag information
        j,k,l = ind2tag[i]
        # only look at the entry with even tag l to avoid duplication
        if l%2 == 0 && j < jmax
            # search for the (j,k,l+1) entry.
            # the (j,k,l) and (j,k,l+1) span the same subspace as (j+1,2*k,l/2) and (j+1,2*k+1,l/2)
            # (j,k,l) and (j,k,l+1) are `frequency-direction`
            # (j,k,l) and (j,k,l+1) are `time-direction`
            if haskey(tag2ind,(j,k,l+1)) # check for degenerate case ((j,k,l+1) doesn't exist)
                freqcos = abs.(dmatrix[:,i]).^costp + abs.(dmatrix[:,tag2ind[(j,k,l+1)]]).^costp
            else
                freqcos = abs.(dmatrix[:,i]).^costp
            end
            if ~haskey(tag2ind,(j+1,2*k,l/2)) # check for degenerate case ((j+1,2*k,l/2) or (j+1,2*k+1,l/2) doesn't exist)
                timecos = abs.(dmatrix[:,tag2ind[(j+1,2*k+1,l/2)]]).^costp
            elseif ~haskey(tag2ind,(j+1,2*k+1,l/2))
                timecos = abs.(dmatrix[:,tag2ind[(j+1,2*k,l/2)]]).^costp
            else
                timecos = abs.(dmatrix[:,tag2ind[(j+1,2*k+1,l/2)]]).^costp + abs.(dmatrix[:,tag2ind[(j+1,2*k,l/2)]]).^costp
            end
            tag_tf[:,s] = timecos.>=freqcos
            dmatrix_new[:,s] = tag_tf[:,s].*freqcos + .~tag_tf[:,s].*timecos
            ind2tag_new[s] = (j,k,l/2)
            tag2ind_new[(j,k,l/2)] = s
            s = s + 1
        end
    end
    dmatrix_new = dmatrix_new[:,1:s-1]
    tag_tf = tag_tf[:,1:s-1]
    return dmatrix_new,tag2ind_new,ind2tag_new,Matrix{UInt8}(tag_tf)
end

"""
    dmatrix_new,tag2ind_new,ind2tag_new,tag_tf = tf_core_2d_row(dmatrix::Matrix{Float64},tag2ind::Dict{Tuple{Int,Int,Int},Int},ind2tag::Dict{Int,Tuple{Int,Int,Int}},jmax::Integer)

Almost the same as tf_core_2d_col. But all the operations on matrix are row-wise.
Due to the column-based matrix feature of Julia language. We are currently only using the tf_core_2d_col here.

### Input Arguments
* `dmatrix`: The cost-functional values of coeffcients
* `tag2ind`: tag2ind[(j,k,l)] is the location of the column of coefficients on level j, region k and tag l
* `ind2tag`: key and item of dicionary `tag2ind` interchanged

### Output Arguments
* `dmatrix_new`: The cost-functional values of coeffcients
* `tag2ind_new`: tag2ind_new[(j,k,l)] is the location of the row of coefficients on level j, region k and tag l
* `ind2tag_new`: key and item of dicionary `tag2ind` interchanged
* `tag_tf`: recording the time-or-frequency information

Copyright 2018 The Regents of the University of California

Implemented by Yiqun Shao (Adviser: Dr. Naoki Saito)
"""
function tf_core_2d_row(dmatrix::Matrix{Float64},tag2ind::Dict{Tuple{Int,Int,Int},Int},ind2tag::Dict{Int,Tuple{Int,Int,Int}},jmax::Integer)
    costp = 1
    p,q = size(dmatrix)
    dmatrix_new = zeros(p,q)
    tag_tf = zeros(Bool,(p,q))
    tag2ind_new = Dict{Tuple{Int,Int,Int},Int}()
    ind2tag_new = Dict{Int,Tuple{Int,Int,Int}}()

    s = 1
    for i = 1:p
        j,k,l = ind2tag[i]
        if l%2 == 0 && j < jmax
            if haskey(tag2ind,(j,k,l+1))
                freqcos = abs.(dmatrix[i,:]).^costp + abs.(dmatrix[tag2ind[(j,k,l+1)],:]).^costp
            else
                freqcos = abs.(dmatrix[i,:]).^costp
            end
            if ~haskey(tag2ind,(j+1,2*k,l/2))
                timecos = abs.(dmatrix[tag2ind[(j+1,2*k+1,l/2)],:]).^costp
            elseif ~haskey(tag2ind,(j+1,2*k+1,l/2))
                timecos = abs.(dmatrix[tag2ind[(j+1,2*k,l/2)],:]).^costp
            else
                timecos = abs.(dmatrix[tag2ind[(j+1,2*k+1,l/2)],:]).^costp + abs.(dmatrix[tag2ind[(j+1,2*k,l/2)],:]).^costp
            end
            tag_tf[s,:] = timecos.>=freqcos
            dmatrix_new[s,:] = tag_tf[s,:].*freqcos + .~tag_tf[s,:].*timecos
            ind2tag_new[s] = (j,k,l/2)
            tag2ind_new[(j,k,l/2)] = s
            s = s + 1
        end
    end
    dmatrix_new = dmatrix_new[1:s-1,:]
    tag_tf = tag_tf[1:s-1,:]
    return dmatrix_new,tag2ind_new,ind2tag_new,Matrix{UInt8}(tag_tf)
end

"""
    Bbasis, infovec = ghwt_tf_bestbasis_2d(dmatrix::Matrix{Float64},GProws::GraphPart,GPcols::GraphPart)

Implementation of 2d time-frequency adapted GHWT method.
Modified from the idea in paper "Image compression with adaptive Haar-Walsh tilings" by Maj Lindberg and Lars F. Villemoes.

### Input Arguments
* `dmatrix`: The ghwt expanding coefficients of matrix on all levels contatenated in a 2d matrix.
* `GProws`: Corresponding to the affinity matrix on rows.
* `GPcols`: Corresponding to the affinity matrix on cols.

### Output Arguments
* `Bbasis`: The coefficients of the best-basis.
* `infovec`: [infovec[i,1],infovec[i,2]] is the location of Bbasis[i] in dmatrix.

Copyright 2018 The Regents of the University of California

Implemented by Yiqun Shao (Adviser: Dr. Naoki Saito)
"""
function ghwt_tf_bestbasis_2d(matrix::Matrix{Float64},GProws::GraphPart,GPcols::GraphPart)
    dmatrix = ghwt_analysis_2d(matrix, GProws, GPcols)
    # Initialization
    tag_r_row = rs_to_region(GProws.rs, GProws.tag)
    tag_r_col = rs_to_region(GPcols.rs, GPcols.tag)
    tag2ind_row,ind2tag_row = tf2d_init(Matrix{Int}(GProws.tag),Matrix{Int}(tag_r_row))
    tag2ind_col,ind2tag_col = tf2d_init(Matrix{Int}(GPcols.tag),Matrix{Int}(tag_r_col))

    jmax_col = size(GPcols.tag,2)
    jmax_row = size(GProws.tag,2)

    # (level,region,tag) information on the columns
    TAG2IND_col = Vector(undef, jmax_col)
    TAG2IND_col[1] = tag2ind_col
    IND2TAG_col = Vector(undef, jmax_col)
    IND2TAG_col[1] = ind2tag_col

    # (level,region,tag) information on the rows
    TAG2IND_row = Vector(undef, jmax_row)
    TAG2IND_row[1] = tag2ind_row
    IND2TAG_row = Vector(undef, jmax_row)
    IND2TAG_row[1] = ind2tag_row

    # recording if time or frequency direction is chosen
    TAG_tf = Matrix(undef, jmax_row,jmax_col)

    # storing the cost functional values
    DMATRIX = Matrix(undef, jmax_row, jmax_col)

    DMATRIX[1,1] = dmatrix

    # comparison in the row direction
    for i=1:(jmax_row - 1)
        DMATRIX[i+1,1], TAG2IND_row[i+1], IND2TAG_row[i+1], TAG_tf[i+1,1] = tf_core_2d_col(Array{Float64,2}(DMATRIX[i,1]'), TAG2IND_row[i], IND2TAG_row[i],jmax_row+1-i)
        TAG_tf[i+1,1] = Array{UInt8,2}((TAG_tf[i+1,1] .+ UInt8(2))')
        DMATRIX[i+1,1] = Array{Float64,2}(DMATRIX[i+1,1]')
    end

    # comparison in the column direction
    for i=1:(jmax_col - 1)
        DMATRIX[1,i+1], TAG2IND_col[i+1], IND2TAG_col[i+1], TAG_tf[1,i+1] = tf_core_2d_col(DMATRIX[1,i],TAG2IND_col[i], IND2TAG_col[i], jmax_col+1-i)
    end

    #DMATRIX[1,1] = 0 # release the memory
    # comparison in both directions
    for i = 2:jmax_row
        for j = 2:jmax_col
            (temp_col,_,_,temp_tf_col) = tf_core_2d_col(DMATRIX[i,j-1], TAG2IND_col[j-1], IND2TAG_col[j-1],jmax_col + 2 - j)
            (temp_row,_,_,temp_tf_row) = tf_core_2d_col(Array{Float64,2}(DMATRIX[i-1,j]'), TAG2IND_row[i-1], IND2TAG_row[i-1],jmax_row + 2 - i)
            temp_row = Array{Float64,2}(temp_row')
            temp_tf_row = Array{UInt8,2}((temp_tf_row .+ UInt8(2))')

            row_or_col = temp_col.<=temp_row
            DMATRIX[i,j] = row_or_col.*temp_col + .~row_or_col.*temp_row
            TAG_tf[i,j] = Matrix{UInt8}(row_or_col).*temp_tf_col + Matrix{UInt8}(.~row_or_col).*temp_tf_row

            #DMATRIX[i-1,j] = 0
        end
        #DMATRIX[i,1] = 0 # release the memory
    end

    #for j = 2: jmax_col
        #DMATRIX[jmax_row, j] = 0 # release the memory
    #end
    #recover the bestbasis
    infovec = Array{Int}([jmax_row; jmax_col; 1; 1])
    for iter = 1:(jmax_row + jmax_col - 2)
        newinfovec = -1*ones(Int,(4,2*size(infovec,2)))
        for h = 1:size(infovec,2)
            m = infovec[1,h]
            n = infovec[2,h]
            p = infovec[3,h]
            q = infovec[4,h]

            #if the bestbasis comes from column direction
            if TAG_tf[m,n][p,q] == 0 #time domain
                j,k,l = IND2TAG_col[n][q]
                if haskey(TAG2IND_col[n-1],(j+1,2*k,l))
                    q1 = TAG2IND_col[n-1][(j+1,2*k,l)]
                    newinfovec[:,2*h-1] = [m;n-1;p;q1]
                end
                if haskey(TAG2IND_col[n-1],(j+1,2*k+1,l))
                    q2 = TAG2IND_col[n-1][(j+1,2*k+1,l)]
                    newinfovec[:,2*h] = [m;n-1;p;q2]
                end
            end

            if TAG_tf[m,n][p,q] == 1 # frequency domain
                j,k,l = IND2TAG_col[n][q]
                if haskey(TAG2IND_col[n-1],(j,k,2*l))
                    q1 = TAG2IND_col[n-1][(j,k,2*l)]
                    newinfovec[:,2*h-1] = [m;n-1;p;q1]
                end
                if haskey(TAG2IND_col[n-1],(j,k,2*l+1))
                    q2 = TAG2IND_col[n-1][(j,k,2*l+1)]
                    newinfovec[:,2*h] = [m;n-1;p;q2]
                end
            end

            #if the bestbasis comes from row direction
            if TAG_tf[m,n][p,q] == 2 # time domain
                j,k,l = IND2TAG_row[m][p]
                if haskey(TAG2IND_row[m-1],(j+1,2*k,l))
                    p1 = TAG2IND_row[m-1][(j+1,2*k,l)]
                    newinfovec[:,2*h-1] = [m-1;n;p1;q]
                end
                if haskey(TAG2IND_row[m-1],(j+1,2*k+1,l))
                    p2 = TAG2IND_row[m-1][(j+1,2*k+1,l)]
                    newinfovec[:,2*h] = [m-1;n;p2;q]
                end
            end

            if TAG_tf[m,n][p,q] == 3 # frequency domain
                j,k,l = IND2TAG_row[m][p]
                if haskey(TAG2IND_row[m-1],(j,k,2*l))
                    p1 = TAG2IND_row[m-1][(j,k,2*l)]
                    newinfovec[:,2*h-1] = [m-1;n;p1;q]
                end
                if haskey(TAG2IND_row[m-1],(j,k,2*l+1))
                    p2 = TAG2IND_row[m-1][(j,k,2*l+1)]
                    newinfovec[:,2*h] = [m-1;n;p2;q]
                end
            end
        end
        infovec = newinfovec[:,newinfovec[1,:].!=-1]
    end
#return Array{Float64,2}(dmatrix[sub2ind(size(dmatrix),infovec[3,:],infovec[4,:])]'),infovec[3:4,:]
return dmatrix[(LinearIndices(size(dmatrix)))[CartesianIndex.(infovec[3,:],infovec[4,:])]],infovec[3:4,:]
end

"""
    ghwt_tf_synthesis_2d_core!(dmatrix::Array{Float64,3},tag::Array{Int,2},rs::Array{Int,2})

Core function of ghwt_tf_synthesis_2d. Synthesize on column direction

### Input Arguments
* `dmatrix`: Same size as dmatrix, but only selected basis vectors are nonzero with expanding coeffcients.
* `tag`: Tag information in GP.
* `rs`: Rs information in GP.

### Output Arguments
* `dmatrix`: Synthesized on column direction.

Copyright 2018 The Regents of the University of California

Implemented by Yiqun Shao (Adviser: Dr. Naoki Saito)
"""
function ghwt_tf_synthesis_2d_core!(dmatrix::Array{Float64,3},tag::Array{Int,2},rs::Array{Int,2})
    jmax = size(tag,2)
    for j = 1:jmax-1
        regioncount = count(!iszero, rs[:,j]) - 1
        for r = 1:regioncount
            # the index that marks the start of the first subregion
            rs1 = rs[r,j]

            # the index that is one after the end of the second subregion
            rs3 = rs[r+1,j]

            #the number of points in the current region
            n = rs3 - rs1

            # proceed if not all coefficients are zero in this region
            if count(!iszero, dmatrix[rs1:rs3-1,j,:]) > 0

                if n == 1
                    ### SCALING COEFFICIENT (n == 1)
                    dmatrix[rs1,j+1,:] = dmatrix[rs1,j,:] + dmatrix[rs1,j+1,:]
                elseif n > 1
                    # the index that marks the start of the second subregion
                    rs2 = rs1+1;
                    while rs2 < rs3 && tag[rs2,j+1] != 0
                        rs2 = rs2+1;
                    end

                    # the parent region is a copy of the subregion
                    if rs2 == rs3
                        dmatrix[rs1:rs3-1,j+1,:] = dmatrix[rs1:rs3-1,j,:] + dmatrix[rs1:rs3-1,j+1,:]

                        # the parent region has 2 child regions
                    else
                        # the number of points in the first subregion
                        n1 = rs2 - rs1

                        # the number of points in the second subregion
                        n2 = rs3 - rs2

                        ### SCALING COEFFICIENTS
                        dmatrix[rs1,j+1,:] = ( sqrt(n1).*dmatrix[rs1,j,:] + sqrt(n2).*dmatrix[rs1+1,j,:] )./sqrt(n) + dmatrix[rs1,j+1,:]
                        dmatrix[rs2,j+1,:] = ( sqrt(n2).*dmatrix[rs1,j,:] - sqrt(n1).*dmatrix[rs1+1,j,:] )./sqrt(n) + dmatrix[rs2,j+1,:]

                        # HAAR-LIKE & WALSH-LIKE COEFFICIENTS

                        # search through the remaining coefficients in each subregion
                        parent = rs1+2
                        child1 = rs1+1
                        child2 = rs2+1
                        while child1 < rs2 || child2 < rs3
                            # subregion 1 has the smaller tag
                            if child2 == rs3 || (tag[child1,j+1] < tag[child2,j+1] && child1 < rs2)
                                dmatrix[child1,j+1,:] = dmatrix[parent,j,:] + dmatrix[child1,j+1,:]
                                child1 = child1+1;
                                parent = parent+1;

                                # subregion 2 has the smaller tag
                            elseif child1 == rs2 || (tag[child2,j+1] < tag[child1,j+1] && child2 < rs3)
                                dmatrix[child2,j+1,:] = dmatrix[parent,j,:] + dmatrix[child2,j+1,:]
                                child2 = child2+1;
                                parent = parent+1;

                                # both subregions have the same tag
                            else
                                dmatrix[child1,j+1,:] = ( dmatrix[parent,j,:] + dmatrix[parent+1,j,:] )/sqrt(2) + dmatrix[child1,j+1,:]
                                dmatrix[child2,j+1,:] = ( dmatrix[parent,j,:] - dmatrix[parent+1,j,:] )/sqrt(2) + dmatrix[child2,j+1,:]
                                child1 = child1+1;
                                child2 = child2+1;
                                parent = parent+2;
                            end
                        end
                    end
                end
            end
        end
    end
end

"""
    matrix = ghwt_tf_synthesis_2d(dmatrix::Matrix{Float64}, GProws::GraphPart, GPcols::GraphPart)

Synthesis the matrix from the coefficients of selected basis vectors.

### Input Arguments
* `dmatrix`: Only selected basis vectors are nonzero with expanding coeffcients.
* `GProws`: Corresponding to the affinity matrix on rows.
* `GPcols`: Corresponding to the affinity matrix on cols.

### Output Arguments
* `matrix`: Synthesized matrix

Copyright 2018 The Regents of the University of California

Implemented by Yiqun Shao (Adviser: Dr. Naoki Saito)
"""
function ghwt_tf_synthesis_2d(dvec::Vector{Float64}, infovec::Matrix{Int}, GProws::GraphPart, GPcols::GraphPart)
    tag_col, rs_col = GPcols.tag, GPcols.rs
    tag_row, rs_row = GProws.tag, GProws.rs
    fcols, jmax_col = size(tag_col);
    frows, jmax_row = size(tag_row);

    dmatrix = zeros(frows*jmax_row, fcols*jmax_col)
    for i in 1:size(infovec,2)
        dmatrix[infovec[1,i],infovec[2,i]] = dvec[i]
    end

    matrix_r = Array{Float64,2}(dmatrix');
    matrix_r = reshape(matrix_r, fcols, jmax_col, frows*jmax_row);
    ghwt_tf_synthesis_2d_core!(matrix_r, Matrix{Int}(tag_col), Matrix{Int}(rs_col));
    matrix_r = matrix_r[:,end,:];
    #matrix_r = reshape(matrix_r, fcols, frows*jmax_row);
    matrix_r = reshape(Array{Float64,2}(matrix_r'),frows,jmax_row, fcols);
    ghwt_tf_synthesis_2d_core!(matrix_r, Matrix{Int}(tag_row), Matrix{Int}(rs_row));
    matrix_r = matrix_r[:,end,:];
    #matrix_r = reshape(matrix_r,frows, fcols);
    matrix = zeros(frows, fcols)
    matrix[GProws.ind,GPcols.ind] = matrix_r
    return matrix
end

"""
    dmatrix, GProws, GPcols = ghwt_tf_init_2d(matrix::Matrix{Float64})

    Partition matrix first to get GProws and GPcols. Then expand matrix
    in two directions to get dmatrix

### Input Arguments
* `matrix::Matrix{Float64}`: an input matrix

### Output Argument
* `GProws::GraphPart`: partitioning using rows as samples
* `GPcols::GraphPart`: partitioning using cols as samples
* `dmatrix::matrix{Float64}`: expansion coefficients of matrix
"""
function ghwt_tf_init_2d(matrix::Matrix{Float64})
  #Partition matrix recuisively using Dhillon's method
  GProws, GPcols = PartitionTreeMatrixDhillon(matrix)
  ghwt_core!(GProws)
  ghwt_core!(GPcols)

  matrix_re = matrix[GProws.ind, GPcols.ind]

  (N, jmax_row) = size(GProws.rs)
  N = N-1

  # expand on the row direction
  (frows, fcols) = size(matrix_re)
  dmatrix = zeros((N,jmax_row,fcols))
  dmatrix[:,jmax_row,:] = matrix_re

  ghwt_core!(GProws, dmatrix)

  dmatrix = reshape(dmatrix, (frows*jmax_row, fcols))'

  # expand on the column direction
  (N, jmax_col) = size(GPcols.rs)
  N = N - 1
  dmatrix2 = zeros(N,jmax_col,size(dmatrix,2))
  dmatrix2[:, jmax_col, :] = dmatrix
  ghwt_core!(GPcols, dmatrix2)

  dmatrix = reshape(dmatrix2, (fcols*jmax_col, frows*jmax_row))'

  return Array{Float64,2}(dmatrix), GProws, GPcols
end

function ghwt_tf_init_2d_Lindberg(matrix::Matrix{Float64})
  #Partition matrix recuisively using Dhillon's method
  m,n = size(matrix)
  GProws = regularhaar(m)
  GPcols = regularhaar(n)
  ghwt_core!(GProws)
  ghwt_core!(GPcols)

  matrix_re = matrix[GProws.ind, GPcols.ind]

  (N, jmax_row) = size(GProws.rs)
  N = N-1

  # expand on the row direction
  (frows, fcols) = size(matrix_re)
  dmatrix = zeros((N,jmax_row,fcols))
  dmatrix[:,jmax_row,:] = matrix_re

  ghwt_core!(GProws, dmatrix)

  dmatrix = reshape(dmatrix, (frows*jmax_row, fcols))'

  # expand on the column direction
  (N, jmax_col) = size(GPcols.rs)
  N = N - 1
  dmatrix2 = zeros(N,jmax_col,size(dmatrix,2))
  dmatrix2[:, jmax_col, :] = dmatrix
  ghwt_core!(GPcols, dmatrix2)

  dmatrix = Array{Float64,2}(reshape(dmatrix2, (fcols*jmax_col, frows*jmax_row))')

  return dmatrix, GProws, GPcols
end

function regularhaar(n::Int)
    if log2(n)%1 != 0
        error("Input number should be some power of 2")
    else
        ind = Array{Int64}(1:n)
        l = Int64(log2(n))
        rs = zeros(Int64,n+1,l+1)
        for j = 1:l+1
            gap = Int64(n/(2^(j-1)))
            for i = 1:(2^(j-1) + 1)
                rs[i,j] = Int64((i-1)*gap + 1)
            end
        end
    end
    return GraphPart(ind, rs)
end

"""
    dmatrix = ghwt_tf_init_2d(matrix::Matrix{Float64}, GProws::GraphPart, GPcols::GraphPart)

    Partition matrix first to get GProws and GPcols. Then expand matrix
    in two directions to get dmatrix

### Input Arguments
* `matrix::Matrix{Float64}`: an input matrix
* `GProws::GraphPart`: partitioning using rows as samples
* `GPcols::GraphPart`: partitioning using cols as samples

### Output Argument

* `dmatrix::matrix{Float64}`: expansion coefficients of matrix
"""
function ghwt_tf_init_2d(matrix::Matrix{Float64}, GProws::GraphPart, GPcols::GraphPart)

  ghwt_core!(GProws)
  ghwt_core!(GPcols)
  matrix_re = matrix[GProws.ind, GPcols.ind]

  (N, jmax_row) = size(GProws.rs)
  N = N-1

  # expand on the row direction
  (frows, fcols) = size(matrix_re)
  dmatrix = zeros((N,jmax_row,fcols))
  dmatrix[:,jmax_row,:] = matrix_re

  ghwt_core!(GProws, dmatrix)

  dmatrix = reshape(dmatrix, (frows*jmax_row, fcols))'

  # expand on the column direction
  (N, jmax_col) = size(GPcols.rs)
  N = N - 1
  dmatrix2 = zeros(N,jmax_col,size(dmatrix,2))
  dmatrix2[:, jmax_col, :] = dmatrix
  ghwt_core!(GPcols, dmatrix2)

  dmatrix = Array{Float64,2}(reshape(dmatrix2, (fcols*jmax_col, frows*jmax_row))')

  return dmatrix
end

function BS2ind(GP::GraphPart, BS::BasisSpec)
    tag = GP.tag
    N = size(tag,1)
    jmax = size(tag,2)

    tag_r = rs_to_region(GP.rs, tag)
    tag2ind, ind2tag = tf2d_init(Array{Int,2}(tag),tag_r)

    ind = fill(0, N)
    if BS.c2f
        for i = 1:N
            j2, j = BS.levlist[i]
            k = tag_r[j2,j]
            l = tag[j2,j]
            ind[i] = tag2ind[(j,k,l)]
        end
    else
        tag_r_f2c = Array{Int,2}(fine2coarse!(GP, dmatrix = Array{Float64,3}(reshape(tag_r,(size(tag_r,1),size(tag_r,2),1))), coefp = true)[:,:,1])
        tag_f2c = GP.tagf2c
        for i = 1:N
            j2, j = BS.levlist[i]
            k = tag_r_f2c[j2,j]
            l = tag_f2c[j2,j]
            ind[i] = tag2ind[(jmax + 1 - j, k, l)]
        end
    end
    return ind
end

function BS2loc(dmatrix::Matrix{Float64}, GProws::GraphPart, GPcols::GraphPart, BSrows::BasisSpec, BScols::BasisSpec)
    indrows = BS2ind(GProws, BSrows)
    indcols = BS2ind(GPcols, BScols)

    m = length(indrows)
    n = length(indcols)

    loc = fill(0, (2,m*n))
    dvec = fill(0., m*n)
    global idx = 1
    for i = 1:n
        for j = 1:m
            loc[1,idx] = indrows[j]
            loc[2,idx] = indcols[i]
            dvec[idx] = dmatrix[indrows[j], indcols[i]]
            global idx += 1
        end
    end
    return dvec, loc
end

# A useful and meaningful set of aliases here
const eghwt_init_2d = ghwt_tf_init_2d
const eghwt_bestbasis_2d = ghwt_tf_bestbasis_2d
const eghwt_synthesis_2d = ghwt_tf_synthesis_2d
const eghwt_init_2d_Lindberg = ghwt_tf_init_2d_Lindberg

end # of module GHWT_tf_2d


text:
#!/usr/bin/julia

# Copyright (C) 2021, 2022 Matheus Fernandes Bigolin <mfrdrbigolin@disroot.org>
# SPDX-License-Identifier: MIT

"Day Ten, Syntax Scoring."

#= This day was considerably easy, maybe because I already worked with lexical
 = analysis before (vide Braskag).
 =#

include("Utils.jl")

using Statistics
using .Utils

const OPENER = ['(', '[', '{', '<']
const CLOSER = [')', ']', '}', '>']
const CORRUPT_ADDEND = [3, 57, 1197, 25137]

"""Calculate both the corruption score and the incompletion score of a `system`
and return them, in this order, inside a tuple."""
function solve(system :: Array{String}) :: Tuple{Int64, Int64}
    corrupt_score = 0
    incomplete_scores = Array{Int64}([])

    for line in system
        order = Array{Char}([])
        is_corrupt = false

        for sep in line
            if sep in OPENER
                push!(order, sep)
            elseif sep != CLOSER[findfirst(==(order[end]), OPENER)]
                corrupt_score += CORRUPT_ADDEND[findfirst(==(sep), CLOSER)]
                is_corrupt = true

                break
            else
                pop!(order)
            end
        end

        #= Here I assumed, based on the enunciation of the problem, that a line
         = must be either corrupt or incomplete.
         =#
        if !is_corrupt
            closers = order[end:-1:1]
            addends = [findfirst(==(sep), OPENER) for sep in closers]

            push!(incomplete_scores, reduce((a,b) -> 5*a + b, addends))
        end
    end

    return (corrupt_score, incomplete_scores |> median)
end

function main()
    usage_and_exit(length(ARGS) != 1)

    input_data = readlines(first(ARGS))

    (S₁, S₂) = input_data |> solve

    @assert S₁ == 268845
    @assert S₂ == 4038824534

    S₁ |> println
    S₂ |> println
end

if abspath(PROGRAM_FILE) == @__FILE__
    main()
end


text:
using Measures: Add, Min, Max, Div, Mul, Neg
using LinearAlgebra

# Measure Constants
# -----------------

const cx = Length{:cx}
const cy = Length{:cy}
const sx = Length{:sx}
const sy = Length{:sy}

*(a::T, b::Type{cx}) where T = x_measure(a)
*(a::T, b::Type{cy}) where T = y_measure(a)
*(a::T, b::Type{sx}) where T = size_x_measure(a)
*(a::T, b::Type{sy}) where T = size_y_measure(a)

# Pixels are not typically used in Compose in preference of absolute
# measurements or measurements relative to parent canvases. So for the
# 'px' constant, we just punt and give something do something vaguely
# reasonable.

const assumed_ppmm = 3.78 # equivalent to 96 DPI
const px = mm/assumed_ppmm

const MeasureOrNumber = Union{Measure, Number}
const XYTupleOrVec = Union{Tuple{MeasureOrNumber,MeasureOrNumber}, Vec}

# Scaling w and h components
# --------------------------

# Compute the length of the given type.
sum_component(::Type{T}, l) where T <: Length = 0.0
sum_component(::Type{T}, l::T) where T <: Length = l.value
sum_component(::Type{T}, l::Add) where T <: Length = sum_component(T, l.a) + sum_component(T, l.b)

# Scale a length component by some factor.
scale_component(::Type{T}, scale, l) where T <: Length = l
scale_component(::Type{T}, scale, l::T) where T <: Length = T(scale * l.value)
scale_component(::Type{T}, scale, l::Add) where T <: Length =
        scale_component(T, scale, l.a) + scale_component(T, scale, l.b)

# Interpretation of bare numbers
# ------------------------------

x_measure(a::Measure) = a
x_measure(a::T) where T = Length{:cx, T}(a)

y_measure(a::Measure) = a
y_measure(a::T) where T = Length{:cy, T}(a)

size_x_measure(a::Measure) = a
size_x_measure(a::T) where T = Length{:sx, T}(a)
size_y_measure(a::Measure) = a
size_y_measure(a::T) where T = Length{:sy, T}(a)

x_measure(a::Vector{T}) where T <: Measure = a
x_measure(a::Vector) = Measure[x_measure(x) for x in a]

y_measure(a::Vector{T}) where T <: Measure = a
y_measure(a::Vector) = Measure[y_measure(y) for y in a]

size_x_measure(a::Vector{T}) where T <: Measure = a
size_x_measure(a::Vector) = Measure[size_x_measure(x) for x in a]
size_y_measure(a::Vector{T}) where T <: Measure = a
size_y_measure(a::Vector) = Measure[size_y_measure(y) for y in a]

size_measure(a::Measure) = a
size_measure(a) = a * mm

x_measure(a::Missing) = x_measure(NaN)
y_measure(a::Missing) = y_measure(NaN)

size_x_measure(a::Missing) = size_x_measure(NaN)
size_y_measure(a::Missing) = size_y_measure(NaN)

# Higher-order measures
# ---------------------

# Compute the union of two bounding boxes.
#
# In other words, given two bounding boxes, return a new bounding box that
# contains both.
#
# Unfortunately this is in general uncomputable without knowing the absolute
# size of the parent canvas which may be passed in via the last two parameters.
# If not passed, this throws an error if they would have been required.
#
function union(a::BoundingBox, b::BoundingBox, units=nothing, parent_abs_width=nothing, parent_abs_height=nothing)
    (a.width == Measure() || a.height == Measure()) && return b
    (b.width == Measure() || b.height == Measure()) && return a
    x0 = min(a.x0, b.x0)
    y0 = min(a.y0, b.y0)
    x1 = max(a.x0 + a.width, b.x0 + b.width)
    y1 = max(a.y0 + a.height, b.y0 + b.height)
    # Check whether we had any problematic computations
    for m in (x0,y0,x1,y1)
        # Pure absolute or pure relative points are fine. When they are mixed,
        # there are problems
        if !isabsolute(m) && m.abs != 0.0
            units == nothing || parent_abs_width == nothing || parent_abs_height == nothing &&
                    error("""Bounding boxes are uncomputable without knowledge of the
                        absolute dimensions of the top canvase due to mixing of relative
                        and absolute coordinates. Either pass the dimension as a parameter
                        or restrict the context to one kind of coordinates.""")
            parent_box = AbsoluteBox(0.0,0.0,parent_abs_width,parent_abs_height)
            abb = union(absolute_units(a,IdentityTransform(),units,parent_box),
                        absolute_units(b,IdentityTransform(),units,parent_box))
            return BoundingBox(Measure(;abs = abb.x0), Measure(;abs = abb.x0),
                               Measure(;abs = abb.width), Measure(;abs = abb.height))
        end
    end
    return BoundingBox(x0, y0, x1 - x0, y1 - y0)
end

function union(a::AbsoluteBox, b::AbsoluteBox)
    (a.width == 0.0 || a.height == 0.0) && return b
    (b.width == 0.0 || b.height == 0.0) && return a
    x0 = min(a.x0, b.x0)
    y0 = min(a.y0, b.y0)
    x1 = max(a.x0 + a.width, b.x0 + b.width)
    y1 = max(a.y0 + a.height, b.y0 + b.height)
    return AbsoluteBox(x0, y0, x1 - x0, y1 - y0)
end

# The same type-signature is used for a box used to assign
# a custom coordinate system to a canvas.

struct UnitBox{S,T,U,V}
    x0::S
    y0::T
    width::U
    height::V

    leftpad::AbsoluteLength
    rightpad::AbsoluteLength
    toppad::AbsoluteLength
    bottompad::AbsoluteLength
end

"""
    UnitBox(x0, y0, width, height; leftpad=0mm, rightpad=0mm, toppad=0mm, bottompad=0mm)

Specifies the coordinate system for a [`context`](@ref), with origin `x0`, `y0`, plus `width`, `height`.
"""
UnitBox(x0, y0, width, height; leftpad=0mm, rightpad=0mm, toppad=0mm, bottompad=0mm) =
        UnitBox{typeof(x0), typeof(y0), typeof(width), typeof(height)}(
            x0, y0, width, height, leftpad, rightpad, toppad, bottompad)

"""
    UnitBox(width, height; leftpad=0mm, rightpad=0mm, toppad=0mm, bottompad=0mm)

Specifies the coordinate system for a [`context`](@ref), with origin `0, 0` plus `width`, `height`.
"""
function UnitBox(width, height; leftpad=0mm, rightpad=0mm, toppad=0mm, bottompad=0mm)
    S, T = typeof(width), typeof(height)
    UnitBox{S,T,S,T}(zero(S), zero(T), width, height, leftpad, rightpad, toppad, bottompad)
end

"""
    UnitBox() = UnitBox(0.0, 0.0, 1.0, 1.0)
"""
UnitBox() = UnitBox(0.0, 0.0, 1.0, 1.0)

# copy with substitution
UnitBox(units::UnitBox;
                 x0=nothing,
                 y0=nothing,
                 width=nothing,
                 height=nothing,
                 leftpad=nothing,
                 rightpad=nothing,
                 toppad=nothing,
                 bottompad=nothing) =
        UnitBox(ifelse(x0 === nothing,     units.x0,     x0),
            ifelse(y0 === nothing,     units.y0,     y0),
            ifelse(width === nothing,  units.width,  width),
            ifelse(height === nothing, units.height, height),
            leftpad   = ifelse(leftpad === nothing,   units.leftpad,   leftpad),
            rightpad  = ifelse(rightpad === nothing,  units.rightpad,  rightpad),
            toppad    = ifelse(toppad === nothing,    units.toppad,    toppad),
            bottompad = ifelse(bottompad === nothing, units.bottompad, bottompad))

Measures.width(units::UnitBox) = units.width
Measures.height(units::UnitBox) = units.height

ispadded(units::UnitBox) = units.leftpad != 0mm || units.rightpad != 0mm ||
        units.toppad != 0mm || units.bottompad != 0mm

isxflipped(units::UnitBox{S, T, U, V}) where {S, T, U, V} = units.width < zero(U)
isyflipped(units::UnitBox{S, T, U, V}) where {S, T, U, V} = units.height < zero(V)
hasunits(::Type, x::Measure) = false
hasunits(::Type{Length{u}}, x::Length{u, T}) where {u, T} = true
hasunits(T::Type, x::Measures.BinaryOp) = hasunits(T, x.a) || hasunits(T, x.b)

# Canvas Transforms
# -----------------

# Transform matrix in absolute coordinates

abstract type Transform end

struct IdentityTransform <: Transform
end

struct MatrixTransform <: Transform
    M::Matrix{Float64}
end

combine(a::IdentityTransform, b::IdentityTransform) = a
combine(a::IdentityTransform, b::MatrixTransform) = b
combine(a::MatrixTransform, b::IdentityTransform) = a
combine(a::MatrixTransform, b::MatrixTransform) = MatrixTransform(a.M * b.M)

# Rotation about a point.
struct Rotation{P <: Vec}
    theta::Float64
    offset::P
end

Rotation(theta::Number, offset::XYTupleOrVec) =
        Rotation(convert(Float64, theta), (x_measure(offset[1]), y_measure(offset[2])))

"""
    Rotation(θ, x, y)

Rotate all forms in context around point `(x,y)` by angle `θ` in radians.
"""
Rotation(theta::Number, offset_x, offset_y) =
        Rotation(convert(Float64, theta), (x_measure(offset_x), y_measure(offset_y)))

"""
    Rotation(θ)

`Rotation(θ)=Rotation(θ, 0.5w, 0.5h)`
"""
Rotation(theta::Number) = Rotation{Vec2}(convert(Float64, theta), (0.5w, 0.5h))
Rotation() = Rotation(0.0, (0.5w, 0.5h))

copy(rot::Rotation) = Rotation(rot)

function convert(::Type{Transform}, rot::Rotation)
    if rot.theta == 0.0
        return IdentityTransform()
    else
        ct = cos(rot.theta)
        st = sin(rot.theta)
        x0 = rot.offset[1] - (ct * rot.offset[1] - st * rot.offset[2])
        y0 = rot.offset[2] - (st * rot.offset[1] + ct * rot.offset[2])
        return MatrixTransform([ct  -st  x0.value
                                st   ct  y0.value
                                0.0 0.0  1.0])
    end
end

# Mirror about a point at a given angle
mutable struct Mirror
    theta::Float64
    point::Vec
end

"""
    Mirror(θ, x, y)

Mirror line passing through point `(x,y)` at angle `θ` (in radians).
"""
Mirror(theta::Number, offset_x, offset_y) = Mirror(convert(Float64, theta), (offset_x, offset_y))
Mirror(theta::Number, offset::XYTupleOrVec) =
        Mirror(convert(Float64, theta), (x_measure(offset[1]), y_measure(offset[2])))

"""
    Mirror(θ)

`Mirror(θ)=Mirror(θ, 0.5w, 0.5h)`
"""
Mirror(theta::Number) = Mirror(convert(Float64, theta), 0.5w, 0.5h)
Mirror() = Mirror(0.0, (0.5w, 0.5h))

# copy constructor
Mirror(mir::Mirror) = Mirror(copy(mir.theta), copy(mir.offset))

function convert(::Type{Transform}, mir::Mirror)
    n = [cos(mir.theta), sin(mir.theta)]
    x0 = mir.point[1]
    y0 = mir.point[2]

    offset = (2I - 2n*n') * [x0.value, y0.value]
    scale  = (2n*n' - I)
    M = vcat(hcat(scale, offset), [0 0 1])

    return MatrixTransform(M)
end

copy(mir::Mirror) = Mirror(mir)

# Resolution
# ----------

resolve(box::AbsoluteBox, units::UnitBox, t::Transform, a::Length) = resolve(box, a)
resolve_position(box::AbsoluteBox, units::UnitBox, t::Transform, a::Length{:cx}) =
        ((a.value - units.x0) / width(units)) * box.a[1]
resolve(box::AbsoluteBox, units::UnitBox, t::Transform, a::Length{:cx}) =
        abs(a.value / width(units)) * box.a[1]
resolve_position(box::AbsoluteBox, units::UnitBox, t::Transform, a::Length{:cy}) =
        ((a.value - units.y0) / height(units)) * box.a[2]
resolve(box::AbsoluteBox, units::UnitBox, t::Transform, a::Length{:cy}) =
        abs(a.value / height(units)) * box.a[2]

resolve(box::AbsoluteBox, units::UnitBox, t::Transform, a::Length{:sx}) =
        a.value / width(units) * box.a[1]

resolve(box::AbsoluteBox, units::UnitBox, t::Transform, a::Length{:sy}) =
        a.value / height(units) * box.a[2]

function resolve(box::AbsoluteBox, units::UnitBox, t::Transform, p::Vec2)
    xy = (resolve_position(box, units, t, p[1]) + box.x0[1],
          resolve_position(box, units, t, p[2]) + box.x0[2])
    return xy
end

function resolve(box::AbsoluteBox, units::UnitBox, t::MatrixTransform, p::Vec2)
    x = resolve_position(box, units, t, p[1]) + box.x0[1]
    y = resolve_position(box, units, t, p[2]) + box.x0[2]
    xy = t.M * [x.value, y.value, 1]
    return (xy[1]mm, xy[2]mm)
end

resolve(box::AbsoluteBox, units::UnitBox, t::Transform, a::BoundingBox) =
        BoundingBox(resolve(box, units, IdentityTransform(), a.x0),
            (resolve(box, units, t, a.a[1]), resolve(box, units, t, a.a[2])))

resolve(box::AbsoluteBox, units::UnitBox, t::Transform, a::Rotation) =
        Rotation(a.theta, resolve(box, units, t, a.offset))

resolve(box::AbsoluteBox, units::UnitBox, t::Transform, a::Mirror) =
        Mirror(a.theta, resolve(box, units, t, a.point))

function resolve(box::AbsoluteBox, units::UnitBox, t::Transform, u::UnitBox)
    if !ispadded(u)
        return u
    else
        leftpad   = resolve(box, units, t, u.leftpad)
        rightpad  = resolve(box, units, t, u.rightpad)
        toppad    = resolve(box, units, t, u.toppad)
        bottompad = resolve(box, units, t, u.bottompad)

        # just give up trying to pad the units if it's impossible
        if leftpad + rightpad >= box.a[1] ||
           toppad + bottompad >= box.a[2]
            return UnitBox(u.x0, u.y0, u.width, u.height)
        end

        width = u.width * (box.a[1] / (box.a[1] - leftpad - rightpad))
        height = u.height * (box.a[2] / (box.a[2] - toppad - bottompad))
        x0 = u.x0 - width * (leftpad / box.a[1])
        y0 = u.y0 - height * (toppad / box.a[2])

        return UnitBox(x0, y0, width, height)
    end
end

# Equivalent to the resolve functions in Measures, but pass through the `units`
# and `transform` parameters.
resolve(box::AbsoluteBox, units::UnitBox, t::Transform, x::Neg) =
        -resolve(box, units, t, x.a)
resolve(box::AbsoluteBox, units::UnitBox, t::Transform, x::Add) =
        resolve(box, units, t, x.a) + resolve(box, units, t, x.b)
resolve(box::AbsoluteBox, units::UnitBox, t::Transform, x::Mul) =
        resolve(box, units, t, x.a) * x.b
resolve(box::AbsoluteBox, units::UnitBox, t::Transform, x::Div) =
        resolve(box, units, t, x.a) / x.b
resolve(box::AbsoluteBox, units::UnitBox, t::Transform, x::Min) =
        min(resolve(box, units, t, x.a), resolve(box, units, t, x.b))
resolve(box::AbsoluteBox, units::UnitBox, t::Transform, x::Max) =
        max(resolve(box, units, t, x.a), resolve(box, units, t, x.b))

resolve_position(box::AbsoluteBox, units::UnitBox, t::Transform, a) = resolve(box, units, t, a)
resolve_position(box::AbsoluteBox, units::UnitBox, t::Transform, op::Add) =
        resolve_position(box, units, t, op.a) + resolve_position(box, units, t, op.b)

# Shear about a point at a given angle
struct Shear
    shear::Float64
    theta::Float64
    point::Vec
end

"""
    Shear(s, θ, x, y)

Shear line passing through point `(x,y)` at angle `θ` (in radians), with shear `s`.
"""
Shear(shear::Number, theta::Number, offset_x, offset_y) =
        Shear(float(shear), float(theta), (offset_x, offset_y))

Shear(shear::Number, theta::Number, offset::XYTupleOrVec) =
        Shear(float(shear), float(theta), (x_measure(offset[1]), y_measure(offset[2])))

"""
    Shear(s, θ)

`Shear(s, θ)=Shear(s, θ, 0.5w, 0.5h)`
"""
Shear(shear::Number, theta::Number) = Shear(float(shear), float(theta), 0.5w, 0.5h)

function convert(::Type{Transform}, shear::Shear)
    x, y = shear.point[1].value, shear.point[2].value
    ct, st = cos(shear.theta), sin(shear.theta)
    M1 = [1.0 0 x; 0 1 y; 0 0 1 ]
    M2 = I + shear.shear .* [ct*st -ct*ct 0; st*st -ct*st 0; 0 0 0.0]
    M3 = [1.0 0 -x; 0 1 -y; 0 0 1 ]
    M = M1*M2*M3
    return  MatrixTransform(M)
end

resolve(box::AbsoluteBox, units::UnitBox, t::Transform, a::Shear) =
            Shear(a.shear, a.theta, resolve(box, units, t, a.point))

# for transforming quadratic to cubic
controlpnt(anchor::XYTupleOrVec, control::XYTupleOrVec) = 0.333*anchor + 0.667*control


text:
using DiffRules
using Test

import SpecialFunctions, NaNMath, LogExpFunctions
import Random
Random.seed!(1)

function finitediff(f, x)
    ϵ = cbrt(eps(typeof(x))) * max(one(typeof(x)), abs(x))
    return (f(x + ϵ) - f(x - ϵ)) / (ϵ + ϵ)
end

@testset "DiffRules" begin
@testset "check rules" begin

non_diffeable_arg_functions = [(:Base, :rem2pi, 2), (:Base, :ldexp, 2), (:Base, :ifelse, 3)]

for (M, f, arity) in DiffRules.diffrules(; filter_modules=nothing)
    (M, f, arity) ∈ non_diffeable_arg_functions && continue
    if arity == 1
        @test DiffRules.hasdiffrule(M, f, 1)
        deriv = DiffRules.diffrule(M, f, :goo)
        modifier = if f in (:asec, :acsc, :asecd, :acscd, :acosh, :acoth)
            1.0
        elseif f === :log1mexp
            -1.0
        elseif f === :log2mexp
            -0.5
        else
            0.0
        end
        @eval begin
            let
                goo = rand() + $modifier
                @test isapprox($deriv, finitediff($M.$f, goo), rtol=0.05)
                # test for 2pi functions
                if "mod2pi" == string($M.$f)
                    goo = 4pi + $modifier
                    @test NaN === $deriv
                end
            end
        end
    elseif arity == 2
        @test DiffRules.hasdiffrule(M, f, 2)
        derivs = DiffRules.diffrule(M, f, :foo, :bar)
        @eval begin
            let
                if "mod" == string($M.$f)
                    foo, bar = rand() + 13, rand() + 5 # make sure x/y is not integer
                else
                    foo, bar = rand(1:10), rand()
                end
                dx, dy = $(derivs[1]), $(derivs[2])
                if !(isnan(dx))
                    @test isapprox(dx, finitediff(z -> $M.$f(z, bar), float(foo)), rtol=0.05)
                end
                if !(isnan(dy))
                    @test isapprox(dy, finitediff(z -> $M.$f(foo, z), bar), rtol=0.05)
                end
            end
        end
    elseif arity == 3
        #=
        @test DiffRules.hasdiffrule(M, f, 3)
        derivs = DiffRules.diffrule(M, f, :foo, :bar, :goo)
        @eval begin
            foo, bar, goo = randn(3)
            dx, dy, dz = $(derivs[1]), $(derivs[2]), $(derivs[3])
            if !(isnan(dx))
                @test isapprox(dx, finitediff(x -> $M.$f(x, bar, goo), foo), rtol=0.05)
            end
            if !(isnan(dy))
                @test isapprox(dy, finitediff(y -> $M.$f(foo, y, goo), bar), rtol=0.05)
            end
            if !(isnan(dz))
                @test isapprox(dz, finitediff(z -> $M.$f(foo, bar, z), goo), rtol=0.05)
            end
        end
        =#
    end
end

# Treat rem2pi separately because of its non-numeric second argument:
derivs = DiffRules.diffrule(:Base, :rem2pi, :x, :y)
for xtype in [:Float64, :BigFloat, :Int64]
    for mode in [:RoundUp, :RoundDown, :RoundToZero, :RoundNearest]
        @eval begin
            let
                x = $xtype(rand(1 : 10))
                y = $mode
                dx, dy = $(derivs[1]), $(derivs[2])
                @test isapprox(dx, finitediff(z -> rem2pi(z, y), float(x)), rtol=0.05)
                @test isnan(dy)
            end
        end
    end
end

# Treat ldexp separately because of its integer second argument:
derivs = DiffRules.diffrule(:Base, :ldexp, :x, :y)
for xtype in [:Float64, :BigFloat]
    for ytype in [:Integer, :UInt64, :Int64]
        @eval begin
            let
                x = rand($xtype)
                y = $ytype(rand(1 : 10))
                dx, dy = $(derivs[1]), $(derivs[2])
                @test isapprox(dx, finitediff(z -> ldexp(z, y), x), rtol=0.05)
                @test isnan(dy)
            end
        end
    end
end

end

    @testset "diffrules" begin
        rules = @test_deprecated(DiffRules.diffrules())
        @test Set(M for (M, _, _) in rules) == Set((:Base, :SpecialFunctions, :NaNMath))

        rules = DiffRules.diffrules(; filter_modules=nothing)
        @test Set(M for (M, _, _) in rules) == Set((:Base, :SpecialFunctions, :NaNMath, :LogExpFunctions))

        rules = DiffRules.diffrules(; filter_modules=(:Base, :LogExpFunctions))
        @test Set(M for (M, _, _) in rules) == Set((:Base, :LogExpFunctions))
    end
end

# Test ifelse separately as first argument is boolean
#=
@test DiffRules.hasdiffrule(:Base, :ifelse, 3)
derivs = DiffRules.diffrule(:Base, :ifelse, :foo, :bar, :goo)
for cond in [true, false]
    @eval begin
        foo = $cond
        bar, gee = randn(2)
        dx, dy, dz = $(derivs[1]), $(derivs[2]), $(derivs[3])
        @test isapprox(dy, finitediff(y -> ifelse(foo, y, goo), bar), rtol=0.05)
        @test isapprox(dz, finitediff(z -> ifelse(foo, bar, z), goo), rtol=0.05)
    end
end
=#


text:
using ToyFHE
using OffsetArrays
using Test

params = BFVParams(
    7, # plaintext modulus
    ; eval_mult_count = 3
)

plain = zero(plaintext_space(params))
plain[0] = 2

# Generate a keypair and a corresponding relin key
kp1 = keygen(params)
ek = keygen(EvalMultKey, kp1.priv)

c1 = encrypt(kp1, plain)

@test decrypt(kp1, c1)[0] == 2
let c1squared = c1*c1
    @test decrypt(kp1, c1squared)[0] == 4
    cswitch = keyswitch(ek, c1squared)
    @test length(cswitch.cs) == 2
    @test decrypt(kp1, cswitch)[0] == 4
    @test decrypt(kp1, cswitch*c1)[0] == 1
end


text:
struct FlameColors
    colors::Vector{RGB{N0f8}}
    colorbg::RGB{N0f8}
    colorfont::RGB{N0f8}
    colorsrt::Vector{RGB{N0f8}}
    colorsgc::Vector{RGB{N0f8}}
end

"""
    fcolor = FlameColors(n::Integer=2;
                         colorbg=colorant"white", colorfont=colorant"black",
                         colorsrt=colorant"crimson", colorsgc=colorant"orange")

Choose a set of colors for rendering a flame graph. There are several special colors:

- `colorbg` is the background color
- `colorfont` is used when annotating stackframes with text
- `colorsrt` highlights [runtime dispatch](https://discourse.julialang.org/t/dynamic-dispatch/6963), typically a costly process
- `colorsgc` highlights garbage-collection events

`n` specifies the number of "other" colors to choose when one of the above is not relevant.
`FlameColors` chooses `2n` colors: the first `n` colors for odd depths in the stacktrace and
the last `n` colors for even depths in the stacktrace. Consequently, different stackframes will typically
be distinguishable from one another by color.

The highlighting can be disabled by passing `nothing` or a zero-element vector
for `colorsrt` or `colorsgc`. When a single color is passed for `colorsrt` or
`colorsgc`, this method generates four variant colors slightly different from
the specified color. `colorsrt` or `colorsgc` can also be specified as multiple
colors with a vector. The first half of the vector is for odd depths and the
second half is for even depths. By using a one-element vector instead of a
single color, the specified color is always used.

While the return value is a `struct`, it is callable and can be used as the
`fcolor` input for `flamepixels`.
"""
function FlameColors(n::Integer=2;
                     colorbg=colorant"white", colorfont=colorant"black",
                     colorsrt=colorant"crimson", colorsgc=colorant"orange")
    seeds = [colorbg, colorfont]
    function make_variations(color)
        color === nothing && return RGB{N0f8}[]
        isa(color, AbstractVector) && return color
        c = convert(LCHab, color)
        a = LCHab(min(c.l+2.5, 100.0), c.c+10.0, c.h-7.5)
        b = LCHab(max(c.l-2.5, 0.0), max(c.c-10.0, 0.0), c.h+7.5)
        RGB{N0f8}.(range(a, stop=b, length=4))
    end
    colorsrt = make_variations(colorsrt)
    colorsgc = make_variations(colorsgc)
    append!(seeds, colorsrt)
    append!(seeds, colorsgc)
    nseeds = length(seeds)
    colors = distinguishable_colors(2n+nseeds, seeds,
                                    transform=c->deuteranopic(c, 0.95),
                                    lchoices=Float64[65, 80],
                                    cchoices=Float64[10, 55],
                                    hchoices=range(10, stop=350, length=18))[nseeds+1:end]
    sort!(colors, by=c->colordiff(c, colorbg))
    return FlameGraphs.FlameColors(colors, colorbg, colorfont, colorsrt, colorsgc)
end

const default_colors = FlameColors()

function (colors::FlameColors)(s::Symbol)
    s === :bg && return colors.colorbg
    s === :font && return colors.colorfont
    throw(ArgumentError("unrecognized color id ", s))
end

function (colors::FlameColors)(nextidx, j::Integer, data)
    idx = nextidx[j]
    nextidx[j] = idx + 1
    if length(colors.colorsrt) > 0 && (data.status & runtime_dispatch) != 0
        colorvec = colors.colorsrt
    elseif length(colors.colorsgc) > 0 && (data.status & gc_event) != 0
        colorvec = colors.colorsgc
    else
        colorvec = colors.colors
    end
    n = length(colorvec)
    m = (n + 1)÷2
    return colorvec[mod1(mod1(idx, m) + iseven(j) * m, n)]
end

"""
    img = flamepixels(g; kwargs...)

Return a flamegraph as a matrix of RGB colors. The first dimension corresponds to cost,
the second dimension to depth in the call stack.

See also [`flametags`](@ref).
"""
flamepixels(g::Node; kwargs...) = flamepixels(default_colors, g; kwargs...)

"""
    img = flamepixels(fcolor, g; costscale=nothing)

Return a flamegraph as a matrix of RGB colors, customizing the color choices.

## fcolor

`fcolor` is a function that returns the
color used for the current item in the call stack.
See [`FlameColors`](@ref) for the default implementation of `fcolor`.

If you provide a custom `fcolor`, it must support the following API:

    colorbg = fcolor(:bg)
    colorfont = fcolor(:font)

must return the background and font colors.

    colornode = fcolor(nextidx::Vector{Int}, j, data::NodeData)

chooses the color for the node represented by `data` (see [`NodeData`](@ref)).
`j` corresponds to depth in the call stack and `nextidx[j]` holds the state for the next
color choice.
In general, if you have a list of colors, `fcolor` should cycle `nextidx[j]` to ensure that
the next call to `fcolor` with this `j` moves on to the next color.
(However, you may not want to increment `nextidx[j]` if you are choosing the color by some
means other than cycling through a list.)

By accessing `data.sf`, you can choose to color individual nodes based on the identity of
the stackframe.

## costscale

`costscale` can be used to limit the size of `img` when profiling collected a large number of stacktraces.
The size of the first dimension of `img` is proportional to the total number of
stacktraces collected during profiling. `costscale` is the constant of proportionality;
for example, setting `costscale=0.2` would mean that `size(img, 1)` would be approximately
1/5 the number of stacktraces collected by the profiler. The default value of `nothing`
imposes an upper bound of approximately 1000 pixels along the first dimension, with
`costscale=1` chosen if the number of samples is less than 1000.
"""
function flamepixels(fcolor, g::Node; costscale=nothing)
    ndata = g.data
    w = length(ndata.span)
    if costscale === nothing
        costscale = w < 10^3 ? 1.0 : 10^3/w
    end
    h = depth(g)
    img = fill(fcolor(:bg), round(Int, w*costscale), h)
    nextidx = fill(1, h)
    img[scale(ndata.span, costscale), 1] .= fcolor(nextidx, 1, ndata)
    return flamepixels!(fcolor, img, g, 2, nextidx, costscale)
end

function flamepixels!(fcolor, img, g, j, nextidx, costscale)
    for c in g
        ndata = c.data
        img[scale(ndata.span, costscale), j] .= fcolor(nextidx, j, ndata)
        flamepixels!(fcolor, img, c, j+1, nextidx, costscale)
    end
    return img
end

"""
    tagimg = flametags(g, img)

From a flame graph `g`, generate an array `tagimg` with the same axes as `img`,
encoding the stackframe represented by each pixel of `img`.

See [`flamepixels`](@ref) to generate `img`.
"""
function flametags(g, img)
    tags = fill(1, axes(img))
    ndata = g.data
    sflist = StackFrame[StackTraces.UNKNOWN]
    sf2tag = Dict{StackFrame,Int}(StackTraces.UNKNOWN=>1)
    tags[:,1] .= tagidx(ndata.sf, sf2tag, sflist)
    costscale = size(img, 1) / length(ndata.span)
    flametags!(tags, g, sf2tag, sflist, 2, costscale)
    return IndirectArray(tags, sflist)
end

function flametags!(tags, parent, sf2tag, sflist, level, costscale)
    for c in parent
        ndata = c.data
        tags[scale(ndata.span, costscale), level] .= tagidx(ndata.sf, sf2tag, sflist)
        flametags!(tags, c, sf2tag, sflist, level+1, costscale)
    end
end

scale(rng::UnitRange, costscale) = max(1, round(Int, costscale*first(rng))):round(Int, costscale*last(rng))

function tagidx(sf, sf2tag, sflist)
    idx = get(sf2tag, sf, length(sflist)+1)
    if idx > length(sflist)
        sf2tag[sf] = idx
        push!(sflist, sf)
    end
    return idx
end


text:
import Base.Iterators:flatten

using CodecLz4;
using ZfpCompression;

width = 14
height = 6
precision = 8

pixels = Float32(100.0) * randn(Float32, width, height)
mask = map(x -> x > -90, pixels)

compressed_pixels = zfp_compress(pixels, precision=precision)
compressed_mask = transcode(LZ4HCCompressor, collect(flatten(UInt8.(mask))))

println("compressed pixels size:", length(compressed_pixels))
println("compressed mask size:", length(compressed_mask))

compressed_mask = lz4_hc_compress(collect(flatten(UInt8.(mask))))
println("compressed mask size:", length(compressed_mask))

decomp = zfp_decompress(compressed_pixels)

display(pixels)
display(decomp)

# 1D spectrum compression
depth = 768
spectrum = Float32(100.0) * randn(Float32, depth)
precision = 24

compressed_spectrum = zfp_compress(spectrum, precision=precision)
println("compressed spectrum size: ", length(compressed_spectrum))

decomp = zfp_decompress(compressed_spectrum)

display(spectrum)
display(decomp)


text:
@testset "IntRangeDict_basic" begin
    x = IntRangeDict{Int, Int}()
    @test x[2] == []

    push!(x[3:5], 1)
    @test x[2] == []
    @test x[3] == [1]
    @test x[4] == [1]
    @test x[5] == [1]
    @test x[6] == []

    push!(x[7:7], 2)
    @test x[5] == [1]
    @test x[6] == []
    @test x[7] == [2]
    @test x[8] == []

    push!(x[1:2], 3)
    @test x[0] == []
    @test x[1] == [3]
    @test x[2] == [3]
    @test x[3] == [1]

    push!(x[1:7], 4)
    @test x[1] == [3, 4]
    @test x[3] == [1, 4]
    @test x[6] == [4]
    @test x[8] == []
end

@testset "IntRangeDict_save_and_load" begin
    x = IntRangeDict{Int, Int}()
    buf = IOBuffer()

    push!(x[4:5], 1)
    push!(x[7:8], 2)
    save(buf, x)

    seekstart(buf)

    y = IntRangeDict{Int, Int}(buf)
    @test y[3] == []
    @test y[5] == [1]
    @test y[6] == []
    @test y[7] == [2]
    @test y[9] == []
end


text:
using Base.BinaryPlatforms

# Our platform augmentation simply inspects the environment
function augment_platform!(p::Platform)
    # If the `flooblecrank` tag is already set, don't auto-detect!
    if haskey(p, "flooblecrank")
        return p
    end

    # If the tag is not set, autodetect through magic (in this case, checking environment variables)
    flooblecrank_status = get(ENV, "FLOOBLECRANK", "disengaged")
    if flooblecrank_status == "engaged"
        p["flooblecrank"] = "engaged"
    else
        p["flooblecrank"] = "disengaged"
    end
    return p
end


text:
"""
Defines a basisfunction which can be called for each event, defining the time-expanded basis kernel
$(TYPEDEF)
$(TYPEDSIGNATURES)
$(FIELDS)
# Examples
```julia-repl
julia>  b = BasisFunction(kernelfunction,["hrf"],range(0,(length(kernelfunction([0, 1]))-1)*TR,step=TR),"hrf_kernel","basis_A",0)
```
"""
struct BasisFunction
    "a design-matrix kernel function used to timeexpand, given a timepoint in 'sample' timeunits"
    kernel::Function
    "vector of names along columns of kernel-output"
    colnames::AbstractVector
    " vector of times along rows of kernel-output (in seconds)"
    times::AbstractVector
    "type of basisfunction (only for bookkeeping)"
    type::String
    "name of the event, random 1:1000 if unspecified"
    name::String
    "by how many samples do we need to shift the event onsets? This number is determined by how many 'negative' timepoints the basisfunction defines"
    shiftOnset::Integer
end

function Base.show(io::IO, obj::BasisFunction)
    println(io, "name: $(obj.name)")
    println(io, "colnames: $(obj.colnames)")
    println(io, "kernel: $(obj.type)")
    println(io, "times: $(obj.times)")
    println(io, "shiftOnset: $(obj.shiftOnset)")
end

"""
$(SIGNATURES)
Generate a sparse FIR basis around the *τ* timevector at sampling rate *sfreq*. This is useful if you cannot make any assumptions on the shape of the event responses. If unrounded events are supplied, they are split between samples. E.g. event-latency = 1.2 will result in a "0.8" and a "0.2" entry.

# Examples
Generate a FIR basis function from -0.1s to 0.3s at 100Hz
```julia-repl
julia>  f = firbasis([-0.1,0.3],100)
```
Evaluate at an event occuring at sample 103.3
```julia-repl
julia>  f(103.3)
```

"""
function firbasis(τ,sfreq,name::String)
    τ = round_times(τ,sfreq)
    times =range(τ[1],stop=τ[2],step=1 ./sfreq)

    kernel=e->firkernel(e,times)
    type = "firkernel"

    shiftOnset = Int64(floor(τ[1] * sfreq))

    return BasisFunction(kernel,times,times,type,name,shiftOnset)
end
# cant multiple dispatch on optional arguments
#firbasis(;τ,sfreq)           = firbasis(τ,sfreq)
firbasis(;τ,sfreq,name)      = firbasis(τ,sfreq,name)
firbasis(τ,sfreq)            = firbasis(τ,sfreq,"basis_"*string(rand(1:10000)))

"""
$(SIGNATURES)
Calculate a sparse firbasis
# Examples

```julia-repl
julia>  f = firkernel(103.3,range(-0.1,step=0.01,stop=0.31))
```
"""
function firkernel(e,times)
    @assert ndims(e) <= 1 #either single onset or a row vector where we will take the first one :)
    if size(e,1) > 1
        # XXX we will soon assume that the second entry would be the duration
        e = Float64(e[1])
    end
    e = [1 .- (e .% 1)  e .% 1]
    e[isapprox.(e,0,atol = 1e-15)] .= 0
    ksize=length(times) # kernelsize

    kernel = spdiagm(ksize+1,ksize,0 => repeat([e[1]],ksize), -1 => repeat([e[2]],ksize))

    return(kernel)

end

"""
$(SIGNATURES)
Generate a Hemodynamic-Response-Functio (HRF) basis with inverse-samplingrate "TR" (=1/FS)

Optional Parameters p:
                                                           defaults
                                                          {seconds}
        p(1) - delay of response (relative to onset)          6
        p(2) - delay of undershoot (relative to onset)       16
        p(3) - dispersion of response                         1
        p(4) - dispersion of undershoot                       1
        p(5) - ratio of response to undershoot                6
        p(6) - onset {seconds}                                0
        p(7) - length of kernel {seconds}                    32

# Examples
Generate a HRF basis function object with Sampling rate 1/TR. And evaluate it at an event occuring at TR 103.3 with duration of 4.1 TRs
```julia-repl
julia>  f = hrfbasis(2.3)
julia>  f(103.3,4.1)

```

"""
function hrfbasis(TR::Float64;parameters= [6. 16. 1. 1. 6. 0. 32.],name::String="basis_"*string(rand(1:10000)))
    # Haemodynamic response function adapted from SPM12b "spm_hrf.m"
    # Parameters:
    #                                                           defaults
    #                                                          {seconds}
    #        p(1) - delay of response (relative to onset)          6
    #        p(2) - delay of undershoot (relative to onset)       16
    #        p(3) - dispersion of response                         1
    #        p(4) - dispersion of undershoot                       1
    #        p(5) - ratio of response to undershoot                6
    #        p(6) - onset {seconds}                                0
    #        p(7) - length of kernel {seconds}                    32
    kernel=e->hrfkernel(e,TR,parameters)
    type = "hrfkernel"
    return BasisFunction(kernel,["hrf"],range(0,(length(kernel([0, 1]))-1)*TR,step=TR),type,name,0)
end

"""
$(SIGNATURES)
Calculate a HRF kernel. Input e can be [onset duration]
# Examples

```julia-repl
julia>  f = hrfkernel(103.3,2.3,[6. 16. 1. 1. 6. 0. 32.])
```
"""
function hrfkernel(e,TR,p)
    @assert ndims(e) <= 1 #either single onset or a row vector where we will take the first one :)

    # code adapted from SPM12b
    mt = 32
    dt  = TR/mt;
    firstElem = Int(round((1-(e[1]%1)).*mt))

    if size(e,1) == 2 && e[2]>0.
        duration = e[2]
    else
        duration = 1. /mt
    end
    duration_mt = Int(ceil(duration .* mt))
    box_mt = [zeros(firstElem)' ones(duration_mt)']'
    #box_mt = box_mt ./ (sum(box_mt))
    # u   = [0:ceil(p(7)/dt)] - p(6)/dt;
    u   = range(0,stop = ceil((p[7]+duration)/dt)) .- p[6]/dt;

    #hrf = spm_Gpdf(u,p(1)/p(3),dt/p(3))
    # Note the inverted scale parameter compared to SPM12.
    g1  = Gamma(p[1] ./ p[3],p[3] ./ dt)
     #spm_Gpdf(u,p(2)/p(4),dt/p(4))
    g2  = Gamma(p[2] ./ p[4],p[4] ./ dt);
    # g1 - g2/p(5);
    hrf = pdf.(g1,u) .- pdf.(g2,u) ./ p[5]

    # hrf = hrf([0:floor(p(7)/RT)]*fMRI_T + 1);
    hrf = hrf ./ sum(hrf);
    hrf = conv(box_mt,hrf)'
    #println(box_mt)
    hrf = hrf[((range(0,stop=Int(floor(p[7] ./ TR + duration)))*mt)) .+ 1];

    return(SparseMatrixCSC(sparse(hrf)))

end


text:
include("main.jl")

using HDF5
using Images
using ImageView
using Plots
using MLJBase

"""
    Function to load flattened toy data to test the implementation.
"""
function load_data()
    X_train = Float64.(h5read("train_catvnoncat.h5", "train_set_x"))
    #X_train = reshape(X_train, :, size(X_train, 4))

    y_train = Int64.(h5read("train_catvnoncat.h5", "train_set_y"))
    #y_train = reshape(y_train, (1, size(y_train, 1)))

    X_test = Float64.(h5read("test_catvnoncat.h5", "test_set_x"))
    #X_test = reshape(X_test, :, size(X_test, 4))

    y_test = Int64.(h5read("test_catvnoncat.h5", "test_set_y"))
    #y_test = reshape(y_test, (1, size(y_test, 1)))

    categories = h5read("test_catvnoncat.h5", "list_classes")

    return (TrainingInputs = X_train, ValidationInputs = X_test,
            TrainTarget = y_train, ValidationTarget = y_test,
            categories = categories)
end

data = load_data()
reshape(data.TrainingInputs, :, size(data.TrainingInputs, 4))

"""
    Make predictions on new data using the trained paramaters
"""
function predict(DMatrix, parameters)
    Ŷ , _  = forward_propagate_model_weights(DMatrix, parameters)
    return Ŷ
end

#url = "https://raw.githubusercontent.com/PyDataBlog/NN-From-Scratch-With-Julia/master/train_catvnoncat.h5";

# Generate fake data
X, y = make_blobs(10_000, 3; centers=2, as_table=false, rng=2020);
X = Matrix(X');
y = reshape(y, (1, size(X, 2)));
f(x) =  x == 2 ? 0 : x
y2 = f.(y);

# Input dimensions
input_dim = size(X, 1);

nn_results = train_network([input_dim, 5, 3, 1], X, y2; η=0.01, epochs=50, seed=1, verbose=true);

max_stars_repo_name: PyDataBlog/ML-From-Scratch-With-Julia

max_issues_repo_name: PyDataBlog/ML-From-Scratch-With-Julia

max_forks_repo_name: PyDataBlog/ML-From-Scratch-With-Julia