text:
using Statistics
using Plots
using BenchmarkTools
using FFTW
using Statistics
using Optim
using Images, FileIO, ImageIO
using Printf
using Revise
using Profile
using LinearAlgebra
using Distributions
using FITSIO

push!(LOAD_PATH, pwd()*"/main")
using DHC_2DUtils
#include("Deriv_Utils_New.jl")
#import Deriv_Utils_New
push!(LOAD_PATH, pwd()*"/scratch_NM")
using Deriv_Utils_New
using Data_Utils
using Visualization
using LossFuncs
#Remove later

#Jacobian Test Funcs
function derivtestS1(Nx)
   eps = 1e-4
   fhash = fink_filter_hash(1, 8, nx=Nx, pc=1, wd=1)
   im = rand(Float64, Nx, Nx).*0.1
   im[6,6]=1.0
   im1 = copy(im)
   im0 = copy(im)
   im1[2,3] += eps/2
   im0[2,3] -= eps/2
   blarg = Deriv_Utils_New.wst_S1_deriv(im, fhash)
   der0=DHC_compute(im0,fhash,doS2=false,norm=false)
   der1=DHC_compute(im1,fhash,doS2=false,norm=false)
   dS = (der1-der0) ./ eps

   diff_old = dS[3:end]-blarg[2,3,:]

   println(dS[3:end])
   println("and")
   println(blarg[2,3,:])
   println("Mean abs | Mean abs frac", mean(abs.(diff_old)), mean(abs.(diff_old./dS[3:end])))
   println("stdev: ",std(diff_old))
   return blarg[2, 3, :], dS[3:end]
end

function derivtestS1S2(Nx) #MUST RUN DHC_COMPUTE WITH NORM=FALSE.
    eps = 1e-4
    fhash = fink_filter_hash(1, 8, nx=Nx, pc=1, wd=1)
    im = rand(Float64, Nx, Nx).*0.1
    im[6,6]=1.0
    im1 = copy(im)
    im0 = copy(im)
    im1[2,3] += eps/2
    im0[2,3] -= eps/2
    dS1dp_existing = Deriv_Utils_New.wst_S1_deriv(im, fhash)
    dS1S2dp = Deriv_Utils_New.wst_S1S2_deriv(im, fhash) #Deriv_Utils_New.wst_S1S2_deriv(im, fhash)

    der0=DHC_compute(im0,fhash,doS2=true,doS20=false, norm=false)
    der1=DHC_compute(im1,fhash,doS2=true,doS20=false, norm=false)
    dSlim = (der1-der0) ./ eps
    Nf = length(fhash["filt_index"])
    lasts1ind = 2+Nf

    dS1lim23 = dSlim[3:lasts1ind]
    dS2lim23 = dSlim[lasts1ind+1:end]
    dS1dp_23 = dS1S2dp[2, 3, 1:Nf]
    dS2dp_23 = dS1S2dp[2, 3, Nf+1:end]
    #=
    println("Checking dS1dp using existing deriv")
    Nf = length(fhash["filt_index"])
    lasts1ind = 2+Nf
    diff = dS1dp_existing[2, 3, :]-dSlim[3:lasts1ind]
    println(dSlim[3:lasts1ind])
    println("and")
    println(dS1dp_existing[2, 3, :])
    println("stdev: ",std(diff))
    println("--------------------------")
    =#
    println("Checking dS1dp using dS1S2")
    derdiff = dS1dp_23 - dS1lim23
    #println(dS1lim23)
    #println("and")
    #println(dS1dp_23)
    txt= @sprintf("Range of S1 der0, der1: Max= %.3E, %.3E Min=%.3E, %.3E ",maximum(der0[3:lasts1ind]),maximum(der1[3:lasts1ind]),minimum(der0[3:lasts1ind]),minimum(der1[3:lasts1ind]))
    print(txt)
    txt = @sprintf("Checking dS1dp using S1S2 deriv, mean(abs(diff/dS1lim)): %.3E", mean(abs.(derdiff./dS1lim23)))
    print(txt)
    println("stdev: ",std(derdiff))
    txt= @sprintf("Range of dS1lim: Max= %.3E Min= %.3E",maximum(abs.(dS1lim23)),minimum(abs.(dS1lim23)))
    print(txt)
    println("--------------------------")
    println("Checking dS2dp using S1S2 deriv")
    Nf = length(fhash["filt_index"])
    println("Shape check",size(dS2dp_23),size(dS2lim23))
    derdiff = dS2dp_23 - dS2lim23 #Column major issue here?
    println("Difference",derdiff)
    txt = @sprintf("Range of S2 der0, der1: Max= %.3E, %.3E Min=%.3E, %.3E \n",maximum(der0[lasts1ind+1:end]),maximum(der1[lasts1ind+1:end]),minimum(der0[lasts1ind+1:end]),minimum(der1[lasts1ind+1:end]))
    print(txt)
    txt = @sprintf("Range of dS2lim: Max=%.3E  Min=%.3E \n",maximum(abs.(dS2lim23)),minimum(abs.(dS2lim23)))
    print(txt)
    txt = @sprintf("Difference between dS2dp and dSlim Mean: %.3E stdev: %.3E mean(abs(derdiff/dS2lim)): %.3E \n",mean(derdiff),std(derdiff), mean(abs.(derdiff./dS2lim23)))
    print(txt)

    println("Examining only those S2 coeffs which are greater than eps=1e-3 for der0")
    eps_mask = findall(der0[lasts1ind+1:end] .> 1E-3)
    der0_s2good = der0[lasts1ind+1:end][eps_mask]
    der1_s2good = der1[lasts1ind+1:end][eps_mask]
    dSlim_s2good = dS2lim23[eps_mask]
    derdiff_good = derdiff[eps_mask]
    print(derdiff_good)
    txt = @sprintf("\nRange of S2 der0, der1: Max= %.3E, %.3E Min=%.3E, %.3E \n",maximum(der0_s2good),maximum(der1_s2good),minimum(der0_s2good),minimum(der1_s2good))
    print(txt)
    txt = @sprintf("Range of dS2lim: Max=%.3E  Min=%.3E \n",maximum(abs.(dSlim_s2good)),minimum(abs.(dSlim_s2good)))
    print(txt)
    txt = @sprintf("Mean: %.3E stdev: %.3E mean(abs(derdiff/dS2lim)): %.3E \n",mean(derdiff_good),std(derdiff_good), mean(abs.(derdiff_good ./dSlim_s2good)))
    print(txt)
    return dS1dp_23, dS1lim23, reshape(dS2dp_23, (Nf, Nf)), reshape(dS2lim23, (Nf, Nf))
end

function derivtestS12(Nx)
    eps = 1e-4
    fhash = fink_filter_hash(1, 8, nx=Nx, pc=1, wd=1)
    im = rand(Float64, Nx, Nx).*0.1
    im[6,6]=1.0
    im1 = copy(im)
    im0 = copy(im)
    im1[2,3] += eps/2
    im0[2,3] -= eps/2
    dS12dp = Deriv_Utils_New.wst_S12_deriv(im, fhash)

    der0=DHC_compute(im0, fhash, doS2=false, doS20=false, doS12=true, norm=false)
    der1=DHC_compute(im1, fhash, doS2=false, doS20=false, doS12=true, norm=false)
    dS = (der1-der0) ./ eps

    Nf = length(fhash["filt_index"])
    i0 = 3+Nf
    blarg = dS12dp[2, 3, :, :]
    diff = dS[i0:end]-reshape(blarg,Nf*Nf)
    println(dS[i0:end])
    println("and")
    println(blarg)
    println("stdev: ",std(diff))
    println(diff)
    println(size(diff), size(dS), size(dS12dp))
    txt = @sprintf("Mean: %.3E stdev: %.3E mean(abs(diff/dS12dp)): %.3E, mean(abs(dS12dp/dS)) : %.3E \n",mean(diff),std(diff), mean(abs.(diff ./dS12dp[2, 3, :, :][:])), mean(abs.(dS12dp[2, 3, :, :][:]/dS[i0:end])))
    print(txt)

    #plot(dS[3:end])
    #plot!(blarg[2,3,:])
    #plot(diff)
    return reshape(dS12dp[2, 3, :, :], (Nf, Nf)), reshape(dS[i0:end], (Nf, Nf))
end

#Deriv sum test funcs
function dS1S2sum_combtest(fhash) #Tested for non-iso
    (Nf, )    = size(fhash["filt_index"])
    Nx        = fhash["npix"]
    im = rand(Nx,Nx)
    mywts = rand(Nf + (Nf*Nf))
    println(size(mywts))

    # Use new faster code
    sum1 = Deriv_Utils_New.wst_S1S2_derivsum_comb(im, fhash, mywts)
    wts1 = reshape(mywts[1:Nf], (Nf, 1))
    dwS1 = reshape(Deriv_Utils_New.wst_S1_deriv(im, fhash), (Nx^2, Nf)) * wts1

    # Compare to established code
    dS1S2dp = Deriv_Utils_New.wst_S1S2_deriv(im, fhash)
    dS1dp = dS1S2dp[:, :, 1:Nf]
    dS2dp = dS1S2dp[:, :, Nf+1:end]
    term1 = reshape(dS1dp, (Nx^2, Nf)) * reshape(mywts[1:Nf], (Nf, 1))
    sum2 = zeros(Float64, Nx, Nx)
    for i=1:Nf*Nf sum2 += (dS2dp[:,:,i].*mywts[i+Nf]) end
    dScomb_brute = term1 + reshape(sum2, (Nx^2, 1))

    println("tot abs mean", mean(abs.(dScomb_brute -sum1)))
    println("tot Stdev: ",std(dScomb_brute-sum1))
    println("Ratio: ", mean(dScomb_brute./sum1))
    println("S1")
    println("Abs mean", mean(abs.(term1 - dwS1)))
    println("Abs mean", std(term1 - dwS1))
    println("Ratio: ", mean(term1./dwS1))
    println("S2 separate")
    sum1 = Deriv_Utils_New.wst_S2_deriv_sum(im, fhash, reshape(mywts[Nf+1:end], (Nf, Nf)))
    sum2 = zeros(Float64, Nx, Nx)
    for i=1:Nf*Nf sum2 += (dS2dp[:,:,i].*mywts[Nf+i]) end
    sum2 = reshape(sum2, (Nx^2, 1))
    println("abs mean", mean(abs.(sum2 -sum1)))
    println("Stdev: ",std(sum1-sum2))
    println("Ratio: ", mean(sum2./sum1))
    return
end

function dS20sum_test(fhash) #MODIFIED from sandbox.jl
    (Nf, )    = size(fhash["filt_index"])
    Nx        = fhash["npix"]
    im = rand(Nx,Nx)
    mywts = rand(Nf*Nf)

    # this is symmetrical, but not because S20 is symmetrical!
    wtgrid = reshape(mywts, Nf, Nf)
    wtvec  = reshape(wtgrid, Nf*Nf)

    # Use new faster code
    sum1 = Deriv_Utils_New.wst_S20_deriv_sum(im, fhash, wtgrid, FFTthreads=1)

    # Compare to established code
    dS20 = reshape(Deriv_Utils_New.wst_S20_deriv(im, fhash, FFTthreads=1),Nx,Nx,Nf*Nf)
    sum2 = zeros(Float64, Nx, Nx)
    for i=1:Nf*Nf sum2 += (dS20[:,:,i].*wtvec[i]) end
    sum2 = reshape(sum2, (Nx^2, 1))
    println("Abs mean", mean(abs.(sum2 - sum1)))
    println("Abs mean", std(sum2 - sum1))
    println("Ratio: ", mean(sum2./sum1))
    return
end

