text:
using PyCall
unshift!(PyVector(pyimport("sys")["path"]), "");
@pyimport Compute_Jacobian_Sioux

numNodes = Compute_Jacobian_Sioux.numNodes;
numLinks = Compute_Jacobian_Sioux.numLinks;
numZones = Compute_Jacobian_Sioux.numZones;
numODpairs = Compute_Jacobian_Sioux.numODpairs;
od_pairs = Compute_Jacobian_Sioux.od_pairs;
link_list_js = Compute_Jacobian_Sioux.link_list_js;
link_length_list = Compute_Jacobian_Sioux.link_length_list;

require("prepare_data_Sioux.jl");
require("extract_data_Sioux.jl");
require("../Julia_files/tap_MSA.jl");
require("../Julia_files/inverseVI.jl");
require("../Julia_files/demands_adjustment_gradi.jl");

require("initia_data.jl");

demandsDiffDict[1] = norm(demandsDicToVec(demandsDict[1]) - demandsDicToVec(demandsDict[0]))/
                     norm(demandsDicToVec(demandsDict[0]));

fcoeffs = [1, 0, 0, 0, .15];

# get observed flow vector (corresponding to ground truth demands and ground truth costs)
tapFlowDicDict[0], tapFlowVecDict[0] = tapMSA(demandsDict[0], fcoeffs);

# get arcs data corresponding to ground truth demands and flows
arcsDict[0] = observFlow("./benchmark_data/Sioux_net.txt", tapFlowDicDict[0]);

deg = 6
c = 3.41
lam = 1.

coeffs_dict[(deg, c, lam, 1)], objInvVIDict[1] = train(lam, deg, c, demandsDict[1], arcsDict[0]);

# get initial flow vector (corresponding to initial demands)
tapFlowDicDict[1], tapFlowVecDict[1] = tapMSA(demandsDict[1], coeffs_dict[(deg, c, lam, 1)]);

demandsVecDict[0] = demandsDicToVec(demandsDict[0]);

objFunDict[1] = objF(demandsVecDict[1], coeffs_dict[(deg, c, lam, 1)]);

# get arcs data corresponding to initial demands and flows
arcsDict[1] = observFlow("./benchmark_data/Sioux_net.txt", tapFlowDicDict[1]);

linkCostDicDict[1] = tapFlowVecToLinkCostDict(tapFlowVecDict[1], coeffs_dict[(deg, c, lam, 1)]);

linkCostDicDict[1]["0"], link_length_list[1]

jacobiSpiessDict[1] = Compute_Jacobian_Sioux.jacobianSpiess(numNodes, numLinks, numODpairs, od_pairs,
                                              link_list_js, [linkCostDicDict[1]["$(i)"] for i=0:numLinks-1]);

# maximum number of iterations
N = 5;

# Armijo rule parameters
rho = 2;
M = 10;

# search direction parameter
epsilon_1 = 0;

# stop criterion parameter
epsilon_2 = 1e-20;

for l = 1:N

    jacobDict[l] = jacobiSpiessDict[l]

    descDirecDict[l] = descDirec(tapFlowVecDict[l], tapFlowVecDict[0], jacobDict[l], numODpairs, numLinks);

    demandsVecDict[l] = demandsDicToVec(demandsDict[l]);

    searchDirecDict[l] = searchDirec(demandsVecDict[l], descDirecDict[l], epsilon_1);

    thetaMaxDict[l] = thetaMax(demandsVecDict[l], searchDirecDict[l]);

    demandsVecDict[l+1] = similar(demandsVecDict[0]);

    demandsVecDict[l+1], objFunDict[l+1] = armijo(objFunDict[l], demandsVecDict[l], coeffs_dict[(deg, c, lam, l)],
    searchDirecDict[l], thetaMaxDict[l], rho, M);

    demandsDict[l+1] = demandsVecToDic(demandsVecDict[l+1]);

    tapFlowDicDict[l+1], tapFlowVecDict[l+1] = tapMSA(demandsDict[l+1], coeffs_dict[(deg, c, lam, l)]);

    arcsDict[l+1] = observFlow("./benchmark_data/Sioux_net.txt", tapFlowDicDict[l+1]);

    coeffs_dict[(deg, c, lam, l+1)], objInvVIDict[l+1] = train(lam, deg, c, demandsDict[l+1], arcsDict[0]);

    linkCostDicDict[l+1] = tapFlowVecToLinkCostDict(tapFlowVecDict[l+1], coeffs_dict[(deg, c, lam, l+1)]);

    jacobiSpiessDict[l+1] = Compute_Jacobian_Sioux.jacobianSpiess(numNodes, numLinks, numODpairs, od_pairs,
                                              link_list_js, [linkCostDicDict[l+1]["$(i)"] for i=0:numLinks-1]);

    demandsDiffDict[l+1] = norm(demandsVecDict[l+1] - demandsVecDict[0]) / norm(demandsVecDict[0]);

    # stopping criterion
    if (objFunDict[l] - objFunDict[l+1]) / objFunDict[1] < epsilon_2
        break
    end

    println("$(l) iterations finished...")
    println("current objective function value is $(objFunDict[l+1])")
    println("previous objective function value is $(objFunDict[l])")
    println("----------------------------")
end

# normalize objective function value
for l = 1:(length(objFunDict))
    norObjFunDict[l] = objFunDict[l] / objFunDict[1];
end

#demandsVecDict

#norObjFunDict

using PyPlot

iterNum = 1:(length(norObjFunDict))

objFun = map(iterNum->norObjFunDict[iterNum], iterNum)

# plot(iterNum, objFun, "s-g", label="True")

plot(iterNum, objFun, "s-b")

# legend(loc="upper right",fancybox="true")

grid("on")

xlim(.8, length(norObjFunDict) + 1 + .2)
ylim(minimum(objFun)-.01, maximum(objFun)+.01)

font1 = ["family"=>"serif", "color"=>"darkred", "weight"=>"normal", "size"=>14]

xlabel("Number of iterations", fontdict=font1)
ylabel("Normalized objective function value (biLev)", fontdict=font1)

savefig("./results/objFun_biLev_Sioux.pdf")

using PyPlot

iterNum = 1:(length(demandsDiffDict))

demandsDiff = map(iterNum->demandsDiffDict[iterNum], iterNum)

# plot(iterNum, objFun, "s-g", label="True")

plot(iterNum, demandsDiff, "o-g")

# legend(loc="upper right",fancybox="true")

grid("on")

xlim(.8, length(demandsDiffDict) + 1 + .2)
ylim(minimum(demandsDiff)-.0005, maximum(demandsDiff)+.0005)

font1 = ["family"=>"serif", "color"=>"darkred", "weight"=>"normal", "size"=>14]

xlabel("Number of iterations", fontdict=font1)
ylabel("Normalized demand difference", fontdict=font1)

savefig("./results/demandsDiff_biLev_Sioux.pdf")

using PyPlot

iterNum = 1:(length(objInvVIDict))

objFun = map(iterNum->objInvVIDict[iterNum], iterNum)

# plot(iterNum, objFun, "s-g", label="True")

plot(iterNum, objFun, "p-r")

# legend(loc="upper right",fancybox="true")

grid("on")

xlim(0.8, length(objInvVIDict) + 1 + .2)
# ylim(0.9, 2.0)

font1 = ["family"=>"serif", "color"=>"darkred", "weight"=>"normal", "size"=>14]

xlabel("Number of iterations", fontdict=font1)
ylabel("Objective function value (invVI)", fontdict=font1)

savefig("./results/objFun_invVI_Sioux.pdf")

#coeffs_dict

using PyPlot

iterN = 1

true_coeffs = [1, 0, 0, 0, .15]
# true_coeffs = [1, .2, .5, .8, .15]

est_coeffs = coeffs_dict[(6, 3.41, 1., iterN)]

xs = linspace(0, 2, 20)
zs_true = map(x->polyEval(true_coeffs, x), xs)

zs = map(x->polyEval(est_coeffs, x), xs)

plot(xs, zs_true, "s-g", label="True")

plot(xs, zs, "^-m", label="deg=6")
legend(loc="upper left",fancybox="true")

grid("on")
xlim(-0.1, 1.6);
ylim(0.9, 2.0);

font1 = ["family"=>"serif","color"=>"darkred","weight"=>"normal","size"=>14]
xlabel("Scaled Flow", fontdict=font1)
ylabel("Scaled Cost", fontdict=font1)

savefig("./results/fitting_Sioux_iter_$(iterN).pdf")

using PyPlot

iterN = 2

true_coeffs = [1, 0, 0, 0, .15]
# true_coeffs = [1, .2, .5, .8, .15]

est_coeffs = coeffs_dict[(6, 3.41, 1., iterN)]

xs = linspace(0, 2, 20)
zs_true = map(x->polyEval(true_coeffs, x), xs)

zs = map(x->polyEval(est_coeffs, x), xs)

plot(xs, zs_true, "s-g", label="True")

plot(xs, zs, "^-m", label="deg=6")
legend(loc="upper left",fancybox="true")

grid("on")
xlim(-0.1, 1.6);
ylim(0.9, 2.0);

font1 = ["family"=>"serif","color"=>"darkred","weight"=>"normal","size"=>14]
xlabel("Scaled Flow", fontdict=font1)
ylabel("Scaled Cost", fontdict=font1)

savefig("./results/fitting_Sioux_iter_$(iterN).pdf")

using PyPlot

iterN = 3

true_coeffs = [1, 0, 0, 0, .15]
# true_coeffs = [1, .2, .5, .8, .15]

est_coeffs = coeffs_dict[(6, 3.41, 1., iterN)]

xs = linspace(0, 2, 20)
zs_true = map(x->polyEval(true_coeffs, x), xs)

zs = map(x->polyEval(est_coeffs, x), xs)

plot(xs, zs_true, "s-g", label="True")

plot(xs, zs, "^-m", label="deg=6")
legend(loc="upper left",fancybox="true")

grid("on")
xlim(-0.1, 1.6);
ylim(0.9, 2.0);

font1 = ["family"=>"serif","color"=>"darkred","weight"=>"normal","size"=>14]
xlabel("Scaled Flow", fontdict=font1)
ylabel("Scaled Cost", fontdict=font1)

savefig("./results/fitting_Sioux_iter_$(iterN).pdf")

using PyPlot

iterN = length(coeffs_dict)

true_coeffs = [1, 0, 0, 0, .15]
# true_coeffs = [1, .2, .5, .8, .15]

est_coeffs = coeffs_dict[(6, 3.41, 1., iterN)]

xs = linspace(0, 2, 20)
zs_true = map(x->polyEval(true_coeffs, x), xs)

zs = map(x->polyEval(est_coeffs, x), xs)

plot(xs, zs_true, "s-g", label="True")

plot(xs, zs, "^-m", label="deg=6")
legend(loc="upper left",fancybox="true")

grid("on")
xlim(-0.1, 1.6);
ylim(0.9, 2.0);

font1 = ["family"=>"serif","color"=>"darkred","weight"=>"normal","size"=>14]
xlabel("Scaled Flow", fontdict=font1)
ylabel("Scaled Cost", fontdict=font1)

savefig("./results/fitting_Sioux_iter_$(iterN).pdf")

outfile = open("./results/norObjFunDict.json", "w")

JSON.print(outfile, norObjFunDict)

close(outfile)

outfile = open("./results/demandsDiffDict.json", "w")

JSON.print(outfile, demandsDiffDict)

close(outfile)

outfile = open("./results/objInvVIDict.json", "w")

JSON.print(outfile, objInvVIDict)

close(outfile)

outfile = open("./results/coeffs_dict.json", "w")

JSON.print(outfile, coeffs_dict)

close(outfile)

outfile = open("./results/demandsVecDict.json", "w")

JSON.print(outfile, demandsVecDict)

close(outfile)

outfile = open("./results/demandsDict.json", "w")

JSON.print(outfile, demandsDict)

close(outfile)

outfile = open("./results/tapFlowDicDict.json", "w")

JSON.print(outfile, tapFlowDicDict)

close(outfile)

outfile = open("./results/tapFlowVecDict.json", "w")

JSON.print(outfile, tapFlowVecDict)

close(outfile)

max_stars_repo_name: jingzbu/InverseVITraffic

max_issues_repo_name: jingzbu/InverseVITraffic

max_forks_repo_name: jingzbu/InverseVITraffic


text:
# Standard sigmoid layer: f(A,x,b)=1/(1+exp(-(A*x+b)))
# Note that the bias (a column vector) gets expanded here to match the dimension of A*x

# To make the CPU better performing, I think one would need to call c code that does inplace computation of Hadamard products. Essentially this is what is happending in the GPU code.

#FsigmoidAXplusBias(A,X,b)=(1./(1+exp(-(A*X+b*ones(1,size(X,2))))),nothing)

#TODO: don't need aux to be a tuple since we only use aux[1]
FsigmoidAXplusBias(A,X,b)=(1./(1+exp(-(A*X+b*ones(1,size(X,2))))),(zeros(size(A,1),size(X,2)),nothing)) # allocate memory for inplace gradients

FsigmoidAXplusBias_inplace(value,aux,A,X,b)=copy!(value,1./(1+exp(-(A*X+b*ones(1,size(X,2))))))

#TODO; Need to do similar inplace computations for the other transfer functions
function DsigmoidAXplusBias(derivativeIDX,f_c,faux_c,grad_c,grad_n,A,X,b)
    copy!(faux_c[1],f_c)
    A_elmult_B_update!(-1.0,f_c,f_c,1.0,faux_c[1])
    A_elmult_B_update!(1.0,grad_c,faux_c[1],0.0,faux_c[1])

    if derivativeIDX==1
        #axpy!(1.0,(grad_c.*f_c.*(1-f_c))*X',grad_n) # this bad since it creates creates temporary storage which will need to be collected
        BLAS.gemm!('N','T',1.0,faux_c[1],X,1.0,grad_n)
    elseif derivativeIDX==2
        #axpy!(1.0,A'*(grad_c.*f_c.*(1-f_c)),grad_n)
        BLAS.gemm!('T','N',1.0,A,faux_c[1],1.0,grad_n)
    elseif derivativeIDX==3
        axpy!(1.0,sum(grad_c.*f_c.*(1-f_c),2),grad_n)
    end
