text:
# Weighted sum and related reduction function

macro code_wsumfuns(AN, fname)
    h = codegen_helper(AN)

    t1 = h.term(1)
    ti = h.term(:i)
    tidx = h.term(:idx)

    fname! = symbol("$(fname)!")
    _fname! = symbol("_$(fname)!")
    aparams = h.contiguous_aparams
    args = h.args

    quote
        global $(fname)
        function ($fname)(w::ContiguousRealArray, $(aparams...))
            n::Int = $(h.inputlen)
            n == length(w) || error("Inconsistent argument dimensions.")

            if n == 0
                0.0
            else
                @inbounds s = $(t1) * w[1]
                for i = 2 : n
                    @inbounds s += ($ti) * w[i]
                end
                s
            end
        end

        global $(_fname!)
        function $(_fname!)(dst::ContiguousArray, m::Int, n::Int, k::Int, w::ContiguousRealArray, $(aparams...))

            if n == 1  # each page has a single row (simply evaluate)
                @inbounds w1 = w[1]
                for i = 1:m*k
                    @inbounds dst[i] = w1 * $(ti)
                end

            elseif m == 1  # each page has a single column
                idx = 0
                for l = 1:k
                    idx += 1
                    @inbounds s = w[1] * ($tidx)
                    for j = 2:n
                        idx += 1
                        @inbounds s += w[j] * ($tidx)
                    end
                    @inbounds dst[l] = s
                end

            elseif k == 1 # only one page
                @inbounds w1 = w[1]
                for i = 1:m
                    @inbounds dst[i] = w1 * $(ti)
                end
                idx = m
                for j = 2:n
                    @inbounds wj = w[j]
                    for i = 1:m
                        idx += 1
                        @inbounds dst[i] += wj * $(tidx)
                    end
                end

            else  # multiple generic pages
                idx = 0
                od = 0
                @inbounds w1 = w[1]
                for l = 1:k
                    for i = 1:m
                        idx += 1
                        @inbounds dst[od+i] = w1 * $(tidx)
                    end
                    for j = 2:n
                        @inbounds wj = w[j]
                        for i = 1:m
                            idx += 1
                            odi = od + i
                            @inbounds dst[odi] += wj * $(tidx)
                        end
                    end
                    od += m
                end
            end

            return dst
        end

        global $(fname!)
        function $(fname!)(dst::ContiguousArray, w::ContiguousRealArray, $(aparams...), dim::Int)
            shp = $(h.inputsize)
            1 <= dim <= length(shp) || error("Invalid value of dim.")
            length(w) == shp[dim] || error("Inconsistent argument dimensions.")
            m = prec_length(shp, dim)
            n = shp[dim]
            k = succ_length(shp, dim)
            $(_fname!)(dst, m, n, k, w, $(args...))
        end

        global $(fname)
        function ($fname)(w::ContiguousRealArray, $(aparams...), dim::Int)
            tt = $(h.termtype)
            shp = $(h.inputsize)
            r = Array(promote_type(tt, eltype(w)), Base.reduced_dims(shp, dim))
            ($fname!)(r, w, $(args...), dim)
        end
    end
end

@code_wsumfuns 0 wsum
@code_wsumfuns 1 wsum
@code_wsumfuns 2 wsum
@code_wsumfuns 3 wsum
@code_wsumfuns (-2) wsumfdiff

# Specialized cases

wsum{T<:BlasFP}(w::Array{T}, x::Array{T}) = Base.BLAS.dot(w, x)

function wsum!{T<:BlasFP}(dst::Array{T}, w::Array{T}, x::Vector{T}, dim::Int)
    dst[1] = wsum(x, w)
    dst
end

function wsum!{T<:BlasFP}(dst::Array{T}, w::Array{T}, x::Matrix{T}, dim::Int)
    if dim == 1
        Base.BLAS.gemv!('T', one(T), x, vec(w), zero(T), vec(dst))
    elseif dim == 2
        Base.BLAS.gemv!('N', one(T), x, vec(w), zero(T), vec(dst))
    else
        error("dim must be either 1 or 2.")
    end
    dst
end

function wsum!{T<:BlasFP}(dst::Array{T}, w::Array{T}, x::Array{T}, dim::Int)
    siz = size(x)
    nd = length(siz)
    rd = siz[dim]  # this ensures 1 <= dim <= nd

    if dim == 1
        rx = reshape(x, siz[1], succ_length(siz, 1))
        Base.BLAS.gemv!('T', one(T), rx, vec(w), zero(T), vec(dst))
    elseif dim == nd
        rx = reshape(x, prec_length(siz, nd), siz[nd])
        Base.BLAS.gemv!('N', one(T), rx, vec(w), zero(T), vec(dst))
    else
        m::Int = prec_length(siz, dim)
        n::Int = siz[dim]
        k::Int = succ_length(siz, dim)
        plen = m * n

        vw = vec(w)
        for l in 1 : k
            sx = pointer_to_array(pointer(x, plen * (l - 1) + 1), (m, n))
            sdst = pointer_to_array(pointer(dst, m * (l - 1) + 1), m)
            Base.BLAS.gemv!('N', one(T), sx, vw, zero(T), sdst)
        end
    end
    dst
end

# Convenience functions

wsumabs(w::ContiguousRealArray, x::ContiguousArray) = wsum(w, AbsFun(), x)
wsumabs!(dst::ContiguousRealArray, w::ContiguousRealArray, x::ContiguousArray, dim::Int) = wsum!(dst, w, AbsFun(), x, dim)
wsumabs(w::ContiguousRealArray, x::ContiguousArray, dim::Int) = wsum(w, AbsFun(), x, dim)

wsumabsdiff(w::ContiguousRealArray, x::ContiguousArray, y::ContiguousArrOrNum) = wsumfdiff(w, AbsFun(), x, y)
wsumabsdiff!(dst::ContiguousRealArray, w::ContiguousRealArray, x::ContiguousArray, y::ContiguousArrOrNum, dim::Int) = wsumfdiff!(dst, w, AbsFun(), x, y, dim)
wsumabsdiff(w::ContiguousRealArray, x::ContiguousArray, y::ContiguousArrOrNum, dim::Int) = wsumfdiff(w, AbsFun(), x, y, dim)

