text:
module SysImageBuilder
using ..DevEnvConfig: devEnvConfig_pkg_path
using ..RuntimeContext
using Dates
using UUIDs
using Pkg

export ImageSpec
export name, path, size, date, exists, delete, build

import Base: repr, size

const JULIA_TARGET = "native"

const VERSION_PackageCompiler = "1.3.0" # must be a string
const ENVPATH_SEPARATOR = @static Sys.iswindows() ? ";" : ":"
const IMAGE_FILENAME = "sys.jsi"

const BOOTSTRAP_LOCAL = true
const BOOTSTRAP_PATH = devEnvConfig_pkg_path()
const BOOTSTRAP_SPEC = """PackageSpec(path = raw"$BOOTSTRAP_PATH")"""

#region ImageSpec

function isvalidname(name::String)
    check(c::AbstractChar) = isletter(c) || isdigit(c) || c == '_' || c == '-'
    return isascii(name) && all(check, name)
end

function checkvalidname(name::String)
    isvalidname(name) || error("Invalid image name '$name'. Only ascii letters, digits, '_', '-' are allowed.")
end

struct ImageSpec
    name::String
    extends:: Union{ImageSpec,Nothing}
    function ImageSpec(name::String)
        checkvalidname(name)
        new(name, nothing)
    end
    function ImageSpec(name::String, extends::ImageSpec)
        checkvalidname(name)
        new(name, extends)
    end
end

function repr(s::ImageSpec)
    # Must be a true repr. Used to generate code.
    # TODO: Test with eval(Meta.parse(repr(a))) == a
    if isnothing(s.extends)
        """ImageSpec("$(s.name)")"""
    else
        """ImageSpec("$(s.name)",$(repr(s.extends)))"""
    end
end

name(s::ImageSpec)   = s.name
path(s::ImageSpec)   = joinpath(RuntimeContext.builddevenv_path(), name(s), IMAGE_FILENAME)
size(s::ImageSpec)   = round(filesize(path(s)) / (1 << 20); digits=0) # 0 MB if not exists
date(s::ImageSpec)   = Dates.unix2datetime(mtime(path(s)))            # epoc if not exists
exists(s::ImageSpec) = isfile(path(s))

#endregion

#region Julia isolated invokation

function julia_invoke_cmd()
    jr = escape_string(RuntimeContext.julia_executable())
    return `$jr --startup-file=no --history-file=no --warn-overwrite=yes --cpu-target $JULIA_TARGET --quiet -t 1`
end

function julia(;
    project::Union{String,Nothing} = nothing,
    source::Union{String,Nothing} = nothing,
    code::Union{String,Nothing} = nothing,
    trace::Union{String,Nothing} = nothing,
    indir::Union{String,Nothing} = nothing,
    image::Union{ImageSpec,Nothing} = nothing,
    optlvl::Int = 2,
    dbglvl::Int = 1,
    debug::Bool = false,
    )
    nrmpath(s) = isnothing(s) ? nothing : realpath(expanduser(s))
    escape(s) = isnothing(s) ? nothing : escape_string(s)
    addarg(cmd::Cmd, arg::Cmd) = Cmd(Cmd(vcat(cmd.exec, arg.exec)), cmd.ignorestatus, cmd.flags, cmd.env, cmd.dir)

    addproj(cmd::Cmd, ::Nothing) = cmd
    addproj(cmd::Cmd, val::String) = addarg(cmd, `--project=$val`)

    addflag(cmd::Cmd, ::String, ::Nothing) = cmd
    addflag(cmd::Cmd, flag::String, val::String) = addarg(cmd, `$flag $val`)
    addflag(cmd::Cmd, flag::String, val::Int) = addarg(cmd, `$flag $val`)

    # FIXME: https://docs.julialang.org/en/v1/manual/environment-variables/#JULIA_LOAD_PATH
    # "JULIA_LOAD_PATH"  => ":" as per the doc. After test, must use the ENVPATH_SEPARATOR of the platform

    dir = isnothing(indir) ? pwd() : nrmpath(indir)
    img = isnothing(image) || isnothing(image.extends) ? nothing : nrmpath(path(image.extends))
    isnothing(img) || isfile(img) || error("Internal error. Image $image depends on $img which does not exists.")

    # We will inherit the env. var. of the exec system call (e.g. all USER/MACHINE scoped one on Windows)
    # We override the variable that may cause issue. TODO: We probably need to add more here...
    # JULIA_PKG_DEVDIR, ...
    env = Dict(
        "JULIA_DEPOT_PATH" => join(RuntimeContext.depotdevenv(), ENVPATH_SEPARATOR),
        "JULIA_LOAD_PATH"  => ENVPATH_SEPARATOR,
    )

    cmd = julia_invoke_cmd()
    cmd = addproj(cmd, nrmpath(project))
    cmd = addflag(cmd, "--sysimage", escape(img))
    cmd = addflag(cmd, "--trace-compile", escape(trace))

    cmd = addflag(cmd, "-O", optlvl)
    cmd = addflag(cmd, "-g", dbglvl)
    cmd = addflag(cmd, "-e", code)
    cmd = addflag(cmd, "--", nrmpath(source))
    cmd = Cmd(cmd; env, dir)

    if !isnothing(code)
        printstyled("Julia execution code\n"; color=:red)
        printstyled(code; color=:blue)
        print('\n')
    end

    printstyled("Julia execution starts...\n"; color=:red)
    !debug || dump(cmd)
    run(cmd)
    printstyled("Julia execution completed.\n"; color=:red)
end

#endregion

#region Project generation

function create_project_lazy(path::String, name::String)
    path_proj = joinpath(path, name)
    file_proj = joinpath(path_proj, "Project.toml")
    isfile(file_proj) && return path_proj

    ver = Base.VERSION
    mkpath(path_proj)
    open(file_proj, "w") do io
        print(io, """
        name = "$(name)"
        uuid = "$(uuid4())"
        authors = [""]
        version = "0.1.0"

        [compat]
        julia = "$(ver.major).$(ver.minor)"
    """)
    end
    path_src = joinpath(path_proj, "src")
    mkpath(path_src)
    open(joinpath(path_src, "$name.jl"), "w") do io
        print(io, """
        module $(name)
        end
        """)
    end
    return path_proj
end

function setup_project_builder_lazy()
    path_root = RuntimeContext.builddevenv_path()
    path_proj = joinpath(path_root, "Builder")
    file_proj = joinpath(path_proj, "Project.toml")
    isfile(file_proj) && return path_proj

    create_project_lazy(path_root, "Builder")
    julia(;project=path_proj, code="""
        using Pkg
        using UUIDs
        compiler = PackageSpec("PackageCompiler", UUID("9b87118b-4619-50d2-8e1e-99f35a4d4d9d"), v"$VERSION_PackageCompiler")
        self = $BOOTSTRAP_SPEC
        Pkg.add(compiler)
        $(BOOTSTRAP_LOCAL ? "Pkg.develop(self)" : "Pkg.add(self)")
        println("Created project \$(Base.ACTIVE_PROJECT[])")
    """)

    return path_proj
end

#endregion

function build_init(s::ImageSpec)
    opt = Base.JLOptions()
    default_jsi = !convert(Bool, opt.image_file_specified)
    current_img = unsafe_string(opt.image_file)
    printstyled("PackageCompiler invoked in environment:\n";color=:yellow)
    printstyled("           PWD: $(pwd())\n"; color=:blue)
    printstyled("ACTIVE_PROJECT: $(Base.ACTIVE_PROJECT[])\n"; color=:blue)
    printstyled("     LOAD_PATH: $(Base.load_path()[1])\n"; color=:blue)
    printstyled("    DEPOT_PATH: $(Base.DEPOT_PATH[1])\n"; color=:blue)
    printstyled("  TARGET IMAGE: $(path(s))\n";   color=:blue)
    printstyled(" CURRENT IMAGE: $current_img\n"; color=:blue)
    printstyled("   DEFAULT JSI: $default_jsi\n"; color=:blue)

    if !isnothing(s.extends) && current_img != path(s.extends)
        error("Internal error. Current image should be $(path(s.extends))")
    end
end

function build_finish(s::ImageSpec)
    printstyled("Build completed: $(size(s)) MB, $(date(s))\n";color=:yellow)
end

function build_stdlib_min(create::Function, s::ImageSpec)
    # If we use [filter_stdlibs=true] we can create the absolute minimal image.
    # - It contains only Core+Base, the loaded_modules at Core,Base,Main
    # - There is a very limited REPL (not color, etc.)
    # - Thre size is about 63M
    # By including the StdLib the size rise to about 100M (without any precomp)

    packages = Symbol[] # none !
    params = (
        sysimage_path=IMAGE_FILENAME,
        project=("./Build"),
        precompile_execution_file = String[],
        precompile_statements_file = String[],
        incremental = false,
        filter_stdlibs = false,
        replace_default = false,
        version = Base.VERSION,
        cpu_target = "native",
    )

    build_init(s)
    if exists(s)
        printstyled("Building Base+Stdlib: image file exists. Doing nothing.\n";color=:light_red)
    else
        printstyled("Building Base+Stdlib from scratch. No precompilation.\n";color=:yellow)
        invoke(create, Tuple{Vector{Symbol}}, packages; params...) # no extra packages
    end
    build_finish(s)
end

function build_stdlib(create::Function, s::ImageSpec)
    packages = Symbol[] # none !
    params = (
        sysimage_path=IMAGE_FILENAME,
        project=("./Build"),
        precompile_execution_file = String[],
        precompile_statements_file = String["stdlib.jl"],
        incremental = true,
        filter_stdlibs = false,
        replace_default = false,
        version = Base.VERSION,
        cpu_target = "native",
    )

    build_init(s)
    if exists(s)
        printstyled("Building Base+Stdlib: image file exists. Doing nothing.\n";color=:light_red)
    else
        printstyled("Building Base+Stdlib from scratch. No precompilation.\n";color=:yellow)
        invoke(create, Tuple{Vector{Symbol}}, packages; params...) # no extra packages
    end
    build_finish(s)
end

function delete(s::ImageSpec)
    build_path = dirname(path(s))
    contains(build_path, RuntimeContext.DEVENV_DIRNAME) || error("delete")
    rm(build_path; recursive=true, force=true)