end

if PROC=="GPU"
    function FsigmoidAXplusBias(A::CudaArray,X::CudaArray,b::CudaArray)
        return(sigmoid(FAXplusBias(A,X,b)[1]),nothing)
    end

    function FsigmoidAXplusBias_inplace(value,aux,A::CudaArray,X::CudaArray,b::CudaArray)
#        println(size(A))
#        println(size(X))
#        println(size(b))
        FAXplusBias_inplace(value,aux,A,X,b)
        sigmoid!(value,value)
    end

    function DsigmoidAXplusBias(derivativeIDX,f_c,faux_c,grad_c,grad_n,A::CudaArray{Float64},X::CudaArray{Float64},b::CudaArray{Float64})
        tmp=CudaArray(Float64,size(grad_c)); fill!(tmp,0.0)
        tx1mx!(grad_c,f_c,tmp)
        if derivativeIDX==1
            gemm!('N','T',1.0,tmp,X,1.0,grad_n)
        elseif derivativeIDX==2
            gemm!('T','N',1.0,A,tmp,1.0,grad_n)
        elseif derivativeIDX==3
            ons=CudaArray(Float64,(size(X,2),1)); fill!(ons,1.0)
            gemm!('N','N',1.0,tmp,ons,1.0,grad_n)
            free(ons)
        end
        free(tmp)
    end

    function DsigmoidAXplusBias(derivativeIDX,f_c,faux_c,grad_c,grad_n,A::CudaArray{Float32},X::CudaArray{Float32},b::CudaArray{Float32})
        tmp=CudaArray(Float32,size(grad_c)); fill!(tmp,0.0)
        tx1mx!(grad_c,f_c,tmp)
        if derivativeIDX==1
            gemm!('N','T',1.0,tmp,X,1.0,grad_n)
        elseif derivativeIDX==2
            gemm!('T','N',1.0,A,tmp,1.0,grad_n)
        elseif derivativeIDX==3
            ons=CudaArray(Float32,(size(X,2),1)); fill!(ons,1.0)
            gemm!('N','N',1.0,tmp,ons,1.0,grad_n)
            free(ons)
        end
        free(tmp)
    end

end

Derivative[FsigmoidAXplusBias]=DsigmoidAXplusBias
Inplace[FsigmoidAXplusBias]=FsigmoidAXplusBias_inplace

sigmoidAXplusBias(A,X,b)=ADnode(FsigmoidAXplusBias,[A X b])
export sigmoidAXplusBias

ext: jl

max_stars_repo_name: davidbarber/AutoDiff

max_issues_repo_name: davidbarber/AutoDiff

max_forks_repo_name: davidbarber/AutoDiff


text:
module TrackError

using ..Util, ..Tagging, ..Condition, ..Traits
import ..Condition
export applytrackerr, condvar, tagerror

# Issue is that type of err value depends on
# gradient method used
# Solutions
# Specialize: condvar/applytrackerr based on type of Omega

"""
`applytrackerr(x, ω, initerr)`

Is `ω` in the domain of `x`?
Returns `(f(ω), v)` where `v = ω ∈ dom(x)`
"""
function applytrackerr(x, ω, initerr)
  ω_ = tagerror(ω, initerr) # zt: surprisingly expensive
  fx = x(ω_)
  (fx = fx, err = ω_.tags.err.val)
end

# tagnotrackerr(ω) = tag(ω, (donttrack = true,))
# applynotrackerr(x, ω) = apl(x, tagnotrackerr(ω))

# "Soft `condvar`: distance from `ω` to the domain of `x`"
# condvarₛ(x, ω, initerr = softtrue()) = applytrackerr(x, ω, initerr).err
# condvarₛ(x) = ω -> condvarₛ(x, ω)

"Is `ω` in the domain of `x`?"
condvar(x, ω, initerr = true) = applytrackerr(x, ω, initerr).err
condvar(x) = ω -> condvar(x, ω)

"Update `err` by conjoining cu`rrent error with `b`"
conjoinerror!(err::Box, b) = err.val &= b

dotrack(ω) = !haskey(ω.tags, :donttrack) && haskey(ω.tags, :err)

function Condition.condf(::trait(Err), ω, x, y)
  dotrack(ω) && conjoinerror!(ω.tags.err, y(ω))
  x(ω)
end

function Condition.cond!(::trait(Err), ω, bool)
  dotrack(ω) && conjoinerror!(ω.tags.err, bool)
  bool
end

"Tag `ω` with `err`"
tagerror(ω, initerr) = tag(ω, (err = Box(initerr),))

# tagerror(ω, initerr::AbstractBool) = tag(ω, (err = Ref(initerr),)) #FIXME UNCOMMENT THIS AND SPECIALISE BELOW TO OMEGA
# "Tag tω with error"
# function tagerror(tω::TaggedΩ, initerr::AbstractBool)
#   # If error already there, use that!
#   haskey(tω.tags, :err) ? tω : tag(tω, (err = Ref{Real}(initerr),))
# end

function combinetags(::Type{Val{:err}}, a, b)
  # @show a
  # @show b
  # @assert false
  @show a[]
  @show a[] & b[]
end

end


text:
"""
    module IJV

Algorithms for sparse matrices stored in COO form.

`IJV` does not depend on any types not in `Core`.
Do not import any symbol from `IJV` into another module. The methods
in `IJV` are incompatible with functions in `Base` with the same
(unqualified) name. Instead, use `IJV` like this

#Example

```julia
Base.getindex(S::MySparseMatrixCOO, i::Integer, j::Integer) = IJV.getindex(S.I, S.J, S.V, i, j)
```
"""
module IJV

import Random
import StatsBase
import SparseArrays
import SparseUtils: findrepeated

spzeros(m::Integer, n::Integer) = spzeros(Float64, m, n)
spzeros(::Type{Tv}, m::Integer, n::Integer) where {Tv} = spzeros(Tv, Int, m, n)
function spzeros(::Type{Tv}, ::Type{Ti}, m::Integer, n::Integer) where {Tv, Ti}
    ((m < 0) || (n < 0)) && throw(ArgumentError("invalid Array dimensions"))
    return (m, n, Ti[], Ti[], Tv[])
end

### sorting and indexing

"""
    issorted(I::AbstractVector, J::AbstractVector)

Return true if the COO indices `I` and `J` are sorted canonically.
"""
issorted(I::AbstractVector, J::AbstractVector) = Base.issorted(zip(J, I))

"""
    getindex(I, J, V, i::Integer, j::Integer; zeroel=zero(eltype(V)))

Return the element at coordinates `(i, j)` in the matrix represented by  `I, J, V`,
if it exists, and `zeroel` otherwise.
"""
function getindex(I, J, V, i::Integer, j::Integer; zeroel=zero(eltype(V)))
    index = findindex(I, J, i, j)
    index == nothing && return zeroel
    return V[index]
end

## This relies on efficient `Union{nothing, Int}`
"""
    findindex(I::AbstractArray, J::AbstractArray, i::Integer, j::Integer)

Return the common linear index into `I` and `J` (and `V`) representing coordinates
`(i, j)`, or `nothing` if there is no value stored for `(i, j)`.
"""
@inline function findindex(I::AbstractArray, J::AbstractArray, i::Integer, j::Integer)
    colrange = searchsorted(J, j)
    isempty(colrange) && return nothing
    rowrange = searchsorted(view(I, colrange), i)
    isempty(rowrange) && return nothing
    return first(rowrange) + first(colrange) - 1
end

function ijvinsert!(ind, arrays = (I, J, V), elements = (i, j, val))
    for (array, element) in zip(arrays, elements)
        insert!(array, ind, element)
    end
    return nothing
end

function ijvsetindex!(ind, arrays = (I, J, V), elements = (i, j, val))
    for (array, element) in zip(arrays, elements)
        array[ind] = element
    end
    return nothing
end

function setindex!(I, J, V, val, i::Integer, j::Integer; zeroel=zero(eltype(V)))
    rowrange = searchsorted(I, i)
    if isempty(rowrange)
        ind = first(rowrange)
        ijvinsert!(ind, (I, J, V), (i, j, val))
    else
        colrange = searchsorted(view(J, rowrange), j)
        if isempty(colrange)
            ijvinsert!(first(colrange), (I, J, V), (i, j, val))
        else
            valind = first(colrange) + first(rowrange) - 1
            ijvsetindex!(valind, (I, J, V), (i, j, val))
        end
    end
    return val
end

@inline sortpermindex(I, J) = sortperm(1:length(I), by = i -> @inbounds (J[i], I[i]))
@inline function ijvsort!(I, J, V)
    p = sortpermindex(I, J)
    permute!(I, p)
    permute!(J, p)
    permute!(V, p)
    return nothing
end

function permutedims!(m, n, I, J, V)
    ijvsort!(J, I, V) # swap I and J before sorting
    return (n, m, J, I, V)
end
permutedims(m, n, I, J, V) = permutedims!(m, n, copy(I), copy(J), copy(V))

function Array(m, n, I, J, V)
    a = zeros(eltype(V), m, n)
    for ind in 1:length(I)
        a[I[ind], J[ind]] = V[ind]
    end
    return a
end

### dropzeros, dropstored!

function dropzeros!(I, J, V::AbstractArray{T}) where T
    mask = V .== zero(T)
    map(x -> deleteat!(x, mask), (I, J, V))
end
dropzeros(I, J, V::AbstractArray) = dropzeros!(copy(I), copy(J), copy(V))

function dropstored!(I, J, V, i::Integer, j::Integer)
    index = findindex(I, J, i, j)
    index == nothing && return nothing
    deleteat!(I, index)
    deleteat!(J, index)
    deleteat!(V, index)
    return nothing
end

### prunecols

## We'd like `x[inds]`, but this appears to be a missing feature in Iterators.flatten
_applyinds(x, inds) = [x[i] for i in inds]

function _make_new_Inds(ranges)
    fillarr = [fill(i, length(ranges[i])) for i in 1:length(ranges)]
    return reduce(append!, fillarr)
end

shiftrange(range, newstart::Integer=1) = range .- (range.start - newstart)

function prunecols(m, n, I, J, V, min_entries; renumber=true)
    inds = findrepeated(J, min_entries)
    flatinds = Iterators.flatten(inds)
    if renumber
        shiftedinds = map(shiftrange, inds)
        new_J = _make_new_Inds(shiftedinds)
        new_n = length(shiftedinds)
    else
        new_J = _applyinds(J, flatinds)
        new_n = n
    end
    return (m, new_n, _applyinds(I, flatinds), new_J, _applyinds(V, flatinds))
end

function prunerows(m, n, I, J, V, min_entries; renumber=true)
    permmat = permutedims!(m, n, I, J, V)
    newmat = prunecols(permmat..., min_entries; renumber=renumber)
    permutedims!(permmat...)
    return permutedims!(newmat...)
end

### rotate

# This is copied from sparsematrix.jl, modified only to change the order of the arguments
function rot180!(m, n, I, J, V)
   @inbounds for i=1:length(I)
        I[i] = m - I[i] + 1
        J[i] = n - J[i] + 1
    end
    return (m, n, I, J, V)
end

function rotr90!(m, n, I, J, V)
    #old col inds are new row inds
   @inbounds for i=1:length(I)
        I[i] = m - I[i] + 1
    end
    return (n, m, J, I, V)
end

function rotl90!(m, n, I, J, V)
    #old row inds are new col inds
  @inbounds for i=1:length(J)
        J[i] = n - J[i] + 1
    end
    return (n, m, J, I, V)
end

### sprand

_empty_IJV(T) = (Int[], Int[], T[])

function sprand(m::Integer, n::Integer, density::AbstractFloat,
                rfn::Function, ::Type{T}=eltype(rfn(rng,1))) where T
    return sprand(Random.GLOBAL_RNG, m, n, density, rfn, eltype)
end

function sprand(rng::Random.AbstractRNG, m::Integer, n::Integer, density::AbstractFloat,
                rfn::Function, ::Type{T}=eltype(rfn(rng,1))) where T
    N = m*n
    if N == 0
        (I, J, V) = _empty_IJV(T)
    elseif N == 1
        if rand(rng) <= density
            (I, J, V) = ([1], [1], rfn(rng, 1))
        else
            (I, J, V) = _empty_IJV(T)
        end
    else
        I, J = SparseArrays.sprand_IJ(rng, m, n, density)
        V = rfn(rng, length(I))
    end
    return (m, n, I, J, V)
end

truebools(r::Random.AbstractRNG, n::Integer) = fill(true, n)
sprand(m::Integer, n::Integer, density::AbstractFloat) = sprand(Random.GLOBAL_RNG,m,n,density)

sprand(r::Random.AbstractRNG, m::Integer, n::Integer, density::AbstractFloat) = sprand(r,m,n,density,rand,Float64)
sprand(r::Random.AbstractRNG, ::Type{T}, m::Integer, n::Integer, density::AbstractFloat) where {T} = sprand(r,m,n,density,(r, i) -> rand(r, T, i), T)
sprand(r::Random.AbstractRNG, ::Type{Bool}, m::Integer, n::Integer, density::AbstractFloat) = sprand(r,m,n,density, truebools, Bool)
sprand(::Type{T}, m::Integer, n::Integer, density::AbstractFloat) where {T} = sprand(Random.GLOBAL_RNG, T, m, n, density)

end # module


text:
export MutablePose2Pose2Gaussian, PackedMutablePose2Pose2Gaussian
# import IncrementalInference: getSample

"""
    $TYPEDEF

Specialized Pose2Pose2 factor type (Gaussian), which allows for rapid accumulation of odometry information as a branch on the factor graph.
"""
mutable struct MutablePose2Pose2Gaussian  <: IIF.AbstractManifoldMinimize
  Z::MvNormal
  timestamp::DateTime
  MutablePose2Pose2Gaussian(Zij::MvNormal=MvNormal(zeros(3),diagm([0.01; 0.01; 0.001].^2)), timestamp::DateTime=now()) = new(Zij, timestamp)
end

DFG.getManifold(::MutablePose2Pose2Gaussian) = Manifolds.SpecialEuclidean(2)

