text:
using PandocFiltersLiveJuliaCode: walk!, Plain, Null, Code, Str, Space
using Test, JSON
import PandocFiltersLiveJuliaCode

PandocFiltersLiveJuliaCode.walk!(x, y) = walk!(x, y, "", Dict{String,Any}())
PandocFiltersLiveJuliaCode.walk!(x) = walk!(x, (t, c, m, f) -> nothing)

@testset "Manual roundtrip" begin
  manual_str = open("MANUAL.json", "r") do f
    read(f, String)
  end
  j_manual = JSON.parse(manual_str)

  doc_out = walk!(j_manual, (t, c, m, f) -> nothing)

  @test doc_out == j_manual == JSON.parse(manual_str)

end

@testset "Paragraph" begin
  para =
    JSON.parse(raw"""{"c":[{"c":"Brief","t":"Str"},{"t":"Space"},{"c":"mathematical","t":"Str"}],"t":"Para"}""")
  j_para = JSON.json(para)
  p_in = copy(para)
  j_test_no_change = JSON.json(walk!(para, (t, c, f, m) -> nothing))

  # Test the json was unchanged
  @test j_para == j_test_no_change
  out = walk!(para, (t, c, f, m) -> t == "Space" ? [] : nothing)
  j_test_no_space = JSON.json(out)

  @test j_test_no_space ==
        raw"""{"c":[{"c":"Brief","t":"Str"},{"c":"mathematical","t":"Str"}],"t":"Para"}"""

  # Check para was mutated
  @test para != p_in
  @test para == out

  # Try doubling the spaces
  j_test_double_space = JSON.json(walk!(
    p_in,
    (t, c, f, m) -> t == "Space" ? [Space(), Space()] : nothing,
  ))
  @test j_test_double_space ==
        raw"""{"c":[{"c":"Brief","t":"Str"},{"t":"Space"},{"t":"Space"},{"c":"mathematical","t":"Str"}],"t":"Para"}"""

end

@testset "Testing walk with (t,c)->1" begin
  action(t, c, format, meta) = 1
  @test walk!(1, action) == 1
  @test walk!(["string1", "string2"], action) == ["string1", "string2"]
  @test walk!([1, 2, 3], action) == [1, 2, 3]
  @test walk!(Dict("x" => 2, "y" => 3), action) == Dict("x" => 2, "y" => 3)
  @test walk!([Dict("t" => 2, "c" => 3)], action) == [1]
end

@testset "Testing walk with t=>x -> t=>y" begin
  action(t, c, format, meta) = (t == "x") ? Dict("t" => "y", "c" => c) : nothing
  dict_x = Dict("t" => "x", "c" => "z")
  dict_y = Dict("t" => "y", "c" => "z")
  dict_z = Dict("t" => "z", "c" => "w")
  @test walk!([dict_x], action, "", Dict()) == [dict_y]
  @test walk!([dict_x, dict_y, dict_z], action) == [dict_y, dict_y, dict_z]
  @test walk!(Dict("t" => "w", "c" => [dict_x]), action) ==
        Dict("t" => "w", "c" => [dict_y])
  @test walk!([dict_z, Dict("t" => "w", "c" => [dict_x])], action) ==
        [dict_z, Dict("t" => "w", "c" => [dict_y])]
end

@testset "Testing Pandoc elements" begin
  @test Plain("Plain text") == Dict("t" => "Plain", "c" => "Plain text")
  @test Null() == Dict("t" => "Null")
  @test Code(["fun"; Any[[], []]], "1+1") ==
        Dict("t" => "Code", "c" => Any[["fun"; Any[[], []]], "1+1"])
end

@testset "Testing walk with string replacement" begin
  action(t, c, format, meta) =
    (t == "Str") ? Dict("t" => "Str", "c" => "This string was $c") : nothing
  @test walk!([Str("abc"), Code("xy", "ab")], action) ==
        [Str("This string was abc"), Code("xy", "ab")]
end


text:
mat = [
  75 00 00 00 00 00 00 00 00 00 00 00 00 00 00;
  95 64 00 00 00 00 00 00 00 00 00 00 00 00 00;
  17 47 82 00 00 00 00 00 00 00 00 00 00 00 00;
  18 35 87 10 00 00 00 00 00 00 00 00 00 00 00;
  20 04 82 47 65 00 00 00 00 00 00 00 00 00 00;
  19 01 23 75 03 34 00 00 00 00 00 00 00 00 00;
  88 02 77 73 07 63 67 00 00 00 00 00 00 00 00;
  99 65 04 28 06 16 70 92 00 00 00 00 00 00 00;
  41 41 26 56 83 40 80 70 33 00 00 00 00 00 00;
  41 48 72 33 47 32 37 16 94 29 00 00 00 00 00;
  53 71 44 65 25 43 91 52 97 51 14 00 00 00 00;
  70 11 33 28 77 73 17 78 39 68 17 57 00 00 00;
  91 71 52 38 17 14 91 43 58 50 27 29 48 00 00;
  63 66 04 68 89 53 67 30 73 16 69 87 40 31 00;
  04 62 98 27 23 09 70 98 73 93 38 53 60 04 23
]

n, m = size(mat)
for i = 2:n
  for j = 1:i
    mat[i,j] += max(mat[i-1,j], j>1 ? mat[i-1,j-1] : 0)
  end
end
println(maximum(mat[end,:]))

ext: jl

max_stars_repo_name: sujimodern/project-euler-jl

max_issues_repo_name: sujimodern/project-euler-jl

max_forks_repo_name: sujimodern/project-euler-jl


text:
using Random
get_image(x::TrainingTest) = x.image
get_action_image(x::TrainingTest) = x.action_image
get_action(x::TrainingTest) = x.action
get_value(x::TrainingTest) = x.value
get_reward(x::TrainingTest) = x.reward
Base.rand(x::TrainingTests) = rand(x.tests)
Base.getindex(x::TrainingTests, idx::Integer) = x.tests[idx]
Base.getindex(x::TrainingTests, idx::Vector{<:Integer}) = x.tests[idx]
Base.length(x::TrainingTests) = length(x.tests)
Base.lastindex(x::TrainingTests) = length(x)
Base.iterate(x::TrainingTests, state=1) = state > length(x) ? nothing : (x[state], state+1)

function initialinferencefn(tests::TrainingTests)
    function fn(n::Network)
        initial_inference(n, tests.images) |> NetworkOutput
    end
end

function recurrentinferencefn(tests::TrainingTests)
    initial_inference_fn = initialinferencefn(tests)
    function fn(n::Network)
        recurrent_inference(
            n,
            hidden_state(
                initial_inference_fn(n)
            ),
            tests.action_images
        ) |> NetworkOutput
    end
end

# @todo - Chopping block
#
# function recurrentinferencefn(tests::TrainingTests, initial_output)
#     function fn(n::Network)
#         recurrent_inference(
#             n,
#             hidden_state(initial_output),
#             tests.action_images
#         ) |> NetworkOutput
#     end
# end

function getvaluemse(tests::TrainingTests, net_output)
    Flux.mse(
        value(net_output),
        tests.actual_value
    )
end

function getrewardmse(tests::TrainingTests, net_output)
    Flux.mse(
        reward(net_output),
        tests.actual_reward
    )
end

function getpolicyfails(tests::Vector{<:TrainingTest}, net_output)
    _policy = policy(net_output)
    greedy_actions = [coord[1] for coord in argmax(_policy, dims=1)[1, :]]
    fails = filter(i->begin
        test = tests[i]
        pred = greedy_actions[i]
        actuals = test.actions
        !in(pred, actuals)
    end, 1:length(tests))
    fails
end

function getpolicyfails(test::TrainingTest, net_output)
    getpolicyfails([test], net_output)
end

function getpolicyacc(tests::TrainingTests, net_output)
    fails = getpolicyfails(tests.tests, net_output)
    total = length(tests.tests)
    acc = PolicyAcc(total, length(fails))
    acc, fails
end


text:
mutable struct LoggablesXMLElement <: MyXMLElement
    els::Vector{<:MyXMLElement}
    already_made::Vector{Bool}
    function LoggablesXMLElement(els::Vector{<:MyXMLElement},
                                 already_made::Vector{Bool})
        if length(els) != length(already_made)
            error("Both arguments must be of the same length")
        end
        return new(els, already_made)
    end
end

function LoggablesXMLElement()
    return LoggablesXMLElement(MyXMLElement[], Bool[])
end

function LoggablesXMLElement(els::Vector{<:MyXMLElement};
                             already_made::Bool = true)
    return LoggablesXMLElement(els, fill(already_made, length(els)))
end

function Base.show(io::IO, lx::LoggablesXMLElement)
    disp = "$(typeof(lx)) with components:"
    for i = 1:length(lx.els)
        disp *= "\n\t$(typeof(lx.els[i])) (already_made = $(lx.already_made[i]))"
    end
    println(io, disp)
end

function make_xml(lg_el::LoggablesXMLElement)
    for el in lg_el.els
        make_xml(el)
    end
end

function join_loggables(lg1::LoggablesXMLElement, lg2::LoggablesXMLElement)
    return LoggablesXMLElement([lg1.els; lg2.els], [lg1.already_made; lg2.already_made])
end

import Base: push!
function push!(lg1::LoggablesXMLElement, lg2::LoggablesXMLElement)
    lg1.els = [lg1.els; lg2.els]
    lg1.already_made = [lg1.already_made; lg2.already_made]
end

function add_loggable(lg_el::LoggablesXMLElement, my_el::MyXMLElement;
                    already_made::Bool = false)
    push!(lg_el.els, my_el)
    push!(lg_el.already_made, already_made)
end

mutable struct MatrixInverseXMLElement <: SimpleXMLElement
    el::XMLOrNothing
    matrixParam_provider::MyXMLElement

    MatrixInverseXMLElement(my_xml::MyXMLElement) = new(nothing, my_xml)
end

function make_xml(mi_el::MatrixInverseXMLElement)
    make_xml(mi_el.matrixParam_provider)
    mp_el = get_matrix_parameter(mi_el.matrixParam_provider)

    el = new_element(bn.MATRIX_INVERSE)
    id = attribute(mp_el, bn.ID)
    set_attribute(el, bn.ID, "inverse.$id")
    add_ref_el(el, mp_el)

    mi_el.el = el
    return el
end

function get_loggables(mi_el::MatrixInverseXMLElement)
    make_xml(mi_el)
    return [mi_el.el]
end

mutable struct CorrelationMatrixXMLElement <: SimpleXMLElement
    el::XMLOrNothing
    matrixParam_provider::MyXMLElement
    invert::Bool

    CorrelationMatrixXMLElement(my_xml::MyXMLElement) = new(nothing, my_xml, false)
    CorrelationMatrixXMLElement(my_xml::MyXMLElement, invert::Bool) = new(nothing, my_xml, invert)
end

function make_xml(cm_el::CorrelationMatrixXMLElement)
    make_xml(cm_el.matrixParam_provider)
    mp_el = get_matrix_parameter(cm_el.matrixParam_provider)

    el = new_element(bn.CORRELATION_MATRIX)
    id = attribute(mp_el, bn.ID)
    nm = "correlation"
    if cm_el.invert
        nm = "$nm.inverted"
    end
    set_attribute(el, bn.ID, "$nm.$id")
    set_attribute(el, bn.INVERT, cm_el.invert)
    add_ref_el(el, mp_el)

    cm_el.el = el
    return el
end

function get_loggables(cm_el::CorrelationMatrixXMLElement)
    make_xml(cm_el)
    return [cm_el.el]
end

mutable struct VarianceProportionXMLElement <: SimpleXMLElement
    el::XMLOrNothing
    td_el::TraitLikelihoodXMLElement
    tm_el::TreeModelXMLElement
    rm_el::RepeatedMeasuresXMLElement
    mbd_el::MBDXMLElement
    id::String
    emp_var::Bool

    function VarianceProportionXMLElement(td::TraitLikelihoodXMLElement,
                                tm::TreeModelXMLElement,
                                rm::RepeatedMeasuresXMLElement,
                                mbd::MBDXMLElement)

        return new(nothing, td, tm, rm, mbd, bn.VARIANCE_PROPORTION_STATISTIC,
                    false)
    end

