text:
##Differentiation and integration

Base.sum{DD}(f::Fun{Laurent{DD}})=fouriersum(domain(f),f.coefficients)
Base.sum{DD}(f::Fun{Fourier{DD}})=fouriersum(domain(f),f.coefficients)

function linesum{DD}(f::Fun{Laurent{DD}})
    d=domain(f)
    if isa(d,Circle)
        sum(setcanonicaldomain(f))*d.radius
    else
        sum(f)
    end
end

linesum{DD<:Circle}(f::Fun{Fourier{DD}})=sum(setcanonicaldomain(f))*d.radius
linesum{DD<:PeriodicInterval}(f::Fun{Fourier{DD}})=sum(f) #TODO: Complex periodic interval

differentiate{DD<:PeriodicInterval}(f::Fun{Taylor{DD}}) = Fun(im*tocanonicalD(f,0)*taylor_diff(f.coefficients),f.space)
differentiate{DD<:PeriodicInterval}(f::Fun{Hardy{false,DD}}) = Fun(im*tocanonicalD(f,0)*hardyfalse_diff(f.coefficients),f.space)
differentiate{DD<:PeriodicInterval}(f::Fun{Laurent{DD}}) = Fun(im*tocanonicalD(f,0)*laurentdiff(f.coefficients),f.space)

differentiate{DD<:PeriodicInterval}(f::Fun{CosSpace{DD}}) = Fun(tocanonicalD(f,0)*cosspacediff(f.coefficients),SinSpace(domain(f)))
differentiate{DD<:PeriodicInterval}(f::Fun{SinSpace{DD}}) = Fun(tocanonicalD(f,0)*sinspacediff(f.coefficients),CosSpace(domain(f)))
differentiate{DD<:PeriodicInterval}(f::Fun{Fourier{DD}}) = Fun(tocanonicalD(f,0)*fourierdiff(f.coefficients),f.space)

differentiate{DD}(f::Fun{Laurent{DD}}) = Derivative(space(f))*f
differentiate{DD}(f::Fun{Fourier{DD}}) = Derivative(space(f))*f

function integrate{D}(f::Fun{Hardy{false,D}})
    if isa(domain(f),Circle) # drop -1 term if zero and try again
        @assert ncoefficients(f)==0 || abs(f.coefficients[1])<100eps()
        integrate(Fun(f,space(f)|(2:∞)))
    else  # Probably periodic itnerval
        Integral(space(f))*f
    end
end

function integrate{D}(f::Fun{Taylor{D}})
    if isa(domain(f),Circle)
        Integral(space(f))*f
    else  # Probably periodic itnerval  drop constant term if zero
        @assert ncoefficients(f)==0 || abs(f.coefficients[1])<100eps()
        Fun(integrate(Fun(f,space(f)|(2:∞))),space(f))
    end
end

function integrate{CS<:CosSpace}(f::Fun{CS})
    if isa(domain(f),Circle)
        error("Integrate not implemented for CosSpace on Circle")
    else  # Probably periodic itnerval, drop constant term if zero
        tol=1E-14 #TODO: smart tolerance.  Here relative is a bit tricky
                  # since this is called by Fourier integrate
        if abs(f.coefficients[1])<tol
            integrate(Fun(f,space(f)|(2:∞)))
        else
            d=domain(f)
            @assert isa(d,PeriodicInterval)
            x=Fun(identity,[first(d),last(d)])
            (f.coefficients[1]*x)⊕integrate(Fun(f,space(f)|(2:∞)))
        end
    end
end

function integrate{SS<:SinSpace}(f::Fun{SS})
    if isa(domain(f),Circle) # drop term containing z^(-1)
        integrate(Fun(f,space(f)|(2:∞)))
    else  # Probably periodic itnerval\
        Integral(space(f))*f
    end
end

#TODO: This is a hack to make sure Fourier maps to Fourier
# we don't have banded differentiate from CosSpace/SinSpace on a circle
for OP in (:differentiate,:integrate)
    @eval begin
        $OP{T,D<:PeriodicInterval}(f::Fun{Fourier{D},T})=$OP(f[2])⊕$OP(f[1])
        $OP{T,D<:Circle}(f::Fun{Fourier{D},T})=$OP(Fun(f,Laurent))
    end
end

fouriersum(d::PeriodicInterval,cfs)=cfs[1].*arclength(d)

function fouriersum{T}(d::Circle,cfs::Vector{T})
    if length(cfs)≥2
        2π*im*cfs[2]*d.radius
    else
        im*zero(T)
    end
end

# O(min(m,n)) Laurent line integral

function linebilinearform{T,D<:Circle}(f::Fun{Laurent{D},T},g::Fun{Laurent{D},T})
    @assert domain(f) == domain(g)
    u,v,mn = f.coefficients,g.coefficients,min(ncoefficients(f),ncoefficients(g))
    if mn > 1
        ret = u[1]*v[1]
        for i=2:2:mn-1
            ret += u[i]*v[i+1] + u[i+1]*v[i]
        end
        return arclength(domain(f))*ret
    elseif mn > 0
        return arclength(domain(f))*u[1]*v[1]
    else
        return zero(T)
    end
end

function bilinearform{T,D<:Circle}(f::Fun{Laurent{D},T},g::Fun{Laurent{D},T})
    @assert domain(f) == domain(g)
    u,v,mn = f.coefficients,g.coefficients,min(ncoefficients(f),ncoefficients(g))
    if mn > 2
        ret = u[1]*v[2] + u[2]*v[1]
        for i=3:2:mn-1
            ret += u[i]*v[i+1] + u[i+1]*v[i]
        end
        return complexlength(domain(f))*ret
    else
        return zero(T)
    end
end


text:
import Distances
include("structs.jl")

function stepped_dynamics(obj::Contraption, settings::PhysicsSettings, Δt::Float64)
    # Name binding for brevity
    bounds = settings.bounds
    n = length(obj.mass)

    # Storing kinematics
    s = deepcopy(obj.position)
    v = deepcopy(obj.velocity)
    a = zeros(2, n)

    # Compute acceleration from forces
    for i in 1:n
        # Gravity
        a[:, i] += [0., -settings.g]

        # Drag
        a[:, i] -= settings.drag * v[:, i]

        # Springs
        # TODO: should connections be memoized?
        for j in 1:n # Don't double-count!
            if obj.springs[i, j] != 0.
                Δs = s[:, j] - s[:, i]
                current_length = Distances.norm(Δs)
                direction = Δs / current_length

                a[:, i] += obj.springs[i,j] / obj.mass[i] *
                     (current_length - obj.rest_length[i,j]) *
                     direction # unit vector of force direction
            end
        end

    end

    # Applying acceleration
    for i in 1:n
        # Updating velocity from forces
        v[:, i] += a[:, i] * Δt

        # Continually attempt to apply velocity until result is in-bounds
        t_remaining = Δt
        while t_remaining > 0.

            # Movement
            proposed_position = s[:, i] + v[:, i] * Δt

            # Collision handling
            if !(bounds.x[1] <= proposed_position[1] <= bounds.x[2]) |
               !(bounds.y[1] <= proposed_position[2] <= bounds.y[2])

                # Determine whether horizontal or vertical collision is first
                # Catch 0 velocity cases
                if v[1, i] == 0.
                    x_or_y = :y

                    Δt_bottom = v[2, i] / (s[2, i] - bounds.y[1])
                    Δt_top    = v[2, i] / (s[2, i] - bounds.y[2])

                    Δt_c = max(Δt_bottom, Δt_top)
                elseif v[2, i] == 0.
                    x_or_y = :x

                    Δt_left  = v[1, i] / (s[1, i] - bounds.x[1])
                    Δt_right = v[1, i] / (s[1, i] - bounds.x[2])

                    Δt_c = max(Δt_left, Δt_right)
                else
                    Δt_left   = v[1, i] / (s[1, i] - bounds.x[1])
                    Δt_right  = v[1, i] / (s[1, i] - bounds.x[2])
                    Δt_bottom = v[2, i] / (s[2, i] - bounds.y[1])
                    Δt_top    = v[2, i] / (s[2, i] - bounds.y[2])

                    # Only one of top/bottom will be positive
                    # Select that one as the collision which occurs in the future
                    # Then select sooner of collisions that will occur
                    Δt_x = max(Δt_left, Δt_right)
                    Δt_y = max(Δt_top, Δt_bottom)

                    x_or_y = ifelse(Δt_x >= Δt_y, :x, :y)
                    Δt_c = min(Δt_x, Δt_y)
                end

                # Update proposal and time remaining in time step
                # Reverse movement to get to point where collision with wall occurs
                proposed_position -= v[:, i]*(Δt - Δt_c)
                t_remaining -= Δt_c

                # Change directions and lose energy
                if x_or_y == :x
                    v[1, i] *= -settings.elasticity
                else
                    v[2, i] *= -settings.elasticity
                end
            else
                s[:, i] = proposed_position
                t_remaining = 0
            end
        end
    end

    return (s, v)

end

ext: jl


text:
# x contains the refs
# P is the index of permutation, where and counts are integers vector of size ngroups + 1, perm is the permutaion of sorted pool, iperm is the inverse of perm
function _group_indexer!(x::Vector, original_P, copy_P, where, counts, lo, hi, ngroups::Integer, perm, iperm)
    # from PooledArrays.jl
    for i in 1:ngroups+1
        where[i] = 0
    end

    @inbounds for i = lo:hi
        where[x[i] + 1] += 1
    end
    counts[1] = 0
    for i in 1:ngroups
        counts[i+1] = where[perm[i] + 1]
    end
    # mark the start of each contiguous group of like-indexed data
    where[1] = 1
    @inbounds for i = 2:ngroups+1
        where[i] = where[i - 1] + counts[i - 1]
    end

    # this is our indexer
    @inbounds for i = lo:hi
        label = iperm[x[i]] + 1
        original_P[where[label] + lo - 1] = copy_P[i]
        where[label] += 1
    end

    # rearrange data
    @inbounds for i in 1:ngroups
        where[i] = 0
    end

    @inbounds for i = lo:hi
        where[x[i]] += 1
    end
    for i in 1:ngroups
        counts[i] = where[perm[i]]
    end

    f_indx = lo
    @inbounds for i = 1:ngroups
        l_indx = f_indx + counts[i] - 1
        val = perm[i]
        for j = f_indx:l_indx
            x[j] = val
        end
        f_indx = l_indx + 1
    end
end


text:
import Base: ==, <, >

"""
A record of a transition and the time.
It's sortable by time. Immutable.
"""
struct NRTransition
	key
	time::Float64
end

function Base.:<(a::NRTransition, b::NRTransition)
	a.time < b.time
end

function Base.:>(a::NRTransition, b::NRTransition)
    a.time > b.time
end

function ==(a::NRTransition, b::NRTransition)
    a.time == b.time
end

"""
Classic First Reaction method
"""
struct FirstReaction
end

function Next(fr::FirstReaction, system, rng)
	least=NRTransition(nothing, Inf)
	Hazards(system, rng) do clock, now, enabled, rng2
	  trial_time = Sample(clock.intensity, now, rng2)
	  @assert(trial_time >= now)
	  if trial_time < least.time
	  	least = NRTransition(clock, trial_time)
	  end
    end
    (least.time, least.key)
end

Observer(fr::FirstReaction) = (hazard, time, updated, rng) -> nothing


text:
using Knet.KnetArrays: KnetArray
param(x::KnetArray; atype=array_type[]) = Param(convert(atype,x))
update!(w::KnetArray, g, p) = gclip_update!(w, g, p)
_optimizers(::KnetArray,otype; o...) = otype(;o...)


text:
include("MLPClassifier.jl")
include("../datasets/Drug-Classification/dataset1.jl")
include("../datasets/Income-Classification/dataset2.jl")
include("../datasets/Wine-Classification/dataset3.jl")
include("../datasets/Drug-Consumism-Classification/dataset4.jl")
include("../datasets/Heart-Disease-Classification/dataset5.jl")

net = MLPClassifier(dataset1_nfeatures, dataset1_labels, [15,5,10])
#net = MLPClassifier(dataset1_nfeatures, dataset1_labels, [5])
x_train, y_train, x_test, y_test = get_dataset1("../datasets/Drug-Classification")
println("Accuracy before training: ", accuracy(zip(x_test, y_test), net))
println("Error rate before training: ", error_rate(zip(x_test, y_test), net))
#fit!(net, x_train, y_train; max_epochs = 1000)
fit!(net, x_train, y_train; max_epochs = 10000)
println("Accuracy after training: ", accuracy(zip(x_test, y_test), net))
println("Error rate after training: ", error_rate(zip(x_test, y_test), net))
#=
-> 5 classes
Accuracy before training: 0.05
Error rate before training: 0.95

Accuracy after training: 0.5
Error rate after training: 0.5
=#

#=
net = MLPClassifier(dataset2_nfeatures, dataset2_labels, [26,13,6])
x_train, y_train, x_test, y_test = get_dataset2("../datasets/Income-Classification")
println("Accuracy before training: ", accuracy(zip(x_test, y_test), net))
println("Error rate before training: ", error_rate(zip(x_test, y_test), net))
#fit!(net, x_train, y_train; max_epochs = 5000)
fit!(net, x_train, y_train; max_epochs = 10000)
println("Accuracy after training: ", accuracy(zip(x_test, y_test), net))
println("Error rate after training: ", error_rate(zip(x_test, y_test), net))
=#
#=
-> 2 classes
Accuracy before training: 0.24888684170121295
Error rate before training: 0.751113158298787

Accuracy after training: 0.7543374788883771
Error rate after training: 0.2456625211116229
=#

