text:
abstract type AbstractLoadFormulation <: AbstractDeviceFormulation end

abstract type AbstractControllablePowerLoadFormulation <: AbstractLoadFormulation end

struct StaticPowerLoad <: AbstractLoadFormulation end

struct InterruptiblePowerLoad <: AbstractControllablePowerLoadFormulation end

struct DispatchablePowerLoad <: AbstractControllablePowerLoadFormulation end

########################### dispatchable load variables ############################################

function activepower_variables(canonical_model::CanonicalModel,
                               devices::IS.FlattenIteratorWrapper{L}) where {L<:PSY.ElectricLoad}
    add_variable(canonical_model,
                 devices,
                 Symbol("P_$(L)"),
                 false,
                 :nodal_balance_active, -1.0;
                 ub_value = x -> PSY.get_maxactivepower(x),
                 lb_value = x -> 0.0)

    return

end

function reactivepower_variables(canonical_model::CanonicalModel,
                                 devices::IS.FlattenIteratorWrapper{L}) where {L<:PSY.ElectricLoad}
    add_variable(canonical_model,
                 devices,
                 Symbol("Q_$(L)"),
                 false,
                 :nodal_balance_reactive, -1.0;
                 ub_value = x -> PSY.get_maxreactivepower(x),
                 lb_value = x -> 0.0)

    return

end

function commitment_variables(canonical_model::CanonicalModel,
                              devices::IS.FlattenIteratorWrapper{L}) where {L<:PSY.ElectricLoad}

    add_variable(canonical_model,
                 devices,
                 Symbol("ON_$(L)"),
                 true)

    return

end

####################################### Reactive Power Constraints ######################################
"""
Reactive Power Constraints on Loads Assume Constant PowerFactor
"""
function reactivepower_constraints(canonical_model::CanonicalModel,
                                   devices::IS.FlattenIteratorWrapper{L},
                                   device_formulation::Type{D},
                                   system_formulation::Type{S}) where {L<:PSY.ElectricLoad,
                                                                       D<:AbstractControllablePowerLoadFormulation,
                                                                       S<:PM.AbstractPowerModel}
    time_steps = model_time_steps(canonical_model)
    key = Symbol("reactive_$(L)")
    canonical_model.constraints[key] = JuMPConstraintArray(undef, (PSY.get_name(d) for d in devices), time_steps)

    for t in time_steps, d in devices
        name = PSY.get_name(d)
        pf = sin(atan((PSY.get_maxreactivepower(d)/PSY.get_maxactivepower(d))))
        canonical_model.constraints[key][PSY.get_name(d), t] = JuMP.@constraint(canonical_model.JuMPmodel,
                        canonical_model.variables[Symbol("Q_$(L)")][name, t] == canonical_model.variables[Symbol("P_$(L)")][name, t]*pf)
    end

    return

end

######################## output constraints without Time Series ###################################
function _get_time_series(devices::IS.FlattenIteratorWrapper{T},
                          time_steps::UnitRange{Int64}) where {T<:PSY.ElectricLoad}

    names = Vector{String}(undef, length(devices))
    series = Vector{Vector{Float64}}(undef, length(devices))

    for (ix, d) in enumerate(devices)
        names[ix] = PSY.get_name(d)
        series[ix] = fill(PSY.get_maxactivepower(d), (time_steps[end]))
    end

    return names, series

end

function activepower_constraints(canonical_model::CanonicalModel,
                                 devices::IS.FlattenIteratorWrapper{L},
                                 device_formulation::Type{DispatchablePowerLoad},
                                 system_formulation::Type{S}) where {L<:PSY.ElectricLoad,
                                                                     S<:PM.AbstractPowerModel}

    time_steps = model_time_steps(canonical_model)

    if model_has_parameters(canonical_model)
        device_timeseries_param_ub(canonical_model,
                                   _get_time_series(devices, time_steps),
                                   Symbol("active_$(L)"),
                                   UpdateRef{L}(Symbol("P_$(L)")),
                                   Symbol("P_$(L)"))
    else
        range_data = [(PSY.get_name(d), (min = 0.0, max = PSY.get_maxactivepower(d))) for d in devices]
        device_range(canonical_model,
                    range_data,
                    Symbol("activerange_$(L)"),
                    Symbol("P_$(L)")
                    )
    end

    return

end

"""
This function works only if the the Param_L <= PSY.get_maxactivepower(g)
"""
function activepower_constraints(canonical_model::CanonicalModel,
                                 devices::IS.FlattenIteratorWrapper{L},
                                 device_formulation::Type{InterruptiblePowerLoad},
                                 system_formulation::Type{S}) where {L<:PSY.ElectricLoad,
                                                          S<:PM.AbstractPowerModel}
    time_steps = model_time_steps(canonical_model)

    if model_has_parameters(canonical_model)
        device_timeseries_ub_bigM(canonical_model,
                                 _get_time_series(devices, time_steps),
                                 Symbol("active_$(L)"),
                                 Symbol("P_$(L)"),
                                 UpdateRef{L}(Symbol("P_$(L)")),
                                 Symbol("ON_$(L)"))
    else
        device_timeseries_ub_bin(canonical_model,
                                _get_time_series(devices, time_steps),
                                Symbol("active_$(L)"),
                                Symbol("P_$(L)"),
                                Symbol("ON_$(L)"))
    end

    return

end

######################### output constraints with Time Series ##############################################

function _get_time_series(forecasts::Vector{PSY.Deterministic{L}}) where {L<:PSY.ElectricLoad}

    names = Vector{String}(undef, length(forecasts))
    ratings = Vector{Float64}(undef, length(forecasts))
    series = Vector{Vector{Float64}}(undef, length(forecasts))

    for (ix, f) in enumerate(forecasts)
        component = PSY.get_component(f)
        names[ix] = PSY.get_name(component)
        series[ix] = values(PSY.get_timeseries(f))
        ratings[ix] = PSY.get_maxactivepower(component)
    end

    return names, ratings, series

end

function activepower_constraints(canonical_model::CanonicalModel,
                                 devices::Vector{PSY.Deterministic{L}},
                                 device_formulation::Type{DispatchablePowerLoad},
                                 system_formulation::Type{S}) where {L<:PSY.ElectricLoad,
                                                                     S<:PM.AbstractPowerModel}

    if model_has_parameters(canonical_model)
        device_timeseries_param_ub(canonical_model,
                                   _get_time_series(devices),
                                   Symbol("active_$(L)"),
                                   UpdateRef{L}(Symbol("P_$(L)")),
                                   Symbol("P_$(L)"))
    else
        device_timeseries_ub(canonical_model,
                            _get_time_series(devices),
                            Symbol("active_$(L)"),
                            Symbol("P_$(L)"))
    end

    return

end

function activepower_constraints(canonical_model::CanonicalModel,
                                 devices::Vector{PSY.Deterministic{L}},
                                 device_formulation::Type{InterruptiblePowerLoad},
                                 system_formulation::Type{S}) where {L<:PSY.ElectricLoad,
                                                                     S<:PM.AbstractPowerModel}

    if model_has_parameters(canonical_model)
        device_timeseries_ub_bigM(canonical_model,
                                 _get_time_series(devices),
                                 Symbol("active_$(L)"),
                                 Symbol("P_$(L)"),
                                 UpdateRef{L}(Symbol("P_$(L)")),
                                 Symbol("ON_$(L)"))
    else
        device_timeseries_ub_bin(canonical_model,
                                _get_time_series(devices),
                                Symbol("active_$(L)"),
                                Symbol("P_$(L)"),
                                Symbol("ON_$(L)"))
    end

    return

end

############################ injection expression with parameters ####################################

########################################### Devices ####################################################

function _nodal_expression_param(canonical_model::CanonicalModel,
                                devices::IS.FlattenIteratorWrapper{L},
                                system_formulation::Type{S}) where {L<:PSY.ElectricLoad,
                                                                    S<:PM.AbstractPowerModel}

    time_steps = model_time_steps(canonical_model)
    ts_data_active = Vector{Tuple{String, Int64, Float64, Vector{Float64}}}(undef, length(devices))
    ts_data_reactive = Vector{Tuple{String, Int64, Float64, Vector{Float64}}}(undef, length(devices))

    for (ix, d) in enumerate(devices)
        bus_number = PSY.get_number(PSY.get_bus(d))
        name = PSY.get_name(d)
        active_power = PSY.get_activepower(d)
        reactive_power = PSY.get_reactivepower(d)
        time_series_vector = ones(time_steps[end])
        ts_data_active[ix] = (name, bus_number, active_power, time_series_vector, -1.0)
        ts_data_reactive[ix] = (name, bus_number, reactive_power, time_series_vector, -1.0)
    end

    include_parameters(canonical_model,
                  ts_data_active,
                  UpdateRef{L}(Symbol("P_$(L)")),
                  :nodal_balance_active)
    include_parameters(canonical_model,
                   ts_data_reactive,
                   UpdateRef{L}(Symbol("Q_$(L)")),
                   :nodal_balance_reactive)

    return

end

function _nodal_expression_param(canonical_model::CanonicalModel,
                                devices::IS.FlattenIteratorWrapper{L},
                                system_formulation::Type{S}) where {L<:PSY.ElectricLoad,
                                                                    S<:PM.AbstractActivePowerModel}

    time_steps = model_time_steps(canonical_model)
    ts_data_active = Vector{Tuple{String, Int64, Float64, Vector{Float64}}}(undef, length(devices))

    for (ix, d) in enumerate(devices)
        bus_number = PSY.get_number(PSY.get_bus(d))
        name = PSY.get_name(d)
        active_power = PSY.get_activepower(d)
        time_series_vector = ones(time_steps[end])
        ts_data_active[ix] = (name, bus_number, active_power, time_series_vector, -1.0)
    end

    include_parameters(canonical_model,
                  ts_data_active,
                  UpdateRef{L}(Symbol("P_$(L)")),
                  :nodal_balance_active)

    return

end

############################################## Time Series ###################################
function _nodal_expression_param(canonical_model::CanonicalModel,
                                forecasts::Vector{PSY.Deterministic{L}},
                                system_formulation::Type{S}) where {L<:PSY.ElectricLoad,
                                                                    S<:PM.AbstractPowerModel}

    time_steps = model_time_steps(canonical_model)

    ts_data_active = Vector{Tuple{String, Int64, Float64, Vector{Float64}}}(undef, length(forecasts))
    ts_data_reactive = Vector{Tuple{String, Int64, Float64, Vector{Float64}}}(undef, length(forecasts))

    for (ix, f) in enumerate(forecasts)
        device = PSY.get_component(f)
        bus_number = PSY.get_number(PSY.get_bus(device))
        name = PSY.get_name(device)
        active_power = PSY.get_maxactivepower(device)
        reactive_power = PSY.get_maxreactivepower(device)
        time_series_vector = values(PSY.get_timeseries(f))
        ts_data_active[ix] = (name, bus_number, active_power, time_series_vector)
        ts_data_reactive[ix] = (name, bus_number, reactive_power, time_series_vector)
    end

    include_parameters(canonical_model,
                    ts_data_active,
                    UpdateRef{L}(Symbol("P_$(L)")),
                    :nodal_balance_active,
                    -1.0)
    include_parameters(canonical_model,
                    ts_data_reactive,
                    UpdateRef{L}(Symbol("Q_$(L)")),
                    :nodal_balance_reactive,
                    -1.0)

    return

end

function _nodal_expression_param(canonical_model::CanonicalModel,
                                forecasts::Vector{PSY.Deterministic{L}},
                                system_formulation::Type{S}) where {L<:PSY.ElectricLoad,
                                                                    S<:PM.AbstractActivePowerModel}

    time_steps = model_time_steps(canonical_model)
    ts_data_active = Vector{Tuple{String, Int64, Float64, Vector{Float64}}}(undef, length(forecasts))

    for (ix, f) in enumerate(forecasts)
        device = PSY.get_component(f)
        bus_number = PSY.get_number(PSY.get_bus(device))
        name = PSY.get_name(device)
        active_power = PSY.get_maxactivepower(device)
        time_series_vector = values(PSY.get_timeseries(f))
        ts_data_active[ix] = (name, bus_number, active_power, time_series_vector)
    end

    include_parameters(canonical_model,
                    ts_data_active,
                    UpdateRef{L}(Symbol("P_$(L)")),
                    :nodal_balance_active,
                    -1.0)

    return

end

############################ injection expression with fixed values ####################################

########################################### Devices ####################################################

function _nodal_expression_fixed(canonical_model::CanonicalModel,
                                devices::IS.FlattenIteratorWrapper{L},
                                system_formulation::Type{S}) where {L<:PSY.ElectricLoad,
                                                                    S<:PM.AbstractPowerModel}

    time_steps = model_time_steps(canonical_model)

    for t in time_steps, d in devices
        bus_number = PSY.get_number(PSY.get_bus(d))
        active_power = PSY.get_activepower(d)
        reactive_power = PSY.get_reactivepower(d)
        _add_to_expression!(canonical_model.expressions[:nodal_balance_active],
                            bus_number,
                            t,
                            -1*active_power);
        _add_to_expression!(canonical_model.expressions[:nodal_balance_reactive],
                            bus_number,
                            t,
                            -1*reactive_power);
    end

    return

end

function _nodal_expression_fixed(canonical_model::CanonicalModel,
                                devices::IS.FlattenIteratorWrapper{L},
                                system_formulation::Type{S}) where {L<:PSY.ElectricLoad,
                                                                    S<:PM.AbstractActivePowerModel}

    time_steps = model_time_steps(canonical_model)

    for t in time_steps, d in devices
        bus_number = PSY.get_number(PSY.get_bus(d))
        active_power = PSY.get_activepower(d)
        _add_to_expression!(canonical_model.expressions[:nodal_balance_active],
                            bus_number,
                            t,
                            -1*active_power);
    end

    return

end

############################################## Time Series ###################################

