text:
using JuMP, GLPK, LinearAlgebra, DataFrames

# Define some input data about the test system
# Maximum power output of generators
g_max = [1000, 1000];
# Minimum power output of generators
g_min = [0, 300];
# Incremental cost of generators
c_g = [50, 100];
# Fixed cost of generators
c_g0 = [1000, 0]
# Incremental cost of wind generators
c_w = 50;
# Total demand
d = 1500;
# Wind forecast
w_f = 200;

# In this cell we create function solve_ed, which solves the economic dispatch problem for a given set of input parameters.
function solve_ed(g_max, g_min, c_g, c_w, d, w_f)
    #Define the economic dispatch (ED) model
    ed = Model(GLPK.Optimizer)

    # Define decision variables
    @variable(ed, 0 <= g[i = 1:2] <= g_max[i]) # power output of generators
    @variable(ed, 0 <= w <= w_f) # wind power injection

    # Define the objective function
    @objective(ed, Min, dot(c_g, g) + c_w * w)

    # Define the constraint on the maximum and minimum power output of each generator
    @constraint(ed, [i = 1:2], g[i] <= g_max[i]) #maximum
    @constraint(ed, [i = 1:2], g[i] >= g_min[i]) #minimum

    # Define the constraint on the wind power injection
    @constraint(ed, w <= w_f)

    # Define the power balance constraint
    @constraint(ed, sum(g) + w == d)

    # Solve statement
    optimize!(ed)

    # return the optimal value of the objective function and its minimizers
    return value.(g), value(w), w_f - value(w), objective_value(ed)
end

# Solve the economic dispatch problem
(g_opt, w_opt, ws_opt, obj) = solve_ed(g_max, g_min, c_g, c_w, d, w_f);

println("\n")
println("Dispatch of Generators: ", g_opt, " MW")
println("Dispatch of Wind: ", w_opt, " MW")
println("Wind spillage: ", w_f - w_opt, " MW")
println("\n")
println("Total cost: ", obj, "\$")

c_g_scale_df = DataFrame(Symbol("Dispatch of Generator 1(MW)") => Float64[],
               Symbol("Dispatch of Generator 2(MW)") => Float64[],
               Symbol("Dispatch of Wind(MW)") => Float64[],
               Symbol("Spillage of Wind(MW)") => Float64[],
               Symbol("Total cost(\$)") => Float64[])
for c_g1_scale = 0.5:0.1:3.0
    c_g_scale = [c_g[1] * c_g1_scale, c_g[2]] # update the incremental cost of the first generator at every iteration
    g_opt, w_opt, ws_opt, obj = solve_ed(g_max, g_min, c_g_scale, c_w, d, w_f) # solve the ed problem with the updated incremental cost
    push!(c_g_scale_df, (g_opt[1], g_opt[2], w_opt, ws_opt, obj))
end

ENV["COLUMNS"]=250 # Helps us display the complete table
c_g_scale_df

function solve_ed_inplace(c_w_scale)
    start = time()
    obj_out = Float64[]
    w_out = Float64[]
    g1_out = Float64[]
    g2_out = Float64[]

    ed = Model(GLPK.Optimizer)

    # Define decision variables
    @variable(ed, 0 <= g[i = 1:2] <= g_max[i]) # power output of generators
    @variable(ed, 0 <= w <= w_f ) # wind power injection

    # Define the objective function
    @objective(ed, Min, dot(c_g, g) + c_w * w)

    # Define the constraint on the maximum and minimum power output of each generator
    @constraint(ed, [i = 1:2], g[i] <= g_max[i]) #maximum
    @constraint(ed, [i = 1:2], g[i] >= g_min[i]) #minimum

    # Define the constraint on the wind power injection
    @constraint(ed, w <= w_f)

    # Define the power balance constraint
    @constraint(ed, sum(g) + w == d)

    optimize!(ed)

    for c_g1_scale = 0.5:0.01:3.0
        @objective(ed, Min, c_g1_scale*c_g[1]*g[1] + c_g[2]*g[2] + c_w_scale*c_w*w)
        optimize!(ed)
        push!(obj_out, objective_value(ed))
        push!(w_out, value(w))
        push!(g1_out, value(g[1]))
        push!(g2_out, value(g[2]))
    end
    elapsed = time() - start
    print(string("elapsed time: ", elapsed, " seconds"))
    return obj_out, w_out, g1_out, g2_out
end

solve_ed_inplace(2.0);

demandscale_df = DataFrame(Symbol("Dispatch of Generators(MW)") => Float64[],
               Symbol("Dispatch of Generator 2(MW)") => Float64[],
               Symbol("Dispatch of Wind(MW)") => Float64[],
               Symbol("Spillage of Wind(MW)") => Float64[],
               Symbol("Total cost(\$)") => Float64[])

for demandscale = 0.2:0.1:1.5
    g_opt,w_opt,ws_opt,obj = solve_ed(g_max, g_min, c_g, c_w, demandscale*d, w_f)

    push!(demandscale_df, (g_opt[1], g_opt[2], w_opt, ws_opt, obj))
end

demandscale_df

# In this cell we introduce binary decision u to the economic dispatch problem (function solve_ed)
function solve_uc(g_max, g_min, c_g, c_w, d, w_f)
    #Define the unit commitment (UC) model
    uc = Model(GLPK.Optimizer)

    # Define decision variables
    @variable(uc, 0 <= g[i=1:2] <= g_max[i]) # power output of generators
    @variable(uc, u[i = 1:2], Bin) # Binary status of generators
    @variable(uc, 0 <= w <= w_f ) # wind power injection

    # Define the objective function
    @objective(uc, Min, dot(c_g, g) + c_w * w)

    # Define the constraint on the maximum and minimum power output of each generator
    @constraint(uc, [i = 1:2], g[i] <= g_max[i]) #maximum
    @constraint(uc, [i = 1:2], g[i] >= g_min[i]) #minimum

    # Define the constraint on the wind power injection
    @constraint(uc, w <= w_f)

    # Define the power balance constraint
        @constraint(uc, sum(g) + w == d)

    # Solve statement
    optimize!(uc)

    status = termination_status(uc)
    if status != MOI.OPTIMAL
        return status, zeros(length(g)), 0.0, 0.0, zeros(length(u)), Inf
    end
    return status, value.(g), value(w), w_f - value(w), value.(u), objective_value(uc)
end

# Solve the economic dispatch problem
status, g_opt, w_opt, ws_opt, u_opt, obj = solve_uc(g_max, g_min, c_g, c_w, d, w_f);

println("\n")
println("Dispatch of Generators: ", g_opt[:], " MW")
println("Commitments of Generators: ", u_opt[:])
println("Dispatch of Wind: ", w_opt, " MW")
println("Wind spillage: ", w_f - w_opt, " MW")
println("\n")
println("Total cost: ", obj, "\$")

uc_df = DataFrame(Symbol("Commitment of Generator 1(MW)") => Float64[],
               Symbol("Commitment of Generator 2(MW)") => Float64[],
               Symbol("Dispatch of Generator 1(MW)") => Float64[],
               Symbol("Dispatch of Generator 2(MW)") => Float64[],
               Symbol("Dispatch of Wind(MW)") => Float64[],
               Symbol("Spillage of Wind(MW)") => Float64[],
               Symbol("Total cost(\$)") => Float64[])

for demandscale = 0.2:0.1:1.5
    status, g_opt, w_opt, ws_opt, u_opt, obj = solve_uc(g_max, g_min, c_g, c_w, demandscale*d, w_f)

    if status == MOI.OPTIMAL
    push!(uc_df, (u_opt[1], u_opt[2], g_opt[1], g_opt[2], w_opt, ws_opt, obj))
    else
        println("Status: $status for demandscale = $demandscale \n")
    end
end

uc_df


text:
using Distributed
using K8sClusterManagers
using K8sClusterManagers: DEFAULT_NAMESPACE, NAMESPACE_FILE
using K8sClusterManagers: DEFAULT_WORKER_CPU, DEFAULT_WORKER_MEMORY
using K8sClusterManagers: create_pod, delete_pod, exec_pod, get_pod, label_pod, pod_status,
    wait_for_running_pod
using LibGit2: LibGit2
using Mocking: Mocking, @patch, apply
using Mustache: Mustache, render
using Random: randstring
using Test
using YAML: YAML
using kubectl_jll: kubectl

Mocking.activate()

include("utils.jl")

@testset "K8sClusterManagers" begin
    include("namespace.jl")
    include("pod.jl")
    include("native_driver.jl")

    # Tests that interact with a real, usually local, Kubernetes cluster
    if parse(Bool, get(ENV, "K8S_CLUSTER_TESTS", "true"))
        @testset "Kubernetes Cluster Tests" begin
            include("cluster.jl")
        end
    end
end

ext: jl


text:
# Type aliases
const AxisRange = Tuple{Float64, Float64}
const AxisTickData  = Tuple{Float64, Tuple{Float64,Float64}, Int}

"""
    Axes(kind, ranges, tickdata, ticklabels, perspective, camera, options)

Return an `Axes` object, containing the graphical specifications of the
coordinate system of a plot.

Axes are determined by their `kind`, which may be `:axes2d` for 2-D plots,
`:axes3d` for 3-D plots, and `:polar` for polar plots. The rest of its fields are:

* **`ranges`**: boundaries of the different axes/scales. They are given as a dictionary
    whose keys are symbols with the name of the axis (`:x`, `:y`, `:z`, `:c`),
    and whose values are tuples with two float values — the minimum and maximum
    values, respectively. The range `(Inf, -Inf)` describes an undefined axis.
* **`tickdata`**: numeric specifications of the "ticks" that are drawn on the axes.
    They are given as a dictionary whose keys are the names of the axis (as for `range`),
    and whose values are tuples that contain for that axis: (1) a `Float64` with
    the "minor" value between consecutive ticks; (2) a 2-tuple of `Float64` with
    the ends of the axis ticks; and (3) an `Int` with the number of minor ticks
    between "major", numbered ticks.
* **`ticklabels`**: transformations between tick values and labels. They are given
    as a dictionary whose keys are the names of the axis, and whose values are
    functions that accept a number as argument, and return a string with the
    text that must be written at major ticks. (This only works for the X and Y axes).
* **`perspective`**: A `Vector{Int}` that contains the "rotation" and "tilt" angles
    that are used to project 3-D axes on the plot plane. (Only for 3-D plots)
* **`camera`**: A `Vector{Float64}` with the camera parameters (camera position,
    view center and "up" vector, only used for 3-D plots).
* **`options`**: A `Dict{Symbol, Int}` with extra options that control the visualization
    of the axes. Currently supported options are:
    + `options[:scale]`, an integer code that defines what axes are must be
        plotted in log scale or reversed (cf. the function `GR.setscale`).
    + `options[:grid] = 0` to hide the plot grid, or any other value to show it.
    + `options[:tickdir]` to determine how the ticks are drawn
        (positive value to draw them inside the plot area, negative value to
        draw them outside, or `0` to hide them).
    + `options[:gr3] ≠ 0` to identify if the axes are a 3-D scene defined for the `gr3` interface.
    + `options[:radians] = 0` to transform angular values to degrees in polar axes.
"""
struct Axes
    kind::Symbol
    ranges::Dict{Symbol, AxisRange}
    tickdata::Dict{Symbol, AxisTickData}
    ticklabels::Dict{Symbol, <:Function}
    perspective::Vector{Int}
    camera::Vector{Float64}
    options::Dict{Symbol, Int}
end

"""
    Axes(kind::Symbol [; kwargs...])

Return an `Axes` object with selected parameters given by keyword arguments.

This constructor only requires the `kind` of the axes (`:axes2d`, `:axes3d` or
`:axespolar`), such that all the other parameters are passed as keyword arguments.
Null or empty values are used by default for the parameters that are not given.
"""
Axes(kind::Symbol; ranges = Dict{Symbol, AxisRange}(),
    tickdata = Dict{Symbol, AxisTickData}(),
    ticklabels = Dict{Symbol, Function}(),
    perspective = Int[],
    camera = Float64[],
    options = Dict{Symbol, Int}()) =
    Axes(kind, ranges, tickdata, ticklabels, perspective, camera, options)