#=
net = MLPClassifier(dataset3_nfeatures, dataset3_labels, [22,10,15,11])
x_train, y_train, x_test, y_test = get_dataset3("../datasets/Wine-Classification")
println("Accuracy before training: ", accuracy(zip(x_test, y_test), net))
println("Error rate before training: ", error_rate(zip(x_test, y_test), net))
#fit!(net, x_train, y_train; max_epochs = 15000)
fit!(net, x_train, y_train; max_epochs = 10000)
println("Accuracy after training: ", accuracy(zip(x_test, y_test), net))
println("Error rate after training: ", error_rate(zip(x_test, y_test), net))
=#
#=
-> 11 classes
Accuracy before training: 0.007692307692307693
Error rate before training: 0.9923076923076923

Accuracy after training: 0.5030769230769231
Error rate after training: 0.4969230769230769
=#

#=
#net = MLPClassifier(dataset4_nfeatures, dataset4_labels, [24,6,12,7])
net = MLPClassifier(dataset4_nfeatures, dataset4_labels, [24,6,12,4,7])
x_train, y_train, x_test, y_test = get_dataset4("../datasets/Drug-Consumism-Classification")
println("Accuracy before training: ", accuracy(zip(x_test, y_test), net))
println("Error rate before training: ", error_rate(zip(x_test, y_test), net))
#fit!(net, x_train, y_train; max_epochs = 20000)
fit!(net, x_train, y_train; max_epochs = 10000)
println("Accuracy after training: ", accuracy(zip(x_test, y_test), net))
println("Error rate after training: ", error_rate(zip(x_test, y_test), net))
=#
#=
-> 7 classes
Accuracy before training: 0.17506631299734748
Error rate before training: 0.8249336870026526

Accuracy after training: 0.3952254641909814
Error rate after training: 0.6047745358090185
=#

#=
net = MLPClassifier(dataset5_nfeatures, dataset5_labels, [26, 13, 7, 5])
x_train, y_train, x_test, y_test = get_dataset5("../datasets/Heart-Disease-Classification")
println("Accuracy before training: ", accuracy(zip(x_test, y_test), net))
println("Error rate before training: ", error_rate(zip(x_test, y_test), net))
fit!(net, x_train, y_train; max_epochs = 10000)
println("Accuracy after training: ", accuracy(zip(x_test, y_test), net))
println("Error rate after training: ", error_rate(zip(x_test, y_test), net))
=#
#=
-> 5 classes
Accuracy before training: 0.03333333333333333
Error rate before training: 0.9666666666666667

Accuracy after training: 0.48333333333333334
Error rate after training: 0.5166666666666667
=#

max_stars_repo_name: Lucas-Muniz/MNEAT

max_issues_repo_name: Lucas-Muniz/MNEAT

max_forks_repo_name: Lucas-Muniz/MNEAT


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 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

Random.seed!(9)

Nx=64
im = readdust(Nx)
true_img = im[:, :, 1]

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

dhc_args = Dict(:doS2=>false, :doS12=>false, :doS20=>true, :apodize=>true) #CHANGE: Change sig for the sfd data since the noise model is super high and the tiny values make sense
white_noise_args = Dict(:loc=>0.0, :sig=>1.0, :Nsam=>1000, :iso=>false, :norm=>false, :smooth=>true, :smoothval=>0.8) #only if you're using noise based covar
optim_settings = Dict([("iterations", 20), ("norm", false)])
recon_settings = Dict([("target_type", "ground truth"), ("covar_type", "white_noise"), ("log", false), ("GaussianLoss", true), ("Invcov_matrix", "Diagonal+Eps"),
  ("optim_settings", optim_settings), ("lambda", 0.0), ("white_noise_args", white_noise_args)]) #Add constraints
#"white_noise_args": args for white_noise sims if GaussianLoss and Covar_type='white_noise'. :loc...

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

recon_img = reconstruction_wrapper(true_img, init, filter_hash, dhc_args, coeff_mask, recon_settings)

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

dhc_args_nonapd = Dict(:doS2=>false, :doS12=>false, :doS20=>true, :apodize=>false)
dhc_args_apd = dhc_args
@benchmark s_apd = DHC_compute_apd(true_img, filter_hash; dhc_args_apd...)
@benchmark s_nonapd = DHC_compute_apd(true_img, filter_hash; dhc_args_nonapd...)

heatmap(true_img)
heatmap(init)

minimum(true_img)

function adaptive_apodizer(input_image::Array{Float64, 2}, dhc_args)#SPEED
        if dhc_args[:apodize]
            apd_image = apodizer(input_image)
        else
            apd_image = input_image
        end
        return apd_image
end

function adaptive_apodizer_fast(input_image::Array{Float64, 2}, dhc_args)#SPEED
        if dhc_args[:apodize]
            input_image .= apodizer(input_image)
        else

        end
        return input_image
end

@benchmark apdval = adaptive_apodizer(true_img, dhc_args_nonapd)
@benchmark apdvalfast = adaptive_apodizer_fast(true_img, dhc_args_nonapd)
#No difference

@benchmark apdimg = apodizer(true_img) #39 μs
@benchmark img = DHC_compute_apd(true_img, filter_hash; dhc_args_nonapd...) #3.475 ms
@benchmark img = DHC_compute_apd(true_img, filter_hash; dhc_args_apd...) #3.595 ms

Profile.clear()
@profile img = DHC_compute_apd(true_img, filter_hash; dhc_args_apd...)
Juno.profiler()

using StatsBase
arr1 = rand(16, 16)
wts = fweights(arr1)

#3/21###########################################################################

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

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

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

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

function dapodizer_test(img)
    eps = 1e-4
    im1 = copy(img)
    im0 = copy(img)
    im1[6, 30] += (eps/2)
    im0[6, 30] -= (eps/2)
    brute = ((apodizer(im1) - apodizer(im0))./ eps)[6, 30]
    Apmat = wind_2d(size(img)[1])
    c = sum(Apmat)
    ct = 1.0 + (1.0/c)
    term2 = (Apmat./c)
    Apmat .*= (ones(size(Apmat)).*ct .- term2)
    clever = Apmat[6, 30]
    println("Brute ", brute)
    println("Clever ", clever)
    println("diff= ", clever - brute)
end

function dwtapdtest(img, filter_hash)
    (Nf, ) = size(filter_hash["filt_index"])
    eps = 1e-4
    im1 = copy(img)
    im0 = copy(img)
    im1[6, 30] += (eps/2)
    im0[6, 30] -= (eps/2)
    brute = ((apodizer(im1) - apodizer(im0))./ eps)[6, 30]
    Apmat = wind_2d(size(img)[1])
    c = sum(Apmat)
    ct = 1.0 + (1.0/c)
    term2 = (Apmat./c)
    Apmat .*= (ones(size(Apmat)).*ct .- term2)
    clever = Apmat[6, 30]
    println("Brute ", brute)
    println("Clever ", clever)
    println("diff= ", clever - brute)

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

    sum1  = Deriv_Utils_New.wst_S20_deriv_sum(apodizer(img), filter_hash, wtgrid, FFTthreads=1)
    Sjac = reshape(Deriv_Utils_New.wst_S20_deriv(apodizer(img), filter_hash, FFTthreads=1), Nx,Nx,Nf*Nf)
    sum2  = zeros(Float64, Nx, Nx)
    for i=1:Nf*Nf sum2 += (Sjac[:,:,i].*wtvec[i]) end #THIS NEEDS TO BE WTGRID
    sum2 = reshape(sum2, (Nx^2, 1))
    println("Abs mean", mean(abs.(sum2 - sum1)))
    println("Abs mean", std(sum2 - sum1))
    println("Ratio: ", mean(sum2./sum1))
end

dapodizer_test(true_img)
dS20sum_test(adaptive_apodizer(true_img), filter_hash)
dwtapdtest(true_img, filter_hash)

#######DEBUG ZONE##############
tonorm = false
s_targ_mean, s_targ_invcov = s2mean, s2icov
function adaptive_apodizer(input_image::Array{Float64, 2})
    if dhc_args[:apodize]
        apd_image = apodizer(input_image)
    else
        apd_image = input_image
    end
    return apd_image