function _nodal_expression_fixed(canonical_model::CanonicalModel,
                                forecasts::Vector{PSY.Deterministic{L}},
                                system_formulation::Type{S}) where {L<:PSY.ElectricLoad,
                                                                    S<:PM.AbstractPowerModel}

    time_steps = model_time_steps(canonical_model)

    for f in forecasts
        device = PSY.get_component(f)
        bus_number = PSY.get_number(PSY.get_bus(device))
        active_power = PSY.get_maxactivepower(device)
        reactive_power = PSY.get_maxreactivepower(device)
        time_series_vector = values(PSY.get_timeseries(f))
        for t in time_steps
            bus_number = PSY.get_number(PSY.get_bus(device))
            _add_to_expression!(canonical_model.expressions[:nodal_balance_active],
                                bus_number,
                                t,
                                -1 * time_series_vector[t] * active_power)
            _add_to_expression!(canonical_model.expressions[:nodal_balance_reactive],
                                bus_number,
                                t,
                                -1 * time_series_vector[t] * reactive_power)
        end
    end

    return

end

function _nodal_expression_fixed(canonical_model::CanonicalModel,
                                forecasts::Vector{PSY.Deterministic{L}},
                                system_formulation::Type{S}) where {L<:PSY.ElectricLoad,
                                                                    S<:PM.AbstractActivePowerModel}

    time_steps = model_time_steps(canonical_model)

    for f in forecasts
        device = PSY.get_component(f)
        bus_number = PSY.get_number(PSY.get_bus(device))
        active_power = PSY.get_maxactivepower(device)
        time_series_vector = values(PSY.get_timeseries(f))
        for t in time_steps
            bus_number = PSY.get_number(PSY.get_bus(device))
            _add_to_expression!(canonical_model.expressions[:nodal_balance_active],
                                bus_number,
                                t,
                                -1 * time_series_vector[t] * active_power)
        end
    end

    return

end

##################################### FormulationControllable Load Cost ######################################

function cost_function(canonical_model::CanonicalModel,
                       devices::IS.FlattenIteratorWrapper{L},
                       device_formulation::Type{DispatchablePowerLoad},
                       system_formulation::Type{S}) where {L<:PSY.ControllableLoad,
                                                           S<:PM.AbstractPowerModel}

    add_to_cost(canonical_model,
                devices,
                Symbol("P_$(L)"),
                :variable,
                -1.0)

    return

end

function cost_function(canonical_model::CanonicalModel,
                       devices::IS.FlattenIteratorWrapper{L},
                       device_formulation::Type{InterruptiblePowerLoad},
                       system_formulation::Type{S}) where {L<:PSY.ControllableLoad,
                                                           S<:PM.AbstractPowerModel}

    add_to_cost(canonical_model,
                devices,
                Symbol("ON_$(L)"),
                :fixed,
                -1.0)

    return

end


text:
# This trajectory is an example for the chemtrajectorys library
# Any copyright is dedicated to the Public Domain.
# http://creativecommons.org/publicdomain/zero/1.0/
#!/usr/bin/env julia
using chemfiles

trajectory = Trajectory("filename.xyz")
frame = read(trajectory)

less_than_five = UInt64[]
pos = positions(frame)

for i in 0:size(frame)
    if pos[1, i] < 5
        push!(less_than_five, i)
    end
end

println("Atoms with x < 5: ")
for i in less_than_five
    println("  - $i")
end

ext: jl


text:
const WASD_KEYS  = (GLFW.KEY_W, GLFW.KEY_A, GLFW.KEY_S, GLFW.KEY_D)
const SHIFT_KEYS = (GLFW.KEY_LEFT_SHIFT, GLFW.KEY_RIGHT_SHIFT)
const CTRL_KEYS  = (GLFW.KEY_LEFT_CONTROL, GLFW.KEY_RIGHT_CONTROL)

@component mutable struct Mouse
    x::Float32
    y::Float32
    dx::Float32
    dy::Float32
    scroll::NTuple{2,Float32}
    dscroll::NTuple{2,Float32}
    button::GLFW.MouseButton
    action::GLFW.Action
end

@component mutable struct Keyboard
    button::GLFW.Key
    modifiers::Int
    action::GLFW.Action
end

struct EventPoller <: System end

Overseer.requested_components(::EventPoller) = (Mouse, Keyboard, Canvas)

function Overseer.update(::EventPoller, m::AbstractLedger)
    c = singleton(m, Canvas)
    mouse = singleton(m, Mouse)
    keyboard = singleton(m, Keyboard)
    pollevents(c)

    x, y = c.cursor_position

    mouse_button    = c.mouse_buttons
    keyboard_button = c.keyboard_buttons

    mouse.dx      = x - mouse.x
    mouse.dy      = y - mouse.y
    mouse.x       = x
    mouse.y       = y
    mouse.button  = mouse_button[1]
    mouse.action  = mouse_button[2]
    mouse.dscroll = c.scroll .- mouse.scroll
    mouse.scroll  = c.scroll

    keyboard.button = keyboard_button[1]
    keyboard.action = keyboard_button[3]
    return keyboard.modifiers = keyboard_button[4]
end

pressed(m::Mouse)     = m.action ∈ (GLFW.PRESS, GLFW.REPEAT)
pressed(k::Keyboard)  = k.action ∈ (GLFW.PRESS, GLFW.REPEAT)
released(m::Mouse)    = m.action == GLFW.RELEASE
released(k::Keyboard) = k.action == GLFW.RELEASE

# Mouse Picking Stuff
@component @with_kw mutable struct Selectable
    color::RGBf0 = RGBf0(0, 0, 0)
    selected::Bool = false
    color_modifier::Float32 = 1.3f0
end

struct IDColorGenerator <: System end

Overseer.requested_components(::IDColorGenerator) = (PolygonGeometry, Selectable)

function Overseer.update(::IDColorGenerator, m::AbstractLedger)
    sel = m[Selectable]
    i = length(sel)
    @inbounds for s in sel
        s.color = RGBf0(((i & 0x000000FF) >> 0) / 255, ((i & 0x0000FF00) >> 8) / 255,
                        ((i & 0x00FF0000) >> 16) / 255)
        i += 1
    end
end

struct MousePicker <: System end

function Overseer.requested_components(::MousePicker)
    return (Selectable, Camera3D, Spatial, UniformColor, Canvas, UpdatedComponents)
end

function Overseer.update(::MousePicker, m::AbstractLedger)
    ucolor             = m[UniformColor]
    bcolor             = m[BufferColor]
    canvas             = singleton(m, Canvas)
    updated_components = singleton(m, UpdatedComponents)
    mouse              = singleton(m, Mouse)
    keyboard           = singleton(m, Keyboard)
    wh                 = size(canvas)

    iofbo = singleton(m, IOTarget).target

    @inbounds if pressed(keyboard) &&
                 keyboard.button ∈ CTRL_KEYS &&
                 mouse.button == GLFW.MOUSE_BUTTON_1 &&
                 pressed(mouse)
        glFlush()
        bind(iofbo)
        glReadBuffer(GL_COLOR_ATTACHMENT1)

        glPixelStorei(GL_UNPACK_ALIGNMENT, 1)
        dat = Ref{RGBf0}()
        # screenspace = (2 * cursor_position[1] / wh[1] - 1,  1-2*cursor_position[2]/wh[2])
        glReadPixels(mouse.x, wh[2] - mouse.y, 1, 1, GL_RGB, GL_FLOAT, dat)
        unbind(iofbo)
        selected_color = dat[]
        for e in @entities_in(m, Selectable)
            mod = e.color_modifier
            if isapprox(e.color, selected_color)
                was_selected = e.selected
                e.selected = true
                if !was_selected
                    if e in ucolor
                        o_c = ucolor[e].color
                        ucolor[e] = UniformColor(RGB(o_c.r * mod, o_c.g * mod, o_c.b * mod))
                        register_update(m, UniformColor)
                    elseif e in bcolor
                        bcolor[e].color .*= mod
                        register_update(m, BufferColor)
                    end
                end
            elseif !isapprox(e.color, selected_color)
                was_not_selected = e.selected
                e.selected = false
                if was_not_selected
                    if e in ucolor
                        o_c = ucolor[e].color
                        ucolor[e] = UniformColor(RGB(o_c.r / mod, o_c.g / mod, o_c.b / mod))
                        register_update(m, UniformColor)
                    elseif e in bcolor
                        bcolor[e].color ./= mod
                        register_update(m, BufferColor)
                    end
                end
            end
        end
    end
end

function mouse_to_screenspace(m::Mouse, c::Canvas)
    return (2 * m.x / size(c, 1) - 1, 1 - 2 * m.y / size(c, 2))
end

struct Ray
    origin::Point3f0
    direction::Vec3f0
end

function Ray(m::Mouse, c::Canvas, camera::Camera3D, eye::Point3f0)
    screenspace = mouse_to_screenspace(m, c)
    ray_clip    = Vec4f0(screenspace..., -1.0, 1.0)
    cam_proj    = camera.proj
    cam_view    = camera.view
    ray_eye     = Vec4f0((inv(cam_proj)*ray_clip)[1:2]..., -1.0, 0.0)
    ray_dir     = normalize(Vec3f0((inv(cam_view)*ray_eye)[1:3]...))
    return Ray(eye, ray_dir)
end

struct Plane
    origin::Point3f0
    normal::Vec3f0
end

function intersect(r::Ray, p::Plane)
    return r.origin +
           r.direction * dot((p.origin - r.origin), p.normal) / (dot(r.direction, p.normal))
end


text:
"""
    PathfinderResult

Container for results of single-path Pathfinder.

# Fields
- `input`: User-provided input object, e.g. either `logp`, `(logp, ∇logp)`, `optim_fun`,
    `optim_prob`, or another object.
- `optimizer`: Optimizer used for maximizing the log-density
- `rng`: Pseudorandom number generator that was used for sampling
- `optim_prob::SciMLBase.OptimizationProblem`: Otimization problem used for
    optimization
- `logp`: Log-density function
- `fit_distribution::Distributions.MvNormal`: ELBO-maximizing multivariate normal
    distribution
- `draws::AbstractMatrix{<:Real}`: draws from multivariate normal with size `(dim, ndraws)`
- `fit_distribution_transformed`: `fit_distribution` transformed to the same space as the
    user-supplied target distribution. This is only different from `fit_distribution` when
    integrating with other packages, and its type depends on the type of `input`.
- `draws_transformed`: `draws` transformed to be draws from `fit_distribution_transformed`.
- `fit_iteration::Int`: Iteration at which ELBO estimate was maximized
- `num_tries::Int`: Number of tries until Pathfinder succeeded
- `optim_solution::SciMLBase.OptimizationSolution`: Solution object of optimization.
- `optim_trace::Pathfinder.OptimizationTrace`: container for optimization trace of points,
    log-density, and gradient. The first point is the initial point.
- `fit_distributions::AbstractVector{Distributions.MvNormal}`: Multivariate normal
    distributions for each point in `optim_trace`, where
    `fit_distributions[fit_iteration + 1] == fit_distribution`
- `elbo_estimates::AbstractVector{<:Pathfinder.ELBOEstimate}`: ELBO estimates for all but
    the first distribution in `fit_distributions`.

# Returns
- [`PathfinderResult`](@ref)
"""
struct PathfinderResult{I,O,R,OP,LP,FD,D,FDT,DT,OS,OT,EE}
    input::I
    optimizer::O
    rng::R
    optim_prob::OP
    logp::LP
    fit_distribution::FD
    draws::D
    fit_distribution_transformed::FDT
    draws_transformed::DT
    fit_iteration::Int
    num_tries::Int
    optim_solution::OS
    optim_trace::OT
    fit_distributions::Vector{FD}
    elbo_estimates::EE
end

function Base.show(io::IO, ::MIME"text/plain", result::PathfinderResult)
    println(io, "Single-path Pathfinder result")
    println(io, "  tries: $(result.num_tries)")
    println(io, "  draws: $(size(result.draws, 2))")
    println(
        io, "  fit iteration: $(result.fit_iteration) (total: $(length(result.optim_trace) - 1))"
    )
    println(io, "  fit ELBO: $(_to_string(result.elbo_estimates[result.fit_iteration]))")
    print(io, "  fit distribution: ", result.fit_distribution)
    return nothing
end

"""
    pathfinder(logp; kwargs...)
    pathfinder(logp, ∇logp; kwargs...)
    pathfinder(fun::SciMLBase::OptimizationFunction; kwargs...)
    pathfinder(prob::SciMLBase::OptimizationProblem; kwargs...)

Find the best multivariate normal approximation encountered while maximizing `logp`.

From an optimization trajectory, Pathfinder constructs a sequence of (multivariate normal)
approximations to the distribution specified by `logp`. The approximation that maximizes the
evidence lower bound (ELBO), or equivalently, minimizes the KL divergence between the
approximation and the true distribution, is returned.

The covariance of the multivariate normal distribution is an inverse Hessian approximation
constructed using at most the previous `history_length` steps.

# Arguments
- `logp`: a callable that computes the log-density of the target distribution.
- `∇logp`: a callable that computes the gradient of `logp`. If not provided, `logp` is
    automatically differentiated using the backend specified in `ad_backend`.
- `fun::SciMLBase.OptimizationFunction`: an optimization function that represents
    `-logp(x)` with its gradient. It must have the necessary features (e.g. a Hessian
    function) for the chosen optimization algorithm. For details, see
    [GalacticOptim.jl: OptimizationFunction](https://galacticoptim.sciml.ai/stable/API/optimization_function/).
- `prob::SciMLBase.OptimizationProblem`: an optimization problem containing a function with
    the same properties as `fun`, as well as an initial point, in which case `init` and
    `dim` are ignored.

# Keywords
- `dim::Int`: dimension of the target distribution. If not provided, `init` or must be.
    Ignored if `init` is provided.
- `init::AbstractVector{<:Real}`: initial point of length `dim` from which to begin
    optimization. If not provided, an initial point of type `Vector{Float64}` and length
    `dim` is created and filled using `init_sampler`.
- `init_scale::Real`: scale factor ``s`` such that the default `init_sampler` samples
    entries uniformly in the range ``[-s, s]``
- `init_sampler`: function with the signature `(rng, x) -> x` that modifies a vector of
    length `dims` in-place to generate an initial point
- `ndraws_elbo::Int=$DEFAULT_NDRAWS_ELBO`: Number of draws used to estimate the ELBO
- `ndraws::Int=ndraws_elbo`: number of approximate draws to return
- `ad_backend=AD.ForwardDiffBackend()`: AbstractDifferentiation.jl AD backend.
- `rng::Random.AbstractRNG`: The random number generator to be used for drawing samples
- `executor::Transducers.Executor=Transducers.SequentialEx()`: Transducers.jl executor that
    determines if and how to perform ELBO computation in parallel. The default
    (`SequentialEx()`) performs no parallelization. If `rng` is known to be thread-safe, and
    the log-density function is known to have no internal state, then
    `Transducers.PreferParallel()` may be used to parallelize log-density evaluation.
    This is generally only faster for expensive log density functions.
- `history_length::Int=$DEFAULT_HISTORY_LENGTH`: Size of the history used to approximate the
    inverse Hessian.
- `optimizer`: Optimizer to be used for constructing trajectory. Can be any optimizer
    compatible with GalacticOptim, so long as it supports callbacks. Defaults to
    `Optim.LBFGS(; m=history_length, linesearch=LineSearches.MoreThuente())`. See
    the [GalacticOptim.jl documentation](https://galacticoptim.sciml.ai/stable) for details.
- `ntries::Int=1_000`: Number of times to try the optimization, restarting if it fails. Before
    every restart, a new initial point is drawn using `init_sampler`.
- `kwargs...` : Remaining keywords are forwarded to
    [`GalacticOptim.solve`](https://galacticoptim.sciml.ai/stable/API/solve).

# Returns
- [`PathfinderResult`](@ref)
"""
function pathfinder end