"""
    Axes(kind, geoms::Array{<:Geometry} [; kwargs...]) where kind

Return an `Axes` object defined by the `kind` of the axes, and a vector of
[`Geometry`](@ref) objects that are meant to be plotted inside the axes, which
are used to calculate the different axis limits, ticks, etc.
Keyword arguments are used to override the default calculations.
"""
function Axes(kind, geoms::Array{<:Geometry}; grid=1, kwargs...)
    # Set limits based on data and scale
    scale = set_scale(; kwargs...)
    ranges = minmax(geoms, scale)
    adjustranges!(ranges; kwargs...)
    ticklabels = Dict{Symbol, Function}()
    perspective = [0, 0]
    camera = zeros(9)
    options = Dict{Symbol, Int}(:scale => scale, :grid => Int(grid))
    # Special cases dependin on axis kind
    if kind == :axes2d
        tickdata = set_ticks(ranges, 5, (:x, :y); kwargs...)
        set_ticklabels!(ticklabels; kwargs...)
    elseif kind == :axes3d
        tickdata = set_ticks(ranges, 2, (:x, :y, :z); kwargs...)
        perspective = [Int(get(kwargs, :rotation, 40)), Int(get(kwargs, :tilt, 70))]
        if get(kwargs, :gr3, false)
            options[:gr3] = 1
            cameradistance = get(kwargs, :cameradistance, 3.0)
            camera = set_camera(cameradistance, perspective...; kwargs...)
        end
    elseif kind == :polar
        tickdata = set_ticks(ranges, 2, (:x, :y); kwargs..., xlog=false, ylog=false, xflip=false, yflip=false)
        options[:scale] = 0
        if !get(kwargs, :radians, true)
            options[:radians] = 0
        end
    else # Not defined
        tickdata = set_ticks(ranges, 0, (:x, :y); kwargs...)
        options[:scale] = 0
        options[:grid] = 1
    end
    haskey(kwargs, :tickdir) && (options[:tickdir] = kwargs[:tickdir])
    Axes(kind, ranges, tickdata, ticklabels, perspective, camera, options)
end

# Calculation of data ranges

"""
    minmax(args...)

Calculate the ranges of a given array of values, a `Geometry` or a collection of `Geometry`
"""
# Nested definitions of minmax: for array, Geometry and Array of Geometries
function minmax(x::AbstractVecOrMat{<:Real}, (min_prev, max_prev), skipnegative)
    isempty(x) && return (float(min_prev), float(max_prev))
    if skipnegative
        x = replace(v -> v > 0 ? v : NaN, x)
    end
    x0, x1 = extrema64(x)
    newmin = min(x0, min_prev)
    newmax = max(x1, max_prev)
    return (newmin, newmax)
end

function minmax(g::Geometry, xminmax, yminmax, zminmax, cminmax, scale)
    # skip negative values if log scales:
    xminmax = minmax(g.x, xminmax, scale & GR.OPTION_X_LOG != 0)
    yminmax = minmax(g.y, yminmax, scale & GR.OPTION_Y_LOG != 0)
    zminmax = minmax(g.z, zminmax, scale & GR.OPTION_Z_LOG != 0)
    if isempty(g.c)
        cminmax = minmax(g.z, cminmax, false)
    else
        cminmax = minmax(g.c, cminmax, false)
    end
    return xminmax, yminmax, zminmax, cminmax
end

function minmax(geoms::Array{<:Geometry}, scale=0)
    # Calculate ranges of given values
    xminmax = yminmax = zminmax = cminmax = extrema64(Float64[])
    for g in geoms
        xminmax, yminmax, zminmax, cminmax = minmax(g, xminmax, yminmax, zminmax, cminmax, scale)
    end
    # Adjust ranges
    (xminmax == (Inf, -Inf)) && (xminmax = (0.0, 1.0))
    (yminmax == (Inf, -Inf)) && (yminmax = (0.0, 1.0))
    xminmax = fix_minmax(xminmax...)
    yminmax = fix_minmax(yminmax...)
    zminmax = fix_minmax(zminmax...)
    # Return values
    return Dict(:x => xminmax, :y => yminmax, :z => zminmax, :c => cminmax)
end

"""
    extrema64(a)

Compute both the minimum and maximum element and return them as a 2-tuple.
This is approximately like `extrema` in the standard library, but always
returns a tuple of `Float64` values, ignores `NaN`, and returns `(Inf, -Inf)`
for inputs that are empty or only contain `NaN`.
"""
function extrema64(a)
    amin =  typemax(Float64)
    amax = -typemax(Float64)
    for el in a
        if !isnan(el)
            if isnan(amin) || el < amin
                amin = Float64(el)
            end
            if isnan(amax) || el > amax
                amax = Float64(el)
            end
        end
    end
    amin, amax
end

"""
    fix_minmax(a, b)

Adjust `a` and `b` to avoid that they coincide.
"""
function fix_minmax(a, b)
    if a == b
        a -= a != 0 ? 0.1 * a : 0.1
        b += b != 0 ? 0.1 * b : 0.1
    end
    a, b
end

# Adjust ranges

"""
    set_limits(limits1, limits2)

Return a tuple of lower and upper limits, using the values given in `limits1`,
or those from `limits2` if either value of `limits1` is a `Nothing`.
"""
set_limits(::Tuple{Nothing, Nothing}, limits2) = float.(limits2)
set_limits((n, max1)::Tuple{Nothing, T}, (min2, max2)) where T = float.((min2, max1))
set_limits((min1, n)::Tuple{T, Nothing}, (min2, max2)) where T = float.((min1, max2))
set_limits(limits1, limits2) = float.(limits1)

"""
    adjustranges!(ranges; kwargs...)

Adjust the pre-calculated ranges of `Axes` — see [`minmax`](@ref),
to make them near to integers or "small decimals". Explicit axes limits or log scales
are considered, through keyword arguments `xlim`, `ylim`, etc. (given as
2-tuples) or `xlog` `ylog`, etc. (given as `Bool`).
"""
function adjustranges!(ranges::Dict{Symbol, AxisRange}; kwargs...)
    for axname in keys(ranges)
        keylim = Symbol(axname, :lim)
        if haskey(kwargs, keylim)
            ranges[axname] = set_limits(kwargs[keylim], ranges[axname])
        else
            keylog = Symbol(axname, :log)
            if !get(kwargs, keylog, false)
                amin, amax = ranges[axname]
                if isfinite(amin) && isfinite(amax)
                    ranges[axname] = GR.adjustlimits(amin, amax)
                end
            end
        end
    end
end

# Set ticks for the different types of axes

"""
    set_ticks(ranges, major_count, coordinates; kwargs...)

Define the tick numeric specifications of a given `Axes`, taking into
account its calculated `ranges` — see [`minmax`](@ref), with a given
number of minor ticks between major ticks (`major_count`). The specifications
are given for the axes defined in `coordinates` as a tuple of symbols
(e.g. `(:x, :y)` for the X and Y axes).
Keyword arguments are used to adjust the tick intervals
and limits if the axes are set in log scale (`xlog`, `ylog`, etc.,
given as `Bool` values), or if they are reversed (`xflip`, `yflip`, etc.,
also given as `Bools`).
"""
function set_ticks(ranges, major_count, coordinates; kwargs...)
    Dict(c => set_axis(c, ranges[c], major_count; kwargs...) for c ∈ coordinates)
end

function set_axis(axname, axrange, major; kwargs...)
    keyticks = Symbol(axname, :ticks)
    if haskey(kwargs, keyticks)
        tick, major = kwargs[keyticks]
    else
        tick = GR.tick(axrange...) / major
    end
    org = get(kwargs, Symbol(axname, :flip), false) ? reverse(axrange) : axrange
    return tick, org, major
end

# Set tick labels - only working for axes2d

"""
    set_ticklabels!(ticklabels; kwargs...)

Modify the given tick label transformation functions for `Axes`, using the
keyword arguments `xticklabels` or `yticklabels` (if they are given).

The keyword arguments may be functions that transform numbers into strings,
or collection of strings that are associated to the sequence of integers `1, 2, ...`.
"""
function set_ticklabels!(ticklabels; kwargs...)
    if haskey(kwargs, :xticklabels)
        ticklabels[:x] = ticklabel_fun(kwargs[:xticklabels])
    end
    if haskey(kwargs, :yticklabels)
        ticklabels[:y] = ticklabel_fun(kwargs[:yticklabels])
    end
    return nothing
end

function ticklabel_fun(f::Function)
    (x, y, svalue, value) -> begin
        GR.textext(x, y, string(f(value)))
        return Int32(0)
    end
end

function ticklabel_fun(labels::AbstractVecOrMat{T}) where T <: AbstractString
    (x, y, svalue, value) -> begin
        pos = findfirst(t->(value≈t), collect(1:length(labels)))
        lab = (pos === nothing) ? "" : labels[pos]
        GR.textext(x, y, lab)
        return Int32(0)
    end
end

# Set scale

"""
    set_scale(; kwargs...)

Set the scale characteristics of the axes (logarithmic or flipped axes),
taking into account the and keyword arguments that determine
which axes are in logarithmic scale (`xlog`, `ylog`, etc., given as `Bool` values)
and which ones are flipped (`xflip`, `yflip`, etc. also given as `Bool`s).

The result is an integer code used by the low level function [`GR.setscale`](@ref).
"""
function set_scale(; kwargs...)
    scale = 0
    get(kwargs, :xlog, false) && (scale |= GR.OPTION_X_LOG)
    get(kwargs, :ylog, false) && (scale |= GR.OPTION_Y_LOG)
    get(kwargs, :zlog, false) && (scale |= GR.OPTION_Z_LOG)
    get(kwargs, :xflip, false) && (scale |= GR.OPTION_FLIP_X)
    get(kwargs, :yflip, false) && (scale |= GR.OPTION_FLIP_Y)
    get(kwargs, :zflip, false) && (scale |= GR.OPTION_FLIP_Z)
    return scale
end

# Camera rotation

function _rotate!(v, angle) # Around the third (vertical) axis
    c = cosd(angle)
    s = sind(angle)
    v[1:2] .= (c*v[1] - s*v[2], s*v[1] + c*v[2])
end

function _tilt!(v, angle) # Around the first (transversal) axis (clockwise)
    c = cosd(angle)
    s = sind(angle)
    v[2:3] .= (c*v[2]  + s*v[3], -s*v[2] + c*v[3])
end

function _focus!(params, target) # Rotate to focus on target point
    oldaxis = normalize([params[4]-params[1], params[5]-params[2], params[6]-params[3]])
    newaxis = normalize([target[1]-params[1], target[2]-params[2], target[3]-params[3]])
    # Use the Rodrigues formula: v2 = v1 + sin(w)û×v1 + (1-cos(w))û×(û×v1)
    f1 = oldaxis × newaxis # sin(w)û
    all(f1 .≈ 0.0) && return nothing
    k = (1/sum(f1.^2))
    f2 = (k - sqrt(k*k-k)) .* f1 # (1-cos(w))û/sin(w)
    vmod = f1 × view(params, 7:9)
    params[4:6] .= target
    params[7:9] .+= vmod .+ f2 × vmod
    return nothing
end

"""
    set_camera(distance, rotation, tilt; focus=zeros(3), twist=0.0)

Return a vector with the 9 camera parameters (camera position, focus and "up" vector),
given the distance from the center of the scene, rotation and tilt angles
in degrees, and optionally (as keyword arguments) one focus point of the
line of sight and the twist angle in degrees.

If `rotation`, `tilt` and `twist` angles are zero, the direction of the line of sight
is the Y axis, and the camera is positioned in the negative direction of that
axis.
"""
function set_camera(distance, rotation, tilt;
    focus = zeros(3), twist = 0.0, kwargs...)

    camera_position = [distance * cosd(tilt) * sind(rotation),
                      -distance * cosd(tilt) * cosd(rotation),
                       distance * sind(tilt)]
    camera_direction = normalize(camera_position .- focus)
    up_vector = [sind(twist), 0, cosd(twist)]
    _tilt!(up_vector, tilt)
    _rotate!(up_vector, rotation)
    parameters = [camera_position..., focus..., up_vector...]
    if any(focus .== 0.0)
        _focus!(parameters, focus)
    end
    return parameters
end

####################
## `draw` methods ##
####################