end
#For the derivative term dA*P/dP
if dhc_args[:apodize]
    #=
    Apmat = wind_2d(Nx)
    c = sum(Apmat)
    ct = 1.0 + (1.0/c)
    term2 = (Apmat./c)
    Apmat = Apmat .* (ct .- term2)
    =#
    Ap = wind_2d(Nx)
    cA = sum(Ap)
    Apflat = reshape(Ap, Nx^2)
    od_nz_idx = findall((x->((x!=0) & (x!=1))), Apflat)
    #tol = 1e-3
    #od_nz_idx2 = findall((x->((abs(x)>tol) & (abs(x-1)>tol))), Apflat)
    avals = Apflat[od_nz_idx]
    dA = zeros(Nx^2, Nx^2)
    dA[od_nz_idx, od_nz_idx] .= ((1.0 .- avals) * avals')./cA
    dA[diagind(dA)] += Apflat

else
    dA = I
end

lambda=0.0
FFTthreads=1
pixmask = falses(size(init))

s_targ_invcov = I
#DEB: input->init
function loss_func20(img_curr::Array{Float64, 2})
    s_curr = DHC_compute_apd(img_curr,  filter_hash, iso=false, norm=tonorm; dhc_args...)[coeff_mask]
    #regterm =  0.5*lambda*sum((adaptive_apodizer(img_curr) - adaptive_apodizer(init)).^2)
    lnlik = ( 0.5 .* (s_curr - s_targ_mean)' * s_targ_invcov * (s_curr - s_targ_mean))
    #println("Lnlik size | Reg Size", size(lnlik), size(regterm))
    neglogloss = lnlik[1] #+ regterm #BUG:??
    return s_curr - s_targ_mean, lnlik, neglogloss
end
#TODO: Need to replace the deriv_sums20 with a deriv_sum wrapper that handles all cases (S12, S20, S2)
function dloss20(storage_grad::Array{Float64, 2}, img_curr::Array{Float64, 2})
    s_curr = DHC_compute_apd(img_curr, filter_hash, norm=tonorm, iso=false; dhc_args...)[coeff_mask]
    diff = s_curr - s_targ_mean
    wt = reshape(convert(Array{Float64, 2}, transpose(diff) * s_targ_invcov), (Nf, Nf))
    apdimg_curr = adaptive_apodizer(img_curr)
    storage_grad .= reshape(Deriv_Utils_New.wst_S20_deriv_sum(apdimg_curr, filter_hash, wt, FFTthreads=FFTthreads)' * dA, (Nx, Nx)) #DEBUG
    #dnll_wrong = reshape(Deriv_Utils_New.wst_S20_deriv_sum(wind_2d(Nx).*img_curr, filter_hash, wt, FFTthreads=FFTthreads), (Nx, Nx)) #DEBUG
    #storage_grad .= Apmat .* (dnll + lambda.*(apdimg_curr - adaptive_apodizer(init)))
    storage_grad[pixmask] .= 0 # better way to do this by taking pixmask as an argument wst_s2_deriv_sum?
    return diff, wt, dA, dnll#, dnll_wrong
end

println("Diff check")
eps = zeros(size(init))
row, col = 24, 18 #convert(Int8, Nx/2), convert(Int8, Nx/2)+3
eps[row, col] = 1e-4
d1, ll1, chisq1 = loss_func20(init+eps./2) #DEB
d0, ll0, chisq0 = loss_func20(init-eps./2) #DEB
d, ll, chisq = loss_func20(init)
brute  = (chisq1-chisq0)/1e-4
#df_brute = DHC_compute(reshape(input, (Nx, Nx)), filter_hash, doS2=true, doS12=false, doS20=false, norm=false)[coeff_mask] - s_targ_mean
clever = zeros(size(init)) #DEB
meanval = mean(adaptive_apodizer(init)) ./ mean(wind_2d(Nx))
ddiff, dwt, dApmat, dnll = dloss20(clever, init) #DEB
println("Chisq Derve Check")
println("Brute:  ",brute)
println("Clever: ",clever[row, col], " Difference: ", brute - clever[row, col], " Mean ", meanval) #DEB

#diffs: same
brute_dnll = ((ll1 - ll0)/1e-4)[1]
dl_dnll = (Apmat .* dnll)[row, col]
brute_dnll - dl_dnll

dl_dnll_wrong = (Apmat .* dnll_wrong)[row, col]
brute_dnll - dl_dnll_wrong

#heatmap(Apmat)
heatmap(Apmat - wind_2d(Nx))

#3-22-23#############################
#Examining dA
Nx=64
Ap = wind_2d(Nx)
cA = sum(Ap)
true_img = readdust(Nx)
filter_hash = fink_filter_hash(1, 8, nx=Nx, pc=1, wd=1, Omega=true)
Apflat = reshape(Ap, Nx^2)
od_nz_idx = findall((x->((x!=0) & (x!=1))), Apflat)
#tol = 1e-3
#od_nz_idx2 = findall((x->((abs(x)>tol) & (abs(x-1)>tol))), Apflat)
avals = Apflat[od_nz_idx]
dA = zeros(Nx^2, Nx^2)
dA[od_nz_idx, od_nz_idx] .= ((1.0 .- avals) * avals')./cA
dA[diagind(dA)] += Apflat
#dA[od_nz_idx, od_nz_idx][diagind(dA[od_nz_idx, od_nz_idx])] .+= avals
count((x->x!=0), dA) / length(dA)

spdA = sparse(dA)

mywts = reshape(rand(Nf*Nf), (Nf, Nf))
dls20 = Deriv_Utils_New.wst_S20_deriv_sum(true_img, filter_hash, mywts, FFTthreads=FFTthreads)'
@benchmark prod = dls20 * dA #21ms mem: 32KiB, allocs=2
@benchmark prod = dls20 * spdA #4ms mem / allocs same
#Conclusion: Sparse is definitely faster (factor of 5)
dA[od_nz_idx, od_nz_idx]
#check
xi = 12
xj = 500
dA[od_nz_idx, od_nz_idx][xj, xi]
dA[od_nz_idx, od_nz_idx][xi, xj]
dA[od_nz_idx, od_nz_idx][xj, xi]
axi = Apflat[od_nz_idx][xi]
axj = Apflat[od_nz_idx][xj]
dA[od_nz_idx, od_nz_idx][xi, xi] - axi*(1 + (1/cA) -(axi/cA))
dA[od_nz_idx, od_nz_idx][xj, xi] - axi*(1-axj)/cA
dA[od_nz_idx, od_nz_idx][xi, xj] - axj*(1-axi)/cA

#Difference STILL EXISTS
Nx=64
Ap = wind_2d(Nx)
cA = sum(Ap)
true_img = readdust(Nx)
filter_hash = fink_filter_hash(1, 8, nx=Nx, pc=1, wd=1, Omega=true)
function dwtapdtest_single(img, filter_hash)
    (Nf, ) = size(filter_hash["filt_index"])
    eps = 1e-4
    im1 = copy(img)
    im0 = copy(img)
    row, col = 20, 39
    flat_ind = Nx*(col-1) + row
    im1[row, col] += (eps/2)
    im0[row, col] -= (eps/2)
    brute = ((apodizer(im1) - apodizer(im0))./ eps)#[row, col]
    #=Apmat = wind_2d(size(img)[1])
    c = sum(Apmat)
    ct = 1.0 + (1.0/c)
    term2 = (Apmat./c)
    Apmat .*= (ones(size(Apmat)).*ct .- term2)
    =#
    Ap = wind_2d(Nx)
    cA = sum(Ap)
    Apflat = reshape(Ap, Nx^2)
    od_nz_idx = findall(!iszero, Apflat) #findall((x->((x!=0) & (x!=1))), Apflat)
    #od_zidx = findall(iszero, Apflat)
    avals = Apflat[od_nz_idx]
    dA = zeros(Nx^2, Nx^2)
    dA[:, od_nz_idx] .= ((1.0 .- Apflat) * avals')./cA
    dA[diagind(dA)] += Apflat
    #tol = 1e-3
    #od_nz_idx2 = findall((x->((abs(x)>tol) & (abs(x-1)>tol))), Apflat)
    #avals = Apflat[od_nz_idx]
    #
    #dA[od_nz_idx, od_nz_idx] .= ((1.0 .- avals) * avals')./cA
    println("Check dA[:, xi]")
    clever = reshape(dA[:, flat_ind], (Nx, Nx))
    #println("Brute ", brute)
    #println("Clever ", clever)
    println("Ratio= ", mean(abs.(clever./brute)))
    println("Mean Abs diff= ", mean(abs.(clever - brute)))

    # this is symmetrical, but not because S20 is symmetrical!
    mywts = rand(Nf, Nf)
    wtgrid = reshape(mywts, Nf, Nf) + reshape(mywts, Nf, Nf)'
    #only one coeff
    #wtgrid = zeros(Nf, Nf)
    #wtgrid[1, 1] = 1.0
    wtvec  = reshape(wtgrid, Nf*Nf)

    sum1  = Deriv_Utils_New.wst_S20_deriv_sum(apodizer(img), filter_hash, wtgrid, FFTthreads=1) #BUG: Is apodizer image doing something weird??
    Sjac = reshape(Deriv_Utils_New.wst_S20_deriv(apodizer(img), filter_hash, FFTthreads=1), Nx,Nx,Nf*Nf)
    sum2  = zeros(Float64, Nx, Nx)
    for i=1:Nf*Nf sum2 += (Sjac[:,:,i].*wtvec[i]) end #THIS NEEDS TO BE WTGRID
    sum2 = reshape(sum2, (Nx^2, 1))
    println("Check wst_deriv_sum vs longer route")
    println("Abs mean", mean(abs.(sum2 - sum1)))
    println("Abs mean", std(sum2 - sum1))
    println("Ratio: ", mean(sum2./sum1))

    println("Check deriv of Weighted Sum{S(A*P)}")
    dll = reshape(sum1' * dA, (Nx, Nx))
    coeffmask = falses(2+Nf+Nf^2)
    coeffmask[Nf+3:end] .= true
    Sh = DHC_compute_apd(im1, filter_hash; doS2=false, doS20=true, doS12=false, apodize=true, norm=false, iso=false)[coeffmask]
    Su =  DHC_compute_apd(im0, filter_hash; doS2=false, doS20=true, doS12=false, apodize=true, norm=false, iso=false)[coeffmask]

    dll_brute = ((wtvec' * Sh) - (wtvec' * Su))./eps
    println("Diff", dll[row, col] - dll_brute)
    println("Ratio", dll[row, col]./dll_brute)
    #dll_brute = (Deriv_Utils_New.wst_S20_deriv_sum(apodizer(im1), filter_hash, wtgrid, FFTthreads=1) - Deriv_Utils_New.wst_S20_deriv_sum(apodizer(im0), filter_hash, wtgrid, FFTthreads=1))./eps
    return clever, brute, dll, dll_brute, dA
end

dAclev, dAbrute, dll_cl, dll_brute, dAjac = dwtapdtest_single(true_img, filter_hash)
heatmap(reshape(dll_cl, (Nx, Nx)))
heatmap(reshape(dll_brute, (Nx, Nx)))
heatmap(apodizer(true_img))

row, col=6, 30
flat_ind = Nx*(col-1) + row
println(dA[1, flat_ind])

maximum(dll_cl)
maximum(dll_brute)

#TODO:
#Replace with fixed mu and look at comparison
#Why is the brute force test and clever SO DIFFERENT?
#shifting>????
#not using coeffmask in one place but using it elsewhere??/
#Revisit window func

Ap = wind_2d(Nx)
cA = sum(Ap)
Apflat = reshape(Ap, Nx^2)
od_nz_idx = findall(!iszero, Apflat) #findall((x->((x!=0) & (x!=1))), Apflat)
#od_zidx = findall(iszero, Apflat)
avals = Apflat[od_nz_idx]
dA = zeros(Nx^2, Nx^2)
dA[:, od_nz_idx] .= ((1.0 .- Apflat) * avals')./cA
dA[diagind(dA)] += Apflat

count((x->x!=0), dA) / length(dA)

spdA = sparse(dA)

mywts = reshape(rand(Nf*Nf), (Nf, Nf))
dls20 = Deriv_Utils_New.wst_S20_deriv_sum(true_img, filter_hash, mywts, FFTthreads=1)'
@benchmark prod = dls20 * dA #22ms mem: 32KiB, allocs=2
@benchmark prod = dls20 * spdA #12ms mem / allocs same

#3/24###########################################################################
#test wrapper
Nx=16
fhash = fink_filter_hash(1, 8, nx=Nx, pc=1, wd=1, Omega=true)
(Nf, ) = size(fhash["filt_index"])
coeff_mask = falses(2+Nf+Nf^2)
coeff_mask[Nf+3:end] .= true
function test_wrapper(Nx, filter_hash, coeff_mask)
    true_img = readdust(Nx)
    swrap = DHC_compute_wrapper(true_img, filter_hash, coeff_mask=coeff_mask)
    sapd = DHC_compute_apd(true_img, filter_hash)[coeff_mask]
    print("Mean Abs Diff ", mean(abs.(swrap - sapd)))
end

test_wrapper(16, fhash, coeff_mask)

#=
filter_hash = fink_filter_hash(1, 8, nx=64, pc=1, wd=1, Omega=true)
Nx=64
true_img = readdust(Nx)
filter_hash = fink_filter_hash(1, 8, nx=Nx, pc=1, wd=1, Omega=true)
DHC_compute(true_img, filter_hash, filter_hash2=filter_hash)

function test(a, b, c=b; d=3)
    return (a/b) -c + d
end

test(3, 4, 2, d=2)
=#

##3-28: Adding Iso and propagating changes everywhere
Nx=16
im = readsfd(Nx)
true_img = im[:, :, 1]
fhash = fink_filter_hash(1, 8, nx=Nx, pc=1, wd=1, Omega=true)
siso = DHC_compute(true_img, fhash, iso=true, norm=false, doS2=true)

#S20+NoIsoChecking that you didnt break old code
fname_save = "scratch_NM/NewWrapper/IsoTests/s20_noiso"
Nx=64
im = readsfd(Nx)
true_img = im[:, :, 1]

noise = rand(Normal(0.0, std(true_img)), Nx, Nx)
init = true_img + noise
#heatmap(init)
#heatmap(true_img)

filter_hash = fink_filter_hash(1, 8, nx=Nx, pc=1, wd=1, Omega=true)
(Nf, ) = size(filter_hash["filt_index"])
coeff_mask = falses(2+Nf+Nf^2)
coeff_mask[Nf+3:end] .= true

dhc_args = Dict(:doS2=>false, :doS12=>false, :doS20=>true, :apodize=>true, :iso=>false) #Iso #CHANGE: Change sig for the sfd data since the noise model is super high and the tiny values make sense
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", 10), ("norm", false), ("minmethod", ConjugateGradient())])
recon_settings = Dict([("target_type", "ground_truth"), ("covar_type", "white_noise"), ("log", false), ("GaussianLoss", true), ("Invcov_matrix", "Diagonal+Eps"),
  ("optim_settings", optim_settings), ("lambda", std(apodizer(true_img)).^(-2)), ("white_noise_args", white_noise_args)]) #Add constraints
regs_true = DHC_compute_wrapper(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=regs_true)
#Check
#s2mean - s_true

recon_settings["s_targ_mean"] = s2mean
recon_settings["s_invcov"] = s2icov
recon_settings["safemode"] = true
recon_settings["fname_save"] = fname_save * ".jld2"
recon_settings["optim_settings"] = optim_settings
resobj, recon_img = reconstruction_wrapper(true_img, init, filter_hash, dhc_args, coeff_mask, recon_settings)

p= plot([t.value for t in Optim.trace(resobj)])
plot!(title="Loss: S20 | Targ = S(GT) | Cov = SFD Covariance", xlabel = "No. Iterations", ylabel = "Loss")
savefig(p, fname_save * "_trace.png")

Visualization.plot_diffscales([apodizer(true_img), apodizer(init), apodizer(recon_img), apodizer(recon_img) - apodizer(true_img)], ["GT", "Init", "Reconstruction", "Residual"], fname=fname_save*".png")
Visualization.plot_synth_QA(apodizer(true_img), apodizer(init), apodizer(recon_img), fname=fname_save*"_6p.png")

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)))

#S20 + Iso check
#S20+NoIsoChecking that you didnt break old code
fname_save = "scratch_NM/NewWrapper/3-29/GTTargWNCov/log_apd_iso"
Nx=64
im = readsfd(Nx)
true_img = im[:, :, 1]

noise = rand(Normal(0.0, std(true_img)), Nx, Nx)
init = true_img + noise
#heatmap(init)
#heatmap(true_img)

filter_hash = fink_filter_hash(1, 8, nx=Nx, pc=1, wd=1, Omega=true)
(S1iso, Nf) = size(filter_hash["S1_iso_mat"])
(S2iso, Nf) = size(filter_hash["S2_iso_mat"])
coeff_mask = falses(2+S1iso+S2iso)
coeff_mask[S1iso+3:end] .= true
#coeff_mask = falses(2+Nf+Nf^2)
#coeff_mask[Nf+3:end] .= true

dhc_args = Dict(:doS2=>false, :doS12=>false, :doS20=>true, :apodize=>true, :iso=>true) #Iso #CHANGE: Change sig for the sfd data since the noise model is super high and the tiny values make sense
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", "ground_truth"), ("covar_type", "white_noise"), ("log", true), ("GaussianLoss", true), ("Invcov_matrix", "Diagonal+Eps"),
  ("optim_settings", optim_settings), ("lambda", std(apodizer(true_img)).^(-2)), ("white_noise_args", white_noise_args)]) #Add constraints
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=regs_true)
#Check
#s2mean - s_true

recon_settings["s_targ_mean"] = s2mean
recon_settings["s_invcov"] = s2icov
recon_settings["safemode"] = false
recon_settings["fname_save"] = fname_save * ".jld2"
recon_settings["optim_settings"] = optim_settings
resobj, recon_img = reconstruction_wrapper(true_img, init, filter_hash, dhc_args, coeff_mask, recon_settings)

p= plot([t.value for t in Optim.trace(resobj)])
plot!(title="Loss: S20 | Targ = S(GT) | Cov = SFD Covariance", xlabel = "No. Iterations", ylabel = "Loss")
savefig(p, fname_save * "_trace.png")

Visualization.plot_diffscales([apodizer(true_img), apodizer(init), apodizer(recon_img), apodizer(recon_img) - apodizer(true_img)], ["GT", "Init", "Reconstruction", "Residual"], fname=fname_save*".png")
Visualization.plot_synth_QA(apodizer(true_img), apodizer(init), apodizer(recon_img), fname=fname_save*"_6p.png")

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)))

#S20 + Iso check
#S20+NoIsoChecking that you didnt break old code
fname_save = "scratch_NM/NewWrapper/3-29/SFDTargSFDCov/log_apd_iso"
Nx=64
im = readsfd(Nx)
true_img = im[:, :, 1]

noise = rand(Normal(0.0, std(true_img)), Nx, Nx)
init = true_img + noise
#heatmap(init)
#heatmap(true_img)

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=>true) #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", true), ("GaussianLoss", true), ("Invcov_matrix", "Diagonal+Eps"),
  ("optim_settings", optim_settings), ("lambda", std(apodizer(log.(true_img))).^(-2)), ("white_noise_args", white_noise_args)]) #Add constraints
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)
#Check
#s2mean - s_true

recon_settings["s_targ_mean"] = s2mean
recon_settings["s_invcov"] = s2icov
recon_settings["safemode"] = false
recon_settings["fname_save"] = fname_save * ".jld2"
recon_settings["optim_settings"] = optim_settings
resobj, recon_img = reconstruction_wrapper(true_img, init, filter_hash, dhc_args, coeff_mask, recon_settings)