end

#  $ jr -J sys.jsi
#  - does not lock sys.jsi
#  - it is possible to delete it will julia is running
#  - but some operation (at least using a package,, precomp check) try to access it
#  => [ Info: Precompiling BenchmarkTools [6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf]
#     ERROR: could not load library "C:\Users\Yoann\.julia\devenvmgr\v1.6-x86_64-w64-mingw32-haswell\stdlib-min\sys.jsi"
#     The specified module could not be found.
#     ERROR: IOError: write: broken pipe (EPIPE)
#

function build(s::ImageSpec)
    build_path = dirname(path(s))

 # FIXME debug helper
 #   ispath(build_path) && error("Image $(name(s)) exists already at $(path(s)). Call delete(..) before build(...) if needed.")
    mkpath(build_path)

    builder = setup_project_builder_lazy()
    create_project_lazy(build_path, "Build")

    call =""
    if name(s) == "stdlib-min"
        call = "build_stdlib_min"
    elseif name(s) == "stdlib"
        call = "build_stdlib"
    else
        error("unknow image")
    end

    julia(project=builder, image=s, indir=dirname(path(s)), code="""
        using PackageCompiler
        using DevEnvConfig.SysImageBuilder
        image = $(repr(s))
        ENV["JULIA_DEBUG"]="all"
        SysImageBuilder.$call(PackageCompiler.create_sysimage, image)
    """)
end

# julia> build(ImageSpec("stdlib-min"))
# Julia execution code
#     using PackageCompiler
#     using DevEnvConfig.SysImageBuilder
#     image = ImageSpec("stdlib-min")
#     SysImageBuilder.build_stdlib_min(PackageCompiler.create_sysimage, image)
# Julia execution starts...
# PackageCompiler invoked in environment:
#            PWD: C:\Users\Yoann\.julia\devenvmgr\v1.6-x86_64-w64-mingw32-haswell\stdlib-min
# ACTIVE_PROJECT: C:\Users\Yoann\.julia\devenvmgr\v1.6-x86_64-w64-mingw32-haswell\Builder
#      LOAD_PATH: C:\Users\Yoann\.julia\devenvmgr\v1.6-x86_64-w64-mingw32-haswell\Builder\Project.toml
#     DEPOT_PATH: C:\Users\Yoann\.julia\devenvmgr\depot
#   TARGET IMAGE: C:\Users\Yoann\.julia\devenvmgr\v1.6-x86_64-w64-mingw32-haswell\stdlib-min\sys.jsi
#  CURRENT IMAGE: C:\Users\Yoann\AppData\Local\Programs\Julia-1.6.1\lib\julia\sys.dll
#    DEFAULT JSI: true
# Building Base+Stdlib from scratch. No precompilation.
#     Updating registry at `C:\Users\Yoann\.julia\devenvmgr\depot\registries\General`
#     Updating git-repo `https://github.com/JuliaRegistries/General.git`
#   No Changes to `C:\Users\Yoann\.julia\devenvmgr\v1.6-x86_64-w64-mingw32-haswell\stdlib-min\Build\Project.toml`
#   No Changes to `C:\Users\Yoann\.julia\devenvmgr\v1.6-x86_64-w64-mingw32-haswell\stdlib-min\Build\Manifest.toml`
# [ Info: PackageCompiler: creating base system image (incremental=false)...
# [ Info: PackageCompiler: creating system image object file, this might take a while...
# Build completed: 99.0 MB, 2021-09-07T14:12:01.292
# Julia execution completed.

end # module SysImageBuilder


text:
@testset "Neighbor search" begin
  @testset "BallSearch" begin
    𝒟 = CartesianGrid((10, 10), T.((-0.5,-0.5)), T.((1.0,1.0)))

    S = BallSearch(𝒟, NormBall(T(1)))
    n = search(P2(0,0), S)
    @test Set(n) == Set([1,2,11])
    n = search(P2(9,0), S)
    @test Set(n) == Set([9,10,20])
    n = search(P2(0,9), S)
    @test Set(n) == Set([91,81,92])
    n = search(P2(9,9), S)
    @test Set(n) == Set([100,99,90])

    S = BallSearch(𝒟, NormBall(T(√2+eps(T))))
    n = search(P2(0,0), S)
    @test Set(n) == Set([1,2,11,12])
    n = search(P2(9,0), S)
    @test Set(n) == Set([9,10,19,20])
    n = search(P2(0,9), S)
    @test Set(n) == Set([81,82,91,92])
    n = search(P2(9,9), S)
    @test Set(n) == Set([89,90,99,100])

    # non MinkowskiMetric example
    𝒟 = CartesianGrid((360,180), T.((0.0,-90.0)), T.((1.0,1.0)))
    S = BallSearch(𝒟, NormBall(T(150), Haversine(T(6371))))
    n = search(P2(0,0), S)
    @test Set(n) == Set([32041, 32400, 32401, 32760])
  end

  @testset "KNearestSearch" begin
    𝒟 = CartesianGrid((10,10), T.((-0.5,-0.5)), T.((1.0,1.0)))
    S = KNearestSearch(𝒟, 3)
    n = search(P2(0,0), S)
    @test Set(n) == Set([1,2,11])
    n = search(P2(9,0), S)
    @test Set(n) == Set([9,10,20])
    n = search(P2(0,9), S)
    @test Set(n) == Set([91,81,92])
    n = search(P2(9,9), S)
    @test Set(n) == Set([100,99,90])
  end

  @testset "KBallSearch" begin
    𝒟 = CartesianGrid((10,10), T.((-0.5,-0.5)), T.((1.0,1.0)))

    s = KBallSearch(𝒟, 10, NormBall(T(100)))
    n = search(P2(5,5), s)
    @test length(n) == 10

    s = KBallSearch(𝒟, 10, Ellipsoid([100., 100.], [0.], convention=TaitBryanIntr))
    n = search(P2(5,5), s)
    @test length(n) == 10

    s = KBallSearch(𝒟, 10, NormBall(T(1)))
    n = search(P2(5,5), s)
    @test length(n) == 5
    @test n[1] == 56

    mask = trues(nelements(𝒟))
    mask[56] = false
    n = search(P2(5,5), s, mask=mask)
    @test length(n) == 4
    n = search(P2(-0.2,-0.2), s)
    @test length(n) == 1
    n = search(P2(-10,-10), s)
    @test length(n) == 0
  end

  @testset "BoundedSearch" begin
    𝒟 = CartesianGrid((10,10), T.((-0.5,-0.5)), T.((1.0,1.0)))
    S1 = BallSearch(𝒟, NormBall(T(5)))
    S2 = KNearestSearch(𝒟, 10)
    B1 = BoundedSearch(S1, 5)
    B2 = BoundedSearch(S2, 5)
    p = centroid(𝒟, rand(1:100))
    n = search(p, B1)
    @test length(n) == 5
    p = centroid(𝒟, rand(1:100))
    n = search(p, B2)
    @test length(n) == 5
  end
end

ext: jl


text:
#   Copyright (c) 2020-2022 Velexi Corporation
#
#   Licensed under the Apache License, Version 2.0 (the "License");
#   you may not use this file except in compliance with the License.
#   You may obtain a copy of the License at
#
#       http://www.apache.org/licenses/LICENSE-2.0
#
#   Unless required by applicable law or agreed to in writing, software
#   distributed under the License is distributed on an "AS IS" BASIS,
#   WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
#   See the License for the specific language governing permissions and
#   limitations under the License.

"""
Pseudoscalar.jl defines the Pseudoscalar type and core methods
"""

# --- Exports

# Types
export Pseudoscalar

# --- Type definitions

"""
    struct Pseudoscalar{T<:AbstractFloat} <: AbstractBlade

Pseudoscalar (an ``(n-1)``-blade) represented with the floating-point precision
of type `T`. The `basis` for a Pseudoscalar is the standard basis for ``R^n``.
The norm and orientation of a Pseudoscalar are encoded in its `value`. The
norm of a Pseudoscalar is equal to `abs(value)` and the orientation of a
Pseudoscalar relative to the standard basis for ``R^n`` is equal to
`sign(value)`.
"""
struct Pseudoscalar{T<:AbstractFloat} <: AbstractBlade{T}
    #=
      Fields
      ------
      * `dim`: the dimension of the space that the blade is embedded in

      * `value`: the value of the pseudoscalar
    =#
    dim::Int
    value::T

    """
    Construct a Pseudoscalar for a geometric algebra in `dim` dimensions having
    the specified `value`.
    """
    Pseudoscalar{T}(dim::Integer, value::Real) where {T<:AbstractFloat} =
        dim <= 0 ?
            throw(ArgumentError("`dim` must be positive")) :
            value == 0 ? zero(Pseudoscalar{T}) : new(dim, T(value))
end

"""
    Pseudoscalar(dim::Integer, value::Real)

Construct a Pseudoscalar for a geometric algebra in `dim` dimensions having the
specified `value`.

When the precision is not specified, the following rules are applied to set
the precision of the Pseudoscalar.

* If `value` is a floating-point value, the precision of the Pseudoscalar is
  inferred from the precision of `value`.

* If `value` is an integer, the precision of the Pseudoscalar defaults to
  `Float64`.
"""
Pseudoscalar(dim::Integer, value::AbstractFloat) =
    Pseudoscalar{typeof(value)}(dim, value)

Pseudoscalar(dim::Integer, value::Integer) = Pseudoscalar(dim, Float64(value))

"""
    Pseudoscalar(B::Pseudoscalar{T};
                 value::Real=value(B)) where {T<:AbstractFloat}

Copy constructor. Construct a Pseudoscalar representing the same space as
`B` having the specified `value`.
"""
Pseudoscalar(B::Pseudoscalar; value::Real=value(B)) =
    Pseudoscalar{typeof(B.value)}(dim(B), value)

# --- Method definitions

"""
    value(B::Pseudoscalar)::AbstractFloat

Return the value of `B` (with the same precision as `B`).
"""
value(B::Pseudoscalar) = B.value

# --- Method definitions for AbstractBlade interface functions

import LinearAlgebra.I

inv(B::Pseudoscalar) =
    mod(grade(B), 4) < 2 ?
        Pseudoscalar(B, value=1 / value(B)) :
        Pseudoscalar(B, value=-1 / value(B))