end

function make_xml(vp_el::VarianceProportionXMLElement)
    make_xml(vp_el.td_el)
    make_xml(vp_el.tm_el)
    make_xml(vp_el.rm_el)
    make_xml(vp_el.mbd_el)

    el = new_element(bn.VARIANCE_PROPORTION_STATISTIC)
    set_attribute(el, bn.ID, vp_el.id)
    set_attribute(el, bn.MATRIX_RATIO, bn.COHERITABILITY)
    set_attribute(el, bn.USE_EMPIRICAL_VARIANCE, vp_el.emp_var)

    add_ref_el(el, vp_el.td_el.el)
    add_ref_el(el, vp_el.tm_el.el)
    add_ref_el(el, vp_el.rm_el.el)
    add_ref_el(el, vp_el.mbd_el.el)

    vp_el.el = el
    return el
end

function get_loggables(vp_el::VarianceProportionXMLElement)
    make_xml(vp_el)
    return [vp_el.el]
end

mutable struct ModelExtensionLoggerXMLElement <: MyXMLElement
    el::XMLOrNothing
    tl_el::TraitLikelihoodXMLElement
    ext_el::MyXMLElement

    ModelExtensionLoggerXMLElement(tl::TraitLikelihoodXMLElement, my::MyXMLElement) =
        new(nothing, tl, my)
end

function make_xml(mel_el::ModelExtensionLoggerXMLElement)
    el = new_element(bn.MODEL_EXTENSION_TRAIT_LOGGER)
    set_attributes(el, [(bn.ID, bn.MODEL_EXTENSION_TRAIT_LOGGER),
                        (bn.DIMENSIONS, bn.MISSING),
                        (bn.TRAIT_NAME, mel_el.tl_el.attrs[bn.TRAIT_NAME])])

    make_xml(mel_el.tl_el)
    make_xml(mel_el.ext_el)

    add_ref_el(el, mel_el.tl_el.el)
    add_ref_el(el, mel_el.ext_el.el)

    mel_el.el = el
    return el
end

function get_loggables(mel_el::ModelExtensionLoggerXMLElement)
    return [mel_el.el]
end

# mutable struct MatrixValidationXMLElement <: MyXMLElement
#     el::XMLOrNothing
#     matrixParam_provider::MyXMLElement
#     true_parameter::AbstractArray{Float64}
#
#     MatrixValidationXMLElement(my_xml::MyXMLElement, x::AbstractArray{Float64, 2})
#         = new(my_xml, x)
# end
#
# function make_xml(mv_el::MatrixValidationXMLElement)
#     make_xml(mv_el.matrixParam_provider)
#
#     el = new_element(bn.MATRIX_VALIDATION)
#     tp_el = new_child(el, bn.TRUE_PARAMETER)
#     add_matrix_parameter(tp_el, true_parameter)
#     ip_el = new_child(el, bn.INFERRED_PARAMETER)
#     mat_param = get_matrix_parameter(mv_el.matrixParam_provider)
#     add_ref_el(ip_el, mat_param)
#


text:
# Do everything inside a let block to not introduce new variables in the REPL
let trixi_root_dir
  # Get Trixi root directory
  trixi_root_dir = dirname(@__DIR__)

  # Set project to current directory
  import Pkg
  Pkg.activate(trixi_root_dir)

  # Try to load Revise
  revise_message = "Revise initialized: changes to Trixi source code are tracked."
  try
    @eval using Revise
  catch
    # Do nothing... it probably means that Revise is not installed
    revise_message = "Revise not found (run `julia -e 'import Pkg; Pkg.add(\"Revise\")'` to install)."
  end

  # Load Trixi
  using Trixi

  # Inform user
  @info """$revise_message
           Project directory set to '$trixi_root_dir'. Adding/removing packages will only affect this project.

           Execute the following line to start a Trixi simulation:

           Trixi.run(default_example())
           """
end


text:
# TODO: higher level abstraction once it's as fast
for T in (:GeometricJacobian, :MomentumMatrix)
    @eval @inline function set_col!(dest::$T, col::Integer, src::$T)
        @framecheck dest.frame src.frame
        @boundscheck size(src, 2) == 1 || throw(ArgumentError())
        @boundscheck col ∈ Base.OneTo(size(dest, 2)) || throw(DimensionMismatch())
        @inbounds begin
            start = LinearIndices(dest.angular)[1, col]
            dest.angular[start] = src.angular[1]
            dest.angular[start + 1] = src.angular[2]
            dest.angular[start + 2] = src.angular[3]
            dest.linear[start] = src.linear[1]
            dest.linear[start + 1] = src.linear[2]
            dest.linear[start + 2] = src.linear[3]
        end
    end
end

## findunique
function findunique(f, A::AbstractArray)
    i = findfirst(f, A)
    i === nothing && error("No results found.")
    findnext(f, A, i + 1) === nothing || error("Multiple results found.")
    @inbounds return A[i]
end

# Cached download
const module_tempdir = joinpath(Base.tempdir(), string(nameof(@__MODULE__)))

function cached_download(url::String, local_file_name::String, cache_dir::String = joinpath(module_tempdir, string(hash(url))))
    if !ispath(cache_dir)
        mkpath(cache_dir)
    end
    full_cache_path = joinpath(cache_dir, local_file_name)
    if !isfile(full_cache_path)
        download(url, full_cache_path)
    end
    full_cache_path
end

## VectorSegment: type of a view of a vector
const VectorSegment{T} = SubArray{T,1,Array{T, 1},Tuple{UnitRange{Int64}},true} # TODO: a bit too specific

quatnorm(quat::UnitQuaternion) = sqrt(quat.w^2 + quat.x^2 + quat.y^2 + quat.z^2)

## Modification count stuff
function modcount end

struct ModificationCountMismatch <: Exception
    msg::String
end
Base.showerror(io::IO, ex::ModificationCountMismatch) = print(io, "ModificationCountMismatch: $(ex.msg)")

macro modcountcheck(a, b)
    quote
        modcount($(esc(a))) == modcount($(esc(b))) || modcount_check_fail($(QuoteNode(a)), $(QuoteNode(b)), $(esc(a)), $(esc(b)))
    end
end

@noinline function modcount_check_fail(asym, bsym, a, b)
    amod = modcount(a)
    bmod = modcount(b)
    msg = "Modification count of '$(string(asym))' ($amod) does not match modification count of '$(string(bsym))' ($bmod)."
    throw(ModificationCountMismatch(msg))
end

# Bounds
"""
$(TYPEDEF)

Bounds is a scalar-like type representing a closed interval from ``lower`` to
``upper``. To indicate that a vector of values falls with some range, use a
``Vector{Bounds{T}}``.
"""
struct Bounds{T}
    lower::T
    upper::T

    function Bounds{T}(lower, upper) where T
        @assert lower <= upper
        new{T}(lower, upper)
    end
    Bounds{T}() where {T} = new{T}(typemin(T), typemax(T))
end

Bounds(lower::T1, upper::T2) where {T1, T2} = Bounds{promote_type(T1, T2)}(lower, upper)

upper(b::Bounds) = b.upper
lower(b::Bounds) = b.lower
Base.:(==)(b1::Bounds, b2::Bounds) = b1.lower == b2.lower && b1.upper == b2.upper
Base.:-(b::Bounds) = Bounds(-b.upper, -b.lower)
Base.show(io::IO, b::Bounds) = print(io, "(", lower(b), ", ", upper(b), ")")
Base.convert(::Type{Bounds{T1}}, b::Bounds{T2}) where {T1, T2} = Bounds{T1}(convert(T1, lower(b)), convert(T1, upper(b)))
Base.broadcastable(b::Bounds) = Ref(b)

"""
$(SIGNATURES)

Return the closest value to ``x`` within the interval described by ``b``.
"""
Base.clamp(x, b::Bounds) = clamp(x, b.lower, b.upper)
Base.intersect(b1::Bounds, b2::Bounds) = Bounds(max(b1.lower, b2.lower), min(b1.upper, b2.upper))

# Create a new Int-backed index type, along with necessary methods to support creating ranges
macro indextype(ID)
    esc(quote
        struct $ID <: Integer
            value::Int
        end
        $ID(id::$ID) = id
        Base.hash(i::$ID, h::UInt) = hash(i.value, h)
        Base.convert(::Type{Int}, i::$ID) = i.value
        Base.Integer(i::$ID) = i.value
        Base.Int(i::$ID) = i.value
        Base.convert(::Type{$ID}, i::Integer) = $ID(i)
        Base.promote_type(::Type{Int}, ::Type{$ID}) = $ID
        Base.promote_type(::Type{$ID}, ::Type{Int}) = $ID
        Base.:<(x::$ID, y::$ID) = x.value < y.value
        Base.:<=(x::$ID, y::$ID) = x.value <= y.value
        Base.:-(x::$ID, y::$ID) = x.value - y.value
        Base.:+(x::$ID, y::Int) = $ID(x.value + y)
    end)
end


text:
# Standard LP interface
importall MathProgBase.SolverInterface

###############################################################################
# Solver objects
export IpoptSolver
immutable IpoptSolver <: AbstractMathProgSolver
    options
end
IpoptSolver(;kwargs...) = IpoptSolver(kwargs)

type IpoptMathProgModel <: AbstractNonlinearModel
    inner::IpoptProblem
    numvar::Int
    numconstr::Int
    warmstart::Vector{Float64}
    options

    function IpoptMathProgModel(;options...)
        model = new()
        model.options   = options
        model.numvar    = 0
        model.numconstr = 0
        model.warmstart = Float64[]
        model
    end
end
NonlinearModel(s::IpoptSolver) = IpoptMathProgModel(;s.options...)
LinearQuadraticModel(s::IpoptSolver) = NonlinearToLPQPBridge(NonlinearModel(s))

###############################################################################
# Begin interface implementation

# generic nonlinear interface
function loadproblem!(m::IpoptMathProgModel, numVar::Integer, numConstr::Integer, x_l, x_u, g_lb, g_ub, sense::Symbol, d::AbstractNLPEvaluator)

    m.numvar = numVar
    features = features_available(d)
    has_hessian = (:Hess in features)
    if has_hessian
        initialize(d, [:Grad, :Jac, :Hess])
        Ihess, Jhess = hesslag_structure(d)
    else
        initialize(d, [:Grad, :Jac])
        Ihess = Int[]
        Jhess = Int[]
    end
    Ijac, Jjac = jac_structure(d)
    @assert length(Ijac) == length(Jjac)
    @assert length(Ihess) == length(Jhess)
    @assert sense == :Min || sense == :Max

    # Objective callback
    scl = (sense == :Min) ? 1.0 : -1.0
    eval_f_cb(x) = scl*eval_f(d,x)

    # Objective gradient callback
    eval_grad_f_cb(x, grad_f) = (eval_grad_f(d, grad_f, x); scale!(grad_f,scl))

    # Constraint value callback
    eval_g_cb(x, g) = eval_g(d, g, x)

    # Jacobian callback
    function eval_jac_g_cb(x, mode, rows, cols, values)
        if mode == :Structure
            for i in 1:length(Ijac)
                rows[i] = Ijac[i]
                cols[i] = Jjac[i]
            end
        else
            eval_jac_g(d, values, x)
        end
    end

    if has_hessian
        # Hessian callback
        function eval_h_cb(x, mode, rows, cols, obj_factor,
            lambda, values)
            if mode == :Structure
                for i in 1:length(Ihess)
                    rows[i] = Ihess[i]
                    cols[i] = Jhess[i]
                end
            else
                obj_factor *= scl
                eval_hesslag(d, values, x, obj_factor, lambda)
            end
        end
    else
        eval_h_cb = nothing
    end

    m.inner = createProblem(numVar, float(x_l), float(x_u), numConstr,
    float(g_lb), float(g_ub), length(Ijac), length(Ihess),
    eval_f_cb, eval_g_cb, eval_grad_f_cb, eval_jac_g_cb,
    eval_h_cb)
    m.inner.sense = sense
    if !has_hessian
        addOption(m.inner, "hessian_approximation", "limited-memory")
    end