p= plot([t.value for t in Optim.trace(resobj)])
plot!(title="Loss: S20 | Targ = S(GT) | Cov = SFD Covariance", xlabel = "No. Iterations", ylabel = "Loss")
savefig(p, fname_save * "_trace.png")

Visualization.plot_diffscales([apodizer(true_img), apodizer(init), apodizer(recon_img), apodizer(recon_img) - apodizer(true_img)], ["GT", "Init", "Reconstruction", "Residual"], fname=fname_save*".png")
Visualization.plot_synth_QA(apodizer(true_img), apodizer(init), apodizer(recon_img), fname=fname_save*"_6p.png")

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)))

fname_save = "scratch_NM/NewWrapper/3-29/GTTargWNCov/log_apd_iso"
Nx=64
im = readsfd(Nx)
true_img = im[:, :, 1]

noise = rand(Normal(0.0, std(true_img)), Nx, Nx)
init = true_img + noise
#heatmap(init)
#heatmap(true_img)

filter_hash = fink_filter_hash(1, 8, nx=Nx, pc=1, wd=1, Omega=true)
(S1iso, Nf) = size(filter_hash["S1_iso_mat"])
(S2iso, Nf) = size(filter_hash["S2_iso_mat"])
coeff_mask = falses(2+S1iso+S2iso)
coeff_mask[S1iso+3:end] .= true
#coeff_mask = falses(2+Nf+Nf^2)
#coeff_mask[Nf+3:end] .= true

dhc_args = Dict(:doS2=>false, :doS12=>false, :doS20=>true, :apodize=>true, :iso=>true) #Iso #CHANGE: Change sig for the sfd data since the noise model is super high and the tiny values make sense
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", "ground_truth"), ("covar_type", "white_noise"), ("log", true), ("GaussianLoss", true), ("Invcov_matrix", "Diagonal+Eps"),
  ("optim_settings", optim_settings), ("lambda", std(apodizer(true_img)).^(-2)), ("white_noise_args", white_noise_args)]) #Add constraints
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=regs_true)
#Check
#s2mean - s_true

recon_settings["s_targ_mean"] = s2mean
recon_settings["s_invcov"] = s2icov
recon_settings["safemode"] = false
recon_settings["fname_save"] = fname_save * ".jld2"
recon_settings["optim_settings"] = optim_settings
resobj, recon_img = reconstruction_wrapper(true_img, init, filter_hash, dhc_args, coeff_mask, recon_settings)

p= plot([t.value for t in Optim.trace(resobj)])
plot!(title="Loss: S20 | Targ = S(GT) | Cov = SFD Covariance", xlabel = "No. Iterations", ylabel = "Loss")
savefig(p, fname_save * "_trace.png")

Visualization.plot_diffscales([apodizer(true_img), apodizer(init), apodizer(recon_img), apodizer(recon_img) - apodizer(true_img)], ["GT", "Init", "Reconstruction", "Residual"], fname=fname_save*".png")
Visualization.plot_synth_QA(apodizer(true_img), apodizer(init), apodizer(recon_img), fname=fname_save*"_6p.png")

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)))

#S20 + Iso check
#S20+NoIsoChecking that you didnt break old code
fname_save = "scratch_NM/NewWrapper/3-29/SFDTargSFDCov/log_apd_iso_exp2"
Nx=64
im = readsfd(Nx)
true_img = im[:, :, 1]

noise = rand(Normal(0.0, std(true_img)), Nx, Nx)
init = true_img + noise
#heatmap(init)
#heatmap(true_img)

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=>true) #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", true), ("GaussianLoss", true), ("Invcov_matrix", "Diagonal+Eps"),
  ("optim_settings", optim_settings), ("lambda", 0.00456), ("white_noise_args", white_noise_args)]) #Add constraints
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)
#Check
#s2mean - s_true

recon_settings["s_targ_mean"] = s2mean
recon_settings["s_invcov"] = s2icov
recon_settings["safemode"] = false
recon_settings["fname_save"] = fname_save * ".jld2"
recon_settings["optim_settings"] = optim_settings
resobj, recon_img = reconstruction_wrapper(true_img, init, filter_hash, dhc_args, coeff_mask, recon_settings)

p= plot([t.value for t in Optim.trace(resobj)])
plot!(title="Loss: S20 | Targ = S(SFD) | Cov = SFD Covariance", xlabel = "No. Iterations", ylabel = "Loss")
savefig(p, fname_save * "_trace.png")

Visualization.plot_diffscales([apodizer(true_img), apodizer(init), apodizer(recon_img), apodizer(recon_img) - apodizer(true_img)], ["GT", "Init", "Reconstruction", "Residual"], fname=fname_save*".png")
Visualization.plot_synth_QA(apodizer(true_img), apodizer(init), apodizer(recon_img), fname=fname_save*"_6p.png")

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)))

##fname_save = "scratch_NM/NewWrapper/3-29/SFDTargSFDCov/log_apd_iso_exp2"
Nx=64
im = readsfd(Nx)
true_img = im[:, :, 1]

noise = rand(Normal(0.0, std(true_img)), Nx, Nx)
init = true_img + noise
#heatmap(init)
#heatmap(true_img)

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=>true) #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", true), ("GaussianLoss", true), ("Invcov_matrix", "Diagonal+Eps"),
  ("optim_settings", optim_settings), ("lambda", 0.00456), ("white_noise_args", white_noise_args)]) #Add constraints
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)
#Check
#s2mean - s_true

recon_settings["s_targ_mean"] = s2mean
recon_settings["s_invcov"] = s2icov
recon_settings["safemode"] = false
recon_settings["fname_save"] = fname_save * ".jld2"
recon_settings["optim_settings"] = optim_settings
resobj, recon_img = reconstruction_wrapper(true_img, init, filter_hash, dhc_args, coeff_mask, recon_settings)

p= plot([t.value for t in Optim.trace(resobj)])
plot!(title="Loss: S20 | Targ = S(SFD) | Cov = SFD Covariance", xlabel = "No. Iterations", ylabel = "Loss")
savefig(p, fname_save * "_trace.png")

Visualization.plot_diffscales([apodizer(true_img), apodizer(init), apodizer(recon_img), apodizer(recon_img) - apodizer(true_img)], ["GT", "Init", "Reconstruction", "Residual"], fname=fname_save*".png")
Visualization.plot_synth_QA(apodizer(true_img), apodizer(init), apodizer(recon_img), fname=fname_save*"_6p.png")

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)))

max_stars_repo_name: andrew-saydjari/DHC

max_issues_repo_name: andrew-saydjari/DHC

max_forks_repo_name: andrew-saydjari/DHC


text:
#: ====================== concat  =====================

concat0(a::IntRange) = a
concat0(a::IntRange,b::IntRange) = first(a):last(b)

##: ===================== helpers =====================

@inline function correct_R!(C::Vector{Block{TR}})  where TR
    for i=1:length(C)
        C[i].R = compute_label(C[i])
    end
    return
end

##: ============== elementary operation ================

function merge_conseq_iden_blocks(
    C::Vector{Block{TR}}
    )::Vector{Block{TR}}  where TR

    if length(C)==0   return []   end
    @assert is_valid_C(C)
    N   = length(C)
    C1  = Block{TR}[]
    T   = copy(C[1])
    pos = 2
    while pos<=N
        Cp = C[pos]
        h  = label(Cp)   #! was hash(Cp)
        if h != label(T)
            T.R = compute_label(T)
            push!(C1, T)
            T = copy(Cp)
            pos += 1
            continue
        else # if h == T.R
            #: merge / accumulate
            T.n += Cp.n
            #@assert last(T.x)+1==first(Cp.x)
            T.x = concat0(T.x, Cp.x)
            pos += 1
            continue
        end
    end
    T.R = compute_label(T)
    push!(C1,T)

    correct_R!(C1)
    @assert is_valid_C(C1)

    return C1
end

function fold_C_by_blks(
    C::Vector{Block{TR}},
    blocks::Vector{IntRange}
    )::Vector{Block{TR}}  where TR

    if length(blocks)==0 || length(C)==0  return C  end
    @assert is_valid_C(C)

    B  = Set(vcat(collect.(blocks)...))
    @assert length(B)==sum(length.(blocks)) "blocks=$(blocks) overlapping !!!" #!! FIXME
    NC = length(C)
    C1 = Block{TR}[]

    i  = 1
    while i<=NC
        if i ∈ B
            ib = findfirst(x->i∈x, blocks)
            @assert i==first(blocks[ib])
            ML = Block(C[blocks[ib]])   #TODO ???
            push!(C1, ML)
            i = last(blocks[ib])+1
            for j ∈ blocks[ib]  delete!(B,j);  end
        else
            push!(C1, C[i])
            i += 1
        end
    end

    correct_R!(C1)
    return merge_conseq_iden_blocks(C1)
    return C1
end

# correct the .x and .n field of M1 according to M
function correct_x_n!(M::Block{TR}, M1::Block{TR})::Block{TR} where TR
    @assert first(M1.x) == first(M.x)
    @assert is_valid_x(M)
    @assert all(is_valid_x.(children(M1)))
    @assert length(M1.x)*M.n == length(M.x)
    M1.n = M1.n * M.n
    M1.x = M.x
    M1.R = compute_label(M1)
    return M1
end

function fold_block(b::Block{TR}, blocks::Vector{IntRange})::Block{TR} where TR
    return (length(blocks)>0
                ? correct_1_1_1(correct_x_n!(b, Block(fold_C_by_blks(children(b), blocks))))
                : b)
end

is_1_1(b) = (b.n==1 && length(children(b))==1)

function correct_1_1_1(b::Block{TR})::Block{TR} where TR
    for i=1:length(children(b))
        if is_1_1(b.C[i])
            b.C[i] = b.C[i].C[1]
        else
            b.C[i] = correct_1_1_1(b.C[i])
        end
    end
    return b
end

function merge_children(b::Block)::Block
    b1 = copy(b)
    b1.C = merge_conseq_iden_blocks(children(b))
    return b1
end

#: ----------- find blocks -----------

find_block(x::Block; block_identifier=MostFreqSimilarSubsq)::Block = (is_single(x) ? x : fold_block(x, block_identifier(label.(children(x)))))

# not used
#find_block_MostFreqSubsq(x::Block) = find_block(x; block_identifier=MostFreqSubsq)

find_block_MostFreqSimilarSubsq(x::Block) = find_block(x; block_identifier=MostFreqSimilarSubsq)

#: ----------- init blocks -----------

function build_block_init(patts::Vector{TPattern})::Block{RTYPE}
    return Block(merge_conseq_iden_blocks([Block(p,i) for (i,p) ∈ enumerate(patts)]))
end

# assuming that each "logical block" is terminated by an empty line
function build_block_init_by_linebreak(patts::Vector{TPattern})::Block{RTYPE}
    Q = empty_TPattern()   #  previous pattern
    S = Stack{Tuple{Int,Block{RTYPE}}}()
    L = 0   #  level
    for (i,p) ∈ enumerate(patts)
        if Q==Int[]  # previous line is empty
            L += 1   # increase level; next line is new block
        end
        #: enstack (level, pattern_i)
        push!(S, (L, Block(p,i)))
        if p==Int[]  # current line is empty
            if length(S)>0
                #: destack
                TMP  = Block{RTYPE}[]
                while length(S)>0 && top(S)[1]==L
                    t = pop!(S)
                    push!(TMP, t[2])
                end
                if length(S)>0
                    L = top(S)[1]
                else
                    L = 0
                end
                #: make Multiline (a tree) and enstack
                push!(S, (L, Block(merge_conseq_iden_blocks(TMP[end:-1:1]))))
            else  # stack is empty
                L = 0
            end
        end
        Q = p
        continue
    end

    # collect elements from stack bottom=1 --> top=n
    A = [last(s) for s in S][end:-1:1]
    return Block(merge_conseq_iden_blocks(A))
end

function typical_blocks(
    t::T;
    M = 2
    ) where {T <: Block}
    nl = []
    f(x) = (is_multi(x) && x.n>M ? (push!(nl,x); 0) : 0)
    DFS(t, f, identity)
    return sort(unique(nl),by=x->(-getfield(x,:n)))
end

max_stars_repo_name: algorithmx/SmartParser

max_issues_repo_name: algorithmx/SmartParser

max_forks_repo_name: algorithmx/SmartParser


text:
module Stat

import Gadfly
using DataFrames
using Compose
using Color

import Gadfly.Scale, Gadfly.Coord, Gadfly.element_aesthetics,
       Gadfly.default_scales, Gadfly.isconcrete, Gadfly.nonzero_length
import Distributions.Uniform, Distributions.kde
import Iterators.chain, Iterators.cycle, Iterators.product, Iterators.partition

include("bincount.jl")

# Apply a series of statistics.
#
# Args:
#   stats: Statistics to apply in order.
#   scales: Scales used by the plot.
#   aes: A Aesthetics instance.
#
# Returns:
#   Nothing, modifies aes.
#
function apply_statistics(stats::Vector{Gadfly.StatisticElement},
                          scales::Dict{Symbol, Gadfly.ScaleElement},
                          coord::Gadfly.CoordinateElement,
                          aes::Gadfly.Aesthetics)
    for stat in stats
        apply_statistic(stat, scales, coord, aes)
    end
    nothing
end

type Nil <: Gadfly.StatisticElement
end

const nil = Nil()

type Identity <: Gadfly.StatisticElement
end

function apply_statistic(stat::Identity,
                         scales::Dict{Symbol, Gadfly.ScaleElement},
                         coord::Gadfly.CoordinateElement,
                         aes::Gadfly.Aesthetics)
    nothing
end

const identity = Identity()

type HistogramStatistic <: Gadfly.StatisticElement
end