function getSample(cf::CalcFactor{<:MutablePose2Pose2Gaussian})
  Xc = rand(cf.factor.Z)

  M = SpecialEuclidean(2)
  X = hat(M, Manifolds.Identity(M), Xc)
  return X
end

"""
    $SIGNATURES

Residual function for MutablePose2Pose2Gaussian.

Related

Pose2Pose2, Pose3Pose3, InertialPose3, DynPose2Pose2, Point2Point2, VelPoint2VelPoint2
"""
function (cf::CalcFactor{<:MutablePose2Pose2Gaussian})(X, p, q)
  M = getManifold(Pose2)
  q̂ = Manifolds.compose(M, p, exp(M, identity_element(M, p), X)) #for groups
  #TODO allocalte for vee! see Manifolds #412, fix for AD
  Xc = zeros(3)
  vee!(M, Xc, q, log(M, q, q̂))
  return Xc
end

## Serialization support

"""
$(TYPEDEF)
"""
mutable struct PackedMutablePose2Pose2Gaussian  <: IIF.PackedInferenceType
  datastr::String
  timestamp::Int64 # serialized in millisecond
  PackedMutablePose2Pose2Gaussian() = new()
  PackedMutablePose2Pose2Gaussian(x::String, ts::Int=datetime2unix(now())*1e3 |> Int) = new(x, ts)
end
function convert(::Type{MutablePose2Pose2Gaussian}, d::PackedMutablePose2Pose2Gaussian)
  return MutablePose2Pose2Gaussian(convert(SamplableBelief, d.datastr), timestamp=unix2datetime(d.timestamp*1e-3))
end
function convert(::Type{PackedMutablePose2Pose2Gaussian}, d::MutablePose2Pose2Gaussian)
  return PackedMutablePose2Pose2Gaussian(convert(PackedSamplableBelief, d.Z), datetime2unix(d.timestamp)*1e3 |> Int)
end


text:
using JuMP
using Clp

DEMANDS = [10011.22115007576, 9993.540857825825, 10029.493305561082, 9995.513032624265, 10040.407110819071, 10054.429367942921, 10059.76519849408, 9987.000868933652, 10035.54111237392, 9979.076532648018, 9973.268512810797, 10001.801162849199, 9977.013028384847, 10043.66957251199, 10052.449661342594, 9967.208355572673, 10000.69477091062, 10030.05030028466, 10004.993099108486, 9990.907712391863, 10048.51255532599, 10006.528026036742, 9995.689091324659, 9990.868630720339, 9997.195144776251, 9971.848670139489, 9945.536501292985, 9948.336004602528, 9997.577323187325, 9970.52023558712, 10069.401333593078, 10043.51650796613, 10047.132683479025, 10072.330774070775, 10031.89101682498, 10020.955889705769, 9934.338643331117, 10014.21570856934, 10000.845689612031, 10002.553802750448, 10029.850884936972, 9948.279191163045, 9971.233600233412, 9921.28110096277, 9969.291956593264, 9941.660187262343, 9970.955683172295, 10121.63150429881, 9890.3148430648, 10055.434949187049]

EXCESS = [5262.38201199916, 5860.755189729461, 4984.212894426032, 5363.8531442275635, 4620.34431447503, 3825.9519209367336, 5343.680950738127, 5233.478472018467, 5123.639554039856, 4635.80130213247, 4505.55692176387, 4218.191006304447, 5501.7375350818, 5206.40743566552, 4848.244095467865, 4922.152296817367, 5405.514169790265, 4156.03729448171, 4824.382398257035, 4000.2715617359813, 4390.76908613548, 4842.974771246174, 5220.6774383278635, 4178.847030718337, 5881.9659847155, 5380.195019617506, 4621.415119581128, 6104.345451062469, 5160.047534233344, 5370.740725813598, 4229.219081632929, 5048.335218378297, 4820.895310902216, 5144.298907980797, 4765.881280652281, 4707.799430473926, 4008.493355905563, 4681.96931929508, 6194.13279529005, 5615.828913696583, 4279.739726326251, 4650.89876777662, 4966.17758135978, 5151.236396106044, 5177.558350548357, 5293.24954830797, 4730.71724368456, 5493.7205140204305, 4524.476628100214, 4795.57008939593]
NUC_COST = 300
COAL_COST = 250
T = length(DEMANDS)

m = Model(Clp.Optimizer)
@variable(m, x[1:T] >= 0)
@variable(m, y[1:T] >= 0)
@variable(m, w[1:T] >= 0)
@variable(m, z[1:T] >= 0)
@objective(m, Min, sum(NUC_COST * y[i] + COAL_COST * x[i] for i = 1:T))

for t = 1:T
    @constraint(m, w[t] == sum(x[t:-1:max(1, t-19)]))
    @constraint(m, z[t] == sum(y[t:-1:max(1, t-14)]))
end

@constraint(m, w + z + EXCESS .≥ DEMANDS)
@constraint(m, 0.8 * z - 0.2 * w .≤ 0.2 * EXCESS)

println(m)

optimize!(m)

#status = termination_status(m)

#println("Solution status: ", status)

#println("Objective value: ", objective_value(m))
println("x = ", value.(x))
println("y = ", value.(y))


text:
@testset "Writer" begin

# NOTE: Not ideal to be testing with dependence on reader,
#       but reader is tested on expected byte representations so assuming reader is good this is a lot easier
@testset "single field $(typeof(x))" for x in [
    Int32(123),
    Int64(123),
    1.25,
    d128"1.25",
    true,
    nothing,
    DateTime(2021, 1, 2, 9, 30),
    BSONTimestamp(1, 2),
    BSONObjectId((
        0x1, 0x2, 0x3, 0x4,
        0x5, 0x6, 0x7, 0x8,
        0x9, 0xA, 0xB, 0xC,
    )),
    "test",
    BSONBinary([0x1, 0x2, 0x3]),
    uuid4(),
    BSONCode("f() = 1;"),
    BSONRegex("test", "abc"),
]
    buf = empty!(fill(0xff, 1000))
    writer = BSONWriter(buf)
    writer["x"] = x
    close(writer)
    BSONReader(buf, StrictBSONValidator())["x"][typeof(x)] == x
end

@testset "document" begin
    buf = empty!(fill(0xff, 1000))
    writer = BSONWriter(buf)
    writer["x"] = w -> begin
        w["a"] = 1
        w["b"] = 2
    end
    close(writer)
    reader = BSONReader(buf, StrictBSONValidator())
    @test reader["x"].type == BSON_TYPE_DOCUMENT
    @test reader["x"]["a"][Int] == 1
    @test reader["x"]["b"][Int] == 2
end

@testset "array" begin
    buf = empty!(fill(0xff, 1000))
    writer = BSONWriter(buf)
    writer["x"] = [1, 2, 3]
    close(writer)
    reader = BSONReader(buf, StrictBSONValidator())
    @test reader["x"].type == BSON_TYPE_ARRAY
    @test reader["x"][Vector{Int}] == [1, 2, 3]
end

@testset "array generator" begin
    buf = empty!(fill(0xff, 1000))
    writer = BSONWriter(buf)
    writer["x"] = (x * x for x in 1:3)
    close(writer)
    reader = BSONReader(buf, StrictBSONValidator())
    @test reader["x"].type == BSON_TYPE_ARRAY
    @test reader["x"][Vector{Int}] == [1, 4, 9]
end

@testset "array generator documents" begin
    buf = empty!(fill(0xff, 1000))
    writer = BSONWriter(buf)
    writer["x"] = (
        w -> begin
            w["a"] = x
            w["b"] = x * x
        end
        for x in 1:3
    )
    close(writer)
    reader = BSONReader(buf, StrictBSONValidator())
    @test reader["x"].type == BSON_TYPE_ARRAY
    @test reader["x"]["0"].type == BSON_TYPE_DOCUMENT
    @test reader["x"]["0"]["a"][Int] == 1
    @test reader["x"]["0"]["b"][Int] == 1
    @test reader["x"]["1"].type == BSON_TYPE_DOCUMENT
    @test reader["x"]["1"]["a"][Int] == 2
    @test reader["x"]["1"]["b"][Int] == 4
    @test reader["x"]["2"].type == BSON_TYPE_DOCUMENT
    @test reader["x"]["2"]["a"][Int] == 3
    @test reader["x"]["2"]["b"][Int] == 9
end

@testset "dict" begin
    buf = empty!(fill(0xff, 1000))
    writer = BSONWriter(buf)
    x = Dict{String, Any}("x" => 1, "y" => Dict{String, Any}("a" => 1, "b" => 2))
    writer[] = x
    close(writer)
    @test BSONReader(buf, StrictBSONValidator())[Any] == x
end

@testset "typed dict" begin
    buf = empty!(fill(0xff, 1000))
    writer = BSONWriter(buf)
    x = Dict{String, Int}("x" => 1, "y" => 2)
    writer[] = x
    close(writer)
    @test BSONReader(buf, StrictBSONValidator())[Any] == x
end

@testset "symbol name" begin
    buf = empty!(fill(0xff, 1000))
    writer = BSONWriter(buf)
    writer[:x] = 123
    close(writer)
    @test BSONReader(buf, StrictBSONValidator())["x"][Any] == 123
end

@testset "pair" begin
    buf = empty!(fill(0xff, 1000))
    writer = BSONWriter(buf)
    writer[] = "x" => 123
    close(writer)
    @test BSONReader(buf, StrictBSONValidator())["x"][Any] == 123
end

@testset "tuple of pairs" begin
    buf = empty!(fill(0xff, 1000))
    writer = BSONWriter(buf)
    writer[] = ("x" => 123, "y" => 1.25)
    close(writer)
    @test BSONReader(buf, StrictBSONValidator())["x"][Any] == 123
    @test BSONReader(buf, StrictBSONValidator())["y"][Any] == 1.25
end

@testset "empty struct" begin
    buf = empty!(fill(0xff, 1000))
    writer = BSONWriter(buf)
    writer["x"] = EmptyStruct()
    close(writer)
    @test BSONReader(buf, StrictBSONValidator())["x"][Any] == Dict{String, Any}()
    @test BSONReader(buf, StrictBSONValidator())["x"][EmptyStruct] == EmptyStruct()
end

end


text:
module Swarm

include("particles.jl")
include("functions.jl")
include("PSO.jl")
include("HPSO.jl")
include("IHPSO.jl")

export Particles
export pso, hpso, ihpso

end # module