end

getsense(m::IpoptMathProgModel) = m.inner.sense
numvar(m::IpoptMathProgModel) = m.numvar
numconstr(m::IpoptMathProgModel) = m.numconstr

numlinconstr(m::IpoptMathProgModel) = 0

numquadconstr(m::IpoptMathProgModel) = 0

function optimize!(m::IpoptMathProgModel)
    copy!(m.inner.x, m.warmstart) # set warmstart
    for (name,value) in m.options
        sname = string(name)
        if match(r"(^resto_)", sname) != nothing
            sname = replace(sname, r"(^resto_)", "resto.")
        end
        addOption(m.inner, sname, value)
    end
    solveProblem(m.inner)
end

function status(m::IpoptMathProgModel)
    # Map all the possible return codes, as enumerated in
    # Ipopt.ApplicationReturnStatus, to the MPB statuses:
    # :Optimal, :Infeasible, :Unbounded, :UserLimit, and :Error
    stat_sym = ApplicationReturnStatus[m.inner.status]
    if  stat_sym == :Solve_Succeeded ||
        stat_sym == :Solved_To_Acceptable_Level
        return :Optimal
    elseif stat_sym == :Infeasible_Problem_Detected
        return :Infeasible
    elseif stat_sym == :Diverging_Iterates
        return :Unbounded
        # Things that are more likely to be fixable by changing
        # a parameter will be treated as UserLimit, although
        # some are error-like too.
    elseif stat_sym == :User_Requested_Stop ||
        stat_sym == :Maximum_Iterations_Exceeded ||
        stat_sym == :Maximum_CpuTime_Exceeded
        return :UserLimit
    else
        # Default is to not mislead user that it worked
        # Includes:
        #   :Search_Direction_Becomes_Too_Small
        #   :Feasible_Point_Found
        #   :Restoration_Failed
        #   :Error_In_Step_Computation
        #   :Not_Enough_Degrees_Of_Freedom
        #   :Invalid_Problem_Definition
        #   :Invalid_Option
        #   :Invalid_Number_Detected
        #   :Unrecoverable_Exception
        #   :NonIpopt_Exception_Thrown
        #   :Insufficient_Memory
        #   :Internal_Error
        warn("Ipopt finished with status $stat_sym")
        return :Error
    end

end
getobjval(m::IpoptMathProgModel) = m.inner.obj_val * (m.inner.sense == :Max ? -1 : +1)
getsolution(m::IpoptMathProgModel) = m.inner.x

function getreducedcosts(m::IpoptMathProgModel)
    sense = m.inner.sense
    redcost = m.inner.mult_x_U - m.inner.mult_x_L
    return sense == :Max ? redcost : -redcost
end
function getconstrduals(m::IpoptMathProgModel)
    v = m.inner.mult_g # return multipliers for all constraints
    return m.inner.sense == :Max ? copy(v) : -v
end

getrawsolver(m::IpoptMathProgModel) = m.inner
setwarmstart!(m::IpoptMathProgModel, x) = (m.warmstart = x)


text:
module HypoTest

import Distributions
import PyPlot
import Statistics

include("histdata.jl")
include("momentestim.jl")
include("disttest.jl")
include("othertest.jl")

struct TestResult
	H0::String
	H1::String
	estim::Float64
	critical::Float64
	alpha::Float64
end

export meanconfid, stdconfid
export histfit
export testhist
export corrp

end #module

max_stars_repo_name: emenems/HypoTest

max_issues_repo_name: emenems/HypoTest

max_forks_repo_name: emenems/HypoTest


text:
include("ConstantSpace.jl")

include("IntervalSpace.jl")
include("PolynomialSpace.jl")
include("PeriodicSpace.jl")

include("Modifier/Modifier.jl")

## Union

# union_rule dictates how to create a space that both spaces can be converted to
# in this case, it means
function union_rule(s1::PiecewiseSpace{S1},s2::PiecewiseSpace{S2}) where {S1<:Tuple{Vararg{PolynomialSpace}},
                    S2<:Tuple{Vararg{PolynomialSpace}}}
    PiecewiseSpace(map(Space,merge(domain(s1),domain(s2)).domains))
end

function union_rule(s1::PiecewiseSpace{S1},s2::PolynomialSpace) where S1<:Tuple{Vararg{PolynomialSpace}}
    PiecewiseSpace(map(Space,merge(domain(s1),domain(s2)).domains))
end

include("Chebyshev/Chebyshev.jl")
include("Ultraspherical/Ultraspherical.jl")
include("Fourier/Fourier.jl")
include("Singularities/Singularities.jl")
include("Jacobi/Jacobi.jl")
include("Hermite/Hermite.jl")
include("Laguerre/Laguerre.jl")
include("CurveSpace.jl")

## Derivative

function invfromcanonicalD(S::Laurent{PeriodicLine{false}})
    d=domain(S)
    @assert d.center==0  && d.L==1.0
    a=Fun(Laurent(),[1.,.5,.5])
end

function invfromcanonicalD(S::LaurentDirichlet{PeriodicLine{false}})
    d=domain(S)
    @assert d.center==0  && d.L==1.0
    a=Fun(Laurent(),[1.,.5,.5])
end


text:
# Helper functions for generating train/test data

using Random

function batch_dataset((xs, ys), batch_size)
    if length(xs) % batch_size != 0
        error("batch_size must evenly divide dataset")
    end
    batch_xs = [hcat(xs[(i-1)*batch_size+1:i*batch_size]...) for i = 1:div(length(xs), batch_size)]
    batch_ys = [hcat(ys[(i-1)*batch_size+1:i*batch_size]...) for i = 1:div(length(xs), batch_size)]
    return (batch_xs, batch_ys)
end

"""
Generate random vectors with constrained magnitudes.

If a_inds (and optionally b_inds) is provided, then
the magnitude constraints will apply to the subset
of indices specified by these vectors.

Returns a single vector of values.
"""
function randmgn(op, input_size, min, max, a_inds, b_inds=nothing)
    if max < min
        error("min should be less than max in randmgn!")
    end
    if op(min, min) >= max
        error("$op(min, min) results in a value greater than max! min: $min, max: $max")
    end
    while true
        sample = rand(input_size) .* 2*max .- max
        a = sum(sample[a_inds])
        aabs = abs(a)
        if b_inds !== nothing
            b = sum(sample[b_inds])
            babs = abs(b)
            yabs = abs(op(a, b))
        else
            yabs = abs(op(a))
        end
        if aabs < max && aabs > min && (b_inds === nothing || (babs < max && babs > min)) && yabs < max && yabs > min
            return sample
        end
    end
end

"""
Generates a dataset for binary operations.
"""
function binary_data(bin_op,
                     num_samples::Integer,
                     input_size::Integer,
                     max::Integer;
                     perm=nothing,
                     min=0,
                    )
    if perm === nothing
        perm = randperm(input_size)
    end
    a_inds = perm[1:div(input_size, 2)]
    b_inds = perm[div(input_size, 2)+1:end]

    data_xs = [randmgn(bin_op, input_size, min, max, a_inds, b_inds) for _ in 1:num_samples]
    data_abs = [(sum(x[a_inds]), sum(x[b_inds])) for x in data_xs]
    data_ys = [bin_op(a, b) for (a, b) in data_abs]

    return (data_xs, data_ys), (perm, a_inds, b_inds, data_abs)
end

"""
Generates a dataset for unary operations.

Unlike the binary_data function, this does not generate
negative values.
"""
function unary_data(op,
                    num_samples::Integer,
                    input_size::Integer,
                    max::Integer,
                    perm=nothing,
                    min=0
                    )
    if perm === nothing
        perm = randperm(input_size)
    end
    a_inds = perm[1:div(input_size, 2)]
    b_inds = perm[div(input_size, 2)+1:end]

    data_xs = [rand(input_size) .* (2*max) for _ in 1:num_samples]
    data_abs = [(sum(x[a_inds]), sum(x[b_inds])) for x in data_xs]
    data_ys = [op(a) for (a, _) in data_abs]

    return (data_xs, data_ys), (perm, a_inds, b_inds, data_abs)
end

max_stars_repo_name: dillondaudert/neural-alu

max_issues_repo_name: dillondaudert/neural-alu

max_forks_repo_name: dillondaudert/neural-alu


text:
using Luxor, Colors

let
    cols = [Luxor.julia_red, Luxor.julia_green, Luxor.julia_purple]
    Drawing(500, 500,  "/tmp/turing.svg")
    origin()
    sethue("grey5")
    circle(O, 245, :fillpreserve)
    circle(O, 245, :clip)

    translate(-80, 100)

    k = 140 # widens the curve
    h = 130 # height of curve
    redcurve = [Point(k * x, -h * exp(-x^2)) for x in   -10:0.01:10]
    greencurve = [Point(0.65k * x, -250 * exp(-x^2)) for x in -10:0.01:10]
    purplecurve = [Point(0.5k * x, -120 * exp(-x^2)) for x in -10:0.01:10]

    # curves need shifting ↔
    polymove!(greencurve, O, Point(100, 0))
    polymove!(purplecurve, O, Point(240, 5))

    # the thickening
    g(x, θ) = rescale(abs(sin(θ) * sin(θ)), 0, 1, 8, 12)

    for (i, thecurve) in enumerate([redcurve, greencurve, purplecurve])
        @layer begin
            pgon = offsetpoly(thecurve, g)
            poly(pgon, :clip)
            for j in 0:0.1:1
                c = convert(HSB, RGB(cols[i]...))
                sethue(HSB(c.h, 0.8, rescale(j, 0, 1, 0.4, 1.0)))
                setline(20-20j)
                poly(thecurve, :stroke)
            end
        end
    end
    clipreset()
    origin()
    setline(10)
    sethue(0.78, 0.88, 0.95)
    circle(O, 240, :stroke)
    finish()
    preview()
end

max_stars_repo_name: radovankavicky/graphics

max_issues_repo_name: radovankavicky/graphics

max_forks_repo_name: cormullion/graphics


text:
BinQ2low = [650  ,   650  ,   650  ,   650  ,   650  ,   650  ,   650  ,   650  ,   800  ,   800  ,   800  ,   800  ,   800  ,   800  ,   800  ,   800  ,   800  ,   950  ,   950  ,   950  ,   950  ,   950  ,   950  ,   950  ,   950  ,   950  ,   1100  ,   1100  ,   1100  ,   1100  ,   1100  ,   1100  ,   1100  ,   1100  ,   1100  ,   1300  ,   1300  ,   1300  ,   1300  ,   1300  ,   1300  ,   1300  ,   1300  ,   1300  ,   1500  ,   1500  ,   1500  ,   1500  ,   1500  ,   1500  ,   1500  ,   1500  ,   1500  ,   1500  ,   1800  ,   1800  ,   1800  ,   1800  ,   1800  ,   1800  ,   1800  ,   1800  ,   1800  ,   1800  ,   2100  ,   2100  ,   2100  ,   2100  ,   2100  ,   2100  ,   2100  ,   2100  ,   2100  ,   2100  ,   2100  ,   2400  ,   2400  ,   2400  ,   2400  ,   2400  ,   2400  ,   2400  ,   2400  ,   2400  ,   2400  ,   2800  ,   2800  ,   2800  ,   2800  ,   2800  ,   2800  ,   2800  ,   2800  ,   2800  ,   2800  ,   3200  ,   3200  ,   3200  ,   3200  ,   3200  ,   3200  ,   3200  ,   3200  ,   3200  ,   3200  ,   3800  ,   3800  ,   3800  ,   3800  ,   3800  ,   3800  ,   3800  ,   3800  ,   3800  ,   3800  ,   3800  ,   4500  ,   4500  ,   4500  ,   4500  ,   4500  ,   4500  ,   4500  ,   4500  ,   4500  ,   4500  ,   4500  ,   4500  ,   6000  ,   6000  ,   6000  ,   6000  ,   6000  ,   6000  ,   6000  ,   6000  ,   6000  ,   6000  ,   6000  ,   6000  ,   8000  ,   8000  ,   8000  ,   8000  ,   8000  ,   8000  ,   8000  ,   8000  ,   8000  ,   8000  ,   8000  ,   11000  ,   11000]

