text:
## Preloads
using Statistics
using FFTW
using Plots
using BenchmarkTools
using Profile
using LinearAlgebra
#using StaticArrays
#using HybridArrays

# put the cwd on the path so Julia can find the module
push!(LOAD_PATH, pwd())
using DHC_2DUtils

##

# create filter bank
@time filter_hash = fink_filter_hash(1,8)

# instantiate image of wavelet, centered
function wavelet_image(h, ind)
    npix   = h["npix"]
    pixind = h["filt_index"][ind]
    pixval = h["filt_value"][ind]
    ψ̃ = zeros(npix,npix)
    ψ̃[pixind] = pixval
    ψ = ifft(ψ̃)
    return fftshift(ψ)
end

print(sum(wavelet_image(filter_hash,49)))

test_img = rand(256,256)

temp2 = speedy_DHC(test_img, filter_list)

@time speedy_DHC(test_img, filter_list)
@benchmark speedy_DHC(test_img, filter_list)
# mean time:        194.513 ms (8.08% GC)

function DHC(image::Array{Float64,2}, filter_set::Array{Float64,4}; norm_on = 1, coeff_12_on =1, coeff_20_on = 1)
    FFTW.set_num_threads(2)

    #Sizing
    (Nx, Ny) = size(image)
    (_,_,J,L) = size(fink_filter_set)

    out_coeff = []

    #Preallocate Coeff Arrays
    S0 = zeros(2)
    S1 = zeros(J,L)
    if coeff_20_on == 1
        S20 = zeros(Float64,J, L, J, L)
    end
    if coeff_12_on == 1
        S12 = zeros(Float64,J, L, J, L)
    end

    if coeff_20_on == 1
        im_rd_0_1 = zeros(Float64,Nx, Ny, J, L)
    end
    if coeff_12_on == 1
        im_fdf_0_1 = zeros(Float64,Nx, Ny, J, L)
    end

    im_fd_0_1 = zeros(ComplexF64,Nx, Ny, J, L)

    Atmp = zeros(ComplexF64,Nx,Ny)
    Btmp = zeros(Float64,Nx,Ny)
    Ctmp = zeros(Float64,Nx,Ny)
    Dtmp = zeros(Float64,Nx,Ny)
    Etmp = zeros(Float64,Nx,Ny)

    ## 0th Order
    if norm_on == 1
        @inbounds S0[1] = mean(image)
        @inbounds norm_im = image.-S0[1]
        @inbounds S0[2] = BLAS.dot(norm_im,norm_im)/(Nx*Ny)
        @inbounds norm_im = norm_im./sqrt(Nx*Ny*S0[2])
    else
        norm_im = image
    end

    append!(out_coeff,S0)

    ## 1st Order
    im_fd_0 = fft(norm_im)

    @views for l = 1:L
        for j = 1:J
            @inbounds im_fd_0_1[:,:,j,l] .= im_fd_0
        end
    end

    #Precompute loop for 2nd Order
    @views for l = 1:L
        for j = 1:J
            @inbounds had!(im_fd_0_1[:,:,j,l],filter_set[:,:,j,l]) #wavelet already in fft domain not shifted
            @inbounds Btmp .= abs.(im_fd_0_1[:,:,j,l])
            @inbounds S1[j,l]+=BLAS.dot(Btmp,Btmp) #normalization choice arb to make order unity
            if coeff_20_on == 1
                @inbounds Atmp .= ifft(im_fd_0_1[:,:,j,l])
                @inbounds im_rd_0_1[:,:,j,l] .= abs.(Atmp)
            end
            if coeff_12_on == 1
                @inbounds im_fdf_0_1[:,:,j,l] .= abs.(im_fd_0_1[:,:,j,l])
                @inbounds im_fdf_0_1[:,:,j,l] .= fftshift(im_fdf_0_1[:,:,j,l])
            end
        end
    end
    append!(out_coeff,S1)

    ## 2nd Order
    @views for l2 = 1:L
        for j2 = 1:J
            if coeff_20_on == 1
                @inbounds copyto!(Btmp, im_rd_0_1[:,:,j2,l2])
            end
            if coeff_12_on == 1
                @inbounds copyto!(Dtmp, im_fdf_0_1[:,:,j2,l2])
            end
            for l1 = 1:L
                for j1  = 1:J
                    if coeff_20_on == 1
                        @inbounds copyto!(Ctmp, im_rd_0_1[:,:,j1,l1])
                        @inbounds S20[j1,l1,j2,l2] += BLAS.dot(Btmp,Ctmp)
                    end
                    if coeff_12_on == 1
                        @inbounds copyto!(Etmp, im_fdf_0_1[:,:,j1,l1])
                        @inbounds S12[j1,l1,j2,l2] += BLAS.dot(Dtmp,Etmp)
                    end
                end
            end
        end
    end
    if coeff_20_on == 1
        append!(out_coeff,S20)
    end
    if coeff_12_on == 1
        append!(out_coeff,S12)
    end
    return out_coeff
end

temp = DHC(test_img,fink_filter_set,coeff_12_on =0, coeff_20_on = 1)

@time DHC(test_img,fink_filter_set,coeff_12_on =0, coeff_20_on = 1)

Profile.clear()
@profile DHC(test_img,fink_filter_set,coeff_12_on =0, coeff_20_on = 1)

Juno.profiler()

BenchmarkTools.DEFAULT_PARAMETERS.samples = 5
BenchmarkTools.DEFAULT_PARAMETERS.seconds = 120
@benchmark DHC(test_img,fink_filter_set,coeff_12_on =0, coeff_20_on = 1)

temp = DHC(test_img,fink_filter_set,coeff_12_on =1, coeff_20_on = 1)

@time DHC(test_img,fink_filter_set,coeff_12_on =1, coeff_20_on = 1)

Profile.clear()
@profile DHC(test_img,fink_filter_set,coeff_12_on =1, coeff_20_on = 1)

Juno.profiler()

BenchmarkTools.DEFAULT_PARAMETERS.samples = 5
BenchmarkTools.DEFAULT_PARAMETERS.seconds = 120
@benchmark DHC(test_img,fink_filter_set,coeff_12_on =1, coeff_20_on = 1)

## No if ands or buts about it...