element_aesthetics(::HistogramStatistic) = [:x]

const histogram = HistogramStatistic()

function apply_statistic(stat::HistogramStatistic,
                         scales::Dict{Symbol, Gadfly.ScaleElement},
                         coord::Gadfly.CoordinateElement,
                         aes::Gadfly.Aesthetics)
    Gadfly.assert_aesthetics_defined("HistogramStatistic", aes, :x)
    d, bincounts = choose_bin_count_1d(aes.x)

    x_min, x_max = min(aes.x), max(aes.x)
    binwidth = (x_max - x_min) / d

    if aes.color === nothing
        aes.x_min = Array(Float64, d)
        aes.x_max = Array(Float64, d)
        aes.y = Array(Float64, d)

        for j in 1:d
            aes.x_min[j] = x_min + (j - 1) * binwidth
            aes.x_max[j] = x_min + j * binwidth
            aes.y[j] = bincounts[j]
        end
    else
        groups = Dict()
        for (x, c) in zip(aes.x, cycle(aes.color))
            if !haskey(groups, c)
                groups[c] = Float64[x]
            else
                push!(groups[c], x)
            end
        end

        aes.x_min = Array(Float64, d * length(groups))
        aes.x_max = Array(Float64, d * length(groups))
        aes.y = Array(Float64, d * length(groups))
        colors = Array(ColorValue, d * length(groups))

        x_min = min(aes.x)
        x_max = max(aes.x)
        stack_height = zeros(Int, d)
        for (i, (c, xs)) in enumerate(groups)
            fill!(bincounts, 0)
            for x in xs
                bin = max(1, min(d, int(ceil((x - x_min) / binwidth))))
                bincounts[bin] += 1
            end
            stack_height += bincounts[1:d]

            for j in 1:d
                idx = (i-1)*d + j
                aes.x_min[idx] = x_min + (j - 1) * binwidth
                aes.x_max[idx] = x_min + j * binwidth
                aes.y[idx] = bincounts[j]
                colors[idx] = c
            end
        end

        y_drawmax = float64(max(stack_height))
        if aes.y_drawmax === nothing || aes.y_drawmax < y_drawmax
            aes.y_drawmax = y_drawmax
        end

        aes.color = PooledDataArray(colors)
    end
end

type DensityStatistic <: Gadfly.StatisticElement
    # Number of points sampled
    n::Int
end

const density = DensityStatistic(300)

element_aesthetics(::DensityStatistic) = [:x, :y]

function apply_statistic(stat::DensityStatistic,
                         scales::Dict{Symbol, Gadfly.ScaleElement},
                         coord::Gadfly.CoordinateElement,
                         aes::Gadfly.Aesthetics)
    Gadfly.assert_aesthetics_defined("DensityStatistic", aes, :x)

    # TODO: handle grouping by color

    f = kde(aes.x, stat.n)
    aes.x = f.x
    aes.y = f.density
end

type RectangularBinStatistic <: Gadfly.StatisticElement
end

element_aesthetics(::RectangularBinStatistic) = [:x, :y, :color]

default_scales(::RectangularBinStatistic) = [Gadfly.Scale.color_gradient]

const rectbin = RectangularBinStatistic()

function apply_statistic(stat::RectangularBinStatistic,
                         scales::Dict{Symbol, Gadfly.ScaleElement},
                         coord::Gadfly.CoordinateElement,
                         aes::Gadfly.Aesthetics)

    dx, dy, bincounts = choose_bin_count_2d(aes.x, aes.y)

    x_min, x_max = min(aes.x), max(aes.x)
    y_min, y_max = min(aes.y), max(aes.y)

    # bin widths
    wx = (x_max - x_min) / dx
    wy = (y_max - y_min) / dy

    n = 0
    for cnt in bincounts
        if cnt > 0
            n += 1
        end
    end

    aes.x_min = Array(Float64, n)
    aes.x_max = Array(Float64, n)
    aes.y_min = Array(Float64, n)
    aes.y_max = Array(Float64, n)

    k = 1
    for ((i, j), cnt) in zip(product(1:dy, 1:dx), bincounts)
        if cnt > 0
            aes.x_min[k] = x_min + (i - 1) * wx
            aes.x_max[k] = x_min + i * wx
            aes.y_min[k] = y_min + (j - 1) * wy
            aes.y_max[k] = y_min + j * wy
            k += 1
        end
    end

    if !has(scales, :color)
        error("RectangularBinStatistic requires a color scale.")
    end
    color_scale = scales[:color]
    if !(typeof(color_scale) <: Scale.ContinuousColorScale)
        error("RectangularBinStatistic requires a continuous color scale.")
    end

    aes.color_key_title = "Count"

    data = Gadfly.Data()
    data.color = Array(Int, n)
    k = 1
    for cnt in bincounts
        if cnt > 0
            data.color[k] = cnt
            k += 1
        end
    end
    Scale.apply_scale(color_scale, [aes], data)
    nothing
end

default_statistic(stat::RectangularBinStatistic) = [Scale.color_gradient]

# Find reasonable places to put tick marks and grid lines.
type TickStatistic <: Gadfly.StatisticElement
    in_vars::Vector{Symbol}
    out_var::String
end

const x_ticks = TickStatistic([:x, :x_min, :x_max, :x_drawmin, :x_drawmax], "x")
const y_ticks = TickStatistic(
    [:y, :y_min, :y_max, :middle, :lower_hinge, :upper_hinge,
     :lower_fence, :upper_fence, :y_drawmin, :y_drawmax], "y")

# Can a numerical value be treated as an integer
is_int_compatable(::Integer) = true
is_int_compatable{T <: FloatingPoint}(x::T) = abs(x) < maxintfloat(T) && float(int(x)) == x
is_int_compatable(::Any) = false

# Apply a tick statistic.
#
# Args:
#   stat: statistic.
#   aes: aesthetics.
#
# Returns:
#   nothing
#
# Modifies:
#   aes
#
function apply_statistic(stat::TickStatistic,
                         scales::Dict{Symbol, Gadfly.ScaleElement},
                         coord::Gadfly.CoordinateElement,
                         aes::Gadfly.Aesthetics)
    in_values = [getfield(aes, var) for var in stat.in_vars]
    in_values = filter(val -> !(val === nothing), in_values)
    in_values = chain(in_values...)
    # TODO: handle the outliers aesthetic

    minval = Inf
    maxval = -Inf
    all_int = true

    n = 0
    for val in in_values
        if val < minval
            minval = float64(val)
        end

        if val > maxval
            maxval = float64(val)
        end

        if !is_int_compatable(val)
            all_int = false
        end

        n += 1
    end

    # Take into account a forced viewport in cartesian coordinates.
    if typeof(coord) == Coord.Cartesian
        if stat.out_var == "x"
            if !is(coord.xmin, nothing)
                minval == min(minval, float64(coord.xmin))
            end
            if !is(coord.xmax, nothing)
                maxval == max(maxval, float64(coord.xmax))
            end
        elseif stat.out_var == "y"
            if !is(coord.ymin, nothing)
                minval == min(minval, float64(coord.ymin))
            end
            if !is(coord.ymax, nothing)
                maxval == min(maxval, float64(coord.ymax))
            end
        end
    end

    # all the input values in order.
    if all_int
        ticks = Set{Float64}()
        union!(ticks, in_values)
        ticks = Float64[t for t in ticks]
        sort!(ticks)

        maxgap = 0
        for (i, j) in partition(ticks, 2, 1)
            if j - i > maxgap
                maxgap = j -i
            end
        end

        if length(ticks) > 20 || maxgap > 1
            grids = ticks = Gadfly.optimize_ticks(minval, maxval)
            if ticks[1] == 0
                ticks[1] = 1
            end
            grids = ticks
        else
            grids = (ticks - 0.5)[2:end]
        end
        viewmin = min(ticks)
        viewmax = max(ticks)
    else
        ticks = Gadfly.optimize_ticks(minval, maxval)
        viewmin = min(ticks)
        viewmax = max(ticks)

        # Extend ticks
        d = ticks[2] - ticks[1]
        lowerticks = ticks - (ticks[end] - ticks[1]) - d
        upperticks = ticks + (ticks[end] - ticks[1]) + d
        grids = ticks = vcat(lowerticks, ticks, upperticks)
    end

    # We use the first label function we find for any of the aesthetics. I'm not
    # positive this is the right thing to do, or would would be.
    labeler = getfield(aes, symbol(string(stat.out_var, "_label")))

    setfield(aes, symbol(string(stat.out_var, "tick")), ticks)
    setfield(aes, symbol(string(stat.out_var, "grid")), grids)
    setfield(aes, symbol(string(stat.out_var, "tick_label")), labeler)

    viewmin_var = symbol(string(stat.out_var, "_viewmin"))
    if getfield(aes, viewmin_var) === nothing ||
       getfield(aes, viewmin_var) > viewmin
        setfield(aes, viewmin_var, viewmin)
    end

    viewmax_var = symbol(string(stat.out_var, "_viewmax"))
    if getfield(aes, viewmax_var) === nothing ||
       getfield(aes, viewmax_var) < viewmax
        setfield(aes, viewmax_var, viewmax)
    end

    nothing
end

type BoxplotStatistic <: Gadfly.StatisticElement
end

element_aesthetics(::BoxplotStatistic) = [:x, :y]

const boxplot = BoxplotStatistic()

function apply_statistic(stat::BoxplotStatistic,
                         scales::Dict{Symbol, Gadfly.ScaleElement},
                         coord::Gadfly.CoordinateElement,
                         aes::Gadfly.Aesthetics)
    Gadfly.assert_aesthetics_defined("BoxplotStatistic", aes, :y)

    groups = Dict()

    aes_x = aes.x === nothing ? [nothing] : aes.x
    aes_color = aes.color === nothing ? [nothing] : aes.color

    for (x, y, c) in zip(cycle(aes_x), aes.y, cycle(aes_color))
        if !haskey(groups, (x, c))
            groups[(x, c)] = Float64[y]
        else
            push!(groups[(x, c)], y)
        end
    end

    m = length(groups)
    aes.middle = Array(Float64, m)
    aes.lower_hinge = Array(Float64, m)
    aes.upper_hinge = Array(Float64, m)
    aes.lower_fence = Array(Float64, m)
    aes.upper_fence = Array(Float64, m)
    aes.outliers = Vector{Float64}[]

    for (i, ((x, c), ys)) in enumerate(groups)
        aes.lower_hinge[i], aes.middle[i], aes.upper_hinge[i] =
                quantile(ys, [0.25, 0.5, 0.75])
        iqr = aes.upper_hinge[i] - aes.lower_hinge[i]
        aes.lower_fence[i] = aes.lower_hinge[i] - 1.5iqr
        aes.upper_fence[i] = aes.upper_hinge[i] + 1.5iqr
        push!(aes.outliers,
             filter(y -> y < aes.lower_fence[i] || y > aes.upper_fence[i], ys))
    end

    if !is(aes.x, nothing)
        aes.x = Int64[x for (x, c) in keys(groups)]
    end

    if !is(aes.color, nothing)
        aes.color = PooledDataArray(ColorValue[c for (x, c) in keys(groups)],
                                    levels(aes.color))
    end

    nothing
end

end # module Stat


text:
################################################################################
#
#  Conversion to parseble Hecke code
#
################################################################################

# This is helpful to construct code for the tests.
function to_hecke(L::AbsLat; target = "L", skip_field = false)
  return to_hecke(stdout, L, target = target, skip_field = skip_field)
end

function to_hecke_string(L::AbsLat; target = "L", skip_field = false)
  b = IOBuffer()
  to_hecke(b, L, target = target, skip_field = skip_field)
  return String(take!(b))
end

function to_hecke(io::IO, L::QuadLat; target = "L", skip_field = false)
  K = nf(base_ring(L))
  if !skip_field
    println(io, "Qx, x = PolynomialRing(FlintQQ, \"x\", cached = false)")
    f = defining_polynomial(K)
    pol = string(f)
    pol = replace(pol, string(var(parent(f))) => "x")
    println(io, "f = ", pol, ";")
    println(io, "K, a = number_field(f)")
  end
  F = gram_matrix(ambient_space(L))
  Fst = "[" * split(string([F[i, j] for i in 1:nrows(F) for j in 1:ncols(F)]), '[')[2]
  Fst = replace(Fst, string(var(K)) => "a")
  println(io, "D = matrix(K, ", nrows(F), ", ", ncols(F), ", ", Fst, ");")
  gens = generators(L)
  Gs = "["
  for i in 1:length(gens)
    g = gens[i]
    Gs = Gs * "[" * split(string(g), "[")[2]
    if i < length(gens)
      Gs = Gs * ", "
    end
  end
  Gs = Gs * "]"
  println(io, "gens = ", Gs)
  println(io, target, " = quadratic_lattice(K, generators = gens, gram_ambient_space = D)")
end

function to_hecke(io::IO, L::HermLat; target = "L", skip_field = skip_field)
  E = nf(base_ring(L))
  K = base_field(E)
  println(io, "Qx, x = PolynomialRing(FlintQQ, \"x\")")
  f = defining_polynomial(K)
  pol = replace(string(f), "//" => "/")
  pol = replace(pol, string(var(parent(f))) => "x")
  println(io, "f = ", pol)
  println(io, "K, a = NumberField(f, \"a\", cached = false)")
  println(io, "Kt, t = PolynomialRing(K, \"t\")")
  f = defining_polynomial(E)
  pol = replace(string(f), "//" => "/")
  pol = replace(pol, string(var(parent(f))) => "t")
  println(io, "g = ", pol, "")
  println(io, "E, b = NumberField(g, \"b\", cached = false)")
  F = gram_matrix(ambient_space(L))
  Fst = "[" * split(string([F[i, j] for i in 1:nrows(F) for j in 1:ncols(F)]), '[')[2]
  Fst = replace(Fst, string(var(K)) => "a")
  Fst = replace(Fst, string(var(E)) => "b")
  println(io, "D = matrix(E, ", nrows(F), ", ", ncols(F), ", ", Fst, ")")

  gens = generators(L)
  Gs = "Vector{$(elem_type(E))}["
  for i in 1:length(gens)
    g = gens[i]
    gst = replace(string(g), string(var(K)) => "a")
    gst = replace(string(g), string(var(E)) => "b")

    Gs = Gs * "map(E, [" * split(gst, "[")[2] * ")"
    if i < length(gens)
      Gs = Gs * ", "
    end
  end
  Gs = Gs * "]"
  println(io, "gens = ", Gs)

  println(io, target, " = hermitian_lattice(E, generators = gens, gram_ambient_space = D)")