grade(B::Pseudoscalar) = B.dim

basis(B::Pseudoscalar) = LinearAlgebra.I

volume(B::Pseudoscalar) = value(B)

# --- Method definitions for AbstractMultivector interface functions

dim(B::Pseudoscalar) = B.dim

-(B::Pseudoscalar) = Pseudoscalar(B, value=-value(B))

Base.reverse(B::Pseudoscalar) =
    mod(grade(B), 4) < 2 ?  B : Pseudoscalar(B, value=-value(B))

dual(B::Pseudoscalar) = Scalar(value(B))

# --- Comparison methods

==(B::Pseudoscalar, C::Pseudoscalar) =
    (dim(B) == dim(C)) && (value(B) == value(C))

==(B::Pseudoscalar, C::Vector) =
    (dim(B) == 1) && (length(C) == 1) && (value(B) == C[1])

==(B::Vector, C::Pseudoscalar) = (C == B)

isapprox(B::Pseudoscalar{T1}, C::Pseudoscalar{T2};
  atol::Real=0,
  rtol::Real=atol>0 ? 0 : max(√eps(T1), √eps(T2))) where {T1<:AbstractFloat,
                                                          T2<:AbstractFloat} =
    (dim(B) == dim(C)) && isapprox(value(B), value(C), atol=atol, rtol=rtol)

isapprox(B::Pseudoscalar{T1}, C::Vector{T2};
  atol::Real=0,
  rtol::Real=atol>0 ? 0 : max(√eps(T1), √eps(T2))) where {T1<:AbstractFloat,
                                                          T2<:AbstractFloat} =
    (dim(B) == 1) && (length(C) == 1) && isapprox(value(B), C[1], atol=atol, rtol=rtol)

isapprox(B::Vector{T1}, C::Pseudoscalar{T2};
  atol::Real=0,
  rtol::Real=atol>0 ? 0 : max(√eps(T1), √eps(T2))) where {T1<:AbstractFloat,
                                                          T2<:AbstractFloat} =
    isapprox(C, B, atol=atol, rtol=rtol)

# --- Utility methods

convert(::Type{T}, B::Pseudoscalar) where {T<:AbstractFloat} =
    T == typeof(value(B)) ? B : Pseudoscalar{T}(dim(B), value(B))


text:
import Pkg
Pkg.add("Documenter")
Pkg.add("Distributions")
Pkg.add("Measurements")

using Documenter, Jeff, Distributions, Measurements

makedocs(sitename="Jeff.jl",
         modules=[Jeff],
         pages = [
             "Home" => "index.md",
             "API" => [
                "reflect" => "reflect.md",
                "data" => "data.md",
                "model" => "model.md",
                "objective" => "objective.md",
             ]
         ])

deploydocs(repo = "github.com/arm61/Jeff.jl.git")


text:
#=
A place for general utilities that don't fit in elsewhere.

=#

# Traverses an expression, replacing variable symbols with the value arrays
function replace_var_symbols_with_values(ex)
    if typeof(ex) == Expr
        if ex.head === :call # function symbols (like :sin in sin(x)) don't need to be processed
            for i=2:length(ex.args)
                ex.args[i] = replace_var_symbols_with_values(ex.args[i]);
            end
        else
            for i=1:length(ex.args)
                ex.args[i] = replace_var_symbols_with_values(ex.args[i]);
            end
        end

    elseif typeof(ex) == Symbol
        # Check if it is a variable and substitute if needed
        var_index = 0;
        for v in variables
            if ex === v.symbol
                var_index = v.index;
            end
        end
        if var_index > 0
            if variables[var_index].type == SCALAR
                newex = :(Finch.variables[$var_index].values[node_index]);
            else
                newex = :(Finch.variables[$var_index].values[:,node_index]);
            end
            ex = newex;
        end
    end

    return ex;
end

max_stars_repo_name: paralab/Finch

max_issues_repo_name: paralab/Finch

max_forks_repo_name: paralab/Finch


text:
const NucQualityDict = Dict{UInt8,Int}

# dimensions: nuc, pos
# where nuc is A,C,G or T
const NucQualityMap = Matrix{NucQualityDict}

# init Dict for each nuc and position
nucqualitymap(len::Int) = NucQualityDict[NucQualityDict() for n=1:4, i=1:len]

const NucQualityAccumulator = Vector{NucQualityMap}

nucqualityaccumulator(lengths::Vector{Int}) = NucQualityMap[nucqualitymap(l) for l in lengths]
nucqualityaccumulator(len::Int) = nucqualityaccumulator([len])

# increase the count of a nucleotide quality entry
function inc!(m::NucQualityMap,refPos::Int,nucInd::Int,q::UInt8)
	d = m[nucInd,refPos] # the dict for this position and nucleotide
	d[q] = get!(d,q,0)+1 # is there a better way to increase a value for a key?
end
function inc!(a::NucQualityAccumulator,refInd::Int,refPos::Int,nucInd::Int,q::UInt8)
	inc!(a[refInd],refPos,nucInd,q)
end

# qualityfilter! replaces all quality values below threshold with 0
function qualityfilter!(d::NucQualityDict, threshold::Integer, scratch::Vector{UInt8}=Vector{UInt8}())
	isempty(d) && return
	resize!(scratch,0)

	# get all that should be filtered
	for k in keys(d)
		if k<threshold
			push!(scratch,k)
		end
	end

	# for each filtered entry
	for k in scratch
		pop!(d,k)
	end
end
function qualityfilter!(m::NucQualityMap, threshold::Integer, scratch::Vector{UInt8}=Vector{UInt8}())
	for d in m
		qualityfilter!(d,threshold,scratch)
	end
end
function qualityfilter!(a::NucQualityAccumulator, threshold::Integer, scratch::Vector{UInt8}=Vector{UInt8}())
	for m in a
		qualityfilter!(m,threshold,scratch)
	end
end

function processread!(nqa::NucQualityAccumulator, r::BamRead)
	refInd = refID(r)+1 # 0 to 1-based index

	for c in r
		o = op(c)

		if o==CIGAR_M || o==CIGAR_X || o==CIGAR_Eq # the only chunks that are aligned to both read and reference
			refPos = refpos(c)
			cSeq = seq(r,c)
			cQual = qual(r,c)

			# for all nucleotides in this read chunk
			for (n,q) in zip(cSeq,cQual)
				nuc = base_index(n)
				nuc != 0 && inc!(nqa,refInd,refPos,nuc,q)
				refPos += 1
			end
		end
	end
end

# strands is one of: :both, :forward, :reverse
function nucqualitycount(bamFile::BamFile; strands=:both, mappingQualityThreshold::Int=30, ignoreChimeric::Bool=true)
	@assert strands∈[:both, :forward, :reverse] "Strand must be one of :both, :forward, :reverse."

	seqs = sequences(bamFile)
	seqLengths = Int[s[2] for s in seqs]

	flagsUnset = BAM_FUNMAP | BAM_FSECONDARY | BAM_FDUP
	flagsSet = 0

	strands == :forward && (flagsUnset |= BAM_FREVERSE)
	strands == :reverse && (flagsSet   |= BAM_FREVERSE)

	nqa = nucqualityaccumulator(seqLengths)
	processbam!(bamFile,nqa, flagsSet=flagsSet, flagsUnset=flagsUnset,
		        mappingQualityThreshold=mappingQualityThreshold,
		        ignoreChimeric=ignoreChimeric)
	nqa
end


text:
"""
	2D
"""

function integralArray( img::AbstractArray{<:Real,2}; type=Float32, fun=(x)->(x) )
	return integralArray!( img, zeros( type, size(img) .+ 1 ), fun=fun )
end

function integralArray!( img::AbstractArray{<:Real,2}, intArr::AbstractArray{<:AbstractFloat,2}; fun=(x)->(x) )

	@inbounds for c in 1:size(img,2), r in 1:size(img,1)
		intArr[1+r,1+c] = fun(img[r,c]) + intArr[1+r-1,1+c] + intArr[1+r,1+c-1] - intArr[1+r-1,1+c-1]
	end
	return intArr
end

function integralArea( intArr::AbstractArray{<:Real,2}, TL, BR )
	TL   = TL .+ 1;
	BR   = BR .+ 1;
	area = intArr[BR[1],BR[2]] - intArr[BR[1],TL[2]] - intArr[TL[1],BR[2]] + intArr[TL[1],TL[2]]
end

# integral array for ZNCC, copied from quickPIV

function integralArrayZNCC( padData::AbstractArray{Complex{T},N} ) where {T<:AbstractFloat,N}
	return integralArrayZNCC!( padData, zeros( T, size(padData).+1 ) );
end

function integralArrayZNCC!( padData::AbstractArray{Complex{T},2}, intArr::Array{T,2} ) where {T<:AbstractFloat}

	h, w = size(intArr) .- 1

	@inbounds for c in 2:w+1, r in 2:h+1
		intArr[r,c] = (padData[r-1,c-1].re)^2 + intArr[r-1,c] + intArr[r,c-1] - intArr[r-1,c-1]
	end

	return intArr
end

function integralArrayZNCC!( padData::AbstractArray{Complex{T},3}, intArr::Array{T,3} ) where {T<:AbstractFloat}

	h, w, d = size(intArr) .- 1;
	@inbounds for z in 2:d+1, c in 2:w+1
		tmp = 0.0;
		for r in 2:h+1
			val2          = real( padData[r-1,c-1,z-1] )^2
			intArr[r,c,z] = val2 + intArr[r,c-1,z] + intArr[r,c,z-1] - intArr[r,c-1,z-1] + tmp;
			tmp          += val2;
	end end

	return intArr
end

# what did IIM mean...?

function integralArrayIIM!( img::AbstractArray{<:Real,2}, intArr::AbstractArray{<:AbstractFloat,3}, z )

	@inbounds for c in 1:size(img,2), r in 1:size(img,1)
		intArr[1+r,1+c,z] = img[r,c] + intArr[1+r-1,1+c,z] + intArr[1+r,1+c-1,z] - intArr[1+r-1,1+c-1,z]
	end
	return intArr
end

function integralAreaIIM( intArr::AbstractArray{<:Real,3}, z,TL, BR )
	TL   = TL .+ 1;
	BR   = BR .+ 1;
	area = intArr[BR[1],BR[2],z] - intArr[BR[1],TL[2],z] - intArr[TL[1],BR[2],z] + intArr[TL[1],TL[2],z]
