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

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 ReconFuncs

##How does full cov in the log-im case make a difference? vs Diag (Getting the high density and the other artefacts)
fname = "scratch_NM/StandardizedExp/Nx64/noisy_stdtrue/SFDTargSFDCov/log_apd_noiso/1000_triu"
gttarget = load(fname*".jld2")
#Calc 1dps log
Nx = size(gttarget["true_img"])[1]
kbins = collect(1:Nx/2.0)

true_img = gttarget["true_img"]
init = gttarget["init"]
recon_img =gttarget["recon"]

println("Mean Abs Res: Init-True = ", mean(abs.(init - true_img)), " Recon-True = ", mean(abs.(recon_img - true_img)))
println("Mean Abs Frac Res", mean(abs.((init - true_img)./true_img)), " Recon-True=", mean(abs.((recon_img - true_img)./true_img)))
println("RMSE Init-True= ", (mean((init - true_img).^2)).^0.5, " Recon-True=", (mean((recon_img - true_img).^2)).^0.5)

true_lps = calc_1dps(apodizer(log.(gttarget["true_img"])), kbins)
recon_lps = calc_1dps(apodizer(log.(gttarget["recon"])), kbins)
init_lps = calc_1dps(apodizer(log.(gttarget["init"])), kbins)
p = plot(log.(kbins), log.(true_lps), label="True")
plot!(log.(kbins), log.(recon_lps), label="Recon")
plot!(log.(kbins), log.(init_lps), label="Init")
plot!(title="DiagCov Noiso: P(k) of Log Image")
xlabel!("lnk")
ylabel!("lnP(k)")

true_ps = calc_1dps(apodizer(gttarget["true_img"]), kbins)
recon_ps = calc_1dps(apodizer(gttarget["recon"]), kbins)
init_ps = calc_1dps(apodizer(gttarget["init"]), kbins)
p = plot(log.(kbins), log.(true_ps), label="True")
plot!(log.(kbins), log.(recon_ps), label="Recon")
plot!(log.(kbins), log.(init_ps), label="Init")
plot!(title="DiagCov Noiso: P(k) of Image")
xlabel!("lnk")
ylabel!("lnP(k)")

##With the Full Cov:
fname = "scratch_NM/StandardizedExp/Nx64/noisy_stdtrue/SFDTargSFDCov/log_apd_noiso/1000_fullcov_triu"
gttarget = load(fname*".jld2")
#Calc 1dps log
Nx = size(gttarget["true_img"])[1]
kbins = collect(1:Nx/2.0)

true_img = gttarget["true_img"]
init = gttarget["init"]
recon_img =gttarget["recon"]

println("Mean Abs Res: Init-True = ", mean(abs.(init - true_img)), " Recon-True = ", mean(abs.(recon_img - true_img)))
println("Mean Abs Frac Res", mean(abs.((init - true_img)./true_img)), " Recon-True=", mean(abs.((recon_img - true_img)./true_img)))
println("RMSE Init-True= ", (mean((init - true_img).^2)).^0.5, " Recon-True=", (mean((recon_img - true_img).^2)).^0.5)

true_lps = calc_1dps(apodizer(log.(gttarget["true_img"])), kbins)
recon_lps = calc_1dps(apodizer(log.(gttarget["recon"])), kbins)
init_lps = calc_1dps(apodizer(log.(gttarget["init"])), kbins)
p = plot(log.(kbins), log.(true_lps), label="True")
plot!(log.(kbins), log.(recon_lps), label="Recon")
plot!(log.(kbins), log.(init_lps), label="Init")
plot!(title="FullCov Noiso: P(k) of Log Image")
xlabel!("lnk")
ylabel!("lnP(k)")

true_ps = calc_1dps(apodizer(gttarget["true_img"]), kbins)
recon_ps = calc_1dps(apodizer(gttarget["recon"]), kbins)
init_ps = calc_1dps(apodizer(gttarget["init"]), kbins)
p = plot(log.(kbins), log.(true_ps), label="True")
plot!(log.(kbins), log.(recon_ps), label="Recon")
plot!(log.(kbins), log.(init_ps), label="Init")
plot!(title="FullCov Noiso: P(k) of Image")
xlabel!("lnk")
ylabel!("lnP(k)")

##Why Don't Reg plots move away?? Actually it ends up closer to the target than the true image.
#Maybe loss in reg doesn't differentiate noisy from not noisy?
fname = "scratch_NM/StandardizedExp/Nx64/noisy_stdtrue/SFDTargSFDCov/reg_apd_noiso/1000_fullcov_triu"
s_targ_mean, s_targ_invcov = gttarget["dict"]["s_targ_mean"], gttarget["dict"]["s_invcov"]
dhc_args = Dict(:doS2=>false, :doS12=>false, :doS20=>true, :apodize=>true, :iso=>false)

altrue = apodizer(gttarget["true_img"])
alinit = apodizer(gttarget["init"])
alrecon = apodizer(gttarget["recon"])