end

################################################################################
#
#  Conversion to Magma
#
################################################################################

function to_magma(L::AbsLat; target = "L")
  return to_magma(stdout, L, target = target)
end

function to_magma_string(L::AbsLat; target = "L")
  b = IOBuffer()
  to_magma(b, L, target = target)
  return String(take!(b))
end

function to_magma(io::IO, L::HermLat; target = "L")
  E = nf(base_ring(L))
  K = base_field(E)
  println(io, "Qx<x> := PolynomialRing(Rationals());")
  f = defining_polynomial(K)
  pol = replace(string(f), "//" => "/")
  pol = replace(pol, string(var(parent(f))) => "x")
  println(io, "f := ", pol, ";")
  println(io, "K<a> := NumberField(f : DoLinearExtension);")
  println(io, "Kt<t> := PolynomialRing(K);")
  f = defining_polynomial(E)
  pol = replace(string(f), "//" => "/")
  pol = replace(pol, string(var(parent(f))) => "t")
  println(io, "g := ", pol, ";")
  println(io, "E<b> := NumberField(g : DoLinearExtension);")
  F = gram_matrix(ambient_space(L))
  Fst = "[" * split(string([F[i, j] for i in 1:nrows(F) for j in 1:ncols(F)]), '[')[2]
  println(io, "F := Matrix(E, ", nrows(F), ", ", ncols(F), ", ", Fst, ");")
  pm = pseudo_matrix(L)
  M = matrix(pm)
  Mst = "[" * split(string([M[i, j] for i in 1:nrows(M) for j in 1:ncols(M)]), '[')[2]
  Mst = replace(Mst, "//" => "/")
  println(io, "M := Matrix(E, ", nrows(M), ", ", ncols(M), ", ", Mst, ");")
  println(io, "OE := MaximalOrder(E);")
  print(io, "C := [ ")
  for (i, I) in enumerate(coefficient_ideals(pm))
    print(io, "ideal< OE | ")
    bas = "[" * split(string(absolute_basis(I)), '[')[2]
    bas = replace(bas, string(var(K)) => "a")
    bas = replace(bas, string(var(E)) => "b")
    bas = replace(bas, "//" => "/")
    if i < length(coefficient_ideals(pm))
      print(io, bas, ">, ")
    else
      println(io, bas, ">];")
    end
  end
  println(io, "M := Module(PseudoMatrix(C, M));")
  println(io, "$target := HermitianLattice(M, F);")
end

function to_magma(io::IO, L::AbsLat{AnticNumberField}; target = "L")
  K = nf(base_ring(L))
  println(io, "Qx<x> := PolynomialRing(Rationals());")
  f = defining_polynomial(K)
  pol = replace(string(f), "//" => "/")
  pol = replace(pol, string(var(parent(f))) => "x")
  println(io, "f := ", pol, ";")
  println(io, "K<a> := NumberField(f : DoLinearExtension);")
  F = gram_matrix(ambient_space(L))
  Fst = "[" * split(string([F[i, j] for i in 1:nrows(F) for j in 1:ncols(F)]), '[')[2]
  Fst = replace(Fst, "//" => "/")
  Fst = replace(Fst, string(var(K)) => "a")
  println(io, "F := Matrix(K, ", nrows(F), ", ", ncols(F), ", ", Fst, ");")
  pm = pseudo_matrix(L)
  M = matrix(pm)
  Mst = "[" * split(string([M[i, j] for i in 1:nrows(M) for j in 1:ncols(M)]), '[')[2]
  Mst = replace(Mst, "//" => "/")
  println(io, "M := Matrix(K, ", nrows(M), ", ", ncols(M), ", ", Mst, ");")
  println(io, "OK := MaximalOrder(K);")
  print(io, "C := [ ")
  for (i, I) in enumerate(coefficient_ideals(pm))
    print(io, "ideal< OK | ")
    bas = "[" * split(string(absolute_basis(I)), '[')[2]
    bas = replace(bas, string(var(K)) => "a")
    bas = replace(bas, "//" => "/")
    if i < length(coefficient_ideals(pm))
      print(io, bas, ">, ")
    else
      println(io, bas, ">];")
    end
  end
  println(io, "M := Module(PseudoMatrix(C, M));")
  if L isa HermLat
    println(io, "$target := HermitianLattice(M, F);")
  else
    println(io, "$target := LatticeModule(M, F);")
  end
end

function var(E::NfRel)
  return E.S
end

################################################################################
#
#  Zlattices
#
################################################################################

function to_hecke(io::IO, L::ZLat; target = "L")
  B = basis_matrix(L)
  G = gram_matrix(ambient_space(L))
  Bst = "[" * split(string([B[i, j] for i in 1:nrows(B) for j in 1:ncols(B)]), '[')[2]
  Gst = "[" * split(string([G[i, j] for i in 1:nrows(G) for j in 1:ncols(G)]), '[')[2]
  println(io, "B = matrix(FlintQQ, ", nrows(B), ", ", ncols(B), " ,", Bst, ");")
  println(io, "G = matrix(FlintQQ, ", nrows(G), ", ", ncols(G), " ,", Gst, ");")
  println(io, target, " = ", "Zlattice(B, gram = G);")
end


text:
# 2D Heat Diffusion
# LBM
#
# Hongyang Zhou, hyzhou@umich.edu 05/01/2019

# Problem Description
# A two dimensional square slab is subjected to the given boundary conditions.
# Initially the slab was at zero temperature. For time>0, the boundary at x = 0
# is subjected to a high temperature of value 1.0 and other boundaries are kept
# as before. The length of the domain is 100 units. Determine the temperature
# distribution in the slab at time = 400 units (s).
# Thermal diffusivity is 0.25 and xl = yl = 1.0.

using PyPlot

function solveHeatDiffusion2D()
# Parameters
nD = 2             # Dimension of the system
Δx = 1.0           # Discretized space interval
Δt = 1.0           # Discretized time interval
α  = 0.25          # Thermal diffusivity
ω  = inv(nD*Δt*α/Δx^2 + 0.5)  # Relaxation coefficient
TotalTime = 400.   # Total simulation time
nT = Int(TotalTime/Δt)        # Total timesteps
nx = 100           # Number of total space intervals
ny = 100
w  = [0.25,0.25,0.25,0.25]    # Lattice weights
nw = size(w,1)

Tleft  = 1.0       # Left boundary temperature
Tright = 0.0       # Right boundary temperature
Ttop   = 0.0       # Top boundary temperature
Tbot   = 0.0       # Bottom boundary temperature

heatflux  = 100.   # Constant heat flux, [W/mK]
k_thermal = 20.    # Thermal conductivity, [W/mK]

BCDirichlet,BCRiemann = 1,2 # Constants for different boundary types
BCtype = 1         # Type of boundary conditions

UseSource = false  # Logical for calculating source terms
if UseSource
   qg = 1.0        # Rate of heat generation per volume
   ρ  = 1.0        # Density
   C  = 1.0        # Heat capacity
   source = qg/(ρ*C) # Source term. Note that α=k/(ρC)
end

# Initializations
f   = zeros(nw,nx,ny)
feq = zeros(nw,nx,ny) # This can be moved inside the main loop!
T   = zeros(1,nx,ny)
t   = 0.0     # Current simulation time

# Advance
while t < TotalTime

   # Calculate equilibrium distribution
   for iw = 1:nw
      feq[iw,:,:] .= w[iw].*T[1,:,:]
   end

   # Collision
   for iw = 1:nw
      f[iw,:,:] .= (1.0-ω).*f[iw,:,:] .+ ω.*feq[iw,:,:]
   end
   if UseSource
      for iw = 1:nw
         f[iw,:,:] .+= Δt*w[iw]*source
      end
   end

   # Streaming
   f[1,2:nx,:]   .= f[1,1:nx-1,:]
   f[2,1:nx-1,:] .= f[2,2:nx,:]
   f[3,:,2:ny]   .= f[3,:,1:ny-1]
   f[4,:,1:ny-1] .= f[4,:,2:ny]

   # Remaining BCs
   if BCtype == BCDirichlet
      f[1,1,:] .= (w[1]+w[2])*Tleft  .- f[2,1,:]
      f[2,nx,:].= (w[1]+w[2])*Tright .- f[1,nx,:]
      f[3,:,1] .= (w[3]+w[4])*Ttop   .- f[4,:,1]
      f[4,:,ny].= (w[3]+w[4])*Tbot   .- f[3,:,ny]
   elseif BCtype == BCRiemann
      # This is not complete, and not flexible enough!
      f[1,1,:] .= f[1,2,:] .- heatflux*Δx/k_thermal
      #f[2,nx,:] = f[2,nx-1,:] .- heatflux*Δx/k_thermal
   end

   # Calculate macroscopic T
   T .= sum(f, dims=1)

   t += Δt
end

# Visualization
x = range(0,stop=nx-1,length=nx)
y = range(0,stop=ny-1,length=ny)

contourf(x,y,T[1,:,:]); colorbar()
title("2D Heat Diffusion in a Plate")
xlabel("x")
ylabel("y")
grid("on")

# 1D cut plot
#plot(0:nx-1,T[1,:,50],linestyle="-",marker="o",label="LBM")
#xlim(0,30)

end


text:
const max_level = 42

const score_array = [
  100,
  300,
  500,
  800
]

const allowed_labels = [
  "clear",
  "combo",
  "soft",
  "hard",
  "glue",
  "hold"
]

const color_list = [
  "blue",
  "cyan",
  "green",
  "orange",
  "purple",
  "red",
  "yellow"
]

const crayon_dict = Dict()

crayon_dict["invisible"] = Crayon()

for cur_color in color_list
  tmp_color = cur_color

  ( tmp_color == "blue" ) &&
    ( tmp_color = "royalblue")

  cur_piece_color = "$(tmp_color)3"
  cur_shadow_color = "$(tmp_color)1"

  if tmp_color == "yellow"
    cur_piece_color = "yellow2"
  end

  crayon_dict[cur_color] = Crayon(
    background = Colors.color_names[cur_piece_color]
  )

  crayon_dict["$(cur_color)_shadow"] = Crayon(
    background = Colors.color_names[cur_shadow_color]
  )
end

crayon_dict["grey"] = Crayon(
  foreground = Colors.color_names["darkgrey"]
)

const odd_offsets = Dict(
  [ 0 , 1 ] => (
    ( -1 , +0 ),
    ( -1 , +1 ),
    ( +0 , -2 ),
    ( -1 , -2 )
  ),
  [ 1 , 2 ] => (
    ( +1 , +0 ),
    ( +1 , -1 ),
    ( +0 , +2 ),
    ( +1 , +2 )
  ),
  [ 2 , 3 ] => (
    ( +1 , +0 ),
    ( +1 , +1 ),
    ( +0 , -2 ),
    ( +1 , -2 )
  ),
  [ 0 , 3 ] => (
    ( -1 , +0 ),
    ( -1 , -1 ),
    ( +0 , +2 ),
    ( -1 , +2 )
  )
)

const even_offsets = Dict(
  [ 0 , 1 ] => (
    ( -2 , +0 ),
    ( +1 , +0 ),
    ( -2 , -1 ),
    ( +1 , +2 )
  ),
  [ 1 , 2 ] => (
    ( -1 , +0 ),
    ( +2 , +0 ),
    ( -1 , +2 ),
    ( +2 , -1 )
  ),
  [ 2 , 3 ] => (
    ( +2 , +0 ),
    ( -1 , +0 ),
    ( +2 , +1 ),
    ( -1 , -2 )
  ),
  [ 0 , 3 ] => (
    ( +1 , +0 ),
    ( -2 , +0 ),
    ( +1 , -2 ),
    ( -2 , +1 )
  )
)


text:
function calc_solution(xs)
    s = Set{Int64}()
    for x in xs
        if x in s
            return x * (2020 - x)
        end
        push!(s, -x + 2020)
    end
end

function main()
    lines = readlines("../resources/input_01.txt")
    xs = map(line -> parse(Int, line), lines)
    result = calc_solution(xs)
    println(result)
end

main()


text:
"""
    plotseis(S[, fmt=FMT, use_name=false, n=N])

Renormalized, time-aligned trace plot of data in S.x using timestamps in S.t.

Keywords:
* fmt=FMT formats x-axis labels using C-language `strftime` format string `FMT`.
If unspecified, the format is determined by when data in `S` start and end.
* use_name=true uses `S.name`, rather than `S.id`, for labels.
* n=N sets the number of X-axis ticks.

"""
function plotseis(S::SeisData; fmt::String="auto", use_name::Bool=false, nxt::Int64=5)
  xmi = typemax(Int64)
  xma = typemin(Int64)

  fig = PyPlot.figure(figsize=[8.0, 6.0], dpi=150)
  ax = PyPlot.axes([0.20, 0.10, 0.72, 0.85])
  for i = 1:S.n
    x = rescaled(S.x[i].-mean(S.x[i]),i)
    if S.fs[i] > 0
      t = t_expand(S.t[i], S.fs[i])
      plot(t, x, linewidth=1.0)
    else
      t = view(S.t[i], :, 2)
      plot(t, x, "o", linewidth=1, markeredgecolor=[0,0,0])
    end
    xmi = min(xmi, t[1])
    xma = max(xma, t[end])
  end

  xfmt(xmi, xma, fmt, nxt)
  PyPlot.setp(gca().get_yticklabels(), fontsize=8.0, color="black", fontweight="bold", family="serif")
  PyPlot.setp(gca().get_xticklabels(), fontsize=8.0, color="black", fontweight="bold", family="serif")
  PyPlot.yticks(1:S.n, map((i) -> replace(i, " " => ""), use_name ? S.name : S.id))
  PyPlot.ylim(0.5, S.n+0.5)
  return fig