# Convex hull using Graham's scan - to calculate the frame of 3d plots
function graham_hull(x, y)
    @assert length(x) == length(y) > 2 "lengths of x and y must be equal and greater than 2"
    lowest = argmin(y)
    dx = x .- x[lowest]
    dy = y .- y[lowest]
    # sort by angle
    inds = sortperm(dx ./ sqrt.(dx.^2 .+ dy.^2), rev=true)
    hullx = x[inds[1:2]]
    hully = y[inds[1:2]]
    for i in inds[3:end]
        xn, yn = x[i], y[i] # next
        turnleft = false
        while !turnleft # try until there is a left turn
            dxn = xn - hullx[end]
            dyn = yn - hully[end]
            dxp = hullx[end] - hullx[end-1]
            dyp = hully[end] - hully[end-1]
            turnleft = (dxp*dyn - dxn*dyp >= 0)
            if turnleft # add the point
                push!(hullx, xn)
                push!(hully, yn)
            else # remove the previous point
                pop!(hullx)
                pop!(hully)
            end
        end
    end
    return hullx, hully
end

# Convex frame of 3d axes
function axes3frame(ax::Axes)
    xcorners = zeros(8)
    ycorners = zeros(8)
    i=0
    for xi=(0,1), yi=(0,1), zi=(0,1)
        i = xi*4 + yi*2 + zi + 1
        wc = GR.wc3towc(xi, yi, zi)
        xcorners[i], ycorners[i] = GR.wctondc(wc[1], wc[2])
    end
    hullx, hully = graham_hull(xcorners, ycorners)
end

function fillaxesbackground(ax)
    GR.savestate()
    GR.selntran(0)
    GR.setfillintstyle(GR.INTSTYLE_SOLID)
    GR.setfillcolorind(0)
    if ax.kind == :axes3d
        GR.fillarea(axes3frame(ax)...)
    else
        GR.fillrect(GR.inqviewport()...)
    end
    GR.selntran(1)
    GR.restorestate()
end

function draw(ax::Axes, background=true)
    # Special draw functions for polar axes and gr3
    ax.kind == :polar && return draw_polaraxes(ax, background)
    if ax.kind == :axes3d
        get(ax.options, :gr3, 0) ≠ 0 && return draw_gr3axes(ax)
        GR.setwindow(0, 1, 0, 1)
        GR.setspace(0, 1, ax.perspective...)
    end
    background && fillaxesbackground(ax)
    # Set the window of data seen
    GR.setwindow(ax.ranges[:x]..., ax.ranges[:y]...)
    # Modify scale (log or flipped axes)
    GR.setscale(ax.options[:scale])
    # Set the specifications of guides (grid and ticks)
    GR.setlinecolorind(1)
    GR.setlinewidth(1)
    ticksize, charheight = _tickcharheight()
    haskey(ax.options, :tickdir) && (ticksize *= ax.options[:tickdir])
    GR.setcharheight(charheight)
    xtick, xorg, majorx = ax.tickdata[:x]
    ytick, yorg, majory = ax.tickdata[:y]
    # Modify if xlog, etc.
    if ax.options[:scale] & GR.OPTION_X_LOG != 0
        xtick = 10
        majorx = 1
    end
    if ax.options[:scale] & GR.OPTION_Y_LOG != 0
        ytick = 10
        majory = 1
    end
    # Branching for different kinds of axes
    if ax.kind == :axes3d
        GR.setspace(ax.ranges[:z]..., ax.perspective...)
        ztick, zorg, majorz = ax.tickdata[:z]
        if ax.options[:scale] & GR.OPTION_Z_LOG != 0
            ztick = 10
            majorz = 1
        end
        # Draw ticks as 2-D if it is the XY plane
        if ax.perspective == [0, 90]
            (ax.options[:grid] ≠ 0) && GR.grid(xtick, ytick, 0, 0, majorx, majory)
            GR.axes(xtick, ytick, xorg[1], yorg[1], majorx, majory, ticksize)
            GR.axes(xtick, ytick, xorg[2], yorg[2], -majorx, -majory, -ticksize)
        else
            if (ax.options[:grid] ≠ 0)
                GR.grid3d(xtick, 0, ztick, xorg[1], yorg[2], zorg[1], 2, 0, 2)
                GR.grid3d(0, ytick, 0, xorg[1], yorg[2], zorg[1], 0, 2, 0)
            end
            GR.axes3d(xtick, 0, ztick, xorg[1], yorg[1], zorg[1], majorx, 0, majorz, -ticksize)
            GR.axes3d(0, ytick, 0, xorg[2], yorg[1], zorg[1], 0, majory, 0, ticksize)
        end
    elseif ax.kind == :axes2d
        (ax.options[:grid] ≠ 0) && GR.grid(xtick, ytick, 0, 0, majorx, majory)
        if !isempty(ax.ticklabels)
            fx = get(ax.ticklabels, :x, ticklabel_fun(identity))
            fy = get(ax.ticklabels, :y, ticklabel_fun(identity))
            GR.axeslbl(xtick, ytick, xorg[1], yorg[1], majorx, majory, ticksize, fx, fy)
        else
            GR.axes(xtick, ytick, xorg[1], yorg[1], majorx, majory, ticksize)
        end
        GR.axes(xtick, ytick, xorg[2], yorg[2], -majorx, -majory, -ticksize)
    end
    return nothing
end

function draw_polaraxes(ax, background=true)
    # Set the window as the unit circle
    GR.setwindow(-1.0, 1.0, -1.0, 1.0)
    # Modify scale (log or flipped axes)
    GR.setscale(ax.options[:scale]) # ??
    # Set the specifications of guides (grid and ticks)
    GR.savestate()
    GR.setlinewidth(1)
    _, charheight = _tickcharheight()
    GR.setcharheight(charheight)
    GR.setlinetype(GR.LINETYPE_SOLID)
    rmax = maximum(abs.(ax.ranges[:y]))
    tick = 0.5 * GR.tick(0.0, rmax)
    if background
        # Fill with background color
        GR.setfillintstyle(GR.INTSTYLE_SOLID)
        GR.setfillcolorind(0)
        GR.fillarc(-1, 1, -1, 1, 0, 359)
    end
    # Draw the arcs and radii
    n = round(Int, rmax / tick + 0.5)
    for i in 0:n
        r = float(i) / n
        if i % 2 == 0
            GR.setlinecolorind(88)
            if i > 0
                GR.drawarc(-r, r, -r, r, 0, 359)
            end
            GR.settextalign(GR.TEXT_HALIGN_LEFT, GR.TEXT_VALIGN_HALF)
            x, y = GR.wctondc(0.05, r)
            rounded = round(i * tick, sigdigits = 12, base = 10)
            GR.text(x, y, string(rounded))
        else
            GR.setlinecolorind(90)
            GR.drawarc(-r, r, -r, r, 0, 359)
        end
    end
    if get(ax.options, :radians, 1) == 0
        labels = (("$θ^o" for θ in 0:45:315)...,)
    else
        labels = ("0", "\\frac{\\pi}{4}", "\\frac{\\pi}{2}",
            "\\frac{3}{4}\\pi", "\\pi", "\\frac{5}{4}\\pi",
            "\\frac{3}{2}\\pi", "\\frac{7}{4}\\pi")
    end
    for i = 0:7
        sinf = sin(i * 0.25π)
        cosf = cos(i * 0.25π)
        GR.polyline([cosf, 0], [sinf, 0])
        GR.settextalign(GR.TEXT_HALIGN_CENTER, GR.TEXT_VALIGN_HALF)
        x, y = GR.wctondc(1.1 * cosf, 1.1 * sinf)
        GR.mathtex(x, y, labels[i + 1])
    end
    GR.restorestate()
    return nothing
end

draw_gr3axes(ax) = GR.gr3.cameralookat(ax.camera...)

"""
    _tickcharheight([vp])

Return the size of the tick marks and the character height,
proportional to the size of the rectangle defined by `vp`
(by default the current viewport).
"""
function _tickcharheight(vp=GR.inqviewport())
    diag = sqrt((vp[2] - vp[1])^2 + (vp[4] - vp[3])^2)
    ticksize = 0.0075 * diag
    charheight = max(0.018 * diag, 0.012)
    ticksize, charheight
end

"""
    charheight([vp])

Return the character height proportional to the size of the
rectangle defined by `vp` (by default the current viewport)
"""
charheight(vp...) = _tickcharheight(vp...)[2]


text:
include("calc_h.jl")
include("calc_i_m.jl")
include("calc_n_20.jl")
include("calc_nu_j.jl")
include("calc_t_k.jl")

include("calc_q_star.jl")
include("calc_p_w.jl")
include("calc_h_para.jl")
include("calc_cost.jl")
include("calc_b_c.jl")

include("solve.jl")
include("update.jl")
include("is_valid.jl")

include("solve_steady.jl")
include("solve_pulsed.jl")

export solve!


text:
using OECD_ICIO
using Test

@testset "OECD_ICIO.jl" begin
    # Write your tests here.
end


text:
using Distributions, StatsBase

function plothomingsearching(todo, level)
    @pgf p = GroupPlot(
                       {
                        group_style =
                        {
                         group_size = "2 by 1", yticklabels_at="edge left"
                        },
                        legend_style={legend_columns=-1},
                        legend_to_name={"CommonLegend"},
                        xmajorgrids,
                        ymajorgrids,
                        # "ytick distance=25",
                        # "xtick distance=25",
                        # scale_only_axis,
                        xlabel = "X (cm)",
                        axis_equal_image,
                        ymin = -180,
                        ymax = 130,
                        xmin = -180,
                        xmax = 180,
                       },
                      );
    for x in values(todo)
        for r in x.runs
            rotate!(r.data)
            center2!(r.data, r.data.originalnest)
        end
    end
    @pgf push!(p, {ylabel = "Y (cm)", xmin = -80, xmax = 80}, Plot({only_marks, color = "black"}, Coordinates([0], [0])));
    push!(p, LegendEntry("nest"))
    for (x, color) in zip(values(todo), ("magenta", "green", "blue", "cyan"))
        r = x.runs[1]
        @pgf push!(p, Plot({no_markers, color = color*"!50"}, Coordinates(Tuple.(r.data.homing))))
        push!(p, LegendEntry(get(x.runs[1].metadata.setup, level, "")))
    end
    for (x, color) in zip(values(todo), ("magenta", "green", "blue", "cyan"))
        for r in x.runs[2:end]
            @pgf push!(p, Plot({no_markers, color = color*"!50"}, Coordinates(Tuple.(r.data.homing))))
        end
    end
    labels = Dict{Tuple{Float64, Float64}, Int}()
    for (x, color) in zip(values(todo), ("magenta", "green", "blue", "cyan"))
        sina, cosa = mean(x.runs) do r
            xdata = first.(r.data.homing)
            xdata .-= mean(xdata)
            ydata = last.(r.data.homing)
            n = length(xdata)
            a, _ = hcat(xdata, ones(n))\ydata
            α = atan(a)
            if α < 0
                α += π
            end
            sina, cosa = sincos(α)
            SVector(sina, cosa)
        end
        α = atan(sina, cosa)
        ad = 90 - round(Int, rad2deg(α))
        a = tan(α)
        μ = mean(xy for r in x.runs for xy in r.data.homing)
        b = μ[2] - a*μ[1]
        m,M = extrema(yy for r in x.runs for (_, yy) in r.data.homing)
        # m -= 10
        # M += 10
        m = (m - b)/a
        M = (M - b)/a
        xl = [m, M]
        yl = a*xl .+ b
        labels[(xl[2], yl[2])] = ad
        @pgf push!(p, Plot({no_markers, color = color*"!75!black", ultra_thick, "->", ">=stealth"}, Coordinates(xl, yl)), )
    end
    for (k,v) in labels
        x,y = k
        anchor = v < 0 ? "south west" : "south east"
        push!(p, "\\node[inner sep = 0, fill=white, fill opacity = 0.5, text opacity=1, anchor = $anchor] at ($x,$y) {\$$v^{\\circ}\$};")
    end
    for x in values(todo)
        for r in x.runs
            rotate!(r.data)
            center2!(r.data, r.data.originalnest)
        end
    end
    @pgf push!(p, {}, Plot({only_marks, color = "black"}, Coordinates([0], [0])));
    push!(p, LegendEntry("nest"))
    for (x, color) in zip(values(todo), ("magenta", "green", "blue", "cyan"))
        r = x.runs[1]
        @pgf push!(p, Plot({no_markers, color = color*"!50"}, Coordinates(Tuple.(r.data.searching))))
        push!(p, LegendEntry(get(x.runs[1].metadata.setup, level, "")))
    end
    for (x, color) in zip(values(todo), ("magenta", "green", "blue", "cyan"))
        for r in x.runs[2:end]
            @pgf push!(p, Plot({no_markers, color = color*"!50"}, Coordinates(Tuple.(r.data.searching))))
        end
    end
    for (x, color) in zip(values(todo), ("magenta", "green", "blue", "cyan"))
        xy = reduce(vcat, r.data.searching for r in x.runs)
        xy2 = reduce(hcat, xy)
        d = fit_mle(MvNormal, xy2)
        μx, μy = mean(d)
        σ = var(d)
        fwhm = 2sqrt(2log(2)).*sqrt.(σ)
        rx, ry = fwhm/2
        t = range(0, 2π, length = 100)
        @pgf push!(p, Plot({fill = color, draw = "none", fill_opacity = 0.25}, Coordinates(μx .+ rx.*cos.(t), μy .+ ry.*sin.(t))), Plot({only_marks, color = color, draw = color*"!80!black"}, Coordinates([μx], [μy])))
    end
    p = TikzPicture(p);
    push!(p, "\\path (group c1r1.north west) -- node[above]{\\ref{CommonLegend}} (group c2r1.north east);");
    p = push_preamble!(TikzDocument(p), raw"""
                       \def\mystrut{\vphantom{hg}}
                       %\def\mystrut{\strut}
                       \pgfplotsset{
                                    legend style={font=\mystrut},
                                    y tick label style={font=\mystrut}
                                   }
                       """
                      )
    pgfsave(string("homing-searching.pdf"), p)
    # pgfsave(string(join(keys(todo), " "), "homing-searching.pdf"), p)