function imgreconS2test(Nx, pixmask; norm=norm, bias=0.0, std_added=nothing, std_sim=nothing, sim_smoothed=false, sim_smoothedval=0.8, input_smoothed=false, coeff_choice="S2", invcov="Diagonal", mask=nothing, lambda=nothing)
    #=
    Experiments handled by this wrapper:
    S2 | S12 | S20
    Create covariance using simulated noise + img true (Choice of smoothing or not smoothing)
    Init that is smoothed or not smoothed
    Choice of Diagonal | Diagonal+Eps | Full | Full+Eps:
        +Eps adds an epsilon to the diagonal of the covariance matrix to reduce Cond No
        Diagonal considers just the standard deviation of the simulated weights and Full adds an epsilon
    Pixmask: Keeps certain pixels fixed to true values in the image (Not tested using this wrapper yet). False=> Floating pixels.
    Coeffmask: Of length 2+Nf+Nf^2. Choice of coefficients to be optimized. True=>Optimized
    =#

    img = readdust(Nx)
    if std_added==nothing
        std_added=std(img)
    else

    end
    if std_sim==nothing
        std_sim =std(img)
    else

    end
    noise = reshape(rand(Normal(bias, std_added), Nx^2), (Nx, Nx))
    init = img+noise
    if input_smoothed init = imfilter(init, Kernel.gaussian(0.8)) end

    if coeff_choice=="S2"
        #=
        Default assumes coefficients to be optimized are: All S1 and S2 coefficients larger than 1e-5
        =#
        epsilon=1e-5
        fhash = fink_filter_hash(1, 8, nx=Nx, pc=1, wd=1, Omega=true)
        (Nf,) = size(fhash["filt_index"])

        #Uniform Added Noise
        #noise = rand(Nx, Nx).*(2*high) .- high
        #init = copy(img)+noise
        #init[1, 2] += 10.0
        #init[23, 5] -= 25.0
        s2targ = DHC_compute(img, fhash, doS2=true, norm=norm, iso=false)
        #Mask
        if mask==nothing
            mask = (s2targ .>= epsilon)
            mask[1]=false
            mask[2]=false
        else
            if size(mask)[1]!=(2+Nf+Nf^2) error("Wrong length mask") end
        end
        sig2, s2cov = Data_Utils.S2_whitenoiseweights(img, fhash, Nsam=10, loc=bias, sig=std_sim, smooth=sim_smoothed, smoothval = sim_smoothedval, coeff_mask=mask) #S2_uniweights(img, fhash, Nsam=10, high=high, iso=false, norm=norm)
        sig2 = sig2[:]
        println("NF=", size(fhash["filt_index"]), "Sel Coeffs", count((i->(i==true)), mask), size(mask), " Size s2targ, type ", typeof(s2targ[mask]), " Size s2w ", typeof(sig2))

        if invcov=="Diagonal"
            s2icov = invert_covmat(sig2)
        elseif invcov=="Diagonal+Eps"
            s2icov = invert_covmat(sig2, 1e-5)
        elseif invcov=="Full+Eps"
            s2icov = invert_covmat(s2cov, 1e-5)
        else#s2icov
            s2icov = invert_covmat(s2cov)
        end
        recon_img = Deriv_Utils_New.image_recon_derivsum(init, fhash, Float64.(s2targ[mask]), s2icov, pixmask, "S2", coeff_mask=mask, optim_settings=Dict([("iterations", 1000), ("norm", norm)]))
        #recon_img = Deriv_Utils_New.img_reconfunc_coeffmask(init, fhash, s2targ[mask], s2icov[mask, mask], pixmask, Dict([("iterations", 500)]), coeff_mask=mask)
    elseif coeff_choice=="S20"
        #=
        Default assumes coefficients to be optimized are: All S20 coefficients
        =#
        fhash = fink_filter_hash(1, 8, nx=Nx, pc=1, wd=1, Omega=true)
        (Nf,) = size(fhash["filt_index"])
        #img = (img .- mean(img))./std(img) #normalized

        #Uniform Added Noise
        #noise = rand(Nx, Nx).*(2*high) .- high
        #init = copy(img)+noise
        #init[1, 2] += 10.0
        #init[23, 5] -= 25.0
        #Mask
        if mask==nothing
            mask = trues(2+Nf+Nf^2)
            mask[1:2+Nf] .= false
        else
            if size(mask)[1]!=(2+Nf+Nf^2) error("Wrong length mask") end
        end

        s2targ = DHC_compute(img, fhash, doS2=false, doS20=true, norm=norm, iso=false)
        println("NF=", size(fhash["filt_index"]))

        sig2, s2cov = Data_Utils.S20_whitenoiseweights(img, fhash, Nsam=10, loc=bias, sig=std_sim, smooth=sim_smoothed, smoothval = sim_smoothedval, coeff_mask=mask, coeff_choice="S20") #S2_uniweights(img, fhash, Nsam=10, high=high, iso=false, norm=norm)
        sig2 = sig2[:]
        println("NF=", size(fhash["filt_index"]), "Sel Coeffs", count((i->(i==true)), mask), size(mask), " Size s2targ, type ", typeof(s2targ[mask]), " Size s2w ", typeof(sig2))

        if invcov=="Diagonal"
            s2icov = invert_covmat(sig2)
        elseif invcov=="Diagonal+Eps"
            s2icov = invert_covmat(sig2, 1e-10)
        elseif invcov=="Full+Eps"
            s2icov = invert_covmat(s2cov, 1e-10)
        else#s2icov
            s2icov = invert_covmat(s2cov)
        end
        if lambda==nothing
            recon_img = Deriv_Utils_New.image_recon_derivsum(init, fhash, Float64.(s2targ[mask]), s2icov, pixmask, "S20", coeff_mask=mask, optim_settings=Dict([("iterations", 1000), ("norm", norm)]))
        else
            recon_img = Deriv_Utils_New.image_recon_derivsum_regularized(init, fhash, Float64.(s2targ[mask]), s2icov, pixmask, "S20", coeff_mask=mask, optim_settings=Dict([("iterations", 1000), ("norm", norm)]), lambda=0.001)
        end

    else
        if coeff_choice!="S12" error("Not implemented") end
        img = readdust(Nx)
        fhash = fink_filter_hash(1, 8, nx=Nx, pc=1, wd=1, Omega=true)
        (Nf,) = size(fhash["filt_index"])

        #Uniform Added Noise
        #noise = rand(Nx, Nx).*(2*high) .- high
        #init = copy(img)+noise
        #init[1, 2] += 10.0
        #init[23, 5] -= 25.0
        s2targ = DHC_compute(img, fhash, doS2=false, doS20=false, doS12=true, norm=norm, iso=false)
        #Mask
        if mask==nothing
            mask = (s2targ .>= 1e-5)
            mask[1:Nf+2] .= false
        else

        end
        sig2, s2cov = Data_Utils.S20_whitenoiseweights(img, fhash, Nsam=10, loc=bias, sig=std_sim, smooth=sim_smoothed, smoothval = sim_smoothedval, coeff_choice="S12", coeff_mask=mask) #S2_uniweights(img, fhash, Nsam=10, high=high, iso=false, norm=norm)
        sig2 = sig2[:]
        println("NF=", size(fhash["filt_index"]), "Sel Coeffs", count((i->(i==true)), mask), size(mask), " Size s2targ, type ", typeof(s2targ[mask]), " Size s2w ", typeof(sig2))

        if invcov=="Diagonal"
            s2icov = invert_covmat(sig2)
        elseif invcov=="Diagonal+Eps"
            s2icov = invert_covmat(sig2, 1e-10)
        elseif invcov=="Full+Eps"
            s2icov = invert_covmat(s2cov, 1e-10)
        else#s2icov
            s2icov = invert_covmat(s2cov)
        end
        recon_img = Deriv_Utils_New.image_recon_derivsum(init, fhash, Float64.(s2targ[mask]), s2icov, pixmask, "S12", coeff_mask=mask, optim_settings=Dict([("iterations", 1000), ("norm", norm)]))

    end
    return img, init, recon_img
end

function logimgreconS2test(Nx, pixmask; norm=norm, std_added=nothing, std_sim=nothing, sim_smoothed=false, sim_smoothedval=0.8, input_smoothed=false, coeff_choice="S2", invcov="Diagonal", mask=nothing)
    #=
    Uses the coefficients of the LOG of the image
    Experiments handled by this wrapper:
    S2 | S12 | S20
    Create covariance using simulated noise + img true (Choice of smoothing or not smoothing)
    Init that is smoothed or not smoothed
    Choice of Diagonal | Diagonal+Eps | Full | Full+Eps:
        +Eps adds an epsilon to the diagonal of the covariance matrix to reduce Cond No
        Diagonal considers just the standard deviation of the simulated weights and Full adds an epsilon
    Pixmask: Keeps certain pixels fixed to true values in the image (Not tested using this wrapper yet). False=> Floating pixels.
    Coeffmask: Of length 2+Nf+Nf^2. Choice of coefficients to be optimized. True=>Optimized
    =#
    img = readdust(Nx)
    oriimg = copy(img)
    logimg = log.(img)
    if std_added==nothing
        std_added = std(oriimg)
    else

    end

    epsilon=1e-5
    fhash = fink_filter_hash(1, 8, nx=Nx, pc=1, wd=1, Omega=true)
    (Nf,) = size(fhash["filt_index"])

    #Std Normal
    noise = reshape(rand(Normal(0.0, std_added), Nx^2), (Nx, Nx))
    init = oriimg+noise

    if std_sim==nothing
        std_sim = std(oriimg)
    else

    end

    if input_smoothed init = imfilter(init, Kernel.gaussian(0.8)) end

    if coeff_choice=="S2"
        s2targ = DHC_compute(logimg, fhash, doS2=true, norm=norm, iso=false)
        #Mask: Default assumes coefficients to be optimized are: All S1 and S2 coefficients larger than 1e-5
        if mask==nothing
            mask = (s2targ .>= epsilon)
            mask[1]=false
            mask[2]=false
        else
            if size(mask)[1]!=(2+Nf+Nf^2) error("Wrong length mask") end
        end
    elseif coeff_choice=="S20"
        s2targ = DHC_compute(logimg, fhash, doS20=true, doS2=false, norm=norm, iso=false)
        #Mask: Default assumes only S20
        if mask==nothing
            mask = trues(2+Nf+Nf^2)
            mask[1:2+Nf] .= false
        else
            if size(mask)[1]!=(2+Nf+Nf^2) error("Wrong length mask") end
        end
    else
        if coeff_choice=="S12"
            s2targ = DHC_compute(logimg, fhash, doS20=false, doS12=true, doS2=false, norm=norm, iso=false)
            #Mask
            if mask==nothing
                mask = (s2targ .>= 1e-5)
                mask[1:Nf+2] .= false
            else
                if size(mask)[1]!=(2+Nf+Nf^2) error("Wrong length mask") end
            end
        end
    end

    sig2, s2cov = whitenoiseweights_forlog(oriimg, fhash, coeff_choice, Nsam=10, loc=0.0, sig=std_sim, smooth=sim_smoothed, smoothval = sim_smoothedval, coeff_mask=mask) #S2_uniweights(img, fhash, Nsam=10, high=high, iso=false, norm=norm)
    sig2 = sig2[:]
    println("NF=", size(fhash["filt_index"]), "Sel Coeffs", count((i->(i==true)), mask), size(mask), " Size s2targ, type ", typeof(s2targ[mask]), " Size s2w ", typeof(sig2))

    if invcov=="Diagonal"
        s2icov = invert_covmat(sig2)
    elseif invcov=="Diagonal+Eps"
        s2icov = invert_covmat(sig2, 1e-5)
    elseif invcov=="Full+Eps"
        s2icov = invert_covmat(s2cov, 1e-5)
    else#s2icov
        s2icov = invert_covmat(s2cov)
    end
    logrecon_img = Deriv_Utils_New.image_recon_derivsum(log.(init), fhash, Float64.(s2targ[mask]), s2icov, pixmask, coeff_choice, coeff_mask=mask, optim_settings=Dict([("iterations", 1000), ("norm", norm)]))
    recon_img = exp.(logrecon_img)
    return oriimg, init, recon_img
end#redundant?

function dbn_experiment(true_img, noisy_init, dbn_covar, pixmask; norm=norm, bias=0.0, std_sim=nothing, sim_smoothed=false, sim_smoothedval=0.8, input_smoothed=false, coeff_choice="S2", invcov="Diagonal", mask=nothing, lambda=nothing)
    #=
    Experiments handled by this wrapper:
    S20
    Regularized vs Non regularized
    Use externally calculated covariance
    Choice of Diagonal | Diagonal+Eps | Full | Full+Eps:
        +Eps adds an epsilon to the diagonal of the covariance matrix to reduce Cond No
        Diagonal considers just the standard deviation of the simulated weights and Full adds an epsilon
    Pixmask: Keeps certain pixels fixed to true values in the image (Not tested using this wrapper yet). False=> Floating pixels.
    Coeffmask: Of length 2+Nf+Nf^2. Choice of coefficients to be optimized. True=>Optimized
    =#

    (Nx, ) = size(true_img)

    #=
    Default assumes coefficients to be optimized are: All S20 coefficients
    =#
    fhash = fink_filter_hash(1, 8, nx=Nx, pc=1, wd=1, Omega=true)
    (Nf,) = size(fhash["filt_index"])
    #img = (img .- mean(img))./std(img) #normalized

    #Uniform Added Noise
    #noise = rand(Nx, Nx).*(2*high) .- high
    #init = copy(img)+noise
    #init[1, 2] += 10.0
    #init[23, 5] -= 25.0
    #Mask
    if mask==nothing
        mask = trues(2+Nf+Nf^2)
        mask[1:2+Nf] .= false
    else
        if size(mask)[1]!=(2+Nf+Nf^2) error("Wrong length mask") end
    end

    s2targ = DHC_compute(img, fhash, doS2=false, doS20=true, norm=false, iso=false)
    println("NF=", size(fhash["filt_index"]))
    println("NF=", size(fhash["filt_index"]), "Sel Coeffs", count((i->(i==true)), mask), size(mask), " Size s2targ, type ", typeof(s2targ[mask]), " Size s2w ", typeof(sig2))

    if invcov=="Diagonal"
        s2icov = invert_covmat(Diagonal(dbn_covar))
    elseif invcov=="Diagonal+Eps"
        s2icov = invert_covmat(Diagonal(dbn_covar), 1e-10)
    elseif invcov=="Full+Eps"
        s2icov = invert_covmat(dbn_covar, 1e-10)
    else#s2icov
        s2icov = invert_covmat(dbn_covar)
    end
    if lambda==nothing
        recon_img = Deriv_Utils_New.image_recon_derivsum(init, fhash, Float64.(s2targ[mask]), s2icov, pixmask, "S20", coeff_mask=mask, optim_settings=Dict([("iterations", 1000), ("norm", norm)]))
    else
        recon_img = Deriv_Utils_New.image_recon_derivsum_regularized(init, fhash, Float64.(s2targ[mask]), s2icov, pixmask, "S20", coeff_mask=mask, optim_settings=Dict([("iterations", 1000), ("norm", norm)]), lambda=0.001)
    end
    return img, init, recon_img
end#redundant?

#Calling functions############################################
ds1code, ds1lim = derivtestS1(16)
ds1code, ds1lim, ds2code, ds2lim = derivtestS1S2(16)
ds12code, ds12lim = derivtestS12(16)
#All above work correctly for norm=false

Nx=16
fhash = fink_filter_hash(1, 8, nx=Nx, pc=1, wd=1)
dS1S2sum_combtest(fhash)

dS20sum_test(fhash)

#Image Recon tests############################################
#S2 with no pixmask, coeff_mask for large selection
img, init, recon = imgreconS2test(64, falses((64, 64)), norm=false, sim_smoothed=false, input_smoothed=false, coeff_choice="S2", invcov="Diagonal")
Visualization.plot_synth_QA(img, init, recon, fink_filter_hash(1, 8, nx=32, pc=1, wd=1, Omega=true), fname="scratch_NM/TestPlots/WhiteNoiseS2tests/S20_passed_chisq.png")
mean(abs.(init - img)), mean(abs.(recon - img))
#Not working
#Compare all quantities with the equivalentimgrecon code that didnt use derivsum and only uses coeff_mask
#That fails with DHC_compute too. Try DHC_compute_old

#Examining S2w with new dhc-compute
Nx=32
img = readdust(Nx)
fhash = fink_filter_hash(1, 8, nx=Nx, pc=1, wd=1, Omega=true)
(Nf,) = size(fhash["filt_index"])
s2w, s2cov = S2_weights(img, fhash, 10, iso=false, norm=true)
mask = s2w .>= 1e-5
mask[1]=false
mask[2]=false
println("Max/Min coeffmasked S2w", maximum(s2w[mask]), minimum(s2w[mask]))
println("Max/Min coeffmasked S2icov", maximum(s2icov[mask, mask]), minimum(s2icov[mask, mask]))

Nx=64
img = readdust(Nx)
fhash = fink_filter_hash(1, 8, nx=Nx, pc=1, wd=1, Omega=true)
(Nf,) = size(fhash["filt_index"])
s2targ = DHC_compute(img, fhash, doS2=false, doS20=false, doS12=true, norm=false, iso=false)
std_sim=std(img)
s2w, s2icov = Data_Utils.S20_whitenoiseweights(img, fhash, Nsam=10, loc=0.0, sig=std_sim, smooth=false, smoothval =0.8, coeff_choice="S12") #S2_uniweights(img, fhash, Nsam=10, high=high, iso=false, norm=norm)
mask = s2targ .>= 1e-5
mask[1:Nf+2] .= false
println("Nf=", size(fhash["filt_index"]))
println("NumSelCoeffs=", count((i->(i==1)), mask))