BinQ2high = [  800  ,   800  ,   800  ,   800  ,   800  ,   800  ,   800  ,   800  ,   950  ,   950  ,   950  ,   950  ,   950  ,   950  ,   950  ,   950  ,   950  ,   1100  ,   1100  ,   1100  ,   1100  ,   1100  ,   1100  ,   1100  ,   1100  ,   1100  ,   1300  ,   1300  ,   1300  ,   1300  ,   1300  ,   1300  ,   1300  ,   1300  ,   1300  ,   1500  ,   1500  ,   1500  ,   1500  ,   1500  ,   1500  ,   1500  ,   1500  ,   1500  ,   1800  ,   1800  ,   1800  ,   1800  ,   1800  ,   1800  ,   1800  ,   1800  ,   1800  ,   1800  ,   2100  ,   2100  ,   2100  ,   2100  ,   2100  ,   2100  ,   2100  ,   2100  ,   2100  ,   2100  ,   2400  ,   2400  ,   2400  ,   2400  ,   2400  ,   2400  ,   2400  ,   2400  ,   2400  ,   2400  ,   2400  ,   2800  ,   2800  ,   2800  ,   2800  ,   2800  ,   2800  ,   2800  ,   2800  ,   2800  ,   2800  ,   3200  ,   3200  ,   3200  ,   3200  ,   3200  ,   3200  ,   3200  ,   3200  ,   3200  ,   3200  ,   3800  ,   3800  ,   3800  ,   3800  ,   3800  ,   3800  ,   3800  ,   3800  ,   3800  ,   3800  ,   4500  ,   4500  ,   4500  ,   4500  ,   4500  ,   4500  ,   4500  ,   4500  ,   4500  ,   4500  ,   4500  ,   6000  ,   6000  ,   6000  ,   6000  ,   6000  ,   6000  ,   6000  ,   6000  ,   6000  ,   6000  ,   6000  ,   6000  ,   8000  ,   8000  ,   8000  ,   8000  ,   8000  ,   8000  ,   8000  ,   8000  ,   8000  ,   8000  ,   8000  ,   8000  ,   11000  ,   11000  ,   11000  ,   11000  ,   11000  ,   11000  ,   11000  ,   11000  ,   11000  ,   11000  ,   11000  ,   20000  ,   20000 ]

Binxlow=[0.054,   0.07,   0.088,   0.11,   0.14,   0.17,   0.21,   0.26,   0.044,   0.059,   0.076,   0.096,   0.12,   0.15,   0.19,   0.23,   0.28,   0.044,   0.059,   0.076,   0.096,   0.12,   0.15,   0.18,   0.22,   0.32,   0.048,   0.063,   0.08,   0.1,   0.13,   0.16,   0.2,   0.24,   0.34,   0.054,   0.07,   0.088,   0.11,   0.14,   0.17,   0.21,   0.25,   0.36,   0.044,   0.059,   0.076,   0.096,   0.12,   0.15,   0.19,   0.23,   0.28,   0.39,   0.054,   0.07,   0.088,   0.11,   0.14,   0.17,   0.21,   0.26,   0.31,   0.43,   0.048,   0.063,   0.08,   0.1,   0.13,   0.16,   0.2,   0.24,   0.29,   0.34,   0.46,   0.076,   0.096,   0.12,   0.15,   0.18,   0.22,   0.27,   0.32,   0.38,   0.5,   0.088,   0.11,   0.14,   0.17,   0.21,   0.25,   0.3,   0.36,   0.42,   0.54,   0.1,   0.13,   0.16,   0.2,   0.24,   0.29,   0.34,   0.4,   0.46,   0.58,   0.096,   0.12,   0.15,   0.19,   0.23,   0.28,   0.33,   0.39,   0.45,   0.51,   0.63,   0.096,   0.12,   0.15,   0.18,   0.22,   0.27,   0.32,   0.38,   0.44,   0.5,   0.56,   0.69,   0.1,   0.13,   0.16,   0.2,   0.24,   0.29,   0.35,   0.41,   0.47,   0.53,   0.59,   0.73,   0.15,   0.19,   0.23,   0.28,   0.33,   0.39,   0.45,   0.51,   0.57,   0.64,   0.78,   0.256,   0.6]

 Binxhigh=[0.07,   0.088,   0.11,   0.14,   0.17,   0.21,   0.26,   1,   0.059,   0.076,   0.096,   0.12,   0.15,   0.19,   0.23,   0.28,   1,   0.059,   0.076,   0.096,   0.12,   0.15,   0.18,   0.22,   0.32,   1,   0.063,   0.08,   0.1,   0.13,   0.16,   0.2,   0.24,   0.34,   1,   0.07,   0.088,   0.11,   0.14,   0.17,   0.21,   0.25,   0.36,   1,   0.059,   0.076,   0.096,   0.12,   0.15,   0.19,   0.23,   0.28,   0.39,   1,   0.07,   0.088,   0.11,   0.14,   0.17,   0.21,   0.26,   0.31,   0.43,   1,   0.063,   0.08,   0.1,   0.13,   0.16,   0.2,   0.24,   0.29,   0.34,   0.46,   1,   0.096,   0.12,   0.15,   0.18,   0.22,   0.27,   0.32,   0.38,   0.5,   1,   0.11,   0.14,   0.17,   0.21,   0.25,   0.3,   0.36,   0.42,   0.54,   1,   0.13,   0.16,   0.2,   0.24,   0.29,   0.34,   0.4,   0.46,   0.58,   1,   0.12,   0.15,   0.19,   0.23,   0.28,   0.33,   0.39,   0.45,   0.51,   0.63,   1,   0.12,   0.15,   0.18,   0.22,   0.27,   0.32,   0.38,   0.44,   0.5,   0.56,   0.69,   1,   0.13,   0.16,   0.2,   0.24,   0.29,   0.35,   0.41,   0.47,   0.53,   0.59,   0.73,   1,   0.19,   0.23,   0.28,   0.33,   0.39,   0.45,   0.51,   0.57,   0.64,   0.78,   1,   0.6,   1]


text:
using Colors, ColorSchemes
using CairoMakie, FileIO
using CSV, DataFrames
using Statistics

imgpath = joinpath("Molenaar", "Figure 4", "Plots")
resultspath = joinpath("Molenaar", "Results")

#x = [res["S"] for res in results_test]
#a = [res["v_catef"] for res in results_test]
#b = [res["v_metef"] for res in results_test]

#fig = Figure()
#ax = Axis(fig[1,1])
#scatter1 = scatter!(ax, x, a)
#line1 =lines!(ax, x, a)
#scatter2 = scatter!(ax, x, b)
#line2 =lines!(ax, x, b)

#ax.xlabel = "Growth rate"
#ax.ylabel = "Fraction of the total substarte flux"
#fig

#FileIO.save(joinpath(imgpath, "Figure_remake.pdf"), fig)

# Remake 2

x = [res["sub_out"] for res in results_test]
a = [res["v_catef"] for res in results_test]
b = [res["v_metef"] for res in results_test]

fig = Figure()
ax = Axis(fig[1, 1])
scatter1 = scatter!(ax, x, a)
line1 = lines!(ax, x, a)
scatter2 = scatter!(ax, x, b)
line2 = lines!(ax, x, b)

ax.xlabel = "Substrate concentration"
ax.ylabel = "Fraction of the total substarte flux"

fig[1, 1] = Legend(
    fig,
    [[scatter1, line1], [scatter2, line2]],
    ["MetEf_flux", "CatEf_flux"],
    tellheight = false,
    tellwidth = false,
    halign = :right,
    valign = :center,
    labelsize = 14,
)
fig

FileIO.save(joinpath(imgpath, "Figure_remake2.pdf"), fig)


text:
# SPD-License-Identifier: MIT

function padutterance(utt, Nmax)
    D, N = size(utt)
    pad = similar(utt, D, Nmax - N)
    pad .= utt[:, end]
    hcat(utt, pad)
end

struct BatchLoader
    h5data
    alifsms
    batchsize
    uttids
end

function BatchLoader(h5data, alifsms, batchsize; shuffledata = false)
    uttids = if shuffledata
        shuffle(collect(keys(h5data)))
    else
        sort(keys(h5data), by = k -> size(h5data[k], 2))
    end
    BatchLoader(h5data, alifsms, batchsize, uttids)
end

Base.length(bl::BatchLoader) = Int(ceil(length(bl.h5data)/bl.batchsize))

function Base.iterate(bl::BatchLoader)
    iterate(bl, 1)
end

function Base.iterate(bl::BatchLoader, state)
    idx = state
    idx <= length(bl) || return nothing

    startidx = (idx-1) * bl.batchsize + 1
    endidx = min(startidx + bl.batchsize - 1, length(bl.h5data))

    utts = [read(bl.h5data[key]) for key in bl.uttids[startidx:endidx]]
    batch_alifsms = [bl.alifsms[key] for key in bl.uttids[startidx:endidx]]
    seqlengths = [size(utt, 2) for utt in utts]

    # Merge the fsms into one.
    batch_alis = union(batch_alifsms...)

    # Pad the utterances' features to get a 3D tensor.
    Nmax = maximum(seqlengths)
    batch_data = cat(padutterance.(utts, Nmax)..., dims = 3)

    (batch_data, batch_alis, seqlengths), idx+1
end

struct DumpBatchLoader
    h5data
    batchsize
    uttids
end

function DumpBatchLoader(h5data, batchsize)
    uttids = collect(keys(h5data))
    DumpBatchLoader(h5data, batchsize, uttids)
end

Base.length(bl::DumpBatchLoader) = Int(ceil(length(bl.h5data)/bl.batchsize))

function Base.iterate(bl::DumpBatchLoader)
    iterate(bl, 1)
end

function Base.iterate(bl::DumpBatchLoader, state)
    idx = state
    idx <= length(bl) || return nothing

    startidx = (idx-1) * bl.batchsize + 1
    endidx = min(startidx + bl.batchsize - 1, length(bl.h5data))

    utts = [read(bl.h5data[key]) for key in bl.uttids[startidx:endidx]]
    seqlengths = [size(utt, 2) for utt in utts]

    # Pad the utterances' features to get a 3D tensor.
    Nmax = maximum(seqlengths)
    batch_data = cat(padutterance.(utts, Nmax)..., dims = 3)

    (batch_data, seqlengths, bl.uttids[startidx:endidx]), idx+1
end

max_stars_repo_name: MartinKocour/SpeechLab

max_issues_repo_name: MartinKocour/SpeechLab

max_forks_repo_name: MartinKocour/SpeechLab


text:
using .GPUCompiler: classify_arguments, BITS_REF, BITS_VALUE

struct VECompilerParams <: AbstractCompilerParams
    device::Int32
    global_hooks::NamedTuple
    loop_vectorizer::Function
end

VECompilerJob = CompilerJob{VECompilerTarget, VECompilerParams}

GPUCompiler.runtime_module(::VECompilerJob) = VectorEngine

# filter out functions from device libs
GPUCompiler.isintrinsic(job::VECompilerJob, fn::String) =
    invoke(GPUCompiler.isintrinsic,
           Tuple{CompilerJob{VECompilerTarget}, typeof(fn)},
           job, fn) ||
    true
    #startswith(fn, "ve")

# For now, completely disable the check for isbitstype arguments
# TODO: find a better way to identify arguments that can be passed on stack
function GPUCompiler.check_invocation(job::VECompilerJob, entry::LLVM.Function)
    return
end

function GPUCompiler.process_entry!(@nospecialize(job::CompilerJob), mod::LLVM.Module,
                                    entry::LLVM.Function)
    ctx = context(mod)

    if haskey(ENV, "JULIA_VE_ARGS") && ENV["JULIA_VE_ARGS"] == "1"
        if job.source.kernel
            args = classify_arguments(job, eltype(llvmtype(entry)))
            for arg in args
                @debug "arg : $arg"
            end
        end
    end

    return entry