wsumsq(w::ContiguousRealArray, x::ContiguousArray) = wsum(w, Abs2Fun(), x)
wsumsq!(dst::ContiguousRealArray, w::ContiguousRealArray, x::ContiguousArray, dim::Int) = wsum!(dst, w, Abs2Fun(), x, dim)
wsumsq(w::ContiguousRealArray, x::ContiguousArray, dim::Int) = wsum(w, Abs2Fun(), x, dim)

wsumsqdiff(w::ContiguousRealArray, x::ContiguousArray, y::ContiguousArrOrNum) = wsumfdiff(w, Abs2Fun(), x, y)
wsumsqdiff!(dst::ContiguousRealArray, w::ContiguousRealArray, x::ContiguousArray, y::ContiguousArrOrNum, dim::Int) = wsumfdiff!(dst, w, Abs2Fun(), x, y, dim)
wsumsqdiff(w::ContiguousRealArray, x::ContiguousArray, y::ContiguousArrOrNum, dim::Int) = wsumfdiff(w, Abs2Fun(), x, y, dim)


text:
abstract type AbstractAlgorithm end

# trait dispatch
RatesCache(::AbstractAlgorithm) = HasRates()

# by default, initialization is the same as updating the rates vector
# algorithms should overload this behavior as needed
@inline function initialize!(algorithm::AbstractAlgorithm, state, model)
  update_jump_rates!(algorithm, state, model)
end

# we assume that `s` is a time increment rather than an absolute time
# as such, the default behavior is to return `t + s` as the next system time
# algorithms such as NRM can overload this behavior
@inline @fastmath function get_new_time(algorithm::AbstractAlgorithm, t, s)
  return t + s
end

abstract type DecoupledSearchAlgorithm <: AbstractAlgorithm end

# delegate the update procedure based on the RatesCacheType
@inline function update_jump_rates!(algorithm::DecoupledSearchAlgorithm, state, model)
  algorithm.total_rate = update_jump_rates!(RatesCache(algorithm), algorithm.rates, algorithm.total_rate, state, model)
end

@inline function update_jump_rates!(algorithm::DecoupledSearchAlgorithm, state, model, j)
  algorithm.total_rate = update_jump_rates!(RatesCache(algorithm), algorithm.rates, algorithm.total_rate, state, model, j)
end

# A DecoupledSearchSampler always generates the next reaction time based on an exponential distribution.
# The next reaction index depends exclusively on the sampler.
@inline function generate_jump(algorithm::DecoupledSearchAlgorithm)
  j = search_jumps(algorithm)
  s = randexp() / algorithm.total_rate

  return (j, s)
end

abstract type CoupledSearchAlgorithm <: AbstractAlgorithm end

# delegate the update procedure based on the RatesCacheType
@inline function update_jump_rates!(algorithm::CoupledSearchAlgorithm, state, model)
  algorithm.total_rate = update_jump_rates!(RatesCache(algorithm), algorithm.rates, algorithm.total_rate, state, model)
end

@inline function update_jump_rates!(algorithm::CoupledSearchAlgorithm, state, model, j)
  algorithm.total_rate = update_rates!(RatesCache(algorithm), algorithm.rates, algorithm.total_rate, state, model, j)
end

@inline function generate_jump(algorithm::CoupledSearchAlgorithm)
  (j, s) = search_jumps(algorithm)

  return (j, s)
end

abstract type UnsafeLeapAlgorithm <: AbstractAlgorithm end

@inline function update_jump_rates!(algorithm::UnsafeLeapAlgorithm, state, model)
  algorithm.total_rate = update_jump_rates!(RatesCache(algorithm), algorithm.rates, algorithm.total_rate, state, model)
end

function initialize!(algorithm::UnsafeLeapAlgorithm, state, model)
  update!(algorithm, state, model)
end

function generate_leap!(v, algorithm::UnsafeLeapAlgorithm, Δt_max)
  # unpack data
  rates = algorithm.rates
  total_rate = algorithm.total_rate
  validate_leap! = algorithm.validate_leap!

  # evaluate the leap formula to compute τ
  s = algorithm.leap_formula(rates, total_rate)
  s = min(s, Δt_max)

  # if there is a valid leap...
  if isfinite(s)
    generate_events!(algorithm, v, rates, s)
    validate_leap!(v, s)
  else
    # otherwise fill the jumps with zeros
    fill!(v, zero(eltype(v)))
  end

  return v, s
end

# abstract type SafeLeapAlgorithm <: AbstractAlgorithm end
#
# @inline function update_jump_rates!(algorithm::SafeLeapAlgorithm, state, model)
#   algorithm.total_rate = update_jump_rates!(RatesCache(algorithm), algorithm.rates, algorithm.total_rate, state, model)
# end
#
# function initialize!(algorithm::SafeLeapAlgorithm, state, model)
#   update!(algorithm, state, model)
# end
#
# function generate_leap!(v, algorithm::SafeLeapAlgorithm, Δt_max)
#   # unpack data
#   rates = algorithm.rates
#   total_rate = algorithm.total_rate
#
#   # evaluate the leap formula to compute τ
#   s = algorithm.leap_formula(rates, total_rate)
#   s = min(s, Δt_max)
#
#   # if there is a valid leap...
#   if isfinite(s)
#     generate_events!(algorithm, v, rates, s)
#   else
#     # otherwise fill the jumps with zeros
#     fill!(v, zero(eltype(v)))
#   end
#
#   return v, s
# end

ext: jl


text:
@muladd function perform_step!(integrator,cache::LambaEMConstantCache)
  @unpack t,dt,uprev,u,W,p,f = integrator
  du1 = integrator.f(uprev,p,t)
  K = @.. uprev + dt * du1

  if is_split_step(integrator.alg)
    L = integrator.g(uprev,p,t+dt)
  else
    L = integrator.g(K,p,t+dt)
  end

  mil_correction = zero(u)
  if !is_diagonal_noise(integrator.sol.prob) || typeof(W.dW) <: Number
      noise = L*W.dW
  else
      noise = L.*W.dW
  end

  u = K+noise

  if integrator.opts.adaptive
    du2 = integrator.f(K,p,t+dt)
    Ed = dt*(du2 - du1)/2

    utilde =  K + noise*integrator.sqdt #L*integrator.sqdt
    ggprime = (integrator.g(utilde,p,t).-L)./(integrator.sqdt)
    if !is_diagonal_noise(integrator.sol.prob) || typeof(W.dW) <: Number
        En = ggprime*(W.dW.^2 .- dt)./2
    else
        En = ggprime.*(W.dW.^2 .- dt)./2
    end
    resids = calculate_residuals(Ed, En, uprev, u, integrator.opts.abstol,
                                 integrator.opts.reltol, integrator.opts.delta,
                                 integrator.opts.internalnorm, t)
    integrator.EEst = integrator.opts.internalnorm(resids, t)
  end

  integrator.u = u