end


text:
"""
    UncertainDataset

Generic dataset containing uncertain values.

## Fields
- **`values::AbstractVector{<:AbstractUncertainValue}`**: The uncertain values.
"""
struct UncertainDataset <: AbstractUncertainValueDataset
    values::AbstractVector{<:AbstractUncertainValue}
end

"""
    ConstrainedUncertainDataset

Generic constrained dataset containing uncertain values.

## Fields
- **`values::AbstractVector{<:AbstractUncertainValue}`**: The uncertain values.
"""
struct ConstrainedUncertainDataset <: AbstractUncertainValueDataset
    values::AbstractVector{<:AbstractUncertainValue}
end

UncertainDataset(uv::T) where {T <:AbstractUncertainValue} =
    UncertainDataset([uv])
ConstrainedUncertainDataset(uv::T) where {T<:AbstractUncertainValue} =
    ConstrainedUncertainDataset([uv])

##########################
# Sorting
#########################

export
UncertainDataset,
ConstrainedUncertainDataset


text:
## Core compute function

function eqws_compute(image::Array{Float64,2}, filter_hash::Dict, filter_hash2::Dict=filter_hash;
    doS2::Bool=true, doS20::Bool=false, norm=true, iso=false, FFTthreads=2, normS1::Bool=false, normS1iso::Bool=false)
    # image        - input for WST
    # filter_hash  - filter hash from fink_filter_hash
    # filter_hash2 - filters for second order.  Default to same as first order.
    # doS2         - compute S2 coeffs
    # doS20        - compute S2 coeffs
    # norm         - scale to mean zero, unit variance
    # iso          - sum over angles to obtain isotropic coeffs

    # Use 2 threads for FFT
    FFTW.set_num_threads(FFTthreads)

    # array sizes
    (Nx, Ny)  = size(image)
    if Nx != Ny error("Input image must be square") end
    (Nf, )    = size(filter_hash["filt_index"])
    if Nf == 0  error("filter hash corrupted") end
    @assert Nx==filter_hash["npix"] "Filter size should match npix"
    @assert Nx==filter_hash2["npix"] "Filter2 size should match npix"
    @assert (normS1 && normS1iso) != 1 "normS1 and normS1iso are mutually exclusive"

    # allocate coeff arrays
    out_coeff = []
    S0  = zeros(Float64, 2)
    S1  = zeros(Float64, Nf)
    if doS2  S2  = zeros(Float64, Nf, Nf) end  # traditional 2nd order
    if doS20 S20 = zeros(Float64, Nf, Nf) end  # real space correlation
    anyM2 = doS2 | doS20
    anyrd = doS2 | doS20             # compute real domain with iFFT

    # allocate image arrays for internal use
    if anyrd im_rd_0_1  = Array{Float64, 3}(undef, Nx, Ny, Nf) end

    ## 0th Order
    S0[1]   = mean(image)
    norm_im = image.-S0[1]
    S0[2]   = sum(norm_im .* norm_im)/(Nx*Ny)
    if norm
        norm_im ./= sqrt(Nx*Ny*S0[2])
    else
        norm_im = copy(image)
    end

    append!(out_coeff,S0[:])

    ## 1st Order
    im_fd_0 = fft(norm_im)  # 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"]

    zarr = zeros(ComplexF64, Nx, Ny)  # temporary array to fill with zvals

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

    ## Main 1st Order and Precompute 2nd Order
    for f = 1:Nf
        S1tot = 0.0
        f_i = f_ind[f]  # CartesianIndex list for filter
        f_v = f_val[f]  # Values for f_i
        # for (ind, val) in zip(f_i, f_v)   # this is slower!
        if length(f_i) > 0
            for i = 1:length(f_i)
                ind       = f_i[i]
                zval      = f_v[i] * im_fd_0[ind]
                S1tot    += abs2(zval)
                zarr[ind] = zval        # filter*image in Fourier domain
            end
            S1[f] = S1tot/(Nx*Ny)  # image power
            if anyrd
                im_rd_0_1[:,:,f] .= abs2.(P*zarr) end
            zarr[f_i] .= 0
        end
    end

    append!(out_coeff, iso ? filter_hash["S1_iso_mat"]*S1 : S1)

    if normS1iso
        S1iso = vec(reshape(filter_hash["S1_iso_mat"]*S1,(1,:))*filter_hash["S1_iso_mat"])
    end

    # we stored the abs()^2, so take sqrt (this is faster to do all at once)
    if anyrd im_rd_0_1 .= sqrt.(im_rd_0_1) end

    Mat2 = filter_hash["S2_iso_mat"]
    if doS2
        f_ind2   = filter_hash2["filt_index"]  # (J, L) array of filters represented as index value pairs
        f_val2   = filter_hash2["filt_value"]

        ## Traditional second order
        for f1 = 1:Nf
            thisim = fft(im_rd_0_1[:,:,f1])  # Could try rfft here
            # Loop over f2 and do second-order convolution
            if normS1
                normS1pwr = S1[f1]
            elseif normS1iso
                normS1pwr = S1iso[f1]
            else
                normS1pwr = 1
            end
            for f2 = 1:Nf
                f_i = f_ind2[f2]  # CartesianIndex list for filter
                f_v = f_val2[f2]  # Values for f_i
                # sum im^2 = sum(|fft|^2/npix)
                S2[f1,f2] = sum(abs2.(f_v .* thisim[f_i]))/(Nx*Ny)/normS1pwr
            end
        end
        append!(out_coeff, iso ? Mat2*S2[:] : S2[:])
    end

    # Real domain 2nd order
    if doS20
        Amat = reshape(im_rd_0_1, Nx*Ny, Nf)
        S20  = Amat' * Amat
        append!(out_coeff, iso ? Mat2*S20[:] : S20[:])
    end

    return out_coeff
end

function eqws_compute_convmap(image::Array{Float64,2}, filter_hash::Dict, filter_hash2::Dict=filter_hash;
    doS2::Bool=true, doS20::Bool=false, norm=true, iso=false, FFTthreads=2, normS1::Bool=false, normS1iso::Bool=false,
    p=2)
    # image        - input for WST
    # filter_hash  - filter hash from fink_filter_hash
    # filter_hash2 - filters for second order.  Default to same as first order.
    # doS2         - compute S2 coeffs
    # doS20        - compute S2 coeffs
    # norm         - scale to mean zero, unit variance
    # iso          - sum over angles to obtain isotropic coeffs

    # Use 2 threads for FFT
    FFTW.set_num_threads(FFTthreads)

    # array sizes
    (Nx, Ny)  = size(image)
    if Nx != Ny error("Input image must be square") end
    (Nf, )    = size(filter_hash["filt_index"])
    if Nf == 0  error("filter hash corrupted") end
    @assert Nx==filter_hash["npix"] "Filter size should match npix"
    @assert Nx==filter_hash2["npix"] "Filter2 size should match npix"
    @assert (normS1 && normS1iso) != 1 "normS1 and normS1iso are mutually exclusive"

    # allocate coeff arrays
    out_coeff = []
    S0  = zeros(Float64, 2)
    S1  = zeros(Float64, Nf)
    if doS2  S2  = zeros(Float64, Nf, Nf) end  # traditional 2nd order
    if doS20 S20 = zeros(Float64, Nf, Nf) end  # real space correlation
    anyM2 = doS2 | doS20

    # allocate image arrays for internal use
    im_rd_0_1  = Array{Float64, 3}(undef, Nx, Ny, Nf)
    im_rd_0_2  = Array{Float64, 4}(undef, Nx, Ny, Nf, Nf)

    ## 0th Order
    S0[1]   = mean(image)
    norm_im = image.-S0[1]
    S0[2]   = sum(abs.(norm_im).^p)/(Nx*Ny)
    if norm
        norm_im ./= (Nx*Ny*S0[2]).^(1/p)
    else
        norm_im = copy(image)
    end

    append!(out_coeff,S0[:])

    ## 1st Order
    im_fd_0 = fft(norm_im)  # 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"]

    zarr = zeros(ComplexF64, Nx, Ny)  # temporary array to fill with zvals

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

    ## Main 1st Order and Precompute 2nd Order
    for f = 1:Nf
        f_i = f_ind[f]  # CartesianIndex list for filter
        f_v = f_val[f]  # Values for f_i
        if length(f_i) > 0
            zarr[f_i] = f_v.*im_fd_0[f_i]
            im_rd_0_1[:,:,f] .= abs.(P*zarr)
            S1[f] = sum(im_rd_0_1[:,:,f].^(p))  # image power
            #zero out the intermediate arrays
            zarr[f_i] .= 0
        end
    end

    #append!(out_coeff, iso ? filter_hash["S1_iso_mat"].*S1 : S1)
    append!(out_coeff, [im_rd_0_1[:,:,i].^(p)  for i in 1:Nf])

    if normS1iso
        S1iso = vec(reshape(filter_hash["S1_iso_mat"]*S1,(1,:))*filter_hash["S1_iso_mat"])
    end

    Mat2 = filter_hash["S2_iso_mat"]
    if doS2
        f_ind2   = filter_hash2["filt_index"]  # (J, L) array of filters represented as index value pairs
        f_val2   = filter_hash2["filt_value"]

        ## Traditional second order
        for f1 = 1:Nf
            thisim = fft(im_rd_0_1[:,:,f1])  # Could try rfft here
            # Loop over f2 and do second-order convolution
            if normS1
                normS1pwr = S1[f1]
            elseif normS1iso
                normS1pwr = S1iso[f1]
            else
                normS1pwr = 1
            end
            for f2 = 1:Nf
                f_i = f_ind2[f2]  # CartesianIndex list for filter
                f_v = f_val2[f2]  # Values for f_i
                zarr[f_i] = f_v .* thisim[f_i]
                im_rd_0_2[:,:,f1,f2] .= abs.(P*zarr).^(p)
                S2[f1,f2] = sum(im_rd_0_2[:,:,f1,f2])/normS1pwr
                #zero out the intermediate arrays
                zarr[f_i] .= 0
            end
        end
        #append!(out_coeff, iso ? Mat2*S2[:] : S2[:])
        append!(out_coeff, [im_rd_0_2[:,:,f1,f2] for f2 in 1:Nf for f1 in 1:Nf])
    end

    # Real domain 2nd order
    if doS20
        Amat = reshape(im_rd_0_1, Nx*Ny, Nf)
        S20  = Amat' * Amat
        append!(out_coeff, iso ? Mat2*S20[:] : S20[:])
    end

    return out_coeff