text:
module Deriv_Utils_New
    using Statistics
    using Plots
    using BenchmarkTools
    using Profile
    using FFTW
    using Statistics
    using Optim
    using Images, FileIO, ImageIO
    using Printf

    #using Profile
    using LinearAlgebra

    # put the cwd on the path so Julia can find the module
    push!(LOAD_PATH, pwd()*"/main")
    using eqws

    #Jacobians of coeffs wrt image #######################################################

    function wst_S1_deriv(image::Array{Float64,2}, filter_hash::Dict; FFTthreads::Int=1)
        #=
        output: Nx, Nx, Nf
        =#
        # Use nthread threads for FFT -- but for Nx<512 nthread=1 is fastest.  Overhead?
        FFTW.set_num_threads(FFTthreads)

        # array sizes
        (Nx, Ny)  = size(image)
        (Nf, )    = size(filter_hash["filt_index"])

        # allocate output array
        dS1dα  = zeros(Float64, Nx, Nx, Nf)

        ## 1st Order
        im_fd = fft(image)

        # unpack filter_hash
        f_ind   = filter_hash["filt_index"]  # (J, L) array of filters represented as index value pairs
        f_val   = filter_hash["filt_value"]
        zarr = zeros(ComplexF64, Nx, Nx)  # temporary array to fill with zvals

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

        # Loop over filters
        for f = 1:Nf
            f_i = f_ind[f]  # CartesianIndex list for filter
            f_v = f_val[f]  # Values for f_i

            zarr[f_i] = (f_v.*f_v) .* im_fd[f_i]
            dS1dα[:,:,f] = 2 .* real.(P*zarr)

            zarr[f_i] .= 0   # reset zarr for next loop
        end
        return dS1dα
    end

    function wst_S20_deriv(image::Array{Float64,2}, filter_hash::Dict; FFTthreads::Int=1)
        #=
        output: Nx, Nx, Nf, Nf
        =#
        # Use nthread threads for FFT -- but for Nx<512 nthread=1 is fastest.  Overhead?
        FFTW.set_num_threads(FFTthreads)

        # array sizes
        (Nx, Ny)  = size(image)
        (Nf, )    = size(filter_hash["filt_index"])

        # allocate output array
        dS20dα  = zeros(Float64, Nx, Nx, Nf, Nf)

        # allocate image arrays for internal use
        im_rdc = Array{ComplexF64, 3}(undef, Nx, Ny, Nf)  # real domain, complex
        im_rd  = Array{Float64, 3}(undef, Nx, Ny, Nf)  # real domain, complex
        im_fd  = fft(image)

        # unpack filter_hash
        f_ind   = filter_hash["filt_index"]  # (J, L) array of filters represented as index value pairs
        f_val   = filter_hash["filt_value"]
        zarr = zeros(ComplexF64, Nx, Nx)  # temporary array to fill with zvals

        # make a FFTW "plan" a complex array, both forward and inverse transform
        P_fft  = plan_fft(im_fd)   # P_fft is an operator,  P_fft*im is fft(im)
        P_ifft = plan_ifft(im_fd)  # P_ifft is an operator, P_ifft*im is ifft(im)

        # Loop over filters
        for f = 1:Nf
            f_i = f_ind[f]  # CartesianIndex list for filter
            f_v = f_val[f]  # Values for f_i

            zarr[f_i] = f_v .* im_fd[f_i]
            Z_λ = P_ifft*zarr  # complex valued ifft of zarr
            zarr[f_i] .= 0   # reset zarr for next loop
            im_rdc[:,:,f] = Z_λ
            im_rd[:,:,f]  = abs.(Z_λ)
        end

        zarr = zeros(ComplexF64, Nx, Nx)  # temporary array to fill with zvals
        for f2 = 1:Nf
            f_i = f_ind[f2]  # CartesianIndex list for filter
            f_v = f_val[f2]  # Values for f_i
            uvec = im_rdc[:,:,f2] ./ im_rd[:,:,f2]
            for f1 = 1:Nf
                temp = P_fft*(im_rd[:,:,f1].*uvec)
                zarr[f_i] = f_v .* temp[f_i]

                Z1dZ2 = real.(P_ifft*zarr)
                #  It is possible to do this with rifft, but it is not much faster...
                #   Z1dZ2 = myrealifft(zarr)
                dS20dα[:,:,f1,f2] += Z1dZ2
                dS20dα[:,:,f2,f1] += Z1dZ2
                zarr[f_i] .= 0   # reset zarr for next loop
            end
        end
        return dS20dα
    end

    function wst_S12_deriv(image::Array{Float64,2}, filter_hash::Dict; FFTthreads::Int=1)
        #=
        output: Nx, Nx, Nf, Nf
        =#
        FFTW.set_num_threads(FFTthreads)

        # array sizes
        (Nx, Ny)  = size(image)
        (Nf, )    = size(filter_hash["filt_index"])

        # unpack filter_hash
        f_ind   = filter_hash["filt_index"]  # (J, L) array of filters represented as index value pairs
        f_val   = filter_hash["filt_value"]

        zarr = zeros(ComplexF64, Nf, Nx, Nx)  # temporary array to fill with zvals
        ## 1st Order
        im_fd_0 = fft(image)  # total power=1.0
        # make a FFTW "plan" for an array of the given size and type
        P = plan_ifft(im_fd_0)  # P is an operator, P*im is ifft(im)
        for f1 = 1:Nf
            f_i1 = f_ind[f1]  # CartesianIndex list for filter1
            f_v1 = f_val[f1]  # Values for f_i1
            zarr[f1, f_i1] = im_fd_0[f_i1] .* f_v1 #check if this notation is fine
        end
        zarrc =conj.(zarr)
        absz = abs.(zarr)

        #Implemented in the more memory intensive / parallelized ops way
        #fhash_indsarr = [hcat((x->[x[1], x[2]]).(fh)...)' for fh in fhash["filt_index"]] #Each elem gives you fhash["filt_index"] in non-CartesianIndices
        #dabsz = #Mem Req: Nf*Nx^2*sp*Nx^2 and sp~2%

        #Less Mem Route
        dabsz_absz = zeros(Float64, Nf, Nf, Nx, Nx)  # Lookup arrays for dS12 terms
        xigrid = reshape(CartesianIndices((1:Nx, 1:Nx)), (1, Nx^2))
        xigrid = hcat((x->[x[1], x[2]]).(xigrid)...) #2*Nx^2
        #Calculate pnz's here and pre-allocate memory

        for f1 = 1:Nf
            f_i1 = f_ind[f1]  # CartesianIndex list for filter1
            f_v1 = f_val[f1]  # Values for f_i1

            for f2 = 1:Nf
                f_i2 = f_ind[f2]  # CartesianIndex list for filter1
                f_v2 = f_val[f2]  # Values for f_i1
                pnz = intersect(f_i1, f_i2) #Fd indices where both filters are nonzero
                if length(pnz)!=0
                    pnz_arr = vcat((x->[x[1], x[2]]').(pnz)...) #pnz*2
                    Φmat =  exp.((-2π*im)/Nx .* ((pnz_arr .- 1) * (xigrid .- 1))) #pnz*Nx^2, BN1
                    f_v1pnz = f_v1[findall(in(pnz), f_i1)]
                    t2 = real.(zarrc[f1, pnz] .* Φmat) #Check, BN3
                    #println(size(f_v1pnz), size(t2), size(absz[f2, pnz]))
                    term = sum(((absz[f2, pnz] ./ absz[f1, pnz]) .* f_v1pnz) .* t2, dims=1) #p*Nx^2 -> 1*Nx^2, BN2
                    dabsz_absz[f1, f2, :, :] = reshape(term, (Nx, Nx))
                end
            end
        end
        dS12 = dabsz_absz + permutedims(dabsz_absz, [2, 1, 3, 4])
        return permutedims(dS12, [3, 4, 1, 2]) #MODIFIED
    end
    #=
    function wst_S12_derivfast(image::Array{Float64,2}, filter_hash::Dict; FFTthreads::Int=1)
        #=
        output: Nx, Nx, Nf, Nf
        =#
        FFTW.set_num_threads(FFTthreads)

        # array sizes
        (Nx, Ny)  = size(image)
        (Nf, )    = size(filter_hash["filt_index"])

        # unpack filter_hash
        f_ind   = filter_hash["filt_index"]  # (J, L) array of filters represented as index value pairs
        f_val   = filter_hash["filt_value"]

        zarr = zeros(ComplexF64, Nf, Nx, Nx)  # temporary array to fill with zvals
        ## 1st Order
        im_fd_0 = fft(image)  # total power=1.0
        # make a FFTW "plan" for an array of the given size and type
        P = plan_ifft(im_fd_0)  # P is an operator, P*im is ifft(im)
        for f1 = 1:Nf
            f_i1 = f_ind[f1]  # CartesianIndex list for filter1
            f_v1 = f_val[f1]  # Values for f_i1
            zarr[f1, f_i1] = im_fd_0[f_i1] .* f_v1 #check if this notation is fine
        end
        zarrc =conj.(zarr)
        absz = abs.(zarr)

        #Implemented in the more memory intensive / parallelized ops way
        #fhash_indsarr = [hcat((x->[x[1], x[2]]).(fh)...)' for fh in fhash["filt_index"]] #Each elem gives you fhash["filt_index"] in non-CartesianIndices
        #dabsz = #Mem Req: Nf*Nx^2*sp*Nx^2 and sp~2%

        #Less Mem Route
        dabsz_absz = zeros(Float64, Nf, Nf, Nx, Nx)  # Lookup arrays for dS12 terms
        xigrid = reshape(CartesianIndices((1:Nx, 1:Nx)), (1, Nx^2))
        xigrid = hcat((x->[x[1], x[2]]).(xigrid)...) #2*Nx^2
        #Calculate pnz's here and pre-allocate memory
        for f1  = 1:Nf
            f_i1 = f_ind[f1]  # CartesianIndex list for filter1
            f_v1 = f_val[f1]  # Values for f_i1
            for f2 = 1:Nf
                f_i2 = f_ind[f2]  # CartesianIndex list for filter1
                f_v2 = f_val[f2]  # Values for f_i1

        for f1 = 1:Nf
            f_i1 = f_ind[f1]  # CartesianIndex list for filter1
            f_v1 = f_val[f1]  # Values for f_i1

            for f2 = 1:Nf
                f_i2 = f_ind[f2]  # CartesianIndex list for filter1
                f_v2 = f_val[f2]  # Values for f_i1
                pnz = intersect(f_i1, f_i2) #Fd indices where both filters are nonzero
                if length(pnz)!=0
                    pnz_arr = vcat((x->[x[1], x[2]]').(pnz)...) #pnz*2
                    Φmat =  exp.((-2π*im)/Nx .* ((pnz_arr .- 1) * (xigrid .- 1))) #pnz*Nx^2, BN1
                    f_v1pnz = f_v1[findall(in(pnz), f_i1)]
                    t2 = real.(zarrc[f1, pnz] .* Φmat) #Check, BN3
                    #println(size(f_v1pnz), size(t2), size(absz[f2, pnz]))
                    term = sum(((absz[f2, pnz] ./ absz[f1, pnz]) .* f_v1pnz) .* t2, dims=1) #p*Nx^2 -> 1*Nx^2, BN2
                    dabsz_absz[f1, f2, :, :] = reshape(term, (Nx, Nx))
                end
            end
        end
        dS12 = dabsz_absz + permutedims(dabsz_absz, [2, 1, 3, 4])
        return permutedims(dS12, [3, 4, 1, 2]) #MODIFIED
    end
    =#
    function wst_S1S2_deriv(image::Array{Float64,2}, filter_hash::Dict; FFTthreads::Int=1)
        #=
        output: Nx, Nx, Nf+Nf^2
        =#
        # Use 2 threads for FFT
        FFTW.set_num_threads(FFTthreads)

        # array sizes
        (Nx, Ny)  = size(image)
        (Nf, )    = size(filter_hash["filt_index"])

        # allocate output array
        dS1dp  = zeros(Float64, Nf, Nx, Nx) #Mod zeros(Float64, Nx, Nx, Nf)
        dS2dp  = zeros(Float64, Nf, Nf, Nx, Nx) #Mod

        # allocate image arrays for internal use
        im_rdc_0_1 = Array{ComplexF64, 3}(undef, Nx, Ny, Nf)  # real domain, complex
        im_rd_0_1  = Array{Float64, 3}(undef, Nx, Ny, Nf)  # real domain, complex

        # Not sure what to do here -- for gradients I don't think we want these
        ## 0th Order
        #S0[1]   = mean(image)
        #norm_im = image.-S0[1]
        #S0[2]   = sum(norm_im .* norm_im)/(Nx*Ny)
        #norm_im ./= sqrt(Nx*Ny*S0[2])

        ## 1st Order
        im_fd_0 = fft(image)  # total power=1.0

        # unpack filter_hash
        f_ind   = filter_hash["filt_index"]  # (J, L) array of filters represented as index value pairs
        f_val   = filter_hash["filt_value"]

        f_ind_rev = [[CartesianIndex(mod1(Nx+2 - ci[1],Nx), mod1(Nx+2 - ci[2],Nx)) for ci in f_Nf] for f_Nf in f_ind]

        #CartesianIndex(17,17) .- f_ind
        # tmp arrays for all tmp[f_i] = arr[f_i] .* f_v opns
        zarr1 = zeros(ComplexF64, Nx, Nx)
        fz_fψ1 = zeros(ComplexF64, Nx, Nx)
        fterm_a = zeros(ComplexF64, Nx, Nx)
        fterm_ct1 = zeros(ComplexF64, Nx, Nx)
        fterm_ct2 = zeros(ComplexF64, Nx, Nx)
        #rterm_bt2 = zeros(ComplexF64, Nx, Nx)
        # make a FFTW "plan" for an array of the given size and type
        P = plan_ifft(im_fd_0)  # P is an operator, P*im is ifft(im)

        # Loop over filters
        for f1 = 1:Nf
            f_i1 = f_ind[f1]  # CartesianIndex list for filter1
            f_v1 = f_val[f1]  # Values for f_i1
            f_i1rev = f_ind_rev[f1]

            zarr1[f_i1] = f_v1 .* im_fd_0[f_i1] #F[z] #BUG: here????
            zarr1_rd = P*zarr1 #for frzi
            fz_fψ1[f_i1] = f_v1 .* zarr1[f_i1]
            fz_fψ1_rd = P*fz_fψ1
            dS1dp[f1, :,:] = 2 .* real.(fz_fψ1_rd) #real.(conv(I_λ, ψ_λ)) #Mod
            #CHECK; that this equals derivS1fast code

            #dS2 loop prep
            #ψ_λ1  = realspace_filter(Nx, f_i1, f_v1) #Slow
            #fcψ_λ1 = fft(conj.(ψ_λ1)) #Slow
            I_λ1 = sqrt.(abs2.(zarr1_rd))
            fI_λ1 = fft(I_λ1)
            rterm_bt1 = zarr1_rd./I_λ1 #Z_λ1/I_λ1_bar.
            rterm_bt2 = conj.(zarr1_rd)./I_λ1
            for f2 = 1:Nf
                f_i2 = f_ind[f2]  # CartesianIndex list for filter2
                f_v2 = f_val[f2]  # Values for f_i2

                fterm_a[f_i2] = fI_λ1[f_i2] .* (f_v2).^2 #fterm_a = F[I_λ1].F[ψ_λ]^2
                rterm_a = P*fterm_a                #Finv[fterm_a]
                fterm_bt1 = fft(rterm_a .* rterm_bt1) #F[Finv[fterm_a].*Z_λ1/I_λ1_bar] for T1
                fterm_bt2 = fft(rterm_a .* rterm_bt2) #F[Finv[fterm_a].*Z_λ1_bar/I_λ1_bar] for T2
                fterm_ct1[f_i1] = fterm_bt1[f_i1] .* f_v1    #fterm_b*F[ψ_λ]
                #fterm_ct2 = fterm_bt2 .* fcψ_λ1             #Slow
                #println(size(fterm_ct2[f_i1rev]),size(fterm_bt2[f_i1rev]),size(f_v1),size(conj.(f_v1)))
                fterm_ct2[f_i1rev] = fterm_bt2[f_i1rev] .* f_v1 #f_v1 is real
                #fterm_ct2slow = fterm_bt2 .* fcψ_λ1
                dS2dp[f1, f2, :, :] = real.(P*(fterm_ct1 + fterm_ct2)) #Mod THIS SHOULD be Nx*f2-1 + f1 or do the reshaping in loss func.

                #Reset f2 specific tmp arrays to 0
                fterm_a[f_i2] .= 0
                fterm_ct1[f_i1] .=0
                fterm_ct2[f_i1rev] .=0

            end
            # reset all reused variables to 0
            zarr1[f_i1] .= 0
            fz_fψ1[f_i1] .= 0

        end
        return permutedims(vcat(dS1dp, reshape(dS2dp, (Nf^2, Nx, Nx))), [2, 3, 1]) #Mod #MODIFIED
    end

    #Chi-Square Derivatives###############################################################
    function wst_S2_deriv_sum(image::Array{Float64,2}, filter_hash::Dict, wt::Array{Float64}; FFTthreads::Int=1)
        #=
        output:
        =#
        FFTW.set_num_threads(FFTthreads)

        # array sizes
        (Nx, Ny)  = size(image)
        (Nf, )    = size(filter_hash["filt_index"])

        if size(wt)!=(Nf, Nf) error("Wt has wrong shape") end

        f_ind   = filter_hash["filt_index"]  # (J, L) array of filters represented as index value pairs
        f_val   = filter_hash["filt_value"]
        f_ind_rev = [[CartesianIndex(mod1(Nx+2 - ci[1],Nx), mod1(Nx+2 - ci[2],Nx)) for ci in f_Nf] for f_Nf in f_ind]

        im_fd_0 = fft(image)
        P = plan_fft(im_fd_0)
        Pi = plan_ifft(im_fd_0)
        zarr1 = zeros(ComplexF64, Nx, Nx)
        Ilam1 = zeros(ComplexF64, Nx, Nx)
        fterm_bt1 = zeros(ComplexF64, Nx, Nx)
        fterm_bt2 = zeros(ComplexF64, Nx, Nx)
        ifterm = zeros(ComplexF64, Nx, Nx)
        dLossdα = zeros(Float64, Nx, Nx)
        for f1=1:Nf
            vlam1 = zeros(ComplexF64, Nx, Nx) #defined for each lambda1, computed by summing over lambda2
            f_i1 = f_ind[f1]
            f_v1 = f_val[f1]
            f_i1rev = f_ind_rev[f1]

            zarr1[f_i1] = im_fd_0[f_i1] .* f_v1
            zarr1_rd = Pi*(zarr1)
            Ilam1 = abs.(zarr1_rd)
            fIlam1 = P*Ilam1

            for f2=1:Nf
                f_i2 = f_ind[f2]
                f_v2 = f_val[f2]
                vlam1[f_i2] += wt[f1, f2] .* f_v2.^2
            end
            #vlam1= sum_λ2 w.ψl1,l2^2
            vlam1 = (Pi*(fIlam1 .* vlam1))./Ilam1 #sum over f2 -> (Nx, Nx) -> real space Vλ1/Iλ1
            fterm_bt1 = P*(vlam1 .* zarr1_rd)
            fterm_bt2 = P*(vlam1 .* conj.(zarr1_rd))
            ifterm[f_i1] = fterm_bt1[f_i1] .* f_v1
            ifterm[f_i1rev] += fterm_bt2[f_i1rev] .* f_v1
            dLossdα += real.(Pi*ifterm)

            #zero out alloc
            zarr1[f_i1] .= 0
            ifterm[f_i1] .=0
            ifterm[f_i1rev] .=0

        end
        return reshape(dLossdα, (Nx^2, 1))

    end

    function wst_S1S2_derivsum_comb(image::Array{Float64,2}, filter_hash::Dict, wt::Array{Float64}; FFTthreads::Int=1)
        FFTW.set_num_threads(FFTthreads)
        # array sizes
        (Nx, Ny)  = size(image)
        (Nf, )    = size(filter_hash["filt_index"])

        if size(wt)!=(Nf + (Nf*Nf),) error("Wt has wrong shape: the argument should contain both S1 and S2 weights") end
        wts1 = reshape(wt[1:Nf], (Nf, 1))
        wts2 = reshape(wt[Nf+1:end], (Nf, Nf))

        dwS1 = reshape(wst_S1_deriv(image, filter_hash), (Nx^2, Nf)) * wts1 #NOTE: Replace this with eqws_2DUtils.wst_S1_deriv
        dwS2 = reshape(wst_S2_deriv_sum(image, filter_hash, wts2), (Nx^2, 1))
        #println(size(dwS1), size(dwS2))
        return dwS1 + dwS2 #Why was another reshape needed here?

    end
    #=
    function wst_S20_deriv_sum_old(image::Array{Float64,2}, filter_hash::Dict, wt::Array{Float64}; FFTthreads::Int=1)
        # Sum over (f1,f2) filter pairs for S20 derivative.  This is much faster
        #   than calling wst_S20_deriv() because the sum can be moved inside the FFT.
        # Use FFTthreads threads for FFT -- but for Nx<512 FFTthreads=1 is fastest.  Overhead?
        # On Cascade Lake box, 4 is good for 2D, 8 or 16 for 3D FFTs
        FFTW.set_num_threads(FFTthreads)

        # array sizes
        (Nx, Ny)  = size(image)
        (Nf, )    = size(filter_hash["filt_index"])

        # allocate image arrays for internal use
        Uvec   = Array{ComplexF64, 3}(undef, Nx, Ny, Nf)  # real domain, complex
        im_rd  = Array{Float64, 3}(undef, Nx, Ny, Nf)  # real domain, complex
        im_fd  = fft(image) #Should this be the apodized image?

        # unpack filter_hash
        f_ind   = filter_hash["filt_index"]  # (J, L) array of filters represented as index value pairs
        f_val   = filter_hash["filt_value"]
        zarr    = zeros(ComplexF64, Nx, Nx)  # temporary array to fill with zvals

        # make a FFTW "plan" a complex array, both forward and inverse transform
        P_fft  = plan_fft(im_fd)   # P_fft is an operator,  P_fft*im is fft(im)
        P_ifft = plan_ifft(im_fd)  # P_ifft is an operator, P_ifft*im is ifft(im)

        # Loop over filters
        for f = 1:Nf
            f_i = f_ind[f]  # CartesianIndex list for filter
            f_v = f_val[f]  # Values for f_i

            zarr[f_i] = f_v .* im_fd[f_i]
            Z_λ = P_ifft*zarr  # complex valued ifft of zarr
            zarr[f_i] .= 0     # reset zarr for next loop
            im_rd[:,:,f] = abs.(Z_λ)
            Uvec[:,:,f]  = Z_λ ./ im_rd[:,:,f]
        end

        zarr = zeros(ComplexF64, Nx, Nx)  # temporary array to fill with zvals
        for f2 = 1:Nf
            f_i = f_ind[f2]  # CartesianIndex list for filter
            f_v = f_val[f2]  # Values for f_i

            Wtot = reshape( reshape(im_rd,Nx*Nx,Nf)*wt[:,f2], Nx, Nx) #SPEED: Preallocate mem for Wtot?
            temp = P_fft*(Wtot.*Uvec[:,:,f2])
            zarr[f_i] .+= f_v .* temp[f_i]
        end
        ΣdS20dα = 2* real.(P_ifft*zarr) #added here

        return reshape(ΣdS20dα, (Nx^2, 1))
    end
    =#

    function wst_S20_deriv_sum(image::Array{Float64,2}, filter_hash::Dict, wt::Array{Float64}; FFTthreads::Int=1)
        # Sum over (f1,f2) filter pairs for S20 derivative.  This is much faster
        #   than calling wst_S20_deriv() because the sum can be moved inside the FFT.
        # Use FFTthreads threads for FFT -- but for Nx<512 FFTthreads=1 is fastest.  Overhead?
        # On Cascade Lake box, 4 is good for 2D, 8 or 16 for 3D FFTs
        FFTW.set_num_threads(FFTthreads)

        # array sizes
        (Nx, Ny)  = size(image)
        (Nf, )    = size(filter_hash["filt_index"])

        # allocate image arrays for internal use
        Uvec   = Array{ComplexF64, 3}(undef, Nx, Ny, Nf)  # real domain, complex
        im_rd  = Array{Float64, 3}(undef, Nx, Ny, Nf)  # real domain, complex
        im_fd  = fft(image) #Should this be the apodized image?

        # unpack filter_hash
        f_ind   = filter_hash["filt_index"]  # (J, L) array of filters represented as index value pairs
        f_val   = filter_hash["filt_value"]
        zarr    = zeros(ComplexF64, Nx, Nx)  # temporary array to fill with zvals

        # make a FFTW "plan" a complex array, both forward and inverse transform
        P_fft  = plan_fft(im_fd)   # P_fft is an operator,  P_fft*im is fft(im)
        P_ifft = plan_ifft(im_fd)  # P_ifft is an operator, P_ifft*im is ifft(im)

        # Loop over filters
        for f = 1:Nf
            f_i = f_ind[f]  # CartesianIndex list for filter
            f_v = f_val[f]  # Values for f_i

            zarr[f_i] = f_v .* im_fd[f_i]
            Z_λ = P_ifft*zarr  # complex valued ifft of zarr
            zarr[f_i] .= 0     # reset zarr for next loop
            im_rd[:,:,f] = abs.(Z_λ)
            Uvec[:,:,f]  = Z_λ ./ im_rd[:,:,f]
        end

        zarr = zeros(ComplexF64, Nx, Nx)  # temporary array to fill with zvals
        for f2 = 1:Nf
            f_i = f_ind[f2]  # CartesianIndex list for filter
            f_v = f_val[f2]  # Values for f_i

            Wtot = reshape( reshape(im_rd,Nx*Nx,Nf)* (wt[:,f2] + wt'[:, f2]), Nx, Nx) #SPEED: Preallocate mem for Wtot?
            temp = P_fft*(Wtot.*Uvec[:,:,f2])
            zarr[f_i] .+= f_v .* temp[f_i]
        end
        ΣdS20dα = real.(P_ifft*zarr) #added here

        return reshape(ΣdS20dα, (Nx^2, 1))
    end

    function wst_S20_deriv_sum_norm(image::Array{Float64,2}, filter_hash::Dict, wt::Array{Float64}, std::Float64; FFTthreads::Int=1)
        # Sum over (f1,f2) filter pairs for S20 derivative.  This is much faster
        #   than calling wst_S20_deriv() because the sum can be moved inside the FFT.
        # Use FFTthreads threads for FFT -- but for Nx<512 FFTthreads=1 is fastest.  Overhead?
        # On Cascade Lake box, 4 is good for 2D, 8 or 16 for 3D FFTs
        FFTW.set_num_threads(FFTthreads)

        # array sizes
        (Nx, Ny)  = size(image)
        (Nf, )    = size(filter_hash["filt_index"])

        # allocate image arrays for internal use
        Uvec   = Array{ComplexF64, 3}(undef, Nx, Ny, Nf)  # real domain, complex
        im_rd  = Array{Float64, 3}(undef, Nx, Ny, Nf)  # real domain, complex
        im_fd  = fft(image) #Should this be the apodized image?

        # unpack filter_hash
        f_ind   = filter_hash["filt_index"]  # (J, L) array of filters represented as index value pairs
        f_val   = filter_hash["filt_value"]
        zarr    = zeros(ComplexF64, Nx, Nx)  # temporary array to fill with zvals

        # make a FFTW "plan" a complex array, both forward and inverse transform
        P_fft  = plan_fft(im_fd)   # P_fft is an operator,  P_fft*im is fft(im)
        P_ifft = plan_ifft(im_fd)  # P_ifft is an operator, P_ifft*im is ifft(im)

        # Loop over filters
        for f = 1:Nf
            f_i = f_ind[f]  # CartesianIndex list for filter
            f_v = f_val[f]  # Values for f_i

            zarr[f_i] = f_v .* im_fd[f_i]
            Z_λ = P_ifft*zarr  # complex valued ifft of zarr
            zarr[f_i] .= 0     # reset zarr for next loop
            im_rd[:,:,f] = abs.(Z_λ)
            Uvec[:,:,f]  = Z_λ ./ im_rd[:,:,f]
        end

        zarr = zeros(ComplexF64, Nx, Nx)  # temporary array to fill with zvals
        for f2 = 1:Nf
            f_i = f_ind[f2]  # CartesianIndex list for filter
            f_v = f_val[f2]  # Values for f_i

            Wtot = reshape(reshape(im_rd,Nx*Nx,Nf)* (wt[:,f2] + wt'[:, f2]), Nx, Nx) #SPEED: Preallocate mem for Wtot?
            temp = P_fft*(Wtot.*Uvec[:,:,f2])
            zarr[f_i] .+= f_v .* temp[f_i]
        end
        ΣdS20dα = real.(P_ifft*zarr).*(1.0/(Nx*std)) #added here

        return reshape(ΣdS20dα, (Nx^2, 1))
    end

    #Image Reconstruction Funcs##########################################################
    #=
    function image_recon_derivsum(input::Array{Float64, 2}, filter_hash::Dict, s_targ_mean::Array{Float64, 1}, s_targ_invcov::Array{Float64, 2}, pixmask::BitArray{2}, coeff_choice;
        FFTthreads::Int=1, optim_settings=Dict([("iterations", 10)]), coeff_mask=nothing) #add iso here and a check that returns error if both coeffmask is not nothing and iso is present.
        #=
        Cases to be handled by this function:
        S2 | S20 | S12
        Select pixels masked or all floating
        Select coefficients optimized wrt

        Cases NOT handled by this function:
        Iso
        Only S1
        No Regularization
        Non-Gaussian loss function
        Using the S0 params (mean and variance) in the loss function
        eqws_compute with norm=true
        Log
        Apodizing
        =#
        (Nf,) = size(filter_hash["filt_index"])
        println("Coeff mask:", (coeff_mask!=nothing))

        (Nx, Ny)  = size(input)
        if Nx != Ny error("Input image must be square") end

        if Nf == 0 error("filter hash corrupted") end
        if coeff_choice=="S2"
            println("S2")
            pixmask = pixmask[:] #flattened: Nx^2s

            if coeff_mask!=nothing
                if length(coeff_mask)!= (2+Nf + Nf^2) error("Wrong dim mask") end
                if size(s_targ_mean)[1]!=count((i->(i==true)), coeff_mask) error("s_targ_mean should only contain coeffs to be optimized") end
                if size(s_targ_invcov)[1]!=count((i->(i==true)), coeff_mask) error("s_targ_invcov should only have coeffs to be optimized") end
            else #No mask: all coeffs (default: All S1 and S2 coeffs will be optim wrt. Mean, var not incl)
                #Currently assuming inputs have mean, var params in case that's something we wanna optimize at some point
                if size(s_targ_mean)[1]!=(Nf+Nf^2) error("s_targ_mean should only contain coeffs to be optimized") end
                if (size(s_targ_invcov)!=(Nf+Nf^2, Nf+Nf^2))  error("s_targ_invcov of wrong size") end #(s_targ_invcov!=I) &
                #At this point both have dims |S1+S2|
                #Create coeff_mask subsetting 3:end
                coeff_mask = fill(true, 2+Nf+Nf^2)
                coeff_mask[1] = false
                coeff_mask[2] = false
            end

            num_freecoeff = count((i->(i==true)), coeff_mask) #Nf+Nf^2 if argument coeff_mask was empty
            num_freecoeffS1 = count((i->(i==true)), coeff_mask[3:Nf+2])
            num_freecoeffS2 = count((i->(i==true)), coeff_mask[Nf+3:end])

            #After this all cases have a coeffmask, and s_targ_mean and s_targ_invcov have the shapes of the coefficients that we want to select.

            #Optimization Hyperparameters
            norm=false
            numitns_dict = get(optim_settings, "iterations", 100)

            function loss_func(img_curr)
                #size(img_curr) must be (Nx^2, 1)
                s_curr = eqws_compute(reshape(img_curr, (Nx, Nx)),  filter_hash, doS2=true, doS12=false, doS20=false, norm=tonorm)[coeff_mask]
                neglogloss = 0.5 .* (s_curr - s_targ_mean)' * s_targ_invcov * (s_curr - s_targ_mean)
                return neglogloss[1]
            end

            function dloss(storage_grad, img_curr)
                #storage_grad, img_curr must be (Nx^2, 1)
                s_curr = eqws_compute(reshape(img_curr, (Nx, Nx)), filter_hash, doS2=true, doS12=false, doS20=false, norm=tonorm)[coeff_mask]
                diff = s_curr - s_targ_mean
                wts1s2 = zeros(Nf+Nf^2)
                wtall = reshape(diff' * s_targ_invcov, num_freecoeff)
                wts1s2[coeff_mask[3:end]] .= wtall
                dterm = wst_S1S2_derivsum_comb(reshape(img_curr, (Nx, Nx)), filter_hash, wts1s2, FFTthreads=FFTthreads)
                storage_grad .= reshape(dterm, (Nx^2, 1))
                storage_grad[pixmask, 1] .= 0 # better way to do this by taking pixmask as an argument wst_s2_deriv_sum?
                return
            end

            #Debugging stuff
            println("Diff check")
            eps = zeros(size(input))
            eps[1, 2] = 1e-4
            chisq1 = loss_func(input+eps./2)
            chisq0 = loss_func(input-eps./2)
            brute  = (chisq1-chisq0)/1e-4
            #df_brute = eqws_compute(reshape(input, (Nx, Nx)), filter_hash, doS2=true, doS12=false, doS20=false, norm=false)[coeff_mask] - s_targ_mean
            clever = reshape(zeros(size(input)), (Nx*Nx, 1))
            _bar = dloss(clever, reshape(input, (Nx^2, 1)))
            #println("dS1S2comb check")
            #wts1s2 = zeros(Nf+Nf^2)
            #wtall = reshape(df' * s_targ_invcov, num_freecoeff)
            #wts1s2[coeff_mask[3:end]] .= wtall
            #dterm = Deriv_Utils.wst_S1S2_derivsum_comb(reshape(input, (Nx, Nx)), filter_hash, wts1s2)
            #dS1S2comb = Deriv_Utils.wst_S1S2_derivfast(reshape(input, (Nx, Nx)), filter_hash)
            #dS1dp = permutedims(dS1S2comb[1:Nf, :, :], [2, 3, 1])
            #dS2dp = permutedims(dS1S2comb[Nf+1:end, :, :], [2, 3, 1])
            #dS2dp = reshape(dS2dp, Nx, Nx, Nf^2)
            #sum2 = zeros(Float64, Nx, Nx)
            #for i=1:Nf*Nf sum2 += (dS2dp[:,:,i].*wts1s2[Nf+i]) end
            #term1 = reshape(dS1dp, (Nx^2, Nf)) * reshape(wts1s2[1:Nf], (Nf, 1))
            #dterm_brute = term1 + reshape(sum2, (Nx^2, 1))
            println("Chisq Derve Check")
            println("Brute:  ",brute)
            println("Clever: ",clever[Nx*(1)+1])

            res = optimize(loss_func, dloss, reshape(input, (Nx^2, 1)), ConjugateGradient(), Optim.Options(iterations = numitns_dict, store_trace = true, show_trace = true))
            result_img = Optim.minimizer(res)
        elseif coeff_choice=="S20"
            println("S20")
            pixmask = pixmask[:] #flattened: Nx^2s

            if coeff_mask!=nothing
                if count((i->i==true), coeff_mask[1:Nf+2])!=0 println("Warning: S1 coeffs being double counted / S0 aren't 0") end #wrap
                if length(coeff_mask)!= (2+Nf + Nf^2) error("Wrong dim mask") end #wrap
                if size(s_targ_mean)[1]!=count((i->(i==true)), coeff_mask) error("s_targ_mean should only contain coeffs to be optimized") end
                if size(s_targ_invcov)[1]!=count((i->(i==true)), coeff_mask) error("s_targ_invcov should only have coeffs to be optimized") end
            else #No mask: all coeffs (default: Only S20 coeffs will be optim wrt.)
                #Currently assuming inputs have mean, var params in case that's something we wanna optimize at some point
                if size(s_targ_mean)[1]!=(Nf^2) error("s_targ_mean should only contain coeffs to be optimized") end
                if (size(s_targ_invcov)!=(Nf^2, Nf^2))  error("s_targ_invcov of wrong size") end #(s_targ_invcov!=I) &
                #At this point both have dims |S1+S2|
                #Create coeff_mask subsetting 3:end
                coeff_mask = fill(true, 2+Nf+Nf^2)
                coeff_mask[1:Nf+2] .= false #Mask out S0 and S1
            end

            num_freecoeff = count((i->(i==true)), coeff_mask) #Nf^2 if argument coeff_mask was empty

            #After this all cases have a coeffmask, and s_targ_mean and s_targ_invcov have the shapes of the coefficients that we want to select.

            #Optimization Hyperparameters
            tonorm = get(optim_settings, "norm", false)
            numitns_dict = get(optim_settings, "iterations", 100)

            function loss_func20(img_curr)
                #size(img_curr) must be (Nx^2, 1)
                s_curr = eqws_compute(reshape(img_curr, (Nx, Nx)),  filter_hash, doS2=false, doS12=false, doS20=true, norm=tonorm)[coeff_mask]
                neglogloss = 0.5 .* (s_curr - s_targ_mean)' * s_targ_invcov * (s_curr - s_targ_mean)
                return neglogloss[1]
            end

            function dloss20(storage_grad, img_curr)
                #storage_grad, img_curr must be (Nx^2, 1)
                s_curr = eqws_compute(reshape(img_curr, (Nx, Nx)), filter_hash, doS2=false, doS12=false, doS20=true, norm=tonorm)[coeff_mask]
                diff = s_curr - s_targ_mean
                wt = reshape(convert(Array{Float64, 2}, transpose(diff) * s_targ_invcov), (Nf, Nf))
                dterm = wst_S20_deriv_sum(reshape(img_curr, (Nx, Nx)), filter_hash, wt, FFTthreads=FFTthreads)
                storage_grad .= reshape(dterm, (Nx^2, 1))
                storage_grad[pixmask, 1] .= 0 # better way to do this by taking pixmask as an argument wst_s2_deriv_sum?
                return
            end

            #Debugging stuff
            println("Diff check")
            eps = zeros(size(input))
            eps[1, 2] = 1e-4
            chisq1 = loss_func20(input+eps./2)
            chisq0 = loss_func20(input-eps./2)
            brute  = (chisq1-chisq0)/1e-4
            #df_brute = eqws_compute(reshape(input, (Nx, Nx)), filter_hash, doS2=true, doS12=false, doS20=false, norm=false)[coeff_mask] - s_targ_mean
            clever = reshape(zeros(size(input)), (Nx*Nx, 1))
            _bar = dloss20(clever, reshape(input, (Nx^2, 1)))
            println("Chisq Derve Check")
            println("Brute:  ",brute)
            println("Clever: ",clever[Nx*(1)+1])

            res = optimize(loss_func20, dloss20, reshape(input, (Nx^2, 1)), ConjugateGradient(), Optim.Options(iterations = numitns_dict, store_trace = true, show_trace = true))
            result_img = Optim.minimizer(res)
        else
             if coeff_choice!="S12" error("Invalid coeff_choice") end
             println("S12")
             pixmask = pixmask[:] #flattened: Nx^2s

             if coeff_mask!=nothing
                 num_freecoeff = count((i->(i==true)), coeff_mask)
                 println("Number of coefficients being optimized= ", num_freecoeff)
                 if count((i->i==true), coeff_mask[1:Nf+2])!=0 println("Warning: S1 coeffs being double counted / S0 aren't 0") end #wrap
                 if length(coeff_mask)!= (2+Nf + Nf^2) error("Wrong dim mask") end #wrap
                 if size(s_targ_mean)[1]!=count((i->(i==true)), coeff_mask) error("s_targ_mean should only contain coeffs to be optimized") end
                 if size(s_targ_invcov)[1]!=count((i->(i==true)), coeff_mask) error("s_targ_invcov should only have coeffs to be optimized") end
             else #No mask: all coeffs (default: Only S20 coeffs will be optim wrt.)
                 error("S12 must be run with the coeff_mask argument")
                 #=
                 if size(s_targ_mean)[1]!=(Nf^2) error("s_targ_mean should only contain coeffs to be optimized") end
                 if (size(s_targ_invcov)!=(Nf^2, Nf^2))  error("s_targ_invcov of wrong size") end #(s_targ_invcov!=I) &
                 #At this point both have dims |S1+S2|
                 #Create coeff_mask subsetting 3:end
                 coeff_mask = fill(true, 2+Nf+Nf^2)
                 coeff_mask[1:Nf+2] .= false #Mask out S0 and S1 =#
             end

             #After this all cases have a coeffmask, and s_targ_mean and s_targ_invcov have the shapes of the coefficients that we want to select.

             #Optimization Hyperparameters
             tonorm = get(optim_settings, "norm", false)
             numitns_dict = get(optim_settings, "iterations", 100)
             coeffmask_S12 = coeff_mask[Nf+3:end]

             function loss_func12(img_curr)
                 #size(img_curr) must be (Nx^2, 1)
                 s_curr = eqws_compute(reshape(img_curr, (Nx, Nx)),  filter_hash, doS2=false, doS12=true, doS20=false, norm=tonorm)[coeff_mask]
                 neglogloss = 0.5 .* (s_curr - s_targ_mean)' * s_targ_invcov * (s_curr - s_targ_mean)
                 return neglogloss[1]
             end

             function dloss12(storage_grad, img_curr)
                 #storage_grad, img_curr must be (Nx^2, 1)
                 s_curr = eqws_compute(reshape(img_curr, (Nx, Nx)), filter_hash, doS2=false, doS12=true, doS20=false, norm=tonorm)[coeff_mask]
                 diff = s_curr - s_targ_mean
                 wt = transpose(diff) * s_targ_invcov
                 storage_grad .= reshape(wst_S12_deriv(reshape(img_curr, (Nx, Nx)), filter_hash, FFTthreads=FFTthreads), (Nx^2, Nf^2))[:, coeffmask_S12] * transpose(wt) #->(nx^2, 1)
                 storage_grad[pixmask, 1] .= 0 # better way to do this by taking pixmask as an argument wst_s2_deriv_sum?
                 return
             end

             #Debugging stuff
             println("Diff check")
             eps = zeros(size(input))
             eps[1, 2] = 1e-4
             chisq1 = loss_func12(input+eps./2)
             chisq0 = loss_func12(input-eps./2)
             brute  = (chisq1-chisq0)/1e-4
             #df_brute = eqws_compute(reshape(input, (Nx, Nx)), filter_hash, doS2=true, doS12=false, doS20=false, norm=false)[coeff_mask] - s_targ_mean
             clever = reshape(zeros(size(input)), (Nx*Nx, 1))
             _bar = dloss12(clever, reshape(input, (Nx^2, 1)))
             println("Chisq Derve Check")
             println("Brute:  ",brute)
             println("Clever: ",clever[Nx*(1)+1])

             res = optimize(loss_func12, dloss12, reshape(input, (Nx^2, 1)), ConjugateGradient(), Optim.Options(iterations = numitns_dict, store_trace = true, show_trace = true))
             result_img = Optim.minimizer(res)

        end
        return reshape(result_img, (Nx, Nx))
    end
    =#
    function image_recon_derivsum_regularized(input::Array{Float64, 2}, filter_hash::Dict, s_targ_mean::Array{Float64, 1}, s_targ_invcov, pixmask::BitArray{2}, dhc_args;
        FFTthreads::Int=1, optim_settings=Dict([("iterations", 10)]), coeff_mask::BitArray{1}=nothing, lambda=0.001,show_trace0::Bool=true) #add iso here and a check that returns error if both coeffmask is not nothing and iso is present.
        #=
        Cases to be handled by this function:
        S20
        Apodized or not
        Iso
        Log or Regular: If log then input, targ and cov should have been constructed using log
        Select pixels masked or all floating--not using for masked currently
        Select coefficients optimized wrt
        Lambda: regularization

        Cases NOT handled by this function:
        Only S1
        Non-Gaussian loss function
        Using the S0 params (mean and variance) in the loss function
        eqws_compute with norm=true
        =#
        (N1iso, Nf) = size(filter_hash["S1_iso_mat"])
        println("Coeff mask:", (coeff_mask!=nothing))

        (Nx, Ny)  = size(input)
        if Nx != Ny error("Input image must be square") end
        if !dhc_args[:doS20] error("Not implemented") end
        if Nf == 0 error("filter hash corrupted") end

        println("S20")
        #pixmask = pixmask[:] #flattened: Nx^2s #DEB
        ori_input = reshape(copy(input), (Nx^2, 1))

        if coeff_mask!=nothing
            if (!dhc_args[:iso]) & (count((i->i==true), coeff_mask[1:Nf+2])!=0) error("Code assumes S1 coeffs and S0 are masked out") end #wrap
            if (dhc_args[:iso]) & (count((i->i==true), coeff_mask[1:N1iso+2])!=0) error("Code assumes S1iso coeffs and S0 are masked out") end
            if size(s_targ_mean)[1]!=count((i->(i==true)), coeff_mask) error("s_targ_mean should only contain coeffs to be optimized") end
            if size(s_targ_invcov)[1]!=count((i->(i==true)), coeff_mask) error("s_targ_invcov should only have coeffs to be optimized") end
        else #No mask: all coeffs (default: Only S20 coeffs will be optim wrt.)
            error("You must supply a coeff_mask")
        end

        #Apodization adjustments
        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 Jacobian dA*P/dP
        if dhc_args[:apodize]
            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
        else
            dA = I
        end

        num_freecoeff = count((i->(i==true)), coeff_mask)

        #After this all cases have a coeffmask, and s_targ_mean and s_targ_invcov have the shapes of the coefficients that we want to select.
        #Optimization Hyperparameters
        tonorm = false
        numitns_dict = get(optim_settings, "iterations", 100)
        minmethod = get(optim_settings, "minmethod", ConjugateGradient())

        if (dhc_args[:doS20]) & (dhc_args[:iso])
            iso2nf2mask = zeros(Int64, Nf^2)
            M20 = filter_hash["S2_iso_mat"]
            for id=1:Nf^2 iso2nf2mask[M20.colptr[id]] = M20.rowval[id] end
            coeff_masks20 = coeff_mask[3+N1iso:end]
        else# (dhc_args[:doS20]) & (!dhc_args[:iso])
            coeff_masks20 = coeff_mask[3+Nf:end]
        end

        function augment_weights(inpwt::Array{Float64})
            if (dhc_args[:doS20]) & (dhc_args[:iso]) #wt (Nc) -> |S2_iso| -> |Nf^2|
                w_s2iso = zeros(size(filter_hash["S2_iso_mat"])[1])
                w_s2iso[coeff_masks20] .= inpwt
                w_nf2 = zeros(Nf^2)
                w_nf2 .= w_s2iso[iso2nf2mask]
                return reshape(w_nf2, (Nf, Nf))
            else#if (dhc_args[:doS20]) & (!dhc_args[:iso])
                w_nf2 = zeros(Nf^2)
                w_nf2[coeff_masks20] .= inpwt
                return reshape(w_nf2, (Nf, Nf))
            end
        end

        function loss_func20(img_curr::Array{Float64, 2})
            s_curr = eqws_compute_wrapper(img_curr, filter_hash, norm=tonorm; dhc_args...)[coeff_mask]
            regterm =  0.5*lambda*sum((adaptive_apodizer(img_curr) - adaptive_apodizer(input)).^2)
            lnlik = ( 0.5 .* (s_curr - s_targ_mean)' * s_targ_invcov * (s_curr - s_targ_mean))
            neglogloss = lnlik[1] + regterm
            return 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 = eqws_compute_wrapper(img_curr, filter_hash, norm=tonorm; dhc_args...)[coeff_mask]
            diff = s_curr - s_targ_mean
            #Add branches here
            wt = augment_weights(convert(Array{Float64,1}, reshape(transpose(diff) * s_targ_invcov, (length(diff),))))
            apdimg_curr = adaptive_apodizer(img_curr)
            storage_grad .= reshape((wst_S20_deriv_sum(apdimg_curr, filter_hash, wt, FFTthreads=FFTthreads)' + reshape(lambda.*(apdimg_curr - adaptive_apodizer(input)), (1, Nx^2))) * dA, (Nx, Nx))
            storage_grad[pixmask] .= 0 # better way to do this by taking pixmask as an argument wst_s2_deriv_sum?
            return
        end

        println("Diff check")
        eps = zeros(size(input))
        row, col = 24, 18 #convert(Int8, Nx/2), convert(Int8, Nx/2)+3
        epsmag = 1e-4
        eps[row, col] = epsmag
        chisq1 = loss_func20(input+eps./2) #DEB
        chisq0 = loss_func20(input-eps./2) #DEB
        brute  = (chisq1-chisq0)/epsmag
        #df_brute = eqws_compute(reshape(input, (Nx, Nx)), filter_hash, doS2=true, doS12=false, doS20=false, norm=false)[coeff_mask] - s_targ_mean
        clever = zeros(size(input)) #DEB
        meanval = mean(adaptive_apodizer(input)) ./ mean(wind_2d(Nx))
        _bar = dloss20(clever, input) #DEB
        println("Chisq Derve Check")
        println("Brute:  ",brute)
        println("Clever: ",clever[row, col], " Difference: ", brute - clever[row, col], " Mean ", meanval) #DEB
        println("Initial Loss ", loss_func20(input))
        res = optimize(loss_func20, dloss20, input, minmethod, Optim.Options(iterations = numitns_dict, store_trace = show_trace0, show_trace = show_trace0))
        result_img = Optim.minimizer(res)
        println("Final Loss ", loss_func20(result_img))
        return res, reshape(result_img, (Nx, Nx))
    end

    function image_recon_derivsum_regularized_S2(input::Array{Float64, 2}, filter_hash::Dict, s_targ_mean::Array{Float64, 1}, s_targ_invcov, pixmask::BitArray{2}, dhc_args;
        FFTthreads::Int=1, optim_settings=Dict([("iterations", 10)]), coeff_mask=nothing, lambda=0.001,show_trace0::Bool=true) #add iso here and a check that returns error if both coeffmask is not nothing and iso is present.
        #=
        Cases to be handled by this function:
        S20
        Apodized or not
        Iso
        Log or Regular: If log then input, targ and cov should have been constructed using log
        Select pixels masked or all floating--not using for masked currently
        Select coefficients optimized wrt
        Lambda: regularization

        Cases NOT handled by this function:
        Only S1
        Non-Gaussian loss function
        Using the S0 params (mean and variance) in the loss function
        eqws_compute with norm=true
        =#
        (N1iso, Nf) = size(filter_hash["S1_iso_mat"])
        println("Coeff mask:", (coeff_mask!=nothing))

        (Nx, Ny)  = size(input)
        if Nx != Ny error("Input image must be square") end
        if !dhc_args[:doS2] error("Not implemented") end
        if Nf == 0 error("filter hash corrupted") end

        println("S2")
        #pixmask = pixmask[:] #flattened: Nx^2s #DEB
        ori_input = reshape(copy(input), (Nx^2, 1))

        if coeff_mask!=nothing
            if (!dhc_args[:iso]) & (count((i->i==true), coeff_mask[1:Nf+2])!=0) error("Code assumes S1 coeffs and S0 are masked out") end #wrap
            if (dhc_args[:iso]) & (count((i->i==true), coeff_mask[1:N1iso+2])!=0) error("Code assumes S1iso coeffs and S0 are masked out") end
            if size(s_targ_mean)[1]!=count((i->(i==true)), coeff_mask) error("s_targ_mean should only contain coeffs to be optimized") end
            if size(s_targ_invcov)[1]!=count((i->(i==true)), coeff_mask) error("s_targ_invcov should only have coeffs to be optimized") end
        else #No mask: all coeffs (default: Only S20 coeffs will be optim wrt.)
            error("You must supply a coeff_mask")
        end

        #Apodization adjustments
        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 Jacobian dA*P/dP
        if dhc_args[:apodize]
            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
        else
            dA = I
        end

        num_freecoeff = count((i->(i==true)), coeff_mask)

        #After this all cases have a coeffmask, and s_targ_mean and s_targ_invcov have the shapes of the coefficients that we want to select.
        #Optimization Hyperparameters
        tonorm = false
        numitns_dict = get(optim_settings, "iterations", 100)
        minmethod = get(optim_settings, "minmethod", ConjugateGradient())

        if (dhc_args[:doS2]) & (dhc_args[:iso])
            iso2nf2mask = zeros(Int64, Nf^2)
            M20 = filter_hash["S2_iso_mat"]
            for id=1:Nf^2 iso2nf2mask[M20.colptr[id]] = M20.rowval[id] end
            coeff_masks20 = coeff_mask[3+N1iso:end]
        else# (dhc_args[:doS20]) & (!dhc_args[:iso])
            coeff_masks20 = coeff_mask[3+Nf:end]
        end

        function augment_weights(inpwt::Array{Float64})
            if (dhc_args[:doS2]) & (dhc_args[:iso]) #wt (Nc) -> |S2_iso| -> |Nf^2|
                w_s2iso = zeros(size(filter_hash["S2_iso_mat"])[1])
                w_s2iso[coeff_masks20] .= inpwt
                w_nf2 = zeros(Nf^2)
                w_nf2 .= w_s2iso[iso2nf2mask]
                return reshape(w_nf2, (Nf, Nf))
            else#if (dhc_args[:doS20]) & (!dhc_args[:iso])
                w_nf2 = zeros(Nf^2)
                w_nf2[coeff_masks20] .= inpwt
                return reshape(w_nf2, (Nf, Nf))
            end
        end

        function loss_func2(img_curr::Array{Float64, 2})
            s_curr = eqws_compute_wrapper(img_curr, filter_hash, norm=tonorm; dhc_args...)[coeff_mask]
            regterm =  0.5*lambda*sum((adaptive_apodizer(img_curr) - adaptive_apodizer(input)).^2)
            lnlik = ( 0.5 .* (s_curr - s_targ_mean)' * s_targ_invcov * (s_curr - s_targ_mean))
            neglogloss = lnlik[1] + regterm
            return neglogloss
        end
        #TODO: Need to replace the deriv_sums20 with a deriv_sum wrapper that handles all cases (S12, S20, S2)

        function dloss2(storage_grad::Array{Float64, 2}, img_curr::Array{Float64, 2})
            s_curr = eqws_compute_wrapper(img_curr, filter_hash, norm=tonorm; dhc_args...)[coeff_mask]
            diff = s_curr - s_targ_mean
            #Add branches here
            wt = augment_weights(convert(Array{Float64,1}, reshape(transpose(diff) * s_targ_invcov, (length(diff),))))
            apdimg_curr = adaptive_apodizer(img_curr)
            storage_grad .= reshape((wst_S2_deriv_sum(apdimg_curr, filter_hash, wt, FFTthreads=FFTthreads)' + reshape(lambda.*(apdimg_curr - adaptive_apodizer(input)), (1, Nx^2))) * dA, (Nx, Nx))
            storage_grad[pixmask] .= 0 # better way to do this by taking pixmask as an argument wst_s2_deriv_sum?
            return
        end

        println("Diff check")
        eps = zeros(size(input))
        row, col = 24, 18 #convert(Int8, Nx/2), convert(Int8, Nx/2)+3
        epsmag = 1e-4
        eps[row, col] = epsmag
        chisq1 = loss_func2(input+eps./2) #DEB
        chisq0 = loss_func2(input-eps./2) #DEB
        brute  = (chisq1-chisq0)/epsmag
        #df_brute = eqws_compute(reshape(input, (Nx, Nx)), filter_hash, doS2=true, doS12=false, doS20=false, norm=false)[coeff_mask] - s_targ_mean
        clever = zeros(size(input)) #DEB
        meanval = mean(adaptive_apodizer(input)) ./ mean(wind_2d(Nx))
        _bar = dloss2(clever, input) #DEB
        println("Chisq Derve Check")
        println("Brute:  ",brute)
        println("Clever: ",clever[row, col], " Difference: ", brute - clever[row, col], " Mean ", meanval) #DEB
        println("Initial Loss ", loss_func2(input))
        res = optimize(loss_func2, dloss2, input, minmethod, Optim.Options(iterations = numitns_dict, store_trace = show_trace0, show_trace = show_trace0))
        result_img = Optim.minimizer(res)
        println("Final Loss ", loss_func2(result_img))
        return res, reshape(result_img, (Nx, Nx))
    end

    function image_recon_derivsum_regularized_transformed_gaussian(input::Array{Float64, 2}, filter_hash::Dict, fs_targ_mean::Array{Float64, 1}, fs_targ_invcov::Array{Float64, 2}, pixmask::BitArray{2}, dhc_args, func, dfunc;
        FFTthreads::Int=1, optim_settings=Dict([("iterations", 10)]), coeff_mask=nothing, lambda=0.001) #add iso here and a check that returns error if both coeffmask is not nothing and iso is present.
        #=
        f(S) ~ G(fs_targ_mean, fs_targ_invcov)
        =#
        (N1iso, Nf) = size(filter_hash["S1_iso_mat"])
        println("Coeff mask:", (coeff_mask!=nothing))

        (Nx, Ny)  = size(input)
        if Nx != Ny error("Input image must be square") end
        if !dhc_args[:doS20] error("Not implemented") end
        if Nf == 0 error("filter hash corrupted") end

        println("S20")
        #pixmask = pixmask[:] #flattened: Nx^2s #DEB
        ori_input = reshape(copy(input), (Nx^2, 1))

        if coeff_mask!=nothing
            if (!dhc_args[:iso]) & (count((i->i==true), coeff_mask[1:Nf+2])!=0) error("Code assumes S1 coeffs and S0 are masked out") end #wrap
            if (dhc_args[:iso]) & (count((i->i==true), coeff_mask[1:N1iso+2])!=0) error("Code assumes S1iso coeffs and S0 are masked out") end
            #if length(coeff_mask)!= (2+Nf + Nf^2) error("Wrong dim mask") end #wrap
            if size(fs_targ_mean)[1]!=count((i->(i==true)), coeff_mask) error("s_targ_mean should only contain coeffs to be optimized") end
            if size(fs_targ_invcov)[1]!=count((i->(i==true)), coeff_mask) error("s_targ_invcov should only have coeffs to be optimized") end
        else #No mask: all coeffs (default: Only S20 coeffs will be optim wrt.)
            error("You must supply a coeff_mask")
        end

        #Apodization adjustments
        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 Jacobian dA*P/dP
        if dhc_args[:apodize]
            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
        else
            dA = I
        end

        num_freecoeff = count((i->(i==true)), coeff_mask)

        #After this all cases have a coeffmask, and s_targ_mean and s_targ_invcov have the shapes of the coefficients that we want to select.
        #Optimization Hyperparameters
        tonorm = false
        numitns_dict = get(optim_settings, "iterations", 100)
        minmethod = get(optim_settings, "minmethod", ConjugateGradient())

        if (dhc_args[:doS20]) & (dhc_args[:iso])
            iso2nf2mask = zeros(Int64, Nf^2)
            M20 = filter_hash["S2_iso_mat"]
            for id=1:Nf^2 iso2nf2mask[M20.colptr[id]] = M20.rowval[id] end
            coeff_masks20 = coeff_mask[3+N1iso:end]
        else# (dhc_args[:doS20]) & (!dhc_args[:iso])
            coeff_masks20 = coeff_mask[3+Nf:end]
        end

        function augment_weights(inpwt::Array{Float64})
            if (dhc_args[:doS20]) & (dhc_args[:iso]) #wt (Nc) -> |S2_iso| -> |Nf^2|
                w_s2iso = zeros(size(filter_hash["S2_iso_mat"])[1])
                w_s2iso[coeff_masks20] .= inpwt
                w_nf2 = zeros(Nf^2)
                w_nf2 .= w_s2iso[iso2nf2mask]
                return reshape(w_nf2, (Nf, Nf))
            else#if (dhc_args[:doS20]) & (!dhc_args[:iso])
                w_nf2 = zeros(Nf^2)
                w_nf2[coeff_masks20] .= inpwt
                return reshape(w_nf2, (Nf, Nf))
            end
        end

        function loss_func20(img_curr::Array{Float64, 2})
            #println(typeof(eqws_compute_wrapper(img_curr, filter_hash, norm=tonorm; dhc_args...)[coeff_mask]))
            s_curr = func(eqws_compute_wrapper(img_curr, filter_hash, norm=tonorm; dhc_args...)[coeff_mask]) #Edit
            regterm =  0.5*lambda*sum((adaptive_apodizer(img_curr) - adaptive_apodizer(input)).^2)
            lnlik = ( 0.5 .* (s_curr - fs_targ_mean)' * fs_targ_invcov * (s_curr - fs_targ_mean))
            neglogloss = lnlik[1] + regterm
            return 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 = eqws_compute_wrapper(img_curr, filter_hash, norm=tonorm; dhc_args...)[coeff_mask]
            diff = func(s_curr) - fs_targ_mean
            #Add branches here
            wt = convert(Array{Float64,1}, reshape(transpose(diff) * fs_targ_invcov, (length(diff),))) .* dfunc(s_curr)
            wtaug = augment_weights(wt)
            #println(size(wtaug))
            apdimg_curr = adaptive_apodizer(img_curr)
            storage_grad .= reshape((wst_S20_deriv_sum(apdimg_curr, filter_hash, wtaug, FFTthreads=FFTthreads)' + reshape(lambda.*(apdimg_curr - adaptive_apodizer(input)), (1, Nx^2))) * dA, (Nx, Nx))
            storage_grad[pixmask] .= 0 # better way to do this by taking pixmask as an argument wst_s2_deriv_sum?
            return
        end

        println("Diff check")
        eps = zeros(size(input))
        row, col = 24, 18 #convert(Int8, Nx/2), convert(Int8, Nx/2)+3
        epsmag = 1e-6
        eps[row, col] = epsmag
        chisq1 = loss_func20(input+eps./2) #DEB
        chisq0 = loss_func20(input-eps./2) #DEB
        brute  = (chisq1-chisq0)/epsmag
        #df_brute = eqws_compute(reshape(input, (Nx, Nx)), filter_hash, doS2=true, doS12=false, doS20=false, norm=false)[coeff_mask] - s_targ_mean
        clever = zeros(size(input)) #DEB
        meanval = mean(adaptive_apodizer(input)) ./ mean(wind_2d(Nx))
        _bar = dloss20(clever, input) #DEB
        println("Chisq Derve Check")
        println("Brute:  ",brute)
        println("Clever: ",clever[row, col], " Difference: ", brute - clever[row, col], " Mean ", meanval) #DEB
        println("Initial Loss ", loss_func20(input))
        res = optimize(loss_func20, dloss20, input, minmethod, Optim.Options(iterations = numitns_dict, store_trace = true, show_trace = true))
        result_img = Optim.minimizer(res)
        println("Final Loss ", loss_func20(result_img))
        return res, reshape(result_img, (Nx, Nx))
    end

    function image_recon_derivsum_separate(input::Array{Float64, 2}, iteration_input::Array{Float64, 2}, filter_hash::Dict, s_targ_mean::Array{Float64, 1}, s_targ_invcov, pixmask::BitArray{2}, dhc_args;
        FFTthreads::Int=1, optim_settings=Dict([("iterations", 10)]), coeff_mask=nothing, lambda=0.001) #add iso here and a check that returns error if both coeffmask is not nothing and iso is present.
        #=
        Cases to be handled by this function:
        S20
        Apodized or not
        Iso
        Log or Regular: If log then input, targ and cov should have been constructed using log
        Select pixels masked or all floating--not using for masked currently
        Select coefficients optimized wrt
        Lambda: regularization

        Cases NOT handled by this function:
        Only S1
        Non-Gaussian loss function
        Using the S0 params (mean and variance) in the loss function
        eqws_compute with norm=true
        =#
        (N1iso, Nf) = size(filter_hash["S1_iso_mat"])
        println("Coeff mask:", (coeff_mask!=nothing))

        (Nx, Ny)  = size(input)
        if Nx != Ny error("Input image must be square") end
        if !dhc_args[:doS20] error("Not implemented") end
        if Nf == 0 error("filter hash corrupted") end

        println("S20")
        #pixmask = pixmask[:] #flattened: Nx^2s #DEB
        ori_input = reshape(copy(input), (Nx^2, 1))

        if coeff_mask!=nothing
            if (!dhc_args[:iso]) & (count((i->i==true), coeff_mask[1:Nf+2])!=0) error("Code assumes S1 coeffs and S0 are masked out") end #wrap
            if (dhc_args[:iso]) & (count((i->i==true), coeff_mask[1:N1iso+2])!=0) error("Code assumes S1iso coeffs and S0 are masked out") end
            if size(s_targ_mean)[1]!=count((i->(i==true)), coeff_mask) error("s_targ_mean should only contain coeffs to be optimized") end
            if size(s_targ_invcov)[1]!=count((i->(i==true)), coeff_mask) error("s_targ_invcov should only have coeffs to be optimized") end
        else #No mask: all coeffs (default: Only S20 coeffs will be optim wrt.)
            error("You must supply a coeff_mask")
        end

        #Apodization adjustments
        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 Jacobian dA*P/dP
        if dhc_args[:apodize]
            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
        else
            dA = I
        end

        num_freecoeff = count((i->(i==true)), coeff_mask)

        #After this all cases have a coeffmask, and s_targ_mean and s_targ_invcov have the shapes of the coefficients that we want to select.
        #Optimization Hyperparameters
        tonorm = false
        numitns_dict = get(optim_settings, "iterations", 100)
        minmethod = get(optim_settings, "minmethod", ConjugateGradient())

        if (dhc_args[:doS20]) & (dhc_args[:iso])
            iso2nf2mask = zeros(Int64, Nf^2)
            M20 = filter_hash["S2_iso_mat"]
            for id=1:Nf^2 iso2nf2mask[M20.colptr[id]] = M20.rowval[id] end
            coeff_masks20 = coeff_mask[3+N1iso:end]
        else# (dhc_args[:doS20]) & (!dhc_args[:iso])
            coeff_masks20 = coeff_mask[3+Nf:end]
        end

        function augment_weights(inpwt::Array{Float64})
            if (dhc_args[:doS20]) & (dhc_args[:iso]) #wt (Nc) -> |S2_iso| -> |Nf^2|
                w_s2iso = zeros(size(filter_hash["S2_iso_mat"])[1])
                w_s2iso[coeff_masks20] .= inpwt
                w_nf2 = zeros(Nf^2)
                w_nf2 .= w_s2iso[iso2nf2mask]
                return reshape(w_nf2, (Nf, Nf))
            else#if (dhc_args[:doS20]) & (!dhc_args[:iso])
                w_nf2 = zeros(Nf^2)
                w_nf2[coeff_masks20] .= inpwt
                return reshape(w_nf2, (Nf, Nf))
            end
        end

        function loss_func20(img_curr::Array{Float64, 2})
            s_curr = eqws_compute_wrapper(img_curr, filter_hash, norm=tonorm; dhc_args...)[coeff_mask]
            regterm =  0.5*lambda*sum((adaptive_apodizer(img_curr) - adaptive_apodizer(input)).^2)
            lnlik = ( 0.5 .* (s_curr - s_targ_mean)' * s_targ_invcov * (s_curr - s_targ_mean))
            neglogloss = lnlik[1] + regterm
            return 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 = eqws_compute_wrapper(img_curr, filter_hash, norm=tonorm; dhc_args...)[coeff_mask]
            diff = s_curr - s_targ_mean
            #Add branches here
            wt = augment_weights(convert(Array{Float64,1}, reshape(transpose(diff) * s_targ_invcov, (length(diff),))))
            apdimg_curr = adaptive_apodizer(img_curr)
            storage_grad .= reshape((wst_S20_deriv_sum(apdimg_curr, filter_hash, wt, FFTthreads=FFTthreads)' + reshape(lambda.*(apdimg_curr - adaptive_apodizer(input)), (1, Nx^2))) * dA, (Nx, Nx))
            storage_grad[pixmask] .= 0 # better way to do this by taking pixmask as an argument wst_s2_deriv_sum?
            return
        end

        println("Diff check")
        eps = zeros(size(input))
        row, col = 24, 18 #convert(Int8, Nx/2), convert(Int8, Nx/2)+3
        epsmag = 1e-4
        eps[row, col] = epsmag
        chisq1 = loss_func20(input+eps./2) #DEB
        chisq0 = loss_func20(input-eps./2) #DEB
        brute  = (chisq1-chisq0)/epsmag
        #df_brute = eqws_compute(reshape(input, (Nx, Nx)), filter_hash, doS2=true, doS12=false, doS20=false, norm=false)[coeff_mask] - s_targ_mean
        clever = zeros(size(input)) #DEB
        meanval = mean(adaptive_apodizer(input)) ./ mean(wind_2d(Nx))
        _bar = dloss20(clever, input) #DEB
        println("Chisq Derve Check")
        println("Brute:  ",brute)
        println("Clever: ",clever[row, col], " Difference: ", brute - clever[row, col], " Mean ", meanval) #DEB
        println("Initial Loss ", loss_func20(input))
        res = optimize(loss_func20, dloss20, iteration_input, minmethod, Optim.Options(iterations = numitns_dict, store_trace = true, show_trace = true))
        result_img = Optim.minimizer(res)
        println("Final Loss ", loss_func20(result_img))
        return res, reshape(result_img, (Nx, Nx))
    end
end

max_stars_repo_name: andrew-saydjari/DHC

max_issues_repo_name: andrew-saydjari/DHC

max_forks_repo_name: andrew-saydjari/DHC


text:
const L = 2*10^4
iclasslabel = ["Deficient", "Perfect", "Abundant"]
iclass = zeros(Int64, 3)
iclass[1] = one(Int64) #by convention 1 is deficient

for n in 2:L
    if isprime(n)
        iclass[1] += 1
    else
        iclass[sign(divisorsum(n)-n)+2] += 1
    end
end

println("Classification of integers from 1 to ", L)
for i in 1:3
    println("   ", iclasslabel[i], ", ", iclass[i])
end

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
using BSON
using Flux
using Plots
using Images

# EPIC COLOR PALLETTE
col = [RGB(i) for i in [0x00ab9f,0x7409b7,0xd2cf00,0xff6f00] ] # mah precious colors

loadModelsGPU = (path;top="model-up.bson",down="model-down.bson") -> (
    filename = path*top;
    up = BSON.load(filename) ;
    model_up = up ;
    filename = path*down;
    down = BSON.load(filename) ;
    model_down = down;
    return model_up,model_down)

function calc_one(image,model_up,model_down)
    image = image |> gpu
    y_up= model_up(image)
    y_down= model_down(image)
    y_up = y_up |> cpu
    y_down = y_down |> cpu
    return y_up,y_down
end

a = loadModelsGPU("/home/freshstart/DiplomaThesisData/NeuralNetSaves_050/")
CNN_C5_up = a[1][:md] |> gpu
CNN_C5_down = a[2][:md] |> gpu

# FATHER FOIL RAE-2822 Loading for testing
FF=Array{Float32}(undef,(32,32,3,1))
#--- image loading-----
beep = "/home/freshstart/DiplomaThesisData/RAE_2822_baseline/"
# beep = "/home/freshstart/DiplomaThesisData/r_-50_50/50%/"
# beep = "/home/freshstart/DiplomaThesisData/r_-50_50/-50%/"
image = load(beep*"RAE_var.png")
image = channelview(image) #CWH not WHC ie. 3X32X32
Channels,Width,Height=size(image[:,:,:])
for C = 1:Channels-1,W = 1:Width,H = 1:Height
    FF[W,H,C,1] = image[C,W,H] #Switching from CWH to WHCN
end
#-----------------------

#----- Cp data loading--

file = open(beep*"Cp_RAE_up_var.dat")
lines=readlines(file);
CP_up=[]
for i in 1:length(lines)
    push!(CP_up,parse(Float32, lines[i]))
end
close(file)

file = open(beep*"Cp_RAE_down_var.dat")
lines=readlines(file);
CP_down=[]
for i in 1:length(lines)
    push!(CP_down,parse(Float32, lines[i]))
end
close(file)

bench = []

@time for i in 1:1000
t=@elapsed up,down=calc_one(FF,CNN_C5_up,CNN_C5_down)
push!(bench,t)
end

sum = 0
for i in bench
    sum+=i
end

println("Average evaluation time = ",sum/length(bench))

up,down=calc_one(FF,CNN_C5_up,CNN_C5_down)
x = [0,0.01,0.02,0.03,0.04,0.05,0.06,0.07,0.08,0.09,0.1,0.15,0.2,0.25,0.3,0.35,0.4,0.41,0.42,0.43,0.44,0.45,0.46,0.47,0.48,0.49,0.5,0.51,0.52,0.53,0.54,0.55,0.56,0.57,0.58,0.59,0.6,0.65,0.7,0.75,0.8,0.85,0.9,0.95,0.96,0.97,0.98,0.99,1]
col = [RGB(i) for i in [0x00ab9f,0x7409b7,0xd2cf00,0xff6f00] ] # mah precious colors
pgfplotsx()
plot([0],[0],c=RGB(0xFFFFFF),title = "Cp distribution between CFD simulation and Neural Net Prediction", yflip = true,label = "Deformation rate = -50%")
plot!(x,up,label = "CNN prediction distribution",c=col[1])
plot!(x,down,label="",c=col[1])
plot!(x,CP_up,label = "CFD simulation distribution",c=col[4])
plot!(x,CP_down,label="",c=col[4])

loss_up = Flux.mse(up,CP_up)
loss_down = Flux.mse(down,CP_down)