end


text:
function anisodiff(im,niter,kappa,lambda,option)
    diff=zeros(Float64,size(im))
    anisodiff!(im,niter,kappa,lambda,option,diff)
    diff
end

function anisodiff!(im::Array{Float64,2}, niter::Int64, kappa::Float64, lambda::Float64, option::Int64,dd=zeros(Float64,size(dd)))
#=
% Arguments:
%         im     - input image
%         niter  - number of iterations.
%         kappa  - conduction coefficient 20-100 ?
%         lambda - max value of .25 for stability
%         option - 1 Perona Malik diffusion equation No 1
%                  2 Perona Malik diffusion equation No 2
%
% Returns:
%         diff   - diffused image.
%
% kappa controls conduction as a function of gradient.  If kappa is low
% small intensity gradients are able to block conduction and hence diffusion
% across step edges.  A large value reduces the influence of intensity
% gradients on conduction.
%
% lambda controls speed of diffusion (you usually want it at a maximum of
% 0.25)
%
% Diffusion equation 1 favours high contrast edges over low contrast ones.
% Diffusion equation 2 favours wide regions over smaller ones.

% Reference:
% P. Perona and J. Malik.
% Scale-space and edge detection using ansotropic diffusion.
% IEEE Transactions on Pattern Analysis and Machine Intelligence,
% 12(7):629-639, July 1990.
%
% Peter Kovesi
% www.peterkovesi.com/matlabfns/
%
% June 2000  original version.
% March 2002 corrected diffusion eqn No 2.
=#

    (rows,cols) = size(dd);
    #rows=480
    #cols=640

    for i = 1:niter
        for j=1:length(im)
            dd[j]=im[j]
        end

        for k=1:cols
            for j=1:rows

                if j==1
                    deltaN = -1 * dd[j,k]
                else
                    @inbounds deltaN = dd[j-1,k] - dd[j,k]
                end

                if j==rows
                    deltaS = -1 * dd[j,k]
                else
                    @inbounds deltaS = dd[j+1,k] - dd[j,k]
                end

                if k==1
                    deltaW = -1 * dd[j,k]
                else
                    @inbounds deltaW = dd[j,k-1] - dd[j,k];
                end

                if k==cols
                    deltaE = -1 * dd[j,k]
                else
                    @inbounds deltaE = dd[j,k+1] - dd[j,k];
                end

                #if option == 1
                    cN = exp(-(deltaN/kappa)*(deltaN/kappa))
                    cS = exp(-(deltaS/kappa)*(deltaS/kappa))
                    cE = exp(-(deltaE/kappa)*(deltaE/kappa))
                    cW = exp(-(deltaW/kappa)*(deltaW/kappa))
                #else
                    #cN = 1/(1 + (deltaN/kappa)*(deltaN/kappa))
                    #cS = 1/(1 + (deltaS/kappa)*(deltaS/kappa))
                    #cE = 1/(1 + (deltaE/kappa)*(deltaE/kappa))
                    #cW = 1/(1 + (deltaW/kappa)*(deltaW/kappa))
                #end

                @inbounds im[j,k] = dd[j,k] + lambda*(cN*deltaN + cS*deltaS + cE*deltaE + cW*deltaW);
            end
        end
    end

    nothing
end

#The clahe method from Images is not optomized for iterations.
#Should make a new method that preallocates a temp array to be reused.
function local_contrast_enhance!(img,out_img)

    for i=1:length(img)
        out_img[i] = img[i]/255
    end

    img2 = Images.clahe(out_img,256,xblocks=30,yblocks=30,clip=15)

    for i=1:length(img2)
        out_img[i] = img2[i]*255
    end

    nothing
end

function local_contrast_enhance(img)

    img2 = img ./ 255

    img2 = Images.clahe(img2,256,xblocks=30,yblocks=30,clip=15) .* 255

    img2
end

function subtract_background(han::Tracker_Handles)

    mydiff = han.current_frame2 .- mean_image(han)
    new_diff = (mydiff - minimum(mydiff))
    new_diff = new_diff ./ maximum(new_diff)

    han.current_frame = round.(UInt8,new_diff .* 255)

    nothing
end

# https://web.stanford.edu/class/cs448f/lectures/2.1/Sharpening.pdf
#=
Filter type, 1 = median, 2 = gaussian
=#
function sharpen_image(img::AbstractArray{T,2},sigma=21,rep = 5,filter_type=1) where T
    if iseven(sigma)
        sigma = sigma + 1
    end
    kern = KernelFactors.gaussian((sigma,sigma))
    img_out = zeros(T,size(img))
    img_temp = zeros(Float64,size(img))
    img_median = zeros(UInt8,size(img))
    if filter_type == 1
        for i=1:length(img)
            img_median[i] = round(UInt8,img[i])
        end
    end

    for i=1:length(img)
        img_temp[i] = img[i]
    end

    for i=1:rep
        if filter_type == 1
            img_temp = 1.5 .* img_temp .- 0.5 .* median_filter(img_median,(sigma,sigma))
        else
            img_temp = 1.5 .* img_temp .- 0.5 .* imfilter(img_temp,kern)
        end
        adjust_contrast(img_temp,0,255)
        if filter_type == 1
            for i=1:length(img)
                img_median[i] = round(UInt8,img_temp[i])
            end
        end
    end
    if eltype(img) <: Integer
        for i=1:length(img)
            img_out[i] = round(T,img_temp[i])
        end
    else
        img_out[:] = img_temp
    end
    img_out
end

#=
Median Filter adapted from Julia code by Tejus Gupta (c) 2018
https://gist.github.com/tejus-gupta/2f29e101e18a58dcf2477991151b3a42

Based on algorithm in 'A Fast Two-Dimensional Median Filtering Algorithm' by Huang, Yang and Tang.
=#
function median_filter(img::Array{UInt8, 2}, window_size::Tuple{Int64,Int64})

    function update_median(median_val, n_pixels_below_median, hist, half_pixels)
        if n_pixels_below_median <= half_pixels
            for val in median_val:2^8
                if n_pixels_below_median + hist[val] > half_pixels
                    median_val = val
                    break
                else
                    n_pixels_below_median += hist[val]
                end
            end
        elseif n_pixels_below_median > half_pixels
            for val in median_val-1:-1:1
                n_pixels_below_median -= hist[val]
                if n_pixels_below_median <= half_pixels
                    median_val = val
                    break
                end
            end
        end
        return median_val, n_pixels_below_median
    end

    result = copy(img)

    if iseven(window_size[1]) || iseven(window_size[2])
        error("window width and height must be odd.")
    end

    half_window_size = window_size.>>1
    img_size = size(img)
    half_pixels = (window_size[1]*window_size[2])>>1

    for j in 1 + half_window_size[2]: img_size[2] - half_window_size[2]
        hist = zeros(Int, 2^8)
        for I in CartesianIndices((1:window_size[1], j-half_window_size[2]:j+half_window_size[2]))
            hist[img[I]+1] += 1
        end

        median_val, n_pixels_below_median = update_median(1, 0, hist, half_pixels)
        result[1+half_window_size[2], j] = UInt8(median_val-1)

        for i in 2 + half_window_size[1]: img_size[1] - half_window_size[1]
            for I in CartesianIndices((i-half_window_size[1]-1, j-half_window_size[2]:j+half_window_size[2]))
                hist[img[I]+1] -= 1
                if img[I] + 1 < median_val
                    n_pixels_below_median -= 1
                end
            end

            for I in CartesianIndices((i+half_window_size[1], j-half_window_size[2]:j+half_window_size[2]))
                hist[img[I]+1] += 1
                if img[I] + 1 < median_val
                    n_pixels_below_median += 1
                end
            end

            median_val, n_pixels_below_median = update_median(median_val, n_pixels_below_median, hist, half_pixels)
            result[i, j] = UInt8(median_val-1)
        end
    end
    result
end

function adjust_contrast(img::AbstractArray{T,2},min_c::Real,max_c::Real) where T

    out_max = 255
    out_min = 0

    out_min_rel = min_c / out_max
    out_max_rel = max_c / out_max
    out_range_f = (max_c - min_c) / (out_max - out_min)

    for i=1:length(img)

        inten = img[i] / 255 #normalize to float of intensity

        inten = inten - out_min_rel # remove bottom
        if inten < 0.0
            inten = 0.0
        end

        inten = inten / out_range_f

        if inten > 1.0
            inten = 1.0
        end

        if eltype(img) <: Integer
            img[i] = round(T,inten * out_max)
        else
            img[i] = inten * out_max
        end
    end

    nothing
end


text:
using Test

@testset "m-sets" begin
  m = 2
  i = MSetInstance(Float64[5, 4, 3], m)
  g = msets_greedy(i)
  d = msets_dp(i)
  l = ! is_travis && msets_lp(i, solver=Gurobi.Optimizer)

  @test i.m == m

  @test value(g) == value(d)
  @test length(g.items) <= m
  @test length(d.items) <= m

  @test g.instance == i
  @test d.instance == i

  if ! is_travis
    @test value(g) == value(l)
    @test length(l.items) <= m
    @test l.instance == i
  end
end

@testset "Budgeted m-sets" begin
  function test_solution_at(s::BudgetedMSetSolution, kv::Dict{Int, Float64})
    for (k, v) in kv
      @test value(s, k) ≈ v
    end
  end

  function test_items_at(s::BudgetedMSetSolution, kv::Dict{Int, Vector{Int}})
    for (k, v) in kv
      @test items(s, k) == v
    end
  end

  @testset "Basic" begin
    m = 2
    i = BudgetedMSetInstance(Float64[5, 4, 3], Int[1, 1, 1], m)
    d = budgeted_msets_dp(i)
    l = ! is_travis && budgeted_msets_lp_all(i, solver=Gurobi.Optimizer)

    @test i.m == m
    @test d.instance == i

    @test d.state[m, 0 + 1, 0 + 1] == 9.0
    @test d.state[m, 0 + 1, 1 + 1] == 9.0
    @test d.state[m, 0 + 1, 2 + 1] == 9.0
    @test d.state[m, 0 + 1, 3 + 1] == -Inf

    @test d.solutions[m, 0, 0] == [1, 2]
    @test d.solutions[m, 0, 1] == [1, 2]
    @test d.solutions[m, 0, 2] == [1, 2]
    @test d.solutions[m, 0, 3] == [-1]

    if ! is_travis
      for i in 0:3
        @test l.state[m, 0 + 1, i + 1] ≈ d.state[m, 0 + 1, i + 1]
        @test l.solutions[m, 0, i] == d.solutions[m, 0, i]
      end
    end

    # Accessors.
    expected_items = Dict{Int, Vector{Int}}(0 => [1, 2], 1 => [1, 2], 2 => [1, 2], 3 => [-1])
    test_items_at(d, expected_items)
    ! is_travis && test_items_at(l, expected_items)

    expected = Dict{Int, Float64}(0 => 9.0, 1 => 9.0, 2 => 9.0, 3 => -Inf)
    test_solution_at(d, expected)
    ! is_travis && test_solution_at(l, expected)
  end

  @testset "Conformity" begin
    # More advanced tests to ensure the algorithm works as expected.

    # 1
    a = 3.612916190062782
    b = 7.225832380125564
    v = Float64[a, a, b, b, b, b, b, b, b, b]
    w = Int[32, 32, 32, 32, 0, 32, 32, 32, 0, 32]
    m = 3

    i = BudgetedMSetInstance(v, w, m)
    d = budgeted_msets_dp(i)
    # TODO: stop the algorithm in this case? Don't waste too much time on this part of the table.
    l = ! is_travis && budgeted_msets_lp_select(i, [0, 96, 97, 320], solver=Gurobi.Optimizer)
    expected = Dict{Int, Float64}(0 => 3 * b, 96 => 3 * b, 96 + 1 => -Inf, 320 => -Inf) # No more solutions after 96.
    test_solution_at(d, expected)
    ! is_travis && test_solution_at(l, expected)

    # 2
    v = [7.840854066284411, 3.9204270331422055, 7.840854066284411, 3.9204270331422055, 7.840854066284411, 7.840854066284411, 7.840854066284411, 7.840854066284411, 15.681708132568822, 5.227236044189607]
    w = [16, 32, 16, 24, 16, 16, 0, 16, 0, 32]
    m = 3

    i = BudgetedMSetInstance(v, w, m)
    d = budgeted_msets_dp(i)
    l = ! is_travis && budgeted_msets_lp_select(i, [0, 4, 20, 24, 25, 32, 33, 40, 41, 48, 49, 72, 73, 96, 97, 280, 319, 320], solver=Gurobi.Optimizer)

    a = 31.363416265137644
    b = 28.74979824304284
    c = 24.829371209900632
    e = 16.988517143616225
    expected = Dict{Int, Float64}(0 => a, 4 => a, 20 => a, 24 => a, 25 => a, 32 => a,
      33 => b, 40 => b, 41 => b, 48 => b,
      49 => c, 72 => e, 73 => e,
      96 => -Inf, 97 => -Inf, 280 => -Inf, 319 => -Inf, 320 => -Inf
    )
    test_solution_at(d, expected)
    ! is_travis && test_solution_at(l, expected)
  end