end

# Considering offset of the img within the integral array

function integralArray_off( img::AbstractArray{<:Real,2}, off=(0,0); type=Float32, fun=(x)->(x) )
	return integralArray_off!( img, zeros( type, size(img) .+ 1 ), off, fun=fun )
end

function integralArray_off!( img::AbstractArray{<:Real,2}, intArr::AbstractArray{<:AbstractFloat,2}, off=(0,0); fun=(x)->(x))

	@inbounds for c in 1:size(img,2), r in 1:size(img,1)
		intArr[1+r+off[1],1+c+off[2]] = fun(img[r,c]) + intArr[1+r-1+off[1],1+c+off[2]] + intArr[1+r+off[1],1+c-1+off[2]] - intArr[1+r-1+off[1],1+c-1+off[2]]
	end
	return intArr
end

# Specific function to compute mean maxima and minima, which requires padded integral arrays

function integralArray_pad( img::AbstractArray{<:Real,2}, pad=(0,0); type=Float32, fun=(x)->(x) )
	return integralArray_pad!( img, zeros( type, size(img) .+ 1 .+ 2 .* pad ), pad, fun=fun )
end

function integralArray_pad!( img::AbstractArray{<:Real,2}, intArr::AbstractArray{<:AbstractFloat,2}, pad; fun=(x)->(x) )

	for c in 1:size(img,2)      # c    = column in input image
		c_ia = 1+c+pad[2];      # c_ia = offset column index in integral array
		for r in 1:size(img,1)  # r    = row in input image
			r_ia = 1+r+pad[1];  # r_ia = offset row index in integral array
			intArr[r_ia,c_ia] = fun(img[r,c]) + intArr[r_ia-1,c_ia] + intArr[r_ia,c_ia-1]- intArr[r_ia-1,c_ia-1]
	end end

	lastr = 1+size(img,1)+pad[1];
	lastc = 1+size(img,2)+pad[2];

	for c in 1:lastc
		val = intArr[lastr,c]
		for r in lastr+1:size(intArr,1)
			intArr[r,c] = val
	end end

	for c in lastc+1:size(intArr,2)
		intArr[:,c] .= intArr[:,c-1]
	end

	return intArr
end

function integralArray_pad( vol::AbstractArray{<:Real,3}, pad=(0,0,0); type=Float32, fun=(x)->(x) )
	return integralArray_pad!( vol, zeros( type, size(vol) .+ 1 .+ 2 .*pad ), pad, fun=fun )
end

function integralArray_pad!( vol::AbstractArray{<:Real,3}, intArr::AbstractArray{<:AbstractFloat,3}, pad; fun=(x)->(x))

	for z in 1:size(vol,3)         # z    = depth index in input volume
		z_ia = z+pad[3]+1;         # z_ia = offset depth index in integral array
		for c in 1:size(vol,2)     # c    = columns in input volume
			c_ia = c+pad[2]+1;     # c_ia = offset column index in integral array
			reuse = 0.0
			for r in 1:size(vol,1) # r    = row in input volume
				r_ia = r+pad[1]+1; # r_ia = offset row index in integral array
				tmp = fun(vol[r,c,z])
				intArr[r_ia,c_ia,z_ia] = tmp + intArr[r_ia,c_ia-1,z_ia] + intArr[r_ia,c_ia,z_ia-1] - intArr[r_ia,c_ia-1,z_ia-1] + reuse;
				reuse += tmp
	end end end

	lastr = 1+size(vol,1)+pad[1];
	lastc = 1+size(vol,2)+pad[2];
	lastz = 1+size(vol,3)+pad[3];

	for z in 1:lastz
		for c in 1:lastc
			val = intArr[lastr,c,z]
			for r in lastr+1:size(intArr,1)
				intArr[r,c,z] = val
		end end

		for c in lastc+1:size(intArr,2)
			intArr[:,c,z] .= intArr[:,c-1,z]
		end
	end

	for z in lastz+1:size(intArr,3)
		intArr[:,:,z] .= intArr[:,:,z-1]
	end

	return intArr
end

# Reversed variant, which is used during normalized cross-correlation
function integralArrayRev!( img::AbstractArray{<:Real,2}, intArr::AbstractArray{<:AbstractFloat,2}, fun=(x)->(x) )
	h, w = size(img)
	@inbounds for c in 2:w+1, r in 2:h+1
		intArr[r,c] = fun(img[h-(r-2),w-(c-2)]) + intArr[r-1,c] + intArr[r,c-1] - intArr[r-1,c-1]
	end
	return intArr
end

function integralArrayRev( img::AbstractArray{<:Real,2}; type=Float32, fun=(x)->(x) )

	return integralArrayRev!( img, zeros( type, size(img) .+ 1 ), fun )
end

"""
	3D: optimized integral array by reusing previous operations
"""
# Below is the unoptimized code, in case I want to check again:
#    for z in 2:d+1, c in 2:w+1, r in 2:h+1
#        arr[r,c,z] = fun(vol[r-1,c-1,z-1]) + arr[r-1,c,z] + arr[r,c-1,z] + arr[r,c,z-1] -
#                     arr[r,c-1,z-1] - arr[r-1,c,z-1] - arr[r-1,c-1,z] + arr[r-1,c-1,z-1]
#    end

function integralArray!( vol::AbstractArray{<:Real,3}, intArr::AbstractArray{<:AbstractFloat,3}, fun=(x)->(x) )

	@inbounds for z in 2:size(vol,3)+1, c in 2:size(vol,2)+1
		reuse = 0.0
		for r in 2:size(vol,1)+1
			tmp = fun(vol[r-1,c-1,z-1])
			intArr[r,c,z] = tmp + intArr[r,c-1,z] + intArr[r,c,z-1] - intArr[r,c-1,z-1] + reuse;
			reuse += tmp
		end
	end
	return intArr
end

function integralArray( vol::AbstractArray{<:Real,3}; type=Float32, fun=(x)->(x) )

	return integralArray!( vol, zeros( type, size(vol) .+ 1 ), fun )
end

function integralArea( intArr::AbstractArray{<:Real,3}, TL, BR )
	TL = TL .+ 1;
	BR = BR .+ 1;
	area  = intArr[BR[1],BR[2],BR[3]] - intArr[TL[1],TL[2],TL[3]]
	area -= intArr[TL[1],BR[2],BR[3]] + intArr[BR[1],TL[2],BR[3]] + intArr[BR[1],BR[2],TL[3]]
	area += intArr[BR[1],TL[2],TL[3]] + intArr[TL[1],BR[2],TL[3]] + intArr[TL[1],TL[2],BR[3]]
    return area
end

# Reversed variant, which is used during normalized cross-correlation
function integralArrayRev!( vol::AbstractArray{<:Real,3}, intArr::AbstractArray{<:AbstractFloat,3}, fun=(x)->(x) )

	h, w, d = size(vol)
	@inbounds for z in 2:d+1, c in 2:w+1
		reuse = 0.0
		for r in 2:h+1
			tmp = fun(vol[h-(r-2),w-(c-2),d-(z-2)])
			intArr[r,c,z] = tmp + intArr[r,c-1,z] + intArr[r,c,z-1] - intArr[r,c-1,z-1] + reuse;
			reuse += tmp
		end
	end
	return intArr
end

function integralArrayRev( vol::AbstractArray{<:Real,3}; type=Float32, fun=(x)->(x) )

	return integralArrayRev!( vol, zeros( type, size(vol) .+ 1 ), fun )
end

"""
	Derived functions
"""
function scan4Simmetry( intArr, rad )

	score = zeros( length( rad+1:size(intArr,1)-rad), length( rad+1:size(intArr,2)-rad ) )

	cont = 1;
	for c in rad+1:size(intArr,2)-rad;     c1, c2 = c - rad, c + rad
		for r in rad+1:size(intArr,1)-rad; r1, r2 = r - rad, r + rad

			Q1  = integralArea( intArr, (r1,c ), ( r-1,c2-1) )
			Q2  = integralArea( intArr, (r1,c1), ( r-1,c -1) )
			Q3  = integralArea( intArr, (r ,c1), (r2-1,c -1) )
			Q4  = integralArea( intArr, (r ,c ), (r2-1,c2-1) )

			acc = Q1 + Q2 + Q3 + Q4;
			ps  = [ Q1, Q2, Q3, Q4 ] ./ acc

			score[ cont ] = 1/( 1 + sum( abs.( ps .- 0.25 ) ) )
			score[ cont ] = isnan( score[cont] ) ? 0.0 : score[cont];
			cont += 1
	end	end

	return score
end

function followCentroids( data, intArr, rad; iters=5 )

	# Goal: At any position, compute centroid, move to centroid
	# recompute centroid and so on until converging

	score = zeros( Int64, size(data) )

	xx = [ x for y in -rad:rad, x in -rad:rad ]
	yy = [ y for y in -rad:rad, x in -rad:rad ]

	h, w = size( data );
	for c in 1:w;
		for r in 1:h;

			pos = r, c;
			#println( pos )
			# Contemplating the case where the position is closer than "rad"
			# to the edges of the images. As we iterate and update the position
			# of the centroid, it becomes necessary. Here is not needed though.
			left   = rad + min(0, pos[1]-1-rad)
			up     = rad + min(0, pos[2]-1-rad)
			right  = min(  h - pos[ 1 ], rad  )
			down   = min(  w - pos[ 2 ], rad  )

			TL     = pos .- ( left , up   ) .- 1
			BR     = pos .+ ( right, down )

			acc    = integralArea( intArr, TL, BR )
			( acc == 0 ) && ( continue; )

			subimg = pos[1]-left:pos[1]+right, pos[2]-up:pos[2]+down
			subxy  =  rad+1-left:rad+1+right ,  rad+1-up:rad+1+down

			xmean  = sum( (data[ subimg... ] ./ acc) .* xx[ subxy... ] )
			ymean  = sum( (data[ subimg... ] ./ acc) .* yy[ subxy... ] )

			for i in 1:iters

				offset = round.( Int64, (ymean,xmean) )

				if offset[1] == 0 && offset[2] == 0
					break
				end

				pos    = pos .+ offset
				#println( "\t", pos )
				left   = rad + min(0, pos[1]-1-rad)
				up     = rad + min(0, pos[2]-1-rad)
				right  = min(  h - pos[ 1 ], rad  )
				down   = min(  w - pos[ 2 ], rad  )

				TL     = pos .- ( left , up   ) .- 1
				BR     = pos .+ ( right, down )

				acc    = integralArea( intArr, TL, BR )
				( acc == 0 ) && ( break; )

				subimg = pos[1]-left:pos[1]+right, pos[2]-up:pos[2]+down
				subxy  =  rad+1-left:rad+1+right ,  rad+1-up:rad+1+down

				xmean  = sum( (data[ subimg... ] ./ acc) .* xx[ subxy... ] )
				ymean  = sum( (data[ subimg... ] ./ acc) .* yy[ subxy... ] )
			end

			pos = pos .+ round.( Int64, (ymean,xmean) )

			score[ pos... ] += 1;
	end	end

	return score