#NOTE: #1 Choosing what the covariance matrix should be for different coefficients
#S2
Nx=64
img = readdust(Nx)
fhash = fink_filter_hash(1, 8, nx=Nx, pc=1, wd=1, Omega=true)
(Nf,) = size(fhash["filt_index"])
s2targ = DHC_compute(img, fhash, doS2=true, doS20=false, doS12=false, norm=false, iso=false)
#Incl only S1, S2 that are larger than 1e-5
mask = s2targ .>= 1e-5
mask[1:2] .= false
std_sim=std(img)
s2w, s2cov = Data_Utils.S2_whitenoiseweights(img, fhash, Nsam=10, loc=0.0, sig=std_sim, smooth=false, smoothval =0.8, coeff_mask=mask) #S2_uniweights(img, fhash, Nsam=10, high=high, iso=false, norm=norm)
Data_Utils.invert_covmat(s2w) #1.41E18, #1.26E-20
Data_Utils.invert_covmat(s2w, 1e-10) #39k, #1.26E-20
Data_Utils.invert_covmat(s2w, 1e-5) #4.92, #1.26E-20 Try this!

Data_Utils.invert_covmat(s2cov) #2e33, #1e6 CAN'T DO THIS!
Data_Utils.invert_covmat(s2cov, 1e-10) #800k
Data_Utils.invert_covmat(s2cov, 1e-5) #9.24, 1e-19 Try this!

#Incl all S1, S2
mask = trues(2+Nf+Nf^2)
mask[1:2] .= false
std_sim=std(img)
s2w, s2cov = Data_Utils.S2_whitenoiseweights(img, fhash, Nsam=10, loc=0.0, sig=std_sim, smooth=false, smoothval =0.8, coeff_mask=mask) #S2_uniweights(img, fhash, Nsam=10, high=high, iso=false, norm=norm)
Data_Utils.invert_covmat(s2w) #1.41E28, #1.26E-20
Data_Utils.invert_covmat(s2w, 1e-5) #6.2, 1-20
Data_Utils.invert_covmat(s2cov) #2e43, 1e11 CAN'T DO THIS!
Data_Utils.invert_covmat(s2cov, 1e-5) #12, 1e-19 Try this!

#S20
Nx=64
img = readdust(Nx)
fhash = fink_filter_hash(1, 8, nx=Nx, pc=1, wd=1, Omega=true)
(Nf,) = size(fhash["filt_index"])
s2targ = DHC_compute(img, fhash, doS2=false, doS20=true, doS12=false, norm=false, iso=false)
#Including all S20
mask=trues(2+Nf+Nf^2)
mask[1:Nf+2] .= false
std_sim=std(img)
s2w, s2cov = Data_Utils.S20_whitenoiseweights(img, fhash, Nsam=10, loc=0.0, sig=std_sim, smooth=false, smoothval =0.8, coeff_choice="S20", coeff_mask=mask) #S2_uniweights(img, fhash, Nsam=10, high=high, iso=false, norm=norm)
Data_Utils.invert_covmat(s2w) #513, #1.26E-20
Data_Utils.invert_covmat(s2w, 1e-5) #6.7, #1.26E-20 Try this!

Data_Utils.invert_covmat(s2cov) #Breaks
Data_Utils.invert_covmat(s2cov, 1e-10) #5e6, 1e-10
Data_Utils.invert_covmat(s2cov, 1e-5) #55, 1e-17 Choosing this!

#Incl only S20 that are larger than 1e-5
mask = s2targ .>= 1e-5
mask[1:Nf+2] .= false
std_sim=std(img)
s2w, s2cov = Data_Utils.S20_whitenoiseweights(img, fhash, Nsam=10, loc=0.0, sig=std_sim, smooth=false, smoothval =0.8, coeff_choice="S20", coeff_mask=mask) #S2_uniweights(img, fhash, Nsam=10, high=high, iso=false, norm=norm)
Data_Utils.invert_covmat(s2w) #361
Data_Utils.invert_covmat(s2w, 1e-5) #6
Data_Utils.invert_covmat(s2cov) #e7, Breaks
Data_Utils.invert_covmat(s2cov, 1e-5) #51

#S12
Nx=64
img = readdust(Nx)
fhash = fink_filter_hash(1, 8, nx=Nx, pc=1, wd=1, Omega=true)
(Nf,) = size(fhash["filt_index"])
s2targ = DHC_compute(img, fhash, doS2=false, doS20=false, doS12=true, norm=false, iso=false)
#Including all S12:
mask=trues(2+Nf+Nf^2)
mask[1:Nf+2] .= false
std_sim=std(img)
s2w, s2cov = Data_Utils.S20_whitenoiseweights(img, fhash, Nsam=10, loc=0.0, sig=std_sim, smooth=false, smoothval =0.8, coeff_choice="S12", coeff_mask=mask) #S2_uniweights(img, fhash, Nsam=10, high=high, iso=false, norm=norm)
Data_Utils.invert_covmat(s2w, 1e-5) #e7
Data_Utils.invert_covmat(s2cov, 1e-5) #e8

#Including only S12 larger than some eps
mask = s2targ .>=1e-5
mask[1:Nf+2] .= false
std_sim=std(img)
s2w, s2cov = Data_Utils.S20_whitenoiseweights(img, fhash, Nsam=10, loc=0.0, sig=std_sim, smooth=false, smoothval =0.8, coeff_choice="S12", coeff_mask=mask) #S2_uniweights(img, fhash, Nsam=10, high=high, iso=false, norm=norm)
Data_Utils.invert_covmat(s2w, 1e-5) #e7
Data_Utils.invert_covmat(s2cov, 1e-5) #e8
#appear to be some coefficients which aren't too small O(0.6) but still have 0 std
findall(s2w .<1e-5)
s2targ[mask][findall(s2w .<1e-5)]
#Exclude all below 1e0
mask = s2targ .>=1 #This does not exclude any of the diagonal elements (S1)
mask[1:Nf+2] .= false
std_sim=std(img)
s2w, s2cov = Data_Utils.S20_whitenoiseweights(img, fhash, Nsam=10, loc=0.0, sig=std_sim, smooth=false, smoothval =0.8, coeff_choice="S12", coeff_mask=mask) #S2_uniweights(img, fhash, Nsam=10, high=high, iso=false, norm=norm)
Data_Utils.invert_covmat(s2w, 1e-5) #100k Try this!
Data_Utils.invert_covmat(s2cov, 1e-5) #e8, Num Err=1e-8

#NOTE: Trying all the variants with reasonable condition numbers
#S2: Diagonal, All S1, S2 coeffs greater than eps=1e-5: 22->7 Works well
img, init, recon = imgreconS2test(64, falses((64, 64)), norm=false, sim_smoothed=false, input_smoothed=false, coeff_choice="S2", invcov="Diagonal+Eps")
Visualization.plot_synth_QA(img, init, recon, fname="scratch_NM/TestPlots/CoeffCombinations_WhiteNoise_Nosmooth/S2_Diag_allgreaterthaneps.png")
mean(abs.(init - img)), mean(abs.(recon - img))

#S2: FullCov, All S1, S2 coeffs greater than eps=1e-5: Breaks
img, init, recon = imgreconS2test(64, falses((64, 64)), norm=false, sim_smoothed=false, input_smoothed=false, coeff_choice="S2", invcov="Full+Eps")
Visualization.plot_synth_QA(img, init, recon, fname="scratch_NM/TestPlots/CoeffCombinations_WhiteNoise_Nosmooth/S2_Cov_allgreaterthaneps.png")
mean(abs.(init - img)), mean(abs.(recon - img))

#S2: Incl All S1+S2, adding an eps to the super small coeffs
#Diagonal mat,  22->8 Works well
Nx=64
img = readdust(Nx)
fhash = fink_filter_hash(1, 8, nx=Nx, pc=1, wd=1, Omega=true)
(Nf,) = size(fhash["filt_index"])
s2targ = DHC_compute(img, fhash, doS2=true, doS20=false, doS12=false, norm=false, iso=false)
mask = trues(2+Nf+Nf^2)
mask[1:2] .= false
img, init, recon = imgreconS2test(64, falses((64, 64)), norm=false, sim_smoothed=false, input_smoothed=false, coeff_choice="S2", invcov="Diagonal+Eps",mask=mask)
Visualization.plot_synth_QA(img, init, recon, fname="scratch_NM/TestPlots/CoeffCombinations_WhiteNoise_Nosmooth/S2_DiagEps_allincl.png")
mean(abs.(init - img)), mean(abs.(recon - img))

#Cov, 22->8 Works well
img, init, recon = imgreconS2test(64, falses((64, 64)), norm=false, sim_smoothed=false, input_smoothed=false, coeff_choice="S2", invcov="Full+Eps",mask=mask)
Visualization.plot_synth_QA(img, init, recon, fname="scratch_NM/TestPlots/CoeffCombinations_WhiteNoise_Nosmooth/S2_CovEps_allincl.png")
mean(abs.(init - img)), mean(abs.(recon - img))

#S20: All S20 only incl: Best so far.
img, init, recon = imgreconS2test(64, falses((64, 64)), norm=false, sim_smoothed=false, input_smoothed=false, coeff_choice="S20", invcov="Diagonal+Eps")
Visualization.plot_synth_QA(img, init, recon, fname="scratch_NM/TestPlots/CoeffCombinations_WhiteNoise_Nosmooth/S20_DiagEps_allincl.png")
mean(abs.(init - img)), mean(abs.(recon - img))

#S20: All S20 only incl: Breaks
img, init, recon = imgreconS2test(64, falses((64, 64)), norm=false, sim_smoothed=false, input_smoothed=false, coeff_choice="S20", invcov="Full+Eps")
Visualization.plot_synth_QA(img, init, recon, fname="scratch_NM/TestPlots/CoeffCombinations_WhiteNoise_Nosmooth/S20_FullEps_allincl.png")
mean(abs.(init - img)), mean(abs.(recon - img))

#S12:

#NOTE: With Smoothing
#S2
#Cov, 22->8 Works well
mask = trues(2+Nf+Nf^2)
mask[1:2] .= false
img, init, recon = imgreconS2test(64, falses((64, 64)), norm=false, sim_smoothed=true, input_smoothed=true, coeff_choice="S2", invcov="Full+Eps",mask=mask)
Visualization.plot_synth_QA(img, init, recon, fname="scratch_NM/TestPlots/CoeffCombinations_WhiteNoise_Smoothed/S2_CovEps_allincl.png")
mean(abs.(init - img)), mean(abs.(recon - img))

#S20: All S20 only incl: Best so far. All S20 incl.
img, init, recon = imgreconS2test(64, falses((64, 64)), norm=false, sim_smoothed=true, input_smoothed=true, coeff_choice="S20", invcov="Diagonal+Eps")
Visualization.plot_synth_QA(img, init, recon, fname="scratch_NM/TestPlots/CoeffCombinations_WhiteNoise_Smoothed/S20_DiagEps_allincl.png")
mean(abs.(init - img)), mean(abs.(recon - img))

#S12: All S12 only incl
Nx=64
img = readdust(Nx)
mean(img)
fhash = fink_filter_hash(1, 8, nx=Nx, pc=1, wd=1, Omega=true)
(Nf,) = size(fhash["filt_index"])
s12targ = DHC_compute(img, fhash, doS2=false, doS20=false, doS12=true, norm=false, iso=false)
mask = s2targ .>=1 #This does not exclude any of the diagonal elements (S1)
mask[1:Nf+2] .= false
img, init, recon = imgreconS2test(64, falses((64, 64)), norm=false, sim_smoothed=true, input_smoothed=true, coeff_choice="S12", invcov="Diagonal+Eps", mask=mask)
Visualization.plot_synth_QA(img, init, recon, fname="scratch_NM/TestPlots/CoeffCombinations_WhiteNoise_Smoothed/S12_DiagEps_allincl.png")
mean(abs.(init - img)), mean(abs.(recon - img))

#NOTE: With the log of the image's coefficients
img = readdust(64)
img, init, recon = logimgreconS2test(64, falses((64, 64)), norm=false, std_added=std(img)/5.0, std_sim=std(img)/5.0, sim_smoothed=false, input_smoothed=false, coeff_choice="S20", invcov="Diagonal+Eps")
Visualization.plot_synth_QA(img, init, recon, fname="scratch_NM/TestPlots/CoeffCombinations_WhiteNoise_Nosmooth/S20_LogCoeff_DiagEps_allincl.png")
mean(abs.(init - img)), mean(abs.(recon - img))

img, init, recon = imgreconS2test(64, falses((64, 64)), norm=false, sim_smoothed=false, std_added=std(img)/5.0, std_sim=std(img)/5.0, input_smoothed=false, coeff_choice="S20", invcov="Diagonal+Eps")
Visualization.plot_synth_QA(img, init, recon, fname="scratch_NM/TestPlots/CoeffCombinations_WhiteNoise_Nosmooth/S20_DiagEps_allincl_comparison-with-log.png")
mean(abs.(init - img)), mean(abs.(recon - img))

#NOTE: Can it correct a bias?
img = readdust(64)
img, init, recon = imgreconS2test(64, falses((64, 64)), norm=false, bias=mean(img)/10.0, sim_smoothed=false, input_smoothed=false, coeff_choice="S20", invcov="Diagonal+Eps")
mean(abs.(init - img)), mean(abs.(recon - img))
Visualization.plot_diffscales([img, init, recon, recon - img], ["Ground Truth", "Init (N(0.1<I_gt>, std(I_gt))", "Reconstruction", "Reconstruction - Ground Truth"], fname="scratch_NM/TestPlots/CoeffCombinations_WhiteNoise_Nosmooth/S20_DiagEps_allincl_BiasedWhiteNoise")

#NOTE: What effect does adding a regularizing term have?
img, init, recon = imgreconS2test(64, falses((64, 64)), norm=false, sim_smoothed=false, input_smoothed=false, coeff_choice="S20", invcov="Diagonal+Eps")
Visualization.plot_synth_QA(img, init, recon, fname="scratch_NM/3-10_Plots/NoReg_1ksamps.png")
mean(abs.(init - img)), mean(abs.(recon - img))