function DHC(image::Array{Float64,2}, filter_set::Array{Float64,4})
    FFTW.set_num_threads(2)

    #Sizing
    (Nx, Ny)  = size(image)
    (_,_,J,L) = size(fink_filter_set)

    out_coeff = []
    #Preallocate Coeff Arrays
    S0  = zeros(2)
    S1  = zeros(J,L)
    S20 = zeros(Float64,J, L, J, L)
    S12 = zeros(Float64,J, L, J, L)
    im_rd_0_1  = zeros(Float64,Nx, Ny, J, L)
    im_fdf_0_1 = zeros(Float64,Nx, Ny, J, L)
    im_fd_0_1  = zeros(ComplexF64,Nx, Ny, J, L)

    Atmp = zeros(ComplexF64,Nx,Ny)
    Btmp = zeros(Float64,Nx,Ny)
    Ctmp = zeros(Float64,Nx,Ny)
    Dtmp = zeros(Float64,Nx,Ny)
    Etmp = zeros(Float64,Nx,Ny)

    ## 0th Order
    @inbounds S0[1] = mean(image)
    @inbounds norm_im = image.-S0[1]
    # make this a vector before computing variance with BLAS.dot
    # (This was a bug!)
    norm_im_vec = reshape(norm_im, Nx*Ny)
    @inbounds S0[2] = BLAS.dot(norm_im_vec,norm_im_vec)/(Nx*Ny)
    @inbounds norm_im = norm_im./sqrt(Nx*Ny*S0[2])

    # Was this intentional to set this back to image?
    norm_im = image

    append!(out_coeff,S0[:])

    ## 1st Order
    im_fd_0 = fft(norm_im)

    @views for l = 1:L
        for j = 1:J
            @inbounds im_fd_0_1[:,:,j,l] .= im_fd_0
        end
    end

    ## Main 1st Order and Precompute 2nd Order
    @views for l = 1:L
        for j = 1:J
            @inbounds had!(im_fd_0_1[:,:,j,l],filter_set[:,:,j,l]) #wavelet already in fft domain not shifted
            @inbounds Btmp .= abs.(im_fd_0_1[:,:,j,l])
            @inbounds im_fdf_0_1[:,:,j,l] .= fftshift(Btmp)
            # I think you are doing a 256x256 matrix multiplication here
            #@inbounds S1[j,l]+=BLAS.dot(Btmp,Btmp) #normalization choice arb to make order unity
            # You mean this:  (and BLAS won't speed this up much)
            S1[j,l] = sum(Btmp.*Btmp)
            @inbounds Atmp .= ifft(im_fd_0_1[:,:,j,l])
            @inbounds im_rd_0_1[:,:,j,l] .= abs.(Atmp)

        end
    end
    append!(out_coeff,S1[:])

    ## 2nd Order
    @views for l2 = 1:L
        for j2 = 1:J
                @inbounds copyto!(Btmp, im_rd_0_1[:,:,j2,l2])
                @inbounds copyto!(Dtmp, im_fdf_0_1[:,:,j2,l2])
            for l1 = 1:L
                for j1  = 1:J
                        @inbounds copyto!(Ctmp, im_rd_0_1[:,:,j1,l1])
                        @inbounds copyto!(Etmp, im_fdf_0_1[:,:,j1,l1])
                        @inbounds S20[j1,l1,j2,l2] += BLAS.dot(Btmp,Ctmp)
                        @inbounds S12[j1,l1,j2,l2] += BLAS.dot(Dtmp,Etmp)
                    end
                end
            end
        end
    append!(out_coeff,S20)
    append!(out_coeff,S12)

    #DPF version of 2nd Order
    Amat    = reshape(im_fdf_0_1, Nx*Ny, J*L)
    FastS12 = reshape(Amat' * Amat, J, L, J, L)
    diffS12 = S12-FastS12
    println("minmax(diffS12)", minimum(diffS12), "   ", maximum(diffS12))

    Amat    = reshape(im_rd_0_1, Nx*Ny, J*L)
    FastS20 = reshape(Amat' * Amat, J, L, J, L)
    diffS20 = S20-FastS20
    println("minmax(diffS20)", minimum(diffS20), "   ", maximum(diffS20))

    return out_coeff
end

temp = DHC(test_img,fink_filter_set)

@time DHC(test_img,fink_filter_set)

Profile.clear()
@profile DHC(test_img,fink_filter_set)

Juno.profiler()

BenchmarkTools.DEFAULT_PARAMETERS.samples = 10000
BenchmarkTools.DEFAULT_PARAMETERS.seconds = 120
@benchmark DHC(test_img,fink_filter_set,coeff_12_on =1, coeff_20_on = 1)

1.682/8258
8/1857

##Seems the order in memory changes performance a bit...
#Let me just check for a faster abs.() by hand and then call this good
#for a bit

function man_abs_dot!(A,B)
    m,n = size(A)
    @assert (m,n) == size(B)
    for j in 1:n
       for i in 1:m
         @inbounds A[i,j] = abs(B[i,j])
       end
    end
    return A
end

function wrapper(A,B)
    A .= abs.(B)
end

@time man_abs_dot!(zeros(256,256),test_img)

@time wrapper(zeros(256,256),test_img)

compare_img = zeros(256,256)

BenchmarkTools.DEFAULT_PARAMETERS.seconds = 10
@benchmark wrapper($compare_img,test_img)

@benchmark man_abs_dot!($compare_img,test_img)

## Nope, not all that much faster... ok. Good night.
function DHC(image::Array{Float64,2}, filter_set::Array{Float64,4})
    FFTW.set_num_threads(2)
    #Sizing
    (Nx, Ny) = size(image)
    (_,_,J,L) = size(fink_filter_set)
    out_coeff = []
    #Preallocate Coeff Arrays
    S0 = zeros(2)
    S1 = zeros(J,L)
    S20 = zeros(Float64,J, L, J, L)
    S12 = zeros(Float64,J, L, J, L)
    im_rd_0_1 = zeros(Float64,Nx, Ny, J, L)
    im_fdf_0_1 = zeros(Float64,Nx, Ny, J, L)
    im_fd_0_1 = zeros(ComplexF64,Nx, Ny, J, L)
    Atmp = zeros(ComplexF64,Nx,Ny)
    Btmp = zeros(Float64,Nx,Ny)
    Ctmp = zeros(Float64,Nx,Ny)
    Dtmp = zeros(Float64,Nx,Ny)
    Etmp = zeros(Float64,Nx,Ny)
    ## 0th Order
    @inbounds S0[1] = mean(image)
    @inbounds norm_im = image.-S0[1]
    @inbounds S0[2] = BLAS.dot(norm_im,norm_im)/(Nx*Ny)
    @inbounds norm_im = norm_im./sqrt(Nx*Ny*S0[2])
    norm_im = image
    append!(out_coeff,S0[:])
    ## 1st Order
    im_fd_0 = fft(norm_im)
    @views for l = 1:L
        for j = 1:J
            @inbounds im_fd_0_1[:,:,j,l] .= im_fd_0
        end
    end
    ## Main 1st Order and Precompute 2nd Order
    @views for l = 1:L
        for j = 1:J
            @inbounds had!(im_fd_0_1[:,:,j,l],filter_set[:,:,j,l]) #wavelet already in fft domain not shifted
            @inbounds Btmp .= abs.(im_fd_0_1[:,:,j,l])
            @inbounds im_fdf_0_1[:,:,j,l] .= fftshift(Btmp)
            @inbounds S1[j,l]+=BLAS.dot(Btmp,Btmp) #normalization choice arb to make order unity
            @inbounds Atmp .= ifft(im_fd_0_1[:,:,j,l])
            @inbounds im_rd_0_1[:,:,j,l] .= abs.(Atmp)
        end
    end
    append!(out_coeff,S1[:])
    ## 2nd Order
    @views for l2 = 1:L
        for j2 = 1:J
                @inbounds copyto!(Btmp, im_rd_0_1[:,:,j2,l2])
                @inbounds copyto!(Dtmp, im_fdf_0_1[:,:,j2,l2])
            for l1 = 1:L
                for j1  = 1:J
                        @inbounds copyto!(Ctmp, im_rd_0_1[:,:,j1,l1])
                        @inbounds copyto!(Etmp, im_fdf_0_1[:,:,j1,l1])
                        @inbounds S20[j1,l1,j2,l2] += BLAS.dot(Btmp,Ctmp)
                        @inbounds S12[j1,l1,j2,l2] += BLAS.dot(Dtmp,Etmp)
                    end
                end
            end
        end
    append!(out_coeff,S20)
    append!(out_coeff,S12)
    return out_coeff
end

## No if ands or buts about it...

# Todo list
# Check if 2 threads really help FFT when computer is busy

#test_img = rand(256,256)
#temp1 = DHC(test_img,fink_filter_set)
temp2 = speedy_DHC(test_img, filter_list)

@time speedy_DHC(test_img, filter_list);

Profile.clear()
@profile speedy_DHC(test_img, filter_list)
Juno.profiler()

function doug_DHC(image::Array{Float64,2}, filter_list)
    # Use 2 threads for FFT
    FFTW.set_num_threads(2)

    # array sizes
    (Nx, Ny)  = size(image)
    (J,L)     = size(filter_list[1])

    # allocate coeff arrays
    out_coeff = []
    S0  = zeros(Float64, 2)
    S1  = zeros(Float64, J, L)
    S20 = zeros(Float64, J, L, J, L)
    S12 = zeros(Float64, J, L, J, L)
    S2  = zeros(Float64, J, L, J, L)

    # allocate image arrays for internal use
    im_fdf_0_1 = zeros(Float64,           Nx, Ny, J, L)   # this must be zeroed!
    im_rd_0_1  = Array{Float64, 4}(undef, Nx, Ny, J, L)

    ## 0th Order
    S0[1]   = mean(image)
    norm_im = image.-S0[1]
    S0[2]   = sum(norm_im .* norm_im)/(Nx*Ny)
    norm_im ./= sqrt(Nx*Ny*S0[2])

    append!(out_coeff,S0[:])

    ## 1st Order
    im_fd_0 = fft(norm_im)

    # unpack filter_list
    f_ind   = filter_list[1]  # (J, L) array of filters represented as index value pairs
    f_val   = filter_list[2]

    zarr = zeros(ComplexF64, Nx, Ny)  # temporary array to fill with zvals

    # make a FFTW "plan" for an array of the given size and type
    P   = plan_ifft(im_fd_0)   # P is an operator, P*im is ifft(im)

    ## Main 1st Order and Precompute 2nd Order
    for l = 1:L
        for j = 1:J
            S1tot = 0.0
            f_i = f_ind[j,l]  # CartesianIndex list for filter
            f_v = f_val[j,l]  # Values for f_i
            # for (ind, val) in zip(f_i, f_v)   # this is slower!
            if length(f_i) > 0
                for i = 1:length(f_i)
                    ind       = f_i[i]
                    zval      = f_v[i] * im_fd_0[ind]
                    S1tot    += abs2(zval)
                    zarr[ind] = zval
                    im_fdf_0_1[ind,j,l] = abs(zval)
                end
                S1[j,l] = S1tot
                im_rd_0_1[:,:,j,l] .= abs2.(P*zarr)  # Inv FFT
                zarr[f_i] .= 0
            end
        end
    end
    append!(out_coeff, S1[:])

    # we stored the abs()^2, so take sqrt (this is faster to do all at once)
    im_rd_0_1 .= sqrt.(im_rd_0_1)

    ## 2nd Order
    Amat = reshape(im_fdf_0_1, Nx*Ny, J*L)
    S12  = reshape(Amat' * Amat, J, L, J, L)
    Amat = reshape(im_rd_0_1, Nx*Ny, J*L)
    S20  = reshape(Amat' * Amat, J, L, J, L)
    append!(out_coeff, S20)
    append!(out_coeff, S12)

    ## Traditional second order
    for l = 1:L
        for j = 1:J
            thisim = fft(im_rd_0_1[:,:,j,l])  # Could try rfft here
            println("j ",j,"  l",l,"  sum(fft):",sum(abs2.(thisim))/65536, "  sum(im): ",sum(abs2.(im_rd_0_1[:,:,j,l])))
            # Loop over (l2,j2) and do second-order convolution
            for l2 = 1:L
                for j2 = 1:J
                    S2tot = 0.0
                    f_i = f_ind[j2,l2]  # CartesianIndex list for filter
                    f_v = f_val[j2,l2]  # Values for f_i

                    S2tot = sum(abs2.(f_v .* thisim[f_i])) # sum im^2 same in real space
                    S2[j,l,j2,l2] = S2tot
                end
            end
        end
    end

    append!(out_coeff, S2)

    return out_coeff
end

Profile.clear()
@profile doug_DHC(test_img, filter_list)
Juno.profiler()

out = doug_DHC(test_rod, filter_list2)
S2 = reshape(out[8259:12354],8,8,8,8)

#---------------

function rod(xcen, ycen, length, pa, fwhm)
    rodimage = zeros(256,256)
    nx = 256
    x=0
    y=0
    sig = fwhm/2.355
    dtor = π/180
    # -------- define a unit vector in direction of rod at position angle pa
    ux = sin(pa*dtor)   #  0 deg is up
    uy = cos(pa*dtor)   # 90 deg to the right

    for i=1:nx
        for j=1:nx
            x=i-nx/2+xcen
            y=j-nx/2+ycen

            # -------- distance parallel and perpendicular to
            dpara =  ux*x + uy*y
            dperp = -uy*x + ux*y

            if abs(dpara)-length <0
                dpara= 0
            end
            dpara = abs(dpara)
            dpara = min(abs(dpara-length),dpara)

            rodimage[i,j] = exp(-(dperp^2+dpara^2)/(2*sig^2))
        end
    end
    rodimage ./= sqrt(sum(rodimage.^2))
    return rodimage
end

@time test_rod = rod(10,10,30,45,4)

function rod_test(filter_hash)

    Nf      = length(filter_hash["filt_value"])
    Nj      = length(filter_hash["j_value"])
    Nl      = length(filter_hash["theta_value"])
    Npix    = filter_hash["npix"]
    jlind   = filter_hash["psi_index"]
    angmax  = 180
    angstep = 2
    nang    = angmax ÷ angstep+1
    energy     = zeros(nang)
    Siso       = zeros(Nj, Nj, Nl, nang)
    Sall       = zeros(Nj, Nj, Nl, Nl, nang)
    Sisoodd    = zeros(Nj, Nj, Nl, nang)
    Sisoeven   = zeros(Nj, Nj, Nl, nang)
    S2all      = zeros(Nf, Nf, nang)
    Ef         = zeros(Nf, nang)  # energy in each filter
    for i = 1:nang
        ang = i*angstep
        rod_image = rod(10,10,30,ang,4)
        # t = speedy_DHC(rod_image, filter_list)
        t   = DHC_compute(rod_image, filter_hash)
        S20_big = reshape(t[2+Nf+1:2+Nf+Nf*Nf], Nf, Nf)
        #S20 = S20_big[1:Nj*Nl, 1:Nj*Nl]
        off2 = 2+Nf+Nf*Nf
        S2 = reshape(t[off2+1:off2+Nf*Nf], Nf, Nf)
        S2all[:,:,i] = S2
        energy[i] = sum(diag(S20_big))
        Ef[:,i]   = diag(S20_big)
        # if filter_hash["pc"]==1 Ef[1:Nj*Nl,i] .*= 2 end
        # Verify sum of diag(S20) = sum(S1)   (= energy)
        println(i,"   ",energy[i],"   ",sum(t[3:2+Nf]))

        # what about j1, j2, Delta L
        #S = reshape(S20,Nj,Nl,Nj,Nl)  #  (j1, l1, j2, l2), symmetric in j and l
        S = S2[jlind,jlind]  # amazing this works...
        println(sum(S))
        for j1 = 1:Nj
            for j2 = 1:Nj
                for l1 = 1:Nl
                    for l2 = 1:Nl
                        DeltaL = mod(l1-l2, Nl)
                        Siso[j1,j2,DeltaL+1,i] += S[j1,l1,j2,l2]
                        # Sall[j1,j2,l1,l2,i] = S[j1,l1,j2,l2]
                        Sall[j1,j2,l1,l2,i] = S[j1,l1,j2,l2]
                        #println(Siso[j1,j2,DeltaL+1,i])
                    end
                end
                for l1 = 1:2:Nl
                    for l2 = 1:2:Nl
                        DeltaL = mod(l1-l2, Nl)
                        Sisoodd[j1,j2,DeltaL+1,i] += S[j1,l1,j2,l2]
                    end
                end
                for l1 = 2:2:Nl
                    for l2 = 2:2:Nl
                        DeltaL = mod(l1-l2, Nl)
                        Sisoeven[j1,j2,DeltaL+1,i] += S[j1,l1,j2,l2]
                    end
                end

            end
        end

    end
    println(std(Siso))
    return (Siso, Ef, Sall, S2all)
end

function plot_energy(Siso, S2, Ef, h)

    # -------- angle array for plots
    nang = size(Siso)[end]
    angstep = 180/(nang-1)
    ang = collect(0:angstep:180)

    # -------- compute total energy
    Etot = sum(Ef,dims=1)[1,:]
    p = plot(ang, Etot, xguide="Angle [deg]", yguide="Total energy", label="S1 energy")

    # -------- compute total energy from S2
    Etot = sum(S2,dims=(1,2))[1,1,:]
    plot!(ang, Etot, label="S2 energy")
    display(p)

    # -------- energy in the ϕ component
    Eϕ = sum(S2[:,end,:],dims=1)[1,:]
    p=plot(ang, Eϕ, xguide="Angle [deg]", yguide="ϕ energy")
    display(p)

    jlind = h["psi_index"]
    S2_5d = S2[jlind, jlind,:]
    pow   = sum(S2_5d, dims=(2,4))[:,1,:,1,:]
    p=plot(reshape(pow[:,:,:],36,nang)',yaxis=(:log),
           ylims=(1e-6,1e-1), xguide="Rod angle",yguide="Ang averged energy")
    display(p)

    return
    val = Siso[j1,j2,1,:]
    p = plot(ang, val/mean(val))

    for Δl = 2:8
        val = Siso[j1,j2,Δl,:]
        plot!(p, ang, val/mean(val))
    end
    display(p)
    return
end

function plot8(Siso, j1, j2)

    nang = size(Siso)[end]
    angstep = 180/(nang-1)
    ang = collect(0:angstep:180)
    val = Siso[j1,j2,1,:]
    p = plot(ang, val, xguide="Angle", label="Δl=0")
    for Δl = 2:8
        val = Siso[j1,j2,Δl,:]
        plot!(p, ang, val, label="Δl=$(Δl-1)")
    end
    display(p)

    return
end

filter_hash = fink_filter_hash(1,8,wd=2)
(Siso, Ef, Sall, S2) = rod_test(filter_hash)
plot_energy(Siso,S2,Ef,filter_hash)

heatmap(Siso[:,:,1,14])
plot(Siso[4,4,1,:]/Siso[4,4,1,1])

plot8(Siso,3,3)

max_stars_repo_name: andrew-saydjari/DHC

max_issues_repo_name: andrew-saydjari/DHC

max_forks_repo_name: andrew-saydjari/DHC


text:
using LinearAlgebra
using GLM
using LsqFit
using Distributions
using StatsBase
using Printf
using Plots

#linear fit wrapper

function lfit(ndf::DataFrame)
    lr = lm(@formula(y_mean ~ x_mean), ndf)
    c = coef(lr)
    return x -> c[1] + c[2]*x, predict(lr), c
end

struct RFit
	fitpar::Vector{<:Number}
	fitstd::Vector{<:Number}
	ci::Vector{Tuple{Number, Number}}
	g::Function
end

"""
    gpol1(ct::Vector{<:Real})::Function
return a degree 1 polynomial function with parameters ct.
"""
function gpol1(ct::Vector{<:Real})::Function
    function f1(z::Real)
        return ct[1] + ct[2] * z
    end
    return f1
end

"""
    gpol2(ct::Vector{<:Real})::Function
return a degree 2 polynomial function with parameters ct.
"""
function gpol2(ct::Vector{<:Real})::Function
    function f2(z::Real)
        return ct[1] + ct[2] * z + ct[3] * z^2
    end
    return f2
end

"""
    gpol2(ct::Vector{<:Real})::Function
return a degree 3 polynomial function with parameters ct.
"""
function gpol3(ct::Vector{<:Real})::Function
    function f3(z::Real)
        return ct[1] + ct[2] * z + ct[3] * z^2 + ct[4] * z^3
    end
    return f3
end

"""
    func1dfit(ffit::Function, x::Vector{<:Real},
              y::Vector{<:Real}, p0::Vector{<:Real},
              lb::Vector{Float64}, ub::Vector{Float64})

Fit a function to the data x, y with start prediction p0
and return coefficients and errors.
"""
function func1dfit(ffit::Function, x::Vector{<:Real},
                   y::Vector{<:Real}, p0::Vector{<:Real},
                   lb::Vector{Float64}, ub::Vector{Float64})
    fq = curve_fit(ffit, x, y, p0, lower=lb, upper=ub)
    cfq = coef(fq)
    @info "coef(fq)" cfq
    sfq = stderror(fq)
    @info "std(fq)" sfq
    @info "margin_of_error (90%)" margin_error(fq, 0.1)
    @info " confidence_interval (90%)" confidence_interval(fq, 0.1)
    return cfq, sfq
end

"""
    func1dfit(ffit::Function, x::Vector{<:Real}, y::Vector{<:Real},
              yerr::Vector{<:Real}, p0::Vector{<:Real},
              lb::Vector{<:Real}, ub::Vector{<:Real})
Fit function ffit to data x,y taking into account the errors
on the y (yerr). Assumes weights of 1/sigma^2 as in standard least squares
and returns fit result.
"""
function func1dfit(ffit::Function, x::Vector{<:Real}, y::Vector{<:Real},
	               yerr::Vector{<:Real}, p0::Vector{<:Real},
                   lb::Vector{<:Real}, ub::Vector{<:Real})
	fq = curve_fit(ffit, x, y, yerr.^-2, p0, lower=lb, upper=ub)
    return fq
end

"""
    fit_pol1(x::Vector{<:Real}, y::Vector{<:Real}, ci=0.1)
Fit a 1st degree polynomial to the data x, y and return
a RFit object with the result.
"""
function fit_pol1(x::Vector{<:Real}, y::Vector{<:Real}, ci=0.1)
    @. pol(x, p) = p[1] + p[2] * x
    p0 = [1.0, 1.0]
    fq = curve_fit(pol, x, y, p0)

    RFit(coef(fq), stderror(fq), confidence_interval(fq, ci),
         gpol1(coef(fq)))
end

"""
    fit_pol2(x::Vector{<:Real}, y::Vector{<:Real}, ci=0.1)
Fit a 2nd degree polynomial to the data x, y and return
a RFit object with the result.
"""
function fit_pol2(x::Vector{<:Real}, y::Vector{<:Real}, ci=0.1)
    @. pol(x, p) = p[1] + p[2] * x + p[3] * x^2
    p0 = [1.0, 1.0, 1.0]
    fq = curve_fit(pol, x, y, p0)
    RFit(coef(fq), stderror(fq), confidence_interval(fq, ci),
         gpol2(coef(fq)))
end

"""
    fit_pol3(x::Vector{<:Real}, y::Vector{<:Real}, ci=0.1)
Fit a 3rd degree polynomial to the data x, y and return
a RFit object with the result.
"""
function fit_pol3(x::Vector{<:Real}, y::Vector{<:Real}, ci=0.1)
    @. pol(x, p) = p[1] + p[2] * x + p[3] * x^2 + p[4] * x^3
    p0 = [1.0, 1.0, 1.0, 1.0]
    fq = curve_fit(pol, x, y, p0)
    RFit(coef(fq), stderror(fq), confidence_interval(fq, ci),
         gpol3(coef(fq)))
end

struct FGauss
	mu ::Vector{<:Real}
	std::Vector{<:Real}
	C  ::Vector{<:Real}
	h  ::Histogram
	X  ::Vector{<:Real}
	Y  ::Vector{<:Real}
	g  ::Vector{Function}
end

"""
    gausg(μ::Real, σ::Real, C::Real)::Function
Return a Gaussian function with parameters μ, σ and normalisation C.
"""
function gausg(μ::Real, σ::Real, C::Real)::Function
	function gausx(x::Real)
		return C * pdf(Normal(μ, σ,), x)
	end
	return gausx
end

"""
    gausg2(μ1::Real, σ1::Real, C1::Real, μ2::Real, σ2::Real, C2::Real)
Return a function for the sum of two Gaussians with parameters
    μ1, σ1 and μ2, σ2 and normalisations C1 and C2.
"""
function gausg2(μ1::Real, σ1::Real, C1::Real, μ2::Real, σ2::Real, C2::Real)::Function
	function gausx(x::Real)
		return C1 * pdf(Normal(μ1, σ1,), x) + C2 * pdf(Normal(μ2, σ2,), x)
	end
	return gausx
end

@. gauss1(x, p)   = p[1] * pdf(Normal(p[2], p[3]), x)
@. gauss2(x, p)   = p[1] * pdf(Normal(p[2], p[3]), x) + p[4] * pdf(Normal(p[5], p[6]), x)
@. gauss2cm(x, p) = p[1] * pdf(Normal(p[2], p[3]), x) + p[4] * pdf(Normal(p[2], p[5]), x)

"""
    gaussfm(mu::Real)::Function
Return a Gaussian Fitting function with a fixed mean mu.
"""
function gaussfm(mu::Real)::Function
	function gauss(x::Vector{<:Real}, p::Vector{<:Real})
		return @. p[1]* pdf(Normal(mu, p[2]), x)
	end
	return gauss
end

"""
    gauss2fm(mu::Real)::Function
Return a function of the sum of two Gaussians with a fixed common mean mu.
"""
function gauss2fm(mu::Real)::Function
	function gauss2(x::Vector{<:Real}, p::Vector{<:Real})
		return @. p[1]* pdf(Normal(mu, p[2]), x) + p[3]* pdf(Normal(mu, p[4]), x)
	end
	return gauss2
end

"""
	fit_gauss(h::Histogram)

Fit a normal distribution to a Histogram and return FGauss object.
"""
function fit_gauss(h::Histogram)
	c = centers(h)
	w = eltype(c).(h.weights)
	@debug "histo"  w c
	mu, sigma = mean_and_std(c, Weights(h.weights); corrected = false)
	@debug "mu, std" mu, sigma

	# fit parameters lb, ub, po are lower, upper bounds and pars
    lb = [         0.0, mu - 100.0 * sigma, sigma / 100.0]
    ub = [100 * sum(w), mu + 100.0 * sigma, 100.0 * sigma]
    p0 = [      sum(w), mu                ,         sigma]
	(CC, μ, σ), _ = func1dfit(gauss1, c, w, p0, lb, ub)
	gx = gausg(μ, σ, CC)
	return FGauss([μ], [σ], [CC], h, c, gx.(c), [gx])

end

"""
	fit_gauss(x::Vector{Float64}, xmin::Float64, xmax::Float64;
              bins::Integer=50, norm=false)

Fit a normal distribution to data x,y and return FGauss object.
"""
fit_gauss(x::Vector{Float64}, xmin::Float64, xmax::Float64;
	      bins::Integer=50, norm=false) = fit_gauss(hist1d(x, bins, xmin, xmax, norm))

"""
	fit_gauss_fm(x::Vector{Float64}, xmin::Float64, xmax::Float64, bins=50, fm=0.0)

Fit a gaussian with a fixed mean fm and return an FGauss object.
"""
function fit_gauss_fm(x::Vector{Float64}, xmin::Float64, xmax::Float64;
	                  bins=50, norm=false, fm=0.0)

	# fit the unbinned distribution
	xx =  in_range(x, xmin, xmax)
	σ  = std(xx)
	@debug "gfit_gauss_fm: σ = $σ"

	# bin distribution
    h = hist1d(xx, bins, xmin, xmax, norm)
    c = centers(h)
    w = eltype(c).(h.weights)
    @debug "histo w and c"  w c

	# fit parameters lb, ub, po are lower, upper bounds and pars
    lb = [            0.0, σ / 10.0]
    ub = [1.0e+6 * sum(w), σ * 10.0]
    p0 = [         sum(w), σ       ]

	g1 = gaussfm(fm)
	(CC, sigma), _ = func1dfit(g1, c, w, p0, lb, ub)

	@debug "CC,  sigma"  CC  sigma
	gx = gausg(fm, sigma, CC)

	return FGauss([fm], [sigma], [CC], h, c, gx.(c), [gx])
end

"""
    plot_fit_gauss(x::Vector{Float64}, xs::String, ys::String,
                   bins::Integer, xmin::Float64, xmax::Float64;
                   xgmin::Float64, xgmax::Float64, gbins::Integer=50)
Fit a Gaussian to the histogram of data x between xmin and xmax with bins
and plot the results.
"""
function plot_fit_gauss(x::Vector{Float64}, xs::String, ys::String,
                        bins::Integer, xmin::Float64, xmax::Float64;
                        xgmin::Float64, xgmax::Float64, gbins::Integer=50)

    h, _ = hist1d(x, xs, bins, xmin, xmax, norm=true)
    fg   = fit_gauss(x, xgmin, xgmax, bins=gbins, norm=true)
    gx   = fg.g[1]
    X    = centers(h)
    Y    = h.weights
    σY   = sqrt.(Y)
    lbl  = @sprintf " μ=%5.1f, σ =%5.1f " fg.mu[1] fg.std[1]
    lbl  = string("gaussian fit:\n", lbl)
    p    = scatter(X, Y, yerror=σY,fmt=:png, shape=:circle, color=:black, label="data", legend=true)
    p    = plot!(p, X, gx.(X), lw=2, label=lbl, legend=true, fmt = :png)
    xlabel!(xs)
    xlabel!(ys)
    return fg, p
end

"""
    fitg1(x::Vector{<:Real}, xs::String, nbins::Int64, xmin::Real, xmax::Real;
	      xgmin::Real, xgmax::Real, fbins::Int64=100, norm::Bool=true,
          fm::Real=0.0, flex_mean=false)

Histograms data x and returns a plot and a Gaussian fit to the histogram.
Option flex_mean allows the mean to float or be fixed to fm.
"""
function fitg1(x::Vector{<:Real}, xs::String, nbins::Int64, xmin::Real, xmax::Real;
	           xgmin::Real, xgmax::Real, fbins::Int64=100, norm::Bool=true,
               fm::Real=0.0, flex_mean=false)

	h, p = hist1d(x, xs, nbins, xmin, xmax, norm=norm, legend=true)
    if flex_mean
        fg = fit_gauss(x, xgmin, xgmax, bins=fbins, norm=norm)
    else
        fg = fit_gauss_fm(x, xgmin, xgmax, bins=fbins, norm=norm, fm=fm)
    end
	gx  = fg.g[1]
	X   = centers(h)
	lbl = @sprintf "σ =%4.1f " fg.std[1]
    p   = plot!(p, X, gx.(X), lw=2, label=lbl, legend=true, fmt = :png)
	xlabel!(xs)
	return fg, p
end

"""
	gfit_gauss2_cmean(y, xmin, xmax, bins, sigmas, cs, cmean=0.0)

Fit a double gaussian (with sigmas -->[sigma1, sigma2] cs -->[c1, c2] )
and a fixed mean, cmean, to data.
"""
function gfit_gauss2_cmean(y::Vector{Float64}, xmin::Float64, xmax::Float64,
	                       bins::Integer, sigmas::Vector{Float64}, cs::Vector{Float64},
						   norm=false, cmean=0.0)

    x = in_range(y, xmin, xmax)
    h = hist1d(x, bins, xmin, xmax, norm)
    c = centers(h)
    w = h.weights
    @debug "histo centers and weights in full region"  w c

    g2 = gauss2fm(cmean)
    # fit parameters lb, ub, po are lower, upper bounds and pars

    lb = [cs[1]/100.0, sigmas[1]/5.0, cs[2]/100.0, sigmas[2]/5.0]
    ub = [cs[1]*100.0, sigmas[1]*5.0, cs[2]*100.0, sigmas[2]*5.0]
    p0_bounds = [cs[1], sigmas[1], cs[2], sigmas[2]]

    @debug "pars" p0 lb ub
    # fit double gaussian
    fq = curve_fit(g2, c, w, p0_bounds, lower=lb, upper=ub)
    C1, sigma1, C2,  sigma2   = coef(fq)
    @debug "C1 sigma1 C2 sigma2" C1 sigma1 C2  sigma2

    #
    gx  = gausg2(cmean, sigma1, C1, cmean, sigma2, C2)
    gx1 = gausg( cmean, sigma1, C1)
    gx2 = gausg( cmean, sigma2, C2)

	return FGauss([cmean, cmean], [sigma1, sigma2], [C1, C2],
			       h, c, gx.(c), [gx, gx1, gx2])
    end

"""
	fit_2gauss_cmean(data, gp, g1p, g2p, cm)

Fit two gaussian with common mean JJ Review!!!
"""
function fit_2gauss_cmean(data, gp, g1p, g2p, cm, norm=false)
    gf1 = fit_gauss_fm(data, g1p.xmin,g1p.xmax, bins=g1p.nbin, norm=norm, fm=cm)
    @debug gf1
    gf2 = fit_gauss_fm(data, g2p.xmin,g2p.xmax,bins=g2p.nbin, norm= norm, fm=cm)
    @debug gf2
    gf = gfit_gauss2_cmean(data, gp.xmin,gp.xmax,gp.nbin,
	                       [gf1.std[1], gf2.std[1]], [gf1.C[1], gf2.C[1]], norm)
    @debug gf
    return gf
end

"""
    fitg2(x, xs, xmin, xmax, xg1min, xg1max, xg2min, xg2max, xgmin, xgmax; bins=100)
Fits 2 gaussians with common mean (0 by default) to vector x. JJ Review!!
"""
function fitg2(x::Vector{<:Real}, xs::String, bins::Int64, xmin::Real, xmax::Real;
	           xg1min::Real, xg1max::Real, xg2min::Real, xg2max::Real, xgmin::Real, xgmax::Real, cm=0.0,
         	   g1bins::Int64=100, g2bins::Int64=100, gbins::Int64=100, norm::Bool=true)

	_, p = hist1d(x, xs, bins, xmin, xmax, norm=norm, legend=true)

    g1p = (xmin = xg1min, xmax = xg1max, nbin=g1bins)
    g2p = (xmin = xg2min, xmax = xg2max, nbin=g2bins)
    gp  = (xmin = xgmin , xmax = xgmax , nbin=gbins )

    fg  = fit_2gauss_cmean(x, gp, g1p, g2p, cm, norm)
	gx  = fg.g[1]
	gx1 = fg.g[2]
	gx2 = fg.g[3]

	lbl = @sprintf "σt =%4.1f mm, σ =%4.1f mm" fg.std[1] fg.std[2]
	st  = @sprintf "σt =%4.1f mm " fg.std[1]
	sf  = @sprintf "σ  =%4.1f mm" fg.std[2]
    p = plot!(p, fg.X, fg.Y, label=lbl, lw=2, fmt = :png)
    p = plot!(p, fg.X, gx1.(fg.X), label=st, lw=1, fmt = :png)
    p = plot!(p, fg.X, gx2.(fg.X), label=sf, lw=1, fmt = :png)
    #xlabel!(xs)
    return fg, p
end

"""
	fit_profile(x1, x2, tx1, ty1, fit="pol1", bins=25)
    fit_profile(df1, c1, c2, tx1, ty1, fit="pol1", bins=25)
Create and fit a profile with pol1 or poli2 functions.
Return fit parameters, fit function and plot
"""
fit_profile(df1::DataFrame, c1::String, c2::String,
            tx1::String, ty1::String, fit="pol1", bins=25;
            ybin_width::Real=0.1, ymin::Real=352.4, ymax::Real=392.4,
            min_proportion::Real=0.0) =
	fit_profile(df1[!,c1], df1[!,c2], tx1, ty1, fit, bins,
                ybin_width=ybin_width, ymin=ymin, ymax=ymax, min_proportion=min_proportion)

function fit_profile(x1::Vector{<:Real}, x2::Vector{<:Real},
	                 tx1::String, ty1::String, fit="pol1", bins=25;
                     ybin_width::Real=0.1, ymin::Real=minimum(x2),
                     ymax::Real=maximum(x2), min_proportion::Real=0.0)

    pdf1, _ = p1df(x1,x2, bins, ybin_width=ybin_width, ymin=ymin, ymax=ymax, min_proportion=min_proportion)

    if fit == "pol1"
        fr = fit_pol1(pdf1.x_mean, pdf1.y_mean)
    elseif fit == "pol2"
        fr = fit_pol2(pdf1.x_mean, pdf1.y_mean)
	elseif fit == "pol3"
		fr = fit_pol3(pdf1.x_mean, pdf1.y_mean)
	else
		println("option not implemented")
		return nothing
    end

    p1 = scatter(pdf1.x_mean, pdf1.y_mean, yerror=pdf1.y_std, fmt=:png,
	             shape=:circle, color=:black, legend=false)
    p1 = plot!(p1, pdf1.x_mean, fr.g.(pdf1.x_mean), fmt=:png)
    xlabel!(tx1)
    ylabel!(ty1)

   return fr, p1
end

ext: jl

max_stars_repo_name: andLaing/ATools

max_issues_repo_name: andLaing/ATools

max_forks_repo_name: andLaing/ATools


text:
function pinky10_overlap()
    dir = joinpath(homedir(), "seungmount/Omni/TracerTasks/pinky/evaluation/170123_pinky10_golden")
    base_fn = "chunk_19585-21632_22657-24704_4003-4258"
    gt_seg_fn = joinpath(dir, string(base_fn, "_proofread_compressed.h5"))
    orig_seg_fn = joinpath(dir, string(base_fn, "_original_compressed.h5"))

    valid_fn = joinpath(dir, string(base_fn, ".omni.segments.csv"))
    semantic_fn = joinpath(dir, string("semantic_", base_fn, ".csv"))
    semantic_list = load_semantic_list(semantic_fn)

    gt = subsample3d(load_segmentation(gt_seg_fn), (8,8,1))
    o = subsample3d(load_segmentation(orig_seg_fn), (8,8,1))

    gt_to_o = create_segID_map(gt, o, semantic_list)
    merges = count_splits(gt_to_o)
    semantic_list = append_splits(semantic_list, merges)

    overlap_list = unique(vcat(values(gt_to_o)...))

    o_to_gt = create_segID_map(o, gt, overlap_list)
    splits = count_splits(o_to_gt)
    semantic_list = append_splits(semantic_list, splits)

    seg_list = copy(semantic_list)
    den_splits = seg_list[seg_list[:,2] .== 3, 3] - 1
    plt[:hist](den_splits,bins=27)
    println("No of dendrite splits: $(sum(den_splits))")
end

function pinky10_spine_estimate()

    cd("seungmount/research/Nick/IARPA_draft/recon_stats/first_pass")

    trunks = map(Int,readdlm("checked_trunks"))
    uncorr_es = load_edges("uncorr_cons_edges.csv")
    corr_es = load_edges("corrected_cons_edges.csv")
    impl_om = readdlm("segment_overlap.csv") #NOTE: this doesn't work yet - will change

    om = om_from_impl( impl_om )
    tm = find_trunk_mapping( om, trunks )

    println("Synapse Recovery")
    println(synapse_recovery(uncorr_es, corr_es, tm))

    cd()
end

function s1_graph_clustering()
    dt = Dates.format(Dates.today(), "yymmdd")
    fn = "20000_random"
    dir = joinpath(homedir(), "seungmount/research/tommy/s1/")
    edges_fn = joinpath(dir, "combined_edges.csv")
    edges = load_edges(edges_fn)
    seg_size_fn = joinpath(dir, "seg_size.csv")
    seg_size = load_seg_sizes(seg_size_fn)
    # filtered_seg_ids = filter_seg_ids_by_id(seg_size, 1000)
    # filtered_seg_ids = filter_seg_ids_by_size(seg_size, 5000)
    filtered_seg_ids = filter_seg_ids_randomly(seg_size, 20000, 5000)
    filtered_fn = joinpath(dir, "$(dt)_filtered_seg_ids_$(fn).csv")
    writedlm(filtered_fn, filtered_seg_ids)
    # d = create_graph_dicts(edges, seg_size[:,1])
    d = create_graph_dicts(edges, filtered_seg_ids)
    seg_to_index, index_to_seg, syn_to_segs, segs_to_syn, syn_coords, syn_size = d
    adj = create_adjacency_matrix(seg_to_index, syn_to_segs, syn_size)
    reweight_adj!(adj) # apply log to the synapse weights
    rand_order = randperm(size(adj,1))
    # view_order(adj, rand_order)
    radj = adj[rand_order,rand_order]
    # svd_order = svd_clustering(adj, 16)
    # view_order(radj, svd_order)

    radj_fn = joinpath(dir, "$(dt)_rand_adj_sparse_$(fn).csv")
    write_sparse(radj_fn, radj)

    # RUN ./GenLouvain/graph_clustering.m in MATLAB
    # (can't figure out how to load MEX files correctly in Julia)
    cluster, louvain_order = louvain_clustering(radj)
    # println("Louvain Q1: $Q1")
    view_order(radj, louvain_order)

    # louvain_fn = joinpath(dir, "$(dt)_louvain_order.csv")
    # radj = read_sparse(radj_fn)
    # louvain_fn = joinpath(dir, "$(dt)_louvain_order_large.csv")
    # louvain_order = readdlm(louvain_fn, Int64)[:]
    # view_order(radj, louvain_order)

    perm_fn = joinpath(dir, "$(dt)_cluster_perm_$(fn).csv")
    perm_tbl = hcat(collect(1:size(radj,1)), rand_order, louvain_order, round(Int64, cluster))
    writedlm(perm_fn, perm_tbl)

    cadj = radj[louvain_order, louvain_order]
    cadj_fn = joinpath(dir, "$(dt)_clustered_adj_sparse_$(fn).csv")
    write_sparse(cadj_fn, cadj)

    pre_to_post, post_to_pre = create_post_pre_dicts(syn_to_segs)
    post_to_index, index_to_post = create_index_dict(sort(keys(post_to_pre)))
    dadj = create_dendrite_adjacency_matrix(pre_to_post, post_to_pre, post_to_index)
    dadj_fn = joinpath(dir, "$(dt)_dendrite_adj_sparse_$(fn).csv")
    write_sparse(dadj_fn, dadj)
    dadj_louvain_order, Q1 = louvain_clustering(dadj)
    println("Louvain Q1: $Q1")

    cdadj = dadj[dadj_louvain_order, dadj_louvain_order]
    cdadj_fn = joinpath(dir, "$(dt)_clustered_dendrite_adj_sparse_$(fn).csv")
    write_sparse(cdadj_fn, cdadj)

    perm = rand_order[louvain_order]
    write_cluster(dir, index_to_seg, cadj, perm, 12050:12390, 12050:12390)
end


text:
module UtilityFuncs

using ..PerFlagFuncs

export Flag
export FlagSet
export substring_decomp_by_index
export compute_level_total
export compute_level_per_flag
export get_remaining_flags

"""
```
Flag(word::String,
     word_boundaries_left::Vector{String},
     word_boundaries_right::Vector{String})
```
A flag that BetweenFlags looks for to denote
the start/stop position of a given "level" or
scope. The word boundaries need only be unique
since every permutation of left and right
word boundaries are taken to determine levels.
```
julia>
using BetweenFlags
# find: ["\\nfunction", " function", ";function"]
start_flag = BetweenFlags.Flag("function",
                               ["\\n", "\\s", ";"],
                               ["\\n", "\\s"])
# find: ["\\nend", " end", ";end"]
stop_flag = BetweenFlags.Flag("end",
                              ["\\n", "\\s", ";"],
                              ["\\n", "\\s", ";"])
```
"""
struct Flag
  word :: String
  word_boundaries_left :: Vector{String}
  word_boundaries_right :: Vector{String}
  trigger :: Vector{String}
  function Flag(word, word_boundaries_left, word_boundaries_right)
    trigger = Vector{String}()
    for left in word_boundaries_left
      for right in word_boundaries_right
        push!(trigger, string(left, word, right))
      end
    end
    return new(word, word_boundaries_left, word_boundaries_right, trigger)
  end
end

"""
```
FlagSet(start::Flag, stop::Flag)
```
A flag set that defines the start and stop of
the substring of interest.

```
julia>
using BetweenFlags
# find: ["\\nfunction", " function", ";function"]
start_flag = BetweenFlags.Flag("function",
                               ["\\n", "\\s", ";"],
                               ["\\n", "\\s"])
# find: ["\\nend", " end", ";end"]
stop_flag = BetweenFlags.Flag("end",
                              ["\\n", "\\s", ";"],
                              ["\\n", "\\s", ";"])
flag_set = FlagSet(start_flag, stop_flag)
```
"""
struct FlagSet
  start :: Flag
  stop  :: Flag
  ID  :: String
  function FlagSet(start::Flag, stop::Flag)
    # return new(start, stop, start.word*stop.word) # needs testing
    return new(start, stop, start.word)
  end
end

function substring_decomp_by_index(s::String,
                                   i_start::Int,
                                   i_end::Int,
                                   flags_start::Vector{String},
                                   flags_stop::Vector{String},
                                   inclusive::Bool = true)
  if inclusive # middle and after depend on length of stop flag (LSTOP)...
    LSTOP = [length(x) for x in flags_stop if occursin(x, s[i_end:end])][1]
    before = s[1:i_start-1]
    middle = s[i_start:i_end+LSTOP-1]
    after = s[i_end+LSTOP:end]
  else # before and middle depend on length of start flag (LSTART)...
    LSTART = [length(x) for x in flags_start if occursin(x, s[1:i_start+length(x)])][1]
    before = s[1:i_start+LSTART-1]
    middle = s[i_start+LSTART:i_end-1]
    after = s[i_end:end]
  end
  return before, middle, after
end

function compute_level_total(s::String,
                             flags_start::Vector{String},
                             flags_stop::Vector{String})
  L_s = length(s)
  level = zeros(Int, L_s)
  D_o = Dict(x => length(x) for x in flags_start)
  D_c = Dict(x => length(x) for x in flags_stop)
  L_c_arr = [v for (k, v) in D_o]
  L_o_arr = [v for (k, v) in D_c]
  L_delim_max = max([max(L_o_arr...), max(L_c_arr...)]...)
  L_delim_min = min([min(L_o_arr...), min(L_c_arr...)]...)
  if L_s>L_delim_min
    for i in 1:L_s-L_delim_max
      if any([s[i:i+D_o[x]-1] == x for x in flags_start])
        level[i:end] .= level[i:end] .+ 1
      elseif any([s[i:i+D_c[x]-1] == x for x in flags_stop])
        level[i+1:end] .= level[i+1:end] .- 1
      end
    end
  end
  return level
end

function compute_level_per_flag(s::String,
                                level_total::Vector{Int},
                                flag_set_all::Vector{FlagSet})
  # Algorithm:
  # 1) Initialize level_total_modified = level_total_total
  # 2) Find the maximum of level_total_modified and
  #    ask "Which key is responsible for the most recent
  #    increase in level_total_modified?" Assign the corresponding
  #    indexes in the dictionary.
  # 3) Set level_total_modified = level_total_modified-1
  #    where these dictionary indexes were set (since
  #    the solution in these locations are now known).
  # 4) Repeat 2-3 until level_total_modified = 0 everywhere
  level_total_modified = copy(level_total)
  N = length(level_total)
  D = Dict(FS.ID => zeros(Int, N) for FS in flag_set_all)
  max_lev = max(level_total...)
  for i in max_lev:-1:1
    L_max = max(level_total_modified...)
    i_maxes = [i for (i, x) in enumerate(level_total_modified) if x==L_max]
    L_maxes = split_by_consecutives(i_maxes)
    for i_max in L_maxes
      for FS in flag_set_all
        cond_any = any([x==s[i_max[1]:i_max[1]+length(x)-1] for x in FS.start.trigger])
        if cond_any
          D[FS.ID][i_max] .+= 1
        end
      end
      level_total_modified[i_max].-= 1
    end
  end
  temp = [D[FS.ID] for FS in flag_set_all]
  level_total_check = sum(temp, dims=1)[1]
  err = [abs(x-y) for (x, y) in zip(level_total_check, level_total)]
  if !all([x<0.01 for x in err])
    error("Error: levels not conservative.")
  end
  return D
end

function get_remaining_flags(s::String,
                             flags_start::Vector{String},
                             flags_stop::Vector{String})::Bool
  same_flags = all([x==y for x in flags_start for y in flags_stop])
  if same_flags
    remaining_flags = any([occursin(y, s) for y in flags_start]) && any([occursin(y, s) for y in flags_stop])
  else
    c_start = sum([count_flags(s, y) for y in flags_start])
    c_stop  = sum([count_flags(s, y) for y in flags_stop])
    if c_start == 0 || c_stop == 0
      remaining_flags = false
    else
      f_start = [findfirst(y, s)[1] for y in flags_start]
      f_stop  = [findfirst(y, s)[1] for y in flags_stop]
      remaining_flags = any([a<b for a in f_start for b in f_stop])
    end
  end
  return remaining_flags
end

end


text:
struct BoxplotGeometry <: Gadfly.GeometryElement
    default_statistic::Gadfly.StatisticElement
    suppress_outliers::Bool
    tag::Symbol
end

BoxplotGeometry(; method=:tukey, suppress_outliers=false, tag=empty_tag) =
    BoxplotGeometry(Gadfly.Stat.boxplot(method=method), suppress_outliers, tag)

"""
    Geom.boxplot[(; method=:tukey, suppress_outliers=false)]

Draw box plots of the `middle`, `lower_hinge`, `upper_hinge`, `lower_fence`,
`upper_fence`, and `outliers` aesthetics.  The categorical `x` aesthetic is
optional.  If `suppress_outliers` is true, don't draw points indicating
outliers.

Alternatively, if the `y` aesthetic is specified instead, the middle, hinges,
fences, and outliers aesthetics will be computed using [`Stat.boxplot`](@ref).

Boxplots will be automatically dodged by specifying a `color` aesthetic different to the `x` aesthetic.
"""
const boxplot = BoxplotGeometry

element_aesthetics(::BoxplotGeometry) = [:x, :color,
                                         :middle,
                                         :upper_fence, :lower_fence,
                                         :upper_hinge, :lower_hinge]  ### and :y, :outliers

default_statistic(geom::BoxplotGeometry) = geom.default_statistic

function render(geom::BoxplotGeometry, theme::Gadfly.Theme, aes::Gadfly.Aesthetics)
    Gadfly.assert_aesthetics_defined("Geom.bar", aes,
                                     :lower_fence, :lower_hinge,
                                     :upper_hinge, :upper_fence,)
    Gadfly.assert_aesthetics_equal_length("Geom.bar", aes,
                                     :lower_fence, :lower_hinge, :middle,
                                     :upper_hinge, :upper_fence, :outliers)

    default_aes = Gadfly.Aesthetics()
    default_aes.color = discretize_make_ia(RGB{Float32}[theme.default_color])
    default_aes.x = Float64[0.5]
    aes = inherit(aes, default_aes)

    n = length(aes.lower_hinge)

    bw = 1w / n - theme.boxplot_spacing # boxplot width
    if length(aes.x) > 1
        xs_sorted = sort(aes.x)
        minspan = minimum([xj - xi for (xi, xj) in zip(xs_sorted[1:end-1],
                                                       xs_sorted[2:end])])
        bw = minspan*cx - theme.boxplot_spacing
    end

    fw = 2/3 * bw # fence width
    xs = Measure[x*cx for x in takestrict(cycle(aes.x), n)]
    cs = takestrict(cycle(aes.color), n)

    # We allow lower_hinge > upper_hinge, and lower_fence > upper_fence. So we
    # need to organize them for drawing here
    lower_hinge = similar(aes.lower_hinge)
    upper_hinge = similar(aes.upper_hinge)
    lower_fence = similar(aes.lower_fence)
    upper_fence = similar(aes.upper_fence)
    for (i, (lh, uh, lf, uf)) in enumerate(zip(aes.lower_hinge, aes.upper_hinge,
                                               aes.lower_fence, aes.upper_fence))
        if uh > lh
            lower_hinge[i] = aes.lower_hinge[i]
            upper_hinge[i] = aes.upper_hinge[i]
        else
            lower_hinge[i] = aes.upper_hinge[i]
            upper_hinge[i] = aes.lower_hinge[i]
        end

        if uf > lf
            lower_fence[i] = aes.lower_fence[i]
            upper_fence[i] = aes.upper_fence[i]
        else
            lower_fence[i] = aes.upper_fence[i]
            upper_fence[i] = aes.lower_fence[i]
        end
    end

    tbox, tlw, tuw, tlf, tuf, to, tm =
        subtags(geom.tag, :box, :lower_whisker, :upper_whisker,
                          :lower_fence, :upper_fence, :outliers, :middle)

    ctx = compose!(
        context(tag=geom.tag),
        fill(collect(cs)),
        stroke(collect(cs)),
        linewidth(theme.line_width),

        # Box
        rectangle(
            [x - bw/2 for x in xs],
            lower_hinge, [bw],
            [uh - lh for (lh, uh) in zip(lower_hinge, upper_hinge)], tbox),

        (
            context(),

             # Whiskers
            Compose.line([[(x, lh), (x, lf)]
                          for (x, lh, lf) in zip(xs, lower_hinge, lower_fence)], tlw),

            Compose.line([[(x, uh), (x, uf)]
                          for (x, uh, uf) in zip(xs, upper_hinge, upper_fence)], tuw),

            # Fences
            Compose.line([[(x - fw/2, lf), (x + fw/2, lf)]
                          for (x, lf) in zip(xs, lower_fence)], tlf),

            Compose.line([[(x - fw/2, uf), (x + fw/2, uf)]
                          for (x, uf) in zip(xs, upper_fence)], tuf),

            stroke(collect(cs))
        ),

        svgclass("geometry"))

    # Outliers
    if !geom.suppress_outliers && aes.outliers != nothing && !isempty(aes.outliers)
        xys = collect(flatten(zip(cycle([x]), ys, cycle([c]))
                             for (x, ys, c) in zip(xs, aes.outliers, cs)))
        compose!(ctx, (context(),
            Shape.circle([x for (x, y, c) in xys], [y for (x, y, c) in xys], [theme.point_size], to), svgclass("marker"),
             stroke([theme.discrete_highlight_color(c) for (x, y, c) in xys]),
             fill([c for (x, y, c) in xys]) ))
    end

    # Middle
    if aes.middle != nothing
        compose!(ctx, (
           context(order=1),
           Compose.line([[(x - fw/2, mid), (x + fw/2, mid)]
                         for (x, mid) in zip(xs, aes.middle)], tm),
           linewidth(theme.middle_width),
           stroke([theme.middle_color(c) for c in cs])
        ))
    end

    return ctx
end


text:
#
println("As) Test of spin-angular coefficients.")

configs = [Configuration("1s^2 2s"), Configuration("1s^2 2p")]

## configs = [Configuration("1s^2 2s^2 2p^6 3s^2 3p^6 3d^10 4s^2 4p^6 4d^10 5s^2 5p^6 4f^4 6s^1") ] ## ,
           ## Configuration("1s^2 2s^2 2p^6 3s^2 3p^6 3d^10 4s^2 4p^6 4d^10 5s^2 5p^6 4f^4 5d^1"),
           ## Configuration("1s^2 2s^2 2p^6 3s^2 3p^6 3d^10 4s^2 4p^6 4d^10 5s^2 5p^6 4f^3 5d^1 6p^1"),
           ## Configuration("1s^2 2s^2 2p^6 3s^2 3p^6 3d^10 4s^2 4p^6 4d^10 5s^2 5p^6 4f^3 6s^1 6p^1"),
           ## Configuration("1s^2 2s^2 2p^6 3s^2 3p^6 3d^10 4s^2 4p^6 4d^10 5s^2 5p^6 4f^4 6d^1"),
           ## Configuration("1s^2 2s^2 2p^6 3s^2 3p^6 3d^10 4s^2 4p^6 4d^10 5s^2 5p^6 4f^3 5d^1 7p^1"),
           ## Configuration("1s^2 2s^2 2p^6 3s^2 3p^6 3d^10 4s^2 4p^6 4d^10 5s^2 5p^6 4f^4 7s^1")]

# Generate a list of relativistic configurations and determine an ordered list of subshells for these configurations
relconfList = ConfigurationR[]
for  conf in configs    wax = Basics.generate("configuration list: relativistic", conf);     append!( relconfList, wax)   end
for  i = 1:length(relconfList)    println(">> include ", relconfList[i])    end
subshellList = Basics.generate("subshells: ordered list for relativistic configurations", relconfList)
Defaults.setDefaults("relativistic subshell list", subshellList; printout=false)

# Generate the relativistic CSF's for the given subshell list
csfList = CsfR[]
for  relconf in relconfList
    newCsfs = Basics.generate("CSF list: from single ConfigurationR", relconf, subshellList);   append!( csfList, newCsfs)
end

sumSA = 0.;     sumFortran = 0.

# using Profile
# Profile.clear()
# @profile if  false
@time if false
    # Calculate angular coefficients for a scalar one- or two-particle operator
    for  leftCsf in csfList
        for rightCsf in csfList
            opa    = SpinAngular.OneParticleOperator(0, plus, true)     ## SpinAngular.TwoParticleOperator(0, plus, true)
            coeffs = SpinAngular.computeCoefficients(opa, leftCsf, rightCsf, subshellList)
            coefft = Basics.compute("angular coefficients: e-e, Ratip2013", leftCsf, rightCsf)
            ## println(">> Angular coeffients for <$leftCsf || O^($(op.rank))  || $rightCsf> = \n $coeffs ")
            if  length(coeffs)    != 0     println("\n>> Angular coeffients from SpinAngular = $coeffs ")         end
            if  length(coefft[1]) != 0     println(  ">> Angular coeffients from Ratip2013   = $(coefft[1]) ")    end
            for coeff in coeffs       global sumSA      = sumSA      + coeff.v     end
            for coeff in coefft[1]    global sumFortran = sumFortran + coeff.T     end
        end
    end
    @show sumSA, sumFortran
    #
elseif  false
    # Calculate angular coefficients for a nonscalar one- particle operator
    for  leftCsf in csfList
        for rightCsf in csfList
            opb    = SpinAngular.OneParticleOperator(1, plus, true)  ## SpinAngular.OneParticleOperator(1, minus, true)
            coeffs = SpinAngular.computeCoefficients(opb, leftCsf, rightCsf, subshellList)
            ## coefft = Basics.compute("angular coefficients: 1-p, Ratip2013", 1, leftCsf, rightCsf)
            coefft = Basics.compute("angular coefficients: 1-p, Grasp92", 0, 1, leftCsf, rightCsf)
            ## println(">> Angular coeffients for <$leftCsf || O^($(op.rank))  || $rightCsf> = \n $coeffs ")
            if  length(coeffs) != 0     println("\n>> Angular coeffients from SpinAngular = $coeffs ")    end
            if  length(coefft) != 0     println(  ">> Angular coeffients from MCT/Grasp92 = $coefft ")    end
            for coeff in coeffs       global sumSA      = sumSA      + coeff.v     end
            for coeff in coefft       global sumFortran = sumFortran + coeff.T     end
        end
    end
    @show sumSA, sumFortran
    #
elseif  false
    # Calculate angular coefficients for a nonscalar one- particle operator
    for  leftCsf in csfList
        for rightCsf in csfList
            opb    = SpinAngular.OneParticleOperator(2, plus, true)  ## SpinAngular.OneParticleOperator(2, plus, true)
            coeffs = SpinAngular.computeCoefficients(opb, leftCsf, rightCsf, subshellList)
            ## coefft = Basics.compute("angular coefficients: 1-p, Ratip2013", 2, leftCsf, rightCsf)
            coefft = Basics.compute("angular coefficients: 1-p, Grasp92", 0, 2, leftCsf, rightCsf)
            ## println(">> Angular coeffients for <$leftCsf || O^($(op.rank))  || $rightCsf> = \n $coeffs ")
            if  length(coeffs) != 0     println("\n>> Angular coeffients from SpinAngular = $coeffs ")    end
            if  length(coefft) != 0     println(  ">> Angular coeffients from MCT/Grasp92 = $coefft ")    end
            for coeff in coeffs       global sumSA      = sumSA      + coeff.v     end
            for coeff in coefft       global sumFortran = sumFortran + coeff.T     end
        end
    end
    @show sumSA, sumFortran
    #
elseif  true
    # Calculate angular coefficients for a scalar two- particle operator
    for  leftCsf in csfList
        for rightCsf in csfList
            opb    = SpinAngular.TwoParticleOperator(0, plus, true)  ## SpinAngular.TwoParticleOperator(0, plus, true)
            coeffs = SpinAngular.computeCoefficients(opb, leftCsf, rightCsf, subshellList)
            coefft = Basics.compute("angular coefficients: e-e, Ratip2013", leftCsf, rightCsf)
            ## println(">> Angular coeffients for <$leftCsf || O^($(op.rank))  || $rightCsf> = \n $coeffs ")
            if  length(coeffs)    != 0     println("\n>> Angular coeffients from SpinAngular = $coeffs ")    end
            if  length(coefft[2]) != 0     println(  ">> Angular coeffients from Ratip2013   = $(coefft[2]) ")    end
            for coeff in coeffs       global sumSA      = sumSA      + coeff.v     end
            for coeff in coefft[2]    global sumFortran = sumFortran + coeff.V     end
        end
    end
    @show sumSA, sumFortran
    #
elseif  true
    N1 = N2 = 0
    # Calculate angular coefficients for a scalar one- or two-particle operator
    op = SpinAngular.TwoParticleOperator(0, plus, true)
    for  leftCsf in csfList
        for rightCsf in csfList
            coeffs = SpinAngular.computeCoefficients(op, leftCsf, rightCsf, subshellList)
            global N2 = N2 + length(coeffs)
        end
    end
    #
    # Calculate angular coefficients for a nonscalar one- particle operator
    op = SpinAngular.OneParticleOperator(0, plus, true)
    for  leftCsf in csfList
        for rightCsf in csfList
            coeffs = SpinAngular.computeCoefficients(op, leftCsf, rightCsf, subshellList)
            global N1 = N1 + length(coeffs)
        end
    end
    @show length(csfList), N1, N2
    println("The Nd II test (Julia) has been completed")
    #
end

# Profile.print(maxdepth=30, mincount=1000)
nothing


text:
export basemeasure
export basemeasure_type

function basemeasure end

basemeasure(μ, x) = basemeasure(μ)

function basemeasure_type end

basemeasure_depth(μ::M) where {M} = basemeasure_depth(M)

basemeasure_depth(::Type{M}) where {M} = static(1) + basemeasure_depth(M)


text:
using Fast12thRoot

using Core.Intrinsics: sqrt_llvm

@noinline function throw_complex_domainerror(f, x)
    throw(DomainError(x, string("$f will only return a complex result if called with a ",
                                "complex argument. Try $f(Complex(x)).")))
end

# naive 12th root that uses built-in functions
# modern cpus e.g. Intel with SSE have built-in hardware support for sqrt
@inline twelfthroot_base(x::Float64) = cbrt(sqrt(sqrt(x)))

@inline function twelfthroot_base2(x::Float64)
    #x < zero(x) && throw(DomainError(x,"x must be non-negative"))
    if x > zero(x)
        return _twelfthroot_base2(x)
    elseif x == zero(Float64)
        return zero(Float64)
    else
        throw_complex_domainerror(:twelfthroot_base, x)
    end
    #x < zero(x) && throw_complex_domainerror(:twelfthroot_base, x)
    #cbrt(sqrt_llvm(sqrt_llvm(x)))
end

# use sqrt_llvm intrinsic as attempt to avoid double checks in sqrt
@inline function _twelfthroot_base2(x::Float64)
    cbrt(sqrt_llvm(sqrt_llvm(x)))
end

function benchtestbase(N::Int)
    sum = zero(Float64)
    srand(123456789)
    for i=1:N
        sum += twelfthroot_base(rand()*1_000_000)
    end
    sum
end

function benchtestbase2(N::Int)
    sum = zero(Float64)
    srand(123456789)
    for i=1:N
        sum += twelfthroot_base2(rand()*1_000_000)
    end
    sum
end

function benchtestfast12throot(N::Int)
    sum = zero(Float64)
    srand(123456789)
    for i=1:N
        sum += fast12throot(rand()*1_000_000)
    end
    sum
end

function benchtestpow(N::Int)
    sum = zero(Float64)
    srand(123456789)
    for i=1:N
        sum += (rand()*1_000_000)^(1/12)
    end
    sum
end

function runbench(N::Int)
    println("compiling")
    benchtestbase(1)
    benchtestbase2(1)
    benchtestfast12throot(1)
    benchtestpow(1)
    println()

    println("benchmarking sum of twelfth power of $(N) random numbers")
    println("---------------------------------------------------")
    print("fast12throot:")
    @time(benchtestfast12throot(N))
    print("x^(1/12):")
    @time(benchtestpow(N))
    print("cbrt(sqrt(sqrt(x))):")
    @time(benchtestbase(N))
    print("cbrt(sqrt_llvm(sqrt_llvm(x))):")
    @time(benchtestbase2(N))

    println()

end

runbench(10_000_000)


text:
function epipole(𝐅::AbstractArray)
    𝐔, 𝐒, 𝐕 = svd(𝐅)
    𝐞 = 𝑛(MVector{3}(𝐕[:,end]))
    Point2DH(𝐞)
end


text:
abstract type Disturbances end

function disturbances(d::Disturbances, x, t)
    @warn "disturbances not defined"
    return nothing
end

"""
    no disturbances
"""
struct NoDisturbances{T} <: Disturbances
    w::Vector{T}
end

function no_disturbances(model::ContactModel)
    NoDisturbances(zeros(model.dim.w))
end

function disturbances(d::NoDisturbances, x, t)
    return d.w
end

"""
    open-loop disturbances
"""
mutable struct OpenLoopDisturbance{W} <: Disturbances
    w::Vector{W} # nominal disturbances
    idx::Int
    cnt::Int
    N_sample::Int
end

function open_loop_disturbances(w, N_sample=Int)
    OpenLoopDisturbance(w, 0, N_sample, N_sample)
end

function disturbances(d::OpenLoopDisturbance, x, t)
    # With N_sample = 3
    # t   = 1   2   3   4   5   6   7   8   9   10
    # idx = 1   1   1  12   2   2  23   3   3   34
    # cnt = 0   1   2  30   1   2  30   1   2   30

    # reset
    if t == 1
        d.idx = 0
        d.cnt = d.N_sample
    end

    if d.cnt == d.N_sample
        d.idx += 1
        d.cnt = 0
    end

    d.cnt += 1

    return d.w[d.idx] ./ d.N_sample
end

"""
    impulse disturbances
"""
mutable struct ImpulseDisturbance{W} <: Disturbances
    w::Vector{W} # nominal disturbances
    idx::Vector{Int}
end

function impulse_disturbances(w, idx)
    ImpulseDisturbance(w, idx)
end

function disturbances(d::ImpulseDisturbance, x, t)
    for (i, w) in enumerate(d.w)
        if d.idx[i] == t
            return w
        end
    end
    return zero(d.w[1])
end

"""
    random disturbances
    Sample random disturbance forces uniformly between 0 and -w_amplitude (one-sided).
    Apply these forces along the x (and y) axes depending on the model considered.
"""
struct RandomDisturbance{W, T} <: Disturbances
    w_amp::Vector{T} # disturbance amplitude
    w::Vector{W}     # disturbances
    t::Vector{T}     # time trajectory
end

function random_disturbances(model::ContactModel, w_amp::Vector{T}, H::Int, h::T) where {T}
    if length(w_amp) == model.dim.w
        w = [rand(model.dim.w) .* w_amp for i=1:H]
    elseif length(w_amp) == 1
        w = [rand(model.dim.w) .* w_amp[1] for i=1:H]
    else
        @warn "w_amp is not of the correct size."
    end
    RandomDisturbance(w_amp, [rand(model.dim.w) .* w_amp[1] for i=1:H], [(t - 1) * h for t = 1:H])
end

function disturbances(d::RandomDisturbance, x, t)
    k = searchsortedlast(d.t, t)
    return d.w[k]
end


text:
template_plugins(::PDTN"default") = PkgTemplates.Plugin[]

function create_template(::PDTN"test", dir, user, interactive)
    return Template(;
        dir=dir,
        user=user,
        interactive=interactive,
        plugins=[
            Git(;
                ignore=String[],
                name="Roger-luo",
                email="rogerluo.rl18@gmail.com",
            )
        ]
    )
end


text:
using Gadfly, DataFrames

x = collect(11:20)
a = x .* 3 .+ 4
b = a .- 1
c = a .- 2
d = a .- 3

df = DataFrame(x=x, a=a, b=b, c=c, d=d, s=repeat(["A"],inner=10))
plot(df, x="x", upper_fence="a", upper_hinge="b", lower_hinge="c", lower_fence="d",
     Geom.boxplot)


text:
# MIT License
#
# Copyright (c) 2020 Joshua E Gentry

# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
#
# The above copyright notice and this permission notice shall be included in all
# copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
# SOFTWARE.
using .Pipeline

@enum ResultRendererType begin
    RENDERER_UNKNOWN
    RENDERER_ALTO
    RENDERER_HOCR
    RENDERER_LSTM_BOX
    RENDERER_PDF
    RENDERER_TEXT
    RENDERER_TSV
    RENDERER_UNLV
    RENDERER_WORD_BOX
end
# =========================================================================================
"""
    mutable struct TessPipeline
        inst::TessInst
        ptr::Ptr{Cvoid}
        tasks::Vector{PipelineTask}
        types::Vector{ResultRendererType}
    end

Allows the client to process multiple images in sequence.

__Values:__

| Name   | Description
| :----- | :----------
| inst   | The TessInst that will be processing the images.
| ptr    | The Tesseract pipeline that will be outputing the data.
| task   | The list of background tasks associated with pipeline.
| types  | The list of renderers used in the pipeline.

See also: [`tess_run_pipeline`](@ref)
"""
mutable struct TessPipeline
    inst::TessInst
    ptr::Ptr{Cvoid}
    tasks::Vector{PipelineTask}
    types::Vector{ResultRendererType}
end

# =========================================================================================
"""
    TessPipeline(
        inst::TessInst
    )

Construct a new instance of the object.  Registers a finalizer to clean up allocated
resources.

| Name   | Description
| :----- | :----------
| inst   | The TessInst that will be processing the images.

See also: [`tess_run_pipeline`](@ref)
"""
function TessPipeline(
            inst::TessInst
        )
    local retval = TessPipeline(inst, C_NULL, Vector{PipelineTask}(), Vector{ResultRendererType}())

    finalizer(retval) do obj
        tess_delete!(obj, false)
    end

    return retval
end

# =========================================================================================
"""
    show(
        io::IO,
        inst::TessPipeline
    )::Nothing

Display summary information about the `TessPipeline` instance.

__Arguments:__

| T | Name  | Default | Description
|---| :---- | :------ | :----------
| R | io    |         | The stream to write the information to.
| R | inst  |         | The `TessPipeline` instance to display info about.
"""
function Base.show(
            io::IO,
            inst::TessPipeline
        )::Nothing
    if inst.ptr == C_NULL
        print(io, "Pipeline is empty.")
    else
        print(io, "Pipeline is populated.")
    end
    nothing
end

# =========================================================================================
"""
    unsafe_convert(
        ::Type{Ptr{Cvoid}},
        p::TessPipeline
    )

Retreieve a pointer to the Tesseract pipeline we are using.

__Arguments:__

| T | Name               | Default  | Description
|:--| :----------------- | :------- | :----------
| R | ::Type{Ptr{Cvoid}} |          | The type we "convert" the object to.
| R | p                  |          | The pipeline we are converting.
"""
Base.unsafe_convert(::Type{Ptr{Cvoid}}, p::TessPipeline) = p.ptr

# =========================================================================================
"""
    tess_delete!(
        pipe::TessPipeline,
        wait::Bool = true
    )

Release all resources associated with the pipeline.

__Arguments:__

| T | Name | Default | Description
|:--| :--- | :------ | :----------
| R | pipe |         | The pipeline to release.
| O | wait | `true`  | Wait for the tasks to complete?

See also: [`TessPipeline`](@ref)
"""
function tess_delete!(
            pipe::TessPipeline,
            wait::Bool = true
        )
    if pipe.ptr != C_NULL
        ccall(
            (:TessDeleteResultRenderer, TESSERACT),
            Cvoid,
            (Ptr{Cvoid},),
            pipe
        )
        pipe.ptr = C_NULL

        for task in pipe.tasks
            pipeline_stop_task(task, wait)
        end
        resize!(pipe.tasks, 0)
        resize!(pipe.types, 0)
    end
end

# =========================================================================================
"""
    pipeline_add_renderer(
        pipe::TessPipeline,
        ptr::Ptr{Cvoid},
        task::PipelineTask,
        type::ResultRendererType
    )

Add a renderer to the pipeline.

__Arguments:__

| T | Name | Default | Description
|:--| :--- | :------ | :----------
| R | pipe |         | The pipeline to add the renderer to.
| R | ptr  |         | The Tesseract renderer to add.
| R | task |         | The task that processes the results from the renderer.
| R | type | unknown | The type of the renderer.

__Details:__

This method actually inserts the renderer at the front of the pipeline.  The pipeline is
really a series of independant extractors so the order doesn't matter.
"""
function pipeline_add_renderer(
            pipe::TessPipeline,
            ptr::Ptr{Cvoid},
            task::PipelineTask,
            type::ResultRendererType = RENDERER_UNKNOWN
        )::Nothing

    # ---------------------------------------------------------------------------------------------
    # If we have a current renderer, add it AFTER this new renderer.
    if pipe.ptr != C_NULL
        ccall(
            (:TessResultRendererInsert, TESSERACT),
            Cvoid,
            (Ptr{Cvoid}, Ptr{Cvoid}),
            ptr,
            pipe.ptr
        )
    end

    # ---------------------------------------------------------------------------------------------
    # Save the pointer in our structure.
    pipe.ptr = ptr

    # ---------------------------------------------------------------------------------------------
    # Add the task to the end of our list.
    push!(pipe.tasks, task)
    push!(pipe.types, type)

    nothing
end

# =========================================================================================
"""
    pipeline_add_image(
        pipe::TessPipeline,
        inst::TessInst
    )

Push the currently recognized image out to the renderers in the pipeline.

__Arguments:__

| T | Name | Default | Description
|:--| :--- | :------ | :----------
| R | pipe |         | The pipeline to push the recognized image to.
| R | inst |         | The tesseract image to pull the text from.

__Details:__

`tess_recognize` must be called before this method is called so the renderers can pick up
the recognized text.
"""
function pipeline_add_image(
            pipe::TessPipeline
        )::Bool

    if pipe.ptr == C_NULL
        @error "Pipline has been freed."
        return false
    end

    if is_valid(pipe.inst) == false
        @error "Instance has been freed."
        return false
    end

    local retval = ccall(
        (:TessResultRendererAddImage, TESSERACT),
        Cint,
        (Ptr{Cvoid}, Ptr{Cvoid}),
        pipe,
        pipe.inst
    )

    return retval == 1
end

# =========================================================================================
"""
    pipeline_start(
        pipe::TessPipeline,
        title::AbstractString
    )::Bool

Notifies the pipeline that we are starting, allowing the renderers to write out header
information before the content.  Returns `false` if Tesseract reports an error.

__Arguments:__

| T | Name  | Default | Description
|:--| :---- | :------ | :----------
| R | pipe  |         | The pipeline to start.
| R | title |         | An optional title for the document.
"""
function pipeline_start(
            pipe::TessPipeline,
            title::AbstractString
        )::Bool

    if pipe.ptr == C_NULL
        @warn "Pipeline has no renderers."
        return false
    end

    local result = @threadcall(
        (:TessResultRendererBeginDocument, TESSERACT),
        Cint,
        (Ptr{Cvoid}, Cstring),
        pipe,
        title
    );

    return result == 1
end

# =========================================================================================
"""
    pipeline_finish(
        pipe::TessPipeline
    )::Bool

Tell the renders to write out any footer information they want at the end of the document.
Returns `false` if Tesseract reports an error.

__Arguments:__

| T | Name  | Default | Description
|:--| :---- | :------ | :----------
| R | pipe  |         | The pipeline to complete.
"""
function pipeline_finish(
            pipe::TessPipeline
        )::Bool
    local retval = false

    if pipe.ptr == C_NULL
        @warn "Pipeline has no renderers."
        return false
    end

    local result = @threadcall(
        (:TessResultRendererEndDocument, TESSERACT),
        Cint,
        (Ptr{Cvoid},),
        pipe
    )

    return result == 1
end

# =========================================================================================
"""
    tess_run_pipeline(
        func::Function,
        pipe::TessPipeline,
        title::AbstractString = ""
    )::Bool

Run the pipeline against a series of images.

__Arguments:__

| T | Name  | Default | Description
|:--| :---- | :------ | :----------
| R | func  |         | The callback to call for the images.
| R | pipe  |         | The pipline we're going to execute.
| O | title | `""`    | Optional title for the document.

__Details:__

This is the primary function of the `TessPipeline` object.  This method is used to combine
output from multiple images into a single document (or multiple documents of different
types).

The argument passed back to the function is a method with the form:

    function add(Pix, Integer)::Bool

The client calls this method with the image and pixels per inch of the image.  The image
will be processed immediately and added to the pipeline's output.

If you the callback returns false then the `tess_run_pipeline()` will return `false`.  If it
doesn't return a boolean or returns `true` then `tess_run_pipeline()` returns `true`
assuming there are no other errors.

This method can only be used "once".  After it's called new output needs to be added to the
pipeline if you want to reuse the [`TessPipeline`](@ref) object.

__Examples:__

```jldoctest
using Tesseract

# Generate some pages to load.
write("page01.tiff", sample_tiff())
write("page02.tiff", sample_tiff())
write("page03.tiff", sample_tiff())

download_languages() # Make sure we have the data files.

instance = TessInst()
pipeline = TessPipeline(instance)

book = tess_pipeline_text(pipeline)

tess_run_pipeline(pipeline, "My First Book") do add
    add(pix_read("page01.tiff"), 72)
    add(pix_read("page02.tiff"), 72)
    add(pix_read("page03.tiff"), 72)
end

for line in split(book[], "\\n")[1:10]
    println(line)
end

# output

No one would have believed in the last years of the

the nineteenth century that this world was being watched
watched keenly and closely by intelligences greater than
than man’s and yet as mortal as his own; that as men busied
busied themselves about their various concerns they were
were scrutinised and studied, perhaps almost as narrowly as
as a man with a microscope might scrutinise the transient
transient creatures that swarm and multiply in a drop of
```

See also: [`TessPipeline`](@ref), [`tess_pipeline_text`](@ref)
"""
function tess_run_pipeline(
            func::Function,
            pipe::TessPipeline,
            title::AbstractString = ""
        )::Bool
    local result = true

    # -------------------------------------------------------------------------------------
    # Make sure there is somethign to do.
    if pipe.ptr == C_NULL
        @error "No renderers attached to the pipeline."
        return false
    end

    # ---------------------------------------------------------------------------------
    # Allow the renderers to write out a header.
    if pipeline_start(pipe, title) == false
        tess_delete!(pipe)
        return false
    end
    # See https://github.com/tesseract-ocr/tesseract/issues/2700
    if RENDERER_ALTO ∈ pipe.types && tess_get_input_name(pipe.inst) === nothing
        tess_set_input_name(pipe.inst, title == "" ? " " : title)
    end

    try
        # ---------------------------------------------------------------------------------
        # This function is called by the client to process an image.
        function process(
                    pix::Pix,
                    res::Integer
                )::Bool
            # -----------------------------------------------------------------------------
            # Set the image.
            tess_image(pipe.inst, pix)
            tess_resolution(pipe.inst, res)

            # -----------------------------------------------------------------------------
            # Look for text.
            if tess_recognize(pipe.inst) == false
                result = false
                return false
            end

            # -----------------------------------------------------------------------------
            # Allow the renderers to process the new text.
            if pipeline_add_image(pipe) == false
                result = false
                return false
            end

            # -----------------------------------------------------------------------------
            # Notify the tasks that there is probably new data to read.
            for task in pipe.tasks
                notify(task.event)
            end

            return true
        end

        # ---------------------------------------------------------------------------------
        # Call the client to provide all the images for processing.
        local retval = func(process)

        if isa(retval, Bool) && retval == false
            result = false
        end

    finally

        # ---------------------------------------------------------------------------------
        # Allow the renderers to finish the documents.
        if pipeline_finish(pipe) == false
            result = false
        end

        # ---------------------------------------------------------------------------------
        # Delete the renderers so they finish writing out their data.
        tess_delete!(pipe)
    end

    return result
end

# =========================================================================================
"""
    tess_run_pipeline(
        pipe::TessPipeline,
        filename::AbstractString;
        retryConfig::AbstractString = "",
        timeout::Integer = Int32(0)
    )::Bool

Run the pipline against a file that contains the list of images to process or a TIFF
file with multiple images.  Returns `false` if the process fails.

__Arguments:__

| T | Name        | Default | Description
|:--| :---------- | :--------- | :----------
| R | pipe        |            | The pipline we're going to execute.
| R | filename    |            | The name of the file to read.
| K | retryConfig | `""`       | A configuration file to load if an image cannot be processed.
| K | timeout     | `Int32(0)` | The maximum time in milliseconds to spend per page.

__Details:__

This method provides a simplified version of the `tess_run_pipeline` function.  This
uses method exposed by the Tesseract library to process a list of files or a single
TIFF file without additional input.

This method can only be used "once".  After it's called new output needs to be added to the
pipeline if you want to reuse the [`TessPipeline`](@ref) object.

__Examples:__

Providing a file with a list of images to process:

```jldoctest
using Tesseract
using Suppressor

# Generate some pages to load.
write("page01.tiff", sample_tiff())
write("page02.tiff", sample_tiff())
write("page03.tiff", sample_tiff())

download_languages() # Make sure we have the data files.

instance = TessInst()
pipeline = TessPipeline(instance)

book = tess_pipeline_text(pipeline)

# Create a text file with the files to process.
open("pages.lst", create = true, write = true) do io
    println(io, "page01.tiff")
    println(io, "page02.tiff")
    println(io, "page03.tiff")
end

@suppress tess_run_pipeline(pipeline, "pages.lst")

for line in split(book[], "\\n")[1:10]
    println(line)
end

# output

No one would have believed in the last years of the

the nineteenth century that this world was being watched
watched keenly and closely by intelligences greater than
than man’s and yet as mortal as his own; that as men busied
busied themselves about their various concerns they were
were scrutinised and studied, perhaps almost as narrowly as
as a man with a microscope might scrutinise the transient
transient creatures that swarm and multiply in a drop of
```

Using a multipage TIFF:

```jldoctest
using Tesseract
using Suppressor

# Generate some pages to load.
pix_write_tiff("book.tiff", sample_pix())
pix_write_tiff("book.tiff", sample_pix(); append = true)
pix_write_tiff("book.tiff", sample_pix(); append = true)

download_languages() # Make sure we have the data files.

instance = TessInst()
pipeline = TessPipeline(instance)

book = tess_pipeline_text(pipeline)

@suppress tess_run_pipeline(pipeline, "book.tiff")

for line in split(book[], "\\n")[1:10]
    println(line)
end

# output

No one would have believed in the last years of the

the nineteenth century that this world was being watched
watched keenly and closely by intelligences greater than
than man’s and yet as mortal as his own; that as men busied
busied themselves about their various concerns they were
were scrutinised and studied, perhaps almost as narrowly as
as a man with a microscope might scrutinise the transient
transient creatures that swarm and multiply in a drop of
```

See also: [`TessPipeline`](@ref), [`tess_pipeline_text`](@ref)
"""
function tess_run_pipeline(
            pipe::TessPipeline,
            filename::AbstractString;
            retryConfig::AbstractString = "",
            timeout::Integer = Int32(0)
        )::Bool
    local result = true

    # -------------------------------------------------------------------------------------
    # Make sure there is something to do.
    if pipe.ptr == C_NULL
        @error "No renderers attached to the pipeline."
        return false
    end

    # -------------------------------------------------------------------------------------
    # Make sure the instanc is still valid.
    if is_valid(pipe.inst) == false
        @error "Instance has been freed."
        return false
    end

    # -------------------------------------------------------------------------------------
    # Make sure the file name is provided.
    if filename == ""
        @error "No filename provided."
        return false
    end

    # -------------------------------------------------------------------------------------
    # Make the call.
    local retval = ccall(
        (:TessBaseAPIProcessPages, TESSERACT),
        Cint,
        (Ptr{Cvoid}, Cstring, Cstring, Cint, Ptr{Cvoid}),
        pipe.inst,
        filename,
        isempty(retryConfig) ? C_NULL : retryConfig,
        timeout,
        pipe
    )

    # ---------------------------------------------------------------------------------
    # Delete the renderers so they finish writing out their data.
    tess_delete!(pipe)

    return result
end


text:
using
  FourierFlows.TracerAdvDiff,
  FFTW

using LinearAlgebra: mul!, ldiv!

# Dictionary of (stepper, nsteps) pairs to test. Each stepper is tested by
# stepping forward nstep times.
steppersteps = Dict([
  ("ForwardEuler", 2000),
  ("FilteredForwardEuler", 2000),
  ("AB3", 400),
  ("FilteredAB3", 400),
  ("RK4", 40),
  ("FilteredRK4", 40),
  ("ETDRK4", 40),
  ("FilteredETDRK4", 40)
])

dualsteppersteps = Dict([
  ("DualRK4", 40),
  ("DualFilteredRK4", 40),
  ("DualETDRK4", 40),
  ("DualFilteredETDRK4", 40),
])

function tracerdiffusionproblem(nx=16, ny=2, Lx=2π, kap=1e-2; nsteps=100, stepper="ForwardEuler")
  k1 = 2π/Lx
  dt = 1e-9 / (kap*k1^2) # resolved?
  prob = TracerAdvDiff.ConstDiffProblem(nx=nx, Lx=Lx, ny=ny, kap=kap, dt=dt, stepper=stepper, steadyflow=true)
  g = prob.grid

  t1 = dt*nsteps
  c0 = @. sin(k1*g.X)
  c1 = @. exp(-kap*k1^2*t1) * sin(k1*g.X) # analytical solution

  set_c!(prob, c0)
  stepforward!(prob, nsteps)
  updatevars!(prob)
  prob, c0, c1, nsteps
end

function tracerdiffusiontest(args...; kwargs...)
  prob, c0, c1, nsteps = tracerdiffusionproblem(args...; kwargs...)
  isapprox(c1, prob.vars.c, rtol=prob.grid.nx^2*nsteps*1e-13)
end

for (stepper, steps) in steppersteps
  @test tracerdiffusiontest(; stepper=stepper, nsteps=steps)
end

#=
"""
    testtwodturbstepforward(n, L, nu, nnu; kwargs...)

Build a twodturb initial value problem and use it to test time-stepping methods.
We integrate a random IC from 0 to t=tf. The amplitude of the initial condition
is kept low (e.g. multiplied by 1e-5) so that nonlinear terms do not come into play.
This way we can compare the final state qh(t=tf) with qh(t=0)*exp(-nu k^nnu tf).
We choose linear drag (nnu=0) so that we can test the energy since in that
case E(t=tf) = E(t=0)*exp(-2nu tf).
"""
function testtwodturbstepforward(n=64, L=2π, nu=1e-2, nnu=0; nsteps=100, stepper="ForwardEuler")
  dt = tf/nsteps
  prob = TwoDTurb.InitialValueProblem(nx=n, Lx=L, ny=n, Ly=L, nu=nu, nnu=nnu, dt=dt, stepper=stepper)
  g, p, v, s = prob.grid, prob.params, prob.vars, prob.state
  E = Diagnostic(energy, prob, nsteps=nsteps)
  diags = [E]

  qih = 1e-5*(rand(g.nkr, g.nl) + im*rand(g.nkr, g.nl)) # initial condition
  qih[1, 1] = 0

  # Remove high wavenumber power from IC
  cutoff = 0.3
  highwavenumbers = @. sqrt((g.Kr*g.dx/π)^2 + (g.Lr*g.dy/π)^2 ) > cutoff
  @. qih[highwavenumbers] = 0

  qi = irfft(qih, g.nx)
  qih = rfft(qi)
  q0h = deepcopy(qih)

  TwoDTurb.set_q!(prob, qi)
  E0 = energy(s, v, g)
  stepforward!(prob, diags, nsteps)

  qfh = @. q0h*exp(-p.nu*g.KKrsq^p.nnu*s.step*dt)
  Ef = @. E0*exp(-2*p.nu*s.step*dt)
  isapprox(qfh, prob.state.sol, rtol=n^2*nsteps*1e-13) &  isapprox(Ef, E[E.count], rtol=n^2*nsteps*1e-13)
end

"""
    testverticallyfourierstepforward(n, L, nu, nnu; kwargs...)

Build a Vertically Fourier Boussinesq initial value problem and use it to test
time-stepping methods.
"""
function testverticallyfourierstepforward(n=64, L=2π, nu=1e-2, nnu=0; nsteps=100, stepper="DualRK4")
  dt = tf/nsteps
  prob = VerticallyFourierBoussinesq.Problem(
    nx=n, Lx=L, ny=n, Ly=L, nu0=nu, nnu0=nnu, nu1=nu, nnu1=nnu, dt=dt, stepper=stepper)
  g, p, v, s = prob.grid, prob.params, prob.vars, prob.state

  Zih = 1e-5*(rand(g.nkr, g.nl) + im*rand(g.nkr, g.nl)) # initial condition
  Zih[1, 1] = 0

  uih = 1e-5*(rand(g.nk, g.nl) + im*rand(g.nk, g.nl)) # initial condition
  vih = 1e-5*(rand(g.nk, g.nl) + im*rand(g.nk, g.nl)) # initial condition
  pih = 1e-5*(rand(g.nk, g.nl) + im*rand(g.nk, g.nl)) # initial condition
  uih[1, 1] = 0
  vih[1, 1] = 0
  pih[1, 1] = 0

  # Remove high wavenumber power from IC
  cutoff = 0.3
  highwavenumbersr = @. sqrt((g.Kr*g.dx/π)^2 + (g.Lr*g.dy/π)^2 ) > cutoff
  highwavenumbersc = @. sqrt((g.K*g.dx/π)^2 + (g.L*g.dy/π)^2 ) > cutoff
  Zih[highwavenumbersr] .= 0
  uih[highwavenumbersc] .= 0
  vih[highwavenumbersc] .= 0
  pih[highwavenumbersc] .= 0

  Zi = irfft(Zih, g.nx)
  Zih = rfft(Zi)
  Z0h = deepcopy(Zih)

  ui = ifft(uih)
  uih = fft(ui)
  u0h = deepcopy(uih)

  vi = ifft(vih)
  vih = fft(vi)
  v0h = deepcopy(vih)

  pi = ifft(vih)
  pih = fft(vi)
  p0h = deepcopy(vih)

  VerticallyFourierBoussinesq.set_Z!(prob, Zi)
  stepforward!(prob, nsteps)

  Zfh = @. Z0h*exp(-p.nu0*g.KKrsq^p.nnu0*s.step*dt)
  isapprox(Zfh, prob.state.solr, rtol=n^2*nsteps*1e-13)
end

# Run the tests
n = 64
tf = 0.2

for (stepper, steps) in steppersteps
  @test testtwodturbstepforward(n; stepper=stepper, nsteps=steps)
end

for (stepper, steps) in dualsteppersteps
  @test testverticallyfourierstepforward(n; stepper=stepper, nsteps=steps)
end
=#


text:
# definition of lyapunov
using DrWatson
@quickactivate "NonlinearDynamicsTextbook"
include(srcdir("style.jl"))
using DynamicalSystems, PyPlot, Random

using LinearAlgebra
lo = Systems.lorenz([20,20,20.0])
X₁ = trajectory(lo, 45)
u₂ = get_state(lo) .+ 1e-6
X₂ = trajectory(lo, 45, u₂)
δ  = norm.(X₂.data .- X₁.data)
λ = lyapunov(lo, 100000)

fig = figure(figsize=(figx/2,figy/2))
ax = gca()
ax.plot(0:0.01:45, log.(δ), c = COLORS[1], label ="\$\\ln(\\delta(t)))\$")
ax.set_yticks(-12:4:4)
ax.set_xlabel("\$t\$", labelpad=-20)
# Lyapunov
ax.plot([0, 15] .+ 4, λ .* [0, 15] .- 13, color = COLORS[2])
ax.text(12, -9, "\$\\lambda_1\$=$(round(λ;digits=2))", color = COLORS[2])
ax.legend(;handlelength = 1.0)
xticks(0:15:45)
tight_layout(pad = 0.25)
subplots_adjust(bottom = 0.24, left = 0.12)
wsave(plotsdir("3", "lyapunov"), fig)

max_stars_repo_name: JuliaDynamics/NonlinearDynamicsTextbook

max_issues_repo_name: JuliaDynamics/NonlinearDynamicsTextbook

max_forks_repo_name: JuliaDynamics/NonlinearDynamicsTextbook


text:
function calculate_isometries(bondenv::AbstractTensorMap, D::Integer)
    # convention:
    # 3. e1u |-<-  -<-| w2u 1.
    # 2. e1d |->-  ->-| w2d 4.
    # ==========
    # -<-|u>-<-
    #

    σ = id(domain(bondenv)[1])

    u, s, v = tsvd(σ, (1,), (2,); trunc=truncdim(D))
    u = TensorMap(randn, codomain(u), domain(u))
    v = TensorMap(randn, codomain(v), domain(v))
    # u = isometry(domain(bondenv)[1], D)
    # v = isometry(codomain(bondenv)[1], D)
    @tensor B1[w2u e1d; e1u w2d] := bondenv[w2u 1; 2 w2d] * (u')[e1d; 1] * u[2; e1u]
    @tensor P1[e1d; w2d] := bondenv[2 1; 2 w2d] * (u')[e1d; 1]
    B1inv = permute(pinv(B1), (2, 4), (3, 1))

    R1 = permute(permute(P1, (), (2, 1)) * B1inv, (2,), (1,))
    u, s, v = tsvd(u * R1, (1,), (2,); trunc=truncdim(D))

    @tensor B2[w2u e1d; e1u w2d] := bondenv[1 e1d; e1u 2] * (v)[w2u; 1] * (v')[2; w2d]
    @tensor P2[e1d; w2d] := bondenv[2 e1d; 2 1] * (v')[1; w2d]
    B2inv = permute(pinv(B2), (2, 4), (3, 1))

    R2 = permute(permute(P2, (), (2, 1)) * B2inv, (2,), (1,))
    u, s, v = tsvd(R2 * v, (1,), (2,); trunc=truncdim(D))

    for i in 1:20
        print(PEPOKit.fidelity(bondenv, u, s, v), "\n")
        @tensor B1[w2u e1d; e1u w2d] = bondenv[w2u 1; 2 w2d] * (u')[e1d; 1] * u[2; e1u]
        @tensor P1[e1d; w2d] = bondenv[2 1; 2 w2d] * (u')[e1d; 1]
        B1inv = permute(pinv(B1), (2, 4), (3, 1))

        R1 = permute(permute(P1, (), (2, 1)) * B1inv, (2,), (1,))
        u, s, v = tsvd(u * R1, (1,), (2,); trunc=truncdim(D))

        @tensor B2[w2u e1d; e1u w2d] = bondenv[1 e1d; e1u 2] * (v)[w2u; 1] * (v')[2; w2d]
        @tensor P2[e1d; w2d] = bondenv[2 e1d; 2 1] * (v')[1; w2d]
        B2inv = permute(pinv(B2), (2, 4), (3, 1))

        R2 = permute(permute(P2, (), (2, 1)) * B2inv, (2,), (1,))
        u, s, v = tsvd(R2 * v, (1,), (2,); trunc=truncdim(D))
    end
    return u, s, v
end

function fidelity(bondenv, u, s, v)
    @tensor trϕϕ =
        bondenv[w2u e1d; e1u w2d] *
        u[e1u; s1u] *
        s[s1u; s2u] *
        v[s2u; w2u] *
        (v')[w2d; s1d] *
        s[s1d; s2d] *
        (u')[s2d; e1d]
    @tensor trρρ = bondenv[1 2; 1 2]
    @tensor trρϕ = bondenv[1 e1d; 1 w2d] * (v')[w2d; s1d] * s[s1d; s2d] * (u')[s2d; e1d]

    return real(trρϕ / sqrt(trϕϕ * trρρ))
end

max_stars_repo_name: jack-dunham/PEPOKit

max_issues_repo_name: jack-dunham/PEPOKit

max_forks_repo_name: jack-dunham/PEPOKit


text:
module PseudoPotentials
using AtomicLevels
import AtomicLevels: spectroscopic_label, AbstractOrbital
using AtomicPotentials
import AtomicPotentials: islocal
using AngularMomentumAlgebra
import AngularMomentumAlgebra: jmⱼ, spin
using WignerSymbols
using HalfIntegers
using PrettyTables

abstract type AbstractPseudoPotential{T} <: AbstractPotential{T} end

islocal(::AbstractPseudoPotential) = false

include("gaussian_expansions.jl")
include("nonrelativistic_pseudopotentials.jl")
include("relativistic_pseudopotentials.jl")
include("parse.jl")
include("misc_pseudopotentials.jl")

export @PP_str, @pc_str,
    AbstractPotential, PseudoPotential, RelativisticPseudoPotential,
    charge, ground_state, spin_orbit_potential

end # module


text:
export MetaBayesTree, BayesTree
export CSMHistoryTuple

## ========================================================================================================================
## Bayes Trees
## ========================================================================================================================

abstract type AbstractBayesTree end

# TODO DEV MetaGraphs bayes tree, will potentially also make a LightBayesTree, CloudBayesTree,
"""
$(TYPEDEF)
Data structure for the Bayes (Junction) tree, which is used for inference and constructed from a given `::AbstractDFG`.
"""
mutable struct MetaBayesTree <: AbstractBayesTree
  bt::MetaDiGraph{Int,Float64}
  btid::Int
  frontals::Dict{Symbol,CliqueId{Int}}
  eliminationOrder::Vector{Symbol}
  buildTime::Float64
end

const BayesTree = MetaBayesTree

MetaBayesTree() = MetaBayesTree(MetaDiGraph{Int,Float64}(), 0, Dict{AbstractString, Int}(), Symbol[], 0.0)

Base.propertynames(x::MetaBayesTree, private::Bool=false) = (:bt, :btid, :cliques, :frontals, :variableOrder, :buildTime)

Base.getproperty(x::MetaBayesTree,f::Symbol) = begin
    if f == :cliques
      @warn "Don't use cliques field directly, use eg. getClique(tree, cliqId)" maxlog=1
      d = Dict{Int,Any}()
      for (k,v) in x.bt.vprops
        d[k] = v[:clique]
      end
      return d
    elseif f == :variableOrder
      @warn "MetaBayesTree.variableOrder has been deprecated in favor of .eliminationOrder"
      return getfield(x,:eliminationOrder)
    else
      getfield(x,f)
    end
  end

function Base.setproperty!(x::MetaBayesTree, f::Symbol, val)
  if f == :cliques
    @warn "`setproperty!(clique)` Don't use cliques field directly, use eg. addClique(tree, cliqId)" maxlog=1
    for (k,v) in val
      set_prop!(x.bt, k, :clique, v)
    end
  else
    setfield!(x,f,val)
  end
end

function getMessageChannels(tree::MetaBayesTree)

  d = Dict{Int, NamedTuple{(:upMsg, :downMsg),Tuple{Channel{LikelihoodMessage},Channel{LikelihoodMessage}}}}()

  for (k,e) in tree.bt.eprops
    d[k.dst] = (upMsg=e[:upMsg], downMsg=e[:downMsg])
  end

  return d

end

function Base.show(io::IO, mbt::MetaBayesTree)
  printstyled(io, "MetaBayesTree\n", color=:blue)
  println(io, "  Nr cliques:  ", length(mbt.cliques))

  # TODO ad dmore stats: max depth, widest point, longest chain, max clique size, average nr children

  nothing
end

Base.show(io::IO, ::MIME"text/plain", mbt::MetaBayesTree) = show(io, mbt)

"""
    $TYPEDEF

Container for upward tree solve / initialization.

DevNotes
- TODO more direct clique access (cliq, parent, children), for multi-process solves
"""
mutable struct CliqStateMachineContainer{BTND, G <: AbstractDFG, InMemG <: InMemoryDFGTypes, BT <: AbstractBayesTree}
  dfg::G
  cliqSubFg::InMemG
  tree::BT
  cliq::TreeClique
  incremental::Bool
  drawtree::Bool
  dodownsolve::Bool
  delay::Bool
  opts::SolverParams
  refactoring::Dict{Symbol, String}
  oldcliqdata::BTND
  logger::SimpleLogger
  cliqId::CliqueId
  algorithm::Symbol
  init_iter::Int
  enableLogging::Bool
  solveKey::Symbol
  _csm_iter::Int
end

#TODO use @NamedTuple if julia compat > 1.5

const CSMHistoryTuple =  NamedTuple{(:timestamp, :id, :f, :csmc), Tuple{DateTime, Int, Function, CliqStateMachineContainer}}
const CSMHistory = Vector{CSMHistoryTuple}

Base.show(io::IO, o::CSMHistoryTuple) = print(io, "$(o[1]), $(o[2]), $(o[3])")

function CliqStateMachineContainer( dfg::G,
                                    cliqSubFg::M,
                                    tree::T,
                                    cliq::TreeClique,
                                    incremental::Bool,
                                    drawtree::Bool,
                                    dodownsolve::Bool,
                                    delay::Bool,
                                    opts::SolverParams,
                                    refactoring::Dict{Symbol,String}=Dict{Symbol,String}(),
                                    oldcliqdata::BTND=BayesTreeNodeData(),
                                    logger::SimpleLogger=SimpleLogger(Base.stdout);
                                    cliqId::CliqueId = cliq.id,
                                    algorithm::Symbol = :default,
                                    init_iter::Int=0,
                                    enableLogging::Bool=true,
                                    solveKey::Symbol=:default,
                                    _csm_iter::Int=0  ) where {BTND, G <: AbstractDFG, M <: InMemoryDFGTypes, T <: AbstractBayesTree}
  #
  CliqStateMachineContainer{BTND, G, M, T}( dfg,
                                            cliqSubFg,
                                            tree,
                                            cliq,
                                            incremental,
                                            drawtree,
                                            dodownsolve,
                                            delay,
                                            opts,
                                            refactoring,
                                            oldcliqdata,
                                            logger,
                                            cliqId,
                                            algorithm,
                                            init_iter,
                                            enableLogging,
                                            solveKey,
                                            _csm_iter )
  #
end

# TODO resolve name conflict
function DFG.compare( cs1::CliqStateMachineContainer{BTND1, T1, InMemG1, BT1},
                      cs2::CliqStateMachineContainer{BTND2, T2, InMemG2, BT2};
                      skip::Vector{Symbol}=Symbol[] ) where {BTND1, T1 <: AbstractDFG, InMemG1 <: InMemoryDFGTypes, BT1 <: AbstractBayesTree, BTND2, T2 <: AbstractDFG, InMemG2 <: InMemoryDFGTypes, BT2 <: AbstractBayesTree}
  #
  BTND1 == BTND2 ? nothing : @warn("oldcliqdata::$BTND1 != oldcliqdata::$BTND2")
  T1 == T2 ? nothing : @warn("dfg::$T1 != dfg::$T2")
  InMemG1 == InMemG2 ? nothing : @warn("cliqSubFg::$InMemG1 != cliqSubFg::$InMemG2")
  BT1 == BT2 ? nothing : @warn("tree::$BQ1 != tree::$BT2")

  TP = true
  @warn "Skipping compare of CSMC.dfg and .cliqSubFg"
  # TP = TP && compare(cs1.dfg,  cs2.dfg)
  # TP = TP && compare(cs1.cliqSubFg,  cs2.cliqSubFg)
  @warn "Skipping compare of CSMC.tree"
  # TP = TP && compare(cs1.tree,  cs2.tree)
  TP = TP && compare(cs1.cliq,  cs2.cliq)
  TP = TP && compare(cs1.cliqId,  cs2.cliqId)
  TP = TP && length(cs1.parentCliq) == length(cs2.parentCliq)
  for i in 1:length(cs1.parentCliq)
    TP = TP && compare(cs1.parentCliq[i],  cs2.parentCliq[i])
  end
  TP = TP && length(cs1.childCliqs) == length(cs2.childCliqs)
  for i in 1:length(cs1.childCliqs)
    TP = TP && compare(cs1.childCliqs[i],  cs2.childCliqs[i])
  end
  TP = TP && compare(cs1.incremental,  cs2.incremental)
  TP = TP && compare(cs1.drawtree,  cs2.drawtree)
  TP = TP && compare(cs1.dodownsolve,  cs2.dodownsolve)
  TP = TP && compare(cs1.delay,  cs2.delay)
  @warn "skipping compare on csmc.opts::SolverParams"
  # TP = TP && compare(cs1.opts,  cs2.opts)
  TP = TP && compare(cs1.refactoring,  cs2.refactoring)
  # TP = TP && compare(cs1.oldcliqdata,  cs2.oldcliqdata)
  # TP = TP && compare(cs1.logger,  cs2.logger)

  return TP
end

#


text:
# minimal definition of a custom kernel
struct MyKernel <: SimpleKernel end

KernelFunctions.kappa(::MyKernel, d2::Real) = exp(-d2)
KernelFunctions.metric(::MyKernel) = SqEuclidean()

@testset "custom" begin
    @test kappa(MyKernel(), 3) == kappa(SqExponentialKernel(), 3)
    @test kernelmatrix(MyKernel(), [1 2; 3 4], [5 6; 7 8]) == kernelmatrix(SqExponentialKernel(), [1 2; 3 4], [5 6; 7 8])
    @test kernelmatrix(MyKernel(), [1 2; 3 4]) == kernelmatrix(SqExponentialKernel(), [1 2; 3 4])
end


text:
# test/runtests.jl

using LinearMapsAA
using Test: @test, @testset, detect_ambiguities

include("multiply.jl")

list = [
"ambiguity"
"identity"
"kron"
"cat"
"getindex"
"setindex"
"block_diag"
"tests"
]

for file in list
    @testset "$file" begin
        include("$file.jl")
    end
end

@testset "ambiguities" begin
    @test length(detect_ambiguities(LinearMapsAA)) == 0
end