end

@muladd function perform_step!(integrator,cache::LambaEMCache)
  @unpack du1,du2,K,tmp,L,gtmp,dW_cache = cache
  @unpack t,dt,uprev,u,W,p,f = integrator

  integrator.f(du1,uprev,p,t)
  @.. K = uprev + dt * du1

  if is_split_step(integrator.alg)
    integrator.g(L,K,p,t+dt)
  else
    integrator.g(L,uprev,p,t+dt)
  end

  if is_diagonal_noise(integrator.sol.prob)
    @.. tmp=L*W.dW
  else
    mul!(tmp,L,W.dW)
  end

  @.. u = K+tmp

  if integrator.opts.adaptive

    if !is_diagonal_noise(integrator.sol.prob)
      g_sized = norm(L,2)
    else
      g_sized = L
    end

    if !is_diagonal_noise(integrator.sol.prob)
      @.. tmp = K + integrator.sqdt * g_sized
      integrator.g(gtmp,tmp,p,t)
      g_sized2 = norm(gtmp,2)
      @.. dW_cache = W.dW.^2 - dt
      diff_tmp = integrator.opts.internalnorm(dW_cache, t)
      En = (g_sized2-g_sized)/(2integrator.sqdt)*diff_tmp
      @.. tmp = En
    else
      @.. tmp = K + integrator.sqdt * L
      integrator.g(gtmp,tmp,p,t)
      @.. tmp = (gtmp-L)/(2integrator.sqdt)*(W.dW.^2 - dt)
    end

    # Ed
    integrator.f(du2,K,p,t+dt)
    @.. tmp += integrator.opts.internalnorm(integrator.opts.delta * dt * (du2 - du1) / 2, t)

    calculate_residuals!(tmp, tmp, uprev, u, integrator.opts.abstol,
                         integrator.opts.reltol, integrator.opts.internalnorm, t)
    integrator.EEst = integrator.opts.internalnorm(tmp,t)
  end
end

@muladd function perform_step!(integrator,cache::LambaEulerHeunConstantCache)
  @unpack t,dt,uprev,u,W,p,f = integrator
  du1 = integrator.f(uprev,p,t)
  K = uprev + dt * du1
  L = integrator.g(uprev,p,t)

  if is_diagonal_noise(integrator.sol.prob)
    noise = L.*W.dW
  else
    noise = L*W.dW
  end
  tmp = K .+ noise
  gtmp2 = (1/2).*(L.+integrator.g(tmp,p,t+dt))
  if is_diagonal_noise(integrator.sol.prob)
    noise2 = gtmp2.*W.dW
  else
    noise2 = gtmp2*W.dW
  end

  u = uprev .+ (dt / 2) .* (du1 .+ integrator.f(tmp,p,t+dt)) .+ noise2

  if integrator.opts.adaptive
    du2 = integrator.f(K,p,t+dt)
    Ed = dt*(du2 - du1)/2

    utilde = uprev + L*integrator.sqdt
    ggprime = (integrator.g(utilde,p,t).-L)./(integrator.sqdt)
    En = ggprime.*(W.dW.^2)./2

    resids = calculate_residuals(Ed, En, uprev, u, integrator.opts.abstol,
                                 integrator.opts.reltol, integrator.opts.delta,
                                 integrator.opts.internalnorm, t)
    integrator.EEst = integrator.opts.internalnorm(resids, t)
  end

  integrator.u = u
end

@muladd function perform_step!(integrator,cache::LambaEulerHeunCache)
  @unpack du1,du2,K,tmp,L,gtmp,dW_cache = cache
  @unpack t,dt,uprev,u,W,p,f = integrator
  integrator.f(du1,uprev,p,t)
  integrator.g(L,uprev,p,t)
  @.. K = uprev + dt * du1

  if is_diagonal_noise(integrator.sol.prob)
    @.. tmp=L*W.dW
  else
    mul!(tmp,L,W.dW)
  end

  @.. tmp = K+tmp

  integrator.f(du2,tmp,p,t+dt)
  integrator.g(gtmp,tmp,p,t+dt)

  if is_diagonal_noise(integrator.sol.prob)
    @.. tmp=(1/2)*W.dW*(L+gtmp)
  else
    @.. gtmp = (1/2)*(L+gtmp)
    mul!(tmp,gtmp,W.dW)
  end

  dto2 = dt*(1/2)
  @.. u = uprev + dto2*(du1+du2) + tmp

  if integrator.opts.adaptive

    if !is_diagonal_noise(integrator.sol.prob)
      g_sized = norm(L,2)
    else
      g_sized = L
    end

    if !is_diagonal_noise(integrator.sol.prob)
      @.. tmp = uprev + integrator.sqdt * g_sized
      integrator.g(gtmp,tmp,p,t)
      g_sized2 = norm(gtmp,2)
      @.. dW_cache = W.dW.^2
      diff_tmp = integrator.opts.internalnorm(dW_cache, t)
      En = (g_sized2-g_sized)/(2integrator.sqdt)*diff_tmp
      @.. tmp = En
    else
      @.. tmp = uprev + integrator.sqdt * L
      integrator.g(gtmp,tmp,p,t)
      @.. tmp = (gtmp-L)/(2integrator.sqdt)*(W.dW.^2)
    end

    # Ed
    integrator.f(du2,K,p,t+dt)
    @.. tmp += integrator.opts.internalnorm(integrator.opts.delta * dt * (du2 - du1) / 2, t)

    calculate_residuals!(tmp, tmp, uprev, u, integrator.opts.abstol,
                         integrator.opts.reltol, integrator.opts.internalnorm, t)
    integrator.EEst = integrator.opts.internalnorm(tmp, t)
  end
end

max_stars_repo_name: JuliaDiffEq/StochasticDiffEq


text:
module Runtests

using Test, PeriodicSystems

@testset "Test PeriodicSystems" begin
# test constructors, basic tools
include("test_psutils.jl")
include("test_conversions.jl")
include("test_pslifting.jl")
include("test_psanalysis.jl")
end

end