function pathfinder(logp; ad_backend=AD.ForwardDiffBackend(), kwargs...)
    return pathfinder(build_optim_function(logp; ad_backend); input=logp, kwargs...)
end
function pathfinder(logp, ∇logp; ad_backend=AD.ForwardDiffBackend(), kwargs...)
    return pathfinder(
        build_optim_function(logp, ∇logp; ad_backend); input=(logp, ∇logp), kwargs...
    )
end
function pathfinder(
    optim_fun::SciMLBase.OptimizationFunction;
    rng=Random.GLOBAL_RNG,
    init=nothing,
    dim::Int=-1,
    init_scale=2,
    init_sampler=UniformSampler(init_scale),
    input=optim_fun,
    kwargs...,
)
    if init !== nothing
        _init = init
        allow_mutating_init = false
    elseif init === nothing && dim > 0
        _init = Vector{Float64}(undef, dim)
        init_sampler(rng, _init)
        allow_mutating_init = true
    else
        throw(ArgumentError("An initial point `init` or dimension `dim` must be provided."))
    end
    prob = build_optim_problem(optim_fun, _init)
    return pathfinder(prob; rng, input, init_sampler, allow_mutating_init, kwargs...)
end
function pathfinder(
    prob::SciMLBase.OptimizationProblem;
    rng::Random.AbstractRNG=Random.GLOBAL_RNG,
    history_length::Int=DEFAULT_HISTORY_LENGTH,
    optimizer=default_optimizer(history_length),
    ndraws_elbo::Int=DEFAULT_NDRAWS_ELBO,
    ndraws::Int=ndraws_elbo,
    input=prob,
    kwargs...,
)
    if prob.f.grad === nothing || prob.f.grad isa Bool
        throw(ArgumentError("optimization function must define a gradient function."))
    end
    logp(x) = -prob.f.f(x, nothing)
    path_result = ProgressLogging.progress(; name="Optimizing") do progress_id
        return _pathfinder_try_until_succeed(
            rng,
            prob,
            logp;
            history_length,
            optimizer,
            progress_id,
            ndraws_elbo,
            kwargs...,
        )
    end
    @unpack (
        itry,
        success,
        optim_solution,
        optim_trace,
        fit_distributions,
        fit_iteration,
        elbo_estimates,
    ) = path_result

    if !success
        ndraws_elbo_actual = 0
        @warn "Pathfinder failed after $itry tries. Increase `ntries`, inspect the model for numerical instability, or provide a more suitable `init_sampler`."
    else
        elbo_estimate_opt = elbo_estimates[fit_iteration]
        ndraws_elbo_actual = ndraws_elbo
    end

    fit_distribution = fit_distributions[fit_iteration + 1]

    # reuse existing draws; draw additional ones if necessary
    draws = if ndraws_elbo_actual == 0
        rand(rng, fit_distribution, ndraws)
    elseif ndraws_elbo_actual < ndraws
        hcat(elbo_estimate_opt.draws, rand(rng, fit_distribution, ndraws - ndraws_elbo_actual))
    else
        elbo_estimate_opt.draws[:, 1:ndraws]
    end

    # placeholders
    fit_distribution_transformed = fit_distribution
    draws_transformed = draws

    return PathfinderResult(
        input,
        optimizer,
        rng,
        optim_solution.prob,
        logp,
        fit_distribution,
        draws,
        fit_distribution_transformed,
        draws_transformed,
        fit_iteration,
        itry,
        optim_solution,
        optim_trace,
        fit_distributions,
        elbo_estimates,
    )
end

function _pathfinder_try_until_succeed(
    rng,
    prob,
    logp;
    ntries::Int=1_000,
    init_scale=2,
    init_sampler=UniformSampler(init_scale),
    allow_mutating_init::Bool=false,
    kwargs...,
)
    itry = 1
    progress_name = "Optimizing (try 1)"
    result = _pathfinder(rng, prob, logp; progress_name, kwargs...)
    _prob = prob
    while !result.success && itry < ntries
        if itry == 1 && !allow_mutating_init
            _prob = deepcopy(prob)
        end
        itry += 1
        init_sampler(rng, _prob.u0)
        progress_name = "Optimizing (try $itry)"
        result = _pathfinder(rng, _prob, logp; progress_name, kwargs...)
    end
    return (; itry, result...)
end

function _pathfinder(
    rng,
    prob,
    logp;
    history_length::Int=DEFAULT_HISTORY_LENGTH,
    optimizer=default_optimizer(history_length),
    ndraws_elbo=DEFAULT_NDRAWS_ELBO,
    executor::Transducers.Executor=Transducers.SequentialEx(),
    kwargs...,
)
    # compute trajectory
    optim_solution, optim_trace = optimize_with_trace(prob, optimizer; kwargs...)
    L = length(optim_trace) - 1
    success = L > 0

    # fit mv-normal distributions to trajectory
    fit_distributions = fit_mvnormals(
        optim_trace.points, optim_trace.gradients; history_length
    )

    # find ELBO-maximizing distribution
    fit_iteration, elbo_estimates = @views maximize_elbo(
        rng, logp, fit_distributions[(begin + 1):end], ndraws_elbo, executor
    )
    if isempty(elbo_estimates)
        success = false
    else
        elbo = elbo_estimates[fit_iteration].value
        success &= !isnan(elbo) & (elbo != -Inf)
    end

    return (;
        success,
        optim_solution,
        optim_trace,
        fit_distributions,
        fit_iteration,
        elbo_estimates,
    )
end

"""
    UniformSampler(scale::Real)

Sampler that in-place modifies an array to be IID uniformly distributed on `[-scale, scale]`
"""
struct UniformSampler{T<:Real}
    scale::T
    function UniformSampler(scale::T) where {T<:Real}
        scale > 0 || throw(DomainError(scale, "scale of uniform sampler must be positive."))
        return new{T}(scale)
    end
end

function (s::UniformSampler)(rng::Random.AbstractRNG, point)
    scale = s.scale
    @. point = rand(rng) * 2scale - scale
    return point
end


text:
#
println("Ag) Test of the QED model corrections to the level structure of atoms and ions.")
#
wa = Atomic.Computation("QED estimates for carbon-like Xe",  Nuclear.Model(54.);
                        configs=[Configuration("1s^2 2s^2 2p^6")],
                        asfSettings=AsfSettings(true, false, "meanDFS", "hydrogenic", "xxx", [1],    40, 1.0e-6, JAC.Subshell[],
                                                true, false, NoneQed(), "yyy", LSjjSettings(true),
                                                false, [1,2,3,4], false, JAC.LevelSymmetry[] )  )

wb = perform(wa)


text:
#cd(@__DIR__)
using Git

cd(dirname(@__FILE__))
Git.run(`submodule init`)
Git.run(`submodule update`)


text:
export EChart

mutable struct EChart

  options::Dict
  width::Int64
  height::Int64
  html::String

EChart(options::Dict)=new(options,800,500,"")
EChart(options::Dict,width::Int64,height::Int64)=new(options,width,height,"")

end

function axys(arr)

    if eltype(arr)<:Number

        return Dict("type"=>"value")

    elseif eltype(arr)<:Date

        return Dict("type"=>"time")

    elseif length(arr)!=0

        return Dict("type"=>"category","data"=>sort(unique(arr)))

    else
      return ""

    end

end

function EChart(;x=[], y=[],symbolSize=[],series=[],names=[],ctype="bar", width=800,height=600,options=Dict())

    leng=(x==[] ? length(y) : length(x))

    if x!=[] && y!=[]
    xAxis=axys(x)
    yAxis=axys(y)
    xAxis=="" ? "" : options["xAxis"]=xAxis
    options["xAxis"]=convert(Dict{Any,Any},options["xAxis"])
    yAxis=="" ? "" : options["yAxis"]=yAxis
    options["yAxis"]=convert(Dict{Any,Any},options["yAxis"])
    end

    if series==[]
    ### Single series
    ser=Dict("Series1"=>[i for i in 1:leng])
    else

        ser=class_series(series)
        options["legend"]=Dict("data"=>collect(keys(ser)))
    end

        options["series"]=[]
        for s in ser

            push!(options["series"],Dict("type"=>ctype,"name"=>s.first,"data"=>treat_series(x=get_pos(x,s.second),
                        y=get_pos(y,s.second))))

            act_ser=length(options["series"])

            ### Additional series attributes
            ### symbolSize
            if symbolSize!=[]
              s_arr=get_pos(symbolSize,s.second)
              options["series"][act_ser]["data"]=convert(Array{Dict{String,Any},1},options["series"][act_ser]["data"])
              for i in 1:length(options["series"][act_ser]["data"])
              options["series"][act_ser]["data"][i]["symbolSize"]=s_arr[i]
              end
            end
            ### names
            if names!=[]
              s_arr=get_pos(names,s.second)
              options["series"][act_ser]["data"]=convert(Array{Dict{String,Any},1},options["series"][act_ser]["data"])
              for i in 1:length(options["series"][act_ser]["data"])
              options["series"][act_ser]["data"][i]["name"]=s_arr[i]
              end
            end

        end

    ### Series Type

    return EChart(options,width,height)

end

function treat_series(;x=[],y=[])

    lx=length(x)
    ly=length(y)

    if lx==0
         return [Dict("value"=>[y[i]]) for i in 1:ly]
    elseif ly==0
        return [Dict("value"=>[x[i]]) for i in 1:lx]
    else
        return [Dict("value"=>[x[i],y[i]]) for i in 1:lx]
    end

end

function class_series(arr)
ret=Dict()
for i in 1:length(arr)

    haskey(ret,arr[i]) ? push!(ret[arr[i]],i) : ret[arr[i]]=[i]

end
return ret
end

function get_pos(arr,pos)
arr_out=[]
  if length(arr)==0
          return arr_out
      else
          for i in pos

         push!(arr_out,arr[i])
          end
        return arr_out
  end

end


text:
println("Loading helpers ...")

#include("errors.jl")
include("eigenvalues.jl")
include("search.jl")
#include("woodbury.jl")


text:
using Statistics
using Plots
using BenchmarkTools
using Profile
using FFTW
using Statistics
using Optim
using Images, FileIO, ImageIO
using Printf
using SparseArrays
using Distributions
using Measures
using ImageFiltering
#using Plotly
#using Weave

using Profile
using LinearAlgebra

push!(LOAD_PATH, pwd()*"/scratch_NM")
using Deriv_Utils
push!(LOAD_PATH, pwd()*"/main")
using DHC_2DUtils#: fink_filter_bank, fink_filter_hash, DHC_compute

#=function img_reconfunc(input, s_targ_mean, s_targ_invcov, mask, coeff_choice)
    #Only for denoising rn. Need to adapt to use mask and coeff_choice

    (Nx, Ny)  = size(input)
    if Nx != Ny error("Input image must be square") end
    filter_hash = fink_filter_hash(1, 8, nx=Nx, pc=1, wd=1)

    (Nf, )    = size(filter_hash["filt_index"])
    if Nf == 0 error("filter hash corrupted") end

    storage_beta = Array{Float64, 2}(undef, 1, Nx^2)

    function loss_func(img_curr)
        s_curr = Deriv_Utils.DHC(img_curr,  filter_hash, doS2=true, doS12=false, doS20=false)
        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)
        s_curr = Deriv_Utils.DHC(img_curr, filter_hash, doS2=true, doS12=false, doS20=false)
        diff = s_curr - s_targ_mean
        diff = diff[3:end]
        dS1S2 = Deriv_Utils.wst_S1S2_derivfast(img_curr, filter_hash)
        term1 = Transpose(s_targ_invcov * diff)
        term2 = reshape(dS1S2, (size(dS1S2)[1], Nx^2))
        mul!(storage_beta, term1, term2)
        storage_grad = reshape(storage_beta, 1, Nx, Nx) #term1: 1x (Nf+Nf^2) | term2: #(Nf+Nf^2) x Nx x Nx
        #storage_grad .=  term1 * term2
    end

    res = optimize(loss_func, dloss, input, BFGS())
    result_img = zeros(Float64, Nx, Nx)
    result_img = Optim.minimizer(res)
    return result_img

end
=#

#function optim_subroutine()