end

function integralAreaZNCC( sumS::Array{T,2}, sumS2::Array{T,2},
                           row, col, ss, si ) where {T<:Real}
    pad = 1;
    row = row + pad;
    col = col + pad;
    r0, r1 = max( 1, row - si[1] ), min( ss[1] + pad, row );
    c0, c1 = max( 1, col - si[2] ), min( ss[2] + pad, col );

    opS  = 0.0
    opS += sumS[ r1, c1 ];
    opS -= sumS[ r1, c0 ];
    opS -= sumS[ r0, c1 ];
    opS += sumS[ r0, c0 ];

    opS2  = 0.0
    opS2 += sumS2[ r1, c1 ];
    opS2 -= sumS2[ r1, c0 ];
    opS2 -= sumS2[ r0, c1 ];
    opS2 += sumS2[ r0, c0 ];

    return opS, opS2
end

max_stars_repo_name: BrucePBon/ImageAnalysis

max_issues_repo_name: BrucePBon/ImageAnalysis

max_forks_repo_name: BrucePBon/ImageAnalysis


text:
import Base.isempty
using LinearAlgebra
export fgmres

function fgmres(A::SparseMatrixCSC{T1,Int},b::Array{T2,1},restrt::Int; kwargs...) where {T1,T2}
	Ax = zeros(promote_type(T1,T2),size(A,1))
	return fgmres(x -> mul!(Ax,A,x,1.0,0.0),b,restrt;kwargs...)
end

fgmres(A,b::Vector,restrt;kwargs...) = fgmres(x -> A*x ,b,restrt;kwargs...)

mutable struct FGMRESmem
	V    			::Array
	Z				::Array
end

function getFGMRESmem(n::Int,flexible::Bool,T::Type,k::Int,nrhs::Int=1)
	if flexible
		return FGMRESmem(zeros(T,n,k*nrhs),zeros(T,n,k*nrhs));
	else
		return FGMRESmem(zeros(T,n,k*nrhs),zeros(T,0));
	end
end

function getEmptyFGMRESmem()
	return FGMRESmem(zeros(0),zeros(0));
end

function isempty(mem::FGMRESmem)
	return size(mem.V,1)==0;
end

function checkMemorySize(mem::FGMRESmem,n::Int,k::Int,TYPE::Type,flexible::Bool,nrhs::Int=1)
	if isempty(mem)
		# warn("Allocating memory in FGMRES")
		mem = getFGMRESmem(n,flexible,TYPE,k,nrhs);
		return mem
	else
		if size(mem.V,2)!= k*nrhs
			error("FGMRES: size of Krylov subspace is different than inner*nrhs");
		end
	end
end

"""
x,flag,err,iter,resvec = fgmres(A,b,restrt,tol=1e-2,maxIter=100,M=1,x=[],out=0)

(flexible) Generalized Minimal residual ( (F)GMRES(m) ) method with restarts applied to A*x = b.

This is the "right preconditioning" version of GMRES: A*M^{-1}Mx = b. See gmres.jl for "left preconditioning".

Input:

A       - function computing A*x
b       - right hand side vector
restrt  - number of iterations between restarts
tol     - error tolerance
maxIter - maximum number of iterations
M       - preconditioner, function computing M\\x
x       - starting guess
out     - flag for output (0 : only errors, 1 : final status, 2: error at each iteration)

Output:

x       - approximate solution
flag    - exit flag (  0 : desired tolerance achieved,
	-1 : maxIter reached without converging
	-9 : right hand side was zero)
	err     - norm of relative residual, i.e., norm(A*x-b)/norm(b)
	iter    - number of iterations
	resvec  - norm of relative residual at each iteration

	preconditioner M(r) must return a copy of a vector. Cannot reuse memory of r.
	"""
function fgmres(A::Function,b::Vector,restrt::Int; tol::Real=1e-2,maxIter::Int=100,M::Function=t->copy(t),x::Vector=[],out::Int=0,storeInterm::Bool=false,flexible::Bool=false,mem::FGMRESmem =  getEmptyFGMRESmem())
  # initialization
  n  = length(b)
  TYPE = eltype(b)
  mem = checkMemorySize(mem,n,restrt,TYPE,flexible);

  if norm(b)==0.0
	return zeros(eltype(b),n),-9,0.0,0,[0.0];
  end

  if Base.isempty(x)
    x = zeros(eltype(b),n)
	r = copy(b);
  elseif norm(x) < eps(real(TYPE))
	r = copy(b);
  else
	r = copy(b);
	r.-=A(x);
  end

  if eltype(b) <: Complex
       x = complex(x)
   end

  if storeInterm
    X = zeros(eltype(b),n,maxIter)	# allocate space for intermediates
  end

  betta = norm(r);
  rnorm0 = norm(b);

  if rnorm0 == 0.0; rnorm0 = 1.0; end

  err = norm( r )/rnorm0
  if err < tol; return x, err; end

  constOne = one(TYPE);
  constZero = zero(TYPE);

  restrt = min(restrt,n-1)
  Z = mem.Z;
  V = mem.V;
  H = zeros(TYPE,restrt+1,restrt)
  xi = zeros(TYPE,restrt+1);
  t = zeros(TYPE,restrt);

  resvec = zeros(restrt*maxIter)
  if out==2
    println(@sprintf("=== fgmres ===\n%4s\t%7s\n","iter","relres"))
  end

  flag = -1

	counter = 0
	w = zeros(TYPE,0);
	iter = 0
	while iter < maxIter

		iter+=1;
		xi[1] = betta;
		BLAS.scal!(n,(1/betta)*constOne,r,1); # w = w./betta
		H[:] .= 0.0;
		t[:] .= 0.0;

		if out==2;; print(@sprintf("%3d\t", iter));end

		for j = 1:restrt

			if j==1
				V[:,j] = r;  # no memory problem with this line....
				z = M(r)
			else
				V[:,j] = w;  # no memory problem with this line....
				z = M(w)
			end
			if flexible
				Z[:,j] = z;
			end
			w = A(z); # w = A'*z;

			counter += 1;

			# ## modified Gram Schmidt:
			# for i=1:j
				# H[i,j] = dot(vec(V[:,i]),w);
				# w = w - H[i,j]*vec(V[:,i]);
			# end

			# Gram Schmidt (much faster than MGS even though zeros are multiplied, does relatively well):
			BLAS.gemv!('C', constOne, V, w,constZero,t);# t = V'*w;
			t[j+1:end] .= 0.0;
			H[1:restrt,j] = t;
			BLAS.gemv!('N', -constOne, V, t,constOne,w);# w = w - V*t;

			betta = norm(w);
			H[j+1,j] = betta;

			BLAS.scal!(n,(1/betta)*constOne,w,1); # w = w*(1/betta);

			# the following 2 lines are equivalent to the 2 next
			# y = H[1:j+1,1:j]\xi[1:j+1];
			# err = norm(H[1:j+1,1:j]*y - xi[1:j+1])/rnorm0
			Q = qr(H[1:j+1,1:j]).Q;
			err = abs(Q[1,end]*xi[1])/rnorm0

			if out==2 print(@sprintf("%1.1e ", err)); end
			resvec[counter] = err;
			if err <= tol
				if flexible
					Z[:,j+1:end] .= 0.0;
				else
					V[:,j+1:end] .= 0.0;
				end
				if out==2; print("\n"); end
				flag = 0; break
			end
		end # end for j to restrt

		y = pinv(H)*xi;

		if flexible
			BLAS.gemv!('N', constOne, Z, y,constZero,w); # w = Z*y  #This is the correction that corresponds to the residual.
		else
			BLAS.gemv!('N', constOne, V, y, constZero, w); #w = V*y
			z = M(w);
			w[:] = z;
		end

		x .+= w;

		if storeInterm; X[:,iter] = x; end

		if out==2; print("\n"); end
		if err <= tol
			flag = 0;
			break
		end
		if iter < maxIter
			r = copy(b);
			r.-=A(x);
			betta = norm(r)
		end
		# if out==2; print(@sprintf("\t %1.1e\n", err)); end
	end #for iter to maxiter
	if out>=0
		if flag==-1
			println(@sprintf("fgmres iterated maxIter (=%d) outer iterations without achieving the desired tolerance. Acheived: %1.2e",maxIter,err))
		elseif flag==0 && out>=1
			println(@sprintf("fgmres achieved desired tolerance at inner iteration %d. Residual norm is %1.2e.",counter,resvec[counter]))
		end
	end

	if storeInterm
		return X[:,1:iter],flag,resvec[counter],iter,resvec[1:counter]
	else
		return x,flag,resvec[counter],iter,resvec[1:counter]
	end
end #fgmres


text:
#=
methods for handling backtest results of strategy objects
=#

struct Backtest
    trades::Dict{String,TS}
    backtest::Dict{String,TS{Float64}}
    optimization::Matrix{Float64}
    function Backtest(trades::Dict{String,TS}=Dict{String,TS}(),
                      backtest::Dict{String,TS{Float64}}=Dict{String,TS{Float64}}(),
                      optimization::Matrix{Float64}=Matrix{Float64}(undef,0,0))
        return new(trades, backtest, optimization)
    end
end

function cum_pnl(results::Backtest)::Float64
    result = 0.0
    @inbounds for val in values(results.backtest)
        pnl::Vector = val[:PNL].values[:]
        result += sum(pnl)
    end
    return result