end

function eqws_compute_p(image::Array{Float64,2}, filter_hash::Dict, filter_hash2::Dict=filter_hash;
    doS2::Bool=true, doS20::Bool=false, norm=true, iso=false, FFTthreads=2, normS1::Bool=false, normS1iso::Bool=false,
    p=2)
    # image        - input for WST
    # filter_hash  - filter hash from fink_filter_hash
    # filter_hash2 - filters for second order.  Default to same as first order.
    # doS2         - compute S2 coeffs
    # doS20        - compute S2 coeffs
    # norm         - scale to mean zero, unit variance
    # iso          - sum over angles to obtain isotropic coeffs

    # Use 2 threads for FFT
    FFTW.set_num_threads(FFTthreads)

    # array sizes
    (Nx, Ny)  = size(image)
    if Nx != Ny error("Input image must be square") end
    (Nf, )    = size(filter_hash["filt_index"])
    if Nf == 0  error("filter hash corrupted") end
    @assert Nx==filter_hash["npix"] "Filter size should match npix"
    @assert Nx==filter_hash2["npix"] "Filter2 size should match npix"
    @assert (normS1 && normS1iso) != 1 "normS1 and normS1iso are mutually exclusive"

    # allocate coeff arrays
    out_coeff = []
    S0  = zeros(Float64, 2)
    S1  = zeros(Float64, Nf)
    if doS2  S2  = zeros(Float64, Nf, Nf) end  # traditional 2nd order
    if doS20 S20 = zeros(Float64, Nf, Nf) end  # real space correlation
    anyM2 = doS2 | doS20

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

    ## 0th Order
    S0[1]   = mean(image)
    norm_im = image.-S0[1]
    S0[2]   = sum(abs.(norm_im).^p)/(Nx*Ny)
    if norm
        norm_im ./= (Nx*Ny*S0[2]).^(1/p)
    else
        norm_im = copy(image)
    end

    append!(out_coeff,S0[:])

    ## 1st Order
    im_fd_0 = fft(norm_im)  # 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"]

    zarr = zeros(ComplexF64, Nx, Ny)  # temporary array to fill with zvals

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

    ## Main 1st Order and Precompute 2nd Order
    for f = 1:Nf
        f_i = f_ind[f]  # CartesianIndex list for filter
        f_v = f_val[f]  # Values for f_i
        if length(f_i) > 0
            zarr[f_i] = f_v.*im_fd_0[f_i]
            im_rd_0_1[:,:,f] .= abs.(P*zarr)
            S1[f] = sum(im_rd_0_1[:,:,f].^(p))  # image power
            #zero out the intermediate arrays
            zarr[f_i] .= 0
        end
    end

    append!(out_coeff, iso ? filter_hash["S1_iso_mat"]*S1 : S1)

    if normS1iso
        S1iso = vec(reshape(filter_hash["S1_iso_mat"]*S1,(1,:))*filter_hash["S1_iso_mat"])
    end

    Mat2 = filter_hash["S2_iso_mat"]
    if doS2
        f_ind2   = filter_hash2["filt_index"]  # (J, L) array of filters represented as index value pairs
        f_val2   = filter_hash2["filt_value"]

        ## Traditional second order
        for f1 = 1:Nf
            thisim = fft(im_rd_0_1[:,:,f1])  # Could try rfft here
            # Loop over f2 and do second-order convolution
            if normS1
                normS1pwr = S1[f1]
            elseif normS1iso
                normS1pwr = S1iso[f1]
            else
                normS1pwr = 1
            end
            for f2 = 1:Nf
                f_i = f_ind2[f2]  # CartesianIndex list for filter
                f_v = f_val2[f2]  # Values for f_i
                zarr[f_i] = f_v .* thisim[f_i]
                S2[f1,f2] = sum(abs.(P*zarr).^(p))/normS1pwr
                #zero out the intermediate arrays
                zarr[f_i] .= 0
            end
        end
        append!(out_coeff, iso ? Mat2*S2[:] : S2[:])
    end

    # Real domain 2nd order
    if doS20
        Amat = reshape(im_rd_0_1, Nx*Ny, Nf)
        S20  = Amat' * Amat
        append!(out_coeff, iso ? Mat2*S20[:] : S20[:])
    end

    return out_coeff
end

function eqws_compute_RGB(image::Array{Float64}, filter_hash::Dict, filter_hash2::Dict=filter_hash;
    doS20::Bool=true, norm=true, iso=false, FFTthreads=2)
    # image        - input for WST
    # filter_hash  - filter hash from fink_filter_hash
    # filter_hash2 - filters for second order.  Default to same as first order.
    # doS2         - compute S2 coeffs
    # doS12        - compute S2 coeffs
    # doS20        - compute S2 coeffs
    # norm         - scale to mean zero, unit variance
    # iso          - sum over angles to obtain isotropic coeffs

    # Use 2 threads for FFT
    FFTW.set_num_threads(FFTthreads)

    # array sizes
    (Nx, Ny, Nc)  = size(image)
    if Nx != Ny error("Input image must be square") end
    (Nf, )    = size(filter_hash["filt_index"])
    if Nf == 0  error("filter hash corrupted") end
    @assert Nx==filter_hash["npix"] "Filter size should match npix"
    @assert Nx==filter_hash2["npix"] "Filter2 size should match npix"

    # allocate coeff arrays
    out_coeff = []
    S0  = zeros(Float64, Nc*2)
    S1  = zeros(Float64, Nc*Nf)

    if doS20 S20 = zeros(Float64, Nf, Nf) end  # real space correlatio
    # anyM2 = doS2 | doS12 | doS20
    anyrd = doS20 #| doS2             # compute real domain with iFFT

    # allocate image arrays for internal use
    mean_im = zeros(Float64,1,1,Nc)
    pwr_im = zeros(Float64,1,1,Nc)
    norm_im = zeros(Float64,Nx,Ny,Nc)
    im_fd_0 = zeros(ComplexF64, Nx, Ny, Nc)
    im_fd_0_sl = zeros(ComplexF64, Nx, Ny)

    if doS20
        Amat1 = zeros(Nx*Ny, Nf)
        Amat2 = zeros(Nx*Ny, Nf)
    end

    if anyrd im_rd_0_1  = Array{Float64, 4}(undef, Nx, Ny, Nf, Nc) end

    ## 0th Order
    mean_im = mean(image, dims=(1,2))
    S0[1:Nc]   = dropdims(mean_im,dims=(1,2))
    norm_im = image.-mean_im
    pwr_im = sum(norm_im .* norm_im,dims=(1,2))
    S0[1+Nc:end]   = dropdims(pwr_im,dims=(1,2))./(Nx*Ny)
    if norm
        norm_im ./= sqrt.(pwr_im)
    else
        norm_im = copy(image)
    end

    append!(out_coeff,S0[:])

    ## 1st Order
    im_fd_0 .= fft(norm_im,(1,2))  # 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"]

    zarr = zeros(ComplexF64, Nx, Ny)  # temporary array to fill with zvals

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

    ## Main 1st Order and Precompute 2nd Order
    for f = 1:Nf
        S1tot = 0.0
        f_i = f_ind[f]  # CartesianIndex list for filter
        f_v = f_val[f]  # Values for f_i
        # for (ind, val) in zip(f_i, f_v)   # this is slower!
        for chan = 1:Nc
            im_fd_0_sl .= im_fd_0[:,:,chan]
            if length(f_i) > 0
                for i = 1:length(f_i)
                    ind       = f_i[i]
                    zval      = f_v[i] * im_fd_0_sl[ind]
                    S1tot    += abs2(zval)
                    zarr[ind] = zval        # filter*image in Fourier domain
                end
                S1[f+(chan-1)*Nf] = S1tot/(Nx*Ny)  # image power
                if anyrd
                    im_rd_0_1[:,:,f,chan] .= abs2.(P*zarr) end
                zarr[f_i] .= 0
            end
        end
    end

    if iso
        S1M = filter_hash["S1_iso_mat"]
        M1 = blockdiag(S1M,S1M,S1M)
    end
    append!(out_coeff, iso ? M1*S1 : S1)

    # we stored the abs()^2, so take sqrt (this is faster to do all at once)
    if anyrd im_rd_0_1 .= sqrt.(im_rd_0_1) end

    # Real domain 2nd order
    if doS20
        for chan1 = 1:Nc
            for chan2 = 1:Nc
                Amat1 = reshape(im_rd_0_1[:,:,:,chan1], Nx*Ny, Nf)
                Amat2 = reshape(im_rd_0_1[:,:,:,chan2], Nx*Ny, Nf)
                S20  = Amat1' * Amat2
                append!(out_coeff, iso ? filter_hash["S2_iso_mat"]*S20[:] : S20[:])
            end
        end
    end

    return out_coeff
end

function eqws_compute_wrapper(image::Array{Float64,2}, filter_hash::Dict;
    doS2::Bool=true, doS20::Bool=false, apodize=false, norm=false, iso=false, FFTthreads=1, filter_hash2::Dict=filter_hash, coeff_mask=nothing)
    Nf = size(filter_hash["filt_index"])[1]
    #Not using coeff_mask here after all. ASSUMES ANY ONE of doS12, doS2 and doS20 are true, when using coeff_mask
    #@assert !iso "Iso not implemented yet"

    if apodize
        ap_img = apodizer(image)
    else
        ap_img = image
    end

    sall = eqws_compute(ap_img, filter_hash, filter_hash2, doS2=doS2, doS20=doS20, norm=norm, iso=iso, FFTthreads=FFTthreads)
    if coeff_mask!=nothing
        @assert length(coeff_mask)==length(sall) "The length of the coeff_mask must match the length of the output coefficients"
        return sall[coeff_mask]
    else
        return sall
    end
end