text:
#loop fix example
using Base64, Profile, Test

function sum_sqrt_n(n)
    r = 0.0 # this is initialized with float instead of int to solve instability
    for i = 1:n
        r = r +sqrt(i)
    end
    return r
end


text:
function spgr1()

    spgr! = SPGRBlochSim(10, 5, deg2rad(20), Val(1))
    s = Spin(1, 1000, 100, 0)
    spgr!(s)

    show(devnull, spgr!)
    show(devnull, "text/plain", spgr!)
    return true

end

function spgr2()

    spgr1! = SPGRBlochSim(10, 5, RF([0.1], 0.01), RFSpoiling(), Val(2), Val(true))
    s1 = Spin(1, 1000, 100, 0)
    s2 = SpinMC(1, (1, 0), (1000, 100), (100, 10), (0, 0), (Inf, Inf))
    M11 = spgr1!(s1)
    M12 = spgr1!(s2)

    spgr2! = SPGRBlochSim(10, 5, RF([0.1], 0.01), RFandGradientSpoiling((0, 0, 0), 0), Val(2), Val(true))
    s1 = Spin(1, 1000, 100, 0)
    s2 = SpinMC(1, (1, 0), (1000, 100), (100, 10), (0, 0), (Inf, Inf))
    M21 = spgr2!(s1)
    M22 = spgr2!(s2)

    return all(M11[i] ≈ M12[i][1] ≈ M21[i] ≈ M22[i][1] for i = 1:2)

end

@testset "SPGR" begin

    @test spgr1()
    @test spgr2()

end


text:
using Statistics
using FFTW

function DHC(image, filter_set; norm_on = 1, coeff_12_on =1, coeff_20_on = 1)
    (Nx, Ny) = size(image)

    (J,L,_,_) = size(fink_filter_set)

    out_coeff = []

    ## Coeff at (1,0)
    #takes care of image normalization
    S0 = zeros(2)
    if norm_on == 1
        S0[1] = mean(image)
        norm_im = image.-S0[1]
        S0[2] = sum(abs2.(norm_im))/(Nx*Ny)
        norm_im = norm_im./sqrt(Nx*Ny*S0[2])
    else
        norm_im = image
    end

    append!(out_coeff,S0)

    ##Coeff at (1,1)
    #store the fft of images and im*filters for later
    im_fd_0 = fft(norm_im)

    S1 = zeros((J,L))
    im_fd_0_1 = Array{ComplexF64,4}(undef, J, L, Nx, Ny)

    for j = 1:J
        for l = 1:L
            im_fd_0_1[j,l,:,:] = im_fd_0.*filter_set[j,l,:,:] #wavelet already in fft domain not shifted
            S1[j,l]=sum(abs2.(im_fd_0_1[j,l,:,:])) #normalization choice arb to make order unity
        end
    end
    append!(out_coeff,S1)

    ##Coeff at (2,0)
    if coeff_20_on == 1
        S20 = Array{Float64,4}(undef, J, J, L, L)

        im_rd_0_1 = Array{ComplexF64,4}(undef, J, L, Nx, Ny)
        for j1 = 1:J
            for l1 = 1:L
                im_rd_0_1[j1,l1,:,:] = ifft(im_fd_0_1[j1,l1,:,:])
            end
        end

        for j1 = 1:J
            for l1 = 1:L
                for l2 = 1:L
                    for j2  = 1:J
                        S20[j1,j2,l1,l2] = sum(
                            abs.(
                                im_rd_0_1[j1,l1,:,:].*
                                conj.(im_rd_0_1[j2,l2,:,:])
                            )
                        )
                    end
                end
            end
        end

        append!(out_coeff,S20)
    end

    ##Coeff at (1,2) level
    if coeff_12_on == 1
        S12 = Array{Float64,4}(undef, J, J, L, L)
        for j1 = 1:J
            for l1 = 1:L
                not_rot_im = fftshift(im_fd_0_1[j1,l1,:,:])
                for l2 = 1:L
                    for j2 = 1:J
                        S12[j1,j2,l1,l2] = sum(
                            abs.(
                                not_rot_im.*
                                conj(fftshift(im_fd_0_1[j2,l2,:,:]))
                            )
                        )
                    end
                end
            end
        end
    append!(out_coeff,S12)
    end
end

using Plots

##
heatmap(fftshift(fink_filter))

##
function finklet(j, l)
    # -------- set filters
    jrad = 7-j
    dθ = π/8        # 8 angular bins hardwired
    θ_l = dθ*l
    # -------- define coordinates
    nx = 256
    xbox = LinRange(-nx/2, nx/2-1 , nx)
    # make a 256x256 grid of X
    sx = xbox' .* ones(nx)
    sy = ones(nx)' .* xbox
    r  = sqrt.((sx).^2 + (sy).^2)
    θ  = mod.(atan.(sy, sx).+π .-θ_l,2*π)
    nozeros = r .> 0
    logr = log2.(r[nozeros])
    r[nozeros] = logr
    # -------- in Fourier plane, envelope of psi_j,l
    mask = (abs.(θ.-π).<= dθ) .& (abs.(r.-jrad) .<= 1)
    # -------- angular part
    ang = cos.((θ.-π).*4)
    # -------- radial part
    rad = cos.((r.-jrad).*π./2)
    psi = mask.*ang.*rad             #mask times angular part times radial part
    return psi
end

function fink_filter_bank(J,L)
    fink_filter = Array{Float64, 4}(undef, J, L, 256, 256)
    for j = 1:J-1
        for l = 1:L
            fink_filter[j,l,:,:]=fftshift(finklet(j-1,l-1))
        end
    end
    return fink_filter
end

fink_filter_set = fink_filter_bank(8,8)

size(fink_filter_set)

##
using BenchmarkTools

DHC(fink_filter,fink_filter_set)

##No FFT
@benchmark DHC(fink_filter,fink_filter_set)

##FFT
@benchmark DHC(fink_filter,fink_filter_set)

##FFT + Zeros
DHC(fink_filter,fink_filter_set)

@benchmark DHC(fink_filter,fink_filter_set)

##FFT S1 Simplest mat mult
DHC(fink_filter,fink_filter_set)

@benchmark DHC(fink_filter,fink_filter_set)

## Need to switch to 8,8

##FFT S1 S20
DHC(fink_filter_set[1,3,:,:],fink_filter_set)

20.8/(8*8*8*8+8*8+2)
8/1857