img, init, recon = imgreconS2test(64, falses((64, 64)), norm=false, sim_smoothed=false, input_smoothed=false, coeff_choice="S20", invcov="Diagonal+Eps", lambda=0.1)
Visualization.plot_synth_QA(img, init, recon, fname="scratch_NM/3-10_Plots/Reg_1ksamps_0-1.png")
mean(abs.(init - img)), mean(abs.(recon - img))

#NOTE: Using the Gaussian covariance of fits
im, covar = mldust(64)
im, covar = mldust(64, logbool=False)
true_img = im[:, :, 1]
Nx=64
noise = rand(Normal(0.0, std(true_img)), Nx, Nx)
init = true_img + noise
filter_hash = fink_filter_hash(1, 8, nx=Nx, pc=1, wd=1, Omega=true)
Ncovsamp = 10000
(Nf, ) = size(filter_hash["filt_index"])
s20_dbn = zeros(Float64, Ncovsamp, 2+Nf+Nf^2)
for idx=1:Ncovsamp
    s20_dbn[idx, :] = DHC_compute(im[:, :, idx], filter_hash, doS2=false, doS20=true, norm=false, iso=false)
end

s_targ_mean = mean(s20_dbn, dims=1)
scov  = (s20_dbn .- s_targ_mean)' * (s20_dbn .- s_targ_mean) ./(Ncovsamp-1)

sfddbn_experiment(im[:, :, 1], init, scov, falses(Nx, Nx), norm=false,) #complete

mean(im[:, :, 1])
#####################################################################################################
#Using NEW Reconstruction wrapper and SFD dbn

img, covar = mldust(64, logbool=false)

true_img = im[:, :, 1]
Nx=64
noise = rand(Normal(0.0, std(true_img)), Nx, Nx)
init = true_img + noise
filter_hash = fink_filter_hash(1, 8, nx=Nx, pc=1, wd=1, Omega=true)
(Nf, ) = size(filter_hash["filt_index"])
coeff_mask = falses(2+Nf+Nf^2)
coeff_mask[Nf+3:end] .= true

dhc_args = Dict(:doS2=>false, :doS12=>false, :doS20=>true, :apodize=>false)
white_noise_args = Dict(:loc=>0.0, :sig=>1.0, :Nsam=>1000, :iso=>false, :norm=>false, :smooth=>false, :smoothval=>0.8) #only if you're using noise based covar
optim_settings = Dict([("iterations", 1000), ("norm", false)])
recon_settings = Dict([("target_type", "ground truth"),("covar_type", "white_noise"), ("log", false), ("GaussianLoss", true), ("Invcov_matrix", "Diagonal+Eps"),
    ("Covar_type", "sfd_dbn"), ("apd", "non_apd"), ("optim_settings", optim_settings), ("lambda", (std(true_img).^(-2))./100), ("white_noise_args", white_nois)]) #Add constraints
#"white_noise_args": args for white_noise sims if GaussianLoss and Covar_type='white_noise'. :loc...

s2mean, s2icov = meancov_generator(true_img, filter_hash, dhc_args, coeff_mask, recon_settings)
recon_settings["s_targ_mean"] = s2mean
recon_settings["s_invcov"] = s2icov

recon_img = reconstruction_wrapper(true_img, init, filter_hash, dhc_args, coeff_mask, recon_settings)
#Visualization.plot_synth_QA(true_img, init, recon_img, fname="scratch_NM/NewWrapper/sfd_dbntest.png")
Visualization.plot_diffscales([true_img, init, recon_img, recon_img - true_img], ["GT", "Init", "Reconstruction", "Residual"], fname="scratch_NM/NewWrapper/sfd_dbntest_std-by100.png")

@benchmark dhc = DHC_compute_apd(true_img, filter_hash; Dict(:doS2=>false, :doS12=>false, :doS20=>true, :apodize=>false)...)
@benchmark dhc_apd = DHC_compute_apd(true_img, filter_hash; Dict(:doS2=>false, :doS12=>false, :doS20=>true, :apodize=>true)...)

fname = "scratch_NM/data/dust10000.fits"
f = FITS(fname, "r")
big = read(f[1])
big64 = imresize(big, 64, 64, 10000)
img_old, covar = mldust(64, logbool=false)
mean(big64[:, :, 1])
bigf64 = imresize(Float64.(big), 64, 64, 10000)

(_,__,Nslice) = size(big)
println(Nslice, " slices")
size(big)
minimum(big64)
minimum(im)
heatmap(big[:, :, 1])
heatmap(im[:, :, 1])
##############################################################################################################################

#Using NEW Reconstruction wrapper and SFD dbn

im, covar = mldust(, logbool=false)
true_img = im[:, :, 1]
Nx=64
noise = rand(Normal(0.0, std(true_img)), Nx, Nx)
init = imfilter(true_img + noise, Kernel.gaussian(0.8))
filter_hash = fink_filter_hash(1, 8, nx=Nx, pc=1, wd=1, Omega=true)
(Nf, ) = size(filter_hash["filt_index"])
coeff_mask = falses(2+Nf+Nf^2)
coeff_mask[Nf+3:end] .= true

dhc_args = Dict(:doS2=>false, :doS12=>false, :doS20=>true, :apodize=>true)
white_noise_args = Dict(:loc=>0.0, :sig=>1.0, :Nsam=>1000, :iso=>false, :norm=>false, :smooth=>false, :smoothval=>0.8) #only if you're using noise based covar
optim_settings = Dict([("iterations", 100), ("norm", false)])
recon_settings = Dict([("target_type", "ground truth"),("covar_type", "white_noise"), ("log", false), ("GaussianLoss", true), ("Invcov_matrix", "Diagonal+Eps"),
 ("apd", "non_apd"), ("optim_settings", optim_settings), ("lambda", 0.0), ("white_noise_args", white_noise_args)]) #Add constraints
#"white_noise_args": args for white_noise sims if GaussianLoss and Covar_type='white_noise'. :loc...

s2mean, s2icov = meancov_generator(true_img, filter_hash, dhc_args, coeff_mask, recon_settings)
recon_settings["s_targ_mean"] = s2mean
recon_settings["s_invcov"] = s2icov

recon_img = reconstruction_wrapper(true_img, init, filter_hash, dhc_args, coeff_mask, recon_settings)
Visualization.plot_diffscales([true_img, init, recon_img, recon_img - true_img], ["GT", "Init", "Reconstruction", "Residual"], fname="scratch_NM/NewWrapper/sfd_gtwhitenoisecov_lam0_smoothinit.png")

#Old image##############################################################################
Nx=64
im = readdust(Nx)
true_img = im

noise = rand(Normal(0.0, std(true_img)), Nx, Nx)
init = imfilter(true_img + noise, Kernel.gaussian(0.8))
filter_hash = fink_filter_hash(1, 8, nx=Nx, pc=1, wd=1, Omega=true)
(Nf, ) = size(filter_hash["filt_index"])
coeff_mask = falses(2+Nf+Nf^2)
coeff_mask[Nf+3:end] .= true

dhc_args = Dict(:doS2=>false, :doS12=>false, :doS20=>true, :apodize=>false)
white_noise_args = Dict(:loc=>0.0, :sig=>1.0, :Nsam=>1000, :iso=>false, :norm=>false, :smooth=>false, :smoothval=>0.8) #only if you're using noise based covar
optim_settings = Dict([("iterations", 1000), ("norm", false)])
recon_settings = Dict([("target_type", "ground truth"),("covar_type", "white_noise"), ("log", false), ("GaussianLoss", true), ("Invcov_matrix", "Diagonal+Eps"),
 ("apd", "non_apd"), ("optim_settings", optim_settings), ("lambda", 0.0), ("white_noise_args", white_noise_args)]) #Add constraints
#"white_noise_args": args for white_noise sims if GaussianLoss and Covar_type='white_noise'. :loc...

s2mean, s2icov = meancov_generator(true_img, filter_hash, dhc_args, coeff_mask, recon_settings)
recon_settings["s_targ_mean"] = s2mean
recon_settings["s_invcov"] = s2icov

recon_img = reconstruction_wrapper(true_img, init, filter_hash, dhc_args, coeff_mask, recon_settings)
Visualization.plot_diffscales([true_img, init, recon_img, recon_img - true_img], ["GT", "Init", "Reconstruction", "Residual"], fname="scratch_NM/NewWrapper/sfd_gtwhitenoisecov_lam0_smoothinit_apd.png")

##############################################################################################################################
using StatsBase

image = readdust(256)
apdimg_fixed = apodizer(image)
apdimg = apodizer_ori(image, 1, 1, 255)
mean(image)
mean(apdimg_fixed)

heatmap(image)
heatmap(apdimg)

image = readdust(64)
apdimg = apodizer(image, 1, 1, 63)

heatmap(image)
heatmap(apdimg)

big64

apdimg_new = apodizer_new(image)
heatmap(apdimg_new)
mean(image)

mean(apdimg_new)
mean(apdimg)

#Apodization func check
image = readdust(64)
Nx=64
Amat = wind_2d(Nx)
datad_w = fweights(Amat);#why semi-col?#BUG: Replace 256 with im_size? #<A>
meanVal = mean(image,datad_w) #<AF>
temp2d_a = (image.-meanVal).*wind_2d(Nx).+meanVal #<AF>
mean(Amat.* (image .- meanVal)) #E-14

##After modifying dS20 derivsum

filter_hash = fink_filter_hash(1, 8, nx=64, Omega=true)

dS20sum_test(filter_hash)

max_stars_repo_name: andrew-saydjari/DHC

max_issues_repo_name: andrew-saydjari/DHC

max_forks_repo_name: andrew-saydjari/DHC


text:
#Demo 6: Smith chart for publication (labels etc removed - will be in caption).
#-------------------------------------------------------------------------------
using InspectDR
using Colors

#==Input
===============================================================================#

#Constants
#-------------------------------------------------------------------------------
black = RGB24(0, 0, 0)
white = RGB24(1, 1, 1)
red = RGB24(1, 0, 0)
green = RGB24(0, 1, 0)
blue = RGB24(0, 0, 1)

#Input data
#-------------------------------------------------------------------------------
f = collect(50e6:50e6:10e9)
C1 = 5e-12; GC1 = 100
C2 = .2e-12; RC2 = 10
L = 1.5e-9; RL = 15

#==Equations
===============================================================================#
Γ(Z; Zref::Real=50.0) = (Z - Zref) ./ (Z + Zref)

#Calculations
#-------------------------------------------------------------------------------
ω = 2π*f; j = im
YC1 = (j*ω*C1)+1/GC1
ZC1 = 1./YC1
ZC2 = 1./(j*ω*C2)+RC2
ZL = j*ω*L+RL

#==Generate plot
===============================================================================#
mplot = InspectDR.Multiplot(title="Smith Chart (Publication Layout)")

plot = add(mplot, InspectDR.Plot2D)
	plot.axes = InspectDR.axes(:smith, :Z, ref=50)
	plot.ext_full = InspectDR.PExtents2D(xmin=-1.1,xmax=1.1,ymin=-1.1,ymax=1.1)
	#plot.layout.legend.enabled=true
	plot.layout.legend.width=110
	lyt = plot.layout

#Zero-out title areas, etc:
	lyt.htitle = 0
	lyt.waxlabel = 0
	lyt.haxlabel = 0
	lyt.wnolabels = 0
	lyt.wticklabel = 0
	lyt.hticklabel = 0

#Make data area sqare when saved:
	lyt.wdata = 500
	lyt.hdata = 500

	wfrm = add(plot, f, Γ(ZC1), id="cap: $(C1/1e-12)pF")
	wfrm.line = line(color=blue, width=3)
	wfrm = add(plot, f, Γ(ZC2), id="cap: $(C2/1e-12)pF")
	wfrm.line = line(color=red, width=3)
	wfrm = add(plot, f, Γ(ZL), id="ind: $(L/1e-9)nH")
	wfrm.line = line(color=green, width=3)

gplot = display(InspectDR.GtkDisplay(), mplot)
InspectDR.write_png("export_pubsmith.png", plot)
#InspectDR.write_svg("export_pubsmith.svg", plot)
#InspectDR.write_eps("export_pubsmith.eps", plot)
#InspectDR.write_pdf("export_pubsmith.pdf", plot)

:DONE

ext: jl


text:
@testset "Sink Display" begin
    fname = string(tempname(), ".wav")
    testbuf = SampleBuf(rand(Float32, 10000, 2) .- 0.5f0, srate)
    # set up a 2-channel Float32 stream
    stream = savestreaming_wav(fname, 2, srate, Float32)
    io = IOBuffer()
    show(io, stream)
    @test replace(String(take!(io)), " " => "") ==
    replace(
    """
    LibSndFile.SndFileSink{Float32, String}
      path: "$fname"
      channels: 2
      samplerate: 44100Hz
      position: 0 of 0 frames
                0.00 of 0.00 seconds"""
    , " " => "")
    write(stream, testbuf)
    show(io, stream)
    @test replace(String(take!(io)), " " => "") ==
    replace(
    """
    LibSndFile.SndFileSink{Float32, String}
      path: "$fname"
      channels: 2
      samplerate: 44100Hz
      position: 10000 of 10000 frames
                0.23 of 0.23 seconds"""
    , " " => "")
end

@testset "Source Display" begin
    fname = string(tempname(), ".wav")
    testbuf = SampleBuf(rand(Float32, 10000, 2) .- 0.5f0, srate)
    save_wav(fname, testbuf)
    # set up a 2-channel Float32 stream
    stream = loadstreaming_wav(fname)
    io = IOBuffer()
    show(io, stream)
    @test replace(String(take!(io)), " " => "") ==
    replace(
    """
    LibSndFile.SndFileSource{Float32, String}
      path: "$fname"
      channels: 2
      samplerate: 44100Hz
      position: 0 of 10000 frames
                0.00 of 0.23 seconds"""
    , " " => "")
    read(stream, 5000)
    show(io, stream)
    @test replace(String(take!(io)), " " => "") ==
    replace(
    """
    LibSndFile.SndFileSource{Float32, String}
      path: "$fname"
      channels: 2
      samplerate: 44100Hz
      position: 5000 of 10000 frames
                0.11 of 0.23 seconds"""
    , " " => "")
end