end


text:
"""
```
dualarray(a::AbstractArray, b::AbstractArray)

```

returns an `Array` that has the size of input `a` and an element type consistent with the element types
of the inputs `a` and `b`. For example, suppose you want to write a function of the form

```
julia> function f(a, b)
           F = similar(a)
           F[1] = a[1] * b[1]
           F[2] = a[2] + b[2]
       end
```

If you were to automatically differentiate `f` with respect to `b`, then `F` would not have the correct element type since its type
will be the same as `a`. Rather than require the user to write

```
julia> F = if eltype(b) <: ForwardDiff.Dual
           similar(a, eltype(b))
       else
           similar(a)
       end
```

the user can use `dualvector` to write

```@meta
DocTestSetup = quote
    import ForwardDiff
    using RiskAdjustedLinearizations, ForwardDiff
end
```

```jldoctest
julia> a = rand(3)
julia> b = ones(ForwardDiff.Dual, 5)
julia> F = RiskAdjustedLinearizations.dualvector(a, b)
3-element Array{ForwardDiff.Dual,1}:
 #undef
 #undef
 #undef
```

```@meta
DocTestSetup = nothing
```

Note that the element types of `a` and `b` must be subtypes of `Real` (or else `ForwardDiff` will not work).
"""
@inline dualarray(a::AbstractArray{<: ForwardDiff.Dual}, b::AbstractArray{<: ForwardDiff.Dual}) = similar(a)
@inline dualarray(a::AbstractArray{<: ForwardDiff.Dual}, b::AbstractArray{<: Real})             = similar(a)
@inline dualarray(a::AbstractArray{<: Real}, b::AbstractArray{<: Real})                         = similar(a)
@inline dualarray(a::AbstractArray{<: Real}, b::AbstractArray{<: ForwardDiff.Dual})             = similar(a, eltype(b))

"""
```
dualvector(a::AbstractVector, b::AbstractVector)
```

has the same behavior of `dualarray` but acts specifically on vectors. This function
is primarily for user convenience.
"""
@inline dualvector(a::AbstractVector{<: ForwardDiff.Dual}, b::AbstractVector{<: ForwardDiff.Dual}) = similar(a)
@inline dualvector(a::AbstractVector{<: ForwardDiff.Dual}, b::AbstractVector{<: Real})             = similar(a)
@inline dualvector(a::AbstractVector{<: Real}, b::AbstractVector{<: Real})                         = similar(a)
@inline dualvector(a::AbstractVector{<: Real}, b::AbstractVector{<: ForwardDiff.Dual})             = similar(a, eltype(b))

# Extend get_tmp(dc::DiffCache, ...) to allow for two arguments
function get_tmp(dc::DiffCache, u1::AbstractArray{T1}, u2::AbstractArray{T2},
                 select::Tuple{Int, Int}) where {T1 <: ForwardDiff.Dual, T2 <: ForwardDiff.Dual}
    if select[1] == 1
        get_tmp(dc, u1)
    elseif select[1] == 2
        get_tmp(dc, u2)
    else
        throw(MethodError("Fourth input argument to get_tmp points to a non-existent cache."))
    end
end
function get_tmp(dc::DiffCache, u1::AbstractArray{T1}, u2::AbstractArray{T2},
                 select::Tuple{Int, Int}) where {T1 <: Number, T2 <: ForwardDiff.Dual}
    if select[1] == 1
        get_tmp(dc, u1)
    elseif select[1] == 2
        get_tmp(dc, u2)
    else
        throw(MethodError("Fourth input argument to get_tmp points to a non-existent cache."))
    end
end
function get_tmp(dc::DiffCache, u1::AbstractArray{T1}, u2::AbstractArray{T2},
                 select::Tuple{Int, Int}) where {T1 <: ForwardDiff.Dual, T2 <: Number}
    if select[1] == 1
        get_tmp(dc, u1)
    elseif select[1] == 2
        get_tmp(dc, u2)
    else
        throw(MethodError("Fourth input argument to get_tmp points to a non-existent cache."))
    end
end

get_tmp(dc::DiffCache, u1::AbstractArray{T1}, u2::AbstractArray{T2}, select::Tuple{Int, Int}) where {T1 <: Number, T2 <: Number} = dc.du

# Extend get_tmp to allow 3 input arguments, only done for the case required by RiskAdjustedLinearizations.jl
function get_tmp(dc::DiffCache, u1::AbstractArray{T1}, u2::AbstractArray{T2},
                 u3::AbstractArray{T3}, select::Tuple{Int, Int}) where {T1 <: ForwardDiff.Dual, T2 <: ForwardDiff.Dual,
                                                                        T3 <: ForwardDiff.Dual}
    if select[1] == 1
        get_tmp(dc, u1)
    elseif select[1] == 2
        get_tmp(dc, u2)
    elseif select[1] == 3
        get_tmp(dc, u3)
    else
        throw(MethodError("Sixth input argument to get_tmp points to a non-existent cache."))
    end
end

get_tmp(dc::DiffCache, u1::AbstractArray{<: Number}, u2::AbstractArray{<: Number}, u3::AbstractArray{<: Number}, select::Tuple{Int, Int}) = dc.du

# Extend get_tmp to allow 4 input arguments, only done for the case required by RiskAdjustedLinearizations.jl
function get_tmp(dc::DiffCache, u1::AbstractArray{T1}, u2::AbstractArray{T2},
                 u3::AbstractArray{T3}, u4::AbstractArray{T4},
                 select::Tuple{Int, Int}) where {T1 <: ForwardDiff.Dual, T2 <: ForwardDiff.Dual,
                                                 T3 <: ForwardDiff.Dual, T4 <: ForwardDiff.Dual}
    if select[1] == 1
        get_tmp(dc, u1)
    elseif select[1] == 2
        get_tmp(dc, u2)
    elseif select[1] == 3
        get_tmp(dc, u3)
    elseif select[1] == 4
        get_tmp(dc, u4)
    else
        throw(MethodError("Sixth input argument to get_tmp points to a non-existent cache."))
    end
end

function get_tmp(dc::DiffCache, u1::AbstractArray{T1}, u2::AbstractArray{T2},
                 u3::AbstractArray{T3}, u4::AbstractArray{T4},
                 select::Tuple{Int, Int}) where {T1 <: Number, T2 <: Number,
                                                 T3 <: Number, T4 <: ForwardDiff.Dual}
    if select[1] == 1
        get_tmp(dc, u1)
    elseif select[1] == 2
        get_tmp(dc, u2)
    elseif select[1] == 3
        get_tmp(dc, u3)
    elseif select[1] == 4
        get_tmp(dc, u4)
    else
        throw(MethodError("Sixth input argument to get_tmp points to a non-existent cache."))
    end
end

get_tmp(dc::DiffCache, u1::AbstractArray{<: Number}, u2::AbstractArray{<: Number}, u3::AbstractArray{<: Number}, u4::AbstractArray{<: Number}, select::Tuple{Int, Int}) = dc.du

"""
```
TwoDiffCache
```

The TwoDiffCache type extends DiffCache from DiffEqBase to permit two Dual Array caches
to permit the case where you need to autodiff w.r.t 2 different lengths of arguments.
For example, suppose we have

```
function f(out, x, y) # in-place function
    out .= x .* y
end

function F(cache::DiffCache, x, y) # wrapper function for f
    f(get_tmp(cache, x), x, y)
    return get_tmp(cache, x)
end

# Instantiate inputs and caches
x = rand(3)
y = rand(3)
cache3 = dualcache(zeros(3), Val{3})
cache6 = dualcache(zeros(3), Val{6})
```

Then the following block of code will work

```
JF3 = (G, x1, y1) -> ForwardDiff.jacobian!(G, z -> F(cache3, z, y1), x1)
JF3(rand(3, 3), x, y)
```

but it may be the case that we also sometimes need to calculate

```
JF6 = (G, x1, y1) -> ForwardDiff.jacobian!(G, z -> F(cache3, z[1:3], z[4:6]), vcat(x1, y1))
JF6(rand(3, 3), x, y)
```

This block of code will error because the chunk size needs to be 6, not 3 here.
Therefore, the correct definition of `JF6` is

```
JF6 = (G, x1, y1) -> ForwardDiff.jacobian!(G, z -> F(cache6, z[1:3], z[4:6]), vcat(x1, y1))
```

Rather than carry around two `DiffCache` objects, it is better to simply add another dual cache
to a single `DiffCache` object. In principle, this code could be generalized to `n` dual caches,
but this would require some additional thought to implement generically.
"""
struct TwoDiffCache{T <: AbstractArray, C1 <: AbstractArray, C2 <: AbstractArray}
    du::T
    dual_du1::C1
    dual_du2::C2
end

function TwoDiffCache(u::AbstractArray{T}, siz, ::Type{Val{chunk_size1}}, ::Type{Val{chunk_size2}}) where {T, chunk_size1, chunk_size2}
    x1 = ArrayInterface.restructure(u, zeros(ForwardDiff.Dual{nothing, T, chunk_size1}, siz...))
    x2 = ArrayInterface.restructure(u, zeros(ForwardDiff.Dual{nothing, T, chunk_size2}, siz...))
    TwoDiffCache(u, x1, x2)
end

twodualcache(u::AbstractArray, N1, N2) = TwoDiffCache(u, size(u), N1, N2)

# get_tmp for AbstractArray cases
function get_tmp(tdc::TwoDiffCache, u1::AbstractArray{T1}, u2::AbstractArray{T2},
                 select::Tuple{Int, Int}) where {T1 <: ForwardDiff.Dual, T2 <: ForwardDiff.Dual}
    if select[1] == 1
        x = reinterpret(T1, select[2] == 1 ? tdc.dual_du1 : tdc.dual_du2)
    elseif select[1] == 2
        x = reinterpret(T2, select[2] == 1 ? tdc.dual_du1 : tdc.dual_du2)
    else
        throw(MethodError("Fourth input argument to get_tmp points to a non-existent cache."))
    end
    return x
end
function get_tmp(tdc::TwoDiffCache, u1::AbstractArray{T1}, u2::AbstractArray{T2},
                 select::Tuple{Int, Int}) where {T1 <: Number, T2 <: ForwardDiff.Dual}
    if select[1] == 1
        x = reinterpret(T1, select[2] == 1 ? tdc.dual_du1 : tdc.dual_du2)
    elseif select[1] == 2
        x = reinterpret(T2, select[2] == 1 ? tdc.dual_du1 : tdc.dual_du2)
    else
        throw(MethodError("Fourth input argument to get_tmp points to a non-existent cache."))
    end
    return x