@benchmark DHC(fink_filter_set[1,3,:,:],fink_filter_set)

##FFT S1 S20 S21
DHC(fink_filter_set[1,3,:,:],fink_filter_set)

@benchmark DHC(fink_filter_set[1,3,:,:],fink_filter_set)

39.348/(2*8*8*8*8+8*8+2)

## Let me benchmark on S1 only first.

@benchmark DHC(fink_filter_set[1,3,:,:],fink_filter_set,coeff_12_on =0, coeff_20_on = 0)

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

function DHC(image, filter_set; norm_on = 1, coeff_12_on =1, coeff_20_on = 1)
    (Nx, Ny) = size(image)

    (J,L,_,_) = size(fink_filter_set)

    out_coeff = []

    ## Coeff at (1,0)
    #takes care of image normalization
    S0 = zeros(2)
    if norm_on == 1
        S0[1] = mean(image)
        norm_im = image.-S0[1]
        S0[2] = sum(abs2.(norm_im))/(Nx*Ny)
        norm_im = norm_im./sqrt(Nx*Ny*S0[2])
    else
        norm_im = image
    end

    append!(out_coeff,S0)

    ##Coeff at (1,1)
    #store the fft of images and im*filters for later
    im_fd_0 = fft(norm_im)

    S1 = zeros((J,L))
    im_fd_0_1 = repeat(im_fd_0,inner=[1,1,J,L])

    for j = 1:J
        for l = 1:L
            had!(im_fd_0_1[:,:,j,l],filter_set[j,l,:,:]) #wavelet already in fft domain not shifted
            S1[j,l]=sum(abs2.(im_fd_0_1[:,:,j,l])) #normalization choice arb to make order unity
        end
    end
    append!(out_coeff,S1)

    ##Coeff at (2,0)
    if coeff_20_on == 1
        S20 = Array{Float64,4}(undef, J, J, L, L)

        im_rd_0_1 = Array{ComplexF64,4}(undef, J, L, Nx, Ny)
        for j1 = 1:J
            for l1 = 1:L
                im_rd_0_1[j1,l1,:,:] = ifft(im_fd_0_1[:,:,j1,l1])
            end
        end

        # for j1 = 1:J
        #     for l1 = 1:L
        #         for l2 = 1:L
        #             for j2  = 1:J
        #                 S20[j1,j2,l1,l2] = sum(
        #                     abs.(
        #                         im_rd_0_1[j1,l1,:,:].*
        #                         conj.(im_rd_0_1[j2,l2,:,:])
        #                     )
        #                 )
        #             end
        #         end
        #     end
        # end

        append!(out_coeff,S20)
    end

    ##Coeff at (1,2) level
    if coeff_12_on == 1
        S12 = Array{Float64,4}(undef, J, J, L, L)
        for j1 = 1:J
            for l1 = 1:L
                not_rot_im = fftshift(im_fd_0_1[j1,l1,:,:])
                for l2 = 1:L
                    for j2 = 1:J
                        S12[j1,j2,l1,l2] = sum(
                            abs.(
                                not_rot_im.*
                                conj(fftshift(im_fd_0_1[j2,l2,:,:]))
                            )
                        )
                    end
                end
            end
        end
    append!(out_coeff,S12)
    end
    return out_coeff
end

temp = repeat(fink_filter_set[1,3,:,:],inner=[1,1,8,8])

@benchmark had!(fink_filter_set[1,3,:,:],fink_filter_set[1,3,:,:])

function wraper(A,B)
    A.*B
end

@benchmark wraper(fink_filter_set[1,3,:,:],fink_filter_set[1,3,:,:])

## Slightly better inplace S1

@benchmark DHC(fink_filter_set[1,3,:,:],fink_filter_set,coeff_12_on =0, coeff_20_on = 0)

## Start checking the IFFT step
out1 = DHC(fink_filter_set[1,3,:,:],fink_filter_set,coeff_12_on =0, coeff_20_on = 1)

@benchmark DHC(fink_filter_set[1,3,:,:],fink_filter_set,coeff_12_on =0, coeff_20_on = 1)

## Quick FFT

using LinearAlgebra

x = rand(Complex{Float64}, 256,256)
FFTW.set_num_threads(1)
p1 = plan_fft(x, flags=FFTW.MEASURE)
FFTW.set_num_threads(2)
p2 = plan_fft(x, flags=FFTW.MEASURE)
y = similar(x)

@time mul!(y, p1, x)

@time mul!(y, p2, x)

## FFT plan
function DHC(image, filter_set; norm_on = 1, coeff_12_on =1, coeff_20_on = 1)
    (Nx, Ny) = size(image)

    (J,L,_,_) = size(fink_filter_set)

    out_coeff = []

    ## Coeff at (1,0)
    #takes care of image normalization
    S0 = zeros(2)
    if norm_on == 1
        S0[1] = mean(image)
        norm_im = image.-S0[1]
        S0[2] = sum(abs2.(norm_im))/(Nx*Ny)
        norm_im = norm_im./sqrt(Nx*Ny*S0[2])
    else
        norm_im = image
    end

    append!(out_coeff,S0)

    ##Coeff at (1,1)
    #store the fft of images and im*filters for later
    im_fd_0 = fft(norm_im)

    S1 = zeros((J,L))
    im_fd_0_1 = repeat(im_fd_0,inner=[1,1,J,L])

    for j = 1:J
        for l = 1:L
            had!(im_fd_0_1[:,:,j,l],filter_set[j,l,:,:]) #wavelet already in fft domain not shifted
            S1[j,l]=sum(abs2.(im_fd_0_1[:,:,j,l])) #normalization choice arb to make order unity
        end
    end
    append!(out_coeff,S1)

    ##Coeff at (2,0)
    if coeff_20_on == 1
        x = rand(Complex{Float64}, 256,256)
        FFTW.set_num_threads(2)
        p2 = plan_ifft(x, flags=FFTW.MEASURE)

        S20 = Array{Float64,4}(undef, J, J, L, L)

        im_rd_0_1 = Array{ComplexF64,4}(undef, J, L, Nx, Ny)
        for j1 = 1:J
            for l1 = 1:L
                im_rd_0_1[j1,l1,:,:] = mul!(im_rd_0_1[j1,l1,:,:], p2, im_fd_0_1[:,:,j1,l1])
            end
        end

        # for j1 = 1:J
        #     for l1 = 1:L
        #         for l2 = 1:L
        #             for j2  = 1:J
        #                 S20[j1,j2,l1,l2] = sum(
        #                     abs.(
        #                         im_rd_0_1[j1,l1,:,:].*
        #                         conj.(im_rd_0_1[j2,l2,:,:])
        #                     )
        #                 )
        #             end
        #         end
        #     end
        # end

        append!(out_coeff,S20)
    end

    ##Coeff at (1,2) level
    if coeff_12_on == 1
        S12 = Array{Float64,4}(undef, J, J, L, L)
        for j1 = 1:J
            for l1 = 1:L
                not_rot_im = fftshift(im_fd_0_1[j1,l1,:,:])
                for l2 = 1:L
                    for j2 = 1:J
                        S12[j1,j2,l1,l2] = sum(
                            abs.(
                                not_rot_im.*
                                conj(fftshift(im_fd_0_1[j2,l2,:,:]))
                            )
                        )
                    end
                end
            end
        end
    append!(out_coeff,S12)
    end
    return out_coeff