##In PROGRESS: Rewriting to use deriv_sums: In snippets_tbdebugged
function img_reconfunc_coeffmask(input, filter_hash, s_targ_mean, s_targ_invcov, pixmask, optim_settings; coeff_mask=nothing, coeff_type="S2")
    #=
    input: Initial image. When trues in pixmask, assumed to contain the correct value in those pixels.
    coeff_mask: Allows you to only optimize wrt a subset of coefficients (eg: j1<=j2)
    pixmask: Which pixels are floating (false) vs fixed (true)
    =#
    println("Coeff mask:", (coeff_mask!=nothing))
    (Nx, Ny)  = size(input)
    if Nx != Ny error("Input image must be square") end
    #filter_hash = fink_filter_hash(1, 8, nx=Nx, pc=1, wd=1, Omega=true)

    (Nf, )    = size(filter_hash["filt_index"])
    if Nf == 0 error("filter hash corrupted") end
    if coeff_type=="S2"
        println("S2")
        pixmask = pixmask[:] #flattened: Nx^2s
        cpvals = copy(input[:])[pixmask] #constrained pix values

        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 ((s_targ_invcov!=I) & size(s_targ_invcov)[1]!=count((i->(i==true)), coeff_mask)) error("s_targ_invcov should only contain coeffs to be optimized") end
            #if (s_targ_invcov!=I)  error("E1") end
            if size(s_targ_invcov)[1]!=count((i->(i==true)), coeff_mask) error("E2") end
        else #No mask: all coeffs (default: All S1 and S2 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]== (2+Nf+Nf^2)
                s_targ_mean = s_targ_mean[3:end]
            else
                error("s_targ_mean of wrong size")
            end
            if (s_targ_invcov!=I) & (size(s_targ_invcov)==(2+Nf+Nf^2, 2+Nf+Nf^2))
                s_targ_invcov = s_targ_invcov[3:end, 3:end]
            else
                error("s_targ_invcov of wrong size")
            end
            #At this point both have dims S1+S2
            #3:end
            coeff_mask = fill(true, 2+Nf+Nf^2)
            coeff_mask[1] = false
            coeff_mask[2] = false
        end
        #After this all cases have a coeffmask
        function loss_func(img_curr)
            #println("Img",size(img_curr))
            s_curr = Deriv_Utils.DHC(reshape(img_curr, (Nx, Nx)),  filter_hash, doS2=true, doS12=false, doS20=false)[coeff_mask]
            neglogloss = 0.5 .* (s_curr - s_targ_mean)' * s_targ_invcov * (s_curr - s_targ_mean)
            #println("In Loss:",size(s_curr))
            #println("NegLogLoss",neglogloss[1])
            return neglogloss[1]
        end

        function dloss(storage_grad, img_curr)
            s_curr = Deriv_Utils.DHC(reshape(img_curr, (Nx, Nx)), filter_hash, doS2=true, doS12=false, doS20=false)[coeff_mask]
            diff = s_curr - s_targ_mean
            #diff = diff
            dS1S2 = Transpose(reshape(Deriv_Utils.wst_S1S2_derivfast(reshape(img_curr, (Nx, Nx)), filter_hash), (:, Nx^2))) #(Nf+Nf^2, Nx, Nx)->(Nx^2, Nf+Nf^2) Cant directly do this without handling S1 coeffs
            dS1S2 = dS1S2[:, coeff_mask[3:end]] #Nx^2 * |SelCoeff|
            dS1S2[pixmask, :] .= 0 #Zeroing out wrt fixed params
            term1 = s_targ_invcov * diff #(Nf+Nf^2) or |SelCoeff| x 1
            #WARNING: Uses the Deriv_Utils version of dS1S2. Need to rewrite if you want to use the DHC_2DUtils one.
            #println("In dLoss:",size(diff), size(term1), size(term2))
            #println(size(term1),size(term2),size(storage_grad))
            mul!(storage_grad, dS1S2, term1) #Nx^2x1
        end
        #Debugging stuff
        println("Chisq derivative 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

        clever = reshape(zeros(size(input)), (Nx*Nx, 1))
        _bar = dloss(clever, reshape(input, (Nx^2, 1)))
        println("Brute:  ",brute)
        println("Clever: ",clever[Nx*(1)+1])

        res = optimize(loss_func, dloss, reshape(input, (Nx^2, 1)), ConjugateGradient(), Optim.Options(iterations = get(optim_settings, "iterations", 100), store_trace = true, show_trace = true))
        result_img = zeros(Float64, Nx, Nx)
        result_img = Optim.minimizer(res)
    elseif coeff_type=="S12"
        println("S12")
        pixmask = pixmask[:] #flattened: Nx^2s
        cpvals = copy(input[:])[pixmask] #constrained pix values

        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("E2") end
            #Add WARNING that flags when the S1 are additionally incl in coeffmask
        else #No mask: all coeffs (default: All S1 and S2 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]== (2+Nf+Nf^2)
                s_targ_mean = s_targ_mean[3:end]
            else
                error("s_targ_mean of wrong size")
            end
            if (s_targ_invcov!=I) & (size(s_targ_invcov)==(2+Nf+Nf^2, 2+Nf+Nf^2))
                s_targ_invcov = s_targ_invcov[3:end, 3:end]
            else
                error("s_targ_invcov of wrong size")
            end
            #At this point both have dims S1+S2
            #3:end
            coeff_mask = fill(true, 2+Nf+Nf^2)
            coeff_mask[1] = false
            coeff_mask[2] = false
        end
        #After this all cases have a coeffmask
        function loss_func12(img_curr)
            #println("Img",size(img_curr))
            s_curr = Deriv_Utils.DHC(reshape(img_curr, (Nx, Nx)),  filter_hash, doS2=false, doS12=true, doS20=false)[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)
            s_curr = Deriv_Utils.DHC(reshape(img_curr, (Nx, Nx)), filter_hash, doS2=false, doS12=true, doS20=false)[coeff_mask]
            diff = s_curr - s_targ_mean
            dS12 = reshape(permutedims(Deriv_Utils.wst_S12_deriv(reshape(img_curr, (Nx, Nx)), filter_hash), [3, 4, 1, 2]), (Nx^2, Nf^2)) #Nf, Nf, Nx, Nx -> (Nx^2, Nf^2)
            println(size(dS12))
            dS12 = dS12[:, coeff_mask[3+Nf:end]] #Nx^2 * |SelCoeff|
            dS12[pixmask, :] .= 0   #Zeroing out wrt fixed pix
            term1 = s_targ_invcov * diff #CHECK: Same as diff^Cinv? Nf^2 x1
            mul!(storage_grad, dS12, term1) #term1: 1x (Nf+Nf^2) | term2: #(Nf+Nf^2) x Nx^2
        end

        #Debugging stuff
        println("Chisq derivative 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

        clever = reshape(zeros(size(input)), (Nx*Nx, 1))
        _bar = dloss12(clever, reshape(input, (Nx^2, 1)))
        println("Brute:  ",brute)
        println("Clever: ",clever[Nx*(1)+1])

        res = optimize(loss_func12, dloss12, reshape(input, (Nx^2, 1)), ConjugateGradient(), Optim.Options(iterations = get(optim_settings, "iterations", 100), store_trace = true, show_trace = true))
        result_img = zeros(Float64, Nx, Nx)
        result_img = Optim.minimizer(res)
    else error("S20 not implemented")
    end
    return reshape(result_img, (Nx, Nx))
end
#Passed chisq test for all floating, and fixed

function img_reconfunc(input, filter_hash, s_targ_mean, s_targ_invcov, pixmask, optim_settings; coeff_type="S20")
    #=
    input: Square.
    s_targ_mean, s_targ_invcov should have the same format as the raw output of the appropriate DHC call. i.e. you need to exclude S0/S1 from s_targ_mean below.
    pixmask: Which pixels are floating (false) vs fixed (true)
    =#
    #Consider making this the wrapper for a subroutine that is constant irrespective of the choice of coeffs, mask etc. SO ths wrapper handles
    #all the choices of S2 vs S20 or only S1 etc--nope probably messier because would need to pass wrapper for appropriate type of DHC
    #
    (Nx, Ny)  = size(input)
    if Nx != Ny error("Input image must be square") end
    filter_hash = fink_filter_hash(1, 8, nx=Nx, pc=1, wd=1, Omega=true)

    (Nf, )    = size(filter_hash["filt_index"])
    if Nf == 0 error("filter hash corrupted") end

    #Handle optimizer settings
    #storage_beta = Array{Float64, 2}(undef, 1, Nx^2)
    #Handle reshaping -- make less flexi and change to asserts to make less messy
    #=if coeff_choice=="S2"
        #Doesn't work currently.
        println("S2")
        if size(s_targ_mean)[1]!=(Nf+Nf^2)
            s_targ_mean = s_targ_mean[3:end]
        end

        if (s_targ_invcov!=I) & (size(s_targ_invcov)!=(Nf+Nf^2, Nf+Nf^2))
            s_targ_invcov = s_targ_invcov[3:end, 3:end]
        end

        function loss_func(img_curr)
            #println("Img",size(img_curr))
            s_curr = Deriv_Utils.DHC(reshape(img_curr, (Nx, Nx)),  filter_hash, doS2=true, doS12=false, doS20=false)[3:end]
            neglogloss = 0.5 .* (s_curr - s_targ_mean)' * s_targ_invcov * (s_curr - s_targ_mean)
            #println("In Loss:",size(s_curr))
            #println("NegLogLoss",neglogloss[1])
            return neglogloss[1]
        end

        function dloss(storage_grad, img_curr)
            s_curr = Deriv_Utils.DHC(reshape(img_curr, (Nx, Nx)), filter_hash, doS2=true, doS12=false, doS20=false)[3:end]
            diff = s_curr - s_targ_mean
            #diff = diff
            dS1S2 = Transpose(reshape(Deriv_Utils.wst_S1S2_derivfast(reshape(img_curr, (Nx, Nx)), filter_hash), (:, Nx^2))) #(Nf+Nf^2, Nx, Nx)->(Nx^2, Nf+Nf^2) Cant directly do this without handling S1 coeffs
            term1 = s_targ_invcov * diff #(Nf+Nf^2) x 1
            #WARNING: Uses the Deriv_Utils version of dS1S2. Need to rewrite if you want to use the DHC_2DUtils one.
            #println("In dLoss:",size(diff), size(term1), size(term2))
            #println(size(term1),size(term2),size(storage_grad))
            mul!(storage_grad, dS1S2, term1) #Nx^2x1
            #storage_grad .=  term1 * term2
        end
        #Debugging stuff
        println("Chisq derivative 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

        clever = reshape(zeros(size(input)), (Nx*Nx, 1))
        _bar = dloss(clever, reshape(input, (Nx^2, 1)))
        println("Brute:  ",brute)
        println("Clever: ",clever[Nx*(1)+1])

        res = optimize(loss_func, dloss, reshape(input, (Nx^2, 1)), ConjugateGradient(), Optim.Options(iterations = get(optim_settings, "iterations", 100), store_trace = true, show_trace = true))
        result_img = zeros(Float64, Nx, Nx)
        result_img = Optim.minimizer(res)

    else #Currently S20
    =#
        #DO NOT RENAME THESE FUNCS TO loss_func and dloss_func: otherwise Julia calls the s20 version of the loss in the S2 case.
    if coeff_type=="S20"
        println("S20")
        pixmask = pixmask[:] #flattened: Nx^2s
        cpvals = copy(input[:])[pixmask] #constrained pix values
        if size(s_targ_mean)[1]== (2+Nf+Nf^2)
            s_targ_mean = s_targ_mean[3+Nf:end]
        else
            error("s_targ_mean of wrong size")
        end
        if (s_targ_invcov!=I) & (size(s_targ_invcov)==(2+Nf+Nf^2, 2+Nf+Nf^2))
            s_targ_invcov = s_targ_invcov[3+Nf:end, 3+Nf:end]
        else
            error("s_targ_invcov of wrong size")
        end
        #Subsequently s_targ_invcov and s_targ_mean have * only * the S20 coeffs

        #println("Shapes:s_targ_mean, s_targ_invcov", size(s_targ_mean), size(s_targ_invcov))
        function loss_func20(img_curr)
            s_curr = Deriv_Utils.DHC(reshape(img_curr, (Nx, Nx)),  filter_hash, doS2=false, doS12=false, doS20=true)[3+Nf:end] #not counting S1 twice
            neglogloss = 0.5 .* (s_curr - s_targ_mean)' * s_targ_invcov * (s_curr - s_targ_mean)
            #println("s_curr",size(s_curr))
            return neglogloss[1]
        end

        function dloss20(storage_grad, img_curr)
            #Diagonal covar: This works!
            #=
            s_curr = Deriv_Utils.DHC(reshape(img_curr, (Nx, Nx)), filter_hash, doS2=false, doS12=false, doS20=true)[3+Nf:end]
            diff = (s_curr - s_targ_mean).* Array(s_targ_invcov[diagind(s_targ_invcov)])
            dS20 = reshape(Deriv_Utils.wst_S20_deriv(reshape(img_curr, (Nx, Nx)), filter_hash), (Nx^2, Nf^2)) #Check (Nx, Nx, Nf, Nf) -> (Nx^2, Nf^2)
            mul!(storage_grad, dS20, diff) #term1: 1x (Nf+Nf^2) | term2: #(Nf+Nf^2) x Nx^2
            =#
            #More general:
            s_curr = Deriv_Utils.DHC(reshape(img_curr, (Nx, Nx)), filter_hash, doS2=false, doS12=false, doS20=true)[3+Nf:end]
            diff = s_curr - s_targ_mean
            dS20 = reshape(Deriv_Utils.wst_S20_deriv(reshape(img_curr, (Nx, Nx)), filter_hash), (Nx^2, Nf^2)) #Check (Nx, Nx, Nf, Nf) -> (Nx^2, Nf^2)
            dS20[pixmask, :] .= 0
            term1 = s_targ_invcov * diff #Nf^2 x1
            mul!(storage_grad, dS20, term1) #term1: 1x (Nf+Nf^2) | term2: #(Nf+Nf^2) x Nx^2

        end
        #
        println("Chisq derivative 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

        clever = reshape(zeros(size(input)), (Nx*Nx, 1))
        _bar = dloss20(clever, reshape(input, (Nx^2, 1)))
        println("Brute:  ",brute)
        println("Clever: ",clever[Nx*(1)+1])

        res = optimize(loss_func20, dloss20, reshape(input, (Nx^2, 1)), ConjugateGradient(), Optim.Options(iterations = get(optim_settings, "iterations", 100), store_trace = true, show_trace = true))
        result_img = zeros(Float64, Nx, Nx)
        result_img = Optim.minimizer(res)
    elseif coeff_type=="S12"
        println("S12")
        pixmask = pixmask[:] #flattened: Nx^2s
        cpvals = copy(input[:])[pixmask] #constrained pix values
        if size(s_targ_mean)[1]== (2+Nf+Nf^2)
            s_targ_mean = s_targ_mean[3+Nf:end]
        else
            error("s_targ_mean of wrong size")
        end
        if (s_targ_invcov!=I) & (size(s_targ_invcov)==(2+Nf+Nf^2, 2+Nf+Nf^2))
            s_targ_invcov = s_targ_invcov[3+Nf:end, 3+Nf:end]
        else
            error("s_targ_invcov of wrong size")
        end
        #Subsequently s_targ_invcov and s_targ_mean have * only * the S20 coeffs

        #println("Shapes:s_targ_mean, s_targ_invcov", size(s_targ_mean), size(s_targ_invcov))
        function loss_func12(img_curr)
            s_curr = Deriv_Utils.DHC(reshape(img_curr, (Nx, Nx)),  filter_hash, doS2=false, doS12=true, doS20=false)[3+Nf:end] #not counting S1 twice
            neglogloss = 0.5 .* (s_curr - s_targ_mean)' * s_targ_invcov * (s_curr - s_targ_mean)
            #println("s_curr",size(s_curr))
            return neglogloss[1]
        end

        function dloss12(storage_grad, img_curr)
            #Diagonal covar: This works!
            #=
            s_curr = Deriv_Utils.DHC(reshape(img_curr, (Nx, Nx)), filter_hash, doS2=false, doS12=false, doS20=true)[3+Nf:end]
            diff = (s_curr - s_targ_mean).* Array(s_targ_invcov[diagind(s_targ_invcov)])
            dS20 = reshape(Deriv_Utils.wst_S20_deriv(reshape(img_curr, (Nx, Nx)), filter_hash), (Nx^2, Nf^2)) #Check (Nx, Nx, Nf, Nf) -> (Nx^2, Nf^2)
            mul!(storage_grad, dS20, diff) #term1: 1x (Nf+Nf^2) | term2: #(Nf+Nf^2) x Nx^2
            =#
            #More general:
            s_curr = Deriv_Utils.DHC(reshape(img_curr, (Nx, Nx)), filter_hash, doS2=false, doS12=true, doS20=false)[3+Nf:end]
            diff = s_curr - s_targ_mean
            dS12 = reshape(permutedims(Deriv_Utils.wst_S12_deriv(reshape(img_curr, (Nx, Nx)), filter_hash), [3, 4, 1, 2]), (Nx^2, Nf^2)) #Nf, Nf, Nx, Nx -> (Nx^2, Nf^2)
            dS12[pixmask, :] .= 0
            term1 = s_targ_invcov * diff #CHECK: Same as diff^Cinv? Nf^2 x1
            mul!(storage_grad, dS12, term1) #term1: 1x (Nf+Nf^2) | term2: #(Nf+Nf^2) x Nx^2

        end
        #
        println("Chisq derivative 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

        clever = reshape(zeros(size(input)), (Nx*Nx, 1))
        _bar = dloss12(clever, reshape(input, (Nx^2, 1)))
        println("Brute:  ",brute)
        println("Clever: ",clever[Nx*(1)+1])

        res = optimize(loss_func12, dloss12, reshape(input, (Nx^2, 1)), ConjugateGradient(), Optim.Options(iterations = get(optim_settings, "iterations", 100), store_trace = true, show_trace = true))
        result_img = zeros(Float64, Nx, Nx)
        result_img = Optim.minimizer(res)

    else error("Not Implemented S2 here")
    end

    return reshape(result_img, (Nx, Ny))

end

function readdust()

    RGBA_img = load(pwd()*"/scratch_DF/t115_clean_small.png")
    img = reinterpret(UInt8,rawview(channelview(RGBA_img)))

    return img[1,:,:]
end

function S2_weights(im, fhash, Nsam=10; iso=iso)

    (Nx, Ny)  = size(im)
    if Nx != Ny error("Input image must be square") end
    (N1iso, Nf)    = size(fhash["S1_iso_mat"])

    # fhash = fink_filter_hash(1, 8, nx=Nx, pc=1, wd=1)

    S2   = DHC_compute(im, fhash, doS2=true, doS20=false, norm=false, iso=iso)
    Ns    = length(S2)
    S2arr = zeros(Float64, Ns, Nsam)
    println("Ns", Ns)
    for j=1:Nsam
        noise = rand(Nx,Nx)
        S2arr[:,j] = DHC_compute(im+noise, fhash, doS2=true, doS20=false, norm=false, iso=iso)
    end
    wt = zeros(Float64, Ns)
    for i=1:Ns
        wt[i] = std(S2arr[i,:])
    end
    msub = S2arr .- mean(S2arr, dims=2)
    cov = (msub * msub')./(Nsam-1)

    return wt, cov
end
####Comparison Functions#########
#Basically sandbox_DF funcs
function S20_weights(im, fhash, Nsam=10; iso=iso)

    (Nx, Ny)  = size(im)
    if Nx != Ny error("Input image must be square") end
    (N1iso, Nf)    = size(fhash["S1_iso_mat"])

    # fhash = fink_filter_hash(1, 8, nx=Nx, pc=1, wd=1)

    S20   = DHC_compute(im, fhash, doS2=false, doS20=true, norm=false, iso=iso)
    Ns    = length(S20)
    S20arr = zeros(Float64, Ns, Nsam)
    for j=1:Nsam
        noise = rand(Nx,Nx)
        S20arr[:,j] = DHC_compute(im+noise, fhash, doS2=false, doS20=true, norm=false, iso=iso)
    end

    wt = zeros(Float64, Ns)
    for i=1:Ns
        wt[i] = std(S20arr[i,:])
    end
    msub = S20arr .- mean(S20arr, dims=2)
    cov = (msub * msub')./(Nsam-1)

    return wt, cov
end

function S12_weights(im, fhash, Nsam=10; iso=iso)

    (Nx, Ny)  = size(im)
    if Nx != Ny error("Input image must be square") end
    (N1iso, Nf)    = size(fhash["S1_iso_mat"])

    # fhash = fink_filter_hash(1, 8, nx=Nx, pc=1, wd=1)

    S12   = DHC_compute(im, fhash, doS2=false, doS20=false, doS12=true, norm=false, iso=iso)
    Ns    = length(S12)
    S12arr = zeros(Float64, Ns, Nsam)
    for j=1:Nsam
        noise = rand(Nx,Nx)
        S12arr[:,j] = DHC_compute(im+noise, fhash, doS2=false, doS20=false, doS12=true, norm=false, iso=iso)
    end

    wt = zeros(Float64, Ns)
    for i=1:Ns
        wt[i] = std(S12arr[i,:])
    end
    msub = S12arr .- mean(S12arr, dims=2)
    cov = (msub * msub')./(Nsam-1)

    return wt, cov
end

function wst_synthS20(im_init, fixmask, S_targ, S20sig; iso=false)
    # fixmask -  0=float in fit   1=fixed

    function wst_synth_chisq(vec_in)

        thisim = copy(im_init)
        thisim[indfloat] = vec_in

        M20 = fhash["S2_iso_mat"]

        S20 = DHC_compute(thisim, fhash, doS2=false, doS20=true, norm=false, iso=iso)
        i0 = 3+(iso ? N1iso : Nf)
        diff  = ((S20[i0:end] - S_targ)./S20sig)

        # should have some kind of weight here
        chisq = diff'*diff
        println(chisq)
        return chisq

    end

    function wst_synth_dchisq(storage, vec_in)

        thisim = copy(im_init)
        thisim[indfloat] = vec_in

        dS20dp = reshape(wst_S20_deriv(thisim, fhash), Nx*Nx, Nf*Nf)
        if iso
            M20 = fhash["S2_iso_mat"]
            dS20dp = dS20dp * M20'
        end
        i0 = 3+(iso ? N1iso : Nf)

        S20arr = (DHC_compute(thisim, fhash, doS2=false, doS20=true, norm=false, iso=iso))[i0:end]

        # put both factors of S20sig in this array to weight
        diff   = (S20arr - S_targ)./(S20sig.^2)
        #println("size of diff", size(diff))
        #println("size of dS20dp", size(reshape(dS20dp, Nx*Nx, Nf*Nf)))
        # dSdp matrix * S1-S_targ is dchisq
        #dchisq_im = (reshape(dS20dp, Nx*Nx, Nf*Nf) * diff).*2
        dchisq_im = (dS20dp * diff).*2
        dchisq = reshape(dchisq_im, Nx, Nx)[indfloat]

        storage .= dchisq
    end

    (Nx, Ny)  = size(im_init)
    if Nx != Ny error("Input image must be square") end
    (N1iso, Nf)    = size(fhash["S1_iso_mat"])

    # index list of pixels to float in fit
    indfloat = findall(fixmask .== 0)

    # initial values for floating pixels
    vec_init = im_init[indfloat]
    println(length(vec_init), " pixels floating")

    eps = zeros(size(vec_init))
    eps[1] = 1e-4
    chisq1 = wst_synth_chisq(vec_init+eps./2)
    chisq0 = wst_synth_chisq(vec_init-eps./2)
    brute  = (chisq1-chisq0)/1e-4

    clever = zeros(size(vec_init))
    _bar = wst_synth_dchisq(clever, vec_init)
    println("Brute:  ",brute)
    println("Clever: ",clever[1])

    # call optimizer
    #res = optimize(wst_synth_chisq, wst_synth_dchisq, copy(vec_init), BFGS())
    res = optimize(wst_synth_chisq, wst_synth_dchisq, copy(vec_init),
        ConjugateGradient(), Optim.Options(iterations=100))

    # copy results into pixels of output image
    im_synth = copy(im_init)
    im_synth[indfloat] = Optim.minimizer(res)
    println(res)

    return im_synth
end

function plot_synth_QA(ImTrue, ImInit, ImSynth, fhash; fname="test2.png")

    # -------- define plot1 to append plot to a list
    function plot1(ps, image; clim=nothing, bin=1.0, fsz=16, label=nothing)
        # ps    - list of plots to append to
        # image - image to heatmap
        # clim  - tuple of color limits (min, max)
        # bin   - bin down by this factor, centered on nx/2+1
        # fsz   - font size
        marg   = 1mm
        nx, ny = size(image)
        nxb    = nx/round(Integer, 2*bin)

        # -------- center on nx/2+1
        i0 = max(1,round(Integer, (nx/2+2)-nxb-1))
        i1 = min(nx,round(Integer, (nx/2)+nxb+1))
        lims = [i0,i1]
        subim = image[i0:i1,i0:i1]
        push!(ps, heatmap(image, aspect_ratio=:equal, clim=clim,
            xlims=lims, ylims=lims, size=(400,400),
            legend=false, xtickfontsize=fsz, ytickfontsize=fsz,#tick_direction=:out,
            rightmargin=marg, leftmargin=marg, topmargin=marg, bottommargin=marg))
        if label != nothing
            annotate!(lims'*[.96,.04],lims'*[.09,.91],text(label,:left,:white,32))
        end
        return
    end

    # -------- initialize array of plots
    ps   = []
    clim  = (0,200)
    clim2 = (0,200).-100

    # -------- 6 panel QA plot
    plot1(ps, ImTrue, clim=clim, label="True")
    plot1(ps, ImSynth, clim=clim, label="Synth")
    plot1(ps, ImInit, clim=clim, label="Init")
    plot1(ps, ImInit-ImTrue, clim=clim2, label="Init-True")
    plot1(ps, ImInit-ImSynth, clim=clim2, label="Init-Synth")
    plot1(ps, ImSynth-ImTrue, clim=clim2, label="Synth-True")

    myplot = plot(ps..., layout=(3,2), size=(1400,2000))
    savefig(myplot, fname)
end
############################################
#Code Checks: S2: Passes chisq test (without coeff_mask)
epsilon=1e-5
dust = Float64.(readdust())
dust = dust[1:256,1:256]

Nx     = 32
doiso  = false
fhash = fink_filter_hash(1, 8, nx=Nx, pc=1, wd=1, Omega=true)
(N1iso, Nf)    = size(fhash["S1_iso_mat"])
i0 = 3+(doiso ? N1iso : Nf)
img    = imresize(dust,(Nx,Nx))
fixmask = rand(Nx,Nx) .< 0.1
#fixmask = falses(Nx, Nx)
#Sanity check: fixmask = trues(Nx, Nx). Was 0.
init = copy(img)
floatind = findall(fixmask .==0)
init[floatind] .+= rand(length(floatind)).*50 .-25

s2targ = DHC_2DUtils.DHC_compute(img, fhash, doS2=true, doS20=false, norm=false, iso=doiso)
s2w, cov = S2_weights(img, fhash, 100, iso=doiso)
s2w[findall(s2w .< epsilon)] .= s2w[findall(s2w .< epsilon)] .+ 1e-5

invcovarS2 = Diagonal(s2w.^(-2))
#Profile.clear()
reconS2cgeps_check = img_reconfuncS2(init, fhash, s2targ,  invcovarS2, fixmask, Dict([("iterations", 10)]))

#Code Checks: S2 (With coeffmask)
epsilon=1e-5
dust = Float64.(readdust())
dust = dust[1:256,1:256]
img = dust[1:16, 1:16]
Nx=16
Omega_bool = true
test = img[1:Nx, 1:Nx]
fhash = fink_filter_hash(1, 8, nx=Nx, pc=1, wd=1, Omega=Omega_bool)
filt64, f64info = fink_filter_bank(1, 8, nx=Nx, wd=1, pc=1, shift=false, Omega=true, safety_on=true)
s2targ = DHC_2DUtils.DHC_compute(img, fhash, doS2=true, doS20=false, norm=false, iso=doiso)
Nf = length(fhash["filt_index"])

#Mask cons
J=log2(Nx) - 2
L = 8
Nf = length(fhash["filt_index"])
#finoffx = Omega_bool ? 2: 1
#j1 = reshape(f64info["J_L"][1:(end-finoffx)], ())
fjvals = f64info["J_L"][:, 1]
outerjf1f2 =Array{Tuple{Int, Int}}(undef, length(fjvals), length(fjvals))
f1f2idx =Array{Tuple{Int, Int}}(undef, length(fjvals), length(fjvals))
for f1=1:Nf
    for f2=1:Nf
        outerjf1f2[f1, f2] = (fjvals[f1], fjvals[f2])
        f1f2idx[f1, f2] = (f1, f2)
    end
end
f1f2idx_flat = f1f2idx[:]
bool_incl = (tu->(tu[1]<=tu[2])).(outerjf1f2)
booltot_flat = cat([false, false], fill(true, Nf), bool_incl[:], dims=1) #false if we don't wanna optimize
largemask = s2targ .> 1e-5

s2w, covs2 = S2_weights(img, fhash, 100, iso=doiso)
#maximum(s2w[booltot_flat]), minimum(s2w[booltot_flat])
#(3604, 0.002) Without booltot, the minimum is 1e-28

invcovarS2 = Diagonal(s2w.^(-2))
invcovarS2 = invcovarS2[booltot_flat, booltot_flat]
img_s2restr = img_reconfunc_coeffmask(init, fhash, s2targ[booltot_flat], invcovarS2, fixmask, Dict([("iterations", 10)]), coeff_mask = booltot_flat,
coeff_type="S2")

#S12
#Code Checks: S12 (With coeffmask)
epsilon=1e-5
dust = Float64.(readdust())
dust = dust[1:256,1:256]
Nx=32
img = dust[1:Nx, 1:Nx]

Omega_bool = true
fhash = fink_filter_hash(1, 8, nx=Nx, pc=1, wd=1, Omega=Omega_bool)
Nf = length(fhash["filt_index"])
#fixmask = rand(Nx,Nx) .< 0.1
noise = reshape(rand(Normal(0, std(img)), Nx^2), (Nx, Nx))
init = img+ noise
init = imfilter(init, Kernel.gaussian(1))
#Mask cons
s12w, covs12 = S12_weights(img, fhash, 100, iso=doiso)
incl= (s12w.!=0)
invcovar = Diagonal(s12w.^(-2))
s12targ = DHC_2DUtils.DHC_compute(img, fhash, doS2=false, doS20=false, doS12=true, norm=false, iso=doiso)
#myfoogens20_100itn = img_reconfunc_old(init, s20targ,  invcovar, fixmask, "S20", Dict([("iterations", 100)]))
incl[1:Nf+2] .= 0 #removing S0 and redundant S1
S12recon = img_reconfunc_coeffmask(init, fhash, s12targ[incl], invcovar[incl, incl], fixmask, Dict([("iterations", 100)]), coeff_mask = incl,
 coeff_type="S12")
heatmap(img)
heatmap(init)
heatmap(S12recon)
#plot_synth_QA(img, init, S12recon, fhash, fname="scratch_NM/TestPlots/s12_100itns.png")
pl12 = plot(
    heatmap(img, title="Ground Truth"),
    heatmap(init, title="GT+N(0, std(I))"),
    heatmap(S12recon,title= "Reconstruction w S12"),
    heatmap(S12recon - img, title="Residual");
    layout=4,
)
plot!(size=(1000, 800))
savefig(pl12, "scratch_NM/TestPlots/s12_100itns_smoothedwhitenoiseinit.png")
#Without smoothing
mean(abs.(S12recon - img)) #2.54
mean(abs.((S12recon - img)./img)) #0.03
#=
invcovarS2 = Diagonal(s2w.^(-2))
invcovarS2 = invcovarS2[booltot_flat, booltot_flat]
img_s2restr = img_reconfunc_coeffmask(init, fhash, s2targ[booltot_flat], invcovarS2, fixmask, Dict([("iterations", 10)]), coeff_mask = booltot_flat,
coeff_type="S2")
=#

#Code Checks: S20. Works
s20w, covs20 = S20_weights(img, fhash, 100, iso=doiso)
invcovar = Diagonal(s20w.^(-2))
s20targ = DHC_2DUtils.DHC_compute(img, fhash, doS2=false, doS20=true, norm=false, iso=doiso)
#myfoogens20_100itn = img_reconfunc_old(init, s20targ,  invcovar, fixmask, "S20", Dict([("iterations", 100)]))

S20recon = img_reconfunc(init, fhash, s20targ, invcovar, fixmask, Dict([("iterations", 10)]), coeff_type="S20")

#Code Checks: S12

####Testing Examples#############
#Experiment 1: all pixels floating,
#Sec1: Don't use this--this was before I fixed the plotting / function bugs--skip to Sec2!!!
#=
dust = readdust()
dust = Array{Float64}(dust[1:16, 1:16])

n = Normal()
noise = reshape(rand(n, 256), (16, 16))
noisyinp = dust+noise

fhash = fink_filter_hash(1, 8, nx=16, pc=1, wd=1)
s_targ = DHC_2DUtils.DHC_compute(dust, fhash, doS2=false, doS20=true)[3+length(fhash["filt_index"]):end]
Nx=16
#noisyinp = zeros(Float64, Nx, Nx)
#noisyinp .= dust
#noisyinp[3, 5] += dust[3, 5] + 8
denoised_s20 = img_reconfunc(noisyinp, s_targ, I, trues(16, 16), "S20", Dict([("iterations", 1000)]))

#Comparing with wst_synthS20
denoised_comp = wst_synthS20(noisyinp, falses((Nx, Nx)), s_targ, 0.2)
#Doug's test code
Nx=16
dust = Float64.(readdust())
dust = dust[1:256,1:256]
img   = imresize(dust,(Nx,Nx))
n = reshape(rand(Normal(0, mean(img)/20.0), Nx^2), (Nx, Nx))
noisyimg = img + n
s_targ = DHC_2DUtils.DHC_compute(img, fhash, doS2=false, doS20=true)[3+length(fhash["filt_index"]):end]
denoised_img = wst_synthS20(noisyimg, falses(Nx,Nx), s_targ, mean(img)/20.0)
=#
#For some reason the above code doesnt work.
#Sec2: This now works!
#Using the exact same example as in sandbox_DFcopy.
dust = Float64.(readdust())
dust = dust[1:256,1:256]

Nx     = 16
doiso  = false
fhash = fink_filter_hash(1, 8, nx=Nx, pc=1, wd=1, Omega=true)
(N1iso, Nf)    = size(fhash["S1_iso_mat"])
i0 = 3+(doiso ? N1iso : Nf)
img    = imresize(dust,(Nx,Nx))
fixmask = rand(Nx,Nx) .< 0.5
#fixmask = falses((Nx, Nx))

init = copy(img)
floatind = findall(fixmask .==0)
init[floatind] .+= rand(length(floatind)).*50 .-25

S20sig, cov = S20_weights(img, fhash, 100, iso=doiso)
#S20sig = S20sig[i0:end]

#2.1. Using wst_synthS20
#foo = wst_synthS20(init, fixmask, S_targ, S20sig, iso=doiso)

#2.2. Using your code.
#S20 ---------------------------------
s20w, covs20 = S20_weights(img, fhash, 100, iso=doiso)
invcovar = Diagonal(s20w.^(-2))
s20targ = DHC_2DUtils.DHC_compute(img, fhash, doS2=false, doS20=true, norm=false, iso=doiso)
#myfoogens20_100itn = img_reconfunc_old(init, s20targ,  invcovar, fixmask, "S20", Dict([("iterations", 100)]))

S20recon = img_reconfuncS20(init, fhash, s20targ, invcovar, fixmask, Dict([("iterations", 10)]))

#Need a new invcovar (S2weights) for doing S2
s2w, covs2 = S2_weights(img, fhash, 100, iso=doiso)
invcovarS2 = Diagonal(s2w.^(-2)) #This doesn't work for S2: use invcovar computed for S20 for S2 as well for the timebeing.
s2targ = DHC_2DUtils.DHC_compute(img, fhash, doS2=true, doS20=false, norm=false, iso=doiso)
println(size(img), size(init), size(s2targ), size(invcovar), size(fixmask))
reconS2cg_ws20 = img_reconfunc(init, s2targ,  invcovar, fixmask, "S2", Dict([("iterations", 100)]))
#Weird bug where the code works using invcovar but doesnt when you use invcovars2
#myfoogen used S20's invcovar
#myfoogenS2 uses the correct incovar
#myfoogens20 uses s20

plot_synth_QA(img, init, myfoogens20, fhash, fname="scratch_NM/TestPlots/s20_100itns.png")
plot_synth_QA(img, init, myfoogenS2, fhash, fname="scratch_NM/TestPlots/s2_100itns.png") #Used LBFGS: slightly better than CG but much slower
plot_synth_QA(img, init, reconS2cg, fhash, fname="scratch_NM/TestPlots/s2_100itns_cg.png")

#Concl: julia> maximum(s2siginvsq), minimum(s2siginvsq)
#(1.418445018807988e55, 1.066753479196218e-7)

#julia> maximum(s20siginvsq), minimum(s20siginvsq)
#(45474.169557552836, 7.544983824162472e-8)
#This is why it doesnt work with S2 weights:

#Sec2.3 Adding an eps to s2w
epsilon=1e-5
dust = Float64.(readdust())
dust = dust[1:256,1:256]

Nx     = 64
doiso  = false
fhash = fink_filter_hash(1, 8, nx=Nx, pc=1, wd=1, Omega=true)
(N1iso, Nf)    = size(fhash["S1_iso_mat"])
i0 = 3+(doiso ? N1iso : Nf)
img    = imresize(dust,(Nx,Nx))
#fixmask = rand(Nx,Nx) .< 0.1
fixmask = falses((Nx, Nx))

init = copy(img)
floatind = findall(fixmask .==0)
init[floatind] .+= rand(length(floatind)).*50 .-25

s2targ = DHC_2DUtils.DHC_compute(img, fhash, doS2=true, doS20=false, norm=false, iso=doiso)
s2w, cov = S2_weights(img, fhash, 100, iso=doiso)
s2w[findall(s2w .< epsilon)] .= s2w[findall(s2w .< epsilon)] .+ 1e-5

invcovarS2 = Diagonal(s2w.^(-2))
reconS2cgeps_100itn = img_reconfuncS2(init, fhash, s2targ,  invcovarS2, fixmask, Dict([("iterations", 100)]))
#reconS2cgeps_new = img_reconfuncS2(init, fhash, s2targ, invcovarS2, fixmask, Dict([("iterations", 100)]))
#Sec2.4: j1<=j2+Using the new S2 specific function
#2.4.1: Without j1<=j2: Check
Nx=64
Omega_bool = true
test = img[1:Nx, 1:Nx]
fhash64 = fink_filter_hash(1, 8, nx=Nx, pc=1, wd=1, Omega=Omega_bool)
filt64, f64info = fink_filter_bank(1, 8, nx=Nx, wd=1, pc=1, shift=false, Omega=true, safety_on=true)
s2targ = DHC_2DUtils.DHC_compute(img, fhash, doS2=true, doS20=false, norm=false, iso=doiso)
Nf = length(fhash64["filt_index"])

#Mask cons
J=log2(Nx) - 2
L = 8
Nf = length(fhash64["filt_index"])
#finoffx = Omega_bool ? 2: 1
#j1 = reshape(f64info["J_L"][1:(end-finoffx)], ())
fjvals = f64info["J_L"][:, 1]
outerjf1f2 =Array{Tuple{Int, Int}}(undef, length(fjvals), length(fjvals))
f1f2idx =Array{Tuple{Int, Int}}(undef, length(fjvals), length(fjvals))
for f1=1:Nf
    for f2=1:Nf
        outerjf1f2[f1, f2] = (fjvals[f1], fjvals[f2])
        f1f2idx[f1, f2] = (f1, f2)
    end
end
f1f2idx_flat = f1f2idx[:]
bool_incl = (tu->(tu[1]<=tu[2])).(outerjf1f2)
booltot_flat = cat([false, false], fill(true, Nf), bool_incl[:], dims=1) #false if we don't wanna optimize
largemask = s2targ .> 1e-5

s2w, covs2 = S2_weights(img, fhash, 100, iso=doiso)
#maximum(s2w[booltot_flat]), minimum(s2w[booltot_flat])
#(3604, 0.002) Without booltot, the minimum is 1e-28

invcovarS2 = Diagonal(s2w.^(-2))
invcovarS2 = invcovarS2[booltot_flat, booltot_flat]
img_s2restr = img_reconfuncS2(init, fhash, s2targ[booltot_flat], invcovarS2, fixmask, Dict([("iterations", 100)]), coeff_mask = booltot_flat)

img_s2restr1k = img_reconfuncS2(init, fhash, s2targ[booltot_flat], invcovarS2, fixmask, Dict([("iterations", 1000)]), coeff_mask = booltot_flat)

#mean(abs.(reconS2cgeps_1000itn - img)), mean(abs.(reshape(img_s2restr1k, (Nx, Nx)) - img)),mean(abs.(myfoogens20_1000itn - img))
#(5.20568411412877, 5.239126791722001, 4.973186001871807)
#mean(abs.(reconS2cgeps_1000itn - img)), mean(abs.(reshape(img_s2restr1k, (Nx, Nx)) - img)),mean(abs.(myfoogens20_1000itn - img))
#(5.20568411412877, 5.239126791722001, 4.973186001871807)

#Experiment 2: Non-informative init: Doesn't work
dust = Float64.(readdust())
dust = dust[1:256,1:256]

Nx     = 64
doiso  = false
fhash = fink_filter_hash(1, 8, nx=Nx, pc=1, wd=1, Omega=true)
(N1iso, Nf)    = size(fhash["S1_iso_mat"])
i0 = 3+(doiso ? N1iso : Nf)
img    = imresize(dust,(Nx,Nx))
#fixmask = rand(Nx,Nx) .< 0.1
fixmask = falses((Nx, Nx))

init = copy(img)
floatind = findall(fixmask .==0)
init[floatind] .+= rand(length(floatind)).*50 .-25

S20sig = S20_weights(img, fhash, 100, iso=doiso)
S20sig = S20sig[i0:end]

invcovar = Diagonal(S20_weights(img, fhash, 100, iso=doiso).^(-2))
s20targ = DHC_2DUtils.DHC_compute(img, fhash, doS2=false, doS20=true, norm=false, iso=doiso)
myfoogens20 = img_reconfunc(fill(mean(img), (Nx, Nx)), s20targ,  invcovar, fixmask, "S20", Dict([("iterations", 1000)]))

#mean(abs.(myfoogens20_500itn - img)), mean(abs.(reconS2cgeps_500itn - img))
#(5.11, 5.82)
#Q2: S2, S20 differences
#mean(s2wsel), median(s2wsel), mean(s2wsel./s2sel), median(s2wsel./s2sel)
#(20.845713040948258, 2.2510762633937498, 0.009551723133768989, 0.008862903344736044)
#mean(s20w), median(s20w), mean(s20w./s20targ), median(s20w./s20targ)
#(137.99458258542106, 70.5409237758934, 0.0016749817200681092, 0.0016540522310717976)
#mean(s2w), median(s2w), mean(s2w./s2targ), median(s2w./s2targ)
#(13.192017376387687, 0.6169549962485834, 0.019359754705979722, 0.009315575046683549)

#Juno.profiler()
#Profile.print(format=:tree, sortedby=:count)

#why is this so slow for some reason?? Takes 80s for 100 itns with all free but now takes 80s for 10itns
#Check
#reconS2cgeps_10pct[fixmask] .== img[fixmask]

#Experiment 3 (Code checks shifted up)

#3.1. Using 64 different parts of the same larger dust img for the 32^2 case
epsilon = 1e-5
doiso=false
heatmap(dust)
Nx=32
LargeNx = 256
numtiles = convert(Int8, LargeNx/Nx)
dbnimg = []
x_end = collect(1:1:256)
x_end = x_end[mod.(x_end, Nx).==0]
for rid in x_end
    for cid in x_end
        push!(dbnimg, dust[rid-Nx+1:rid, cid-Nx+1:cid])
    end
end
dbnimg = (x->reshape(x, (Nx, Nx))).(dbnimg)
fhash = fink_filter_hash(1, 8, nx=Nx, pc=1, wd=1, Omega=true)
Nf = length(fhash["filt_index"])
s2dbn = (x->DHC_2DUtils.DHC_compute(x, fhash, doS2=true, doS20=false, norm=false, iso=doiso)).(dbnimg[2:end])
s2dbn = hcat(s2dbn...)
s2mean = mean(s2dbn, dims=2)
s2sig = std(s2dbn, dims=2)
s2sig[findall(s2sig.<epsilon)] .+= epsilon
s2sig = s2sig[:]
invcovar = Diagonal(s2sig.^(-2))

heatmap(dbnimg[1], title="Ground Truth")

fixmask = rand(Nx,Nx) .< 0.5
#fixmask = falses(Nx, Nx)
#Sanity check: fixmask = trues(Nx, Nx). Was 0.
init = copy(dbnimg[34])
floatind = findall(fixmask .==0)
init[fixmask] .= dbnimg[1][fixmask]
#heatmap(init, title="Init10pct")

dbns2_50pct = img_reconfuncS2(init, fhash, s2mean,  invcovar, fixmask, Dict([("iterations", 100)]))
maximum(s2sig), minimum(s2sig)

plot_synth_QA(dbnimg[1], init, dbns2_10pct, fhash, fname="scratch_NM/TestPlots/dbns2_100itns_50pct.png")
heatmap(init, title="Init 10pct")
heatmap(dbns2_50pct, title="Dbn10pct Reconstruction")

plot(
    heatmap(dbnimg[1], title="Ground Truth"),
    heatmap(init, title="Init 50pct"),
    heatmap(dbns2_50pct,title= "Reconstruction w S2dbn"),
    heatmap(dbns2_50pct - dbnimg[1], title="Residual");
    layout=4
)
#heatmap(dbns2_50pct, title="Reconstruction gn 50% pixels")

strue = DHC_2DUtils.DHC_compute(dbnimg[1], fhash, doS2=true, doS20=false, norm=false, iso=doiso)
sres = DHC_2DUtils.DHC_compute(dbns2_50pct, fhash, doS2=true, doS20=false, norm=false, iso=doiso)
strue[1], strue[2], sres[1], sres[2]
#gets the variance wrong -- maybe constrain this too?
sres[3:3+Nf] - strue[3:3+Nf], mean((sres[3:3+Nf] - strue[3:3+Nf])./strue[3:3+Nf])
mean((s2mean[3:3+Nf] - strue[3:3+Nf])./strue[3:3+Nf])
heatmap(reshape(log.(strue[3+Nf:end]), (Nf, Nf)), title="true S2")
heatmap(reshape(log.(sres[3+Nf:end]), (Nf, Nf)), title="result S2")
heatmap(reshape(log.(abs.((sres[3+Nf:end] - strue[3+Nf:end])./strue[3+Nf:end])), (Nf, Nf)), title="Log Abs (S2res - S2true) / S2true")

#Plotting: Scatter
trace1 = Plotly.scatter(;x=reshape(img, (Nx^2)), y=reshape(myfoogens20_500itn, (Nx^2)), mode="markers", name="S20 Reconstruction")
trace2 = Plotly.scatter(;x=reshape(img, (Nx^2)), y=reshape(reconS2cgeps_500itn, (Nx^2)), mode="markers", name="S2 Reconstruction (w eps)")
data = [trace1, trace2]
layout = Layout(;title="Reconstruction Value vs True Value")
response = Plotly.plot(data, layout)

#plot_url = response["url"]
#Experiment 4: With S12
#Code check above
#Now, with log of image
dust = Float64.(readdust())
dust = dust[1:256,1:256]
Nx=32
img = dust[1:Nx, 1:Nx]
oriimg = copy(img)
img = log.(img)

Omega_bool = true
fhash = fink_filter_hash(1, 8, nx=Nx, pc=1, wd=1, Omega=Omega_bool)
Nf = length(fhash["filt_index"])
#fixmask = rand(Nx,Nx) .< 0.1
noise = reshape(rand(Normal(0, std(oriimg)), Nx^2), (Nx, Nx))
init = oriimg+ noise
init = imfilter(init, Kernel.gaussian(0.8))
init = log.(init)
#Mask cons
s12w, covs12 = S12_weights(img, fhash, 100, iso=doiso)
incl= (s12w.!=0)
invcovar = Diagonal(s12w.^(-2))
s12targ = DHC_2DUtils.DHC_compute(img, fhash, doS2=false, doS20=false, doS12=true, norm=false, iso=doiso)
#myfoogens20_100itn = img_reconfunc_old(init, s20targ,  invcovar, fixmask, "S20", Dict([("iterations", 100)]))
incl[1:Nf+2] .= 0 #removing S0 and redundant S1
S12recon = img_reconfunc_coeffmask(init, fhash, s12targ[incl], invcovar[incl, incl], fixmask, Dict([("iterations", 100)]), coeff_mask = incl,
 coeff_type="S12")
heatmap(img)
heatmap(init)
heatmap(S12recon)
#plot_synth_QA(img, init, S12recon, fhash, fname="scratch_NM/TestPlots/s12_100itns.png")
pl12 = plot(
    heatmap(oriimg, title="I_GT: Ground Truth"),
    heatmap(exp.(init), title="I_GT+N(0, std(I_GT))"),
    heatmap(exp.(S12recon),title= "Reconstruction w S12"),
    heatmap(exp.(S12recon) - oriimg, title="Residual");
    layout=4,
)
plot!(size=(1000, 800))
savefig(pl12, "scratch_NM/TestPlots/s12_100itns_logIgd_smoothgk0-8.png")
#Without smoothing
mean(abs.(exp.(S12recon) - oriimg)) #2.54
mean(abs.((exp.(S12recon) - oriimg)./oriimg)) #0.03

####S2#########
epsilon=1e-5
dust = Float64.(readdust())
dust = dust[1:256,1:256]
img = dust[1:Nx, 1:Nx]
doiso  = false
fhash = fink_filter_hash(1, 8, nx=Nx, pc=1, wd=1, Omega=true)
(N1iso, Nf)    = size(fhash["S1_iso_mat"])
i0 = 3+(doiso ? N1iso : Nf)
img    = imresize(dust,(Nx,Nx))
fixmask = rand(Nx,Nx) .< 0.1
fixmask = falses(Nx, Nx)
#Sanity check: fixmask = trues(Nx, Nx). Was 0.
noise = reshape(rand(Normal(0, std(oriimg)), Nx^2), (Nx, Nx))
init = oriimg+ noise

s2targ = DHC_2DUtils.DHC_compute(img, fhash, doS2=true, doS20=false, norm=false, iso=doiso)
s2w, cov = S2_weights(img, fhash, 100, iso=doiso)
s2w[findall(s2w .< epsilon)] .= s2w[findall(s2w .< epsilon)] .+ 1e-5

invcovarS2 = Diagonal(s2w.^(-2))
#Profile.clear()
reconS2cgeps_check = img_reconfunc_coeffmask(init, fhash, s2targ,  invcovarS2, fixmask, Dict([("iterations", 100)]))
pl12 = plot(
    heatmap(img, title="I_GT: Ground Truth"),
    heatmap(init, title="I_GT+N(0, std(I_GT))"),
    heatmap(reconS2cgeps_check,title= "Reconstruction w S12"),
    heatmap(reconS2cgeps_check - img, title="Residual");
    layout=4,
)
plot!(size=(1000, 800))
savefig(pl12, "scratch_NM/TestPlots/s2_100itns_logIgd_smoothgk0-8.png")
#Without smoothing
mean(abs.(exp.(S12recon) - oriimg)) #2.54
mean(abs.((exp.(S12recon) - oriimg)./oriimg))

###Old DEBUG SECTION: Isolated funcs for debugging
#filename = normpath("scratch_NM", "test.jmd")
#weave(filename, out_path = :pwd)
#=
Nf = length(fhash["filt_index"])
Nx=16
storage_grad = zeros(Float64, 2+Nf+(Nf^2), Nx, Nx)

function dloss(storage_grad, img_curr, filter_hash, s_targ_mean, s_targ_invcov)
    s_curr = Deriv_Utils.DHC(img_curr, filter_hash, doS2=true, doS12=false, doS20=false)
    diff = s_curr - s_targ_mean
    diff = diff[3:end]
    term1 = Transpose(s_targ_invcov * diff)  #1x (Nf+Nf^2)
    storage_beta = Array{Float64, 2}(undef, size(term1)[1], Nx^2)
    dS1S2 = Deriv_Utils.wst_S1S2_derivfast(img_curr, filter_hash)
    term2 = reshape(dS1S2, (size(dS1S2)[1], Nx^2)) #(Nf+Nf^2) x Nx x Nx
    storage_grad = reshape(mul!(storage_beta, term1, term2), (size(term1)[1], Nx, Nx))
    #storage_grad .=  term1 * term2
end

function loss_func(img_curr, filter_hash, s_targ_mean, s_targ_invcov)
    s_curr = Deriv_Utils.DHC(img_curr,  filter_hash, doS2=true, doS12=false, doS20=false)
    neglogloss = 0.5 .* (s_curr - s_targ_mean)' * s_targ_invcov * (s_curr - s_targ_mean)
    return neglogloss
end

loss_func(dust, fhash, s_targ, I)  #Tested: with input=0, with white noise default=10^7
dloss(storage_grad, noisyinp, fhash, s_targ, I)
=#
#=
#Checking Loss Func only
s_targ = Deriv_Utils.DHC(dust, fhash, doS2=true)
s_targ_mean = s_targ[3:end]
s_targ_invcov = I
function loss_func(img_curr, filter_hash, s_targ_mean, s_targ_invcov)
    s_curr = Deriv_Utils.DHC(reshape(img_curr, (Nx, Nx)),  filter_hash, doS2=true, doS12=false, doS20=false)[3:end]
    neglogloss = 0.5 .* (s_curr - s_targ_mean)' * s_targ_invcov * (s_curr - s_targ_mean)
    println("In Loss:",size(s_curr))
    #println("NegLogLoss",neglogloss[1])
    return neglogloss[1]
end

loss_func(reshape(dust, (1, Nx^2)), fhash, s_targ_mean, s_targ_invcov)
=#

max_stars_repo_name: andrew-saydjari/DHC

max_issues_repo_name: andrew-saydjari/DHC

max_forks_repo_name: andrew-saydjari/DHC


text:
using Statistics
using Serialization
using Printf

gctime(stat) = stat.total_time

const BENCH = isempty(ARGS) ? "benches/append/append.jl" : ARGS[1]
const RUNS = isempty(ARGS) ? 10 : parse(Int, ARGS[2])
const JULIAVER = Base.julia_cmd()[1]

path = split(BENCH, "/")
file = pop!(path)
path = join(path, "/")

dir = joinpath(@__DIR__, path)
cd(dir)
@printf("Running Benchmark path = %s file = %s times = %d\n", path, file, RUNS)
@show BENCH
result = open(deserialize, `$JULIAVER --project=. $file $RUNS SERIALIZE`)
(value,times,stats)= result
@printf("run time: %0.0fms min, %0.0fms max %0.0fms median\n",
   minimum(result.times) / 1_000_000,
   maximum(result.times) / 1_000_000,
   median(result.times)  / 1_000_000)

time = map(gctime, stats)
@printf("gc time: %0.0fms min, %0.0fms max, %0.0fms median\n",
   minimum(time)/ 1_000_000,
   maximum(time)/ 1_000_000,
   median(time) / 1_000_000)

max_stars_repo_name: d-netto/GCBenchmarks

max_issues_repo_name: d-netto/GCBenchmarks

max_forks_repo_name: d-netto/GCBenchmarks


text:
function wsplint!(n, x, w)

    if n <= 9
        for i in 1:n
          f[:] .= 0.0
          f[i] = 1.0
          w[i] = splint(n,x,f)
        end
        return
    end

    f = zeros(9)
    for i in 1:4
        f[:] .= 0.0
        f[i] = 1.0
        w[i] = splint(9, x, f)
    end

    f[:] .= 0.0
    f[5] = 1.0
    for i in 5:n-4
      @views w[i] = splint(9, x[i-4:end], f)
    end

    for i in 1:4
      f[:] .= 0.0
      f[i+5] = 1.0
      @views w[n-4+i] = splint(9, x[n-8:end], f)
    end
    return
end


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct ApplicationGatewaySslPolicy <: SwaggerModel
    disabledSslProtocols::Any # spec type: Union{ Nothing, Vector{ProtocolsEnum} } # spec name: disabledSslProtocols
    policyType::Any # spec type: Union{ Nothing, String } # spec name: policyType
    policyName::Any # spec type: Union{ Nothing, PolicyNameEnum } # spec name: policyName
    cipherSuites::Any # spec type: Union{ Nothing, Vector{CipherSuitesEnum} } # spec name: cipherSuites
    minProtocolVersion::Any # spec type: Union{ Nothing, ProtocolsEnum } # spec name: minProtocolVersion

    function ApplicationGatewaySslPolicy(;disabledSslProtocols=nothing, policyType=nothing, policyName=nothing, cipherSuites=nothing, minProtocolVersion=nothing)
        o = new()
        validate_property(ApplicationGatewaySslPolicy, Symbol("disabledSslProtocols"), disabledSslProtocols)
        setfield!(o, Symbol("disabledSslProtocols"), disabledSslProtocols)
        validate_property(ApplicationGatewaySslPolicy, Symbol("policyType"), policyType)
        setfield!(o, Symbol("policyType"), policyType)
        validate_property(ApplicationGatewaySslPolicy, Symbol("policyName"), policyName)
        setfield!(o, Symbol("policyName"), policyName)
        validate_property(ApplicationGatewaySslPolicy, Symbol("cipherSuites"), cipherSuites)
        setfield!(o, Symbol("cipherSuites"), cipherSuites)
        validate_property(ApplicationGatewaySslPolicy, Symbol("minProtocolVersion"), minProtocolVersion)
        setfield!(o, Symbol("minProtocolVersion"), minProtocolVersion)
        o
    end
end # type ApplicationGatewaySslPolicy

const _property_map_ApplicationGatewaySslPolicy = Dict{Symbol,Symbol}(Symbol("disabledSslProtocols")=>Symbol("disabledSslProtocols"), Symbol("policyType")=>Symbol("policyType"), Symbol("policyName")=>Symbol("policyName"), Symbol("cipherSuites")=>Symbol("cipherSuites"), Symbol("minProtocolVersion")=>Symbol("minProtocolVersion"))
const _property_types_ApplicationGatewaySslPolicy = Dict{Symbol,String}(Symbol("disabledSslProtocols")=>"Vector{ProtocolsEnum}", Symbol("policyType")=>"String", Symbol("policyName")=>"PolicyNameEnum", Symbol("cipherSuites")=>"Vector{CipherSuitesEnum}", Symbol("minProtocolVersion")=>"ProtocolsEnum")
Base.propertynames(::Type{ ApplicationGatewaySslPolicy }) = collect(keys(_property_map_ApplicationGatewaySslPolicy))
Swagger.property_type(::Type{ ApplicationGatewaySslPolicy }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_ApplicationGatewaySslPolicy[name]))}
Swagger.field_name(::Type{ ApplicationGatewaySslPolicy }, property_name::Symbol) =  _property_map_ApplicationGatewaySslPolicy[property_name]

const _allowed_ApplicationGatewaySslPolicy_policyType = ["Predefined", "Custom"]

function check_required(o::ApplicationGatewaySslPolicy)
    true
end

function validate_property(::Type{ ApplicationGatewaySslPolicy }, name::Symbol, val)
    if name === Symbol("policyType")
        Swagger.validate_param(name, "ApplicationGatewaySslPolicy", :enum, val, _allowed_ApplicationGatewaySslPolicy_policyType)
    end
end


text:
function fastcap2{T<:Geometry}(geoms::Vector{T};kwargs...)
    # Write list file and shape files to temporary directory
    lstfile = writefilestodir(geoms)
    # Run FastCap2
    cmd = `cscript $fastcap2_vbs -l$lstfile $[string("-",a[1],a[2]) for a in kwargs] //nologo`
    res = runandparse_ffs(cmd)
    # Delete temp directory and files
    rm(dirname(lstfile),recursive=true,force=true)
    # Return result
    return res
end

function fastercap{T<:Geometry}(geoms::Vector{T};kwargs...)
    # Write full hierarchal list file
    lstfile = writefullfile(geoms)
    # Run FasterCap
    cmd = `cscript $fastercap_vbs $lstfile $[string("-",a[1],a[2]) for a in kwargs] //nologo`
    res = runandparse_ffs(cmd)
    # Delete temp file
    rm(lstfile)
    # Return result
    return res
end

function fastcapwr{T<:Geometry}(geoms::Vector{T};kwargs...)
    # Write full hierarchal list file
    lstfile = writefullfile(geoms)
    # Run fastcapwr
    cmd = `fastcap $[string("-",a[1],a[2]) for a in kwargs] -l$lstfile`
    res = runandparse_wr(cmd)
    # Delete temp file
    rm(lstfile)
    # Return result
    return res
end

function linedrawing{T<:Geometry}(geoms::Vector{T},filename::String;kwargs...)
    # Write full hierarchal list file
    lstfile = writefullfile(geoms)
    # Run fastcapwr in line drawing mode
    cmd = `fastcap -m $[string("-",a[1],a[2]) for a in kwargs] -l$lstfile`
    run(cmd)
    # Run ps2pdf
    psfile = lstfile[1:end-3]*"ps"
    run(`$ps2pdf -dEPSCrop $psfile $filename`)
    # Delete temp files
    rm(lstfile)
    rm(psfile)
end

function chargedensity{T<:Geometry}(geoms::Vector{T},filename;kwargs...)
    # Write full hierarchal list file
    lstfile = writefullfile(geoms)
    # Run fastcapwr in charge density mode
    cmd = `fastcap -q $[string("-",a[1],a[2]) for a in kwargs] -l$lstfile`
    run(cmd)
    # Run ps2pdf
    psfile = lstfile[1:end-3]*"ps"
    run(`$ps2pdf -dEPSCrop $psfile $filename`)
    # Delete temporary files
    rm(lstfile)
    rm(psfile)
end

max_stars_repo_name: jebej/FastCap

max_issues_repo_name: jebej/FastCap

max_forks_repo_name: jebej/FastCap


text:
using Test
using Mongoc,MongocUtils

struct Tvalue
   id::String
end

struct TestM
    id::String
    i::Int64
    a::Array{Dict,1}
    vt::Vector
    t::Tvalue
end

t1=TestM("aa", 1, [Dict("aa"=>1),Dict("c"=>"tt")],[Tvalue("tt")], Tvalue("tt"))
bt=Mongoc.BSON(t1)

@test bt["id"]=="aa"
@test bt["i"]==1
@test bt["a"]==[Dict("aa" => 1), Dict("c" => "tt")]
@test bt["vt"][1]["id"]=="tt"
@test bt["t"]["id"]=="tt"

## Second run after method definition
t1=TestM("aa", 1, [Dict("aa"=>1),Dict("c"=>"tt")],[Tvalue("tt")], Tvalue("tt"))
bt=Mongoc.BSON(t1)

@test bt["id"]=="aa"
@test bt["i"]==1
@test bt["a"]==[Dict("aa" => 1), Dict("c" => "tt")]
@test bt["vt"][1]["id"]=="tt"
@test bt["t"]["id"]=="tt"

module md

abstract type abs end
@enum Color red blue green

struct nest_st
    s::String
    nb::Number
    u::Union{String,Number}
    v::Vector
    z
end

struct st<:abs
    s::String
    nb::Number
    u::Union{String,Number}
    e::Color
    v::Vector
    ns::nest_st
    z
end

struct std<:abs
    d1::Dict
    d2::Dict
    d3::Dict{String,Int64}
    d4::Dict{Int64,Int64}
    d5::Dict{Int64,nest_st}
    d6
end

end
dict_complex=Dict("a"=>"a",1=>1,2=>md.nest_st("T",100.10,200,[23,"aa"],"Any"))

inst=md.st("T",123,123.123,md.red,[md.nest_st("T",100.10,200,[],"Any"),Dict("a"=>100,"b"=>"test")],md.nest_st("T",100.10,200,[23,"aa"],"Any"),Dict("a"=>100))
bson=Mongoc.BSON(inst)
inst_st=as_struct(md.st,bson)
@test inst_st.s=="T"
@test inst_st.u==123.123
@test inst_st.v[1].s=="T"
@test inst_st.v[2]==Dict("a"=>100,"b"=>"test")
@test inst_st.ns.s=="T"
@test inst_st.ns.v==[23,"aa"]
@test inst_st.z==Dict("a"=>100)
@test inst_st.e==md.red

## Construction with Abstract Type
inst_st=as_struct(md.abs,bson)
@test inst_st.ns.s=="T"
@test inst_st.u==123.123
@test inst_st.v[1].s=="T"
@test inst_st.v[2]==Dict("a"=>100,"b"=>"test")
@test inst_st.ns.s=="T"
@test inst_st.ns.v==[23,"aa"]
@test inst_st.z==Dict("a"=>100)
@test inst_st.e==md.red

## Dicts
inst=md.std(Dict("ff"=>"ff"),Dict(1=>md.nest_st("T",100.10,200,[],"Any")),Dict("a"=>1),Dict(1=>1),Dict(1=>md.nest_st("T",100.10,200,[],"Any")),dict_complex)
bson=Mongoc.BSON(inst)

inst_st=as_struct(md.std,bson)
@test inst_st.d1==Dict("ff"=>"ff")
@test inst_st.d2[1] isa md.nest_st
@test inst_st.d3["a"]==1
@test inst_st.d4[1]==1
@test inst_st.d5[1] isa md.nest_st
@test inst_st.d6[1]==1
@test inst_st.d6[2] isa md.nest_st

inst_st=as_struct(md.abs,bson)
@test inst_st.d1==Dict("ff"=>"ff")
@test inst_st.d2[1] isa md.nest_st
@test inst_st.d3["a"]==1
@test inst_st.d4[1]==1
@test inst_st.d5[1] isa md.nest_st
@test inst_st.d6[1]==1
@test inst_st.d6[2] isa md.nest_st

## Complex, Dict with non string keys and datatypes
mutable struct complexStruct{T}
    a::T
    b::String
    d::Dict
    v::Vector{Int64}
    t::Type
    ta
end

p1=complexStruct{Int64}(1,"t",Dict("1"=>md.nest_st("T",100.10,200,[],"Any")),[1,2,3],Float64,Number)

bs=Mongoc.BSON(p1)
inst_st=as_struct(complexStruct{Int64},bs)
@test inst_st.t==Float64
@test inst_st.d["1"] isa md.nest_st
@test inst_st.v==[1,2,3]
@test inst_st.ta==Number

p1=complexStruct{Int64}(1,"t",Dict(1=>md.nest_st("T",100.10,200,[],"Any")),[1,2,3],md.st,md.st)
bs=Mongoc.BSON(p1)
inst_st=as_struct(complexStruct{Int64},bs)
@test inst_st.t==md.st
@test inst_st.v==[1,2,3]
@test inst_st.d[1] isa md.nest_st
@test inst_st.ta==md.st

## Dict
bson=Mongoc.BSON(Dict("a"=>md.nest_st("T",100.10,200,[],"Any"),"b"=>100,1=>1,2=>md.nest_st("T",100.10,200,[],"Any"),3=>Dict(2=>1),4=>Dict(1=>md.nest_st("T",100.10,200,[],"Any")),:a=>1,5=>:b))
d=as_struct(Dict,bson)
@test d[2] isa Main.md.nest_st
@test d["a"] isa Main.md.nest_st
@test d[3][2]==1
@test d[1]==1
@test d[4][1] isa Main.md.nest_st
@test d[:a]==1
@test d[5]==:b

## deep dict
bson=Mongoc.BSON(Dict(1=>Tvalue("T"),"a"=>1,2=>Dict(1=>1),:a=>"11",:b=>Dict(:b=>Dict(1=>Tvalue("T"))),Tvalue("T")=>"a",6=>[Dict(1=>1),12,"a",Tvalue("T")],
        "6"=>[Dict(1=>1),12,"a",Tvalue("T")]))
d=as_struct(Dict,bson)
@test d[1].id=="T"
@test d["a"]==1
@test d[2][1]==1
@test d[:a]=="11"
@test d[:b][:b][1]==Tvalue("T")
@test d[Tvalue("T")]=="a"
@test d[6][1][1]==1
@test d[6][2]==12
@test d[6][3]=="a"
@test d[6][4]==Tvalue("T")
@test d["6"][1][1]==1
@test d["6"][2]==12
@test d["6"][3]=="a"
@test d["6"][4]==Tvalue("T")

## Symbols and nested arrays
struct TestS
    a::Vector{Vector{Any}}
    s::Symbol
end

t1=TestS([[1,2,3.4,"aaa",:s]],:A)
bson=Mongoc.BSON(t1)
d=as_struct(TestS,bson)

@test d.a[1][1]==1
@test d.a[1][3]==3.4
@test d.a[1][4]=="aaa"
@test d.a[1][5]==:s
@test d.s==:A