end
function get_tmp(tdc::TwoDiffCache, u1::AbstractArray{T1}, u2::AbstractArray{T2},
                 select::Tuple{Int, Int}) where {T1 <: ForwardDiff.Dual, T2 <: Number}
    if select[1] == 1
        x = reinterpret(T1, select[2] == 1 ? tdc.dual_du1 : tdc.dual_du2)
    elseif select[1] == 2
        x = reinterpret(T2, select[2] == 1 ? tdc.dual_du1 : tdc.dual_du2)
    else
        throw(MethodError("Fourth input argument to get_tmp points to a non-existent cache."))
    end
    return x
end

# get_tmp for no Dual cases
get_tmp(tdc::TwoDiffCache, u1::AbstractArray, u2::AbstractArray, select::Tuple{Int, Int}) = tdc.du

# Extend get_tmp to allow 4 input arguments, only done for the case required by RiskAdjustedLinearizations.jl
function get_tmp(tdc::TwoDiffCache, u1::AbstractArray{T1}, u2::AbstractArray{T2},
                 u3::AbstractArray{T3}, select::Tuple{Int, Int}) where {T1 <: ForwardDiff.Dual, T2 <: ForwardDiff.Dual,
                                                 T3 <: ForwardDiff.Dual}
    if select[1] == 1
        x = reinterpret(T1, select[2] == 1 ? tdc.dual_du1 : tdc.dual_du2)
    elseif select[1] == 2
        x = reinterpret(T2, select[2] == 1 ? tdc.dual_du1 : tdc.dual_du2)
    elseif select[1] == 3
        x = reinterpret(T3, select[2] == 1 ? tdc.dual_du1 : tdc.dual_du2)
    else
        throw(MethodError("Sixth input argument to get_tmp points to a non-existent cache."))
    end
end

get_tmp(tdc::TwoDiffCache, u1::AbstractArray{<: Number}, u2::AbstractArray{<: Number}, u3::AbstractArray{<: Number}, select::Tuple{Int, Int}) = tdc.du

# Extend get_tmp to allow 4 input arguments, only done for the case required by RiskAdjustedLinearizations.jl
function get_tmp(tdc::TwoDiffCache, u1::AbstractArray{T1}, u2::AbstractArray{T2},
                 u3::AbstractArray{T3}, u4::AbstractArray{T4},
                 select::Tuple{Int, Int}) where {T1 <: ForwardDiff.Dual, T2 <: ForwardDiff.Dual,
                                                 T3 <: ForwardDiff.Dual, T4 <: ForwardDiff.Dual}
    if select[1] == 1
        x = reinterpret(T1, select[2] == 1 ? tdc.dual_du1 : tdc.dual_du2)
    elseif select[1] == 2
        x = reinterpret(T2, select[2] == 1 ? tdc.dual_du1 : tdc.dual_du2)
    elseif select[1] == 3
        x = reinterpret(T3, select[2] == 1 ? tdc.dual_du1 : tdc.dual_du2)
    elseif select[1] == 4
        x = reinterpret(T4, select[2] == 1 ? tdc.dual_du1 : tdc.dual_du2)
    else
        throw(MethodError("Sixth input argument to get_tmp points to a non-existent cache."))
    end
end

function get_tmp(tdc::TwoDiffCache, u1::AbstractArray{T1}, u2::AbstractArray{T2},
                 u3::AbstractArray{T3}, u4::AbstractArray{T4},
                 select::Tuple{Int, Int}) where {T1 <: Number, T2 <: Number,
                                                 T3 <: Number, T4 <: ForwardDiff.Dual}
    if select[1] == 1
        x = reinterpret(T1, select[2] == 1 ? tdc.dual_du1 : tdc.dual_du2)
    elseif select[1] == 2
        x = reinterpret(T2, select[2] == 1 ? tdc.dual_du1 : tdc.dual_du2)
    elseif select[1] == 3
        x = reinterpret(T3, select[2] == 1 ? tdc.dual_du1 : tdc.dual_du2)
    elseif select[1] == 4
        x = reinterpret(T4, select[2] == 1 ? tdc.dual_du1 : tdc.dual_du2)
    else
        throw(MethodError("Sixth input argument to get_tmp points to a non-existent cache."))
    end
end

get_tmp(tdc::TwoDiffCache, u1::AbstractArray{<: Number}, u2::AbstractArray{<: Number}, u3::AbstractArray{<: Number}, u4::AbstractArray{<: Number}, select::Tuple{Int, Int}) = tdc.du

"""
```
ThreeDiffCache
```

The ThreeDiffCache type extends DiffCache from DiffEqBase to permit three Dual Array caches.
"""
struct ThreeDiffCache{T <: AbstractArray, C1 <: AbstractArray, C2 <: AbstractArray, C3 <: AbstractArray}
    du::T
    dual_du1::C1
    dual_du2::C2
    dual_du3::C3
end

function ThreeDiffCache(u::AbstractArray{T}, siz, ::Type{Val{chunk_size1}},
                        ::Type{Val{chunk_size2}}, ::Type{Val{chunk_size3}}) where {T, chunk_size1, chunk_size2, chunk_size3}
    x1 = ArrayInterface.restructure(u, zeros(ForwardDiff.Dual{nothing, T, chunk_size1}, siz...))
    x2 = ArrayInterface.restructure(u, zeros(ForwardDiff.Dual{nothing, T, chunk_size2}, siz...))
    x3 = ArrayInterface.restructure(u, zeros(ForwardDiff.Dual{nothing, T, chunk_size3}, siz...))
    ThreeDiffCache(u, x1, x2, x3)
end

threedualcache(u::AbstractArray, N1, N2, N3) = ThreeDiffCache(u, size(u), N1, N2, N3)

# get_tmp for both AbstractArray
function get_tmp(tdc::ThreeDiffCache, u1::AbstractArray{T1}, u2::AbstractArray{T2},
                 select::Tuple{Int, Int}) where {T1 <: ForwardDiff.Dual, T2 <: ForwardDiff.Dual}
    dual_du = if select[2] == 1
        tdc.dual_du1
    elseif select[2] == 2
        tdc.dual_du2
    elseif select[2] == 3
        tdc.dual_du3
    else
        throw(MethodError("Fourth input argument to get_tmp points to a non-existent cache"))
    end

    if select[1] == 1
        x = reinterpret(T1, dual_du)
    elseif select[1] == 2
        x = reinterpret(T2, dual_du)
    else
        throw(MethodError("Fourth input argument to get_tmp points to a non-existent cache."))
    end
    return x
end
function get_tmp(tdc::ThreeDiffCache, u1::AbstractArray{T1}, u2::AbstractArray{T2},
                 select::Tuple{Int, Int}) where {T1 <: Number, T2 <: ForwardDiff.Dual}
    dual_du = if select[2] == 1
        tdc.dual_du1
    elseif select[2] == 2
        tdc.dual_du2
    elseif select[2] == 3
        tdc.dual_du3
    else
        throw(MethodError("Fourth input argument to get_tmp points to a non-existent cache"))
    end

    if select[1] == 1
        x = reinterpret(T1, dual_du)
    elseif select[1] == 2
        x = reinterpret(T2, dual_du)
    else
        throw(MethodError("Fourth input argument to get_tmp points to a non-existent cache."))
    end
    return x
end
function get_tmp(tdc::ThreeDiffCache, u1::AbstractArray{T1}, u2::AbstractArray{T2},
                 select::Tuple{Int, Int}) where {T1 <: ForwardDiff.Dual, T2 <: Number}
    dual_du = if select[2] == 1
        tdc.dual_du1
    elseif select[2] == 2
        tdc.dual_du2
    elseif select[2] == 3
        tdc.dual_du3
    else
        throw(MethodError("Fourth input argument to get_tmp points to a non-existent cache"))
    end

    if select[1] == 1
        x = reinterpret(T1, dual_du)
    elseif select[1] == 2
        x = reinterpret(T2, dual_du)
    else
        throw(MethodError("Fourth input argument to get_tmp points to a non-existent cache."))
    end
    return x
end

# get_tmp for no Dual cases
get_tmp(tdc::ThreeDiffCache, u1::AbstractArray, u2::AbstractArray, select::Tuple{Int, Int}) = tdc.du

# Extend get_tmp to allow 3 input arguments, only done for the case required by RiskAdjustedLinearizations.jl
function get_tmp(tdc::ThreeDiffCache, u1::AbstractArray{T1}, u2::AbstractArray{T2},
                 u3::AbstractArray{T3}, select::Tuple{Int, Int}) where {T1 <: ForwardDiff.Dual, T2 <: ForwardDiff.Dual,
                                                                        T3 <: ForwardDiff.Dual}
    dual_du = if select[2] == 1
        tdc.dual_du1
    elseif select[2] == 2
        tdc.dual_du2
    elseif select[2] == 3
        tdc.dual_du3
    else
        throw(MethodError("Fourth input argument to get_tmp points to a non-existent cache"))
    end

    if select[1] == 1
        x = reinterpret(T1, dual_du)
    elseif select[1] == 2
        x = reinterpret(T2, dual_du)
    elseif select[1] == 3
        x = reinterpret(T3, dual_du)
    else
        throw(MethodError("Sixth input argument to get_tmp points to a non-existent cache."))
    end
end

get_tmp(tdc::ThreeDiffCache, u1::AbstractArray{<: Number}, u2::AbstractArray{<: Number}, u3::AbstractArray{<: Number}, select::Tuple{Int, Int}) = tdc.du

# Extend get_tmp to allow 4 input arguments, only done for the case required by RiskAdjustedLinearizations.jl
function get_tmp(tdc::ThreeDiffCache, u1::AbstractArray{T1}, u2::AbstractArray{T2},
                 u3::AbstractArray{T3}, u4::AbstractArray{T4},
                 select::Tuple{Int, Int}) where {T1 <: ForwardDiff.Dual, T2 <: ForwardDiff.Dual,
                                                 T3 <: ForwardDiff.Dual, T4 <: ForwardDiff.Dual}
    dual_du = if select[2] == 1
        tdc.dual_du1
    elseif select[2] == 2
        tdc.dual_du2
    elseif select[2] == 3
        tdc.dual_du3
    else
        throw(MethodError("Fourth input argument to get_tmp points to a non-existent cache"))
    end

    if select[1] == 1
        x = reinterpret(T1, dual_du)
    elseif select[1] == 2
        x = reinterpret(T2, dual_du)
    elseif select[1] == 3
        x = reinterpret(T3, dual_du)
    elseif select[1] == 4
        x = reinterpret(T4, dual_du)
    else
        throw(MethodError("Sixth input argument to get_tmp points to a non-existent cache."))
    end
end

function get_tmp(tdc::ThreeDiffCache, u1::AbstractArray{T1}, u2::AbstractArray{T2},
                 u3::AbstractArray{T3}, u4::AbstractArray{T4},
                 select::Tuple{Int, Int}) where {T1 <: Number, T2 <: Number,
                                                 T3 <: Number, T4 <: ForwardDiff.Dual}
    dual_du = if select[2] == 1
        tdc.dual_du1
    elseif select[2] == 2
        tdc.dual_du2
    elseif select[2] == 3
        tdc.dual_du3
    else
        throw(MethodError("Fourth input argument to get_tmp points to a non-existent cache"))
    end

    if select[1] == 1
        x = reinterpret(T1, dual_du)
    elseif select[1] == 2
        x = reinterpret(T2, dual_du)
    elseif select[1] == 3
        x = reinterpret(T3, dual_du)
    elseif select[1] == 4
        x = reinterpret(T4, dual_du)
    else
        throw(MethodError("Sixth input argument to get_tmp points to a non-existent cache."))
    end
end

get_tmp(tdc::ThreeDiffCache, u1::AbstractArray{<: Number}, u2::AbstractArray{<: Number}, u3::AbstractArray{<: Number}, u4::AbstractArray{<: Number}, select::Tuple{Int, Int}) = tdc.du


text:
include("../include.jl")