end


text:
module TestFirstOrder2

using TinyModia
using DifferentialEquations
using ModiaPlot

# using RuntimeGeneratedFunctions
# RuntimeGeneratedFunctions.init(@__MODULE__)

inputSignal(t) = sin(t)

FirstOrder1 = Model(
    T = 0.2,
    x = Var(init=0.3),
    equations = :[u = inputSignal(time/u"s"),
                  T * der(x) + x = u,
                  y = 2*x]
)

FirstOrder2 = FirstOrder1 | Map(T = 0.3, x = Var(init=0.6))

firstOrder = @instantiateModel(FirstOrder2, logCode=true)

simulate!(firstOrder, Tsit5(), stopTime = 10, merge = Map(T = 0.4, x = 0.9),
          log=true, logParameters=true, logStates=true,
          requiredFinalStates = [-0.17964872595554535])

plot(firstOrder, [("u", "x"), "der(x)", "y"])

end


text:
#gaussian_process.jl

"""
Returns the values of the mean and varience as well as the distribution (this is returned for ease of iteration)
after adding a single point [x,f]=D or a set of points D=[x1,f1, x2]. Note this is a very inefficient way of
computing the gaussian_process as it requires generating the variance covarience matrix for all possible
selection values. Instead use the cholesky version found below.
Themean and varience functions are as described in A Tutorial on Bayesian Optimization of Expensive Cost
Functions, with Application to Active User Modeling and  Hierarchical Reinforcement Learning Eric Brochu,
Vlad M. Cora and Nando de Freitas December 14, 2010.

Arguments
---------
Ker
    Kernal function
D
    Dataset as tuples of (x,y), where y is f(x)+noise
noise
    noise value (often) 1e-6
xrange
    range of x points to choose (x axis)

"""

function gaussian_process(ker,D,noise,xrange)
    K_ss=cov_gen(ker,xrange,xrange)+eye(length(xrange))*noise
    y=map(x->x[2],D)   # these are our y noisy functions
    x=map(x->x[1],D) #These are our x training points
    K = cov_gen(ker,x,x)+noise*eye(length(x))
    K_s=cov_gen(ker,x,xrange)  #K_s is as in nandos code but in the paper called little k
    Inv_K=inv(K)
    µ=transpose(K_s)*Inv_K*y
    sigma2=diag(K_ss-transpose(K_s)*Inv_K*K_s)
    sigma=sqrt(sigma2)
    return (µ,sigma,D)
end

#Another temporary one:

#gaussian_process_chol.jl

"""
Computes the gaussian process of a given set of points xrange with given noise
Arguments
---------
Ker
    Kernal function
D
    Dataset as tuples of (x,y), where y is f(x)+noise
noise
    noise value (often) 1e-6
xrange
    range of x points to choose (x axis)

"""

function gaussian_process_chol(ker,D,noise,xrange)

    y=map(x->x[2],D)   # these are our y noisy functions
    x=map(x->x[1],D) #These are our x training points
    K = cov_gen(ker,x,x)+noise*eye(length(x))
    L=ctranspose(chol(K))
    # temp=cov_gen(ker,x,xrange)
    Lk=\(L,cov_gen(ker,x,xrange))

    mu=[dot(Lk[:,i],(\(L,y))) for i=1:size(Lk)[2] ]
    K_diag=diag_cov_gen(ker,xrange)+ones(length(xrange))*noise
    s2=K_diag-[ sum( (Lk.*Lk)[:,i] ) for i=1:size(Lk)[2] ]
    sigma=sqrt(s2)
    return (mu,sigma,D)
end


text:
## Load data for testing
using Plots
using SparseArrays
using LinearAlgebra
using TextParse
using DelimitedFiles
using PyCall
using ColorSchemes