end

## checking the IFFT step

out2 = DHC(fink_filter_set[1,3,:,:],fink_filter_set,coeff_12_on =0, coeff_20_on = 1)

@benchmark DHC(fink_filter_set[1,3,:,:],fink_filter_set,coeff_12_on =0, coeff_20_on = 1)

size(fink_filter_set)

A = out2-out1

minimum(out2)
maximum(out2)
minimum(out1)
maximum(out1)

## Ok, something is wrong with the fft preallocation. Seems like we are getting overflows.
# out 1 also seems to have some big numbers though... since I got a factor of 10 speed up using 2 cores..
# seems paralleizing the native version might be the way to go

function DHC(image, filter_set; norm_on = 1, coeff_12_on =1, coeff_20_on = 1)
    FFTW.set_num_threads(2)
    (Nx, Ny) = size(image)

    (J,L,_,_) = size(fink_filter_set)

    out_coeff = []

    ## Coeff at (1,0)
    #takes care of image normalization
    S0 = zeros(2)
    if norm_on == 1
        S0[1] = mean(image)
        norm_im = image.-S0[1]
        S0[2] = sum(abs2.(norm_im))/(Nx*Ny)
        norm_im = norm_im./sqrt(Nx*Ny*S0[2])
    else
        norm_im = image
    end

    append!(out_coeff,S0)

    ##Coeff at (1,1)
    #store the fft of images and im*filters for later
    im_fd_0 = fft(norm_im)

    S1 = zeros((J,L))
    im_fd_0_1 = repeat(im_fd_0,inner=[1,1,J,L])
    if coeff_20_on == 1
        im_rd_0_1 = Array{ComplexF64,4}(undef, J, L, Nx, Ny)
    end

    @views for j = 1:J
        for l = 1:L
            had!(im_fd_0_1[:,:,j,l],filter_set[j,l,:,:]) #wavelet already in fft domain not shifted
            S1[j,l]+=sum(abs2.(im_fd_0_1[:,:,j,l])) #normalization choice arb to make order unity
            if coeff_20_on == 1
                im_rd_0_1[j1,l1,:,:] = ifft(im_fd_0_1[:,:,j1,l1])
            end
        end
    end
    append!(out_coeff,S1)

    ##Coeff at (2,0)
    if coeff_20_on == 1
        S20 = Array{Float64,4}(undef, J, J, L, L)

        # for j1 = 1:J
        #     for l1 = 1:L
        #         for l2 = 1:L
        #             for j2  = 1:J
        #                 S20[j1,j2,l1,l2] = sum(
        #                     abs.(
        #                         im_rd_0_1[j1,l1,:,:].*
        #                         conj.(im_rd_0_1[j2,l2,:,:])
        #                     )
        #                 )
        #             end
        #         end
        #     end
        # end

        append!(out_coeff,S20)
    end

    ##Coeff at (1,2) level
    if coeff_12_on == 1
        S12 = Array{Float64,4}(undef, J, J, L, L)
        for j1 = 1:J
            for l1 = 1:L
                not_rot_im = fftshift(im_fd_0_1[j1,l1,:,:])
                for l2 = 1:L
                    for j2 = 1:J
                        S12[j1,j2,l1,l2] = sum(
                            abs.(
                                not_rot_im.*
                                conj(fftshift(im_fd_0_1[j2,l2,:,:]))
                            )
                        )
                    end
                end
            end
        end
    append!(out_coeff,S12)
    end
    return out_coeff
end

## 2 core native FFT

out = DHC(fink_filter_set[1,3,:,:],fink_filter_set,coeff_12_on =0, coeff_20_on = 1)

@benchmark DHC(fink_filter_set[1,3,:,:],fink_filter_set,coeff_12_on =0, coeff_20_on = 1)

function DHC(image, filter_set; norm_on = 1, coeff_12_on =1, coeff_20_on = 1)
    (Nx, Ny) = size(image)

    (J,L,_,_) = size(fink_filter_set)

    out_coeff = []

    ## Coeff at (1,0)
    #takes care of image normalization
    S0 = zeros(2)
    if norm_on == 1
        S0[1] = mean(image)
        norm_im = image.-S0[1]
        S0[2] = sum(abs2.(norm_im))/(Nx*Ny)
        norm_im = norm_im./sqrt(Nx*Ny*S0[2])
    else
        norm_im = image
    end

    append!(out_coeff,S0)

    ##Coeff at (1,1)
    #store the fft of images and im*filters for later
    im_fd_0 = fft(norm_im)

    S1 = zeros((J,L))
    im_fd_0_1 = repeat(im_fd_0,inner=[1,1,J,L])

    for j = 1:J
        for l = 1:L
            had!(im_fd_0_1[:,:,j,l],filter_set[j,l,:,:]) #wavelet already in fft domain not shifted
            S1[j,l]=sum(abs2.(im_fd_0_1[:,:,j,l])) #normalization choice arb to make order unity
        end
    end
    append!(out_coeff,S1)

    ##Coeff at (2,0)
    if coeff_20_on == 1
        S20 = Array{Float64,4}(undef, J, J, L, L)

        FFTW.set_num_threads(1)
        im_rd_0_1 = Array{ComplexF64,4}(undef, J, L, Nx, Ny)
        for j1 = 1:J
            for l1 = 1:L
                im_rd_0_1[j1,l1,:,:] = ifft(im_fd_0_1[:,:,j1,l1])
            end
        end

        # for j1 = 1:J
        #     for l1 = 1:L
        #         for l2 = 1:L
        #             for j2  = 1:J
        #                 S20[j1,j2,l1,l2] = sum(
        #                     abs.(
        #                         im_rd_0_1[j1,l1,:,:].*
        #                         conj.(im_rd_0_1[j2,l2,:,:])
        #                     )
        #                 )
        #             end
        #         end
        #     end
        # end

        append!(out_coeff,S20)
    end

    ##Coeff at (1,2) level
    if coeff_12_on == 1
        S12 = Array{Float64,4}(undef, J, J, L, L)
        for j1 = 1:J
            for l1 = 1:L
                not_rot_im = fftshift(im_fd_0_1[j1,l1,:,:])
                for l2 = 1:L
                    for j2 = 1:J
                        S12[j1,j2,l1,l2] = sum(
                            abs.(
                                not_rot_im.*
                                conj(fftshift(im_fd_0_1[j2,l2,:,:]))
                            )
                        )
                    end
                end
            end
        end
    append!(out_coeff,S12)
    end
    return out_coeff