summary_netlib = Dict{String, Dict{String,problem_summary}}()
#summary_netlib["netlib/well-behaved"] = load("../results/netlib/well-behaved/summary.jld", "summary")
#summary_netlib["netlib/fast-primal"] = load("results/netlib/fast-primal/summary.jld", "summary")
#summary_netlib["netlib/ipopt-no-perturb"] = convert_JuMP(load("../results/netlib/ipopt-no-perturb/summary.jld", "summary"))
#summary_netlib["netlib/ipopt-perturb"] = convert_JuMP(load("../results/netlib/ipopt-perturb/summary.jld", "summary"))

summary_netlib["netlib/well-behaved-2"] = load("../results/netlib/well-behaved-2/summary.jld", "summary")
summary_netlib["netlib/ipopt-mehrotra"] = convert_JuMP(load("../results/netlib/ipopt-mehrotra/summary.jld", "summary"))

for (method_name, sum_data) in summary_netlib
    println(method_name, "=" ,tot_failures(sum_data), " out of ", length(sum_data))
end

summary_netlib = overlap(summary_netlib)
problem_list, method_list = get_lists(summary_netlib)
all_opt_problems = get_all_optimal(summary_netlib, problem_list)
#restrict_problems(summary_netlib,all_opt_problems)

function compute_dual_maxes(summary, folder_name)
    dual_maxes = Dict{String,Dict{String,Float64}}() #Array{Float64,1}
    for (method_name, sum_data) in summary
      dual_maxes[method_name] = Dict{String,Float64}()
      for (problem_name, info) in sum_data
          file = "$folder_name/$method_name/jld/$problem_name.jld"
          hist = load(file, "history")
          dual_hist = compute_dual_hist(hist)
          dual_maxes[method_name][problem_name] = maximum(dual_hist)
      end
    end

    return dual_maxes
end

linestyle_dict = Dict(
  "netlib/well-behaved" => "--",
  "netlib/well-behaved-2" => "--",
  "netlib/ipopt-no-perturb" => "-",
  "netlib/fast-primal" => "-.",
  "netlib/ipopt-perturb" => ":",
  "netlib/ipopt-mehrotra" => ":"
)

label_dic = Dict(
  "netlib/well-behaved" => "Well behaved IPM",
  "netlib/ipopt-no-perturb" => "IPOPT w/o perturb",
  "netlib/fast-primal" => "IPM reduces primal fast",
  "netlib/ipopt-perturb" => "IPOPT w. perturb",
  "netlib/well-behaved-2" => "Well behaved IPM",
  "netlib/ipopt-mehrotra" => "IPOPT mehrotra"
)

dual_maxes_all = compute_dual_maxes(summary_netlib, "../results")

dual_maxes_restricted = Dict{String,Array{Float64,1}}()
for method_name in method_list
  dual_maxes_restricted[method_name] = Array{Float64,1}()
  for problem_name in all_opt_problems
   push!(dual_maxes_restricted[method_name], dual_maxes_all[method_name][problem_name])
  end
end

using PyPlot

min_y = 1.0

y_vals = collect(1:length(all_opt_problems)) / length(all_opt_problems)
for (method_name, val) in dual_maxes_restricted
  semilogy(y_vals, sort(dual_maxes_restricted[method_name]),linestyle=linestyle_dict[method_name], label=label_dic[method_name], color="black", basey=10)
end

xlabel("proportion of problems")
ylabel("maximum dual value over all iterates")

legend()

#df = DataFrame(problem_name = problem_list, feasible11 = status11, problem_name2 = names2, feasible21 = status21 )

# IPOPT 11/29 failures
# us 4/29 failures

con_vio = Dict{String,Array{Float64,1}}()
for folder_name in method_list
  con_vio[folder_name] = Array{Float64,1}()
  for problem_name in problem_list
      results = load("results/$folder_name/jld/$problem_name.jld", "history")
      push!(con_vio[folder_name], results[end].con_vio)
  end
end

ratios = Dict()
y_vals = collect(1:length(con_vio["netlib"])) / length(con_vio["netlib"])
for (method_name, val) in dual_maxes
  #ratios[method_name] = dual_maxes[method_name] ./ best_duals
  semilogx(sort(con_vio[method_name]), y_vals, label=method_name, basex=10)
  #min_y = min(min_y,sum(ratios[method_name] .== 1.0) / length(best_duals))
end

xlabel("con violation")
ylabel("proportion of problems")

legend()

####################################
### PLOT ONE PARTICULAR PROBLEM
####################################

field_list = Dict("max dual variable" => :y_norm, "primal feasibility" => :con_vio, "dual feasibility" => :norm_grad_lag, "complementarity" => :comp);
line_style_list = Dict("max dual variable" => "-", "primal feasibility" => "--", "dual feasibility" => "-.", "complementarity" => ":");

method_name = "well-behaved"
problem_name = "ADLITTLE"
file_wb = "results/netlib/$method_name/jld/$problem_name.jld"
hist_wb = load(file_wb, "history")
hist_by_it_wb = summarize_by_iteration(hist_wb)
data_wb = aggregate_summary_by_it_for_plots(hist_by_it_wb, field_list)

method_name = "ipopt-no-perturb"
file_ipnp = "results/netlib/$method_name/jld/$problem_name.jld"
hist_ipnp = load(file_ipnp, "history")
hist_by_it_ipnp = summarize_by_iteration(hist_ipnp)
data_ipnp = aggregate_summary_by_it_for_plots(hist_by_it_ipnp, field_list)

method_name = "ipopt-perturb"
file_ipp = "results/netlib/$method_name/jld/$problem_name.jld"
hist_ipp = load(file_ipp, "history")
hist_by_it_ipp = summarize_by_iteration(hist_ipp)
data_ipp = aggregate_summary_by_it_for_plots(hist_by_it_ipp, field_list)

##
## DO THE PLOTTING
##

function plot_single_example_compare_iterates(data, field_list, line_style_list)
  for (label, field) in field_list
    semilogy(data[field],label=label,linestyle=line_style_list[label], color="black")
  end
end

ylims = [1e-9, 1e12]

#plt[:figlegend]( ["-"], ["test"], loc = (0.5, 0), ncol=5 )
#plt[:figlegend]( ["-"], ["test"], loc = (0.0, 0), ncol=1 )

subplot(131)
plot_single_example_compare_iterates(data_wb, field_list, line_style_list)

title("A well-behaved IPM")
ylabel("L-infinity norm")
xlabel("iterations")
ax = gca()
ax[:set_ylim](ylims)

subplot(132)
plot_single_example_compare_iterates(data_ipnp, field_list, line_style_list)

title("IPOPT w/o perturbation")
#ylabel("L-infinity norm")
xlabel("iterations")
ax = gca()
ax[:set_ylim](ylims)
ax[:set_yticklabels]([])

subplot(133)
plot_single_example_compare_iterates(data_ipp, field_list, line_style_list)

title("IPOPT w. perturbation")
xlabel("iterations")
ax = gca()
ax[:set_ylim](ylims)
ax[:set_yticklabels]([])

legend()
ax[:legend]()[:draggable]()

#suptitle("Comparsion of algorithms on $problem_name")
#ttl = ax[:title]
#ttl[:set_position]([.5, 1.05])
#plt[:tight_layout]()

#ax[:set_position]([.5, 1.05])

#########################
## PLOT ITERATION COUNTS
#########################

using PyPlot

its, best, ratios, times = compute_its_etc(summary_netlib,MAX_IT=100)
plot_iteration_ratios(its, best, ratios)


text:
using Test, Omniscape, Circuitscape

### Unit tests for components
## source target matching
source_subset = [1.0 1 1 1 1 1 1 1 1 1; # 6x10 array
                   1 1 1 1 1 1 1 1 1 1;
                   1 1 1 1 1 1 1 1 1 1;
                   1 1 1 1 1 1 1 1 1 1;
                   1 1 1 1 1 1 1 1 1 1;
                   1 1 1 1 1 1 1 1 1 1]

con1pres =   [1.0 4 3 2 7 6 9 2 3 0;
                4 2 6 5 3 9 4 5 3 3;
                1 1 5 2 7 7 8 8 9 1;
                8 7 6 5 5 5 0 0 1 1;
                3 5 4 4 4 8 5 0 1 8;
                1 2 8 2 1 6 5 5 8 7]

con1fut =   [2.0 1 3 2 7 6 7 1 3 0;
               4 2 4 5 3 9 4 5 5 3;
               3 8 3 2 6 3 8 5 9 1;
               8 7 6 5 1 5 0 0 1 1;
               3 1 4 5 4 8 5 7 1 8;
               5 2 3 2 1 6 4 5 3 7]
x = 5
y = 4
con1_lower = -1.0
con1_upper = 1.0
Omniscape.source_target_match!(source_subset,
                     1,
                     con1pres,
                     con1fut,
                     Array{Float64, 2}(undef, 1, 1),
                     Array{Float64, 2}(undef, 1, 1),
                     "within",
                     "within",
                     con1_lower, # lower bound to be within for comparisome
                     con1_upper,
                     0.0,
                     0.0,
                     y,
                     y,
                     x,
                     x,
                     1,
                     6,
                     1,
                     10
                     )
target_val = con1fut[y, x]
# Make sure no present day vals in con1 are outside of range
@test sum((con1pres[source_subset .== 1] .< (target_val + con1_lower)) .|
    (con1pres[source_subset .== 1] .> (target_val + con1_upper))) == 0
@info "conditional connectivity tests passed"

## Check that targets are IDed properly
sources_raw = Circuitscape.read_raster("input/source.asc", Float64)[1]
int_arguments = Dict{String, Int64}()
int_arguments["block_size"] = 7
int_arguments["block_radius"] = 3 # must be (size - 1) / 2
int_arguments["nrows"] = size(sources_raw)[1]
int_arguments["ncols"] = size(sources_raw)[2]
targets = Omniscape.get_targets(sources_raw,
                                int_arguments,
                                Float64)

n_targets = floor(int_arguments["nrows"] / int_arguments["block_size"]) *
    floor(int_arguments["ncols"] / int_arguments["block_size"])

# Correct number of targets
@test size(targets)[1] == n_targets # would be less than if any of the sources
                                    # had 0 strength, but sources.asc does
                                    # not have 0's

# correct source strength for a target
block_sources = sources_raw[Int(targets[1,2] - int_arguments["block_radius"]):Int(targets[1,2] + int_arguments["block_radius"]),
                            Int(targets[1,1] - int_arguments["block_radius"]):Int(targets[1,1] + int_arguments["block_radius"])]
@test targets[1,3] ≈ sum(block_sources)
@info "target tests passed"

### Syntax tests for run_omniscape()
l, f, p = run_omniscape("input/config4.ini")
g = run_omniscape("input/config5.ini")
h = run_omniscape("input/config6.ini")
a, b, c = run_omniscape("input/config.ini")
a1, b1, c1 = run_omniscape("input/config_32bit.ini")
q, e = run_omniscape("input/config3.ini")
d = run_omniscape("input/config2.ini")
d_1 = run_omniscape("input/config2.ini")
d_2 = run_omniscape("input/config2.ini")

@test typeof(f) == Array{Float64,2}
@test typeof(g) == Array{Float64,2}
@test typeof(h) == Array{Float64,2}
@test typeof(a) == Array{Float64,2}
@test typeof(b) == Array{Float64,2}
@test typeof(c) == Array{Float64,2}
@test typeof(a1) == Array{Float32,2}
@test typeof(b1) == Array{Float32,2}
@test typeof(c1) == Array{Float32,2}
@test typeof(d) == Array{Float64,2}
@test typeof(e) == Array{Float64,2}
@test a ≈ d #parallel and serial produce same result

# Single and double produce similar results
@test a ≈ a1
@test b ≈ b1
@test c ≈ c1

# Test error throws
@info "Testing error throws"
@test run_omniscape("input/config7.ini") == nothing
@test run_omniscape("input/bad_config.ini") == nothing

GC.gc()

rm("test1", recursive = true)
rm("test1_32bit", recursive = true)
rm("test2", recursive = true)
rm("test2_1", recursive = true)
rm("test2_2", recursive = true)
rm("test3", recursive = true)
rm("test4", recursive = true)
rm("test5", recursive = true)
rm("test6", recursive = true)

@info "run_omniscape tests passed"