function eqws_compute_gpu_fake(image::Array{Float64,2}, filter_hash::Dict, filter_hash2::Dict=filter_hash;
    doS2::Bool=true, doS12::Bool=false, doS20::Bool=false, norm=true, iso=false,batch_mode::Bool=false,opt_memory::Bool=false,max_memory::Int64=0)

    @assert !iso "Isotropic is not implemented on GPU"
    @assert !doS20 "S20 is not supported on GPU"
    # image        - input for WST
    # filter_hash  - filter hash from fink_filter_hash
    # filter_hash2 - filters for second order.  Default to same as first order.
    # doS2         - compute S2 coeffs Similar to Mallat(2014), Abs in real space
    # doS12        - compute S2 coeffs Abs in Fourier space
    # doS20        - compute S2 coeffs Square in real space
    # norm         - scale to mean zero, unit variance
    # iso          - sum over angles to obtain isotropic coeffs

    # array sizes
    (Nx, Ny)  = size(image)
    if Nx != Ny error("Input image must be square") end
    (Nf, )    = size(filter_hash["filt_vals"])
    if Nf == 0  error("filter hash corrupted") end
    @assert Nx==filter_hash["npix"] "Filter size should match npix"
    @assert Nx==filter_hash2["npix"] "Filter2 size should match npix"

    # allocate coeff arrays
    out_coeff = []
    S0  = zeros(Float64, 2)
    S1  = zeros(Float64, Nf)
    if doS2  S2  = zeros(Float64, Nf, Nf) end  # traditional 2nd order
    if doS12 S12 = zeros(Float64, Nf, Nf) end  # Fourier correlation
    if doS20 S20 = zeros(Float64, Nf, Nf) end  # real space correlation
    anyM2 = doS2 | doS12 | doS20
    anyrd = doS2 | doS20             # compute real domain with iFFT

    # allocate image arrays for internal use
    if doS12 im_fdf_0_1 = zeros(Float64,           Nx, Ny) end   # this must be zeroed!
    if anyrd im_rd_0_1  = Array{Float64, 2}(undef, Nx, Ny) end

    image= image
    ## 0th Order
    S0[1]   = mean(image)
    norm_im = image.-S0[1]
    S0[2]   = sum(norm_im .* norm_im)/(Nx*Ny)
    if norm
        norm_im./= sqrt(Nx*Ny*S0[2])
    else
        norm_im = image
    end

    ## 1st Order
    #REview: We fft shift!
    im_fd_0 = FFTW.fftshift(FFTW.fft(norm_im))  # total power=1.0

    # unpack filter_hash
    f_mms   = filter_hash["filt_mms"]  # (J, L) array of filters represented as index value pairs
    f_vals   = filter_hash["filt_vals"]

    zarr = zeros(ComplexF64, Nx, Ny)  # temporary array to fill with zvals

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

    ## Main 1st Order and Precompute 2nd Order
    for f1 = 1:Nf
        f_mm1 = f_mms[:,:,f1]  # Min,Max
        f_val1 = f_vals[f1]  # values

        #Review: Removes length check-> empty filter case???, I think this should be prevented beforehands
        zv=f_val1.*im_fd_0[f_mm1[1,1]:f_mm1[1,2],f_mm1[2,1]:f_mm1[2,2]] #dense multiplication for the filter region
        zarr[f_mm1[1,1]:f_mm1[1,2],f_mm1[2,1]:f_mm1[2,2]]=zv

        if doS12 im_fdf_0_1[f_mm1[1,1]:f_mm1[1,2],f_mm1[2,1]:f_mm1[2,2]] .= abs.(zv) end

        S1[f1] = sum(abs2.(zv))/(Nx*Ny)  # image power
        if anyrd im_rd_0_1.= abs2.(P*FFTW.ifftshift(zarr)) end
        zarr[f_mm1[1,1]:f_mm1[1,2],f_mm1[2,1]:f_mm1[2,2]] .= 0

        im_rd_0_1 .= sqrt.(im_rd_0_1)
        im_fdf_1_1 = FFTW.fftshift(FFTW.fft(im_rd_0_1))
        for f2 = 1:Nf
            f_mm2 = f_mms[:,:,f2]
            f_val2 = f_vals[f2]

            S2[f1,f2] = sum(abs2.(f_val2 .* im_fdf_1_1[f_mm2[1,1]:f_mm2[1,2],f_mm2[2,1]:f_mm2[2,2]]))/(Nx*Ny)
            if !doS12
                continue
            end

            if f2>f1
                S12[f1,f2]=NaN
                continue
            end

            if f_mm1[1,1]>f_mm2[1,2] || f_mm2[1,1]>f_mm1[1,2] || f_mm1[2,1]>f_mm2[2,2] || f_mm2[2,1]>f_mm1[2,2]
                S12[f1,f2]=NaN
                continue
            end

            mx=max(f_mm1[1,1],f_mm2[1,1])
            my=max(f_mm1[2,1],f_mm2[2,1])
            Mx=min(f_mm1[1,2],f_mm2[1,2])
            My=min(f_mm1[2,2],f_mm2[2,2])
            sfx=Mx-mx
            sfy=My-my
            st1x=mx-f_mm1[1,1]+1
            st1y=my-f_mm1[2,1]+1
            st2x=mx-f_mm2[1,1]+1
            st2y=my-f_mm2[2,1]+1

            S12[f1,f2]=sum(f_val1[st1x:st1x+sfx,st1y:st1y+sfy].*f_val2[st2x:st2x+sfx,st2y:st2y+sfy].*abs2.(im_fd_0[mx:Mx,my:My]))
        end
    end
    append!(out_coeff,S0)
    append!(out_coeff, iso ? filter_hash["S1_iso_mat"]*S1 : S1)
    if doS2 append!(out_coeff, iso ? filter_hash["S1_iso_mat"]*S2 : S2) end
    if doS12 append!(out_coeff, iso ? filter_hash["S1_iso_mat"]*S12 : S12) end
    return out_coeff
end

function eqws_compute_gpu(image::Array{Float64,2}, filter_hash::Dict, filter_hash2::Dict=filter_hash;
    doS2::Bool=true, doS12::Bool=false, doS20::Bool=false, norm=true, iso=false,batch_mode::Bool=false,opt_memory::Bool=false,max_memory::Int64=0)

    @assert !iso "Isotropic is not implemented on GPU"
    @assert !doS20 "S20 is not supported on GPU"
    # image        - input for WST
    # filter_hash  - filter hash from fink_filter_hash
    # filter_hash2 - filters for second order.  Default to same as first order.
    # doS2         - compute S2 coeffs Similar to Mallat(2014), Abs in real space
    # doS12        - compute S2 coeffs Abs in Fourier space
    # doS20        - compute S2 coeffs Square in real space
    # norm         - scale to mean zero, unit variance
    # iso          - sum over angles to obtain isotropic coeffs

    # array sizes
    (Nx, Ny)  = size(image)
    if Nx != Ny error("Input image must be square") end
    (Nf, )    = size(filter_hash["filt_vals"])
    if Nf == 0  error("filter hash corrupted") end
    @assert Nx==filter_hash["npix"] "Filter size should match npix"
    @assert Nx==filter_hash2["npix"] "Filter2 size should match npix"

    # allocate coeff arrays
    out_coeff = []
    S0  = zeros(Float64, 2)
    S1  = zeros(Float64, Nf)
    if doS2  S2  = zeros(Float64, Nf, Nf) end  # traditional 2nd order
    if doS12 S12 = zeros(Float64, Nf, Nf) end  # Fourier correlation
    if doS20 S20 = zeros(Float64, Nf, Nf) end  # real space correlation
    anyM2 = doS2 | doS12 | doS20
    anyrd = doS2 | doS20             # compute real domain with iFFT

    # allocate image arrays for internal use
    if doS12 im_fdf_0_1 = CUDA.zeros(Float64, Nx, Ny) end   # this must be zeroed!
    if anyrd im_rd_0_1  = CUDA.CuArray(Array{Float64, 2}(undef, Nx, Ny)) end

    image= CUDA.CuArray(image)
    ## 0th Order
    S0[1]   = mean(image)
    norm_im = image.-S0[1]
    S0[2]   = sum(norm_im .* norm_im)/(Nx*Ny)
    if norm
        norm_im./= sqrt(Nx*Ny*S0[2])
    else
        norm_im = image
    end

    ## 1st Order
    #REview: We fft shift!
    im_fd_0 = CUDA.CUFFT.fftshift(CUDA.CUFFT.fft(norm_im))  # total power=1.0

    # unpack filter_hash
    f_mms   = filter_hash["filt_mms"]  # (J, L) array of filters represented as index value pairs
    f_vals   = filter_hash["filt_vals"]

    zarr = CUDA.zeros(ComplexF64, Nx, Ny)  # temporary array to fill with zvals

    # make a FFTW "plan" for an array of the given size and type
    if anyrd
        P = CUDA.CUFFT.plan_ifft(im_fd_0) end  # P is an operator, P*im is ifft(im)

    ## Main 1st Order and Precompute 2nd Order
    for f1 = 1:Nf
        f_mm1 = f_mms[:,:,f1]  # Min,Max
        f_val1 = f_vals[f1]  # values

        #Review: Removes length check-> empty filter case???, I think this should be prevented beforehands
        zv=f_val1.*im_fd_0[f_mm1[1,1]:f_mm1[1,2],f_mm1[2,1]:f_mm1[2,2]] #dense multiplication for the filter region
        zarr[f_mm1[1,1]:f_mm1[1,2],f_mm1[2,1]:f_mm1[2,2]]=zv

        if doS12 im_fdf_0_1[f_mm1[1,1]:f_mm1[1,2],f_mm1[2,1]:f_mm1[2,2]] .= abs.(zv) end

        S1[f1] = sum(abs2.(zv))/(Nx*Ny)  # image power
        if anyrd im_rd_0_1.= abs2.(P*CUDA.CUFFT.ifftshift(zarr)) end
        zarr[f_mm1[1,1]:f_mm1[1,2],f_mm1[2,1]:f_mm1[2,2]] .= 0

        im_rd_0_1 .= sqrt.(im_rd_0_1)
        im_fdf_1_1 = CUDA.CUFFT.fftshift(CUDA.CUFFT.fft(im_rd_0_1))
        for f2 = 1:Nf
            f_mm2 = f_mms[:,:,f2]
            f_val2 = f_vals[f2]

            S2[f1,f2] = sum(abs2.(f_val2 .* im_fdf_1_1[f_mm2[1,1]:f_mm2[1,2],f_mm2[2,1]:f_mm2[2,2]]))/(Nx*Ny)
            if !doS12
                continue
            end

            if f2>f1
                S12[f1,f2]=NaN
                continue
            end

            if f_mm1[1,1]>f_mm2[1,2] || f_mm2[1,1]>f_mm1[1,2] || f_mm1[2,1]>f_mm2[2,2] || f_mm2[2,1]>f_mm1[2,2]
                S12[f1,f2]=NaN
                continue
            end

            mx=max(f_mm1[1,1],f_mm2[1,1])
            my=max(f_mm1[2,1],f_mm2[2,1])
            Mx=min(f_mm1[1,2],f_mm2[1,2])
            My=min(f_mm1[2,2],f_mm2[2,2])
            sfx=Mx-mx
            sfy=My-my
            st1x=mx-f_mm1[1,1]+1
            st1y=my-f_mm1[2,1]+1
            st2x=mx-f_mm2[1,1]+1
            st2y=my-f_mm2[2,1]+1

            S12[f1,f2]=sum(f_val1[st1x:st1x+sfx,st1y:st1y+sfy].*f_val2[st2x:st2x+sfx,st2y:st2y+sfy].*abs2.(im_fd_0[mx:Mx,my:My]))
        end
    end
    append!(out_coeff,S0)
    append!(out_coeff, iso ? filter_hash["S1_iso_mat"]*S1 : S1)
    if doS2 append!(out_coeff, iso ? filter_hash["S1_iso_mat"]*S2 : S2) end
    if doS12 append!(out_coeff, iso ? filter_hash["S1_iso_mat"]*S12 : S12) end
    return out_coeff
end