end

out = DHC(fink_filter_set[1,3,:,:],fink_filter_set,coeff_12_on =0, coeff_20_on = 1)

@benchmark DHC(fink_filter_set[1,3,:,:],fink_filter_set,coeff_12_on =0, coeff_20_on = 1)

FFTW.set_num_threads(1)
@time ifft(fink_filter_set[1,3,:,:])
FFTW.set_num_threads(2)
@time ifft(fink_filter_set[1,3,:,:])

##
using Profile

DHC(fink_filter_set[1,3,:,:],fink_filter_set,coeff_12_on =0, coeff_20_on = 1)
@time DHC(fink_filter_set[1,3,:,:],fink_filter_set,coeff_12_on =0, coeff_20_on = 1)

Profile.clear()
@profile DHC(fink_filter_set[1,3,:,:],fink_filter_set,coeff_12_on =0, coeff_20_on = 1)

Juno.profiler()

##
using Profile

function DHC(image, filter_set; norm_on = 1, coeff_12_on =1, coeff_20_on = 1)
    FFTW.set_num_threads(2)
    (Nx, Ny) = size(image)

    (J,L,_,_) = size(fink_filter_set)

    out_coeff = []

    ## Coeff at (1,0)
    #takes care of image normalization
    S0 = zeros(2)
    if norm_on == 1
        S0[1] = mean(image)
        norm_im = image.-S0[1]
        S0[2] = sum(abs2.(norm_im))/(Nx*Ny)
        norm_im = norm_im./sqrt(Nx*Ny*S0[2])
    else
        norm_im = image
    end

    append!(out_coeff,S0)

    ##Coeff at (1,1)
    #store the fft of images and im*filters for later
    im_fd_0 = fft(norm_im)

    S1 = zeros((J,L))
    im_fd_0_1 = repeat(im_fd_0,inner=[1,1,J,L])
    if coeff_20_on == 1
        im_rd_0_1 = Array{ComplexF64,4}(undef, J, L, Nx, Ny)
    end

    @views for j = 1:J
        for l = 1:L
            had!(im_fd_0_1[:,:,j,l],filter_set[j,l,:,:]) #wavelet already in fft domain not shifted
            S1[j,l]+=sum(abs2.(im_fd_0_1[:,:,j,l])) #normalization choice arb to make order unity
            if coeff_20_on == 1
                im_rd_0_1[j,l,:,:] = ifft(im_fd_0_1[:,:,j,l])
            end
        end
    end
    append!(out_coeff,S1)

    ##Coeff at (2,0)
    if coeff_20_on == 1
        S20 = Array{Float64,4}(undef, J, J, L, L)

        # for j1 = 1:J
        #     for l1 = 1:L
        #         for l2 = 1:L
        #             for j2  = 1:J
        #                 S20[j1,j2,l1,l2] = sum(
        #                     abs.(
        #                         im_rd_0_1[j1,l1,:,:].*
        #                         conj.(im_rd_0_1[j2,l2,:,:])
        #                     )
        #                 )
        #             end
        #         end
        #     end
        # end

        append!(out_coeff,S20)
    end

    ##Coeff at (1,2) level
    if coeff_12_on == 1
        S12 = Array{Float64,4}(undef, J, J, L, L)
        for j1 = 1:J
            for l1 = 1:L
                not_rot_im = fftshift(im_fd_0_1[j1,l1,:,:])
                for l2 = 1:L
                    for j2 = 1:J
                        S12[j1,j2,l1,l2] = sum(
                            abs.(
                                not_rot_im.*
                                conj(fftshift(im_fd_0_1[j2,l2,:,:]))
                            )
                        )
                    end
                end
            end
        end
    append!(out_coeff,S12)
    end
    return out_coeff
end

DHC(fink_filter_set[1,3,:,:],fink_filter_set,coeff_12_on =0, coeff_20_on = 1)
@time DHC(fink_filter_set[1,3,:,:],fink_filter_set,coeff_12_on =0, coeff_20_on = 1)

Profile.clear()
@profile DHC(fink_filter_set[1,3,:,:],fink_filter_set,coeff_12_on =0, coeff_20_on = 1)

Juno.profiler()

@benchmark DHC(fink_filter_set[1,3,:,:],fink_filter_set,coeff_12_on =0, coeff_20_on = 1)