end

function cum_pnl(results::Dict{String, TS{Float64}})::Float64
    result = 0.0
    @inbounds for val in values(results)
        pnl::Vector = val[:PNL].values[:]
        result += sum(pnl)
    end
    return result
end

max_stars_repo_name: lithiferous/Strathreads

max_issues_repo_name: lithiferous/Strathreads

max_forks_repo_name: lithiferous/Strathreads


text:
################################################################################
#                              OPTIMIZER MODEL BASICS
################################################################################
"""
    optimizer_model(model::InfiniteModel)::JuMP.Model

Return the JuMP model stored in `model` that is used to solve it.

**Example**
```julia-repl
julia> opt_model = optimizer_model(model)
A JuMP Model
Feasibility problem with:
Variables: 0
Model mode: AUTOMATIC
CachingOptimizer state: NO_OPTIMIZER
Solver name: No optimizer attached.
```
"""
optimizer_model(model::InfiniteModel)::JuMP.Model = model.optimizer_model

"""
    JuMP.bridge_constraints(model::InfiniteModel)::Bool

Extend `JuMP.bridge_constraints` to return if an infinite model `model`
has an optimizer model where the optimizer is set and unsupported constraints
are automatically bridged to equivalent supported constraints when an
appropriate transformation is available.

**Example**
```julia-repl
julia> bridge_constraints(model)
false
```
"""
function JuMP.bridge_constraints(model::InfiniteModel)::Bool
    return JuMP.bridge_constraints(optimizer_model(model))
end

"""
    JuMP.add_bridge(model::InfiniteModel,
                    BridgeType::Type{<:MOI.Bridges.AbstractBridge})

Extend `JuMP.add_bridge` to add `BridgeType` to the list of bridges that can
be used by the optimizer model to transform unsupported constraints into an
equivalent formulation using only constraints supported by the optimizer.
"""
function JuMP.add_bridge(model::InfiniteModel,
                         BridgeType::Type{<:MOI.Bridges.AbstractBridge})
    JuMP.add_bridge(optimizer_model(model), BridgeType)
    return
end

"""
    optimizer_model_ready(model::InfiniteModel)::Bool

Return `Bool` if the optimizer model is up to date with `model`.

**Example**
```julia-repl
julia> optimizer_model_ready(model)
false
```
"""
optimizer_model_ready(model::InfiniteModel)::Bool = model.ready_to_optimize

"""
    set_optimizer_model_ready(model::InfiniteModel, status::Bool)

Set the status of the optimizer model to whether it is up to date or not. Note
is more intended as an internal function, but is useful for extensions.

**Example**
```julia-repl
julia> set_optimizer_model_ready(model, true)

julia> optimizer_model_ready(model)
true
```
"""
function set_optimizer_model_ready(model::InfiniteModel, status::Bool)
     model.ready_to_optimize = status
     return
end

"""
    add_infinite_model_optimizer(opt_model::JuMP.Model, inf_model::InfiniteModel)

Parse the current optimizer and its attributes associated with `inf_model` and load
them into `opt_model`. This is intended to be used as an internal method
for [`set_optimizer_model`](@ref).
"""
function add_infinite_model_optimizer(opt_model::JuMP.Model,
                                      inf_model::InfiniteModel)
    if !isa(inf_model.optimizer_constructor, Nothing)
        bridge_constrs = JuMP.bridge_constraints(inf_model)
        JuMP.set_optimizer(opt_model, inf_model.optimizer_constructor,
                           add_bridges =  bridge_constrs)
    end
    # parse the attributes (this is a hacky workaround)
    for (attr, val) in JuMP.backend(inf_model).model_cache.optattr
        MOI.set(opt_model, attr, val)
    end
    return
end

"""
    set_optimizer_model(inf_model::InfiniteModel, opt_model::JuMP.Model;
                        inherit_optimizer::Bool = true)

Specify the JuMP model that is used to solve `inf_model`. This is intended for
internal use and extensions. Note that `opt_model` should contain extension
data to allow it to map to `inf_model` in a manner similar to
[`TranscriptionModel`](@ref). `inherit_optimizer` indicates whether
[`add_infinite_model_optimizer`](@ref) should be invoked on the new optimizer
mode to inherit the optimizer constuctor and attributes currently stored in
`inf_model`.

**Example**
```julia-repl
julia> set_optimizer_model(model, TranscriptionModel())

julia> optimizer_model(model)
A JuMP Model
Feasibility problem with:
Variables: 0
Model mode: AUTOMATIC
CachingOptimizer state: NO_OPTIMIZER
Solver name: No optimizer attached.
```
"""
function set_optimizer_model(
    inf_model::InfiniteModel,
    opt_model::JuMP.Model;
    inherit_optimizer::Bool = true
    )
    if inherit_optimizer
        add_infinite_model_optimizer(opt_model, inf_model)
    end
    inf_model.optimizer_model = opt_model
    set_optimizer_model_ready(inf_model, false)
    return
end

"""
    optimizer_model_key(model::JuMP.Model)::Any

Return the extension key used in the optimizer model `model`. Errors if
`model.ext` contains more than one key. This is intended for internal
use and extensions. For extensions this is used to dispatch to the appropriate
optmizer model functions such as extensions to [`build_optimizer_model!`](@ref).
This is intended as an internal method. See
[`optimizer_model_key`](@ref optimizer_model_key(::InfiniteModel))
for the public method
"""
function optimizer_model_key(model::JuMP.Model)
    if length(model.ext) != 1
        error("Optimizer models should have 1 and only 1 extension key of the " *
              "form `Model.ext[:my_ext_key] = MyExtData`.")
    end
    return first(keys(model.ext))
end

"""
    optimizer_model_key(model::InfiniteModel)::Any

Return the extension key used in the optimizer model of `model`. Errors if
`optimizer_model.ext` contains more than one key. This is intended for internal
use and extensions. For extensions this is used to dispatch to the appropriate
optmizer model functions such as extensions to [`build_optimizer_model!`](@ref).

**Example**
```julia-repl
julia> optimizer_model_key(model)
:TransData
```
"""
function optimizer_model_key(model::InfiniteModel)::Any
    return optimizer_model_key(optimizer_model(model))
end

################################################################################
#                         OPTIMIZER METHOD EXTENSIONS
################################################################################
"""
    JuMP.set_optimizer(model::InfiniteModel,
                       [optimizer_constructor;
                       add_bridges::Bool = true])

Extend `JuMP.set_optimizer` to set optimizer of infinite models.
Specifically, the optimizer of the optimizer model is modified.

**Example**
```julia-repl
julia> set_optimizer(model, Clp.Optimizer)

julia> optimizer_model(model)
A JuMP Model
Feasibility problem with:
Variables: 0
Model mode: AUTOMATIC
CachingOptimizer state: EMPTY_OPTIMIZER
Solver name: SolverName() attribute not implemented by the optimizer.
```
"""
function JuMP.set_optimizer(
    model::InfiniteModel,
    optimizer_constructor;
    add_bridges::Bool = true
    )
    JuMP.set_optimizer(optimizer_model(model), optimizer_constructor,
                       add_bridges = add_bridges)
    _set_optimizer_constructor(model, optimizer_constructor)
    return
end

"""
    JuMP.set_silent(model::InfiniteModel)

Extend `JuMP.set_silent` for infinite models to take precedence over any other
attribute controlling verbosity and requires the solver to produce no output.

**Example**
```julia-repl
julia> set_silent(model)
true
```
"""
function JuMP.set_silent(model::InfiniteModel)
    return JuMP.set_silent(optimizer_model(model))
end

"""
    JuMP.unset_silent(model::InfiniteModel)

Extend `JuMP.unset_silent` for infinite models to neutralize the effect of the
`set_silent` function and let the solver attributes control the verbosity.

**Example**
```julia-repl
julia> unset_silent(model)
false
```
"""
function JuMP.unset_silent(model::InfiniteModel)
    return JuMP.unset_silent(optimizer_model(model))
end

"""
    JuMP.set_time_limit_sec(model::InfiniteModel, limit)

Extend `set_time_limit_sec` to set the time limit (in seconds) of the solver.
Can be unset using `unset_time_limit_sec` or with `limit` set to `nothing`.

**Example**
```julia-repl
julia> set_time_limit_sec(model, 100)
100
```
"""
function JuMP.set_time_limit_sec(model::InfiniteModel, limit)
    return JuMP.set_time_limit_sec(optimizer_model(model), limit)
end

"""
    JuMP.unset_time_limit_sec(model::InfiniteModel)

Extend `unset_time_limit_sec` to unset the time limit of the solver. Can be set
using `set_time_limit_sec`.

**Example**
```julia-repl
julia> unset_time_limit_sec(model)
```
"""
function JuMP.unset_time_limit_sec(model::InfiniteModel)
    return JuMP.unset_time_limit_sec(optimizer_model(model))
end

"""
    JuMP.time_limit_sec(model::InfiniteModel)

Extend `time_limit_sec` to get the time limit (in seconds) of the solve used by
the optimizer model (`nothing` if unset). Can be set using `set_time_limit_sec`.

**Example**
```julia-repl
julia> time_limit_sec(model)
100
```
"""
function JuMP.time_limit_sec(model::InfiniteModel)
    return JuMP.time_limit_sec(optimizer_model(model))
end

"""
    JuMP.set_optimizer_attribute(model::InfiniteModel, name::String, value)

Extend `set_optimizer_attribute` to specify a solver-specific attribute
identified by `name` to `value`.

**Example**
```julia-repl
julia> set_optimizer_attribute(model, "SolverSpecificAttributeName", true)
true
```
"""
function JuMP.set_optimizer_attribute(model::InfiniteModel, name::String, value)
    return JuMP.set_optimizer_attribute(optimizer_model(model), name, value)
end

"""
    JuMP.set_optimizer_attribute(model::InfiniteModel,
                                 attr::MOI.AbstractOptimizerAttribute,
                                 value)

Extend `set_optimizer_attribute` to set the solver-specific attribute `attr` in
`model` to `value`.

**Example**
```julia-repl
julia> set_optimizer_attribute(model, MOI.Silent(), true)
true
```
"""
function JuMP.set_optimizer_attribute(
    model::InfiniteModel,
    attr::MOI.AbstractOptimizerAttribute,
    value
    )
    return MOI.set(optimizer_model(model), attr, value)