function eqws_compute_batch_gpu(image::Union{Array{Float64,3},Array{Float32,3}}, filter_hash::Dict, filter_hash2::Dict=filter_hash;
    doS2::Bool=true, doS20::Bool=false, norm=true, iso=false, normS1::Bool=false, normS1iso::Bool=false,prec=Float64)
    # image        - input for WST
    # filter_hash  - filter hash from fink_filter_hash
    # filter_hash2 - filters for second order.  Default to same as first order.
    # doS2         - compute S2 coeffs
    # doS20        - compute S2 coeffs
    # norm         - scale to mean zero, unit variance
    # iso          - sum over angles to obtain isotropic coeffs

    # array sizes
    (Nx, Ny, Nb)  = size(image)
    if Nx != Ny error("Input image must be square") end
    Nf    = filter_hash["Nf"]
    if Nf == 0  error("filter hash corrupted") end
    @assert Nx==filter_hash["npix"] "Filter size should match npix"
    @assert Nx==filter_hash2["npix"] "Filter2 size should match npix"
    @assert (normS1 && normS1iso) != 1 "normS1 and normS1iso are mutually exclusive"
    if prec == Float64
        cprec = ComplexF64
    elseif prec == Float32
        cprec = ComplexF32
    else
        error("precision is not acceptable, please choose Float32 or Float64")
    end

    # allocate coeff arrays
    anyrd = doS2 | doS20             # compute real domain with iFFT

    if iso
        ind=1
    else
        ind=2
    end

    if !anyrd
        S1s = size(filter_hash["S1_iso_mat"])[ind]
        out_coeff = zeros(S1s+2,Nb)
    elseif doS2 ⊻ doS20
        S1s = size(filter_hash["S1_iso_mat"])[ind]
        S2s = size(filter_hash["S2_iso_mat"])[ind]
        out_coeff = zeros(S1s + S2s + 2,Nb)
    else
        S1s = size(filter_hash["S1_iso_mat"])[ind]
        S2s = size(filter_hash["S2_iso_mat"])[ind]
        out_coeff = zeros(S1s + 2*S2s + 2,Nb)
    end

    S0  = zeros(prec, 2, Nb)
    S1  = zeros(prec, Nf, Nb)
    if doS2  S2  = zeros(prec, Nf, Nf, Nb) end  # traditional 2nd order
    if doS20 S20 = zeros(prec, Nf, Nf, Nb) end  # real space correlation

    # allocate image arrays for internal use
    if anyrd im_rd_0_1  = CUDA.CuArray{prec}(undef, Nx, Ny, Nf, Nb) end

    #image handling
    if typeof(image[1]) != prec
        image = convert(Array{prec,3},image)
    end

    ## 0th Order
    μ = mean(image, dims=(1,2))
    S0[1,:] = dropdims(μ, dims=(1,2))
    norm_im = image.-μ
    σ       = sum(norm_im .* norm_im, dims=(1,2))/(Nx*Ny)
    S0[2,:] = dropdims(σ, dims=(1,2))
    if norm
        norm_im ./= sqrt.(Nx*Ny*σ)
    else
        norm_im = image
    end

    out_coeff[1:2,:] .= S0

    image_in = CUDA.CuArray(norm_im)

    if doS2
        P2 = CUDA.CUFFT.plan_fft(image_in,[1,2])
        im_fd_0 = P2*image_in
    else
        im_fd_0 = CUDA.CUFFT.fft(image_in,[1,2])  # total power=1.0
    end

    CUDA.unsafe_free!(image_in)

    # unpack filter_hash
    filts   = filter_hash["gpu_filt"]  # dense gpu array

    zarr = CUDA.zeros(cprec, Nx, Ny, Nb)  # temporary array to fill with zvals

    # make a FFTW "plan" for an array of the given size and type
    if anyrd
        P1 = CUDA.CUFFT.plan_ifft(im_fd_0,[1,2]) # P is an operator, P*im is ifft(im)
    end

    ## Main 1st Order and Precompute 2nd Order
    for f = 1:Nf
        zarr = filts[:,:,f].*im_fd_0

        if anyrd
            im_rd_0_1[:,:,f,:] = abs2.(P1*zarr)
            S1[f,:] = sum(im_rd_0_1[:,:,f,:],dims=(1,2))  # image power
        else
            S1[f,:] = sum(abs2.(zarr),dims=(1,2))/(Nx*Ny)
        end
    end

    CUDA.unsafe_free!(zarr)

    if iso
        out_coeff[3:3+S1s-1,:] .= filter_hash["S1_iso_mat"]*S1
    else
        out_coeff[3:3+S1s-1,:] .= S1
    end

    if normS1iso
        S1iso = vec(reshape(filter_hash["S1_iso_mat"]*S1,(1,:))*filter_hash["S1_iso_mat"]) #FIXME CuVec?
    end

    # we stored the abs()^2, so take sqrt (this is faster to do all at once)
    if anyrd im_rd_0_1 = sqrt.(im_rd_0_1) end

    if (doS2 | doS20) & iso
        Mat2 = filter_hash["S2_iso_mat"]
    end

    # Real domain 2nd order
    if doS20
        Amat = reshape(im_rd_0_1, Nx*Ny, Nf, Nb)
        S20  = permutedims(Amat,(2,1,3)) ⊠ Amat
        if iso
            out_coeff[3+S1s:3+S1s+S2s-1,:] .= Mat2*reshape(S20,Nf*Nf,Nb)
        else
            out_coeff[3+S1s:3+S1s+S2s-1,:] .= reshape(S20,Nf*Nf,Nb)
        end
    end

    if doS2
        ## Traditional second order
        for f1 = 1:Nf
            thisim = P2*im_rd_0_1[:,:,f1,:]
            # Loop over f2 and do second-order convolution
            if normS1
                normS1pwr = S1[f1,:]
            elseif normS1iso
                normS1pwr = S1iso[f1,:]
            else
                normS1pwr = 1
            end
            for f2 = 1:Nf
                # sum im^2 = sum(|fft|^2/npix)
                S2[f1,f2,:] .= sum(abs2.(filts[:,:,f2].*thisim))/(Nx*Ny)/normS1pwr
            end
        end

        if iso
            out_coeff[end-S2s+1:end,:] .= Mat2*reshape(S2,Nf*Nf,Nb)
        else
            out_coeff[end-S2s+1:end,:] .= reshape(S2,Nf*Nf,Nb)
        end
    end

    return out_coeff
end

function eqws_compute_3d(image::Array{Float64,3}, filter_hash; FFTthreads=2, iso=false)
    # Use 2 threads for FFT
    FFTW.set_num_threads(FFTthreads)

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

    # allocate coeff arrays
    out_coeff = []
    S0  = zeros(Float64, 2)
    S1  = zeros(Float64, Nf)
    S20 = zeros(Float64, Nf, Nf)
    S12 = zeros(Float64, Nf, Nf)
    S2  = zeros(Float64, Nf, Nf)  # traditional 2nd order

    # allocate image arrays for internal use
    im_fdf_0_1 = zeros(Float64,           Nx, Ny, Nz, Nf)   # this must be zeroed!
    im_rd_0_1  = Array{Float64, 4}(undef, Nx, Ny, Nz, Nf)

    ## 0th Order
    S0[1]   = mean(image)
    norm_im = image.-S0[1]
    S0[2]   = sum(norm_im .* norm_im)/(Nx*Ny*Nz)
    norm_im ./= sqrt(Nx*Ny*Nz*S0[2])

    append!(out_coeff,S0[:])

    ## 1st Order
    im_fd_0 = fft(norm_im)  # 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_ind2   = filter_hash2["filt_index"]  # (J, L) array of filters represented as index value pairs
    #f_val2   = filter_hash2["filt_value"]

    zarr = zeros(ComplexF64, Nx, Ny, Nz)  # temporary array to fill with zvals

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

    ## Main 1st Order and Precompute 2nd Order
    for f = 1:Nf
        S1tot = 0.0
        f_i = f_ind[f]  # CartesianIndex list for filter
        f_v = f_val[f]  # Values for f_i
        # for (ind, val) in zip(f_i, f_v)   # this is slower!
        if length(f_i) > 0
            for i = 1:length(f_i)
                ind       = f_i[i]
                zval      = f_v[i] * im_fd_0[ind]
                S1tot    += abs2(zval)
                zarr[ind] = zval
                im_fdf_0_1[ind,f] = abs(zval)
            end
            S1[f] = S1tot/(Nx*Ny*Nz)  # image power ###why do we need to normalize again?
            im_rd_0_1[:,:,:,f] .= abs2.(P*zarr)
            zarr[f_i] .= 0
        end
    end
    append!(out_coeff, iso ? filter_hash["S1_iso_mat"]*S1 : S1)

    # we stored the abs()^2, so take sqrt (this is faster to do all at once)
    im_rd_0_1 .= sqrt.(im_rd_0_1)

    ## 2nd Order
    #Amat = reshape(im_fdf_0_1, Nx*Ny, Nf)
    #S12  = Amat' * Amat
    #Amat = reshape(im_rd_0_1, Nx*Ny, Nf)
    #S20  = Amat' * Amat

    #append!(out_coeff, S20)
    #append!(out_coeff, S12)

    ## Traditional second order
    if iso
        Mat2 = filter_hash["S2_iso_mat"]
    end
    for f1 = 1:Nf
        thisim = fft(im_rd_0_1[:,:,:,f1])  # Could try rfft here
        # println("  f1",f1,"  sum(fft):",sum(abs2.(thisim))/Nx^2, "  sum(im): ",sum(abs2.(im_rd_0_1[:,:,f1])))
        # Loop over f2 and do second-order convolution
        for f2 = 1:Nf
            f_i = f_ind[f2]  # CartesianIndex list for filter
            f_v = f_val[f2]  # Values for f_i
            # sum im^2 = sum(|fft|^2/npix)
            S2[f1,f2] = sum(abs2.(f_v .* thisim[f_i]))/(Nx*Ny*Nz)
        end
    end
    append!(out_coeff, iso ? Mat2*S2[:] : S2[:])

    return out_coeff
end

function eqws_compute_3d_gpu(image::Array{Float64,3}, filter_hash;doS2::Bool=true, doS12::Bool=false, norm=true)
    CUDA.reclaim()
    # array sizes
    (Nx, Ny, Nz)  = size(image)
    (Nf, )    = size(filter_hash["filt_vals"])

    # allocate coeff arrays
    out_coeff = []
    S0  = zeros(Float64, 2)
    S1  = zeros(Float64, Nf)
    #S20 = zeros(Float64, Nf, Nf)
    if doS12 S12 = zeros(Float64, Nf, Nf) end
    if doS2 S2  = zeros(Float64, Nf, Nf)  end # traditional 2nd order

    # allocate image arrays for internal use
    im_fdf_0_1 = CUDA.CuArray(Array{Float64, 3}(undef, Nx, Ny, Nz))
    im_rd_0_1  = CUDA.CuArray(Array{Float64, 3}(undef, Nx, Ny, Nz))

    image=CUDA.CuArray(image)

    ## 0th Order
    S0[1]   = mean(image)
    norm_im = image.-S0[1]
    S0[2]   = sum(norm_im .* norm_im)/(Nx*Ny*Nz)
    if norm
        norm_im ./= sqrt(Nx*Ny*Nz*S0[2])
    else
        norm_im=image
    end

    ## 1st Order
    im_fd_0 = CUDA.CUFFT.fftshift(CUDA.CUFFT.fft(norm_im))  # total power=1.0

    # unpack filter_hash
    f_mms   = filter_hash["filt_mms"]  # (J, L) array of filters represented as index value pairs
    f_vals   = filter_hash["filt_vals"]

    #f_ind2   = filter_hash2["filt_index"]  # (J, L) array of filters represented as index value pairs
    #f_val2   = filter_hash2["filt_value"]

    zarr = CUDA.zeros(ComplexF64, Nx, Ny, Nz)  # temporary array to fill with zvals

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

    #depth first search
    for f1 = 1:Nf
        im_fdf_0_1.=0.
        f_mm1 = f_mms[:,:,f1]  # CartesianIndex list for filter
        f_val1 = f_vals[f1]  # Values for f_i
        # for (ind, val) in zip(f_i, f_v)   # this is slower!

        zv=f_val1.*im_fd_0[f_mm1[1,1]:f_mm1[1,2],f_mm1[2,1]:f_mm1[2,2],f_mm1[3,1]:f_mm1[3,2]]
        S1[f1] = sum(abs2.(zv))/(Nx*Ny*Nz)

        zarr[f_mm1[1,1]:f_mm1[1,2],f_mm1[2,1]:f_mm1[2,2],f_mm1[3,1]:f_mm1[3,2]]=zv

        im_fdf_0_1[f_mm1[1,1]:f_mm1[1,2],f_mm1[2,1]:f_mm1[2,2],f_mm1[3,1]:f_mm1[3,2]] .= abs.(zv)

        im_rd_0_1=abs2.(P*CUDA.CUFFT.ifftshift(zarr))

        zarr[f_mm1[1,1]:f_mm1[1,2],f_mm1[2,1]:f_mm1[2,2],f_mm1[3,1]:f_mm1[3,2]] .= 0.

        im_rd_0_1 .= sqrt.(im_rd_0_1)
        im_fdf_1_1 = CUDA.CUFFT.fftshift(CUDA.CUFFT.fft(im_rd_0_1))
        for f2 = 1:Nf
            f_mm2 = f_mms[:,:,f2]
            f_val2 = f_vals[f2]

            S2[f1,f2] = sum(abs2.(f_val2 .* im_fdf_1_1[f_mm2[1,1]:f_mm2[1,2],f_mm2[2,1]:f_mm2[2,2],f_mm2[3,1]:f_mm2[3,2]] ))/(Nx*Ny*Nz)

            if !doS12
                continue
            end

            if f2>f1
                S12[f1,f2]=NaN
                continue
            end

            if f_mm1[1,1]>f_mm2[1,2] || f_mm2[1,1]>f_mm1[1,2] || f_mm1[2,1]>f_mm2[2,2] || f_mm2[2,1]>f_mm1[2,2] || f_mm1[3,1]>f_mm2[3,2] || f_mm2[3,1]>f_mm1[3,2]
                S12[f1,f2]=NaN
                continue
            end

            mx=max(f_mm1[1,1],f_mm2[1,1])
            my=max(f_mm1[2,1],f_mm2[2,1])
            mz=max(f_mm1[3,1],f_mm2[3,1])
            Mx=min(f_mm1[1,2],f_mm2[1,2])
            My=min(f_mm1[2,2],f_mm2[2,2])
            Mz=min(f_mm1[3,2],f_mm2[3,2])
            sfx=Mx-mx
            sfy=My-my
            sfz=Mz-mz
            st1x=mx-f_mm1[1,1]+1
            st1y=my-f_mm1[2,1]+1
            st1z=mz-f_mm1[3,1]+1
            st2x=mx-f_mm2[1,1]+1
            st2y=my-f_mm2[2,1]+1
            st2z=mz-f_mm2[3,1]+1

            S12[f1,f2]=sum(f_val1[st1x:st1x+sfx,st1y:st1y+sfy,st1z:st1z+sfz].*f_val2[st2x:st2x+sfx,st2y:st2y+sfy,st2z:st2z+sfz].*abs2.(im_fd_0[mx:Mx,my:My,mz:Mz]))
            CUDA.reclaim()
        end
        CUDA.reclaim()
    end
    append!(out_coeff,S0)
    append!(out_coeff, S1)
    if doS2 append!(out_coeff, S2) end
    if doS12 append!(out_coeff, S12) end
    return out_coeff