coeff_mask = gttarget["coeff_mask"]
s_true = DHC_compute_wrapper(gttarget["true_img"], filter_hash, norm=false; dhc_args...)[coeff_mask]
s_init = DHC_compute_wrapper(gttarget["init"], filter_hash, norm=false; dhc_args...)[coeff_mask]
s_recon = DHC_compute_wrapper(gttarget["recon"], filter_hash, norm=false; dhc_args...)[coeff_mask]
ltruecheck = ( 0.5 .* (s_true - s_targ_mean)' * s_targ_invcov * (s_true - s_targ_mean))[1]
l2true = 0.5*gttarget["dict"]["lambda"]*sum((alinit - altrue).^2)
#7.82, 2.32

l1init = ( 0.5 .* (s_init - s_targ_mean)' * s_targ_invcov * (s_init - s_targ_mean))[1]
l2init = 0.5*gttarget["dict"]["lambda"]*sum((alinit - alinit).^2)
#11.2, 0

l1recon = ( 0.5 .* (s_recon - s_targ_mean)' * s_targ_invcov * (s_recon - s_targ_mean))[1]
l2recon = 0.5*gttarget["dict"]["lambda"]*sum((alrecon - alinit).^2)
#7.61, 3.57
ltruecheck + l2true
l1init + l2init
l1recon + l2recon

##Log
fname = "scratch_NM/StandardizedExp/Nx64/noisy_stdtrue/SFDTargSFDCov/log_apd_noiso/1000_fullcov_triu"
s_targ_mean, s_targ_invcov = gttarget["dict"]["s_targ_mean"], gttarget["dict"]["s_invcov"]
dhc_args = Dict(:doS2=>false, :doS12=>false, :doS20=>true, :apodize=>true, :iso=>false)

altrue = apodizer(log.(gttarget["true_img"]))
alinit = apodizer(log.(gttarget["init"]))
alrecon = apodizer(log.(gttarget["recon"]))

coeff_mask = gttarget["coeff_mask"]
s_true = DHC_compute_wrapper(log.(gttarget["true_img"]), filter_hash, norm=false; dhc_args...)[coeff_mask]
s_init = DHC_compute_wrapper(log.(gttarget["init"]), filter_hash, norm=false; dhc_args...)[coeff_mask]
s_recon = DHC_compute_wrapper(log.(gttarget["recon"]), filter_hash, norm=false; dhc_args...)[coeff_mask]
ltruecheck = ( 0.5 .* (s_true - s_targ_mean)' * s_targ_invcov * (s_true - s_targ_mean))[1]
l2true = 0.5*gttarget["dict"]["lambda"]*sum((alinit - altrue).^2)
#446, 84503

l1init = ( 0.5 .* (s_init - s_targ_mean)' * s_targ_invcov * (s_init - s_targ_mean))[1]
l2init = 0.5*gttarget["dict"]["lambda"]*sum((alinit - alinit).^2)
#6e9, 0

l1recon = ( 0.5 .* (s_recon - s_targ_mean)' * s_targ_invcov * (s_recon - s_targ_mean))[1]
l2recon = 0.5*gttarget["dict"]["lambda"]*sum((alrecon - alinit).^2)
#2.5e3, 6e4

ztrue_recon = (s_recon - s_true) .* diag(s_targ_invcov)
ztrue_target = (s_true - s_targ_mean) .* diag(s_targ_invcov)
zinit_target = (s_init - s_targ_mean) .* diag(s_targ_invcov)

scatter(log.(abs.(zinit_target)), log.(abs.(ztrue_recon)))
xlabel!("ln Abs Z score: Coeffs Target-Init")
ylabel!("ln Abs Z score: Coeffs Recon-True")

zscoremat = zeros(34, 34)
zscoremat[reshape(gttarget["coeff_mask"][3+34:end], (34, 34))] .= log.(abs.(ztrue_target))
heatmap(zscoremat, title="Im 1000: Ln Abs Z-Score b/w True & Target")

gttarget["fhash"]

zscoremat
#basically the J=3 scales that are problematic

filename = normpath("scratch_NM/NewWrapper/Weave", "4-7.jmd")
weave(filename, out_path = :pwd)
using Pkg
Pkg.add("IJulia")

using IJulia
IJulia.notebook(dir="/home/nayantara/Desktop/GalacticDustProject/WST/DHC/scratch_NM/NewWrapper/4-7/")

##ANALYZING Z SCORES AND COVARIANCES
#Log coeff, log im
Nx=64
dbnimg = readsfd(Nx, logbool=false)
fname = "scratch_NM/StandardizedExp/Nx64/noisy_stdtrue/SFDTargSFDCov/reg_apd_noiso/LogCoeff/1000_0-01_full_triu"
gttarget = load(fname*".jld2")
dhc_args = gttarget["dhc_args"]
coeff_mask = gttarget["coeff_mask"]
fhash = gttarget["fhash"]
dbn1 = log.(get_dbn_coeffs(dbnimg[:, :, 1:5000], fhash, dhc_args, coeff_mask = coeff_mask))
dbn2 = log.(get_dbn_coeffs(dbnimg[:, :, 5001:end], fhash, dhc_args, coeff_mask= coeff_mask))

#Mean dia, cov
Ncovsamp = 5000
s_targ_mean = mean(dbn1, dims=1)
scov  = (dbn1 .- s_targ_mean)' * (dbn1 .- s_targ_mean) ./(Ncovsamp-1)
smean1, sdiag1, scov1 = s_targ_mean, diag(scov), scov
s_targ_mean = mean(dbn2, dims=1)
scov  = (dbn2 .- s_targ_mean)' * (dbn2 .- s_targ_mean) ./(Ncovsamp-1)
smean2, sdiag2, scov2 = s_targ_mean, diag(scov), scov

scov1inv = invert_covmat(scov1, 1e-8)
scov2inv = invert_covmat(scov2, 1e-8)

diff = dbn1 .- reshape(smean1, (1, 595))
chisq1 = sum((diff * scov1inv) .* diff, dims=2)

diff = dbn2 .- reshape(smean2, (1, 595))
chisq2 = sum((diff * scov2inv) .* diff, dims=2)
zscores = (dbn1 .- s_targ_mean).* reshape(diag(scov1inv).^0.5, (1, 595))
p = histogram(chisq1[:], label="Dbn1")
histogram!(sum(zscores, dims=2)./595, label="Z_Dbn1")
histogram!(chisq2[:], label="Dbn2")
chisqsamp = rand(Distributions.Chisq(595), 5000)
histogram!(chisqsamp, label="ChiSq")

fname = "scratch_NM/StandardizedExp/Nx64/noisy_stdtrue/SFDTargSFDCov/reg_apd_noiso/LogCoeff/1000_0-01_full_triu"
gttarget = load(fname*".jld2")
filter_hash = gttarget["fhash"]
s_targ_mean, s_targ_invcov = gttarget["dict"]["fs_targ_mean"], gttarget["dict"]["fs_invcov"]
dhc_args = Dict(:doS2=>false, :doS12=>false, :doS20=>true, :apodize=>true, :iso=>false)

altrue = apodizer(gttarget["true_img"])
alinit = apodizer(gttarget["init"])
alrecon = apodizer(gttarget["recon"])

coeff_mask = gttarget["coeff_mask"]
s_true = DHC_compute_wrapper(gttarget["true_img"], filter_hash, norm=false; dhc_args...)[coeff_mask])
s_init = DHC_compute_wrapper(gttarget["init"], filter_hash, norm=false; dhc_args...)[coeff_mask]
s_recon = DHC_compute_wrapper(gttarget["recon"], filter_hash, norm=false; dhc_args...)[coeff_mask]

struesq = zeros(34, 34)
struesq[reshape(coeff_mask[37:end], (34, 34))] .= s_true
sinitsq = zeros(34, 34)
sinitsq[reshape(coeff_mask[37:end], (34, 34))] .= s_init
sreconsq = zeros(34, 34)
sreconsq[reshape(coeff_mask[37:end], (34, 34))] .= s_recon
stargsq = zeros(34, 34)
stargsq[reshape(coeff_mask[37:end], (34, 34))] .= exp.(gttarget["dict"]["fs_targ_mean"])

s_true = log.(s_true)
s_init  = log.(s_init)
s_recon = log.(s_recon)
ltruecheck = ( 0.5 .* (s_true - s_targ_mean)' * s_targ_invcov * (s_true - s_targ_mean))[1]
l2true = 0.5*gttarget["dict"]["lambda"]*sum((alinit - altrue).^2)
println("True", ltruecheck)
println(l2true)

l1init = ( 0.5 .* (s_init - s_targ_mean)' * s_targ_invcov * (s_init - s_targ_mean))[1]
l2init = 0.5*gttarget["dict"]["lambda"]*sum((alinit - alinit).^2)
println("Init", l1init)
println(l2init)
println("From true", ( 0.5 .* (s_init - s_true)' * s_targ_invcov * (s_init - s_true))[1])

l1recon = ( 0.5 .* (s_recon - s_targ_mean)' * s_targ_invcov * (s_recon - s_targ_mean))[1]
l2recon = 0.5*gttarget["dict"]["lambda"]*sum((alrecon - alinit).^2)
println("Recon", l1recon)
println(l2recon)
println("From true", ( 0.5 .* (s_recon - s_true)' * s_targ_invcov * (s_recon - s_true))[1])

##Z-Score for mv gaussians test case
covmat = reshape([30, -15, -1, -15, 21, 3, -1, 3, 2], (3, 3))
mvg = Distributions.MultivariateNormal([-100, 3, 0.05], covmat)
samps = rand(mvg, 100)
samean = mean(samps, dims=2)
sacov = (samps .- samean) * (samps .- samean)'
sacov = sacov ./(99)
zscores = (samps .- samean) ./ (diag(sacov).^0.5)
#these do look like z scores

##Are the larger J value coefficients close to the true??
Nx=64
dbnimg = readsfd(Nx, logbool=false)
fname = "scratch_NM/StandardizedExp/Nx64/noisy_stdtrue/SFDTargSFDCov/reg_apd_noiso/LogCoeff/1000_0-01_full_triu"
gttarget = load(fname*".jld2")
dhc_args = gttarget["dhc_args"]
coeff_mask = gttarget["coeff_mask"]
fhash = gttarget["fhash"]
dbn = get_dbn_coeffs(dbnimg, fhash, dhc_args, coeff_mask = nothing)

s_true = DHC_compute_wrapper(gttarget["true_img"], fhash, norm=false; dhc_args...)
s_init = DHC_compute_wrapper(gttarget["init"], fhash, norm=false; dhc_args...)
s_recon = DHC_compute_wrapper(gttarget["recon"], fhash, norm=false; dhc_args...)
ffb, ffinfo = fink_filter_bank(1, 8, nx=64, Omega=true)

function J_hashindices(J_values)
    jindlist = []
    for jval in J_values
        push!(jindlist, findall(fhash["J_L"][:, 1].==jval))
    end
    return vcat(jindlist'...)
end

function J_S1indices(J_values)
    #Assumes this is applied to an object of length 2+Nf+Nf^2 or 2+Nf
    return J_hashindices(J_values) .+ 2
end

function distribution_percentiles(dbn, idx, perc)
    idxperc = zeros(size(idx))
    idxperc .= (x-> percentile(dbn[:, x], perc)).(idx)
    return idxperc
end

Jfiltidx = J_hashindices([0, 1, 2, 3])
JS1ind = J_S1indices([0, 1, 2, 3])
p25 = distribution_percentiles(dbn, JS1ind, 25)
p75 = distribution_percentiles(dbn, JS1ind, 75)
strueselect = (x->s_true[x]).(JS1ind)
sinitselect = (x->s_init[x]).(JS1ind)
sreconselect = (x->s_recon[x]).(JS1ind)
truebool = (strueselect .> p25) .& (strueselect .<p75)
initbool = (sinitselect .> p25) .& (sinitselect .<p75)
reconbool = (sreconselect .> p25) .& (sreconselect .<p75)

frac_devn_init_true = (sinitselect .- strueselect)./strueselect
frac_devn_recon_true = (sreconselect .- strueselect)./strueselect
dbn_mean = mean(dbn, dims=1)[:]
meanselect = (x->dbn_mean[x]).(JS1ind)
frac_devn_sfd_true = (meanselect .- strueselect)./strueselect
frac_devn_75_25 = (p75 .- p25)./p25

dbn_median = distribution_percentiles(dbn, JS1ind, 50)
frac_devn_medsfd_true = (dbn_median .- strueselect)./strueselect
clim = (minimum(frac_devn_medsfd_true), maximum(frac_devn_medsfd_true))
heatmap(frac_devn_init_true, title="Frac Devn Init-True", xlabel="L", ylabel="J", xticks=collect(1:8))
heatmap(frac_devn_sfd_true, title="Frac Devn SFDMean-True", xlabel="L", ylabel="J", xticks=collect(1:8))
heatmap(frac_devn_75_25, title="Frac Devn 75-25", xlabel="L", ylabel="J", xticks=collect(1:8))
heatmap(frac_devn_medsfd_true, title="Frac Devn SFDMedian-True", xlabel="L", ylabel="J", xticks=collect(1:8))
frac_devn_mean_med = (meanselect .- dbn_median)./dbn_median
heatmap(frac_devn_mean_med, title="Frac Devn Mean-Median", xlabel="L", ylabel="J", xticks=collect(1:8))

heatmap(frac_devn_medsfd_true, title="Frac Devn SFDMedian-True", xlabel="L", ylabel="J", xticks=collect(1:8), clim=clim)
heatmap(frac_devn_init_true, title="Frac Devn Init-True", xlabel="L", ylabel="J", xticks=collect(1:8), clim=clim)

##Does mean-subtraction make a difference?
meanim = mean(dbnimg, dims=[1, 2])
dbnimg_ms = dbnimg .- meanim
fname = "scratch_NM/StandardizedExp/Nx64/noisy_stdtrue/SFDTargSFDCov/reg_apd_noiso/LogCoeff/1000_0-01_full_triu"
gttarget = load(fname*".jld2")
dhc_args = gttarget["dhc_args"]
coeff_mask = gttarget["coeff_mask"]
fhash = gttarget["fhash"]
dbn = get_dbn_coeffs(dbnimg_ms, fhash, dhc_args, coeff_mask = nothing)

s_true = DHC_compute_wrapper(gttarget["true_img"] .- mean(gttarget["true_img"]), fhash, norm=false; dhc_args...)
s_init = DHC_compute_wrapper(gttarget["init"] .- mean(gttarget["init"]), fhash, norm=false; dhc_args...)
s_recon = DHC_compute_wrapper(gttarget["recon"] .- mean(gttarget["recon"]), fhash, norm=false; dhc_args...)
ffb, ffinfo = fink_filter_bank(1, 8, nx=64, Omega=true)

Jfiltidx = J_hashindices([0, 1, 2, 3])
JS1ind = J_S1indices([0, 1, 2, 3])
p25 = distribution_percentiles(dbn, JS1ind, 25)
p75 = distribution_percentiles(dbn, JS1ind, 75)
strueselect = (x->s_true[x]).(JS1ind)
sinitselect = (x->s_init[x]).(JS1ind)
sreconselect = (x->s_recon[x]).(JS1ind)
truebool = (strueselect .> p25) .& (strueselect .<p75)
initbool = (sinitselect .> p25) .& (sinitselect .<p75)
reconbool = (sreconselect .> p25) .& (sreconselect .<p75)

frac_devn_init_true = (sinitselect .- strueselect)./strueselect
frac_devn_recon_true = (sreconselect .- strueselect)./strueselect
dbn_mean = mean(dbn, dims=1)[:]
meanselect = (x->dbn_mean[x]).(JS1ind)
frac_devn_sfd_true = (meanselect .- strueselect)./strueselect
frac_devn_75_25 = (p75 .- p25)./p25

dbn_median = distribution_percentiles(dbn, JS1ind, 50)
frac_devn_medsfd_true = (dbn_median .- strueselect)./strueselect
heatmap(frac_devn_init_true, title="Frac Devn Init-True MeanSub", xlabel="L", ylabel="J", xticks=collect(1:8))
heatmap(frac_devn_sfd_true, title="Frac Devn SFDMean-True MeanSub", xlabel="L", ylabel="J", xticks=collect(1:8))
heatmap(frac_devn_75_25, title="Frac Devn 75-25 MeanSub", xlabel="L", ylabel="J", xticks=collect(1:8))
heatmap(frac_devn_medsfd_true, title="Frac Devn SFDMedian-True MeanSub", xlabel="L", ylabel="J", xticks=collect(1:8))
frac_devn_mean_med = (meanselect .- dbn_median)./dbn_median
heatmap(frac_devn_mean_med, title="Frac Devn Mean-Median MeanSub", xlabel="L", ylabel="J", xticks=collect(1:8))

##Making sure the log of coefficients didn't have to do with the results
Nx=64
dbnimg = readsfd(Nx, logbool=false)
fname = "scratch_NM/StandardizedExp/Nx64/noisy_stdtrue/SFDTargSFDCov/reg_apd_noiso/LambdaVary/1000_tunedlam_0-01_fullcov_triu"
gttarget = load(fname*".jld2")
dhc_args = gttarget["dhc_args"]
coeff_mask = gttarget["coeff_mask"]
fhash = gttarget["fhash"]
dbn = get_dbn_coeffs(dbnimg, fhash, dhc_args, coeff_mask = nothing)

s_true = DHC_compute_wrapper(gttarget["true_img"], fhash, norm=false; dhc_args...)
s_init = DHC_compute_wrapper(gttarget["init"], fhash, norm=false; dhc_args...)
s_recon = DHC_compute_wrapper(gttarget["recon"], fhash, norm=false; dhc_args...)
ffb, ffinfo = fink_filter_bank(1, 8, nx=64, Omega=true)

Jfiltidx = J_hashindices([0, 1, 2, 3])
JS1ind = J_S1indices([0, 1, 2, 3])
p25 = distribution_percentiles(dbn, JS1ind, 25)
p75 = distribution_percentiles(dbn, JS1ind, 75)
strueselect = (x->s_true[x]).(JS1ind)
sinitselect = (x->s_init[x]).(JS1ind)
sreconselect = (x->s_recon[x]).(JS1ind)
truebool = (strueselect .> p25) .& (strueselect .<p75)
initbool = (sinitselect .> p25) .& (sinitselect .<p75)
reconbool = (sreconselect .> p25) .& (sreconselect .<p75)

frac_devn_init_true = (sinitselect .- strueselect)./strueselect
frac_devn_recon_true = (sreconselect .- strueselect)./strueselect
dbn_mean = mean(dbn, dims=1)[:]
meanselect = (x->dbn_mean[x]).(JS1ind)
frac_devn_sfd_true = (meanselect .- strueselect)./strueselect
frac_devn_75_25 = (p75 .- p25)./p25

dbn_median = distribution_percentiles(dbn, JS1ind, 50)
frac_devn_medsfd_true = (dbn_median .- strueselect)./strueselect
p= heatmap(frac_devn_init_true, title="Frac Devn Init-True", xlabel="L", ylabel="J", xticks=collect(1:8))
savefig(p, "scratch_NM/NewWrapper/4-10/fdevn-init-true.png")
p=heatmap(frac_devn_sfd_true, title="Frac Devn SFDMean-True", xlabel="L", ylabel="J", xticks=collect(1:8))
savefig(p, "scratch_NM/NewWrapper/4-10/fdevn-sfdmean-true.png")
p= heatmap(frac_devn_75_25, title="Frac Devn 75-25", xlabel="L", ylabel="J", xticks=collect(1:8))
p= heatmap(frac_devn_medsfd_true, title="Frac Devn SFDMedian-True", xlabel="L", ylabel="J", xticks=collect(1:8))
savefig(p, "scratch_NM/NewWrapper/4-10/fdevn-sfdmedian-true.png")
frac_devn_mean_med = (meanselect .- dbn_median)./dbn_median
heatmap(frac_devn_mean_med, title="Frac Devn Mean-Median", xlabel="L", ylabel="J", xticks=collect(1:8))

p99 = distribution_percentiles(dbn, JS1ind, 90)
meanselect .> p99

##WHEN you used the log
fname = "scratch_NM/StandardizedExp/Nx64/noisy_stdtrue/SFDTargSFDCov/reg_apd_noiso/LambdaVary/1000_tunedlam_0-01_fullcov_triu"
gttarget = load(fname*".jld2")
dhc_args = gttarget["dhc_args"]
coeff_mask = gttarget["coeff_mask"]
fhash = gttarget["fhash"]
dbn = log.(get_dbn_coeffs(dbnimg, fhash, dhc_args, coeff_mask = nothing))

s_true = log.(DHC_compute_wrapper(gttarget["true_img"], fhash, norm=false; dhc_args...))
s_init = log.(DHC_compute_wrapper(gttarget["init"], fhash, norm=false; dhc_args...))
s_recon = log.(DHC_compute_wrapper(gttarget["recon"], fhash, norm=false; dhc_args...))
ffb, ffinfo = fink_filter_bank(1, 8, nx=64, Omega=true)

Jfiltidx = J_hashindices([0, 1, 2, 3])
JS1ind = J_S1indices([0, 1, 2, 3])
p25 = distribution_percentiles(dbn, JS1ind, 25)
p75 = distribution_percentiles(dbn, JS1ind, 75)
strueselect = (x->s_true[x]).(JS1ind)
sinitselect = (x->s_init[x]).(JS1ind)
sreconselect = (x->s_recon[x]).(JS1ind)
truebool = (strueselect .> p25) .& (strueselect .<p75)
initbool = (sinitselect .> p25) .& (sinitselect .<p75)
reconbool = (sreconselect .> p25) .& (sreconselect .<p75)

frac_devn_init_true = (sinitselect .- strueselect)./abs.(strueselect)
frac_devn_recon_true = (sreconselect .- strueselect)./abs.(strueselect)
dbn_mean = mean(dbn, dims=1)[:]
meanselect = (x->dbn_mean[x]).(JS1ind)
frac_devn_sfd_true = (meanselect .- strueselect)./abs.(strueselect)
frac_devn_75_25 = (p75 .- p25)./abs.(p25)

dbn_median = distribution_percentiles(dbn, JS1ind, 50)
frac_devn_medsfd_true = (dbn_median .- strueselect)./abs.(strueselect)
p= heatmap(frac_devn_init_true, title="Frac Devn Init-True", xlabel="L", ylabel="J", xticks=collect(1:8))
savefig(p, "scratch_NM/NewWrapper/4-10/LogCoeff/fdevn-init-true.png")
p=heatmap(frac_devn_sfd_true, title="Frac Devn SFDMean-True", xlabel="L", ylabel="J", xticks=collect(1:8))
savefig(p, "scratch_NM/NewWrapper/4-10/LogCoeff/fdevn-sfdmean-true.png")
p= heatmap(frac_devn_75_25, title="Frac Devn 75-25", xlabel="L", ylabel="J", xticks=collect(1:8))
p= heatmap(frac_devn_medsfd_true, title="Frac Devn SFDMedian-True", xlabel="L", ylabel="J", xticks=collect(1:8))
savefig(p, "scratch_NM/NewWrapper/4-10/LogCoeff/fdevn-sfdmedian-true.png")
frac_devn_mean_med = (meanselect .- dbn_median)./abs.(dbn_median)
heatmap(frac_devn_mean_med, title="Frac Devn Mean-Median", xlabel="L", ylabel="J", xticks=collect(1:8))

p99 = distribution_percentiles(dbn, JS1ind, 90)
meanselect .> p99

##Z-scores
Nx=64
dbnimg = readsfd(Nx, logbool=false)
fname = "scratch_NM/StandardizedExp/Nx64/noisy_stdtrue/SFDTargSFDCov/reg_apd_noiso/LambdaVary/1000_tunedlam_0-01_fullcov_triu"
gttarget = load(fname*".jld2")
dhc_args = gttarget["dhc_args"]
coeff_mask = gttarget["coeff_mask"]
fhash = gttarget["fhash"]
dbn = get_dbn_coeffs(dbnimg, fhash, dhc_args, coeff_mask = nothing)

sigdbn = std(dbn, dims=1)

s_true = DHC_compute_wrapper(gttarget["true_img"], fhash, norm=false; dhc_args...)
s_init = DHC_compute_wrapper(gttarget["init"], fhash, norm=false; dhc_args...)
s_recon = DHC_compute_wrapper(gttarget["recon"], fhash, norm=false; dhc_args...)
ffb, ffinfo = fink_filter_bank(1, 8, nx=64, Omega=true)

function J_hashindices(J_values)
    jindlist = []
    for jval in J_values
        push!(jindlist, findall(fhash["J_L"][:, 1].==jval))
    end
    return vcat(jindlist'...)
end

function J_S1indices(J_values)
    #Assumes this is applied to an object of length 2+Nf+Nf^2 or 2+Nf
    return J_hashindices(J_values) .+ 2
end

function zscores(a, b, std, idx)
    return (a[idx] .- b[idx])./std[idx]
end

Jfiltidx = J_hashindices([0, 1, 2, 3])
JS1ind = J_S1indices([0, 1, 2, 3])
strueselect = (x->s_true[x]).(JS1ind)
sinitselect = (x->s_init[x]).(JS1ind)
sreconselect = (x->s_recon[x]).(JS1ind)
stdselect = (x->sigdbn[x]).(JS1ind)

z_init_true = (sinitselect .- strueselect)./stdselect
z_recon_true = (sreconselect .- strueselect)./stdselect
dbn_mean = mean(dbn, dims=1)[:]
meanselect = (x->dbn_mean[x]).(JS1ind)
z_sfdmean_true = (meanselect .- strueselect)./stdselect

function distribution_percentiles(dbn, idx, perc)
    idxperc = zeros(size(idx))
    idxperc .= (x-> percentile(dbn[:, x], perc)).(idx)
    return idxperc
end

dbn_median = distribution_percentiles(dbn, JS1ind, 50)
z_sfdmed_true = (dbn_median .- strueselect)./stdselect

heatmap(z_init_true, title="Z Init-True", xlabel="L", ylabel="J", xticks=collect(1:8))
heatmap(z_sfdmean_true, title="Z SFDMean-True", xlabel="L", ylabel="J", xticks=collect(1:8))
heatmap(z_sfdmed_true, title="Z SFDMedian-True", xlabel="L", ylabel="J", xticks=collect(1:8))
z_mean_med = (meanselect .- dbn_median)./stdselect
heatmap(z_mean_med, title="Z Mean-Median", xlabel="L", ylabel="J", xticks=collect(1:8))

clim = (minimum(z_sfdmed_true), maximum(z_sfdmed_true))
heatmap(z_sfdmed_true, title="Z SFDMedian-True", xlabel="L", ylabel="J", xticks=collect(1:8), clim=clim)
heatmap(z_init_true, title="Z Init-True", xlabel="L", ylabel="J", xticks=collect(1:8), clim=clim)

badz_init_true = z_init_true./stdselect

##Log Coeff-Zscores

logdbn = log.(dbn)
sigdbn = std(logdbn, dims=1)

s_true = log.(DHC_compute_wrapper(gttarget["true_img"], fhash, norm=false; dhc_args...))
s_init = log.(DHC_compute_wrapper(gttarget["init"], fhash, norm=false; dhc_args...))
s_recon = log.(DHC_compute_wrapper(gttarget["recon"], fhash, norm=false; dhc_args...))
ffb, ffinfo = fink_filter_bank(1, 8, nx=64, Omega=true)

Jfiltidx = J_hashindices([0, 1, 2, 3])
JS1ind = J_S1indices([0, 1, 2, 3])
strueselect = (x->s_true[x]).(JS1ind)
sinitselect = (x->s_init[x]).(JS1ind)
sreconselect = (x->s_recon[x]).(JS1ind)
stdselect = (x->sigdbn[x]).(JS1ind)

z_init_true = (sinitselect .- strueselect)./stdselect
z_recon_true = (sreconselect .- strueselect)./stdselect
dbn_mean = mean(logdbn, dims=1)[:]
meanselect = (x->dbn_mean[x]).(JS1ind)
z_sfdmean_true = (meanselect .- strueselect)./stdselect

function distribution_percentiles(dbn, idx, perc)
    idxperc = zeros(size(idx))
    idxperc .= (x-> percentile(dbn[:, x], perc)).(idx)
    return idxperc
end

dbn_median = distribution_percentiles(logdbn, JS1ind, 50)
z_sfdmed_true = (dbn_median .- strueselect)./stdselect

heatmap(z_init_true, title="Z Init-True", xlabel="L", ylabel="J", xticks=collect(1:8))
heatmap(z_sfdmean_true, title="Z SFDMean-True", xlabel="L", ylabel="J", xticks=collect(1:8))
heatmap(z_sfdmed_true, title="Z SFDMedian-True", xlabel="L", ylabel="J", xticks=collect(1:8))
z_mean_med = (meanselect .- dbn_median)./stdselect
heatmap(z_mean_med, title="Z Mean-Median", xlabel="L", ylabel="J", xticks=collect(1:8))

clim = (minimum(z_sfdmed_true), maximum(z_sfdmed_true))
heatmap(z_sfdmed_true, title="Z SFDMedian-True", xlabel="L", ylabel="J", xticks=collect(1:8), clim=clim)
heatmap(z_init_true, title="Z Init-True", xlabel="L", ylabel="J", xticks=collect(1:8), clim=clim)

badz_init_true = z_init_true./stdselect

fnamelog = "scratch_NM/StandardizedExp/Nx64/noisy_stdtrue/SFDTargSFDCov/reg_apd_noiso/LogCoeff/1000_0-01_full_triu"
gtlog = load(fnamelog*".jld2")
s_invcov = gtlog["dict"]["fs_invcov"]
diag(s_invcov)
sigdbn[coeff_mask].^(-2)

##SOME DISCREPANCY BETWEEN THE FSTARG_INVCOV AND THE ONE YOU CALCULATE HERE FOR BOTH S_TARG_INVCOV AND LOGS_TARG_INVCOV:
#Resolved in check_invcov. Numerical instability issues.

Nx=64
dbnimg = readsfd(Nx, logbool=false)
fname = "scratch_NM/StandardizedExp/Nx64/noisy_stdtrue/SFDTargSFDCov/reg_apd_noiso/LambdaVary/1000_tunedlam_0-01_fullcov_triu"
gttarget = load(fname*".jld2")
dhc_args = gttarget["dhc_args"]
coeff_mask = gtlog["coeff_mask"]
fhash = gttarget["fhash"]
dbn = get_dbn_coeffs(dbnimg, fhash, dhc_args, coeff_mask = nothing)

sigdbn = std(dbn, dims=1)
saved = gttarget["dict"]["s_invcov"]
diag(saved)
sigdbn[coeff_mask].^(-2)
(dbn[:, coeff_mask]' * dbn[:, coeff_mask]) ./9999 #roughly consistent with sigdbn.^2

#INCONSISTENT WITH DIAG(SAVED)
datfile = "scratch_NM/StandardizedExp/Nx64/Data_10000.jld2" #Replace w SLURM array
loaddf = load(datfile)
true_img = loaddf["true_img"]
init = loaddf["init"]

filter_hash = fink_filter_hash(1, 8, nx=Nx, pc=1, wd=1, Omega=true)
(S1iso, Nf) = size(filter_hash["S1_iso_mat"])
(S2iso, Nfsq) = size(filter_hash["S2_iso_mat"])

dhc_args = Dict(:doS2=>false, :doS12=>false, :doS20=>true, :apodize=>true, :iso=>false) #Iso #CHANGE: Change sig for the sfd data since the noise model is super high and the tiny values make sense
if dhc_args[:iso]
    coeff_mask = falses(2+S1iso+S2iso)
    coeff_mask[S1iso+3:end] .= true
else #Not iso
    coeff_mask = falses(2+Nf+Nf^2)
    coeff_mask[Nf+3:end] .= true
end

white_noise_args = Dict(:loc=>0.0, :sig=>std(true_img), :Nsam=>1000, :norm=>false, :smooth=>false, :smoothval=>0.8) #Iso #only if you're using noise based covar
optim_settings = Dict([("iterations", 1000), ("norm", false), ("minmethod", ConjugateGradient())])
recon_settings = Dict([("target_type", "sfd_dbn"), ("covar_type", "sfd_dbn"), ("log", false), ("GaussianLoss", true), ("Invcov_matrix", "Diagonal+Eps"),
  ("optim_settings", optim_settings), ("white_noise_args", white_noise_args)]) #Add constraints

if recon_settings["log"] & dhc_args[:apodize]
    lval = std(apodizer(log.(true_img))).^(-2)
elseif !recon_settings["log"] & dhc_args[:apodize]
    lval = std(apodizer(true_img)).^(-2)
elseif recon_settings["log"] & !dhc_args[:apodize]
    lval = std(log.(true_img)).^(-2)
else
    lval = std(true_img).^(-2)
end

recon_settings["lambda"] = lval
coeff_mask = gtlog["coeff_mask"]
regs_true = DHC_compute_wrapper(log.(true_img), filter_hash; norm=false, dhc_args...)[coeff_mask]
s2mean, s2icov = meancov_generator(true_img, filter_hash, dhc_args, coeff_mask, recon_settings, safety=nothing)
diag(s2icov)
##THIS Is consistent but the saved s_targ isn't???? HOW

max_stars_repo_name: andrew-saydjari/DHC

max_issues_repo_name: andrew-saydjari/DHC

max_forks_repo_name: andrew-saydjari/DHC


text:
function InitialBurrau(T=Float64)

    Gm = [parse(BigFloat,"5"), parse(BigFloat,"4"),parse(BigFloat,"3")]

    N=length(Gm)
    x1 = parse(BigFloat,"1")
    y1 = parse(BigFloat,"-1")
    z1 = parse(BigFloat,"0")
    x2 = parse(BigFloat,"-2")
    y2 = parse(BigFloat,"-1")
    z2 = parse(BigFloat,"0")
    x3 = parse(BigFloat,"1")
    y3 = parse(BigFloat,"3")
    z3 = parse(BigFloat,"0")

    q =  [x1,y1,z1,x2,y2,z2,x3,y3,z3]
    v = zeros(size(q))
    q0 = reshape(q,3,:)
    v0 = reshape(v,3,:)
    u0 = Array{T}(undef,2,3,N)
    u0[1,:,:] = convert.(T,v0)
    u0[2,:,:] = convert.(T,q0)

    return u0, convert.(T,Gm)

end

ext: jl


text:
function test_rmse(truth::AbstractVector, prediction::AbstractVector)
    @test measure(RMSE(), truth, prediction) == 3.0
end

function test_mae(truth::AbstractVector, prediction::AbstractVector)
    @test measure(MAE(), truth, prediction) == 3.0
end

println("-- Testing accuracy metrics")

truth = [1, 2, 3]
prediction = [4, 5, 6]

test_rmse(truth, prediction)
test_mae(truth, prediction)


text:
using Test

@time begin
  include("interpolation_test.jl")
  include("RSwM1_test.jl")
  include("RSwM2_test.jl")
  include("RSwM3_test.jl")
  include("noise_wrapper.jl")
  include("noise_function.jl")
  include("VBT_test.jl")
  include("noise_grid.jl")
  include("noise_approximation.jl")
  include("sde_noise_wrapper.jl")
  include("multi_dim.jl")
  include("geometric_bm.jl")
  include("compoundpoisson.jl")
  include("ornstein.jl")
  include("bridge_test.jl")
  include("sde_adaptivedistribution_tests.jl")
  include("reversal_test.jl")
  include("two_processes.jl")
  include("extraction_test.jl")
  include("restart_test.jl")
end


text:
# ranges.jl

module RangesTest

using Test
using Dates
using QDates

@testset "ranges" begin

let T=QDates.QDate
    local f, l, dr, dr1
    f1 = T(1914); l1 = T(1913,12,30)
    f2 = T(2014); l2 = T(2014)
    f3 = T(1970); l3 = T(2020)
    # f4 = typemin(T); l4 = typemax(T)

    for P in (Dates.Day, Dates.Month, Dates.Year)
        for pos_step in (P(1),P(2),P(50))
            # empty range
            let dr = f1:pos_step:l1
                @test length(dr) == 0
                @test isempty(dr)
                @test first(dr) == f1
                @test last(dr) < f1
                @test length([i for i in dr]) == 0
                @test_throws ArgumentError minimum(dr)
                @test_throws ArgumentError maximum(dr)
                @test_throws BoundsError dr[1]
                @test findall(in(dr),dr) == Int64[]
                @test [dr;] == T[]
                @test isempty(reverse(dr))
                @test length(reverse(dr)) == 0
                @test first(reverse(dr)) < f1
                @test last(reverse(dr)) >= f1
                @test issorted(dr)
                @test sortperm(dr) == 1:1:0
                @test !(f1 in dr)
                @test !(l1 in dr)
                @test !(f1-pos_step in dr)
                @test !(l1+pos_step in dr)
            end

            # for (f,l) in ((f2,l2),(f3,l3),(f4,l4))
            for (f,l) in ((f2,l2),(f3,l3))
                let dr = f:pos_step:l, len = length(dr)
                    @test len > 0
                    @test typeof(len) <: Int64
                    @test !isempty(dr)
                    @test first(dr) == f
                    @test last(dr) <= l
                    @test minimum(dr) == first(dr)
                    @test maximum(dr) == last(dr)
                    @test dr[1] == f
                    @test dr[end] <= l
                    @static if VERSION < v"1.1.0-DEV.480"
                        @test iterate(dr) == (first(dr), 1)
                    else
                        @test iterate(dr) == (first(dr), (length(dr), 1))
                    end

                    if len < 10000
                        dr1 = [i for i in dr]
                        @test length(dr1) == len
                        @test findall(in(dr),dr) == [1:len;]
                        @test length([dr;]) == len
                    end
                    @test !isempty(reverse(dr))
                    @test length(reverse(dr)) == len
                    @test last(reverse(dr)) == f
                    @test issorted(dr)
                    @test f in dr
                end

            end
        end
        for neg_step in (P(-1),P(-2),P(-50))
            # empty range
            let dr = l1:neg_step:f1
                @test length(dr) == 0
                @test isempty(dr)
                @test first(dr) == l1
                @test last(dr) > l1
                @test length([i for i in dr]) == 0
                @test_throws ArgumentError minimum(dr)
                @test_throws ArgumentError maximum(dr)
                @test_throws BoundsError dr[1]
                @test findall(in(dr),dr) == Int64[]
                @test [dr;] == T[]
                @test isempty(reverse(dr))
                @test length(reverse(dr)) == 0
                @test first(reverse(dr)) > l1
                @test last(reverse(dr)) <= l1
                # @test !issorted(dr)
                # @test sortperm(dr) == 0:-1:1
                @test !(l1 in dr)
                @test !(l1 in dr)
                @test !(l1-neg_step in dr)
                @test !(l1+neg_step in dr)
            end

            # for (f,l) in ((f2,l2),(f3,l3),(f4,l4))
            for (f,l) in ((f2,l2),(f3,l3))
                let dr = l:neg_step:f, len = length(dr)
                    @test len > 0
                    @test typeof(len) <: Int64
                    @test !isempty(dr)
                    @test first(dr) == l
                    @test last(dr) >= f
                    @test minimum(dr) == last(dr)
                    @test maximum(dr) == first(dr)
                    @test dr[1] == l
                    @test dr[end] >= f
                    @static if VERSION < v"1.1.0-DEV.480"
                        @test iterate(dr) == (first(dr), 1)
                    else
                        @test iterate(dr) == (first(dr), (length(dr), 1))
                    end

                    if len < 10000
                        dr1 = [i for i in dr]
                        @test length(dr1) == len
                        @test findall(in(dr),dr) == [1:len;]
                        @test length([dr;]) == len
                    end
                    @test !isempty(reverse(dr))
                    @test length(reverse(dr)) == len
                    # @test !issorted(dr)
                    @test l in dr
                end
            end
        end
    end
end

# All the range representations we want to test
# Date ranges
dr  = QDates.QDate(2014,1,1):QDates.QDate(2014,2,1)
dr1 = QDates.QDate(2014,1,1):QDates.QDate(2014,1,1)
dr2 = QDates.QDate(2014,1,1):QDates.QDate(2013,2,1) # empty range
dr3 = QDates.QDate(2014,1,1):Dates.Day(-1):QDates.QDate(2013,1,1) # negative step
# Big ranges
dr4 = QDates.QDate(500):QDates.QDate(2000,1,1)
dr9 = typemin(QDates.QDate):typemax(QDates.QDate)
# Non-default steps
dr10 = typemax(QDates.QDate):Dates.Day(-1):typemin(QDates.QDate)
dr12 = typemin(QDates.QDate):Dates.Month(1):typemax(QDates.QDate)
dr13 = typemin(QDates.QDate):Dates.Year(1):typemax(QDates.QDate)
dr15 = typemin(QDates.QDate):Dates.Month(100):typemax(QDates.QDate)
dr16 = typemin(QDates.QDate):Dates.Year(1000):typemax(QDates.QDate)
dr20 = typemin(QDates.QDate):Dates.Day(2):typemax(QDates.QDate)

drs = Any[dr,dr1,dr2,dr3,dr4,dr9,dr10,
          dr12,dr13,dr15,dr16,dr20]

@test map(length,drs) == map(x->size(x)[1],drs)
@test all(x->findall(in(x),x) == [1:length(x);], drs[1:4])
@test isempty(dr2)
@test all(x->reverse(x) == last(x):-step(x):first(x),drs)
@test all(x->minimum(x) == (step(x) < zero(step(x)) ? last(x) : first(x)),drs[4:end])
@test all(x->maximum(x) == (step(x) < zero(step(x)) ? first(x) : last(x)),drs[4:end])
@test all(drs[1:3]) do dd
    for (i,d) in enumerate(dd)
        @test d == (first(dd) + Dates.Day(i-1))
    end
    true
end
# @test_throws MethodError dr + 1
a = QDates.QDate(2014,1,1)
b = QDates.QDate(2014,2,1)
@test map!(x->x+Dates.Day(1),Array{QDates.QDate}(undef, 30),dr) == [(a+Dates.Day(1)):(b+Dates.Day(1));]
@test map(x->x+Dates.Day(1),dr) == [(a+Dates.Day(1)):(b+Dates.Day(1));]

@test map(x->a in x,drs[1:4]) == [true,true,false,true]
@test a in dr
@test b in dr
@test QDates.QDate(2014,1,3) in dr
@test QDates.QDate(2014,1,15) in dr
@test QDates.QDate(2014,1,26) in dr
@test !(QDates.QDate(2013,1,1) in dr)

@test all(x->sort(x) == (step(x) < zero(step(x)) ? reverse(x) : x),drs)
@test all(x->step(x) < zero(step(x)) ? issorted(reverse(x)) : issorted(x),drs)

@test length(b:Dates.Day(-1):a) == 30
@test length(b:a) == 0
@test length(b:Dates.Day(1):a) == 0
@test length(a:Dates.Day(2):b) == 15
@test last(a:Dates.Day(2):b) == QDates.QDate(2014,1,29)
@test length(a:Dates.Day(7):b) == 5
@test last(a:Dates.Day(7):b) == QDates.QDate(2014,1,29)
@test length(a:Dates.Day(30):b) == 1
@test last(a:Dates.Day(30):b) == QDates.QDate(2014,1,1)
@test (a:b)[1] == QDates.QDate(2014,1,1)
@test (a:b)[2] == QDates.QDate(2014,1,2)
@test (a:b)[7] == QDates.QDate(2014,1,7)
@test (a:b)[end] == b
@test first(a:QDates.QDate(2099,1,1)) == a
@test first(a:typemax(QDates.QDate)) == a
@test first(typemin(QDates.QDate):typemax(QDates.QDate)) == typemin(QDates.QDate)

# Non-default step sizes
@test length(typemin(QDates.QDate):Dates.Month(1):typemax(QDates.QDate)) == QDates.QREF.LAST_RECORD  # == 21719
@test length(typemin(QDates.QDate):Dates.Year(1):typemax(QDates.QDate)) == QDates.QREF.LAST_YEAR - QDates.QREF.FIRST_YEAR + 1  # == 1756

c = QDates.QDate(2014,6,1)
@test length(a:Dates.Month(1):c) == 6
@test [a:Dates.Month(1):c;] == [a + Dates.Month(1)*i for i in 0:5]
@test [a:Dates.Month(2):QDates.QDate(2014,1,2);] == [a]
@test [c:Dates.Month(-1):a;] == reverse([a:Dates.Month(1):c;])
@test length(a:Dates.Month(1):QDates.QDate(2014,10,1)) == 11

d = QDates.QDate(2020,1,1)
@test length(a:Dates.Year(1):d) == 7
@test first(a:Dates.Year(1):d) == a
@test last(a:Dates.Year(1):d) == d
@test length(a:Dates.Month(12):d) == 7
@test first(a:Dates.Month(12):d) == a
@test last(a:Dates.Month(12):d) == QDates.QDate(2019,11,1)
@test length(a:Dates.Day(365):d) == 6
@test first(a:Dates.Day(365):d) == a
@test last(a:Dates.Day(365):d) == QDates.QDate(2018,12,25)

@test length(a:Dates.Year(1):QDates.QDate(2020,2,1)) == 7
@test length(a:Dates.Year(1):QDates.QDate(2020,6,1)) == 7
@test length(a:Dates.Year(1):QDates.QDate(2020,11,1)) == 7
@test length(a:Dates.Year(1):QDates.QDate(2020,12,30)) == 7
@test length(a:Dates.Year(1):QDates.QDate(2021,1,1)) == 8
@test length(QDates.QDate(2000):Dates.Year(-10):QDates.QDate(1900)) == 11
@test length(QDates.QDate(2000,6,23):Dates.Year(-10):QDates.QDate(1900,2,28)) == 11
@test length(QDates.QDate(2000,1,1):Dates.Year(1):QDates.QDate(2000,2,1)) == 1

# All leap months/years in 20th century
@test length(filter(QDates.isleapmonth, QDates.QDate(1901):Dates.Month(1):QDates.QDate(2000))) ==
      length(filter(QDates.isleapyear, QDates.QDate(1901):Dates.Year(1):QDates.QDate(2000))) == 36

end

end


text:
const always_skip = (
    :x, :y, :z,
    :camera,
    :scale, :rotation, :offset
)


text:
using Documenter, Omniscape

const formats = Any[
    Documenter.HTML(
    	assets = ["assets/custom.css"],
    ),
]

makedocs(
	format = formats,
    modules = [Omniscape],
    authors = "Vincent A. Landau",
    sitename = "Omniscape.jl",
    pages = ["Home" => "index.md",
    		 "Usage" => "usage.md",],
)

deploydocs(
    repo = "github.com/Circuitscape/Omniscape.jl.git"
)


text:
"""

# Example

```

```
"""
function GLPolyhedron(V::Lar.Points,
	EV::Lar.ChainOp,FE::Lar.ChainOp,CF::Lar.ChainOp)
	# TODO
end

"""

# Example

```

```
"""
function GLColoring(V::Lar.Points, FV::Lar.Cells, colors::GLColors) # single 3-cell

	# mesh building

	if size(points,2) == 2
		points = [points zeros(size(points,1),1)]
	end
	for line in lines
		  p2,p1 = points[line[1],:], points[line[2],:]
		  t=p2-p1;  n=LinearAlgebra.normalize([-t[2];+t[1];t[3]])

		  p1 = convert(GL.Point3d, p1)
		  p2 = convert(GL.Point3d, p2)
		  n  = convert(GL.Point3d,  n)

		  append!(vertices,p1); #append!(normals,n)
		  append!(vertices,p2); #append!(normals,n)
	end
	ret=GL.GLMesh(GL.GL_LINES)

        vertices,normals = GL.lar4mesh(points,triangles)

        ret=GL.GLMesh(GL.GL_TRIANGLES)
        ret.vertices = GL.GLVertexBuffer(vertices)
		ret.normals  = GL.GLVertexBuffer(normals)
		ret.colors  = GL.GLVertexBuffer(colors)
        return ret
end


text:
module ReferenceFree

export compute_reads, save_reads, load_reads

#%%

using XAM
using BioSequences
using ProgressMeter
using Serialization

include("utils/Utils.jl")
using .Utils

include("pmdtools/PMDTools.jl")
using .PMDTools

#%%

struct Reads
    sequences_original::Vector{LongDNASeq}
    sequences_aligned::Vector{LongDNASeq}
    PMDs::Vector{Float64}
end

# this is used to handle a call to `print`
function Base.show(io::IO, r::Reads)
    println("""Instance of Reads with fields:""")
    println("""   "sequences_original" """)
    println("""   "sequences_aligned" """)
    println("""   "PMDs" \n """)
    println(io, "sequences_original = ", r.sequences_original[1:5], "\n")
    println(io, "sequences_aligned = ", r.sequences_aligned[1:5], "\n")
    println(io, "PMDs = ", r.PMDs[1:5], "\n")
end

function Base.show(io::IO, ::MIME"text/plain", r::Reads)
    print(io, "Instance of Reads(sequences_original, sequences_aligned, PMDs)")
end

function save_reads(filename::String, reads::Reads)
    named_tuple = (
        sequences_original = reads.sequences_original,
        sequences_aligned = reads.sequences_aligned,
        PMDs = reads.PMDs,
    )
    serialize(filename, named_tuple)
end

function load_reads(filename::String)
    object = deserialize(filename)
    # for (k, v) in pairs(objects)
    #     @eval $k = $v
    # end
    return Reads(object.sequences_original, object.sequences_aligned, object.PMDs)
end

#%%

function compute_reads(filename; max_reads = -1)::Reads

    # Open a BAM file.
    reader = open(BAM.Reader, filename)

    sequences_original = LongDNASeq[]
    sequences_aligned = LongDNASeq[]
    PMDs = Float64[]

    counter = 0

    progress = ProgressUnknown("Reading $filename:", spinner = true)

    record = BAM.Record()
    while !eof(reader)
        empty!(record)
        read!(reader, record)
        ProgressMeter.next!(progress)

        # `record` is a BAM.Record object.
        if BAM.ismapped(record)
            seq_original, seq_aligned, PMD = compute(record)
            push!(sequences_original, seq_original)
            push!(sequences_aligned, seq_aligned)
            push!(PMDs, PMD)
            counter += 1
            if (max_reads > 0) & (counter >= max_reads)
                break
            end
        end

    end

    ProgressMeter.finish!(progress)
    close(reader)

    return Reads(sequences_original, sequences_aligned, PMDs)
end

#%%

function load_damaged_record(filename)
    "Load record with many errors"

    # Open a BAM file.
    reader = open(BAM.Reader, filename)

    record = BAM.Record()
    while !eof(reader)
        empty!(record)
        read!(reader, record)

        # `record` is a BAM.Record object.
        if BAM.ismapped(record)
            if length(BAM.cigar(record)) > 10 | length(record["MD"]::String) > 10
                break
            end
        end
    end

    # println(counter)
    close(reader)

    return record
end

#%%

end

max_stars_repo_name: ChristianMichelsen/reference-free

max_issues_repo_name: ChristianMichelsen/reference-free

max_forks_repo_name: ChristianMichelsen/reference-free


text:
using HDF5
using LCIO
using Glob
using StaticArrays
using LinearAlgebra

function getPrimaryParticle(collection, pdg::Int=22)
    for particle in collection
	if length(getParents(particle)) == 0 && abs(getPDG(particle)) == pdg
	    return particle
        end
    end
    return nothing
end

function getExpectedSeed(primary)
        p = SVector{3, Float64}(getMomentum(primary))
        mag = sqrt(sum(p.^2))
        p ./ mag .* 1.265 # ECAL inner radius
end

function computeCenterOfGravity(collection, expected)
    cellsForParticle = Vector{Vector{Float64}}()
    decode = CellIDDecoder(collection)
    for (cell_idx, hit) in enumerate(collection)
        layer = decode(hit)["layer"]
        e = getEnergy(hit)
        if layer > 20
            e *= 2
        end
        pos = SVector{3, Float64}(getPosition(hit))
        angle = dot(pos, expected) / sqrt(dot(pos, pos)) / sqrt(dot(expected, expected))
        if angle < 0
            continue
        end
        push!(cellsForParticle, [e, pos[1], pos[2], pos[3]])
    end
    eSum = 0.0
    cog = [0, 0, 0]
    for cell in cellsForParticle
        cog += cell[1] * [cell[2], cell[3], cell[4]]
        eSum += cell[1]
    end
    eSum, cog / eSum
end

# read the files for calibration
h5open("SiD_CaloCalibTraining_allElectrons_2calib.h5", "w") do file
for (file_idx, fileName) in enumerate(ARGS)
    LCIO.open(fileName) do reader
        cellsForParticle = zeros(5, length(reader))
        for (evt_idx, event) in enumerate(reader)
            # sum up the uncalibrated ECalHits
            # this needs to be sorted by layer, so we need a decoder
            eCalBarrelHits = getCollection(event, "ECalBarrelHits")
    	    primary = getPrimaryParticle(getCollection(event, "MCParticle"), 11)
            expected = getExpectedSeed(primary)
            trueE = getEnergy(primary)
            # nHits = length(EcalBarrelHits)
            # cellsForParticle = zeros(5, nHits)
            # for (cell_idx, hit) in enumerate(EcalBarrelHits)
            #     pos = SVector{3, Float64}(getPosition(hit))
            #     angle = dot(pos, expected) / sqrt(dot(pos, pos)) / sqrt(dot(expected, expected))
            #     if angle < 0
            #         continue
            #     end
            #     # calibrate the hits in the later layers with a higher number,
            #     # because they are behind thicker tungsten slabs
            #     e = getEnergy(hit)
            #     # rawHit = getRawHit(hit)
        	# 	#t = getTime(hit)
            #     #ts = 0#getTimeStamp(rawHit)
            #     #adc = 0#getAmplitude(rawHit)
            #     cellsForParticle[:, cell_idx] = [e, pos[1], pos[2], pos[3], layer]
            #     #push!(cellsForParticle, [e, pos[1], pos[2], pos[3], layer])
            # end
            (e, pos) = computeCenterOfGravity(eCalBarrelHits, expected)
            cellsForParticle[:, evt_idx] = [e, pos[1], pos[2], pos[3], trueE]
        end
        write(file, "file_$(file_idx)", cellsForParticle)
    end
end
end

max_stars_repo_name: jstrube/LCTutorials

max_issues_repo_name: jstrube/LCTutorials

max_forks_repo_name: jstrube/LCTutorials


text:
module ChanVese
using DistanceTransforms
using Images
using ImageFiltering

include("initializers.jl")
include("classical_chan_vese.jl")

export
    checkerboard,
    disk
    classical_chan_vese

end


text:
# TODO: writeup documentation

# Holder type for forward and backward plans, region, scalars etc..
# =================================================================
# Tf ≡ T_forward_arg
# Ti ≡ T_inverse_arg

struct FFTplan{Tf<:FFTN, d, Ti<:FFTN, Tsf<:Number, Tsi<:Number, FT<:Plan, IT<:Plan}
	unscaled_forward_transform::FT
	unscaled_inverse_transform::IT
	scale_forward::Tsf
	scale_inverse::Tsi
    sz_forward_arg::NTuple{d,Int}
    sz_inverse_arg::NTuple{d,Int}
	region::NTuple{d,Bool}

	function FFTplan{Tf,d}(uft::FT,uit::IT,sf::Tsf,si::Tsi,szf,szi,r) where {Tf<:FFTN,d,FT,IT,Tsf,Tsi}
		Ti = Complex{real(Tf)}
		return new{Tf,d,Ti,Tsf,Tsi,FT,IT}(uft,uit,sf,si,szf,szi,r)
	end
end

struct AdjointFFTplan{Tf<:FFTN, d, Ti<:FFTN, Tsf<:Number, Tsi<:Number, FT<:Plan, IT<:Plan}
	p::FFTplan{Tf, d, Ti, Tsf, Tsi, FT, IT}
end

# Constructors
# --------------------------------
# SizeInt and RegionBool are type wrappers to allow
# @generated plan

struct SizeInt{sz} end #e.g. (512,1024,2,4)#

struct RegionBool{rg} end #e.g. (false,true,false,true)#

@generated function XFields.plan(
		::Type{Tf},
		::Type{SizeInt{sz_forward_arg}},
		::Type{RegionBool{region}},
		scale_forward::Number
	) where {Tf<:FFTN, sz_forward_arg, region}

	d          = length(sz_forward_arg)
	region_tp  = tuple(findall(region)...)
	X          = Array{Tf,d}(undef, sz_forward_arg...)
	# TODO: figure out how to pass a global variable that sets this flag
	flags      = FFTW.ESTIMATE
	# flags      = FFTW.MEASURE
	# flags      = FFTW.PATIENT #fixme: bus error
	timelim    = -1.0 # -> no time limit

	if Tf <: FFTR

		T_backward_arg = Complex{Tf}
		sz_inverse_arg = tuple(FFTW.rfft_output_size(X, region_tp)...)
		Y = Array{T_backward_arg,d}(undef, sz_inverse_arg...)
		unscaled_forward_transform = FFTW.rFFTWPlan{Tf, -1,false,d}(X, Y, region_tp, flags, timelim)
		unscaled_inverse_transform = FFTW.rFFTWPlan{T_backward_arg, 1,false,d}(Y, X, region_tp, flags, timelim)

	elseif Tf <: FFTC

		sz_inverse_arg = sz_forward_arg
		Y = Array{Tf,d}(undef, sz_inverse_arg...)
		unscaled_forward_transform = FFTW.cFFTWPlan{Tf,-1,false,d}(X, Y, region_tp, flags, timelim)
		unscaled_inverse_transform = FFTW.cFFTWPlan{Tf, 1,false,d}(Y, X, region_tp, flags, timelim)

	end

	ifft_normalization = FFTW.normalization(real(Tf), sz_forward_arg, region_tp)

	return quote
 	       $(Expr(:meta, :inline))
 	       FFTplan{$Tf, $d}(
			$unscaled_forward_transform,
			$unscaled_inverse_transform,
			scale_forward,
			$ifft_normalization / scale_forward,
		    $sz_forward_arg,
		    $sz_inverse_arg,
			$region,
		)
 	end
end

function Base.adjoint(p::FFTplan)
	return AdjointFFT(p)
end

function Base.adjoint(p::AdjointFFTplan)
	return p.p
end

# Define how these plan holders operate
# -------------------------------

function LinearAlgebra.mul!(y::Array{Ti,d}, p::FFTplan{Tf,d,Ti}, x::Array{Tf,d}) where {d,Tf,Ti}
	mul!(y, p.unscaled_forward_transform, x)
	## rmul!(y, p.scale_forward)
	@inbounds y .*= p.scale_forward
	return y
end

function LinearAlgebra.ldiv!(x::Array{Tf,d}, p::FFTplan{Tf,d,Ti}, y::Array{Ti,d}) where {d,Tf,Ti}
	mul!(x, p.unscaled_inverse_transform, y)
	## rmul!(x, p.scale_inverse)
	@inbounds x .*= p.scale_inverse
	return x
end

function Base.:*(p::FFTplan{Tf,d}, x::Array{Tf,d}) where {d,Tf}
	## return LinearAlgebra.rmul!(p.unscaled_forward_transform * x, p.scale_forward)
	rtn = p.unscaled_forward_transform * x
	@inbounds rtn .*= p.scale_forward
	return rtn
end

function Base.:\(p::FFTplan{Tf,d,Ti}, y::Array{Ti,d}) where {d,Tf,Ti}
	## return LinearAlgebra.rmul!(p.unscaled_inverse_transform * y, p.scale_inverse)
	rtn   = p.unscaled_inverse_transform * y
	@inbounds rtn .*= p.scale_inverse
	return rtn
end

# the adjoint * x is the unscaled inverse transform but with forward scaling

function Base.:*(p::AdjointFFTplan{Tf,d,Ti}, x::Array{Ti,d}) where {d,Tf,Ti}
	## return LinearAlgebra.rmul!(p.p.unscaled_inverse_transform * x, p.p.scale_forward)
	rtn = p.p.unscaled_inverse_transform * x
	@inbounds rtn .*= p.p.scale_forward
end

function Base.:\(p::AdjointFFTplan{Tf,d,Ti}, y::Array{Tf,d}) where {d,Tf,Ti}
	## return LinearAlgebra.rmul!(p.p.unscaled_forward_transform * y, p.p.scale_inverse)
	rtn   = p.p.unscaled_forward_transform * y
	@inbounds rtn .*= p.p.scale_inverse
	return rtn
end

# Extracting/converting to a real plan and/or a complex plan
# -------------------------------

function Base.real(p::FFTplan{Tf,d,Ti,Sf,Si,FT,IT}) where {Tf,d,Ti,Sf,Si,FT,IT}
	Tf′ = real(Tf)
	G    = NTuple{sum(p.region),Int}
	FT′ = FFTW.rFFTWPlan{Tf′,-1, false, d, G}
	IT′ = FFTW.rFFTWPlan{Ti , 1, false, d, G}
	return plan(
		Tf′,
		SizeInt{p.sz_forward_arg},
		RegionBool{p.region},
		p.scale_forward
	)::FFTplan{Tf′,d,Ti,Sf,Si,FT′,IT′}
end

function Base.complex(p::FFTplan{Tf,d,Ti,Sf,Si,FT,IT}) where {Tf,d,Ti,Sf,Si,FT,IT}
	Tf′ = Complex{real(Tf)}
	G   = NTuple{sum(p.region),Int}
	FT′ = FFTW.cFFTWPlan{Tf′,-1, false, d, G}
	IT′ = FFTW.cFFTWPlan{Ti , 1, false, d, G}
	return plan(
		Tf′,
		SizeInt{p.sz_forward_arg},
		RegionBool{p.region},
		p.scale_forward
	)::FFTplan{Tf′,d,Ti,Sf,Si,FT′,IT′}
end


text:
include("ExpandingBranch.jl")
include("IntervalMap.jl")
include("CircleMap.jl")
include("ComposedMaps.jl")
include("conversion.jl")
include("conjugacy.jl")
include("examples.jl")


text:
using BandwidthBenchmark

bwbench(; verbose=true)
flush(stdout)

max_stars_repo_name: carstenbauer/BandwidthBenchmark

max_issues_repo_name: carstenbauer/BandwidthBenchmark

max_forks_repo_name: carstenbauer/BandwidthBenchmark


text:
"""
    logSequencingSolution(bestKnownSolution::Vector{Int}, bestKnownValue::T;firstLine::Union{String, Nothing}=nothing,fileName::Union{Nothing,String}=nothing) where T<:Real

Print a solution to the terminal

# Arguments
- `bestKnownSolution::Vector{Int}`: The solution as an ordered list of nodes
- `bestKnownValue::T`: The cost of the solution
- `firstLine::Union{String, Nothing}`:An optional preface to the solution
- `fileName::Union{Nothing,String}`: An optional file to write the log statement to
"""
function logSequencingSolution(bestKnownSolution::Vector{Int}, bestKnownValue::T;firstLine::Union{String, Nothing}=nothing,fileName::Union{Nothing,String}=nothing) where T<:Real
    if !isnothing(firstLine)
        println(firstLine)
    end
    print("Sequence: ",)
    for i in 1:length(bestKnownSolution)
        print(bestKnownSolution[i], " ")
    end
    if isinteger(bestKnownValue)
        println("\nCost: ", Int(bestKnownValue))
    else
        println("\nCost: ", bestKnownValue)
    end
    if !isnothing(fileName)
        f = open(fileName,"a")
            if !isnothing(firstLine)
                write(f, string(firstLine, "\n"))
            end
            write(f, "Sequence: ")
            for i in 1:length(bestKnownSolution)
                write(f, string(bestKnownSolution[i], " "))
            end
            if isinteger(bestKnownValue)
                write(f, string("\nCost: ", Int(bestKnownValue),"\n"))
            else
                write(f, string("\nCost: ", bestKnownValue,"\n"))
            end
        close(f)
    end
end

"""
    logOptimalityGap(lower::T, upper::U;fileName::Union{Nothing,String}=nothing)where{T<:Real,U<:Real}

Print an optimality gap to the terminal

# Arguments
- `lower::T`: The smaller number
- `upper::U`: The larger number
- `fileName::Union{Nothing,String}`: An optional file to write the log statement to
"""
function logOptimalityGap(lower::T, upper::U;fileName::Union{Nothing,String}=nothing)where{T<:Real,U<:Real}
    println("Bounds: [",lower,",",upper,"]")
    println("Optimality Gap: ",round((upper-lower)/upper,digits=4)*100,"%")
    if !isnothing(fileName)
        f = open(fileName,"a")
            write(f, string("Bounds: [",lower,",",upper,"]", "\n"))
            write(f, string("Optimality Gap: ",round((upper-lower)/upper,digits=4)*100,"%", "\n"))
            write(f, string("Data: ",lower,",",upper,",",round((upper-lower)/upper,digits=4)*100,"\n"))
        close(f)
    end
end

"""
    logDDWidths(dd::Vector{Vector{RelaxedSequencingNode}})

Print the width of each layer of a dd to the terminal

# Arguments
- `dd::Vector{Vector{RelaxedSequencingNode}}`: The dd
"""
function logDDWidths(dd::Vector{Vector{RelaxedSequencingNode}})
    print("DD Layer Widths: [")
    if !isempty(dd)
        print(length(dd[1]))
        for i in 2:length(dd)
            print(",",length(dd[i]))
        end
    end
    println("]")
end

"""
    logRuntime(startTime::U, additionalTimeElapsed::V;fileName::Union{Nothing,String}=nothing)where{U<:Real, V<:Real}

logs the current runtime

# Arguments
- `startTime::U`: the time the program started (in seconds since the epoch)
- `additionalTimeElapsed::V`: additional time used (in seconds)
- `fileName::Union{Nothing,String}`: An optional file to write the log statement to
"""
function logRuntime(startTime::U, additionalTimeElapsed::V;fileName::Union{Nothing,String}=nothing)where{U<:Real, V<:Real}
    println(round(time() - startTime + additionalTimeElapsed, digits=2), " seconds")
    if !isnothing(fileName)
        f = open(fileName,"a")
            write(f, string(round(time() - startTime + additionalTimeElapsed, digits=2), " seconds", "\n"))
        close(f)
    end
end

"""
    logQueueLength(qLength::U;fileName::Union{Nothing,String}=nothing)where{U<:Int}

Log the length of the queue

# Arguments
- `qLength::U`: the queue length
- `fileName::Union{Nothing,String}`: An optional file to write the log statement to
"""
function logQueueLength(qLength::U;fileName::Union{Nothing,String}=nothing)where{U<:Int}
    println("Nodes In Queue: ", qLength)
    if !isnothing(fileName)
        f = open(fileName,"a")
            write(f, string("Nodes In Queue: ",qLength, "\n"))
        close(f)
    end
end


text:
using Test, SpecialFunctions
import TruncatedNormal as TruncNorm

@test TruncNorm.tnmom2c(0, -Inf, Inf) == 1
@test TruncNorm.tnmom2c(0, -Inf, 0) == 1
@test TruncNorm.tnmom2c(0, 0, +Inf) == 1

@test TruncNorm.tnmom2c(0, +Inf, +Inf) == Inf
@test TruncNorm.tnmom2c(0, -Inf, -Inf) == Inf

@test isnan(TruncNorm.tnmom2c(0, 1, 0))
@test isnan(TruncNorm.tnmom2c(0, +Inf, 0))
@test isnan(TruncNorm.tnmom2c(0, 0, -Inf))

@test isnan(TruncNorm.tnmom2c(0, NaN, NaN))
@test isnan(TruncNorm.tnmom2c(0, 0, NaN))
@test isnan(TruncNorm.tnmom2c(0, NaN, 0))

@test TruncNorm.tnmom2c(10, 0, 10) ≈ 85.0423087839426928824034533318
@test TruncNorm.tnmom2c(1, -3, 3) ≈ 1.97333692466254147658812248690
@test TruncNorm.tnmom2c(1, 3, 4) ≈ 5.15893137101968604368598313652
@test TruncNorm.tnmom2c(-1, 3, 4) ≈ 18.2007485133797706592985387254
@test TruncNorm.tnmom2c(2, 3, 4) ≈ 1.63802279983964373587970534208
@test TruncNorm.tnmom2c(-2, 3, 4) ≈ 27.7216570845598129671048165198

for x = -10:10
    @test TruncNorm.tnmom2c(0, x, +Inf) ≈ 1 + x * √(2/π) / erfcx(+x/√2)
    @test TruncNorm.tnmom2c(0, -Inf, x) ≈ 1 - x * √(2/π) / erfcx(-x/√2)
end

@test TruncNorm.tnmom2c(0, -Inf, Inf, 0, 1) == 1

for x = -10:10
    @test TruncNorm.tnmom2c(0, x, +Inf, 0, 1) ≈ TruncNorm.tnmom2(x, +Inf)
    @test TruncNorm.tnmom2c(0, -Inf, x, 0, 1) ≈ TruncNorm.tnmom2(-Inf, x)
end

@test TruncNorm.tnmom2c(0, -1e200, 1e200) ≈ 1
@test TruncNorm.tnmom2c(0, -1e6, -999000) ≈ 9.9800100000199999999999799599399200002e11
@test TruncNorm.tnmom2c(0, +1e6, +Inf) ≈ 1.00000000000199999999999800000e12
@test TruncNorm.tnmom2c(0, -Inf, -1e6) ≈ 1.00000000000199999999999800000e12

@test TruncNorm.tnmom2c(4, -10, 20, 3, 3//2) ≈ 3.24999999999999956053727803062


text:
export moment, moment_value, monomial, base_function

struct Moment{T, PT <: AbstractPolynomialLike} <: AbstractMoment{T}
    α::T
    x::PT
end

"""
    moment(α, m::AbstractPolynomialLike)

Creates the moment of the base function `m` of value `α`.
"""
moment(α, m::AbstractPolynomialLike) = Moment(α, m)

"""
    moment_value(m::AbstractMomentLike)

Returns the value of the moment `m`.

## Examples

Calling `moment_value(moment(3.1, x*y^2))` should return `3.1`.
"""
moment_value(m::Moment) = m.α

"""
    monomial(m::AbstractMomentLike)

Returns the monomial of the moment `m`.

## Examples

Calling `monomial(moment(3.1, x*y^2))` should return `x*y^2`.
"""
MP.monomial(m::Moment) = m.x

"""
    base_function(m::AbstractMomentLike)

Returns the base_function of the moment `m`.

## Examples

Calling `base_function(moment(3.1, x*y^2))` should return `x*y^2`.
"""
base_function(m::Moment) = m.x

for f in [:variables, :nvariables, :exponents, :powers]
    @eval begin
        MP.$f(m::AbstractMomentLike) = $f(base_function(m))
    end
end

MP.degree(m::AbstractMomentLike) = MP.maxdegree(base_function(m))


text:
ignore = false  # when re-generating

include("text-align.jl")
# comparison failing but looks ok
#ignore || compare("text-align.png")

include("ndy.jl")
ignore || compare("ndy.png")

paint(PNG("text-rotate.png", 100, 100),
      Axes(negative=true),
      Font(size=0.5)) do
    text("abcdefg"; align=4, rotate=30)
end
ignore || compare("text-rotate.png")


text:
using BenchmarkTools
using BSON
using Benchmarks

# Benchmark parameters

N = 512
n_threads = min.(2 .^ (0:10), Sys.CPU_THREADS) |> unique

# Run and collect benchmarks

print_system_info()

for t in n_threads
    @info "Benchmarking multithreading (N=$N, threads=$t)..."
    julia = Base.julia_cmd()
    run(`$julia -t $t --project benchmark_multithreading_single.jl $N`)
end

suite = BenchmarkGroup(["size", "threads"])
for t in n_threads
    suite[(N, t)] = BSON.load("multithreading_benchmark_$t.bson")[:trial]
end

# Summarize benchmarks

df = benchmarks_dataframe(suite)
sort!(df, :threads)
benchmarks_pretty_table(df, title="Multithreading benchmarks")

suite_Δ = speedups_suite(suite, base_case=(N, 1))
df_Δ = speedups_dataframe(suite_Δ)
sort!(df_Δ, :threads)
benchmarks_pretty_table(df_Δ, title="Multithreading speedup")


text:
"$(SIGNATURES) Remove bounds constraints from constraint set"
function remove_bounds!(conSet::ConstraintSet)
    bnds = filter(is_bound, conSet.constraints)
    n,m = conSet.n, conSet.m
    filter!(x->!is_bound(x), conSet.constraints)
    num_constraints!(conSet)  # re-calculate number of constraints after removing bounds
	return bnds
end

@inline remove_goals!(conSet::ConstraintSet) = remove_constraint_type!(conSet, GoalConstraint)

"$(SIGNATURES) Remove a type of constraint from constraint set"
function remove_constraint_type!(conSet::ConstraintSet, ::Type{Con}) where Con <: AbstractConstraint
	goals = filter(x->x.con isa Con, conSet.constraints)
	filter!(x->!(x.con isa Con), conSet.constraints)
    num_constraints!(conSet)  # re-calculate number of constraints after removing goals
	return goals
end

"$(SIGNATURES) Remove bounds from constraint set and return them as vectors"
function get_bounds(conSet::ConstraintSet)
    N = length(conSet.p)

	bnds = remove_bounds!(conSet)
	n,m = size(conSet)
    inds = [bnd.inds for bnd in bnds]

    # Make sure the bounds don't overlap
    if !isempty(bnds)
        Ltotal = mapreduce(length,+,inds)  # total number of knot points covered by bounds
        Linter = length(reduce(union,inds))  # number of knot points with a bound
        @assert Ltotal == Linter
    end

    # Get primal bounds
    zL = [-ones(n+m)*Inf for k = 1:N]
    zU = [ ones(n+m)*Inf for k = 1:N]

    for k = 1:N
        for bnd in bnds
            if k ∈ bnd.inds
                zL[k] = bnd.con.z_min
                zU[k] = bnd.con.z_max
            end
        end
    end

	# Set bounds on goal constraints
	goals = remove_goals!(conSet)
	for goal in goals
		for (i,k) in enumerate(goal.inds)
			zL[k][goal.con.inds] = goal.con.xf
			zU[k][goal.con.inds] = goal.con.xf
		end
	end
	zU = vcat(zU...)
	zL = vcat(zL...)

    # Get bounds for constraints
    p = conSet.p
    NP = sum(p)
    gL = -Inf*ones(NP)
    gU =  Inf*ones(NP)
    cinds = gen_con_inds(conSet)
    for k = 1:N
        for (i,con) in enumerate(conSet.constraints)
            if k ∈ con.inds
                j = _index(con,k)
                gL[cinds[i][j]] = lower_bound(con)
                gU[cinds[i][j]] = upper_bound(con)
            end
        end
    end

    convertInf!(zU)
    convertInf!(zL)
    convertInf!(gU)
    convertInf!(gL)
    return zU, zL, gU, gL
end

"```julia
add_dynamics_constraints!(prob::Problem)
```
Add dynamics constraints to the constraint set"
function add_dynamics_constraints!(prob::Problem{Q}, integration=Q) where Q
	n,m = size(prob)
    conSet = prob.constraints

    # Implicit dynamics
    dyn_con = ConstraintVals(DynamicsConstraint{integration}(prob.model, prob.N), 1:prob.N-1)
    add_constraint!(conSet, dyn_con, 1)

    # Initial condition
    init_con = ConstraintVals( GoalConstraint(prob.x0), 1:1)
    add_constraint!(conSet, init_con, 1)

    return nothing
end

"$(SIGNATURES) Generate the indices into the concatenated constraint vector for each constraint.
Determines the bandedness of the Jacobian"
function gen_con_inds(conSet::ConstraintSet, structure=:by_knotpoint)
	n,m = size(conSet)
    N = length(conSet.p)
    numcon = length(conSet.constraints)
    conLen = length.(conSet.constraints)

    cons = [[@SVector ones(Int,length(con)) for i in eachindex(con.inds)] for con in conSet.constraints]

    # Dynamics and general constraints
    idx = 0
	if structure == :by_constraint
	    for (i,con) in enumerate(conSet.constraints)
			for (j,k) in enumerate(con.inds)
				cons[i][_index(con,k)] = idx .+ (1:conLen[i])
				idx += conLen[i]
	        end
	    end
	elseif structure == :by_knotpoint
		for k = 1:N
			for (i,con) in enumerate(conSet.constraints)
				if k in con.inds
					j = _index(con,k)
					cons[i][j] = idx .+ (1:conLen[i])
					idx += conLen[i]
				end
			end
		end
	end
    return cons
end

"$(SIGNATURES)
Get the constraint Jacobian structure as a sparse array, and fill in the linear indices
used for filling a vector of the non-zero elements of the Jacobian"
function constraint_jacobian_structure(solver::DirectSolver,
		structure=:by_knopoint)
    n,m,N = size(solver)
    conSet = get_constraints(solver)
    idx = 0.0
    linds = jacobian_linear_inds(solver)

    NN = num_primals(solver)
    NP = num_duals(solver)
    D = spzeros(Int,NP,NN)

    # Number of elements in each block
    blk_len = map(con->length(con.∇c[1]), conSet.constraints)

    # Number of knot points for each constraint
    con_len = map(con->length(con.∇c), conSet.constraints)

    # Linear indices
	if structure == :by_constraint
	    for (i,con) in enumerate(conSet.constraints)
	        for (j,k) in enumerate(con.inds)
	            inds = idx .+ (1:blk_len[i])
	            linds[i][j] = inds
	            con.∇c[j] = inds
	            idx += blk_len[i]
	        end
	    end
	elseif structure == :by_knotpoint
		for k = 1:N
			for (i,con) in enumerate(conSet.constraints)
				if k in con.inds
					inds = idx .+ (1:blk_len[i])
					j = _index(con,k)
					linds[i][j] = inds
					con.∇c[j] = inds
					idx += blk_len[i]
				end
			end
		end
	end

    copy_jacobians!(D, solver)
    return D
end


text:
@testset "eachobs" begin
    for (i,x) in enumerate(eachobs(X))
        @test x == X[:,i]
    end

    for (i,x) in enumerate(eachobs(X, buffer=true))
        @test x == X[:,i]
    end

    b = zeros(size(X, 1))
    for (i,x) in enumerate(eachobs(X, buffer=b))
        @test x == X[:,i]
    end
    @test b == X[:,end]

    @testset "batched" begin
        for (i, x) in enumerate(eachobs(X, batchsize=2, partial=true))
            if i != 8
                @test size(x) == (4,2)
                @test x == X[:,2i-1:2i]
            else
                @test size(x) == (4,1)
                @test x == X[:,2i-1:2i-1]
            end
        end

        for (i, x) in enumerate(eachobs(X, batchsize=2, buffer=true, partial=false))
            @test size(x) == (4,2)
            @test x == X[:,2i-1:2i]
        end

        b = zeros(4, 2)
        for (i, x) in enumerate(eachobs(X, batchsize=2, buffer=b, partial=false))
            @test size(x) == (4,2)
            @test x == X[:,2i-1:2i]
        end
    end
end


text:
# MCL (Markov CLustering algorithm)

"""
The output of [`mcl`](@ref) function.

# Fields
 - `mcl_adj::AbstractMatrix`: the final MCL adjacency matrix
   (equilibrium state matrix if the algorithm converged), empty if
   `save_final_matrix` option is disabled
 - `assignments::Vector{Int}`: indices of the points clusters.
   `assignments[i]` is the index of the cluster for the ``i``-th point
    (``0`` if unassigned)
 - `counts::Vector{Int}`: the ``k``-length vector of cluster sizes
 - `nunassigned::Int`: the number of standalone points not assigned to any
   cluster
 - `iterations::Int`: the number of elapsed iterations
 - `rel_Δ::Float64`: the final relative Δ
 - `converged::Bool`: whether the method converged
"""
struct MCLResult <: ClusteringResult
    mcl_adj::AbstractMatrix     # final MCL adjacency matrix (equilibrium state matrix if converged)
    assignments::Vector{Int}    # point-to-cluster assignments (n)
    counts::Vector{Int}         # number of points assigned to each cluster (k)
    nunassigned::Int            # number of single elements not assigned to any cluster
    iterations::Int             # number of elapsed iterations
    rel_Δ::Float64              # final relative Δ
    converged::Bool             # whether the procedure converged
end

# Extract clusters from the final (equilibrium) MCL matrix
# Return the tuple: cluster indices for each element, cluster sizes,
# the number of unassigned (0 cluster index) elements (if `allow_singles` is on)
# `zero_tol` is a minimal value to consider as an element-to-cluster assignment
function _mcl_clusters(mcl_adj::AbstractMatrix, allow_singles::Bool, zero_tol::Float64 = 1E-20)
    # remove rows containing only zero elements and convert into a mask of nonzero elements
    el2clu_mask = mcl_adj[dropdims(sum(mcl_adj, dims=2), dims=2) .> zero_tol, :] .> zero_tol

    # assign cluster indexes to each node
    # cluster index is the index of the first TRUE in a given column
    _ms = mapslices(el_mask->isempty(el_mask) ? 0 : argmax(el_mask), el2clu_mask, dims=1)
    clu_ixs = dropdims(_ms, dims=1)
    clu_sizes = zeros(Int, size(el2clu_mask, 1))
    unassigned_count = 0
    @inbounds for clu_ix in clu_ixs
        (clu_ix > 0) && (clu_sizes[clu_ix] += 1)
    end
    if !allow_singles
        # collapse all size 1 clusters into one with index 0
        @inbounds for i in eachindex(clu_ixs)
            clu_ix = clu_ixs[i]
            if clu_ix > 0 && clu_sizes[clu_ix] == 1
                clu_ixs[i] = 0
                clu_sizes[clu_ix] = 0
                unassigned_count += 1
            end
        end
    else
        unassigned_count = 0
    end

    # recode clusters numbers to be in 1:N range (or 0:N if there's collapsed cluster)
    clu_id_map = zeros(Int, length(clu_sizes))
    next_clu_ix = 0
    @inbounds for i in eachindex(clu_ixs)
        old_clu_ix = clu_ixs[i]
        if old_clu_ix > 0
            new_clu_ix = clu_id_map[old_clu_ix]
            clu_ixs[i] = new_clu_ix == 0 ?
                         clu_id_map[old_clu_ix] = (next_clu_ix += 1) :
                         new_clu_ix
        end
    end
    old_clu_sizes = clu_sizes
    clu_sizes = zeros(Int, next_clu_ix)
    for (old_clu_ix, new_clu_ix) in enumerate(clu_id_map)
        if new_clu_ix > 0
            clu_sizes[new_clu_ix] = old_clu_sizes[old_clu_ix]
        end
    end

    clu_ixs, clu_sizes, unassigned_count
end

# adjacency matrix expansion (matrix-wise raising to a given power) kernel
function _mcl_expand(src::AbstractMatrix, expansion::Number)
    try
        return src^expansion
    catch ex # FIXME: remove this when functionality become available in the standard library
        if isa(ex, MethodError)
            throw(ArgumentError("MCL expansion of $(typeof(src)) with expansion=$expansion not supported"))
        else
            rethrow()
        end
    end
end

# integral power inflation (single matrix element)
_mcl_el_inflate(el::Number, inflation::Integer) = el^inflation

# non-integral power inflation (single matrix element)
_mcl_el_inflate(el::Number, inflation::Number) = real((el+0im)^inflation)

# adjacency matrix inflation (element-wise raising to a given power) kernel
function _mcl_inflate!(dest::AbstractMatrix, src::AbstractMatrix, inflation::Number)
    @inbounds for i in eachindex(src)
        dest[i] = _mcl_el_inflate(src[i], inflation)
    end
end

# adjacency matrix pruning
function _mcl_prune!(mtx::AbstractMatrix, prune_tol::Number)
    for i in 1:size(mtx,2)
        c = view(mtx, :, i)
        θ = mean(c)*prune_tol
        @inbounds @simd for j in eachindex(c)
            c[j] = ifelse(c[j] >= θ, c[j], 0.0)
        end
    end
    issparse(mtx) && dropzeros!(mtx)
    return mtx
end

"""
    mcl(adj::AbstractMatrix; [kwargs...])

Perform MCL (Markov Cluster Algorithm) clustering using ``n×n``
adjacency (points similarity) matrix `adj`.

Returns [`MCLResult`](@ref) object.

# Algorithm Options
 - `add_loops::Bool` (enabled by default): whether the edges of weight 1.0
   from the node to itself should be appended to the graph
 - `expansion::Number` (defaults to 2): MCL *expansion* constant
 - `inflation::Number` (defaults to 2): MCL *inflation* constant
 - `save_final_matrix::Bool` (disabled by default): whether to save the final
   equilibrium state in the `mcl_adj` field of the result; could provide useful
   diagnostic if the method doesn't converge
 - `prune_tol::Number`: pruning threshold
 - `display::Symbol` (defaults to `:none`): `:none` for no output or `:verbose`
   for diagnostic messages
 - `max_iter`, `tol`: see [common options](@ref common_options)

### References
> Stijn van Dongen, *"Graph clustering by flow simulation"*, 2001

> [Original MCL implementation](http://micans.org/mcl).
"""
function mcl(adj::AbstractMatrix{T};
             add_loops::Bool = true,
             expansion::Number = 2, inflation::Number = 2,
             save_final_matrix::Bool = false,
             allow_singles::Bool = true,
             max_iter::Integer = 100, tol::Number=1.0e-5,
             prune_tol::Number=1.0e-5, display::Symbol=:none) where T<:Real
    m, n = size(adj)
    m == n || throw(DimensionMismatch("Square adjacency matrix expected"))

    if add_loops
        @inbounds for i in 1:size(adj, 1)
            adj[i, i] = 1.0
        end
    end

    # initialize the MCL adjacency matrix by normalized `adj` weights
    mcl_adj = copy(adj)
    # normalize in columns
    rmul!(mcl_adj, Diagonal(map(x -> x != 0.0 ?  1.0/x : x, dropdims(sum(mcl_adj, dims=1), dims=1))))
    mcl_norm = norm(mcl_adj)
    if !isfinite(mcl_norm)
        throw(OverflowError("The norm of the input adjacency matrix is not finite"))
    end
    next_mcl_adj = similar(mcl_adj)

    # do MCL iterations
    if display != :none
        @info("Starting MCL iterations...")
    end
    niter = 0
    converged = false
    rel_delta = NaN
    while !converged && niter < max_iter
        expanded = _mcl_expand(mcl_adj, expansion)
        _mcl_inflate!(next_mcl_adj, expanded, inflation)
        _mcl_prune!(next_mcl_adj, prune_tol)

        # normalize in columns
        rmul!(next_mcl_adj, Diagonal(map(x -> x != 0.0 ? 1.0/x : x,
                                         dropdims(sum(next_mcl_adj, dims=1), dims=1))))

        next_mcl_norm = norm(next_mcl_adj)
        if !isfinite(next_mcl_norm)
            @warn("MCL adjacency matrix norm is not finite")
            break
        end
        rel_delta = euclidean(next_mcl_adj, mcl_adj)/mcl_norm
        (display == :verbose) && @info("MCL iter. #$niter: rel.Δ=", rel_delta)
        (converged = rel_delta <= tol) && break
        # update (swap) MCL adjacency
        niter += 1
        mcl_adj, next_mcl_adj = next_mcl_adj, mcl_adj
        mcl_norm = next_mcl_norm
        (mcl_norm < tol) && break # matrix is zero
    end

    if display != :none
        if converged
            info("MCL converged after $niter iteration(s)")
        else
            warn("MCL didn't converge after $niter iteration(s)")
        end
    end

    (display == :verbose) && @info("Generating MCL clusters...")
    el2clu, clu_sizes, nunassigned = _mcl_clusters(mcl_adj, allow_singles,
                                                   tol/length(mcl_adj))

    return MCLResult(save_final_matrix ? mcl_adj : similar(mcl_adj, (0,0)),
              el2clu, clu_sizes, nunassigned, niter, rel_delta, converged)
end


text:
# Initially adapted from PiecewiseAffineTransforms.jl
# https://github.com/dfdx/PiecewiseAffineTransforms.jl

global MESHWARP_POLY2SOURCE_MASK = zeros(Bool, 0, 0);

"""
`MESHWARP` - Apply piecewise affine transform to image using bilinear interpolation

    warped_img, [bb.i, bb.j] = meshwarp(img, src, dst, trigs, offset, interp)

* `img`: 2D array, image (todo: extend to Image type)
* `src`: Nx2 array, positions of mesh nodes in `img` (global space)
* `dst`: Nx2 array, positions of mesh nodes in `warped_img` (global space)
* `trigs`: Nx3 array, each row contains indices of nodes of a triangle
* `offset`: 2-element array, position of img[1,1] in global space (optional - default is [0,0])
* `interp`: bool determining whether to use bilinear interpolation or not
    (optional - default is true)
* `warped_img`: with pixel values the same type as original image (for Int type, pixel values are rounded)
* `warped_offset`: 2-element array, position of warped_img[1,1] in global space

The bounding box of the warped image is defined as the smallest
integer-valued rectangle that contains the `dst` mesh.

This means that `warped_offset` is constrained to be integer-valued,
though `offset` is allowed to have floating point values.  The integer
constraint removes the need for further interpolation in any
subsequent fusing of multiple warped tiles.

See definitions in IMWARP documentation for further help.

"""
function meshwarp{T}(img::SharedArray{T},
                    src::Matrix{Float64}, dst::Matrix{Float64},
                    trigs::Matrix{Int64}, offset=[0,0], interp=true)

  bb = snap_bb(find_mesh_bb(dst))
  warped_img = SharedArray(T, bb.h+1, bb.w+1)
  warped_offset = [bb.i, bb.j];

  Us = Array{Any}(size(trigs, 1));
  Vs = Array{Any}(size(trigs, 1));
  Ms = Array{Array{Float64, 2}}(size(trigs, 1));

  @everywhere gc();

@fastmath @inbounds for t in 1:size(trigs, 1);
    #println("$t / $(size(trigs, 1))")
    tr = squeeze(trigs[t, :], 1)
    # coordinates of the source triangle
    X = src[tr, 1]
    Y = src[tr, 2]
    # coordinates of the destination triangle
    U = dst[tr, 1]
    V = dst[tr, 2]
    # Create matrix to transform from destination triangle to source image
    src_tri = [X Y ones(3,1)]
    dst_tri = [U V ones(3,1)]
    M = dst_tri \ src_tri # dst_tri * M = src_tri
    # Create list of coordinates in warped image that represent this
    # triangle (coordinates in the global space).
    Us[t] = U;
    Vs[t] = V;
    Ms[t] = M;
  end

  function proc_range(idx, arr::Array)
	worker_procs = setdiff(procs(), myid());
	nchunks = length(worker_procs);
	if nchunks == 0 return 1:length(arr); end
	if idx == myid() return 1:0; end
	splits = [round(Int64, s) for s in linspace(0, length(arr), nchunks + 1)];
	return splits[findfirst(worker_procs, idx)]+1:splits[findfirst(worker_procs, idx) + 1]
  end

#=@sync @fastmath @inbounds for t in 1:size(trigs, 1);
    @async remotecall_wait(procs()[rem(t, nprocs()-1)+2], calculate_pixels_in_trig!, Us[t], Vs[t], Ms[t], img, offset, warped_img, warped_offset)
end=#

@sync for p in procs()
      	t = proc_range(p, Us);
      	@async @fastmath @inbounds remotecall_wait(p, calculate_pixels_in_trig_chunk!, Us[t], Vs[t], Ms[t], img, offset, warped_img, warped_offset)
end

  return copy(sdata(warped_img)), [bb.i, bb.j]
end

function calculate_pixels_in_trig_chunk!(Us, Vs, Ms, img, offset, warped_img, warped_offset)
 @simd for i in 1:length(Us)
	calculate_pixels_in_trig!(Us[i], Vs[i], Ms[i], img, offset, warped_img, warped_offset)
  end
  global MESHWARP_POLY2SOURCE_MASK = zeros(Bool, 0, 0);
end

function calculate_pixels_in_trig!(U, V, M, img, offset, warped_img, warped_offset)
    us, vs = poly2source!(U, V)
    @simd for ind in 1:length(us)
        # Convert warped coordinate to pixel space
    @fastmath @inbounds begin
        i, j = us[ind]-warped_offset[1]+1, vs[ind]-warped_offset[2]+1
        # Use warped coordinate in global space for transform
        u, v = us[ind], vs[ind]
        # x, y = M * [u, v, 1]
        @fastmath @inbounds x, y = M[1,1]*u + M[2,1]*v + M[3,1], M[1,2]*u + M[2,2]*v + M[3,2]
        # Convert original image coordinate to pixel space
        x, y = x-offset[1]+1, y-offset[2]+1
        fx, fy = floor(Int64, x), floor(Int64, y)
        wx, wy = x-fx, y-fy
        if 1 <= fx <= size(img, 1)-1 && 1 <= fy <= size(img, 2)-1
          # Expansion of p = [1-wy wy] * img[fy:fy+1, fx:fx+1] * [1-wx; wx]
          @fastmath @inbounds p1 = ((1-wy)*img[fx,fy] + wy*img[fx,fy+1])
	  @fastmath @inbounds p2 = ((1-wy)*img[fx+1,fy] + wy*img[fx+1,fy+1])
	  @fastmath p1 = p1 * (1-wx);
	  @fastmath p2 = p2 * (wx);
          writepixel(warped_img,i,j,p1+p2)
        end
      end
    end
  end

"""
`POLY2SOURCE` - Run fillpoly and findn over the basic bounding box of a given triangle
"""
function poly2source(pts_i, pts_j)
  # Find bb of vertices (vertices in global space)
  top, bottom = floor(Int64,minimum(pts_i)), ceil(Int64,maximum(pts_i))
  left, right = floor(Int64,minimum(pts_j)), ceil(Int64,maximum(pts_j))
  # Create image based on number of pixels in bb that will identify triangle
  mask = zeros(Bool, bottom-top+1, right-left+1)
  # Convert vertices into pixel space and fill the mask to identify triangle
  fillpoly!(mask, pts_j-left+1, pts_i-top+1, true)
  # Create list of pixel coordinates that are contained by the triangle
  us, vs = findn(mask)
  # Convert that list of pixel coordinates back into global space
  us += top-1
  vs += left-1
  return us, vs
end

function poly2source!(pts_i, pts_j)
  # Find bb of vertices (vertices in global space)
  top, bottom = floor(Int64,minimum(pts_i)), ceil(Int64,maximum(pts_i))
  left, right = floor(Int64,minimum(pts_j)), ceil(Int64,maximum(pts_j))
  # Create image based on number of pixels in bb that will identify triangle
  if size(MESHWARP_POLY2SOURCE_MASK::Array{Bool, 2})[1] < bottom - top + 1 || size(MESHWARP_POLY2SOURCE_MASK)[2] < right - left + 1
    global MESHWARP_POLY2SOURCE_MASK = zeros(Bool, bottom-top+1, right-left+1)
  end
  (MESHWARP_POLY2SOURCE_MASK::Array{Bool, 2})[:] = false;
  #mask = zeros(Bool, bottom-top+1, right-left+1)
  # Convert vertices into pixel space and fill the mask to identify triangle
  for i in 1:length(pts_i)
    @fastmath @inbounds pts_i[i] = pts_i[i] + 1 - top;
  end
  for i in 1:length(pts_j)
    @fastmath @inbounds pts_j[i] = pts_j[i] + 1 - left;
  end
  #fillpoly!(mask, pts_j-left+1, pts_i-top+1, true)
   fillpoly!(MESHWARP_POLY2SOURCE_MASK::Array{Bool, 2}, pts_j, pts_i, true)
  # Create list of pixel coordinates that are contained by the triangle
  us, vs = findn(MESHWARP_POLY2SOURCE_MASK::Array{Bool, 2})
  # Convert that list of pixel coordinates back into global space
  for i in 1:length(us)
    @fastmath @inbounds us[i] = us[i] - 1 + top;
  end
  for i in 1:length(vs)
    @fastmath @inbounds vs[i] = vs[i] - 1 + left;
  end
#  us += top-1
#  vs += left-1
  return us, vs
end

"""
`FILLPOLY!` - Fill pixels contained inside a polygon

    M = fillpoly!(M, px, py, value)

* `M`: 2D array, image
* `px`: 1D array, x-components of polygon vertices
* `py`: 1D array, y-components of polygon vertices
* `value`: fill value

Features:

* Intended to work for nonconvex polygons as well as convex polygons, except that currently the corner cases (described below) fail in non-convex poligons since it doesn't add its own pair.
* The vertices should be listed sequentially in px and py. Clockwise/counterclockwise ordering doesn't matter.
* It seems that the polygon is closed automatically if it isn't already, i.e., the last vertex in px, py is not equal to the first.  Does the code work if the input polygon is already closed?
* The code treats the "edge case," where an edge of the polygon lies exactly on a grid line.

Bugs:

* The original code (https://github.com/dfdx/PiecewiseAffineTransforms.jl/blob/master/src/polyline.jl) used implicit conversion to Int64 (presumably rounding).  Tommy/Shang replaced this by ceil.  This might produce inconsistent results, with the same pixel belonging to more than one triangle.
* The "corner case" where a grid line intersects a single vertex of the polygon does not appear to be properly treated.  The corner case is nongeneric if px and py are floats.  But the corner case could be common if px and py are ints, which seems encouraged by the parametric typing.
"""
function fillpoly!{T,P<:Number}(M::Matrix{T}, px::Vector{P}, py::Vector{P}, value::T)
  @assert length(py) == length(px)
  left, right = floor(Int64,minimum(px)), ceil(Int64,maximum(px))
  # Scan poly from left to right
  for x=left:right     # loop over grid lines
    ys = Set{Int64}()
    m = length(px)
    for n=1:length(px)  # loop over edges (m,1),(1,2),(2,3),...,(m-1,m)
      # grid line intersects edge in one of two ways
      if (px[n] <= x <= px[m]) || (px[m] <= x <= px[n])
        if px[n] == px[m]  # intersection is entire edge
          push!(ys, ceil(Int64, py[n]))
          push!(ys, ceil(Int64, py[m]))
        else # intersection is point
            y = py[n] + (x-px[n]) * (py[m]-py[n])/(px[m]-px[n])
	    # deal with rounding error
	    if ceil(Int64, y) > size(M, 1)
	      push!(ys, floor(Int64, y))
	    else
              push!(ys, ceil(Int64, y))
	    end
        end
      end
      m = n
    end
    # generically, two intersections for a convex polygon
    # generically, even number of intersections for a nonconvex polygon
    ys = sort([y for y in ys])  # sort the intersection points
    # if odd number of intersection points, add duplicate end point
    if length(ys) % 2 == 1
      push!(ys, ys[end])
    end
    # Place value in matrix at all the y's between min and max for given x
    for n=1:2:length(ys)
      @simd for y in ys[n]:ys[n+1]
      @inbounds M[y, x] = value
      end
    end
  end
  return M
end

function meshwarp{T}(img::Array{T},
                    src::Matrix{Float64}, dst::Matrix{Float64},
                    trigs::Matrix{Int64}, offset=[0,0], interp=true)
  img = convert(SharedArray, img)
  meshwarp(img, src, dst, trigs, offset, interp)
end


text:
module AlgorithmicAllocation

using FillArrays
using IRTools
using LinearAlgebra
using MacroTools

using IRTools: @dynamo, IR, @code_ir, xcall, insertafter!, meta, blocks, varargs!
using IRTools.Inner: argnames!, update!
using MacroTools: prewalk, postwalk

function preallocate(::typeof(*), A::Matrix{<:Real}, B::VecOrMat{<:Real})
    C = A * B
    mem = deepcopy(C)
    function preallocated_mul(A::Matrix{<:Real}, B::VecOrMat{<:Real})
        return mul!(mem, A, B)
    end
    return C, preallocated_mul
end

function preallocate(::typeof(*), A::Matrix, B::Adjoint{<:Any, <:Matrix})
    C = A * B
    mem = deepcopy(C)
    function preallocated_mul(A::Matrix, B::Adjoint{<:Any, <:Matrix})
        return mul!(mem, A, B)
    end
    return C, preallocated_mul
end

function preallocate(::typeof(+), A::Array, B::Array)
    C = A + B
    mem = deepcopy(C)
    function preallocated_add(A::Array, B::Array)
        C .= A .+ B
        return C
    end
    return C, preallocated_add
end

function preallocate(::typeof(+), A::AbstractArray, B::Zeros)
    C = deepcopy(A)
    function preallocated_add_zero(A::Array, B::Zeros)
        C .= A
        return C
    end
    return A + B, preallocated_add_zero
end

preallocate(::Type{<:Adjoint}, A) = A', x -> x'

function preallocate(::typeof(Core.tuple), args...)
    return Core.tuple(args...), Core.tuple
end

# A really simple data structure that can be used with a pre-allocated function.
struct Preallocated{signature, T<:Tuple}
    preallocated_functions::T
end

function Preallocated{signature}(preallocated_functions::T) where {signature, T<:Tuple}
    return Preallocated{signature, T}(preallocated_functions)
end

function preallocated(signature, preallocated_functions::T) where {T}
    return Preallocated{signature, T}(preallocated_functions)
end

@generated function preallocate(f, args...)
    signature = Tuple{f, args...}
    metadata = meta(signature)
    ir = IR(metadata)
    preallocations = []

    # Modify all expressions to instead call preallocate.
    new_arg = IRTools.argument!(ir; at=1)
    for (x, st) in only(blocks(ir))
        isexpr(st.expr, :call) || continue

        # Insert a call to preallocate just before the function is evaluated.
        new_x = insert!(ir, x, xcall(Main, :preallocate, st.expr.args...))

        # Replace the original expression with the first element of the output of
        # preallocate, which is defined to be the same as just evaluating the function.
        ir[x] = xcall(:getindex, new_x, 1)

        # Insert an expression immediately after containing the preallocated function.
        new_preallocation = insertafter!(ir, x, xcall(:getindex, new_x, 2))

        # Store the location of the preallocated function for later.
        push!(preallocations, new_preallocation)
    end

    # Construct a `Preallocated` struct that runs the preallocated version of the function.
    preallocated_data = push!(ir, xcall(:tuple, preallocations...))
    new_allocation = push!(
        ir, xcall(Main, :preallocated, signature, preallocated_data),
    )

    # Set the function to return a tuple comprising the usual return value and the
    # Preallocated data structure.
    old_return = IRTools.returnvalue(only(IRTools.branches(ir)))
    new_return = push!(ir, xcall(:tuple, old_return, new_allocation))
    IRTools.return!(ir, new_return)

    # Hack the metadata to make everything. I don't understand this bit at all.
    argnames!(metadata, Symbol("#self#"), :f, :args)
    ir = varargs!(metadata, ir, 2)
    return update!(metadata.code, ir)
end

(p::Preallocated)(args...) = _preallocated(p::Preallocated, args...)

@generated function _preallocated(p::Preallocated{signature}, args...) where {signature}
    metadata = meta(signature)
    ir = IR(metadata)

    # Add an extra argument to the IR that is p.
    p_arg = IRTools.argument!(ir; at=2)

    # Pull the preallocated functions out of p.
    preallocated_functions = pushfirst!(ir, xcall(:getfield, p_arg, :(:preallocated_functions)))

    n = 0
    for (x, st) in only(blocks(ir))
        isexpr(st.expr, :call) || continue
        n += 1
        n == 1 && continue

        # Insert a call to preallocate just before the function is evaluated to get the
        # preallocated version of the function.
        new_x = insert!(ir, x, xcall(:getindex, preallocated_functions, n - 1))

        # preallocated_call = insert!(ir, x, xcall(:getindex, ))
        # Replace the original call with the pre-allocated call.
        ir[x] = xcall(new_x, st.expr.args[2:end]...)
    end

    # Again, hack the metadata to make everything. I don't understand this bit at all.
    argnames!(metadata, Symbol("#self#"), :p, :args)
    ir = varargs!(metadata, ir, 2)
    return update!(metadata.code, ir)
end

end


text:
####### VARYING SLOPES #######

#### Import dataframes of simulation results
# df = CSV.read("output_dfs/full_slope.csv", DataFrame)
# df2 = CSV.read("output_dfs/full_asymptote.csv", DataFrame)
# df3 = CSV.read("output_dfs/full_inflection.csv", DataFrame)

# In order to visualize our results as heatmaps (or a contour plots), we first split df by target and turn into 2D arrays:

# subset df by target rate
dfy = filter(:target => target -> target == "y", df)
dfB0 = filter(:target => target -> target == "B0", df)
dfr = filter(:target => target -> target == "r", df)
dfx = filter(:target => target -> target == "x", df)

# empty 2D arrays
DR_array_persistence_y = fill(NaN, length(doses), length(slopes)) #preallocate a 2d array for persistence values target rate = y
DR_array_persistence_B0 = fill(NaN, length(doses), length(slopes))
DR_array_persistence_r = fill(NaN, length(doses), length(slopes))
DR_array_persistence_x = fill(NaN, length(doses), length(slopes))

# Colour palette
palette2 = [RGB(190/255, 242/255, 170/255), RGB(149/255, 188/255, 219/255), RGB(30/255, 80/255, 150/255)]

#### rate = y ####
# fill 2D array for rate = y
for (j,s) in enumerate(slopes)
	for (i,d) in enumerate(doses)
		tmp = dfy[(dfy.slope .== s) .& (dfy.dose .== d),:]
		DR_array_persistence_y[i,j] = tmp.sp_survival[1]
	end
end

# plot for rate = y
p_sl_y = heatmap(string.(log2.(slopes)), string.(doses), DR_array_persistence_y
    , title  = "Persistence (Target = y)", classes = 3, color = palette2
    , xlabel = "log2(Slope)", ylabel ="Dose"
    )

#### rate = B0 ####
# fill 2D array for rate = B0
for (j,s) in enumerate(slopes)
	for (i,d) in enumerate(doses)
		tmp = dfB0[(dfB0.slope .== s) .& (dfB0.dose .== d),:]
		DR_array_persistence_B0[i,j] = tmp.sp_survival[1]
	end
end

# plot for rate = B0
p_sl_B0 = heatmap(string.(log2.(slopes)), string.(doses), DR_array_persistence_B0
    , title  = "Persistence (Target = B0)", classes = 3, color = palette2
    , xlabel = "log2(Slope)", ylabel ="Dose"
    )

#### rate = r ####
# fill 2D array for rate = r
for (j,s) in enumerate(slopes)
	for (i,d) in enumerate(doses)
		tmp = dfr[(dfr.slope .== s) .& (dfr.dose .== d),:]
		DR_array_persistence_r[i,j] = tmp.sp_survival[1]
	end
end

# plot for rate = r
p_sl_r = heatmap(string.(log2.(slopes)), string.(doses), DR_array_persistence_r
    , title  = "Persistence (Target = r)", classes = 3, color = palette2
    , xlabel = "log2(Slope)", ylabel ="Dose"
    )

#### rate = x ####
# fill 2D array for rate = x
for (j,s) in enumerate(slopes)
	for (i,d) in enumerate(doses)
		tmp = dfx[(dfx.slope .== s) .& (dfx.dose .== d),:]
		DR_array_persistence_x[i,j] = tmp.sp_survival[1]
	end
end

# plot for rate = x
p_sl_x = heatmap(string.(log2.(slopes)), string.(doses), DR_array_persistence_x
    , title  = "Persistence (Target = x)", classes = 3, color = palette2
    , xlabel = "log2(Slope)", ylabel ="Dose"
    )

l2 = @layout [a; b]
l4 = @layout [a b; c d]
plot(p_sl_y, p_sl_B0, p_sl_r, p_sl_x, layout = l4)

####### VARYING ASYMPTOTES #######
# subset df2 by target rate
df2y = filter(:target => target -> target == "y", df2)
df2B0 = filter(:target => target -> target == "B0", df2)
df2r = filter(:target => target -> target == "r", df2)
df2x = filter(:target => target -> target == "x", df2)

# empty 2D arrays
DR_array_persistence_y2 = fill(NaN, length(doses), length(asymptotes)) #preallocate a 2d array for persistence values target rate = y
DR_array_persistence_B02 = fill(NaN, length(doses), length(asymptotes))
DR_array_persistence_r2 = fill(NaN, length(doses), length(asymptotes))
DR_array_persistence_x2 = fill(NaN, length(doses), length(asymptotes))

#### rate = y ####
# fill 2D array for rate = y
for (j,s) in enumerate(asymptotes)
	for (i,d) in enumerate(doses)
		tmp = df2y[(df2y.asymptote .== s) .& (df2y.dose .== d),:]
		DR_array_persistence_y2[i,j] = tmp.sp_survival[1]
	end
end

# plot for rate = y
p_as_y = heatmap(string.(asymptotes), string.(doses), DR_array_persistence_y2
    , title  = "Persistence (Target = y)", classes = 3, color = palette2
    , xlabel = "Asymptote", ylabel ="Dose"
    )

#### rate = B0 ####
# fill 2D array for rate = B0
for (j,s) in enumerate(asymptotes)
	for (i,d) in enumerate(doses)
		tmp = df2B0[(df2B0.asymptote .== s) .& (df2B0.dose .== d),:]
		DR_array_persistence_B02[i,j] = tmp.sp_survival[1]
	end
end

# plot for rate = B0
p_as_B0 = heatmap(string.(asymptotes), string.(doses), DR_array_persistence_B02
    , title  = "Persistence (Target = B0)", classes = 3, color = palette2
    , xlabel = "Asymptote", ylabel ="Dose"
    )

#### rate = r ####
# fill 2D array for rate = r
for (j,s) in enumerate(asymptotes)
	for (i,d) in enumerate(doses)
		tmp = df2r[(df2r.asymptote .== s) .& (df2r.dose .== d),:]
		DR_array_persistence_r2[i,j] = tmp.sp_survival[1]
	end
end

# plot for rate = r
p_as_r = heatmap(string.(asymptotes), string.(doses), DR_array_persistence_r2
    , title  = "Persistence (Target = r)", classes = 3, color = palette2
    , xlabel = "Asymptote", ylabel ="Dose"
    )

#### rate = x ####
# fill 2D array for rate = x
for (j,s) in enumerate(asymptotes)
	for (i,d) in enumerate(doses)
		tmp = df2x[(df2x.asymptote .== s) .& (df2x.dose .== d),:]
		DR_array_persistence_x2[i,j] = tmp.sp_survival[1]
	end
end

# plot for rate = x
p_as_x = heatmap(string.(asymptotes), string.(doses), DR_array_persistence_x2
    , title  = "Persistence (Target = x)", classes = 3, color = palette2
    , xlabel = "Asymptote", ylabel ="Dose"
    )

####### VARYING INFLECTION POINTS #######
# subset df3 by target rate
df3y = filter(:target => target -> target == "y", df3)
df3B0 = filter(:target => target -> target == "B0", df3)
df3r = filter(:target => target -> target == "r", df3)
df3x = filter(:target => target -> target == "x", df3)

# empty 2D arrays
DR_array_persistence_y3 = fill(NaN, length(doses), length(inflections)) #preallocate a 2d array for persistence values target rate = y
DR_array_persistence_B03 = fill(NaN, length(doses), length(inflections))
DR_array_persistence_r3 = fill(NaN, length(doses), length(inflections))
DR_array_persistence_x3 = fill(NaN, length(doses), length(inflections))

#### rate = y ####
# fill 2D array for rate = y
for (j,s) in enumerate(inflections)
	for (i,d) in enumerate(doses)
		tmp = df3y[(df3y.xmid .== s) .& (df3y.dose .== d),:]
		DR_array_persistence_y3[i,j] = tmp.sp_survival[1]
	end
end

# plot for rate = y
p_xmid_y = heatmap(string.(inflections), string.(doses), DR_array_persistence_y3
    , title  = "Persistence (Target = y)", classes = 3, color = palette2
    , xlabel = "Inflection", ylabel ="Dose"
    )

#### rate = B0 ####
# fill 2D array for rate = B0
for (j,s) in enumerate(inflections)
	for (i,d) in enumerate(doses)
		tmp = df3B0[(df3B0.xmid .== s) .& (df3B0.dose .== d),:]
		DR_array_persistence_B03[i,j] = tmp.sp_survival[1]
	end
end

# plot for rate = B0
p_xmid_B0 = heatmap(string.(inflections), string.(doses), DR_array_persistence_B03
    , title  = "Persistence (Target = B0)", classes = 3, color = palette2
    , xlabel = "Inflection", ylabel ="Dose"
    )

#### rate = r ####
# fill 2D array for rate = r
for (j,s) in enumerate(inflections)
	for (i,d) in enumerate(doses)
		tmp = df3r[(df3r.xmid .== s) .& (df3r.dose .== d),:]
		DR_array_persistence_r3[i,j] = tmp.sp_survival[1]
	end
end

# plot for rate = r
p_xmid_r = heatmap(string.(inflections), string.(doses), DR_array_persistence_r3
    , title  = "Persistence (Target = r)", classes = 3, color = palette2
    , xlabel = "Inflection", ylabel ="Dose"
    )
#### rate = x ####
# fill 2D array for rate = x
for (j,s) in enumerate(inflections)
	for (i,d) in enumerate(doses)
		tmp = df3x[(df3x.xmid .== s) .& (df3x.dose .== d),:]
		DR_array_persistence_x3[i,j] = tmp.sp_survival[1]
	end
end

# plot for rate = x
p_xmid_x = heatmap(string.(inflections), string.(doses), DR_array_persistence_x3
    , title  = "Persistence (Target = x)", classes = 3, color = palette2
    , xlabel = "Inflection", ylabel ="Dose"
    )

max_stars_repo_name: hanamayall/simple-feeding-DR

max_issues_repo_name: hanamayall/simple-feeding-DR

max_forks_repo_name: hanamayall/simple-feeding-DR