text:
function test_bijector_reals(
    b::Bijector{0},
    x_true::Real,
    y_true::Real,
    logjac_true::Real;
    isequal = true,
    tol = 1e-6
)
    ib = @inferred inv(b)
    y = @inferred b(x_true)
    logjac = @inferred logabsdetjac(b, x_true)
    ilogjac = @inferred logabsdetjac(ib, y_true)
    res = @inferred forward(b, x_true)

    # If `isequal` is false, then we use the computed `y`,
    # but if it's true, we use the true `y`.
    ires = isequal ? @inferred(forward(inv(b), y_true)) : @inferred(forward(inv(b), y))

    # Always want the following to hold
    @test ires.rv ≈ x_true atol=tol
    @test ires.logabsdetjac ≈ -logjac atol=tol

    if isequal
        @test y ≈ y_true atol=tol                      # forward
        @test (@inferred ib(y_true)) ≈ x_true atol=tol # inverse
        @test logjac ≈ logjac_true                     # logjac forward
        @test res.rv ≈ y_true atol=tol                 # forward using `forward`
        @test res.logabsdetjac ≈ logjac_true atol=tol  # logjac using `forward`
    else
        @test y ≠ y_true                          # forward
        @test (@inferred ib(y)) ≈ x_true atol=tol # inverse
        @test logjac ≠ logjac_true                # logjac forward
        @test res.rv ≠ y_true                     # forward using `forward`
        @test res.logabsdetjac ≠ logjac_true      # logjac using `forward`
    end
end

function test_bijector_arrays(
    b::Bijector,
    xs_true::AbstractArray{<:Real},
    ys_true::AbstractArray{<:Real},
    logjacs_true::Union{Real, AbstractArray{<:Real}};
    isequal = true,
    tol = 1e-6
)
    ib = @inferred inv(b)
    ys = @inferred b(xs_true)
    logjacs = @inferred logabsdetjac(b, xs_true)
    res = @inferred forward(b, xs_true)
    # If `isequal` is false, then we use the computed `y`,
    # but if it's true, we use the true `y`.
    ires = isequal ? @inferred(forward(inv(b), ys_true)) : @inferred(forward(inv(b), ys))

    # always want the following to hold
    @test ys isa typeof(ys_true)
    @test logjacs isa typeof(logjacs_true)
    @test mean(abs, ires.rv - xs_true) ≤ tol
    @test mean(abs, ires.logabsdetjac + logjacs) ≤ tol

    if isequal
        @test mean(abs, ys - ys_true) ≤ tol                     # forward
        @test mean(abs, (ib(ys_true)) - xs_true) ≤ tol          # inverse
        @test mean(abs, logjacs - logjacs_true) ≤ tol           # logjac forward
        @test mean(abs, res.rv - ys_true) ≤ tol                 # forward using `forward`
        @test mean(abs, res.logabsdetjac - logjacs_true) ≤ tol  # logjac `forward`
        @test mean(abs, ires.logabsdetjac + logjacs_true) ≤ tol # inverse logjac `forward`
    else
        # Don't want the following to be equal to their "true" values
        @test mean(abs, ys - ys_true) > tol           # forward
        @test mean(abs, logjacs - logjacs_true) > tol # logjac forward
        @test mean(abs, res.rv - ys_true) > tol       # forward using `forward`

        # Still want the following to be equal to the COMPUTED values
        @test mean(abs, ib(ys) - xs_true) ≤ tol           # inverse
        @test mean(abs, res.logabsdetjac - logjacs) ≤ tol # logjac forward using `forward`
    end
end

"""
    test_bijector(b::Bijector, xs::Array; kwargs...)
    test_bijector(b::Bijector, xs::Array, ys::Array, logjacs::Array; kwargs...)

Tests the bijector `b` on the inputs `xs` against the, optionally, provided `ys`
and `logjacs`.

If `ys` and `logjacs` are NOT provided, `isequal` will be set to `false` and
`ys` and `logjacs` will be set to `zeros`. These `ys` and `logjacs` will be
treated as "counter-examples", i.e. values NOT to match.

# Arguments
- `b::Bijector`: the bijector to test
- `xs`: inputs (has to be several!!!)(has to be several, i.e. a batch!!!) to test
- `ys`: outputs (has to be several, i.e. a batch!!!) to test against
- `logjacs`: `logabsdetjac` outputs (has to be several!!!)(has to be several, i.e.
  a batch!!!) to test against

# Keywords
- `isequal = true`: if `false`, it will be assumed that the given values are
  provided as "counter-examples" in the sense that the inputs `xs` should NOT map
  to the given outputs. This is useful in cases where one might not know the expected
  output, but still wants to test that the evaluation, etc. works.
  This is set to `true` by default if `ys` and `logjacs` are not provided.
- `tol = 1e-6`: the absolute tolerance used for the checks. This is also used to check
  arrays where we check that the L1-norm is sufficiently small.
"""
function test_bijector(b::Bijector{0}, xs::AbstractVector{<:Real})
    return test_bijector(b, xs, zeros(length(xs)), zeros(length(xs)); isequal = false)
end

function test_bijector(b::Bijector{1}, xs::AbstractMatrix{<:Real})
    return test_bijector(b, xs, zeros(size(xs)), zeros(size(xs, 2)); isequal = false)
end

function test_bijector(
    b::Bijector{0},
    xs_true::AbstractVector{<:Real},
    ys_true::AbstractVector{<:Real},
    logjacs_true::AbstractVector{<:Real};
    kwargs...
)
    ib = inv(b)

    # Batch
    test_bijector_arrays(b, xs_true, ys_true, logjacs_true; kwargs...)

    # Test `logabsdetjac` against jacobians
    test_logabsdetjac(b, xs_true)
    test_logabsdetjac(b, ys_true)

    for (x_true, y_true, logjac_true) in zip(xs_true, ys_true, logjacs_true)
        test_bijector_reals(b, x_true, y_true, logjac_true; kwargs...)

        # Test AD
        if isclosedform(b)
            test_ad(x -> b(first(x)), [x_true, ])
        end

        if isclosedform(ib)
            y = b(x_true)
            test_ad(x -> ib(first(x)), [y, ])
        end

        test_ad(x -> logabsdetjac(b, first(x)), [x_true, ])
    end
end

function test_bijector(
    b::Bijector{1},
    xs_true::AbstractMatrix{<:Real},
    ys_true::AbstractMatrix{<:Real},
    logjacs_true::AbstractVector{<:Real};
    kwargs...
)
    ib = inv(b)

    # Batch
    test_bijector_arrays(b, xs_true, ys_true, logjacs_true; kwargs...)

    # Test `logabsdetjac` against jacobians
    test_logabsdetjac(b, xs_true)
    test_logabsdetjac(b, ys_true)

    for (x_true, y_true, logjac_true) in zip(eachcol(xs_true), eachcol(ys_true), logjacs_true)
        # HACK: collect to avoid dealing with sub-arrays and thus allowing us to compare the
        # type of the computed output to the "true" output.
        test_bijector_arrays(b, collect(x_true), collect(y_true), logjac_true; kwargs...)

        # Test AD
        if isclosedform(b)
            test_ad(x -> sum(b(x)), collect(x_true))
        end
        if isclosedform(ib)
            y = b(x_true)
            test_ad(x -> sum(ib(x)), y)
        end

        test_ad(x -> logabsdetjac(b, x), x_true)
    end
end

function test_logabsdetjac(b::Bijector{1}, xs::AbstractMatrix; tol=1e-6)
    if isclosedform(b)
        logjac_ad = [logabsdet(ForwardDiff.jacobian(b, x))[1] for x in eachcol(xs)]
        @test mean(logabsdetjac(b, xs) - logjac_ad) ≤ tol
    end
end

function test_logabsdetjac(b::Bijector{0}, xs::AbstractVector; tol=1e-6)
    if isclosedform(b)
        logjac_ad = [log(abs(ForwardDiff.derivative(b, x))) for x in xs]
        @test mean(logabsdetjac(b, xs) - logjac_ad) ≤ tol
    end
end


text:
### A Pluto.jl notebook ###
# v0.19.8

using Markdown
using InteractiveUtils

# This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error).
macro bind(def, element)
    quote
        local iv = try Base.loaded_modules[Base.PkgId(Base.UUID("6e696c72-6542-2067-7265-42206c756150"), "AbstractPlutoDingetjes")].Bonds.initial_value catch; b -> missing; end
        local el = $(esc(element))
        global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : iv(el)
        el
    end
end

# ╔═╡ 2efa05e7-b130-4a95-9eab-d8bb9d16a29b
# ╠═╡ show_logs = false
begin
	using Logging
	global_logger(NullLogger())
	display("")
end

# ╔═╡ c3e429e4-e7e9-4db6-852c-906630f909a4
# ╠═╡ show_logs = false
#Set-up packages
begin

	using DataFrames, HTTP, CSV, Dates, Plots, PlutoUI, Printf, LaTeXStrings, HypertextLiteral, XLSX

	gr();
	Plots.GRBackend()

	#Define html elements
	nbsp = html"&nbsp" #non-breaking space
	vspace = html"""

"""
	br = html"
"

	#Sets the width of cells, caps the cell width by 90% of screen width
	#(setting overwritten by cell below)
	# @bind screenWidth @htl("""
	#

	#
	#

	# """)

	# cellWidth= min(1000, screenWidth*0.9)
	# @htl("""
	#
	# """)

	#Sets the width of the cells
	#begin
	#	html"""


text:
using FunctionalCollections
using Random
_pvec(x::PersistentVector) = x
_pvec(x::AbstractArray) = pvec(x)

const newid = let rng = MersenneTwister()
    function (prefix)
        Base.depwarn("newid(prefix) is deprecated. Scopes don't need this anymore. To get the scope's id, use scopeid(scope)", :newid)
        string(prefix, '-', uuid4(rng))
    end
end
# b can be array of pairs / kwargs etc.
function recmerge!(a, b, f=recmerge!)
    for (k, v) in b
        if isa(v, AbstractDict) && haskey(a, k) && isa(a[k], AbstractDict)
            a[k] = f(a[k], v)
        else
            a[k] = v
        end
    end
    a
end

recmerge(a, b) = recmerge!(Dict{Any, Any}(a), b, recmerge)

function kebab2camel(str)
    strs = split(str, "-")
    titlefrom2 = titlecase.(strs[2:end])
    join([strs[1]; titlefrom2], "")
end

function camel2kebab(str)
    lowstr = lowercase(str)
    lowercaseidxs = Vector{Char}(lowstr) .== Vector{Char}(str)
    sb = IOBuffer()
    for (i,c) in enumerate(lowercase(str))
        !lowercaseidxs[i] && print(sb, '-')
        print(sb, c)
    end
    String(take!(sb))
end

"""
Escape characters for a "safe" representation of JSON.

In particular, we escape '/' characters to avoid the presence of "</" (and
especially "</script>") which cause the browser to break out of the current
<script /> tag.
"""
function escape_json(s::String)
    # Replace all "/" with "\/"'s.
    # This prevents the browser from interpreting "</" as a close tag; since
    # everything within the string is JSON, any appearances of "/" should be
    # within strings and when the JSON is parsed, the "\/"'s will be interpreted
    # as just normal "/"'s.
    return replace(s, "/" => "\\/")
end

escape_json(x::Any) = escape_json(JSON.json(x))


text:
using Plots: plot, savefig, xlabel!, ylabel!

function plotresults(mean_rewards, eps_rewards)
    # Plotting rewards per episode
    x_axis = 1:length(mean_rewards)
    @assert length(eps_rewards) == length(mean_rewards) # Just to be safe...
    y_axis = hcat(eps_rewards, mean_rewards)
    p = plot(x_axis, y_axis,
        title="Monte Carlo Policy Gradient Method on CartPole-v0",
        label=["episodic reward", "10-episode average of reward"],
        lw=2, legend=:topleft)
    xlabel!("Episode #")
    ylabel!("Rewards")
    savefig(p, "Reward.png")
end


text:
# utility
# -------

const ERROR_COLOR = :light_red
const NOERROR_COLOR = :light_green
const RAIL_COLORS = (
    # preserve yellow for future performance linting
    45, # light blue
    123, # light cyan
    150, # ???
    215, # orange
    231, # white
)
const N_RAILS = length(RAIL_COLORS)
const LEFT_ROOF  = "═════ "
const RIGHT_ROOF = " ═════"
const HEADER_COLOR = :reverse
const ERROR_SIG_COLOR = :bold
const TYPE_ANNOTATION_COLOR = :light_cyan

pluralize(n::Integer, one::AbstractString, more::AbstractString = string(one, 's')) =
    string(n, ' ', isone(n) ? one : more)

printlnstyled(args...; kwarg...) = printstyled(args..., '\n'; kwarg...)

function print_rails(io, depth)
    for i = 1:depth
        color = RAIL_COLORS[i%N_RAILS+1]
        printstyled(io, '│'; color)
    end
end

function with_bufferring(f, args...)
    buf = IOBuffer()
    io = IOContext(buf, args...)
    f(io)
    return String(take!(buf))
end

# we may need something like this for stdlibs as well ?
expandbasepath(filename) =
    return normpath(Sys.BINDIR::String, Base.DATAROOTDIR, "julia", "base", filename)
function tofullpath(filename::AbstractString)
    path = abspath(filename)
    return isfile(path) ? path : expandbasepath(filename)
end

# toplevel
# --------

function print_reports(io::IO,
                       reports::AbstractVector{<:ToplevelErrorReport},
                       @nospecialize(postprocess = identity);
                       print_toplevel_sucess::Bool = false,
                       color::Bool = get(io, :color, false),
                       fullpath::Bool = false,
                       # dummy kwargs so that kwargs for other functions can be passed on
                       __kwargs...)
    if isempty(reports)
        if print_toplevel_sucess
            printlnstyled(io, "No toplevel errors !"; color = NOERROR_COLOR)
        end
        return false
    end

    arg = :color => color
    color = ERROR_COLOR

    with_bufferring(arg) do io
        s = string(pluralize(length(reports), "toplevel error"), " found")
        printlnstyled(io, LEFT_ROOF, s, RIGHT_ROOF; color = HEADER_COLOR)

        rail = with_bufferring(arg) do io
            printstyled(io, "│ "; color)
        end

        for report in reports
            s = string("┌ @ ", (fullpath ? tofullpath : identity)(string(report.file)), ':', report.line, ' ')
            printlnstyled(io, s; color)

            errlines = with_bufferring(arg) do io
                print_report(io, report)
            end |> strip
            println(io, join(string.(rail, split(errlines, '\n')), '\n'))

            s = string('└', '─'^(length(s)-1))
            printlnstyled(io, s; color)
        end
    end |> postprocess |> Fix1(print, io)

    return true
end

# don't show stacktrace for syntax errors
print_report(io, report::SyntaxErrorReport) = showerror(io, report.err)
function print_report(io, report::RecursiveIncludeErrorReport)
    printstyled(io, "ERROR: "; bold = true, color = ERROR_COLOR)
    println(io, "recursive `include` call detected:")
    println(io, " ⚈ duplicated file: ", report.duplicated_file)
    println(io, " ⚈  included files: ", join(report.files, ' '))