end

function eqws_compute_3d_gpu_fake(image::Array{Float64,3}, filter_hash;doS2::Bool=true, doS12::Bool=false, norm=true)
    # array sizes
    (Nx, Ny, Nz)  = size(image)
    (Nf, )    = size(filter_hash["filt_vals"])

    # allocate coeff arrays
    out_coeff = []
    S0  = zeros(Float64, 2)
    S1  = zeros(Float64, Nf)
    #S20 = zeros(Float64, Nf, Nf)
    if doS12 S12 = zeros(Float64, Nf, Nf) end
    if doS2 S2  = zeros(Float64, Nf, Nf)  end # traditional 2nd order

    # allocate image arrays for internal use
    im_fdf_0_1 = Array{Float64, 3}(undef, Nx, Ny, Nz)
    im_rd_0_1  = Array{Float64, 3}(undef, Nx, Ny, Nz)

    ## 0th Order
    S0[1]   = mean(image)
    norm_im = image.-S0[1]
    S0[2]   = sum(norm_im .* norm_im)/(Nx*Ny*Nz)
    if norm
        norm_im ./= sqrt(Nx*Ny*Nz*S0[2])
    else
        norm_im=image
    end

    ## 1st Order
    im_fd_0 = FFTW.fftshift(FFTW.fft(norm_im))  # total power=1.0

    # unpack filter_hash
    f_mms   = filter_hash["filt_mms"]  # (J, L) array of filters represented as index value pairs
    f_vals   = filter_hash["filt_vals"]

    #f_ind2   = filter_hash2["filt_index"]  # (J, L) array of filters represented as index value pairs
    #f_val2   = filter_hash2["filt_value"]

    zarr = zeros(ComplexF64, Nx, Ny, Nz)  # temporary array to fill with zvals

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

    #depth first search
    for f1 = 1:Nf
        im_fdf_0_1.=0.
        f_mm1 = f_mms[:,:,f1]  # CartesianIndex list for filter
        f_val1 = f_vals[f1]  # Values for f_i
        # for (ind, val) in zip(f_i, f_v)   # this is slower!

        zv=f_val1.*im_fd_0[f_mm1[1,1]:f_mm1[1,2],f_mm1[2,1]:f_mm1[2,2],f_mm1[3,1]:f_mm1[3,2]]
        S1[f1] = sum(abs2.(zv))/(Nx*Ny*Nz)

        zarr[f_mm1[1,1]:f_mm1[1,2],f_mm1[2,1]:f_mm1[2,2],f_mm1[3,1]:f_mm1[3,2]]=zv

        im_fdf_0_1[f_mm1[1,1]:f_mm1[1,2],f_mm1[2,1]:f_mm1[2,2],f_mm1[3,1]:f_mm1[3,2]] .= abs.(zv)

        im_rd_0_1=abs2.(P*FFTW.ifftshift(zarr))

        zarr[f_mm1[1,1]:f_mm1[1,2],f_mm1[2,1]:f_mm1[2,2],f_mm1[3,1]:f_mm1[3,2]] .= 0.

        im_rd_0_1 .= sqrt.(im_rd_0_1)
        im_fdf_1_1 = FFTW.fftshift(FFTW.fft(im_rd_0_1))
        for f2 = 1:Nf
            f_mm2 = f_mms[:,:,f2]
            f_val2 = f_vals[f2]

            S2[f1,f2] = sum(abs2.(f_val2 .* im_fdf_1_1[f_mm2[1,1]:f_mm2[1,2],f_mm2[2,1]:f_mm2[2,2],f_mm2[3,1]:f_mm2[3,2]] ))/(Nx*Ny*Nz)

            if !doS12
                continue
            end

            if f2>f1
                S12[f1,f2]=NaN
                continue
            end

            if f_mm1[1,1]>f_mm2[1,2] || f_mm2[1,1]>f_mm1[1,2] || f_mm1[2,1]>f_mm2[2,2] || f_mm2[2,1]>f_mm1[2,2] || f_mm1[3,1]>f_mm2[3,2] || f_mm2[3,1]>f_mm1[3,2]
                S12[f1,f2]=NaN
                continue
            end

            mx=max(f_mm1[1,1],f_mm2[1,1])
            my=max(f_mm1[2,1],f_mm2[2,1])
            mz=max(f_mm1[3,1],f_mm2[3,1])
            Mx=min(f_mm1[1,2],f_mm2[1,2])
            My=min(f_mm1[2,2],f_mm2[2,2])
            Mz=min(f_mm1[3,2],f_mm2[3,2])
            sfx=Mx-mx
            sfy=My-my
            sfz=Mz-mz
            st1x=mx-f_mm1[1,1]+1
            st1y=my-f_mm1[2,1]+1
            st1z=mz-f_mm1[3,1]+1
            st2x=mx-f_mm2[1,1]+1
            st2y=my-f_mm2[2,1]+1
            st2z=mz-f_mm2[3,1]+1

            S12[f1,f2]=sum(f_val1[st1x:st1x+sfx,st1y:st1y+sfy,st1z:st1z+sfz].*f_val2[st2x:st2x+sfx,st2y:st2y+sfy,st2z:st2z+sfz].*abs2.(im_fd_0[mx:Mx,my:My,mz:Mz]))
        end
    end
    append!(out_coeff,S0)
    append!(out_coeff, S1)
    if doS2 append!(out_coeff, S2) end
    if doS12 append!(out_coeff, S12) end
    return out_coeff
end

max_stars_repo_name: andrew-saydjari/DHC

max_issues_repo_name: andrew-saydjari/DHC

max_forks_repo_name: andrew-saydjari/DHC


text:
function entropy( img )
	sum = 0.0
	@inbounds @simd for idx in 1:length(img)
		sum += img[idx]
	end

	entrpy = 0.0
	@inbounds @simd for idx in 1:length(img)
		entrpy += (img[idx] == 0) ? 0.0 : log( 2, img[idx]/sum )*img[idx]/sum
	end

	return -entrpy
end


text:
module CLIWrappers

using TSML
using TSML.TSMLTypes
using TSML.TSMLTransformers
using TSML.Utils

using TSML.DataReaders
using TSML.DataWriters
using TSML.Statifiers
using TSML.Monotonicers

using Dates
using DataFrames
using CSV

export tsmlrun

const COMMONARG = Dict(:dateformat=>"dd/mm/yyyy HH:MM",
                        :dateinterval=>Dates.Hour(1),
                        :processmissing=>true
                       )

function tsmlrun(inputname::AbstractString,outputname::AbstractString="",datefmt::AbstractString="dd/mm/yyyy HH:MM",otype::AbstractString="table")
    if otype == "table"
        imputedoutput(inputname,outputname,datefmt)
    elseif otype == "stat"
        imputedstat(inputname,outputname,datefmt)
    else
        error("wrong output type")
    end
end

function rawstat(inputname::AbstractString,outputname::AbstractString="",datefmt::AbstractString="dd/mm/yyyy HH:MM")
    isfile(inputname) || error("input file name does not exist")
    csvreader = DataReader(Dict(:filename=>inputname,:dateformat=>datefmt))
    stfier = Statifier()
    lpipe = Pipeline(Dict(
        :transformers => [csvreader,stfier]
      )
    )
    fit!(lpipe)
    res=transform!(lpipe)
    if outputname != ""
        res |> CSV.write(outputname)
    end
    return res
end

function aggregatedstat(inputname::AbstractString,outputname::AbstractString="",datefmt::AbstractString="dd/mm/yyyy HH:MM")
    isfile(inputname) || error("input file name does not exist")
    csvreader = DataReader(Dict(:filename=>inputname,:dateformat=>datefmt))
    valgator = DateValgator(COMMONARG)
    stfier = Statifier()
    lpipe = Pipeline(Dict(
        :transformers => [csvreader,valgator,stfier]
      )
    )
    fit!(lpipe)
    res=transform!(lpipe)
    if outputname != ""
        res |> CSV.write(outputname)
    end
    return res
end

function aggregatedoutput(inputname::AbstractString,outputname::AbstractString="",datefmt::AbstractString="dd/mm/yyyy HH:MM")
    isfile(inputname) || error("input file name does not exist")
    csvreader = DataReader(Dict(:filename=>inputname,:dateformat=>datefmt))
    valgator = DateValgator(COMMONARG)
    lpipe = Pipeline(Dict(
        :transformers => [csvreader,valgator]
      )
    )
    fit!(lpipe)
    res=transform!(lpipe)
    if outputname != ""
        res |> CSV.write(outputname)
    end
    return res
end

function imputedstat(inputname::AbstractString,outputname::AbstractString="",datefmt::AbstractString="dd/mm/yyyy HH:MM")
    isfile(inputname) || error("input file name does not exist")
    csvreader = DataReader(Dict(:filename=>inputname,:dateformat=>datefmt))
    valgator = DateValgator(COMMONARG)
    valnner = DateValNNer(COMMONARG)
    mononicer = Monotonicer()
    stfier = Statifier()
    lpipe = Pipeline(Dict(
        :transformers => [csvreader,valgator,valnner,mononicer,stfier]
      )
    )
    fit!(lpipe)
    res=transform!(lpipe)
    if outputname != ""
        res |> CSV.write(outputname)
    end
    return res
end

function imputedoutput(inputname::AbstractString,outputname::AbstractString="",datefmt::AbstractString="dd/mm/yyyy HH:MM")
    isfile(inputname) || error("input file name does not exist")
    csvreader = DataReader(Dict(:filename=>inputname,:dateformat=>datefmt))
    valgator = DateValgator(COMMONARG)
    valnner = DateValNNer(COMMONARG)
    mononicer = Monotonicer()
    stfier = Statifier()
    lpipe = Pipeline(Dict(
        :transformers => [csvreader,valgator,valnner,mononicer]
      )
    )
    fit!(lpipe)
    res=transform!(lpipe)
    if outputname != ""
        res |> CSV.write(outputname)
    end
    return res
end

end


text:
# rounding.jl

module RoundingTests

using Test
using Dates
using QDates

# Basic rounding tests
@testset "Basic rounding" begin

qdt = QDates.QDate(2016, 2, 28)
@test floor(qdt, Dates.Year) == QDates.QDate(2016)
@test floor(qdt, Dates.Year(5)) == QDates.QDate(2015)
@test floor(qdt, Dates.Year(10)) == QDates.QDate(2010)
@test floor(qdt, Dates.Month) == QDates.QDate(2016, 2)
@test floor(qdt, Dates.Month(6)) == QDates.QDate(2015, 10)
@test floor(qdt, Dates.Day) == qdt
@test floor(qdt, Dates.Day(3)) == QDates.QDate(2016, 2, 26)
@test ceil(qdt, Dates.Year) == QDates.QDate(2017)
@test ceil(qdt, Dates.Year(5)) == QDates.QDate(2020)
@test ceil(qdt, Dates.Month) == QDates.QDate(2016, 3)
@test ceil(qdt, Dates.Month(6)) == QDates.QDate(2016, 4)
@test ceil(qdt, Dates.Day) == qdt
@test ceil(qdt, Dates.Day(3)) == QDates.QDate(2016, 2, 29)
@test round(qdt, Dates.Year) == QDates.QDate(2016)
@test round(qdt, Dates.Month) == QDates.QDate(2016, 3)

end

end