##
function DHC(image, filter_set; norm_on = 1, coeff_12_on =1, coeff_20_on = 1)
    FFTW.set_num_threads(2)
    (Nx, Ny) = size(image)

    (J,L,_,_) = size(fink_filter_set)

    out_coeff = []

    ## Coeff at (1,0)
    #takes care of image normalization
    S0 = zeros(2)
    if norm_on == 1
        S0[1] = mean(image)
        norm_im = image.-S0[1]
        S0[2] = sum(abs2.(norm_im))/(Nx*Ny)
        norm_im = norm_im./sqrt(Nx*Ny*S0[2])
    else
        norm_im = image
    end

    append!(out_coeff,S0)

    ##Coeff at (1,1)
    #store the fft of images and im*filters for later
    im_fd_0 = fft(norm_im)

    S1 = zeros((J,L))
    im_fd_0_1 = repeat(im_fd_0,inner=[1,1,J,L])
    if coeff_20_on == 1
        im_rd_0_1 = Array{ComplexF64,4}(undef, J, L, Nx, Ny)
    end

    p2 = plan_ifft(im_fd_0_1[:,:,1,1], flags=FFTW.MEASURE)

    for j = 1:J
        for l = 1:L
            had!(im_fd_0_1[:,:,j,l],filter_set[j,l,:,:]) #wavelet already in fft domain not shifted
            S1[j,l]+=sum(abs2.(im_fd_0_1[:,:,j,l])) #normalization choice arb to make order unity
            if coeff_20_on == 1
                mul!(im_rd_0_1[j,l,:,:], p2, im_fd_0_1[:,:,j,l])
            end
        end
    end
    append!(out_coeff,S1)

    ##Coeff at (2,0)
    if coeff_20_on == 1
        S20 = Array{Float64,4}(undef, J, J, L, L)

        # for j1 = 1:J
        #     for l1 = 1:L
        #         for l2 = 1:L
        #             for j2  = 1:J
        #                 S20[j1,j2,l1,l2] = sum(
        #                     abs.(
        #                         im_rd_0_1[j1,l1,:,:].*
        #                         conj.(im_rd_0_1[j2,l2,:,:])
        #                     )
        #                 )
        #             end
        #         end
        #     end
        # end

        append!(out_coeff,S20)
    end

    ##Coeff at (1,2) level
    if coeff_12_on == 1
        S12 = Array{Float64,4}(undef, J, J, L, L)
        for j1 = 1:J
            for l1 = 1:L
                not_rot_im .= fftshift(im_fd_0_1[j1,l1,:,:])
                for l2 = 1:L
                    for j2 = 1:J
                        S12[j1,j2,l1,l2] = sum(
                            abs.(
                                not_rot_im.*
                                conj(fftshift(im_fd_0_1[j2,l2,:,:]))
                            )
                        )
                    end
                end
            end
        end
    append!(out_coeff,S12)
    end
    return out_coeff
end

DHC(fink_filter_set[1,3,:,:],fink_filter_set,coeff_12_on =0, coeff_20_on = 1)
@time DHC(fink_filter_set[1,3,:,:],fink_filter_set,coeff_12_on =0, coeff_20_on = 1)

Profile.clear()
@profile DHC(fink_filter_set[1,3,:,:],fink_filter_set,coeff_12_on =0, coeff_20_on = 1)

Juno.profiler()

@benchmark DHC(fink_filter_set[1,3,:,:],fink_filter_set,coeff_12_on =0, coeff_20_on = 1)

##

function DHC(image, filter_set; norm_on = 1, coeff_12_on =1, coeff_20_on = 1)
    FFTW.set_num_threads(2)
    (Nx, Ny) = size(image)

    (J,L,_,_) = size(fink_filter_set)

    out_coeff = []

    ## Coeff at (1,0)
    #takes care of image normalization
    S0 = zeros(2)
    if norm_on == 1
        S0[1] = mean(image)
        norm_im = image.-S0[1]
        S0[2] = sum(abs2.(norm_im))/(Nx*Ny)
        norm_im = norm_im./sqrt(Nx*Ny*S0[2])
    else
        norm_im = image
    end

    append!(out_coeff,S0)

    ##Coeff at (1,1)
    #store the fft of images and im*filters for later
    im_fd_0 = fft(norm_im)

    S1 = zeros((J,L))
    im_fd_0_1 = Array{ComplexF64,4}(undef, J, L, Nx, Ny)
    @views for j = 1:J
        for l = 1:L
            im_fd_0_1[j,l,:,:] .= im_fd_0
        end
    end

    if coeff_20_on == 1
        im_rd_0_1 = Array{ComplexF64,4}(undef, J, L, Nx, Ny)
    end

    @views for j = 1:J
        for l = 1:L
            had!(im_fd_0_1[j,l,:,:],filter_set[j,l,:,:]) #wavelet already in fft domain not shifted
            S1[j,l]+=sum(abs2.(im_fd_0_1[j,l,:,:])) #normalization choice arb to make order unity
            if coeff_20_on == 1
                im_rd_0_1[j,l,:,:] .= ifft(im_fd_0_1[j,l,:,:])
            end
        end
    end
    append!(out_coeff,S1)

    ##Coeff at (2,0)
    if coeff_20_on == 1
        S20 = Array{Float64,4}(undef, J, J, L, L)

        # for j1 = 1:J
        #     for l1 = 1:L
        #         for l2 = 1:L
        #             for j2  = 1:J
        #                 S20[j1,j2,l1,l2] = sum(
        #                     abs.(
        #                         im_rd_0_1[j1,l1,:,:].*
        #                         conj.(im_rd_0_1[j2,l2,:,:])
        #                     )
        #                 )
        #             end
        #         end
        #     end
        # end

        append!(out_coeff,S20)
    end

    ##Coeff at (1,2) level
    if coeff_12_on == 1
        S12 = Array{Float64,4}(undef, J, J, L, L)
        for j1 = 1:J
            for l1 = 1:L
                not_rot_im = fftshift(im_fd_0_1[j1,l1,:,:])
                for l2 = 1:L
                    for j2 = 1:J
                        S12[j1,j2,l1,l2] = sum(
                            abs.(
                                not_rot_im.*
                                conj(fftshift(im_fd_0_1[j2,l2,:,:]))
                            )
                        )
                    end
                end
            end
        end
    append!(out_coeff,S12)
    end
    return out_coeff
end

DHC(fink_filter_set[1,3,:,:],fink_filter_set,coeff_12_on =0, coeff_20_on = 1)
@time DHC(fink_filter_set[1,3,:,:],fink_filter_set,coeff_12_on =0, coeff_20_on = 1)

Profile.clear()
@profile DHC(fink_filter_set[1,3,:,:],fink_filter_set,coeff_12_on =0, coeff_20_on = 1)

Juno.profiler()

@benchmark DHC(fink_filter_set[1,3,:,:],fink_filter_set,coeff_12_on =0, coeff_20_on = 1)

## testing my own repeat
J=8
L=8

test_img = fink_filter_set[1,3,:,:]

test = Array{ComplexF64,4}(undef, J, L, 256, 256)
@views for j = 1:J
    for l = 1:L
        test[j,l,:,:] = test_img
    end
end

max_stars_repo_name: andrew-saydjari/DHC

max_issues_repo_name: andrew-saydjari/DHC

max_forks_repo_name: andrew-saydjari/DHC