end
# TODO: add context information, i.e. during macroexpansion, defining something
print_report(io, report::ActualErrorWrapped) = showerror(io, report.err, report.st)

# inference
# ---------

function print_reports(io::IO,
                       reports::AbstractVector{<:InferenceErrorReport},
                       @nospecialize(postprocess = identity);
                       print_inference_sucess::Bool = true,
                       color::Bool = get(io, :color, false),
                       fullpath::Bool = false,
                       annotate_types::Bool = false,
                       # dummy kwargs so that kwargs for other functions can be passed on
                       __kwargs...)
    reports = uniquify_reports(reports)

    if isempty(reports)
        if print_inference_sucess
            printlnstyled(io, "No errors !"; color = NOERROR_COLOR)
        end
        return false
    end

    with_bufferring(:color => color) do io
        s = string(pluralize(length(reports), "possible error"), " found")
        printlnstyled(io, LEFT_ROOF, s, RIGHT_ROOF; color = HEADER_COLOR)

        # don't duplicated virtual stack frames for reports from the same toplevel frame
        toplevel_linfo_hash = hash(:dummy)
        wrote_linfos = Set{UInt64}()
        for report in reports
            new_toplevel_linfo_hash = hash(first(report.st))
            if toplevel_linfo_hash != new_toplevel_linfo_hash
                toplevel_linfo_hash = new_toplevel_linfo_hash
                wrote_linfos = Set{UInt64}()
            end
            print_report(io, report, wrote_linfos; fullpath, annotate_types)
        end
    end |> postprocess |> Fix1(print, io)

    return true
end

# this is a dirty fix for duplicated reports
# TODO: analysis can efficiently avoids collecting duplicated reports ?
function uniquify_reports(reports)
    return unique(reports) do report; return #= uniquify keys =# (
        # report itself
        report.msg,
        report.sig,
        # entry point
        first(report.st),
        # error point
        last(report.st),
    ); end
end

# traverse abstract call stack, print frames
function print_report(io, report::InferenceErrorReport, wrote_linfos, depth = 1; kwargs...)
    if length(report.st) == depth # error here
        return print_error_frame(io, report, depth; kwargs...)
    end

    frame = report.st[depth]

    # cache current frame info
    linfo_hash = hash(frame)
    should_print = linfo_hash ∉ wrote_linfos
    push!(wrote_linfos, linfo_hash)

    # print current frame and go into deeper
    should_print && print_frame(io, frame, depth, false; kwargs...)
    print_report(io, report, wrote_linfos, depth + 1; kwargs...)
end

function print_error_frame(io, report, depth; kwargs...)
    frame = report.st[depth]

    len = print_frame(io, frame, depth, true; kwargs...)
    print_rails(io, depth-1)
    print_error_report(io, report)

    print_rails(io, depth-1)
    printlnstyled(io, '└', '─'^len; color = ERROR_COLOR)
end

function print_frame(io, frame, depth, is_err;
                     fullpath = false,
                     kwargs...)
    print_rails(io, depth-1)

    color = is_err ? ERROR_COLOR : RAIL_COLORS[(depth)%N_RAILS+1]
    s = string("┌ @ ", (fullpath ? tofullpath : identity)(string(frame.file)), ':', frame.line)
    printstyled(io, s, ' '; color)
    print_signature(io, frame.sig; kwargs...)

    return length(s) # the length of frame info string
end

function print_signature(io, sig; kwargs...)
    for a in sig
        _print_signature(io, a; kwargs...)
    end
    println(io)
end
_print_signature(io, a::Union{AbstractChar,AbstractString};
                 annotate_types = false,
                 kwargs...) =
    printstyled(io, a; kwargs...)
function _print_signature(io, @nospecialize(typ); annotate_types = false, kwargs...)
    annotate_types || return

    printstyled(io, "::", string(typ); color = TYPE_ANNOTATION_COLOR, kwargs...)
end

# default error report printer
function print_error_report(io, report::InferenceErrorReport)
    printstyled(io, "│ ", report.msg, ": "; color = ERROR_COLOR)
    print_signature(io, report.sig;
                    annotate_types = true, # always annotate types for errored signatures
                    bold = true,
                    )
end
# those don't need explicit signatures
print_error_report(io, report::NoFieldErrorReport)       = printlnstyled(io, "│ ", report.msg; color = ERROR_COLOR)
print_error_report(io, report::LocalUndefVarErrorReport) = printlnstyled(io, "│ ", report.msg; color = ERROR_COLOR)
print_error_report(io, report::UndefKeywordErrorReport)  = printlnstyled(io, "│ ", report.msg; color = ERROR_COLOR)


text:
# ------------------------------------------------------------------
# Licensed under the MIT License. See LICENSE in the project root.
# ------------------------------------------------------------------

# -------------
# DOMAIN VIEWS
# -------------

"""
    DomainView(domain, indices)

A partial view of a `domain` containing only the elements at `indices`.
"""
struct DomainView{Dim,T,D<:Domain{Dim,T},I} <: Domain{Dim,T}
  domain::D
  inds::I
end

# -----------------
# DOMAIN INTERFACE
# -----------------

element(v::DomainView, ind::Int) = element(v.domain, v.inds[ind])

nelements(v::DomainView) = length(v.inds)

centroid(v::DomainView, ind::Int) = centroid(v.domain, v.inds[ind])

# -----------
# IO METHODS
# -----------

function Base.show(io::IO, v::DomainView)
  domain = getfield(v, :domain)
  nelms  = length(getfield(v, :inds))
  print(io, "$nelms View{$domain}")
end

# -----------
# DATA VIEWS
# -----------

# helper function for table view
viewtable(table, rows) =
  view(Tables.rows(table) |> collect, rows) |> Tables.columns

"""
    DataView(data, inds)

Return a view of `data` at indices `inds`.
"""
struct DataView{D<:Data,I} <: Data
  data::D
  inds::I
end

# ---------------
# DATA INTERFACE
# ---------------

function domain(v::DataView)
  data = getfield(v, :data)
  inds = getfield(v, :inds)
  view(domain(data), inds)
end

function values(v::DataView, rank=nothing)
  data = getfield(v, :data)
  inds = getfield(v, :inds)
  R = paramdim(domain(data))
  r = isnothing(rank) ? R : rank
  𝒯 = values(data, r)
  r == R ? viewtable(𝒯, inds) : nothing
end

function constructor(::Type{DataView{D,I}}) where {D<:Data,I}
  function ctor(domain, values)
    data = constructor(D)(domain, values)
    inds = 1:nelements(domain)
    DataView(data, inds)
  end
end

# specialize methods for performance
==(v₁::DataView, v₂::DataView) =
  getfield(v₁, :data) == getfield(v₂, :data) &&
  getfield(v₁, :inds) == getfield(v₂, :inds)

# -----------
# IO METHODS
# -----------

function Base.show(io::IO, v::DataView)
  data  = getfield(v, :data)
  nelms = length(getfield(v, :inds))
  print(io, "$nelms View{$data}")
end


text:
export getEdgeLength, getFaceArea, getCellVolume

function getEdgeLength(mesh::TetraMesh)
# get edge length

if isempty(mesh.L)

  ne = mesh.ne
  L = zeros(ne)

  for i=1:ne
    i1 = mesh.edges[i,1]
    i2 = mesh.edges[i,2]

    L[i] = getEdgeLength(mesh.Points[i1,:],  mesh.Points[i2,:])
  end

  mesh.L = L

end

return mesh.L
end

function getFaceArea(mesh::TetraMesh)
# get face size

if isempty(mesh.F)
  nf = mesh.nf
  F = zeros(nf)
  for i=1:nf
    i1 = mesh.faces[i,1]
    i2 = mesh.faces[i,2]
    i3 = mesh.faces[i,3]

    F[i] = getFaceArea(mesh.Points[i1,:], mesh.Points[i2,:], mesh.Points[i3,:])
  end

mesh.F = F
end

return mesh.F
end

function getCellVolume(mesh::TetraMesh)
# get volume

if isempty(mesh.V)

  nc = mesh.nc
  T  = mesh.Tetras
  V = zeros(nc)

  for i=1:nc
    i1 = T[i,1]
    i2 = T[i,2]
    i3 = T[i,3]
    i4 = T[i,4]

    V[i] = getCellVolume(mesh.Points[i1,:], mesh.Points[i2,:], mesh.Points[i3,:], mesh.Points[i4,:])
  end

  mesh.V = V

end

return mesh.V
end

function getEdgeLength(P1::Array{Float64,1}, P2::Array{Float64,1})
	return sqrt((P1[1] - P2[1])^2 + (P1[2] - P2[2])^2 + (P1[3] - P2[3])^2)
end

function getFaceArea(P1::Array{Float64,1}, P2::Array{Float64,1}, P3::Array{Float64,1})
	return norm(cross(P3 - P1, P3 - P2))/2
end

function getCellVolume(P1::Array{Float64,1}, P2::Array{Float64,1}, P3::Array{Float64,1}, P4::Array{Float64,1})
	return 1/6.0*abs(dot(P4 - P3, cross(P4 - P1, P4 - P2)))
end

max_stars_repo_name: JuliaInv/JTetra

max_issues_repo_name: JuliaInv/JTetra

max_forks_repo_name: JuliaInv/JTetra


text:
module AcronymGenerator

export Acronym, AcronymFamily, ScoreParams, ScoredAcronymFamily
export generate_acronyms, acronyms, score, best_acronym, generate_prefix_map

include("util.jl")
include("acronym.jl")
include("score.jl")
include("findbest.jl")
include("display.jl")

end


text:
module TORA

using NeuralNetworkAnalysis, LaTeXStrings
using NeuralNetworkAnalysis: UniformAdditivePostprocessing
import Plots, DifferentialEquations
using Plots: plot, plot!, xlims!, ylims!, savefig, font, Measures.mm

# problem

@taylorize function TORA!(dx, x, p, t)
    x₁, x₂, x₃, x₄, u = x

    aux = 0.1 * sin(x₃)
    dx[1] = x₂
    dx[2] = -x₁ + aux
    dx[3] = x₄
    dx[4] = u
    dx[5] = zero(u)
    return dx
end

controller = read_nnet_sherlock(@modelpath "" "controller_TORA_sherlock")