function read_graph(filename)
    src, dst, val = csvread(filename;
        spacedelim=true, colparsers=[Int,Int,Float64], header_exists=false)[1]
    src .+= 1
    dst .+= 1
    n = maximum(src)
    n = max(maximum(dst),n)
    A = sparse(src, dst, val, n, n)
    A = max.(A,A')
    return findnz(triu(A,1))[1:2]
end

src, dst = read_graph("../../dataset/lawlor-spectra-k32.edgelist.gz")
# x, y = csvread("$../../dataset/lawlor-spectra-k32.edgelist";
#         spacedelim=true, colparsers=[Float64,Float64],
#         header_exists=false)[1]
V = readdlm("../../dataset/lawlor-spectra-k32.coords")
V[:,1] *= -10
V[:,2] *= -4

function plotperm(v)
  return invperm(sortperm(v))
end
# function myscatter(x,y,color;kwargs...)
#   p = sortperm(color)
#   scatter(x[p],y[p], marker_z = color[p],
#     label="", alpha=0.1, markerstrokewidth=0,
#     colorbar=false, framestyle=:none; kwargs...)
# end

##
pfile = "embeddings/flow_3_bfs_delta_0.1.p.gz"

pickle = pyimport("pickle")
gzip = pyimport("gzip")
f = gzip.open(pfile)
data= pickle.load(f)
f.close()
function _pycsr_to_sparse(rowptr, colinds, vals, m, n)
  rowinds = zeros(eltype(colinds),0)
  for i in 1:length(rowptr)-1
    for j in rowptr[i]:rowptr[i+1]-1
      push!(rowinds, i)
    end
  end
  return sparse(rowinds, colinds, vals, m, n)
end
function pycsr_to_sparse(pycsr)
  rowptr = pycsr.indptr
  colinds = pycsr.indices
  vals = pycsr.data
  colinds .+= 1
  m, n = data.shape
  # put in a fucntion so we can do type inference
  return _pycsr_to_sparse(rowptr, colinds, vals, m, n )
end

X = pycsr_to_sparse(data)

flipexp = findall(X[517173,:] .> 0) # this is at the opposite end
for f in flipexp
  X[:,f] = 1.0 .-X[:,f]
end

X = Matrix(X) # convert to Matrix
X ./= 10 # convert to 0, 1
subset = vec((sum(X,dims=2) .> 0))
subsetnodes = findall(subset)
subset_c = vec((sum(X,dims=2) .<= 0))
subsetnodes_c = findall(subset_c)
U,sig,Vt = svd(X.-1)

new_coords = zeros(size(X,1),2)
new_coords[subsetnodes,1] = plotperm(-1 * U[subsetnodes,1])
new_coords[subsetnodes,2] = plotperm(-1 * U[subsetnodes,2])
i = argmax(new_coords[subsetnodes,1])
new_coords[subsetnodes_c,1] .= new_coords[subsetnodes[i],1] + 10000
new_coords[subsetnodes_c,2] .= new_coords[subsetnodes[i],2] + 10000

x,y = new_coords[:,1],new_coords[:,2]

## Form a random subset of edges
using Random
Random.seed!(0)
p = randperm(length(src))[1:10^5] # random set of 1M edges
srcsub = Int.(src[p])
dstsub = Int.(dst[p])

## Hexbin graph plots
include("hexbinplots_orig.jl")
function shapecoords(hx::HexBinPlots.HexHistogram)
  h,vh = HexBinPlots.make_shapes(hx)
  vmax = maximum(vh)
  xsh = Vector{Float64}()
  ysh = Vector{Float64}()
  for k in eachindex(h)
      append!(xsh,h[k].x)
      push!(xsh,h[k].x[1])
      push!(xsh,NaN)
      append!(ysh,h[k].y)
      push!(ysh,h[k].y[1])
      push!(ysh,NaN)
  end
  return (xsh, ysh, vh)
end
using Interpolations
function control_point(xi, xj, yi, yj, dist_from_mid)
    xmid = 0.5 * (xi+xj)
    ymid = 0.5 * (yi+yj)
    # get the angle of y relative to x
    theta = atan((yj-yi) / (xj-xi)) + 0.5pi
    # dist = sqrt((xj-xi)^2 + (yj-yi)^2)
    # dist_from_mid = curvature_scalar * 0.5dist
    # now we have polar coords, we can compute the position, adding to the midpoint
    (xmid + dist_from_mid * cos(theta),
     ymid + dist_from_mid * sin(theta))
end
function _add_edge_points!(xlines,ylines,xsi,ysi,xdi,ydi,curve_amount,interplen)
  xpt, ypt = control_point(xsi,xdi,ysi,ydi,curve_amount)
  xpts = [xsi, xpt, xdi]
  ypts = [ysi, ypt, ydi]
  t = range(0, stop=1, length=3)
  A = hcat(xpts, ypts)
  itp = scale(interpolate(A, BSpline(Cubic(Natural(OnGrid())))), t, 1:2)
  tfine = range(0, 1, length=interplen)
  for t in tfine
    push!(xlines, itp(t, 1))
  end
  for t in tfine
    push!(ylines, itp(t, 2))
  end
end

# Mutating version from HexBinPlots
function _xy2counts!(counts::Dict{(Tuple{Int, Int}), Int},
                    x::AbstractArray,y::AbstractArray,
                    xsize::Real,ysize::Real,x0,y0)
  @inbounds for i in eachindex(x)
      h = convert(HexagonOffsetOddR,
                  cube_round(x[i] - x0,y[i] - y0,xsize, ysize))
      idx = (h.q, h.r)
      counts[idx] = 1 + get(counts,idx,0)
  end
  counts
end
using Statistics
using Hexagons
# Generate counts for a graph edge plot histogram
# hhparam = (xsize,ysize,x0,y0) are the hex-histogrma params
function graphedgehist(src,dst,x,y,ptscale,hhparams,curvature = 0.05)
  xlines = similar(x,0)
  ylines = similar(y,0)
  sizehint!(xlines,max(length(src)+1000,10000000)) # we'll process them in batches this big
  sizehint!(ylines,max(length(src)+1000,10000000))

  # this is the counts of the edge vectors...
  edgecounts = Dict{(Tuple{Int, Int}), Int}()

  # now, run through each edges...
  for ei in eachindex(src)
    si = src[ei]
    di = dst[ei]
    xsi = x[si]
    ysi = y[si]
    xdi = x[di]
    ydi = y[di]

    if xsi == xdi && ysi == ydi
      continue
    end

    curve_amount = curvature*sign((di-si))

    edgelen = sqrt((xdi-xsi)^2 + (ydi-ysi)^2)
    interplen = max(ceil(Int,2*edgelen/ptscale),2)

    _add_edge_points!(xlines, ylines, xsi, ysi, xdi, ydi, curve_amount, interplen)

    if length(xlines) >= length(src)
      # compress all to edgecounts
      println("Compressing counts... $(ei) / $(length(src))")
      _xy2counts!(edgecounts, xlines, ylines, hhparams...)
      resize!(xlines,0)
      resize!(ylines,0)
    end
  end
  # handle final update
  return _xy2counts!(edgecounts, xlines, ylines, hhparams...)
end

function hexbingraphplot(src,dst,x,y;nbins=500)
  xbins = nbins
  ybins = nbins

  # compute params for hex histogram bins
  xmin,xmax = extrema(x)
  ymin,ymax = extrema(y)
  xspan, yspan = xmax - xmin, ymax - ymin
  xsize, ysize = xspan / xbins, yspan / ybins
  x0, y0 = xmin - xspan / 2,ymin - yspan / 2
  hhparams = (xsize, ysize, x0, y0) # save them as a tuple

  edgescale = 6
  hsize = min((xmax-xmin)/(xbins*edgescale),(ymax-ymin)/(ybins*edgescale))

  # need to compute hhparams

  # get coords for edges
  counts = graphedgehist(src,dst,x,y,hsize,hhparams)
  # do this for the nodes...
  #counts = HexBinPlots.xy2counts_(x,y,hhparams...)

  xh,yh,vh = HexBinPlots.counts2xy_(counts,hhparams...)
  hedges = HexBinPlots.HexHistogram(xh,yh,vh,xsize,ysize,false)

  #hx = fit(HexBinPlots.HexHistogram,x,y,nbins)
  xsh,ysh,vh = shapecoords(hedges)

  #plot(xsh, ysh, seriestype=:shape, alpha=log10.(vh)/maximum(log10.(vh)), fill_z=vh, legend=false,linealpha=0)
  return hedges
end
@time hedges = hexbingraphplot(src,dst,x,y,nbins=1000)

using ColorSchemes
xsh,ysh,vh = shapecoords(hedges)
val_range = (minimum(V[:,1]),0.008)
colors = [ColorSchemes.get(ColorSchemes.CMRmap,i,val_range) for i in V[:,1]]
# using JLD
# save("hexbin_plot_flow.jld","xsh",xsh,"ysh",ysh,"vh",vh,"x",x,"y",y,"colors",colors,"subsetnodes_c",subsetnodes_c)

# ##
# xsh,ysh,vh = shapecoords(hedges)

sizes = 3*ones(length(x))
sizes[subsetnodes_c] .= 8.0

axform = (log10.(vh).^4)/maximum(log10.(vh).^4)
avals = axform/maximum(axform)
plot(xsh, ysh, seriestype=:shape, alpha=avals, fillcolor=:darkblue, legend=false,linealpha=0, framestyle=:none,background=nothing)
scatter!(x,y,markercolor=colors,alpha=0.03,markersize=sizes,legend=false,colorbar=false,markerstrokewidth=0,background=nothing)
##
plot!(dpi=400,size=(1200,1200)) # size is inches times 100, so this should be about 3600x3600
png("hexbin_plot_flow_high_res.png")

axform = (log10.(vh).^4)/maximum(log10.(vh).^4)
avals = axform/maximum(axform)
plot(xsh, ysh, seriestype=:shape, alpha=avals, fillcolor=:darkblue, legend=false,linealpha=0, framestyle=:none,background=nothing)
scatter!(x,y,markercolor=colors,alpha=0.03,markersize=sizes,legend=false,colorbar=false,markerstrokewidth=0,background=nothing)
##
plot!(dpi=150,size=(1200,1200)) # size is inches times 100, so this should be about 3600x3600
png("hexbin_plot_flow_low_res.png")

max_stars_repo_name: dgleich/flowpaper-code

max_issues_repo_name: dgleich/flowpaper-code

max_forks_repo_name: dgleich/flowpaper-code


text:
function solve!(cur_reactor::AbstractReactor)
  isa(cur_reactor.T_bar, SymEngine.Basic) && return cur_reactor

  cur_reactor.is_good || return cur_reactor

  cur_I_P_list = solve(cur_reactor)
  cur_reactor.is_solved = true

  cur_reactor.I_P = isempty(cur_I_P_list) ? NaN : first(cur_I_P_list)
  cur_reactor.is_good = false

  valid_reactors = []
  valid_costs = []

  for cur_I_P in cur_I_P_list
    work_reactor = deepcopy(cur_reactor)

    work_reactor.I_P = cur_I_P
    work_reactor.is_good = true

    update!(work_reactor)
    work_reactor.is_valid || continue

    push!(valid_reactors, work_reactor)
    push!(valid_costs, work_reactor.cost)
  end

  if !isempty(valid_reactors)
    cur_reactor.I_P = valid_reactors[indmin(valid_costs)].I_P
    cur_reactor.is_good = true
  else
    cur_reactor.is_valid = false
  end

  update!(cur_reactor)

  cur_reactor
end

function solve(cur_reactor::AbstractReactor)
  @assert !isa(cur_reactor.T_bar, SymEngine.Basic)

  cur_reactor.T_bar <= 0 && return []

  tmp_reactor = deepcopy(cur_reactor)

  tmp_reactor.sigma_v = calc_sigma_v(tmp_reactor)
  isnan(tmp_reactor.sigma_v) && return []

  cur_equation = calc_I_P(tmp_reactor)

  isa(cur_equation, SymEngine.Basic) || return [cur_equation]

  secondary_constraint = getfield(Fussy, Symbol("$(tmp_reactor.constraint)_equation"))(tmp_reactor)

  cur_dict = equation_set_dict(EquationSet(tmp_reactor, secondary_constraint))

  any(is_nan, values(cur_dict)) && return []

  @assert all(is_finite, values(cur_dict))

  cur_equation = subs(
    cur_equation, cur_dict...
  )

  cur_equation /= symbols(:I_P)

  cur_equation -= 1.0

  _solve(tmp_reactor, cur_equation)
end

function _solve(cur_reactor::AbstractReactor, cur_equation::SymEngine.Basic)
  cur_atol = 10 * eps()

  cur_I_P_list = find_zeros(
    (tmp_value) -> real(float(cur_equation(tmp_value))),
    min_I_P, max_I_P, rtol=1e-2
  )

  isempty(cur_I_P_list) && return []

  good_reactors = Vector{AbstractReactor}()

  for (cur_index, cur_I_P) in enumerate(cur_I_P_list)
    isapprox(0.0, cur_I_P, atol=cur_atol) && continue

    tmp_reactor = deepcopy(cur_reactor)

    tmp_reactor.I_P = cur_I_P

    secondary_constraint = getfield(Fussy, Symbol("$(tmp_reactor.constraint)_equation"))(tmp_reactor)

    cur_equation_set = EquationSet(tmp_reactor, secondary_constraint)

    try

      tmp_reactor.B_0 = convert(Real, cur_equation_set.B_0)
      tmp_reactor.R_0 = convert(Real, cur_equation_set.R_0)

      tmp_reactor.n_bar = convert(Real, calc_n_bar(tmp_reactor))

    catch cur_error

      is_caught_error = isa(cur_error, InexactError)
      is_caught_error |= isa(cur_error, DomainError)
      is_caught_error || rethrow(cur_error)

      continue

    end

    ( 0 <= f_BS(tmp_reactor) <= 1 ) || continue
    ( 0 <= f_CD(tmp_reactor) <= 1 ) || continue

    ( f_IN(tmp_reactor) > -0.01 ) || continue
    ( f_IN(tmp_reactor) <= 1.00 ) || continue

    f_total = f_BS(tmp_reactor) + f_CD(tmp_reactor) + f_IN(tmp_reactor)

    isapprox(1.0, f_total, atol=cur_atol) || continue

    push!(good_reactors, tmp_reactor)
  end

  isempty(good_reactors) && return []

  real_root_count = length(cur_I_P_list)

  @assert real_root_count <= max_roots

  cur_I_P_list = map(tmp_reactor -> tmp_reactor.I_P, good_reactors)

  cur_I_P_list
end


text:
module Turing

##############
# Dependency #
########################################################################
# NOTE: when using anything from external packages,                    #
#       let's keep the practice of explictly writing Package.something #
#       to indicate that's not implemented inside Turing.jl            #
########################################################################

using Requires, Reexport, ForwardDiff
using DistributionsAD, Bijectors, StatsFuns, SpecialFunctions
using Statistics, LinearAlgebra
using Markdown, Libtask
@reexport using Distributions, MCMCChains, Libtask
using Tracker: Tracker

import Base: ~, ==, convert, hash, promote_rule, rand, getindex, setindex!
import DynamicPPL: getspace

const PROGRESS = Ref(true)
function turnprogress(switch::Bool)
    @info "[Turing]: progress logging is $(switch ? "enabled" : "disabled") globally"
    PROGRESS[] = switch
end

const DEBUG = Bool(parse(Int, get(ENV, "DEBUG_TURING", "0")))

# Random probability measures.
include("stdlib/distributions.jl")
include("stdlib/RandomMeasures.jl")
include("utilities/Utilities.jl")
using .Utilities
include("core/Core.jl")
using .Core
include("inference/Inference.jl")  # inference algorithms
using .Inference
include("variational/VariationalInference.jl")
using .Variational

# TODO: re-design `sample` interface in MCMCChains, which unify CmdStan and Turing.
#   Related: https://github.com/TuringLang/Turing.jl/issues/746
#@init @require CmdStan="593b3428-ca2f-500c-ae53-031589ec8ddd" @eval begin
#     @eval Utilities begin
#         using ..Turing.CmdStan: CmdStan, Adapt, Hmc
#         using ..Turing: HMC, HMCDA, NUTS
#         include("utilities/stan-interface.jl")
#     end
# end

@init @require DynamicHMC="bbc10e6e-7c05-544b-b16e-64fede858acb" @eval Inference begin
    import ..DynamicHMC

    if isdefined(DynamicHMC, :mcmc_with_warmup)
        using ..DynamicHMC: mcmc_with_warmup
        include("contrib/inference/dynamichmc.jl")
    else
        error("Please update DynamicHMC, v1.x is no longer supported")
    end
end

###########
# Exports #
###########

# Turing essentials - modelling macros and inference algorithms
export  @model,                 # modelling
        @varname,
        DynamicPPL,

        MH,                     # classic sampling
        RWMH,
        ESS,
        Gibbs,

        HMC,                    # Hamiltonian-like sampling
        SGLD,
        SGHMC,
        HMCDA,
        NUTS,
        DynamicNUTS,
        ANUTS,

        IS,                     # particle-based sampling
        SMC,
        CSMC,
        PG,

        vi,                    # variational inference
        ADVI,

        sample,                 # inference
        psample,
        setchunksize,
        resume,
        @logprob_str,
        @prob_str,

        auto_tune_chunk_size!,  # helper
        setadbackend,
        setadsafe,

        turnprogress,           # debugging

        Flat,
        FlatPos,
        BinomialLogit,
        VecBinomialLogit,
        OrderedLogistic,
        LogPoisson,
        NamedDist,
        filldist,
        arraydist

# Reexports
using AbstractMCMC: sample, psample
export sample, psample

end


text:
"""
    abstract type Grid{Dc,Dp} <: Triangulation{Dc,Dp}

Abstract type that represents conforming triangulations, whose cell-wise
nodal coordinates are defined with a vector of nodal coordinates, plus
a cell-wise vector of node ids.

The interface of `Grid` is defined by overloading the
methods in `Triangulation` plus the following ones:

- [`get_node_coordinates(trian::Grid)`](@ref)
- [`get_cell_nodes(trian::Grid)`](@ref)

From these two methods a default implementation of [`get_cell_coordinates(trian::Triangulation)`](@ref)
is available.

The `Grid`  interface has the following traits

- [`OrientationStyle(::Type{<:Grid})`](@ref)
- [`RegularityStyle(::Type{<:Grid})`](@ref)

The interface of `Grid` is tested with

- [`test_grid`](@ref)

"""
abstract type Grid{Dc,Dp} <: Triangulation{Dc,Dp} end

# Traits

"""
    OrientationStyle(::Type{<:Grid}) -> Val{Bool}
    OrientationStyle(::Grid) -> Val{Bool}

`Val{true}()` if has oriented faces, `Val{false}()` otherwise (default).
"""
OrientationStyle(a::Grid) = OrientationStyle(typeof(a))
OrientationStyle(::Type{<:Grid}) = Val{false}()

"""
    RegularityStyle(::Type{<:Grid}) -> Val{Bool}
    RegularityStyle(::Grid) -> Val{Bool}

`Val{true}()` if no hanging-nodes (refault), `Val{false}()` otherwise.
"""
RegularityStyle(a::Grid) = RegularityStyle(typeof(a))
RegularityStyle(::Type{<:Grid}) = Val{true}()

# Interface

"""
    get_node_coordinates(trian::Grid) -> AbstractArray{<:Point{Dp}}
"""
function get_node_coordinates(trian::Grid)
  @abstractmethod
end

"""
    get_cell_nodes(trian::Grid)
"""
function get_cell_nodes(trian::Grid)
  @abstractmethod
end

"""
    test_grid(trian::Grid)
"""
function test_grid(trian::Grid)
  test_triangulation(trian)
  nodes_coords = get_node_coordinates(trian)
  @test isa(nodes_coords,AbstractArray{<:Point})
  cell_nodes = get_cell_nodes(trian)
  @test isa(cell_nodes,AbstractArray{<:AbstractArray{<:Integer}})
  @test num_nodes(trian) == length(nodes_coords)
  @test isa(is_oriented(trian),Bool)
  @test isa(is_regular(trian),Bool)
  @test OrientationStyle(trian) in (Val{false}(), Val{true}())
  @test RegularityStyle(trian) in (Val{false}(), Val{true}())
end

# Methods from triangulation

function get_cell_coordinates(trian::Grid)
  node_to_coords = get_node_coordinates(trian)
  cell_to_nodes = get_cell_nodes(trian)
  LocalToGlobalArray(cell_to_nodes,node_to_coords)
end

# Some API

"""
    num_nodes(trian::Grid) -> Int
"""
num_nodes(trian::Grid) = length(get_node_coordinates(trian))

"""
    is_oriented(::Type{<:Grid}) -> Bool
    is_oriented(a::Grid) -> Bool
"""
is_oriented(a::Grid) = _is_oriented(OrientationStyle(a))
is_oriented(a::Type{<:Grid}) = _is_oriented(OrientationStyle(a))

"""
    is_regular(::Type{<:Grid}) -> Bool
    is_regular(a::Grid) -> Bool
"""
is_regular(a::Grid) = _is_regular(RegularityStyle(a))
is_regular(a::Type{<:Grid}) = _is_regular(RegularityStyle(a))

"""
    Grid(reffe::LagrangianRefFE)
"""
function Grid(reffe::LagrangianRefFE)
  UnstructuredGrid(reffe)
end

"""
    compute_linear_grid(reffe::LagrangianRefFE)
"""
function compute_linear_grid(reffe::LagrangianRefFE)
  @notimplementedif ! is_n_cube(get_polytope(reffe)) "linear grid only implemented for n-cubes at the moment"
  if get_order(reffe) == 0
    D = num_cell_dims(reffe)
    partition = tfill(1,Val{D}())
  else
    partition = get_orders(reffe)
  end
  pmin, pmax = get_bounding_box(get_polytope(reffe))
  desc = CartesianDescriptor(pmin,pmax,partition)
  CartesianGrid(desc)
end

"""
    compute_reference_grid(p::LagrangianRefFE, nelems::Integer)
"""
function compute_reference_grid(reffe::LagrangianRefFE, nelems::Integer)
  p = get_polytope(reffe)
  r = LagrangianRefFE(Float64,p,nelems)
  compute_linear_grid(r)
end

"""
    Grid(::Type{ReferenceFE{d}},p::Polytope) where d
"""
function Grid(::Type{ReferenceFE{d}},p::Polytope) where d
  UnstructuredGrid(ReferenceFE{d},p)
end

"""
    GridTopology(grid::Grid)
    GridTopology(grid::Grid, cell_to_vertices::Table, vertex_to_node::Vector)
"""
function GridTopology(grid::Grid)
  _grid = UnstructuredGrid(grid)
  UnstructuredGridTopology(_grid)
end

function GridTopology(grid::Grid, cell_to_vertices::Table, vertex_to_node::AbstractVector)
  _grid = UnstructuredGrid(grid)
  UnstructuredGridTopology(_grid,cell_to_vertices,vertex_to_node)
end

"""
    simplexify(grid::Grid)
"""
function simplexify(grid::Grid)
  simplexify(UnstructuredGrid(grid))
end