end

"""
    JuMP.set_optimizer_attributes(model::InfiniteModel, pairs::Pair...)

Extend `set_optimizer_attributes` to set multiple solver attributes given a
list of `attribute => value` pairs. Calls
`set_optimizer_attribute(model, attribute, value)` for each pair.

**Example**
```julia-repl
julia> model = Model(Ipopt.Optimizer);

julia> set_optimizer_attributes(model, "tol" => 1e-4, "max_iter" => 100)
```
is equivalent to:
```julia-repl
julia> set_optimizer_attribute(model, "tol", 1e-4);

julia> set_optimizer_attribute(model, "max_iter", 100);
```
"""
function JuMP.set_optimizer_attributes(model::InfiniteModel, pairs::Pair...)
    for (name, value) in pairs
        JuMP.set_optimizer_attribute(model, name, value)
    end
    return
end

"""
    JuMP.get_optimizer_attribute(model::InfiniteModel, name::String)

Extend `get_optimizer_attribute` to return the value associated with the
solver-specific attribute named `name`.

**Example**
```julia-repl
julia> get_optimizer_attribute(model, "tol")
0.0001
````
"""
function JuMP.get_optimizer_attribute(model::InfiniteModel, name::String)
    return JuMP.get_optimizer_attribute(optimizer_model(model), name)
end

"""
    JuMP.get_optimizer_attribute(model::InfiniteModel,
                                 attr::MOI.AbstractOptimizerAttribute)

Extend `get_optimizer_attribute` to return the value of the solver-specific
attribute `attr` in `model`.

**Example**
```julia-repl
julia> get_optimizer_attribute(model, MOI.Silent())
true
````
"""
function JuMP.get_optimizer_attribute(
    model::InfiniteModel,
    attr::MOI.AbstractOptimizerAttribute
    )
    return MOI.get(optimizer_model(model), attr)
end

"""
    JuMP.solver_name(model::InfiniteModel)

Extend `solver_name` to return the name of the solver being used if there is an
optimizer selected and it has a name attribute. Otherwise, an error is thrown.

**Example**
```julia-repl
julia> solver_name(model)
"Gurobi"
```
"""
function JuMP.solver_name(model::InfiniteModel)
    return JuMP.solver_name(optimizer_model(model))
end

"""
    JuMP.backend(model::InfiniteModel)

Extend `backend` to return the `MathOptInterface` backend associated with the
optimizer model. Note this will be empty if the optimizer model has not been
build yet.

**Example**
```julia-repl
julia> moi_model = backend(model);
```
"""
function JuMP.backend(model::InfiniteModel)
    return JuMP.backend(optimizer_model(model))
end

"""
    JuMP.mode(model::InfiniteModel)

Extend `mode` to return the `MathOptInterface` mode the optimizer model is in.

**Example**
```julia-repl
julia> mode(model)
AUTOMATIC::ModelMode = 0
```
"""
function JuMP.mode(model::InfiniteModel)
    return JuMP.mode(optimizer_model(model))
end

"""
    JuMP.result_count(model::InfiniteModel)

Extend `result_count` to return the number of results available to query after a
call to `optimize!`.

**Example**
```julia-repla
julia> result_count(model)
1
```
"""
function JuMP.result_count(model::InfiniteModel)::Int
    return MOI.get(optimizer_model(model), MOI.ResultCount())
end

################################################################################
#                       OPTIMIZER MODEL BUILD METHODS
################################################################################
"""
    build_optimizer_model!(model::InfiniteModel, key::Val{ext_key_name};
                           [kwargs...])

Build the optimizer model stored in `model` such that it can be
treated as a normal JuMP model, where the `Model.ext` field contains a key
that points to a datastructure that appropriately maps the data between the
two models. The key argument should be be typed to `Val{ext_key_name}`. This
should also use [`clear_optimizer_model_build!`](@ref) to empty the out the current
optimizer model. Ultimately, [`set_optimizer_model`](@ref) should be called
to insert the build optimizer model into `model` and [`set_optimizer_model_ready`](@ref)
should be used to update the optimizer model's status.
"""
function build_optimizer_model! end

"""
    clear_optimizer_model_build!(model::JuMP.Model)::JuMP.Model

Empty the optimizer model using appropriate calls of `Base.empty!`. This
effectively resets `model` except the optimizer, its attributes, and an an emptied
optimizer model data struct are maintained. This is intended as an internal
method for use by [`build_optimizer_model!`](@ref).
"""
function clear_optimizer_model_build!(model::JuMP.Model)::JuMP.Model
    MOI.empty!(model.moi_backend)
    empty!(model.shapes)
    model.nlp_data = nothing
    empty!(model.obj_dict)
    key = optimizer_model_key(model)
    data_type = typeof(model.ext[key])
    model.ext[key] = data_type() # assumes a constructor has been implemented.
    return model
end

"""
    clear_optimizer_model_build!(model::InfiniteModel)::JuMP.Model

Empty the optimizer model using appropriate calls of `Base.empty!`. This
effectively resets `model.optimizer_model` except the optimizer, its attributes,
and an an emptied optimizer model data struct are maintained. This is intended
as an internal method for use by [`build_optimizer_model!`](@ref).
"""
function clear_optimizer_model_build!(model::InfiniteModel)::JuMP.Model
    return clear_optimizer_model_build!(optimizer_model(model))
end

"""
    build_optimizer_model!(model::InfiniteModel; [kwargs...])

Build the optimizer model stored in `model` such that it can be
treated as a normal JuMP model. Specifically, translate the variables and
constraints stored in `model` into ones that are stored in the optimizer model
and can be solved. This is provided generally to accomodate extensions that use
custom optimizer model types in accordance with [`optimizer_model_key`](@ref).
However, it may be useful in certain applications when the user desires to
force a build without calling `optimize!`.
Extensions will need to implement their own version of the function
`build_optimizer_model!(model::InfiniteModel, key::Val{ext_key_name}; kwargs...)`.

**Example**
```julia-repl
julia> build_optimizer_model!(model)

julia> optimizer_model_ready(model)
true
```
"""
function build_optimizer_model!(model::InfiniteModel; kwargs...)
    if num_parameters(model, InfiniteParameter) == 0
        @warn("Finite models (i.e., `InfiniteModel`s with no infinite " *
              "parameters) should be modeled directly via a `Model` in JuMP.jl.")
    end
    key = optimizer_model_key(model)
    build_optimizer_model!(model, Val(key); kwargs...)
    return
end

################################################################################
#                  OPTIMIZER MODEL MAPPING METHODS (VARIABLES)
################################################################################
"""
    optimizer_model_variable(vref::GeneralVariableRef, key::Val{ext_key_name};
                             [kwargs...])

Return the reformulation variable(s) stored in the optimizer model that correspond
to `vref`. This needs to be defined for extensions that implement a custom
optimizer model type. Principally, this is accomplished by typed the `key`
argument to `Val{ext_key_name}`. Keyword arguments can be added as needed.
"""
function optimizer_model_variable end

# Fallback for unextended keys
function optimizer_model_variable(vref::GeneralVariableRef, key; kwargs...)
    error("`optimizer_model_variable` not implemented for optimizer model " *
          "key `$(typeof(key).parameters[1])`.")
end

"""
    optimizer_model_variable(vref::GeneralVariableRef;
                             [label::Type{<:AbstractSupportLabel} = PublicLabel,
                             ndarray::Bool = false,
                             kwargs...])

Return the reformulation variable(s) stored in the optimizer model that correspond
to `vref`. Also errors if no such variable can be found in
the optimizer model.

The keyword arugments `label` and `ndarray` are what `TranscriptionOpt` employ
and `kwargs` denote extra ones that user extensions may employ in accordance with
their implementation of [`optimizer_model_variable`](@ref). Errors if such an
extension has not been written.

By default only the variables associated with public supports are returned, the
full set can be accessed via `label = All`. Moreover, infinite variables are
returned as a list corresponding to their supports. However, a n-dimensional array
can be obtained via `ndarray = true` which is handy when the variable has multiple
infinite parameter dependencies. The corresponding supports are obtained via
`supports` using the same keyword arguments.

**Example**
```julia-repl
julia> optimizer_model_variable(x) # infinite variable
2-element Array{VariableRef,1}:
 x(support: 1)
 x(support: 2)

julia> optimizer_model_variable(z) # finite variable
z
```
"""
function optimizer_model_variable(vref::GeneralVariableRef; kwargs...)
    key = optimizer_model_key(JuMP.owner_model(vref))
    return optimizer_model_variable(vref, Val(key); kwargs...)
end

"""
    variable_supports(optimizer_model::JuMP.Model, vref,
                      key::Val{ext_key_name};
                      [kwargs...])::Vector

Return the supports associated with the mappings of `vref` in `optimizer_model`.
This dispatches off of `key` which permits optimizer model extensions. This
should throw an error if `vref` is not associated with the variable mappings
stored in `optimizer_model`. Keyword arguments can be added as needed. Note that
no extension is necessary for point or finite variables.
"""
function variable_supports end

# fallback for unextended keys
function variable_supports(optimizer_model::JuMP.Model, vref, key; kwargs...)
    error("`variable_supports` not implemented for optimizer model key " *
          "`$(typeof(key).parameters[1])` and/or variable type $(typeof(vref)).")
end

# FiniteRef
function variable_supports(optimizer_model::JuMP.Model, vref::FiniteRef,
                           key; kwargs...)
    return ()
end

"""
    supports(vref::DecisionVariableRef;
             [label::Type{<:AbstractSupportLabel} = PublicLabel,
             ndarray::Bool = false,
             kwargs...])

Return the supports associated with `vref` in the optimizer
model. Errors if [`InfiniteOpt.variable_supports`](@ref) has not been extended for the
optimizer model type or if `vref` is not be reformulated in the optimizer model.

The keyword arugments `label` and `ndarray` are what `TranscriptionOpt` employ
and `kwargs` denote extra ones that user extensions may employ in accordance with
their implementation of `variable_supports`. Errors if such an
extension has not been written.

By default only the public supports are returned, the
full set can be accessed via `label = All`. Moreover, the supports of infinite
variables are returned as a list. However, a n-dimensional array
can be obtained via `ndarray = true` which is handy when the variable has multiple
infinite parameter dependencies.

**Example**
```julia-repl
julia> supports(vref)
2-element Array{Tuple{Float64},1}:
 (0.0,)
 (1.0,)
```
"""
function supports(
    vref::Union{DecisionVariableRef, MeasureRef, ParameterFunctionRef};
    kwargs...
    )
    model = optimizer_model(JuMP.owner_model(vref))
    key = optimizer_model_key(JuMP.owner_model(vref))
    return variable_supports(model, vref, Val(key); kwargs...)