X₀ = Hyperrectangle(low=[0.6, -0.7, -0.4, 0.5], high=[0.7, -0.6, -0.3, 0.6])
U = ZeroSet(1)
vars_idx = Dict(:state_vars=>1:4, :control_vars=>5)
ivp = @ivp(x' = TORA!(x), dim: 5, x(0) ∈ X₀ × U)
control_postprocessing = UniformAdditivePostprocessing(-10.0)

period = 1.0
k = 20
T = k * period
T_warmup = 2 * period

prob = ControlledPlant(ivp, controller, vars_idx, period;
                       postprocessing=control_postprocessing)

safe_states = cartesian_product(BallInf(zeros(4), 2.0), Universe(1))
predicate = X -> X ⊆ safe_states

# simulation

println("simulation")
res = @timed simulate(prob, T=T, trajectories=10, include_vertices=true)
sim = res.value
print_timed(res)

# reachability analysis

alg = TMJets(abstol=1e-10, orderT=8, orderQ=3)
alg_nn = Ai2()
splitter = BoxSplitter([4, 4, 3, 5])

function benchmark(; T=T, silent::Bool=false)
    silent || println("flowpipe construction")
    res = @timed solve(prob, T=T, alg_nn=alg_nn, alg=alg, splitter=splitter)
    sol = res.value
    silent || print_timed(res)

    silent || println("property checking")
    solh = overapproximate(sol, Hyperrectangle)
    res = @timed predicate(solh)
    silent || print_timed(res)
    if res.value
        silent || println("The property is satisfied.")
    else
        silent || println("The property may be violated.")
    end

    return solh
end

benchmark(T=T_warmup, silent=true)
res = @timed benchmark()
sol = res.value

# plotting

vars = (1, 2)
fig = plot(xlab=L"x_1", ylab=L"x_2",
           legendfontsize=15,
           tickfont=font(25, "Times"),
           guidefontsize=25,
           xguidefont=font(25, "Times"),
           yguidefont=font(25, "Times"),
           bottom_margin=0mm, left_margin=0mm, right_margin=3mm, top_margin=2mm,
           size=(500, 450))
plot!(project(X₀, vars), c=:blue, lab=L"X_0", alpha=0.6)
plot!(fig, sol, vars=vars, color=:orange, lab="")
plot!(project(X₀, vars), c=:blue, lab="", alpha=0.6)
plot_simulation!(fig, sim; vars=vars, color=:black, lab="")
xlims!(-1, 1)
ylims!(-1.5, 1)
savefig("tora_x1x2_juliareach.png")

vars = (3, 4)
fig = plot(xlab=L"x_3", ylab=L"x_4",
           legendfontsize=15,
           tickfont=font(25, "Times"),
           guidefontsize=25,
           xguidefont=font(25, "Times"),
           yguidefont=font(25, "Times"),
           bottom_margin=0mm, left_margin=0mm, right_margin=0mm, top_margin=2mm,
           size=(500, 450))
# plot!(project(X₀, vars), c=:blue, lab=L"X_0", alpha=0.6)
plot!(fig, sol, vars=vars, color=:orange, lab="")
plot!(project(X₀, vars), c=:blue, lab="", alpha=0.6)
plot_simulation!(fig, sim; vars=vars, color=:black, lab="")
xlims!(-2, 2)
ylims!(-2, 2)
savefig("tora_x3x4_juliareach.png")

end
nothing


text:
using MySQL_Octo
using Base.Test

adapter = SQL.connect(host="127.0.0.1", user="test", pass="test", db="mysqltest", reconnect=true)
@test !isa(adapter.handle, Void)
u = table_from(adapter, "users")
@test repr(u) == "FROM users"
@test repr(u[*]) == "SELECT * FROM users"
@test repr(u[u.id]) == "SELECT id FROM users"

SQL.arrange(DefaultArrangement, STDOUT, u[*](:where => and(u.id == 1, u.name == "John")))

# custom arrangement
function SQL.arrange(io::IO, query::Queryable)
    write(io, QueryBuilders.colored_repr(query), "\n")
    for row in SQL.iterate(query)
        write(io, string(map(x -> isa(x, Nullable) ? x.value : x, row), "\n"))
    end
end

Octo.Text.colors[:keyword] = :bold

SQL.arrange(STDOUT, u[*](:where => and(u.id == 1, u.name == "John")))

Octo.Text.colors[:keyword] = :normal

SQL.arrange(STDOUT, u[u.name](:where => and(u.id == 1, u.name == "John")))

SQL.arrange(STDOUT, u[u.name, u.id](:where => and(u.id == 1, u.name == "John"), :limit => 2))


text:
using Discord
include("base.jl")

function reserve(client::Client,m::Message,boss_number::Int,damage::Int;over=false,attack=false)
    global management_dict
    channel_id_str = string(m.channel_id)
    if !haskey(management_dict,channel_id_str)
        return reply(client,m,"初期設定が完了しておりません．`.set`コマンドにて初期設定を完了させてください")
    end

    if !haskey(management_dict[channel_id_str],"remain_attack_channel")
        return reply(client,m,"残凸把握板の登録が完了しておりません．`.set`コマンドにて初期設定を完了させてください")
    end

    if  !haskey(management_dict[channel_id_str],"output_channel")
        return reply(client,m,"予約確認板の登録が完了しておりません．`.set`コマンドにて初期設定を完了させてください")
    end

    if  !haskey(management_dict[channel_id_str],"output_message")
        return reply(client,m,"先に`.battle`コマンドを使用して下さい")
    end

    if !over
        #.attackコマンドを使用した場合
        if attack
            nothing
        elseif management_dict[channel_id_str]["reserve"][string(boss_number)]["plan_remain_hp"] == 0
            if boss_number == management_dict[channel_id_str]["boss_number"]
                return reply(client,m,"予約を無視して戦闘を行う場合には`.attack [予定ダメージ]`を送信して下さい")
            else
                return reply(client,m,"既に予約が一杯です")
            end
        end
    end

    reserve_dict = Dict(
        "name" =>  get_display_name(m),
        "id" => m.author.id,
        "damage" => damage,
        "attack" => attack,
        "over" => over
    )

    if attack
        push!(management_dict[channel_id_str]["reserve"][string(management_dict[channel_id_str]["boss_number"])]["reserve_list"],reserve_dict)
    else
        push!(management_dict[channel_id_str]["reserve"][string(boss_number)]["reserve_list"],reserve_dict)
    end

    management_dict[channel_id_str]["reserve"][string(boss_number)]["plan_remain_hp"] = calc_plan_remain_hp(management_dict[channel_id_str]["reserve"][string(boss_number)])
    output_message_content = create_output_message(management_dict[channel_id_str])
    edit_message(client,management_dict[channel_id_str]["output_channel"],management_dict[channel_id_str]["output_message"],content=output_message_content)

    dt_now =now(tz"Asia/Tokyo") - Hour(5)
    dt_now_str = "$(month(dt_now))月$(day(dt_now))日"
    if dt_now_str != management_dict[channel_id_str]["today"]#日付が変わっていれば残凸状況を修正する
        management_dict[channel_id_str]["today"] = dt_now_str
        management_dict[channel_id_str]["remain_attack"] =initiate_remain_attack_dict(client,management_dict[channel_id_str])
        remain_attack_message_content = make_remain_attack_message_content(management_dict[channel_id_str]["remain_attack"],management_dict[channel_id_str]["today"])
        remain_attack_message = fetchval(create_message(client,management_dict[channel_id_str]["remain_attack_channel"],content=remain_attack_message_content))

        management_dict[channel_id_str]["remain_attack_message"] = remain_attack_message.id
    end

    save_management_dict(management_dict)
    create(client, Reaction, m, '👍')

end

max_stars_repo_name: hibibol/clan-battle

max_issues_repo_name: hibibol/clan-battle

max_forks_repo_name: hibibol/clan-battle


text:
using NQCDynamics
using Test
using RecursiveArrayTools: ArrayPartition
using ComponentArrays: ComponentVector

atoms = Atoms([:H])
model = NQCModels.Harmonic()
sim = Simulation(atoms, model)

positions = [randn(1, length(atoms)) for i=1:10]
velocities = [randn(1, length(atoms)) for i=1:10]
distribution = DynamicalDistribution(positions, velocities, (1, 1))
u0 = rand(distribution)
tspan = (0.0, 10.0)

@testset "run_ensemble" for reduction in (:append, :sum, :mean)
    out = run_ensemble(sim, tspan, distribution; reduction,
        output=(:position, :velocity, :hamiltonian), dt=1, trajectories=10)
    if (reduction === :sum) || (reduction === :mean)
        @test out.t == 0.0:10.0
        @test out.position isa Vector{<:AbstractMatrix}
        @test out.velocity isa Vector{<:AbstractMatrix}
        @test out.hamiltonian isa Vector{<:Number}
    else
        @test all(x -> x == 0.0:10.0, (traj.t for traj in out))
        @test all(x -> x isa Vector{<:AbstractMatrix}, (traj.position for traj in out))
        @test all(x -> x isa Vector{<:AbstractMatrix}, (traj.velocity for traj in out))
        @test all(x -> x isa Vector{<:Number}, (traj.hamiltonian for traj in out))
    end
end

@testset "run_ensemble" begin
    out = run_ensemble(sim, tspan, distribution;
        output=Ensembles.OutputFinal(), dt=1, trajectories=10, reduction=:append)
    @test out isa Vector{<:ArrayPartition}
end

@testset "run_ensemble, reduction=$reduction" for reduction ∈ (:sum, :mean)
    out = run_ensemble(sim, tspan, distribution;
        output=Ensembles.OutputFinal(), dt=1, trajectories=10, reduction=reduction, u_init=u0)
    @test out isa ComponentVector
end


text:
# ==================================
# Backward approximation
# ==================================

"""
    Backward{EM, SO, SI, IT, BT} <: AbstractApproximationModel

Backward approximation model.

### Fields

- `exp`     -- exponentiation method
- `setops`  -- set opertaions method
- `sih`     -- symmetric interval hull
- `inv`     -- (optional, default: `false`) if `true`, assume that the state matrix
               is invertible and use its inverse in the `Φ` functions
- `backend` -- (optional, default: `nothing`) used if the algorithm needs to apply
               concrete polyhedral computations

### Algorithm

The transformations are:

- ``Φ ← \\exp(Aδ)``,
- ``Ω_0 ← CH(\\mathcal{X}_0, Φ\\mathcal{X}_0 ⊕ δU(0) ⊕ E_ψ(U(0), δ) ⊕ E^-(\\mathcal{X}_0, δ))``,
- ``V(k) ← δU(k) ⊕ E_ψ(U(k), δ)``.

Here we allow ``U`` to be a sequence of time varying non-deterministic input sets.

For the definition of the sets ``E_ψ`` and ``E^-`` see [[FRE11]](@ref).
The `Forward` method uses ``E^+``.
"""
struct Backward{EM, SO, SI, IT, BT} <: AbstractApproximationModel
    exp::EM
    setops::SO
    sih::SI
    inv::IT
    backend::BT
end

# convenience constructor using symbols
function Backward(; exp=BaseExp, setops=:lazy, sih=:concrete, inv=false, backend=nothing)
    return Backward(_alias(exp), _alias(setops), Val(sih), Val(inv), backend)
end

function Base.show(io::IO, alg::Backward)
    print(io, "`Backward` approximation model with: \n")
    print(io, "    - exponentiation method: $(alg.exp) \n")
    print(io, "    - set operations method: $(alg.setops)\n")
    print(io, "    - symmetric interval hull method: $(alg.sih)\n")
    print(io, "    - invertibility assumption: $(alg.inv)")
    print(io, "    - polyhedral computations backend: $(alg.backend)")
end

Base.show(io::IO, m::MIME"text/plain", alg::Backward) = print(io, alg)

# TODO: add corresponding `discrete` methods <<<<<


text:
# leyends
# TODO fix this...
let
    for (title, colors) in [
            ("exp", exp_colors),
            ("iders", ider_colors),
            ("method", method_colors)
        ]
    p = plot(; framestyle = :none)
        scolors = sort(collect(colors); by = (p) -> string(first(p)))
        for (id, color) in scolors
            scatter!(p, [0], [0];
                thickness_scaling = 1,
                color, ms = 8, label = string(id),
                legendfontsize=10,
                # size = [300, 900],
                # legend = :left
            )
        end
        mysavefig(p, "$(title)_color_legend")
    end
end


text:
using Colors

# Test data from the DIN 6176 specification
const testdata = [
((50,  10,  10), (61.43,   9.70,   3.76)),
((50,  50,  50), (61.43,  28.64,  11.11)),
((50, -10,  10), (61.43,  -5.57,   7.03)),
((50, -50,  50), (61.43, -17.22,  21.75)),
((50, -10, -10), (61.43,  -9.70,  -3.76)),
((50, -50, -50), (61.43, -28.64, -11.11)),
((50,  10, -10), (61.43,   5.57,  -7.03)),
((50,  50, -50), (61.43,  17.22, -21.75)),
(( 0,   0,   0), ( 0,      0,      0)),
((100,  0,   0), (100,     0,      0))]

# A high error threshold has been chosen because converting from DIN99
# to CIELAB with only two decimal places of accuracy yields fairly inaccurate
# results due to the exponentiation.
const conveps = 0.05
const diffeps = 0.01
for (i, (a, b)) in enumerate(testdata)
    converted = convert(DIN99, Lab(a...))
    test = DIN99(b...)

    @assert (abs(converted.l - test.l) < conveps)
    @assert (abs(converted.a - test.a) < conveps)
    @assert (abs(converted.b - test.b) < conveps)

    converted = convert(Lab, DIN99(b...))
    test = Lab(a...)

    @assert (abs(converted.l - test.l) < conveps)
    @assert (abs(converted.a - test.a) < conveps)
    @assert (abs(converted.b - test.b) < conveps)

    # This is not a real test of the color difference metric, but at least
    # makes sure it isn't doing anything really crazy.
    metric = DE_DIN99()
    @assert (abs(colordiff(convert(DIN99, Lab(a...)), DIN99(b...), metric)) < diffeps)

end


text:
struct MatrixConvolution{T} <: AbstractArray{T,3}
    arr::Array{T,3}
end

function Base.:+(x::MatrixConvolution, y::MatrixConvolution)

    A = x.arr
    B = y.arr

    T = promote_type(eltype(A), eltype(B))
    @assert axes(A)[1:2] == axes(B)[1:2]

    kmax = max(size(A,3), size(B,3))
    C = zeros(T,size(A)[1:2]...,kmax)

    for k in axes(A,3)
        C[:,:,k] += A[:,:,k]
    end

    for k in axes(B,3)
        C[:,:,k] += B[:,:,k]
    end

    return MatrixConvolution(C)
end

Base.size(x::MatrixConvolution) = size(x.arr)
Base.getindex(x::MatrixConvolution, i::Int) = x.arr[i]
Base.IndexStyle(x::MatrixConvolution) = IndexLinear()

function Base.:*(x::MatrixConvolution, y::Matrix)

    A = x.arr

    T = promote_type(eltype(A), eltype(y))
    @assert axes(A,2) == axes(y,1)

    C = zeros(T, size(A)...)
    for k in axes(C,3)
        C[:,:,k] = A[:,:,k]*y
    end

    return MatrixConvolution(C)
end

function Base.:*(y::Matrix, x::MatrixConvolution)

    A = x.arr

    T = promote_type(eltype(A), eltype(y))
    @assert axes(A,2) == axes(y,1)

    C = zeros(T, size(A)...)
    for k in axes(C,3)
        C[:,:,k] = y*A[:,:,k]
    end

    return MatrixConvolution(C)
end

Base.:*(a::Number, x::MatrixConvolution) = MatrixConvolution(a * x.arr)
Base.:*(x::MatrixConvolution, a::Number) = MatrixConvolution(x.arr * a)
Base.:/(x::MatrixConvolution, a::Number) = MatrixConvolution(x.arr / a)

function convolve(Z::Array,x,j,k0)
    M,N,K = size(Z)
    y = similar(Z,M)
    fill!(y,0)
    for k ∈ k0 : min(j,K)
        i = j - k + 1
        y += Z[:,:,k] * x[:,i]
    end
    return y
end

convolve(x::MatrixConvolution, y::Matrix, i, j) = convolve(x.arr, y, i, j)


text:
"""
A Julia package for generating auditory stimuli.

"""
module AuditoryStimuli

using DSP
using SampledSignals
using LinearAlgebra
using Random
using Logging
using FFTW
using Unitful
using Plots
using UnitfulRecipes
using Parameters

using Unitful: AbstractQuantity, AffineQuantity, DimensionlessQuantity
import SampledSignals: nchannels, samplerate, unsafe_read!
import Plots: plot

export  bandpass_noise,
        bandpass_filter,
        amplitude_modulate,
        ITD_modulate,
        set_RMS,
        ramp_on,
        ramp_off,
        set_ITD,
        PlotSpectroTemporal,
        NoiseSource,
        CorrelatedNoiseSource,
        SinusoidSource,
        DummySampleSink,
        Amplification,
        AmplitudeModulation,
        TimeDelay,
        samplerate,
        modify,
        plot,
        interaural_coherence,
        plot_cross_correlation

# #########################################
# #####  Signal Generation
# #########################################

include("SignalGenerators/NoiseSource.jl")
include("SignalGenerators/CorrelatedNoiseSource.jl")
include("SignalGenerators/SinusoidSource.jl")
include("SignalGenerators/DummySampleSink.jl")

# #########################################
# #####  Signal Modifiers
# #########################################

include("SignalModifiers/Amplification.jl")
include("SignalModifiers/BandpassFilter.jl")
include("SignalModifiers/Modulation.jl")
include("SignalModifiers/TimeDelay.jl")

# #########################################
# #####  Signal Metrics
# #########################################

include("SignalMetrics/InterauralCoherence.jl")

# #########################################
# #####  Signal Plotting
# #########################################

include("Plotting.jl")

"""
    bandpass_noise(number_samples, number_channels, lower_bound, upper_bound, sample_rate; filter_order=14)

Generates band pass noise with specified upper and lower bounds using a butterworth filter.
"""
function bandpass_noise(number_samples::Int, number_channels::Int, lower_bound::Number, upper_bound::Number, sample_rate::Number; filter_order::Int = 14)
    bandpass_filter(randn(number_samples, number_channels), lower_bound, upper_bound, sample_rate, filter_order=filter_order)
end

# #########################################
# #####  Signal Modifiers
# #########################################

"""
    bandpass_filter(AbstractArray, lower_bound, upper_bound, sample_rate; filter_order=14)
    bandpass_filter(SampledSignal, lower_bound, upper_bound;              filter_order=14)

Signal will be filtered with bandpass butterworth filter between 'lower_bound' and `upper_bound` with filter of `filter_order`.
"""
function bandpass_filter(x::AbstractArray, lower_bound::Number, upper_bound::Number, sample_rate::Number; filter_order::Int = 14)

    responsetype = Bandpass(lower_bound, upper_bound; fs=sample_rate)
    designmethod = Butterworth(filter_order)
    filt(digitalfilter(responsetype, designmethod), x)
end

function bandpass_filter(x::SampledSignals.SampleBuf, lower_bound::typeof(1u"Hz"), upper_bound::typeof(1u"Hz"); filter_order::Int = 14)
    x.data = bandpass_filter(x.data, ustrip(lower_bound), ustrip(upper_bound), x.samplerate, filter_order=filter_order)
    x
end

"""
    amplitude_modulate(data, modulation_frequency, sample_rate; phase=π)

Amplitude modulates the signal

See [wikipedia](https://en.wikipedia.org/wiki/Amplitude_modulation)
"""
function amplitude_modulate(x::AbstractArray, modulation_frequency::Number, sample_rate::Number; phase::Number = π)

    t = 1:size(x, 1)
    t = t ./ sample_rate

    fits = mod(maximum(t), (1/modulation_frequency))
    if !(isapprox(fits, 0, atol = 1e-5) || isapprox(fits, 1/modulation_frequency, atol = 1e-5)  )
        @warn("Not a complete modulation")
    end
        # println(maximum(t))
    # println(1/modulation_frequency)
    # println(mod(maximum(t), (1/modulation_frequency)))

    M = 1 .* cos.(2 * π * modulation_frequency * t .+ phase)
    (1 .+ M) .* x;
end

function amplitude_modulate(x::SampledSignals.SampleBuf, modulation_frequency::typeof(1u"Hz"); phase::Number = π)
    amplitude_modulate(x, modulation_frequency * 1.0, phase=phase)

end

function amplitude_modulate(x::SampledSignals.SampleBuf, modulation_frequency::typeof(1.0u"Hz"); phase::Number = π)
    x.data = amplitude_modulate(x.data, ustrip(modulation_frequency), x.samplerate, phase=phase)
    x
end

"""
    ITD_modulate(data, modulation_frequency, ITD_1, ITD_2, samplerate)

Modulate an applied ITD

"""
function ITD_modulate(x::AbstractArray, modulation_frequency::Number, ITD_1::Int, ITD_2::Int, sample_rate)
    @warn "Unvalidated code"

    t = 1:size(x, 1)
    t = t ./ sample_rate

    Ti = 1 / modulation_frequency
    switches = round.(Int, collect(0:Ti:maximum(t))*sample_rate)
    switches_starts = switches[1:1:end]
    switches_stops  = switches[2:1:end]

    switch_samples = switches_stops[1] - switches_starts[1]

    for idx = 1:2:length(switches_starts)-1
        x[switches_starts[idx]+1:switches_stops[idx], 1]  = x[switches_starts[idx]+1+ITD_1:switches_stops[idx]+ITD_1, 1]  .* tukey(switch_samples, 0.01)
        x[switches_stops[idx]+1:switches_stops[idx+1], 1] = x[switches_stops[idx]+1+ITD_2:switches_stops[idx+1]+ITD_2, 1] .* tukey(switch_samples, 0.01)
    end
    return x
end

"""
    set_RMS(data, desired_rms)

Modify rms of signal to desired value

"""
function set_RMS(data::AbstractArray, desired_rms::Number)

    data / (rms(data) / desired_rms)

end

"""
    ramp_on(data, number_samples)

Apply a linear ramp to start of signal

"""
function ramp_on(data::AbstractArray, number_samples::Int)

    data[1:number_samples, :] = LinRange(0, 1, number_samples) .* data[1:number_samples, :]
    return data
end

"""
    ramp_off(data, number_samples)

Apply a linear ramp to end of signal

"""
function ramp_off(data::AbstractArray, number_samples::Int)

    data[end-number_samples+1:end, :] = LinRange(1, 0, number_samples) .* data[end-number_samples+1:end, :]
    return data
end

"""
    set_ITD(data, number_samples)

Introduce an ITD of number_samples

"""
function set_ITD(data::AbstractArray, number_samples::Int)

    abs_number_samples = abs(number_samples)

    if number_samples > 0
        data[:, 1] = [zeros(abs_number_samples, 1); data[1:end - abs_number_samples, 1]]
    elseif number_samples < 0
        data[:, 2] = [zeros(abs_number_samples, 1); data[1:end - abs_number_samples, 2]]
    end

    return data
end

end # module


text:
printstyled("  ungap!\n", color=:light_green)
printstyled("    subsample negative time gaps (Issue 29)\n", color=:light_green)

gf1 = path * "/SampleFiles/SEED/CIRIO__BHE___2017101.mseed"
opf = path * "/SampleFiles/SEED/obspy.dat"

# read merge/ungap target file
S = read_data(gf1)
C1 = deepcopy(S[1])
i1 = C1.t[2,1]
merge!(S)
ungap!(S)
C2 = S[1]
i2 = C2.t[2,1]

# read ObsPy merge output to compare
io = open(opf, "r")
X2 = Array{Int32,1}(undef, length(C2.x))
read!(io, X2)
X = map(Float32, X2)

# these should be approximately equal
@test isapprox(C2.x, X)


text:
############################################################
## joNumber
############################################################

export joNumber

# type definition
"""
joNumber type

A number type to use for jo operations with number

# TYPE PARAMETERS
- DDT::DataType : domain DataType
- RDT::DataType : range DataType

# FIELDS
- ddt::DDT : number to use when acting on vector to return domain vector
- rdt::RDT : number to use when acting on vector to return range vector

"""
struct joNumber{DDT<:Number,RDT<:Number}
    ddt::DDT
    rdt::RDT
end


text:
function CategoricalPool{S, T <: Integer, R <: Integer}(index::Vector{S},
                                           invindex::Dict{S, T},
                                           order::Vector{R},
                                           ordered::Bool=false)
    invindex = convert(Dict{S, R}, invindex)
    CategoricalPool{S, R, CategoricalValue{S, R}}(index, invindex, order, ordered)
end

@compat (::Type{CategoricalPool{T, R}}){T, R}(ordered::Bool=false) =
    CategoricalPool(T[], Dict{T, R}(), R[], ordered)
@compat (::Type{CategoricalPool{T}}){T}(ordered::Bool=false) =
    CategoricalPool(T[], Dict{T, DefaultRefType}(), DefaultRefType[], ordered)

@compat function (::Type{CategoricalPool{T, R}}){T, R}(index::Vector,
                                                       ordered::Bool=false)
    invindex = buildinvindex(index, R)
    order = Vector{R}(1:length(index))
    CategoricalPool(index, invindex, order, ordered)
end

function CategoricalPool(index::Vector, ordered::Bool=false)
    invindex = buildinvindex(index)
    order = Vector{DefaultRefType}(1:length(index))
    return CategoricalPool(index, invindex, order, ordered)
end

function CategoricalPool{S, R <: Integer}(invindex::Dict{S, R},
                                          ordered::Bool=false)
    index = buildindex(invindex)
    order = Vector{DefaultRefType}(1:length(index))
    return CategoricalPool(index, invindex, order, ordered)
end

# TODO: Add tests for this
function CategoricalPool{S, R <: Integer}(index::Vector{S},
                                          invindex::Dict{S, R},
                                          ordered::Bool=false)
    order = Vector{DefaultRefType}(1:length(index))
    return CategoricalPool(index, invindex, order, ordered)
end

function CategoricalPool{T}(index::Vector{T},
                            levels::Vector{T},
                            ordered::Bool=false)
    invindex = buildinvindex(index)
    order = buildorder(invindex, levels)
    return CategoricalPool(index, invindex, order, ordered)
end

function CategoricalPool{S, R <: Integer}(invindex::Dict{S, R},
                                          levels::Vector{S},
                                          ordered::Bool=false)
    index = buildindex(invindex)
    order = buildorder(invindex, levels)
    return CategoricalPool(index, invindex, order, ordered)
end

Base.convert(::Type{CategoricalPool}, pool::CategoricalPool) = pool
Base.convert{T}(::Type{CategoricalPool{T}}, pool::CategoricalPool{T}) = pool
Base.convert{T, R <: Integer}(::Type{CategoricalPool{T, R}}, pool::CategoricalPool{T, R}) = pool

Base.convert{S, T, R <: Integer}(::Type{CategoricalPool{S}}, pool::CategoricalPool{T, R}) =
    convert(CategoricalPool{S, R}, pool)

function Base.convert{S, R <: Integer}(::Type{CategoricalPool{S, R}}, pool::CategoricalPool)
    if length(levels(pool)) > typemax(R)
        throw(LevelsException{S, R}(levels(pool)[typemax(R)+1:end]))
    end

    indexS = convert(Vector{S}, pool.index)
    invindexS = convert(Dict{S, R}, pool.invindex)
    order = convert(Vector{R}, pool.order)
    return CategoricalPool(indexS, invindexS, order, pool.ordered)
end

function Base.show{T, R}(io::IO, pool::CategoricalPool{T, R})
    @printf(io, "%s{%s,%s}([%s])", typeof(pool).name, T, R,
                join(map(repr, levels(pool)), ","))

    pool.ordered && print(io, " with ordered levels")
end

Base.length(pool::CategoricalPool) = length(pool.index)

Base.getindex(pool::CategoricalPool, i::Integer) = pool.valindex[i]
Base.get(pool::CategoricalPool, level::Any) = pool.invindex[level]
Base.get(pool::CategoricalPool, level::Any, default::Any) = get(pool.invindex, level, default)

function Base.get!{T, R, V}(pool::CategoricalPool{T, R, V}, level)
    get!(pool.invindex, level) do
        x = convert(T, level)
        n = length(pool)
        if n >= typemax(R)
            throw(LevelsException{T, R}([level]))
        end

        i = R(n + 1)
        push!(pool.index, x)
        push!(pool.order, i)
        push!(pool.levels, x)
        push!(pool.valindex, V(i, pool))
        i
    end
end

Base.push!(pool::CategoricalPool, level) = (get!(pool, level); pool)

# TODO: optimize for multiple additions
function Base.append!(pool::CategoricalPool, levels)
    for level in levels
        push!(pool, level)
    end
    return pool
end

function Base.delete!{S, R, V}(pool::CategoricalPool{S, R, V}, levels...)
    for level in levels
        levelS = convert(S, level)
        if haskey(pool.invindex, levelS)
            ind = pool.invindex[levelS]
            delete!(pool.invindex, levelS)
            splice!(pool.index, ind)
            ord = splice!(pool.order, ind)
            splice!(pool.levels, ord)
            splice!(pool.valindex, ind)
            for i in ind:length(pool)
                pool.invindex[pool.index[i]] -= 1
                pool.valindex[i] = V(i, pool)
            end
            for i in 1:length(pool)
                pool.order[i] > ord && (pool.order[i] -= 1)
            end
        end
    end
    return pool
end

function levels!{S, R, V}(pool::CategoricalPool{S, R, V}, newlevels::Vector)
    levs = convert(Vector{S}, newlevels)
    if !allunique(levs)
        throw(ArgumentError(string("duplicated levels found in levs: ",
                                   join(unique(filter(x->sum(levs.==x)>1, levs)), ", "))))
    end

    n = length(levs)

    if n > typemax(R)
        throw(LevelsException{S, R}(setdiff(levs, levels(pool))[typemax(R)-length(levels(pool))+1:end]))
    end

    # No deletions: can preserve position of existing levels
    if issubset(pool.index, levs)
        append!(pool, setdiff(levs, pool.index))
    else
        empty!(pool.invindex)
        resize!(pool.index, n)
        resize!(pool.valindex, n)
        resize!(pool.order, n)
        resize!(pool.levels, n)
        for i in 1:n
            v = levs[i]
            pool.index[i] = v
            pool.invindex[v] = i
            pool.valindex[i] = V(i, pool)
        end
    end

    buildorder!(pool.order, pool.invindex, levs)
    for (i, x) in enumerate(pool.order)
        pool.levels[x] = pool.index[i]
    end
    return pool
end

index(pool::CategoricalPool) = pool.index
levels(pool::CategoricalPool) = pool.levels
order(pool::CategoricalPool) = pool.order

isordered(pool::CategoricalPool) = pool.ordered
ordered!(pool::CategoricalPool, ordered) = (pool.ordered = ordered; pool)

# LevelsException
function Base.showerror{T, R}(io::IO, err::LevelsException{T, R})
    levs = join(repr.(err.levels), ", ", " and ")
    print(io, "cannot store level(s) $levs since reference type $R can only hold $(typemax(R)) levels. Use the decompress function to make room for more levels.")
end


text:
using MAT: matopen
using NPZ: npzread

scale_factor_act = 1e5
scale_factor_den = 1e-2
scale_factor_dose = 1e6
T = Float32

path = "../data/petvp6/"
file = matopen(path*"petvp6-mumap-208.mat")
mumap = read(file, "mumap208")
mumap = T.(mumap)
close(file)

file = matopen(path*"petvp6-psf-208.mat")
psfs = read(file, "psf_208")
psfs = T.(psfs)
close(file)

file = matopen(path*"petvp6-density-128.mat")
density = read(file, "density")
density = max.(density, 1) # minimum density = 1
density = T.(density * scale_factor_den) # keep things in similar range
close(file)

file = matopen(path*"petvp6-doseGT-128.mat")
doseGT = read(file, "doseGT")
doseGT = T.(doseGT * scale_factor_dose)
close(file)

# file = matopen(path*"doseDVK_128.mat")
# doseDVK = read(file, "DVK")
# doseDVK = T.(doseDVK * scale_factor_dose)
# close(file)

file = matopen(path*"petvp6-xinit-128.mat")
xinit = read(file, "xinit")
xinit = T.(xinit / sum(xinit) * scale_factor_act)
close(file)

file = matopen(path*"petvp6-recon-16th.mat")
xrecon16 = read(file, "x")
xrecon16 = T.(xrecon16 / sum(xrecon16) * scale_factor_act)
close(file)

file = matopen(path*"petvp6-xtrue-128.mat")
xtrue = read(file, "xtrue")
xtrue = T.(xtrue / sum(xtrue) * scale_factor_act)
close(file)

file = matopen(path*"petvp6-yi-128.mat")
yi = read(file, "yi")
yi = T.(yi)
close(file)

file = matopen(path*"petvp6-ri-128.mat")
ri = read(file, "ri")
ri = T.(ri)
close(file)

fullkernel = npzread("../data/fullkernel.npy") # kernels are in 257x257x161