end

################################################################################
#                 OPTIMIZER MODEL MAPPING METHODS (EXPRESSIONS)
################################################################################
"""
    optimizer_model_expression(expr, key::Val{ext_key_name}; [kwargs...])

Return the reformulation expression(s) stored in the optimizer model that correspond
to `expr`. This needs to be defined for extensions that implement a custom
optimizer model type. Principally, this is accomplished by typed the `key`
argument to `Val{ext_key_name}`. Keyword arguments can be added as needed.
Note that if `expr` is a `GeneralVariableRef` this just dispatches to
`optimizer_model_variable`.
"""
function optimizer_model_expression end

# Fallback for unextended keys
function optimizer_model_expression(expr, key; kwargs...)
    error("`optimizer_model_expression` not defined for optimizer model " *
          "key `$(typeof(key).parameters[1])` and expression type " *
          "`$(typeof(expr))`.")
end

# Define for variable reference expressions
function optimizer_model_expression(expr::GeneralVariableRef, key; kwargs...)
    return optimizer_model_variable(expr, key; kwargs...)
end

"""
    optimizer_model_expression(expr::JuMP.AbstractJuMPScalar;
                               [label::Type{<:AbstractSupportLabel} = PublicLabel,
                               ndarray::Bool = false,
                               kwargs...])

Return the reformulation expression(s) stored in the optimizer model that correspond
to `expr`. Also errors if no such expression can be found in
the optimizer model (meaning one or more of the underlying variables have not
been transcribed).

The keyword arugments `label` and `ndarray` are what `TranscriptionOpt` employ
and `kwargs` denote extra ones that user extensions may employ in accordance with
their implementation of [`optimizer_model_expression`](@ref). Errors if such an
extension has not been written.

By default only the expressions associated with public supports are returned, the
full set can be accessed via `label = All`. Moreover, infinite expressions are
returned as a list corresponding to their supports. However, a n-dimensional array
can be obtained via `ndarray = true` which is handy when the expression has multiple
infinite parameter dependencies. The corresponding supports are obtained via
`supports` using the same keyword arguments.

**Example**
```julia-repl
julia> optimizer_model_expression(my_expr) # finite expression
x(support: 1) - y
```
"""
function optimizer_model_expression(expr::JuMP.AbstractJuMPScalar; kwargs...)
    model = _model_from_expr(expr)
    if isnothing(model)
        return zero(JuMP.AffExpr) + JuMP.constant(expr)
    else
        key = optimizer_model_key(model)
        return optimizer_model_expression(expr, Val(key); kwargs...)
    end
end

"""
    expression_supports(optimizer_model::JuMP.Model, expr,
                        key::Val{ext_key_name}; [kwargs...])

Return the supports associated with the mappings of `expr` in `optimizer_model`.
This dispatches off of `key` which permits optimizer model extensions. This
should throw an error if `expr` is not associated with the variable mappings
stored in `optimizer_model`. Keyword arguments can be added as needed. Note that
if `expr` is a `GeneralVariableRef` this just dispatches to `variable_supports`.
"""
function expression_supports end

# fallback for unextended keys
function expression_supports(optimizer_model::JuMP.Model, expr, key; kwargs...)
  error("`constraint_supports` not implemented for optimizer model key " *
        "`$(typeof(key).parameters[1])` and/or expressions of type " *
        "`$(typeof(expr))`.")
end

# Variable reference expressions
function expression_supports(model::JuMP.Model, vref::GeneralVariableRef, key;
                             kwargs...)
    return variable_supports(model, dispatch_variable_ref(vref), key; kwargs...)
end

"""
    supports(expr::JuMP.AbstractJuMPScalar;
             [label::Type{<:AbstractSupportLabel} = PublicLabel,
             ndarray::Bool = false,
             kwargs...])

Return the support associated with `expr`. Errors if `expr` is
not associated with the constraint mappings stored in `optimizer_model`.

The keyword arugments `label` and `ndarray` are what `TranscriptionOpt` employ
and `kwargs` denote extra ones that user extensions may employ in accordance with
their implementation of `expression_supports`. Errors if such an
extension has not been written.

By default only the public supports are returned, the
full set can be accessed via `label = All`. Moreover, the supports of infinite
expressions are returned as a list. However, a n-dimensional array
can be obtained via `ndarray = true` which is handy when the expression has multiple
infinite parameter dependencies.

**Example**
```julia-repl
julia> supports(cref)
2-element Array{Tuple{Float64},1}:
 (0.0,)
 (1.0,)
```
"""
function supports(expr::JuMP.AbstractJuMPScalar; kwargs...)
    model = _model_from_expr(expr)
    if isnothing(model)
        return ()
    else
        key = optimizer_model_key(model)
        opt_model = optimizer_model(model)
        return expression_supports(opt_model, expr, Val(key); kwargs...)
    end
end

################################################################################
#                OPTIMIZER MODEL MAPPING METHODS (CONSTRAINTS)
################################################################################
"""
    optimizer_model_constraint(cref::InfOptConstraintRef,
                               key::Val{ext_key_name}; [kwargs...])

Return the reformulation constraint(s) stored in the optimizer model that correspond
to `cref`. This needs to be defined for extensions that implement a custom
optimizer model type. Principally, this is accomplished by typed the `key`
argument to `Val{ext_key_name}`. Keyword arguments can be added as needed.
"""
function optimizer_model_constraint end

# Fallback for unextended keys
function optimizer_model_constraint(
    cref::InfOptConstraintRef,
    key;
    kwargs...
    )
    error("`optimizer_model_constraint` not implemented for optimizer model " *
          "key `$(typeof(key).parameters[1])`.")
end

"""
    optimizer_model_constraint(cref::InfOptConstraintRef;
                               [label::Type{<:AbstractSupportLabel} = PublicLabel,
                               ndarray::Bool = false,
                               kwargs...])

Return the reformulation constraint(s) stored in the optimizer model that correspond
to `cref`. Errors if no such constraint can be found in
the optimizer model.

The keyword arugments `label` and `ndarray` are what `TranscriptionOpt` employ
and `kwargs` denote extra ones that user extensions may employ in accordance with
their implementation of [`optimizer_model_constraint`](@ref). Errors if such an
extension has not been written.

By default only the constraints associated with public supports are returned, the
full set can be accessed via `label = All`. Moreover, infinite constraints are
returned as a list corresponding to their supports. However, a n-dimensional array
can be obtained via `ndarray = true` which is handy when the constraint has multiple
infinite parameter dependencies. The corresponding supports are obtained via
`supports` using the same keyword arguments.

**Example**
```julia-repl
julia> optimizer_model_constraint(c1) # finite constraint
c1 : x(support: 1) - y <= 3.0
```
"""
function optimizer_model_constraint(
    cref::InfOptConstraintRef;
    kwargs...
    )
    key = optimizer_model_key(JuMP.owner_model(cref))
    return optimizer_model_constraint(cref, Val(key); kwargs...)
end

"""
    constraint_supports(optimizer_model::JuMP.Model,
                        cref::InfOptConstraintRef,
                        key::Val{ext_key_name}; [kwargs...])

Return the supports associated with the mappings of `cref` in `optimizer_model`.
This dispatches off of `key` which permits optimizer model extensions. This
should throw an error if `cref` is not associated with the variable mappings
stored in `optimizer_model`. Keyword arguments can be added as needed.
"""
function constraint_supports end

# fallback for unextended keys
function constraint_supports(optimizer_model::JuMP.Model,
                             cref::InfOptConstraintRef,
                             key; kwargs...)
  error("`constraint_supports` not implemented for optimizer model key " *
        "`$(typeof(key).parameters[1])`.")
end

"""
    supports(cref::InfOptConstraintRef;
             [label::Type{<:AbstractSupportLabel} = PublicLabel,
             ndarray::Bool = false,
             kwargs...])

Return the support associated with `cref`. Errors if `cref` is
not associated with the constraint mappings stored in `optimizer_model`.

The keyword arugments `label` and `ndarray` are what `TranscriptionOpt` employ
and `kwargs` denote extra ones that user extensions may employ in accordance with
their implementation of `constraint_supports`. Errors if such an
extension has not been written.

By default only the public supports are returned, the
full set can be accessed via `label = All`. Moreover, the supports of infinite
constraints are returned as a list. However, a n-dimensional array
can be obtained via `ndarray = true` which is handy when the constraint has multiple
infinite parameter dependencies.

**Example**
```julia-repl
julia> supports(cref)
2-element Array{Tuple{Float64},1}:
 (0.0,)
 (1.0,)
```
"""
function supports(cref::InfOptConstraintRef; kwargs...)
    model = optimizer_model(JuMP.owner_model(cref))
    key = optimizer_model_key(JuMP.owner_model(cref))
    return constraint_supports(model, cref, Val(key); kwargs...)
end

################################################################################
#                             OPTIMIZATION METHODS
################################################################################
# Extend the solve error function
JuMP.solve(model::InfiniteModel) = JuMP.solve(optimizer_model(model))

"""
    JuMP.optimize!(model::InfiniteModel; [kwargs...])

Extend `JuMP.optimize!` to optimize infinite models using the internal
optimizer model. Will call [`build_optimizer_model!`](@ref) if the optimizer
model isn't up to date. The `kwargs` correspond to keyword arguments passed to
[`build_optimizer_model!`](@ref) if any are defined.

**Example**
```julia-repl
julia> optimize!(model)

julia> has_values(model)
true
```
"""
function JuMP.optimize!(model::InfiniteModel; kwargs...)
    if !optimizer_model_ready(model)
        build_optimizer_model!(model; kwargs...)
    end
    JuMP.optimize!(optimizer_model(model))
    return
end