text:
function Array2DRowRealMatrix()
    return Array2DRowRealMatrix(())
end

function Array2DRowRealMatrix(arg0::Vector{Vector{jdouble}})
    return Array2DRowRealMatrix((Vector{Vector{jdouble}},), arg0)
end

function Array2DRowRealMatrix(arg0::Vector{Vector{jdouble}}, arg1::jboolean)
    return Array2DRowRealMatrix((Vector{Vector{jdouble}}, jboolean), arg0, arg1)
end

function Array2DRowRealMatrix(arg0::Vector{jdouble})
    return Array2DRowRealMatrix((Vector{jdouble},), arg0)
end

function Array2DRowRealMatrix(arg0::jint, arg1::jint)
    return Array2DRowRealMatrix((jint, jint), arg0, arg1)
end

function add(obj::Array2DRowRealMatrix, arg0::Array2DRowRealMatrix)
    return jcall(obj, "add", Array2DRowRealMatrix, (Array2DRowRealMatrix,), arg0)
end

function add_to_entry(obj::Array2DRowRealMatrix, arg0::jint, arg1::jint, arg2::jdouble)
    return jcall(obj, "addToEntry", void, (jint, jint, jdouble), arg0, arg1, arg2)
end

function copy(obj::Array2DRowRealMatrix)
    return jcall(obj, "copy", RealMatrix, ())
end

function create_matrix(obj::Array2DRowRealMatrix, arg0::jint, arg1::jint)
    return jcall(obj, "createMatrix", RealMatrix, (jint, jint), arg0, arg1)
end

function get_column_dimension(obj::Array2DRowRealMatrix)
    return jcall(obj, "getColumnDimension", jint, ())
end

function get_data(obj::Array2DRowRealMatrix)
    return jcall(obj, "getData", Vector{Vector{jdouble}}, ())
end

function get_data_ref(obj::Array2DRowRealMatrix)
    return jcall(obj, "getDataRef", Vector{Vector{jdouble}}, ())
end

function get_entry(obj::Array2DRowRealMatrix, arg0::jint, arg1::jint)
    return jcall(obj, "getEntry", jdouble, (jint, jint), arg0, arg1)
end

function get_row(obj::Array2DRowRealMatrix, arg0::jint)
    return jcall(obj, "getRow", Vector{jdouble}, (jint,), arg0)
end

function get_row_dimension(obj::Array2DRowRealMatrix)
    return jcall(obj, "getRowDimension", jint, ())
end

function get_sub_matrix(obj::Array2DRowRealMatrix, arg0::jint, arg1::jint, arg2::jint, arg3::jint)
    return jcall(obj, "getSubMatrix", RealMatrix, (jint, jint, jint, jint), arg0, arg1, arg2, arg3)
end

function kronecker_product(obj::Array2DRowRealMatrix, arg0::RealMatrix)
    return jcall(obj, "kroneckerProduct", RealMatrix, (RealMatrix,), arg0)
end

function multiply(obj::Array2DRowRealMatrix, arg0::Array2DRowRealMatrix)
    return jcall(obj, "multiply", Array2DRowRealMatrix, (Array2DRowRealMatrix,), arg0)
end

function multiply_entry(obj::Array2DRowRealMatrix, arg0::jint, arg1::jint, arg2::jdouble)
    return jcall(obj, "multiplyEntry", void, (jint, jint, jdouble), arg0, arg1, arg2)
end

function multiply_transposed(obj::Array2DRowRealMatrix, arg0::Array2DRowRealMatrix)
    return jcall(obj, "multiplyTransposed", RealMatrix, (Array2DRowRealMatrix,), arg0)
end

function multiply_transposed(obj::Array2DRowRealMatrix, arg0::RealMatrix)
    return jcall(obj, "multiplyTransposed", RealMatrix, (RealMatrix,), arg0)
end

function operate(obj::Array2DRowRealMatrix, arg0::Vector{jdouble})
    return jcall(obj, "operate", Vector{jdouble}, (Vector{jdouble},), arg0)
end

function pre_multiply(obj::Array2DRowRealMatrix, arg0::Vector{jdouble})
    return jcall(obj, "preMultiply", Vector{jdouble}, (Vector{jdouble},), arg0)
end

function set_entry(obj::Array2DRowRealMatrix, arg0::jint, arg1::jint, arg2::jdouble)
    return jcall(obj, "setEntry", void, (jint, jint, jdouble), arg0, arg1, arg2)
end

function set_row(obj::Array2DRowRealMatrix, arg0::jint, arg1::Vector{jdouble})
    return jcall(obj, "setRow", void, (jint, Vector{jdouble}), arg0, arg1)
end

function set_sub_matrix(obj::Array2DRowRealMatrix, arg0::Vector{Vector{jdouble}}, arg1::jint, arg2::jint)
    return jcall(obj, "setSubMatrix", void, (Vector{Vector{jdouble}}, jint, jint), arg0, arg1, arg2)
end

function stack(obj::Array2DRowRealMatrix)
    return jcall(obj, "stack", RealMatrix, ())
end

function subtract(obj::Array2DRowRealMatrix, arg0::Array2DRowRealMatrix)
    return jcall(obj, "subtract", Array2DRowRealMatrix, (Array2DRowRealMatrix,), arg0)
end

function transpose_multiply(obj::Array2DRowRealMatrix, arg0::Array2DRowRealMatrix)
    return jcall(obj, "transposeMultiply", RealMatrix, (Array2DRowRealMatrix,), arg0)
end

function transpose_multiply(obj::Array2DRowRealMatrix, arg0::RealMatrix)
    return jcall(obj, "transposeMultiply", RealMatrix, (RealMatrix,), arg0)
end

function unstack_square(obj::Array2DRowRealMatrix)
    return jcall(obj, "unstackSquare", RealMatrix, ())
end

function walk_in_column_order(obj::Array2DRowRealMatrix, arg0::RealMatrixChangingVisitor)
    return jcall(obj, "walkInColumnOrder", jdouble, (RealMatrixChangingVisitor,), arg0)
end

function walk_in_column_order(obj::Array2DRowRealMatrix, arg0::RealMatrixChangingVisitor, arg1::jint, arg2::jint, arg3::jint, arg4::jint)
    return jcall(obj, "walkInColumnOrder", jdouble, (RealMatrixChangingVisitor, jint, jint, jint, jint), arg0, arg1, arg2, arg3, arg4)
end

function walk_in_column_order(obj::Array2DRowRealMatrix, arg0::RealMatrixPreservingVisitor)
    return jcall(obj, "walkInColumnOrder", jdouble, (RealMatrixPreservingVisitor,), arg0)
end

function walk_in_column_order(obj::Array2DRowRealMatrix, arg0::RealMatrixPreservingVisitor, arg1::jint, arg2::jint, arg3::jint, arg4::jint)
    return jcall(obj, "walkInColumnOrder", jdouble, (RealMatrixPreservingVisitor, jint, jint, jint, jint), arg0, arg1, arg2, arg3, arg4)
end

function walk_in_row_order(obj::Array2DRowRealMatrix, arg0::RealMatrixChangingVisitor)
    return jcall(obj, "walkInRowOrder", jdouble, (RealMatrixChangingVisitor,), arg0)
end

function walk_in_row_order(obj::Array2DRowRealMatrix, arg0::RealMatrixChangingVisitor, arg1::jint, arg2::jint, arg3::jint, arg4::jint)
    return jcall(obj, "walkInRowOrder", jdouble, (RealMatrixChangingVisitor, jint, jint, jint, jint), arg0, arg1, arg2, arg3, arg4)
end

function walk_in_row_order(obj::Array2DRowRealMatrix, arg0::RealMatrixPreservingVisitor)
    return jcall(obj, "walkInRowOrder", jdouble, (RealMatrixPreservingVisitor,), arg0)
end

function walk_in_row_order(obj::Array2DRowRealMatrix, arg0::RealMatrixPreservingVisitor, arg1::jint, arg2::jint, arg3::jint, arg4::jint)
    return jcall(obj, "walkInRowOrder", jdouble, (RealMatrixPreservingVisitor, jint, jint, jint, jint), arg0, arg1, arg2, arg3, arg4)
end


text:
## Display graph with weigthed edges
Lwidth=decisionTree.Edges.Weight*19+1
Lwidth=10*Lwidth/maximum(Lwidth)
graphPlot=plot(decisionTree,"Layout','layered','NodeLabel",nodeLabel,...
    "LineWidth",Lwidth)
toHL=find(decisionTree.Nodes.Sum~=0)
for HL=1:length(toHL)
highlight[graphPlot,toHL[HL],"Marker','s','NodeColor','r','MarkerSize",decisionTree.Nodes.Sum[toHL[HL]]+5]
end
toHL=find(decisionTree.Nodes.Sum==maximum(decisionTree.Nodes.Sum))
highlight[graphPlot,toHL,"NodeColor','y"]

ext: jl

max_stars_repo_name: DesignEngrLab/ColorGraph

max_issues_repo_name: DesignEngrLab/ColorGraph

max_forks_repo_name: DesignEngrLab/ColorGraph


text:
using ErrorTypes
using Test

@testset "Result" begin
    # instantiation
    @test Result{Bool, String}(Ok(true)) isa Result{Bool,String}
    @test Result{Bool, String}(Ok(false)) isa Result{Bool,String}
    @test Result{Int, Int}(Err(15)) isa Result{Int, Int}
    @test_throws MethodError Err{Int}(1.0)
    @test_throws MethodError Ok{Bool}(5)

    # conversion
    @test convert(Result{String, Int}, Ok("foo")) isa Result{String, Int}
    @test convert(Result{String, Int}, Ok("foo")).x isa Ok{String}
    @test convert(Result{String, Int}, Err(15)) isa Result{String, Int}
    @test convert(Result{String, Int}, Err(15)).x isa Err{Int}
    @test convert(Result{Union{Int, UInt}, Int}, Ok(15)) isa Result{Union{Int, UInt}, Int}
    @test convert(Result{Union{String, Int}, Int}, Ok("foo")).x isa Ok{Union{String, Int}}

    # convert a result directly to another if all the subtypes match
    x = Ok{Bool}(true)
    @test convert(Result{Bool, Bool}, x).x isa Ok{Bool}
    @test convert(Result{Integer, AbstractString}, Err{String}("boo")) isa Result{Integer, AbstractString}
    @test convert(Result{Signed, String}, Ok{Int32}(Int32(9))) isa Result{Signed, String}
    @test_throws MethodError convert(Result{Dict, Array}, Ok{Int}(1))

    # Errors if converting from different result type
    @test_throws MethodError convert(Result{UInt, Int}, Ok{Int}(1))
    @test_throws MethodError convert(Result{UInt, Int}, Err{Int}(UInt(1)))
    @test_throws MethodError convert(Result{Dict, Integer}, Ok{Bool}(true))
    @test_throws MethodError convert(Result{Signed, String}, Err{AbstractArray}([]))

    @test_throws MethodError convert(Result{Int, String}, Ok("foo"))
    @test_throws MethodError convert(Result{Int, Int}, Err("foo"))
    @test_throws MethodError convert(Result{String, Int}, Err("foo"))
    @test_throws MethodError convert(Result{Int, String}, Ok([1,2,3]))

    # is_error
    @test is_error(convert(Result{Int, Bool}, Err(false)))
    @test is_error(Result{String, AbstractArray}(Err([1,2])))
    @test !is_error(convert(Result{Bool, AbstractFloat}, Ok(true)))
    @test !is_error(Result{Int, Dict}(Ok(11)))

    # unwrap
    @test unwrap(Result{Int, String}(Ok(9))) === 9
    @test unwrap(convert(Result{Dict, Dict}, Ok(Dict(1=>5)))) == Dict(1=>5)
    @test_throws ErrorException unwrap(Result{Int, String}(Err("foo")))
    @test_throws ErrorException unwrap(Result{Bool, Bool}(Err(true)))

    # unwrap_or
    @test unwrap_or(Result{Int, Float64}(Ok(5)), 11) === 5
    @test unwrap_or(Result{Nothing, Nothing}(Ok(nothing)), 1) === nothing
    @test unwrap_or(Result{String, Dict}(Err(Dict())), "hello") == "hello"
    @test unwrap_or(Result{Int, Int}(Err(1)), 1 + 1) === 2

    # unwrap_error_or
    @test unwrap_error_or(Result{Integer, AbstractFloat}(Err(1.1)), 1.2) === 1.1
    @test unwrap_error_or(Result{UInt8, UInt8}(Err(0x00)), Dict()) === 0x00
    @test unwrap_error_or(Result{Dict, String}(Ok(Dict())), missing) === missing
    @test unwrap_error_or(Result{Int, UInt}(Ok(55)), UInt(1)) === UInt(1)

    # expect
    @test expect(Result{Nothing, Int}(Ok(nothing)), "foo") === nothing
    @test expect(convert(Result{Int,Int}, Ok(17)), " ") === 17
    @test_throws ErrorException expect(Result{Int, UInt}(Err(UInt(11))), "bar")
    @test_throws ErrorException expect(Result{Bool, Bool}(Err(true)), "foo")

    # unwrap_error
    @test unwrap_error(Result{String, UInt}(Err(UInt(19)))) == UInt(19)
    @test unwrap_error(Result{AbstractDict, String}(Err("bar"))) == "bar"
    @test_throws ErrorException unwrap_error(Result{Int, Int}(Ok(1)))
    @test_throws ErrorException unwrap_error(Result{String, Int}(Ok("bar")))

    # expect_err
    @test expect_error(Result{Vector, AbstractString}(Err("x")), "foo") == "x"
    @test expect_error(Result{Int, UInt8}(Err(0xa4)), "mistake!") == 0xa4
    @test_throws ErrorException expect_error(Result{Int, UInt8}(Ok(15)), "foobar")
    @test_throws ErrorException expect_error(Result{Vector, AbstractString}(Ok([])), "xx")

    # and_then
    @test and_then(x -> x + 1, Float64, Result{Float64, String}(Ok(5.0))) === Result{Float64, String}(Ok(6.0))
    @test and_then(x -> abs(x) % UInt8, UInt8, Result{Int, Dict}(Ok(-55))) === Result{UInt8, Dict}(Ok(UInt8(55)))
    @test and_then(x -> x + 1, Bool, Result{String, Int}(Err(1))) === Result{Bool, Int}(Err(1))
    @test and_then(x -> "foo", Dict, Result{Int, Int}(Err(1))) === Result{Dict, Int}(Err(1))

    # map_or
    @test map_or(x -> x + 1, Result{Float64, String}(Ok(3.0)), 3.0) === 4.0
    @test map_or(x -> x + 1, Result{Float64, String}(Err("foo")), 3.0) === 3.0
    @test map_or(x -> 0, Result{Int, Int}(Ok(5)), 1) === 0
    @test map_or(x -> 0, Result{Int, Int}(Err(19)), 1) === 1
end

@testset "Option" begin
    # instantiation
    @test none(Int) isa Option{Int}
    @test !(none(Float32) isa Option{AbstractFloat})
    @test !(none(Int) isa Err)

    # Conversion from none
    @test convert(Option{Int}, none) isa Option{Int}
    @test convert(Option{Some{Nothing}}, none) isa Option{Some{Nothing}}
    @test_throws MethodError convert(Int, none)
    @test_throws MethodError convert(Option{String}, "foo")

    # Convert from Option
    @test convert(Option{Integer}, some(1)) isa Option{Integer}
    @test convert(Option{AbstractString}, none(String)) isa Option{AbstractString}
    @test convert(Option{Union{String, Dict}}, some("foo")) isa Option{Union{String, Dict}}
    @test convert(Option{Float64}, none(Float64)) isa Option{Float64}
    @test_throws MethodError convert(Option{String}, some(1))
    @test_throws MethodError convert(Option{String}, none(Int))

    # flatten
    @test flatten(some(some(11))) === some(11)
    @test flatten(some(some("foo"))) == some("foo")
    @test flatten(some(none(Float64))) == none(Float64)
    @test flatten(none(Option{Dict{Int, Int}})) == none(Dict{Int, Int})
    @test flatten(some(some(none(Int)))) == some(none(Int))
    @test_throws MethodError flatten(some(Result{Int, UInt}(Err(UInt(1)))))
    @test_throws MethodError flatten(Result{Option{Int}, Missing}(Ok(some(5))))

    # base
    @test base(some(11)) == Some(11)
    @test base(none(AbstractString)) === nothing
    @test base(some(some(0xf1))) == Some(some(0xf1))
    @test_throws MethodError base(Result{String, Vector}(Err([1,2])))
    @test_throws MethodError base(Result{Int, Int}(Ok(1)))
end

@testset "Result/Option conversion" begin
    # Construct Option from Result
    @test Option(Result{String, Int}(Ok("foo"))) == some("foo")
    @test Option(Result{String, Int}(Err(15))) == none(String)

    # Construct Result from Option
    @test Result(none(String), 0x01) == Result{String, UInt8}(Err(0x01))
    @test Result(some(0x1234), []) == Result{UInt16, Vector{Any}}(Ok(0x1234))
    @test Result(some(1), 2) == Result{Int, Int}(Ok(1))
end

@testset "@?" begin
    # Generic Option usage
    function testf1(x)::Option{typeof(x)}
        iszero(x) ? none : some(x + one(x))
    end

    function testf2(x)::Option{String}
        some(string(@? testf1(x)))
    end

    @test is_error(testf2(0))
    @test unwrap(testf2(1)) == "2"

    # Generic Result usage
    function testf3(x)::Result{Int, Union{Nothing, Missing}}
        if iszero(x)
            return Err(nothing)
        elseif isone(x)
            return Err(missing)
        else
            return Ok(x + 1)
        end
    end

    function testf4(x)::Result{Float64, Union{Nothing, Missing}}
        y = @? testf3(x)
        Ok(y + 0.5)
    end

    @test is_error(testf4(0))
    @test is_error(testf4(1))
    @test !is_error(testf4(3))
    @test unwrap(testf4(5)) === 6.5

    @test_throws TypeError @? 1 + 1

    # Convert from Result to Option
    testf5(x::Int)::Result{Int, Int} = iszero(x) ? Err(x) : Ok(x + 1)
    testf6(x)::Option{Int} = some(@?(testf5(x)) + 2)

    @test testf6(0) == none(Int)
    @test testf6(1) == some(4)
end

@testset "@unwrap_or" begin
    # @unwrap_or
    not_zero(x::Int)::Option{Int} = iszero(x) ? none : some(1 + x)

    function my_sum(f, it)
        sum = 0
        for i in it
            sum += @unwrap_or f(i) continue
        end
        sum
    end

    @test my_sum(not_zero, [1,2,0,3,0,1]) == 11

    function not_zero_two(x::Int)::Result{Float64, String}
        iszero(x) && return Err("Zero")
        isnan(x) && return Err("NaN")
        Ok(1 / x)
    end

    @test my_sum(not_zero_two, [1,2,0,2,0,1]) === 3.0

    @test_throws TypeError @unwrap_or (1 + 1) "foo"

    # @unwrap_error_or
    @test (@unwrap_error_or none(Integer) missing) === nothing
    @test (@unwrap_error_or some(55) missing) === missing
    @test (@unwrap_error_or Result{UInt8, String}(Err("Someerr")) missing) == "Someerr"
end

# I only have this to get 100% code coverage and make sure printing doesn't error
# I don't actually really test what it printed
@testset "printing" begin
    function test_printing(x)
        io = IOBuffer()
        show(io, x)
        return length(take!(io))
    end
    @test test_printing(Result{Int, Bool}(Err(false))) > 5
    @test test_printing(some(1)) > 5
    @test test_printing(none(AbstractString)) > 5
    @test test_printing(some(1)) < test_printing(Option{Integer}(some(1)))
end


text:
#=
Response surfaces implementantion, following:
"A Taxonomy of Global Optimization Methods Based on Response Surfaces"
by DONALD R. JONES
=#
mutable struct RadialBasis{F,X,Y,B,C} <: AbstractSurrogate
    phi::F
    dim_poly::Int
    x::X
    y::Y
    bounds::B
    coeff::C
end

"""
Calculates current estimate of array value 'val' with respect to RadialBasis object.
"""
function (rad::RadialBasis)(val)
    n = length(rad.x)
    d = length(rad.x[1])
    q = rad.dim_poly
    central_point = zeros(eltype(rad.x[1]), d)
    sum = zero(eltype(rad.x[1]))
    @inbounds for i = 1:d
        central_point[i] = (rad.bounds[i][1]+rad.bounds[i][2])/2
        sum += (rad.bounds[i][2]-rad.bounds[i][1])/2
    end
    half_diameter_domain = sum/d
    approx = zero(eltype(rad.x[1]))
    for i = 1:n
        approx = approx + rad.coeff[i]*rad.phi(val .- rad.x[i])
    end
    for i = n+1:n+q
        approx = approx + rad.coeff[i]*centralized_monomial(val,n+1-i,half_diameter_domain,central_point)
    end
    return approx
end

"""
Calculates current estimate of value 'val' with respect to RadialBasis object.
"""
function (rad::RadialBasis)(val::Number)
    approx = zero(eltype(rad.x[1]))
    n = length(rad.x)
    q = rad.dim_poly
    Chebyshev(x,k) = cos(k*acos(-1 + 2/(rad.bounds[2]-rad.bounds[1])*(x.-rad.bounds[1])))
    for i = 1:n
        approx = approx + rad.coeff[i]*rad.phi(val .- rad.x[i])
    end
    for i = n+1:n+q
        approx = approx + rad.coeff[i]*Chebyshev(val,n+1-i)
    end
    return approx
end

#=
linear_basis_function = Basis(z->norm(z), 1)
cubic_basis_function = Basis(z->norm(z)^3, 2)
thinplate_basis_function = Basis(z->norm(z)^2*log(norm(z)),2)
function multiquadric_basis_function(lambda)
    return Basis(z->sqrt(norm(z)^2 + lambda^2),1)
end
=#

"""
    RadialBasis(x,y,a::Number,b::Number,phi::Function,q::Int)

Constructor for RadialBasis surrogate
- (x,y): sampled points
- '(a,b)': interval of interest
-'phi': radial basis of choice
-'q': number of polynomial elements
"""
function RadialBasis(x,y,a::Number,b::Number,phi::Function,q::Int)
    coeff = _calc_coeffs(x,y,a,b,phi,q)
    RadialBasis(phi,q,x,y,(a,b),coeff)
end

function _calc_coeffs(x,y,a,b,phi,q)
    Chebyshev(x,k) = cos(k*acos(-1 + 2/(b-a)*(x-a)))
    n = length(x)
    size = n+q
    D = zeros(eltype(x[1]), size, size)
    d = zeros(eltype(x[1]),size)

    @inbounds for i = 1:n
        d[i] =  y[i]
        for j = 1:n
            D[i,j] = phi(x[i] .- x[j])
        end
        if i < n + 1
            for k = n+1:size
                D[i,k] = Chebyshev(x[i],k)
            end
        end
    end
    Sym = Symmetric(D,:U)
    coeff = Sym\d
end

"""
RadialBasis(x,y,bounds,phi::Function,q::Int)

Constructor for RadialBasis surrogate

- (x,y): sampled points
- bounds: region of interest of the form [[a,b],[c,d],...,[w,z]]
- phi: radial basis of choice
- q: number of polynomial elements
"""
function RadialBasis(x,y,bounds,phi::Function,q::Int)
    coeff = _calc_coeffs(x,y,bounds,phi,q)
    RadialBasis(phi,q,x,y,bounds,coeff)
end

function _calc_coeffs(x,y,bounds,phi,q)
    n = length(x)
    d = length(x[1])
    central_point = zeros(eltype(x[1]), d)
    sum = zero(eltype(x[1]))
    @inbounds for i = 1:d
        central_point[i] = (bounds[i][1]+bounds[i][2])/2
        sum += (bounds[i][2]-bounds[i][1])/2
    end
    half_diameter_domain = sum/d

    size = n+q
    D = zeros(eltype(x[1]), size, size)
    d = zeros(eltype(x[1]),size)

    @inbounds for i = 1:n
        d[i] =  y[i]
        for j = 1:n
            D[i,j] = phi(x[i] .- x[j])
        end
        if i < n + 1
            for k = n+1:size
                D[i,k] = centralized_monomial(x[i],k,half_diameter_domain,central_point)
            end
        end
    end
    Sym = Symmetric(D,:U)
    coeff = Sym\d
end

"""
    centralized_monomial(vect,alpha,half_diameter_domain,central_point)

Returns the value at point vect[] of the alpha degree monomial centralized.

#Arguments:
-'vect': vector of points i.e [x,y,...,w]
-'alpha': degree
-'half_diameter_domain': half diameter of the domain
-'central_point': central point in the domain
"""
function centralized_monomial(vect,alpha,half_diameter_domain,central_point)
    mul = 1
    @inbounds for i = 1:length(vect)
        mul *= vect[i]
    end
    return ((mul-norm(central_point))/(half_diameter_domain))^alpha
end

"""
    add_point!(rad::RadialBasis,new_x,new_y)

Add new samples x and y and updates the coefficients. Return the new object radial.
"""
function add_point!(rad::RadialBasis,new_x,new_y)
    if (length(new_x) == 1 && length(new_x[1]) == 1) || ( length(new_x) > 1 && length(new_x[1]) == 1 && length(rad.bounds[1])>1)
        push!(rad.x,new_x)
        push!(rad.y,new_y)
        if length(rad.bounds[1]) == 1
            rad.coeff = _calc_coeffs(rad.x,rad.y,rad.bounds[1],rad.bounds[2],rad.phi,rad.dim_poly)
        else
            rad.coeff = _calc_coeffs(rad.x,rad.y,rad.bounds,rad.phi,rad.dim_poly)
        end
    else
        append!(rad.x,new_x)
        append!(rad.y,new_y)
        if length(rad.bounds[1]) == 1
            rad.coeff = _calc_coeffs(rad.x,rad.y,rad.bounds[1],rad.bounds[2],rad.phi,rad.dim_poly)
        else
            rad.coeff = _calc_coeffs(rad.x,rad.y,rad.bounds,rad.phi,rad.dim_poly)
        end
    end
    nothing
end


text:
#########################################################################
# Reading
#########################################################################
"""
    ncread(file, var [, selection]; kwargs...) → A
Load the variable `var` from the `file` and convert it into a [`ClimArray`](@ref)
with proper dimension mapping and also containing the variable attributes as a dictionary.
Dimension attributes are also given to the dimensions of `A`, if any exist.
See keywords below for specifications for unstructured grids.

`file` can be a string to a `.nc` file. Or, it can be an
`NCDataset`, which allows you to lazily combine different
`.nc` data (typically split by time), e.g.
```julia
using Glob # for getting all files
alldata = glob("toa_fluxes_*.nc")
file = NCDataset(alldata; aggdim = "time")
A = ClimArray(file, "tow_sw_all")
```

`var` is a `String` denoting which variable to load.
For `.nc` data containing groups `var` can also be a tuple `("group_name", "var_name")`
that loads a specific variable from a specific group.
In this case, the attributes of both the group and the CF-variable are attributed to
the created [`ClimArray`](@ref).

Optionally you can provide a `selection` for selecting a smaller part of the full array.
The `selection` must be a tuple of indices that compose the selection and you must specify
exactly as many ranges as the dimensions of the array and in the correct order.
For example, if `var` corresponds
to an array with three dimensions, such syntaxes are possible:
```julia
(:, :, 1:3)
(1:5:100, 1:1, [1,5,6])
```
The function [`ncsize`](@ref) can be useful for `selection`.

See also [`ncdetails`](@ref), [`nckeys`](@ref) and [`ncwrite`](@ref).

## Smart loading
The following things make loading data with `ncread` smarter than directly trying to use
NCDatasets.jl and then convert to some kind of dimensional container.
1. Data are directly transformed into `ClimArray`, conserving metadata and dimension names.
1. If there are no missing values in the data (according to CF standards), the
   returned array is automatically converted to a concrete type (i.e. `Union{Float32, Missing}`
   becomes `Float32`).
1. Dimensions that are ranges (i.e. sampled with constant step size) are automatically
   transformed to a standard Julia `Range` type (which makes sub-selecting faster).
1. Automatically deducing whether the spatial information is in an orthogonal
   grid or not, and creating a single `Coord` dimension if not.

## Keywords
* `name` optionally rename loaded array.
* `grid = nothing` optionally specify whether the underlying grid is `grid = OrthogonalSpace()`
  or `grid = CoordinateSpace()`, see [Types of spatial information](@ref).
  If `nothing`, we try to deduce automatically based on
  the names of dimensions and other keys of the `NCDataset`.
* `lon, lat`. These two keywords are useful in unstructured grid data where the grid
  information is provided in a *separate .nc file*. What we need is the user to
  provide vectors of the central longitude and central latitude of each grid point.
  This is done e.g. by
  ```julia
  ds = NCDataset("path/to/grid.nc");
  lon = Array(ds["clon"]);
  lat = Array(ds["clat"]);
  ```
  If `lon, lat` are given, `grid` is automatically assumed `CoordinateSpace()`.
"""
function ncread(path::Union{String, Vector{String}}, args...; kwargs...)
    NCDataset(path) do ds
        data = ncread(ds, args...; kwargs...)
        return data
    end
end

# TODO: Allow reading multiple variables at once. This has the performance benefit
# of not re-creating dimensions all the time.

function ncread(ds::NCDatasets.AbstractDataset, var, selection = nothing;
        name = var2name(var), grid = nothing, lon = nothing, lat = nothing,
    )
    if lon isa Vector && lat isa Vector
        gridtype = CoordinateSpace()
    elseif isnothing(grid)
        gridtype = autodetect_grid(ds)
    else
        gridtype = grid
    end

    if gridtype == CoordinateSpace()
        return ncread_unstructured(ds, var, name, lon, lat, selection)
    else
        return ncread_lonlat(ds, var, name, selection)
    end
end

function autodetect_grid(ds)
    if haskey(ds, "reduced_points") || haskey(ds, "clon") ||
        any(x -> x ∈ ds.dim, POSSIBLE_CELL_NAMES)
        # Common cases of coordinate spaces in NetCDF
        return CoordinateSpace()
    elseif haskey(ds, "lat") && length(size(ds["lat"])) > 1
        # Uncommon case of storing lon/lat as matrices, used in some weird ocean grids
        return CoordinateSpace()
    else
        return OrthogonalSpace()
    end
end

var2name(var) = string(var)
var2name(var::Tuple) = join(var, "_")

#########################################################################
# Reading: OrthogonalSpace and main reading functionality
#########################################################################
# Notice that this function properly loads even without any spatial coordinate
function ncread_lonlat(ds::NCDatasets.AbstractDataset, var, name, selection)
    cfvar = var isa String ? ds[var] : ds.group[var[1]][var[2]]
    sel = isnothing(selection) ? selecteverything(cfvar) : selection
    @assert length(sel) == length(size(cfvar))
    attrib = get_attributes_from_var(ds, cfvar, var)
    A = cfvar[sel...]
    dnames = Tuple(NCDatasets.dimnames(cfvar))
    if !any(ismissing, A)
        A = nomissing(A)
    end
    dimensions = create_dims(ds, dnames, sel)
    data = ClimArray(A, dimensions; name = Symbol(name), attrib = attrib)
    return data
end
selecteverything(cfvar) = map(i -> 1:i, size(cfvar))

get_attributes_from_var(ds, cfvar, var::String) = Dict(cfvar.attrib)
function get_attributes_from_var(ds, cfvar, var::Tuple)
    d1 = Dict(cfvar.attrib)
    d2 = Dict(ds.group[var[1]].attrib)
    return merge!(d2, d1)
end

"""
    create_dims(ds::NCDatasets.AbstractDataset, dnames, cfvar, sel = selecteverything(A))
Create a tuple of `Dimension`s from the `dnames` (tuple of strings).
"""
function create_dims(ds::NCDatasets.AbstractDataset, dnames, sel = selecteverything(A))
    true_dims = to_proper_dimensions(dnames)
    dim_values = extract_dim_values(ds, dnames, sel)
    # Some stupid datasets return a union{Missing} type for dimensions.
    # this is of course nonsense, a dimension cannot have "missing" values.
    if any(d -> Missing <: eltype(d), dim_values)
        dim_values = nomissing.(dim_values)
    end
    optimal_values = vector2range.(dim_values)
    attribs = [
        (haskey(ds, d) && ds[d].attrib isa NCDatasets.BaseAttributes) ?
            Dict(ds[d].attrib) :
            Dict()
        for d in dnames
    ]
    out = Dimension[]
    for i in 1:length(true_dims)
        push!(out, true_dims[i](optimal_values[i]; metadata = attribs[i]))
    end
    return (out...,)
end

function to_proper_dimensions(dnames)
    r = []
    for n in dnames
        if haskey(COMMONNAMES, n)
            push!(r, COMMONNAMES[n])
        else
            @warn """
            Dimension name "$n" not in common names.
            Please consider opening an issue/PR on GitHub proposing "$n" in common names!
            Making a generic dimension named `Dim{:$n}` for now...
            """
            push!(r, Dim{Symbol(n)})
        end
    end
    return (r...,)
end
export Dim # for generic dimensions this must be exported

function extract_dim_values(ds::NCDatasets.AbstractDataset, dnames, sel)
    dim_values = AbstractVector[]
    for (i, n) in enumerate(dnames)
        if haskey(ds, n)
            push!(dim_values, vec(ds[n][sel[i]]))
        else
            @warn """
            Dimension named "$n" does not have values in the dataset.
            Using the range `1:size(cfvar, i)` as the dimension values instead,
            where `i` is the dimension index.
            """
            push!(dim_values, sel[i])
        end
    end
    return dim_values
end

#########################################################################
# Reading: CoordinateSpace
#########################################################################
export SVector

function ncread_unstructured(
        ds::NCDatasets.AbstractDataset, var::String, name, lon, lat, selection
    )
    cfvar = var isa String ? ds[var] : ds.group[var[1]][var[2]]
    sel = isnothing(selection) ? selecteverything(cfvar) : selection
    attrib = get_attributes_from_var(ds, cfvar, var)

    # Here we generate the longitudes and latitudes based on whether we have
    # reduced points or not, and obtain the name of the coord dimension in the `ds`
    if isnothing(lon)
        lonlat, original_grid_dim = load_coordinate_points(ds)
    else
        lonlat = [SVector(lo, la) for (lo, la) in zip(lon, lat)]
        original_grid_dim = intersect(NCDatasets.dimnames(cfvar), POSSIBLE_CELL_NAMES)[1]
    end

    alldims = Any[NCDatasets.dimnames(cfvar)...]

    # Set up the remaining dimensions of the dataset
    if original_grid_dim isa Tuple
        # stupid case where `lon` data are `Matrix`

        # TODO: I'm not sure I've set up this to work correctly in the case
        # of `selection` given by user... But the `Matrix` case is so rare
        # I honestly don't care at the moment.

        A = cfvar[sel...]
        sizes = [size(A)...]
        # First, find where lon, lat dimensions are positioned
        is = findall(x -> x ∈ original_grid_dim, alldims)
        @assert length(is) == 2
        remainingdims = deleteat!(copy(alldims), is)
        actualdims = Any[create_dims(ds, remainingdims, A)...]
        # Then, reshape A accordingly
        sizes[is[1]] = sizes[is[1]]*sizes[is[2]]
        deleteat!(sizes, is[2])
        A = reshape(A, sizes...)
        i = is[1]
    else
        @assert original_grid_dim ∈ alldims
        i = findfirst(x -> x == original_grid_dim, alldims)
        remainingdims = deleteat!(copy(alldims), i)
        A = cfvar[sel...]
        remainingsel = Tuple(deleteat!([sel...], i))
        actualdims = Any[create_dims(ds, remainingdims, remainingsel)...]
    end

    # Make coordinate dimension
    lonlat = lonlat[sel[i]]
    si = sortperm(lonlat; by = reverse)
    coords = Coord(lonlat, (Lon, Lat))
    insert!(actualdims, i, coords)

    # Create array, sort by latitude and remove missings
    X = ClimArray(A, Tuple(actualdims))
    X = X[Coord(si)]
    if !any(ismissing, X)
        X = nomissing(X)
    end
    return ClimArray(X; name = Symbol(name), attrib)
end

function has_unstructured_key(ds)
    any(x -> haskey(ds.dim, x), POSSIBLE_CELL_NAMES) ||
    haskey(ds, "lon") || haskey(ds, "clon")
end

function load_coordinate_points(ds)
    if haskey(ds, "reduced_points")
        lonlat = reduced_grid_to_points(ds["lat"], ds["reduced_points"])
        original_grid_dim = "rgrid" # Specific to CDO Gaussian grid
    elseif haskey(ds, "lon") && length(size(ds["lon"])) == 2
        # This is the case of having a non-orthogonal grid, but
        # still saving the lon/lat information as matrices whose dimensions are lon/lat
        # (this is a rather stupid way to format things, unfortunately)
        lons = ds["lon"] |> Matrix |> vec
        lats = ds["lat"] |> Matrix |> vec
        original_grid_dim = ("lon", "lat")
        lonlat = [SVector(lo, la) for (lo, la) in zip(lons, lats)]
    elseif has_unstructured_key(ds)
        if haskey(ds, "lon")
            lons = ds["lon"] |> Vector
            lats = ds["lat"] |> Vector
            original_grid_dim = NCDatasets.dimnames(ds["lon"])[1]
        elseif haskey(ds, "clon")
            lons = ds["clon"] |> Vector
            lats = ds["clat"] |> Vector
            original_grid_dim = NCDatasets.dimnames(ds["clon"])[1]
        else
            error("""
            We didn't find key `"lon"` or `"clon"` that represents the longitude of each
            polygon in a non-orthogonal grid.
            """)
        end
        lonlat = [SVector(lo, la) for (lo, la) in zip(lons, lats)]
    else
        error("""
        We couldn't automatically identify the lon/lat values of cell centers.
        Please provide explicitly keywords `lon, lat` in `ncread`.
        """)
    end
    lonlat = convert_to_degrees(lonlat, ds)
    return lonlat, original_grid_dim
end

function reduced_grid_to_points(lat, reduced_points)
    lonlat = SVector{2, Float32}[]
    for (i, θ) in enumerate(lat)
        n = reduced_points[i]
        dλ = Float32(360/n)
        for j in 0:n-1
            push!(lonlat, SVector(0 + dλ*j, θ))
        end
    end
    return lonlat
end

function convert_to_degrees(lonlat, ds)
    !(haskey(ds, "clat") || haskey(ds, "lat")) && return convert_to_degrees(lonlat)
    x = haskey(ds, "clat") ? ds["clat"] : ds["lat"]
    if get(x.attrib, "units", nothing) == "radian" || !any(ll -> abs(ll[2]) > π/2, lonlat)
        lonlat = [SVector(lo*180/π, la*180/π) for (lo, la) in lonlat]
    end
    return lonlat
end
function convert_to_degrees(lonlat)
    if !any(ll -> abs(ll[2]) > π/2, lonlat)
        lonlat = [SVector(lo*180/π, la*180/π) for (lo, la) in lonlat]
    end
    return lonlat
end


text:
export ExperienceBasedSamplingModel, sample

using Random
import StatsBase: sample

"""
    ExperienceBasedSamplingModel

Randomly generate a transition of (s, a, r, t, s′) based on previous experiences
in each sampling.
"""
Base.@kwdef mutable struct ExperienceBasedSamplingModel{R} <: AbstractEnvironmentModel
    experiences::Dict{Any,Dict{Any,Any}} = Dict{Any,Dict{Any,Any}}()
    sample_count::Int = 0
    rng::R = Random.GLOBAL_RNG
end

function RLBase.update!(
    m::ExperienceBasedSamplingModel,
    t::AbstractTrajectory,
    ::AbstractPolicy,
    ::AbstractEnv,
    ::Union{PreActStage,PostEpisodeStage},
)
    if length(t[:terminal]) > 0
        transition = (
            t[:state][end-1],
            t[:action][end-1],
            t[:reward][end],
            t[:terminal][end],
            t[:state][end],
        )
        update!(m, transition)
    end
end

function RLBase.update!(m::ExperienceBasedSamplingModel, transition::Tuple)
    s, a, r, d, s′ = transition
    if haskey(m.experiences, s)
        m.experiences[s][a] = (reward = r, terminal = d, nextstate = s′)
    else
        m.experiences[s] = Dict(a => (reward = r, terminal = d, nextstate = s′))
    end
end

sample(model::ExperienceBasedSamplingModel) = sample(model.rng, model)

function sample(rng::AbstractRNG, model::ExperienceBasedSamplingModel)
    if length(model.experiences) > 0
        s = rand(rng, keys(model.experiences))
        a = rand(rng, keys(model.experiences[s]))
        model.sample_count += 1
        s, a, model.experiences[s][a]...
    else
        nothing
    end
end


text: const global A=1


text:
#!/usr/bin/env bash
    #=
    exec julia --project="$(realpath $(dirname $0))/" "${BASH_SOURCE[0]}" "$@" -e "include(popfirst!(ARGS))" \
    "${BASH_SOURCE[0]}" "$@"
    =#

println("\033[1;34m===>\033[0;38m\033[1;38m\tLoading required libraries (it will take a moment to precompile if it is your first time doing this)...\033[0;38m")

include(joinpath(dirname(dirname(@__FILE__)), "src", "FaceDetection.jl"))

using .FaceDetection
const FD = FaceDetection
using Images: imresize
using StatsPlots  # StatsPlots required for box plots # plot boxplot @layout :origin savefig
using CSV: write
using DataFrames: DataFrame
using HypothesisTests: UnequalVarianceTTest
using Serialization: deserialize

println("...done")

function main(;
	smart_choose_feats::Bool=false,
	scale::Bool=false,
	scale_to::Tuple=(200, 200))

	include("constants.jl")
	include("main_data.jl")
    pos_testing_path = "/Users/jakeireland/projects/FaceDetection.jl/data/lizzie-testset/faces/"
	neg_testing_path = "/Users/jakeireland/projects/FaceDetection.jl/data/lizzie-testset/nonfaces/"

    # read classifiers from file
	classifiers = deserialize("/Users/jakeireland/Desktop/classifiers_577_577_fixed")

    notify_user("Calculating test face scores and constructing dataset...")

    faces_scores = Vector{Real}(undef, length(filtered_ls(pos_testing_path)))
    non_faces_scores = Vector{Real}(undef, length(filtered_ls(neg_testing_path)))

    faces_scores[:] .= [sum([FD.get_faceness(c, load_image(face, scale=scale, scale_to=scale_to)) for c in classifiers]) for face in filtered_ls(pos_testing_path)]
	non_faces_scores[:] .= [sum([FD.get_faceness(c, load_image(non_face, scale=scale, scale_to=scale_to)) for c in classifiers]) for non_face in filtered_ls(neg_testing_path)]

	face_names = filtered_ls(pos_testing_path)
	non_face_names = filtered_ls(neg_testing_path)

    # filling in the dataset with missing to easily write to csv
    df_faces = faces_scores
    df_non_faces = non_faces_scores
    if length(faces_scores) < length(non_faces_scores)
        to_add = length(filtered_ls(neg_testing_path)) - length(filtered_ls(pos_testing_path))
        df_faces = vcat(df_faces, Matrix{Union{Float64, Missing}}(undef, to_add, 1))
        face_names = vcat(face_names, Matrix{Union{Float64, Missing}}(undef, to_add, 1))
    elseif length(faces_scores) > length(non_faces_scores)
        to_add = length(filtered_ls(pos_testing_path)) - length(filtered_ls(neg_testing_path))
        df_non_faces = vcat(df_non_faces, Matrix{Union{Float64, Missing}}(undef, to_add, 1))
        non_face_names = vcat(non_face_names, Matrix{Union{Float64, Missing}}(undef, to_add, 1))
    end

    # write score data
    write(joinpath(dirname(dirname(@__FILE__)), "data", "faceness-scores.csv"), DataFrame(hcat(face_names, df_faces, non_face_names, df_non_faces)), writeheader=false)

    println("...done.\n")

    notify_user("Computing differences in scores between faces and non-faces...")

    welch_t = UnequalVarianceTTest(faces_scores, non_faces_scores)

    println("...done.  $welch_t\n")

    notify_user("Constructing box plot with said dataset...")

    gr() # set plot backend
    theme(:solarized)
    plot = StatsPlots.plot(
                    StatsPlots.boxplot(faces_scores, xaxis=false, label = false),
                    StatsPlots.boxplot(non_faces_scores, xaxis=false, label = false),
                    title = ["Scores of Faces" "Scores of Non-Faces"],
                    # label = ["faces" "non-faces"],
                    fontfamily = font("Times"),
                    layout = @layout([a b]),
                    # fillcolor = [:blue, :orange],
                    link = :y,
                    # framestyle = [:origin :origin]
                )

    # save plot
    StatsPlots.savefig(plot, joinpath(dirname(dirname(@__FILE__)), "figs", "scores.pdf"))

    println("...done.  Plot created at ", joinpath(dirname(dirname(@__FILE__)), "figs", "scores.pdf"), "\n")
end

@time main(smart_choose_feats=true, scale=true, scale_to=(577, 577))


text:
########################################################################################################
########################################################################################################

"""
"""
function form_sigma_block!(term::ClusteredTerm1B,
                            cluster_ops::Vector{ClusterOps},
                            fock_bra::FockConfig, bra::TuckerConfig,
                            fock_ket::FockConfig, ket::TuckerConfig,
                            bra_coeffs::Array, ket_coeffs::Array)
#={{{=#
    #display(term)
    #println(bra, ket)

    c1 = term.clusters[1]
    length(fock_bra) == length(fock_ket) || throw(Exception)
    length(bra) == length(ket) || throw(Exception)
    n_clusters = length(bra)
    #
    # make sure inactive clusters are diagonal
    for ci in 1:n_clusters
        ci != c1.idx || continue

        fock_bra[ci] == fock_ket[ci] || throw(Exception)
        bra[ci] == ket[ci] || return 0.0
    end

    #
    # make sure active clusters are correct transitions
    fock_bra[c1.idx] == fock_ket[c1.idx] .+ term.delta[1] || throw(Exception)

    op = cluster_ops[c1.idx][term.ops[1]][(fock_bra[c1.idx],fock_ket[c1.idx])][bra[c1.idx],ket[c1.idx]]

    # now transpose state vectors and multiply, also, try without transposing to compare
    indices = collect(1:n_clusters+1)
    indices[c1.idx] = 0
    perm,_ = bubble_sort(indices)

    length(size(ket_coeffs)) == n_clusters + 1 || error(" tensors should be folded")

    n_roots = last(size(ket_coeffs))
    ket_coeffs2 = permutedims(ket_coeffs,perm)
    bra_coeffs2 = permutedims(bra_coeffs,perm)

    dim1 = size(ket_coeffs2)
    ket_coeffs2 = reshape(ket_coeffs2, dim1[1], prod(dim1[2:end]))

    dim2 = size(bra_coeffs2)
    bra_coeffs2 = reshape(bra_coeffs2, dim2[1], prod(dim2[2:end]))

    bra_coeffs2 .+= op * ket_coeffs2
#    if bra==ket
#        display(op)
#        display(bra_coeffs2)
#        display(ket_coeffs2)
#    end

    ket_coeffs2 = reshape(ket_coeffs2, dim1)
    bra_coeffs2 = reshape(bra_coeffs2, dim2)

    # now untranspose
    perm,_ = bubble_sort(perm)
    ket_coeffs2 = permutedims(ket_coeffs2,perm)
    bra_coeffs2 = permutedims(bra_coeffs2,perm)

    bra_coeffs .= bra_coeffs2
    return
#=}}}=#
end
"""
"""
function form_sigma_block!(term::ClusteredTerm2B,
                            cluster_ops::Vector{ClusterOps},
                            fock_bra::FockConfig, bra::TuckerConfig,
                            fock_ket::FockConfig, ket::TuckerConfig,
                            bra_coeffs::Array, ket_coeffs::Array)
#={{{=#
    #display(term)
    #println(bra, ket)

    c1 = term.clusters[1]
    c2 = term.clusters[2]
    length(fock_bra) == length(fock_ket) || throw(Exception)
    length(bra) == length(ket) || throw(Exception)
    n_clusters = length(bra)
    #
    # make sure inactive clusters are diagonal
    for ci in 1:n_clusters
        ci != c1.idx || continue
        ci != c2.idx || continue

        fock_bra[ci] == fock_ket[ci] || throw(Exception)
        bra[ci] == ket[ci] || return 0.0
    end

    #display(fock_bra)
    #display(fock_ket)
    #display(term.delta)
    #display(term)
    #
    # make sure active clusters are correct transitions
    fock_bra[c1.idx] == fock_ket[c1.idx] .+ term.delta[1] || throw(Exception)
    fock_bra[c2.idx] == fock_ket[c2.idx] .+ term.delta[2] || throw(Exception)

    #
    # determine sign from rearranging clusters if odd number of operators
    state_sign = compute_terms_state_sign(term, fock_ket)

    #
    # op[IK,JL] = <I|p'|J> h(pq) <K|q|L>
#    display(term)
#    display(fock_bra)
#    display(fock_ket)
    gamma1 = cluster_ops[c1.idx][term.ops[1]][(fock_bra[c1.idx],fock_ket[c1.idx])][:,bra[c1.idx],ket[c1.idx]]
    gamma2 = cluster_ops[c2.idx][term.ops[2]][(fock_bra[c2.idx],fock_ket[c2.idx])][:,bra[c2.idx],ket[c2.idx]]

    op = Array{Float64}[]
    #display(size(term.ints))
    #display(size(gamma1))
    #display(size(gamma2))
    cache_key = (fock_bra[c1.idx], fock_bra[c2.idx], fock_ket[c1.idx], fock_ket[c2.idx], bra[c1.idx], bra[c2.idx], ket[c1.idx], ket[c2.idx])
    if haskey(term.cache, cache_key)
        op = term.cache[cache_key]
    else
        @tensor begin
            op[q,J,I] := term.ints[p,q] * gamma1[p,I,J]
            op[J,L,I,K] := op[q,J,I] * gamma2[q,K,L]
        end
        term.cache[cache_key] = op
    end

    # possibly cache some of these integrals

    # now transpose state vectors and multiply, also, try without transposing to compare
    indices = collect(1:n_clusters+1)
    indices[c1.idx] = 0
    indices[c2.idx] = 0
    perm,_ = bubble_sort(indices)

    length(size(ket_coeffs)) == n_clusters + 1 || error(" tensors should be folded")

    n_roots = last(size(ket_coeffs))
    ket_coeffs2 = permutedims(ket_coeffs, perm)
    bra_coeffs2 = permutedims(bra_coeffs, perm)

    dim1 = size(ket_coeffs2)
    ket_coeffs2 = reshape(ket_coeffs2, dim1[1]*dim1[2], prod(dim1[3:end]))

    dim2 = size(bra_coeffs2)
    bra_coeffs2 = reshape(bra_coeffs2, dim2[1]*dim2[2], prod(dim2[3:end]))

    op = reshape(op, prod(size(op)[1:2]),prod(size(op)[3:4]))

#    println()
#    display((c1.idx, c2.idx))
#    display(perm')
#    display(size(op))
#    display(size(ket_coeffs))
#    display(size(permutedims(ket_coeffs,perm)))
#    display(size(ket_coeffs2))
    if state_sign == 1
        bra_coeffs2 .+= op' * ket_coeffs2
    elseif state_sign == -1
        bra_coeffs2 .-= op' * ket_coeffs2
    else
        error()
    end

    ket_coeffs2 = reshape(ket_coeffs2, dim1)
    bra_coeffs2 = reshape(bra_coeffs2, dim2)

    # now untranspose
    perm,_ = bubble_sort(perm)
    ket_coeffs2 = permutedims(ket_coeffs2,perm)
    bra_coeffs2 = permutedims(bra_coeffs2,perm)

    bra_coeffs .= bra_coeffs2
    return
#=}}}=#
end
"""
"""
function form_sigma_block!(term::ClusteredTerm3B,
                            cluster_ops::Vector{ClusterOps},
                            fock_bra::FockConfig, bra::TuckerConfig,
                            fock_ket::FockConfig, ket::TuckerConfig,
                            bra_coeffs::Array, ket_coeffs::Array)
#={{{=#
    #display(term)
    #println(bra, ket)

    c1 = term.clusters[1]
    c2 = term.clusters[2]
    c3 = term.clusters[3]
    length(fock_bra) == length(fock_ket) || throw(Exception)
    length(bra) == length(ket) || throw(Exception)
    n_clusters = length(bra)
    #
    # make sure inactive clusters are diagonal
    for ci in 1:n_clusters
        ci != c1.idx || continue
        ci != c2.idx || continue
        ci != c3.idx || continue

        fock_bra[ci] == fock_ket[ci] || throw(Exception)
        bra[ci] == ket[ci] || return 0.0
    end

    #
    # make sure active clusters are correct transitions
    fock_bra[c1.idx] == fock_ket[c1.idx] .+ term.delta[1] || throw(Exception)
    fock_bra[c2.idx] == fock_ket[c2.idx] .+ term.delta[2] || throw(Exception)
    fock_bra[c3.idx] == fock_ket[c3.idx] .+ term.delta[3] || throw(Exception)

    #
    # determine sign from rearranging clusters if odd number of operators
    state_sign = compute_terms_state_sign(term, fock_ket)

    #
    # op[IKM,JLN] = <I|p'|J> h(pqr) <K|q|L> <M|r|N>
    gamma1 = cluster_ops[c1.idx][term.ops[1]][(fock_bra[c1.idx],fock_ket[c1.idx])][:,bra[c1.idx],ket[c1.idx]]
    gamma2 = cluster_ops[c2.idx][term.ops[2]][(fock_bra[c2.idx],fock_ket[c2.idx])][:,bra[c2.idx],ket[c2.idx]]
    gamma3 = cluster_ops[c3.idx][term.ops[3]][(fock_bra[c3.idx],fock_ket[c3.idx])][:,bra[c3.idx],ket[c3.idx]]

    op = Array{Float64}[]
    #@tensor begin
    #    op[J,L,N,I,K,M] := term.ints[p,q,r] * gamma1[p,I,J] * gamma2[q,K,L] * gamma3[r,M,N]
    #end
    cache_key = (fock_bra[c1.idx], fock_bra[c2.idx], fock_bra[c3.idx],
                 fock_ket[c1.idx], fock_ket[c2.idx], fock_ket[c3.idx],
                 bra[c1.idx], bra[c2.idx], bra[c3.idx],
                 ket[c1.idx], ket[c2.idx], ket[c3.idx])

#    if haskey(term.cache, cache_key)
#        op = term.cache[cache_key]
#    else
#        @tensor begin
#            op[q,r,I,J] := term.ints[p,q,r] * gamma1[p,I,J]
#            op[r,I,J,K,L] := op[q,r,I,J] * gamma2[q,K,L]
#            op[J,L,N,I,K,M] := op[r,I,J,K,L] * gamma3[r,M,N]
#        end
#        term.cache[cache_key] = op
#
#        #tucker_decompose(op)
#        # compress this
#        opsize = size(op)
#        op = reshape(op, prod(size(op)[1:3]), prod(size(op)[4:6]))
#        F = svd(op)
#        #display(F.S)
#        cut = 0
#        for si in 1:length(F.S)
#            if F.S[si] < 1e-4
#                F.S[si] = 0
#                cut += 1
#            end
#        end
#        #if cut > 0
#        #    display((length(F.S), cut))
#        #end
#        #op = F.U * Diagonal(F.S) * F.Vt
#        op = reshape(op,opsize)
#        core, factors = tucker_decompose(op, thresh=-1, verbose=0)
#        op = tucker_recompose(core, factors)
#    end

    @tensor begin
        op[q,r,I,J] := term.ints[p,q,r] * gamma1[p,I,J]
        op[r,I,J,K,L] := op[q,r,I,J] * gamma2[q,K,L]
        op[J,L,N,I,K,M] := op[r,I,J,K,L] * gamma3[r,M,N]
    end

    # now transpose state vectors and multiply, also, try without transposing to compare
    indices = collect(1:n_clusters+1)
    indices[c1.idx] = 0
    indices[c2.idx] = 0
    indices[c3.idx] = 0
    perm,_ = bubble_sort(indices)

    length(size(ket_coeffs)) == n_clusters + 1 || error(" tensors should be folded")

    n_roots = last(size(ket_coeffs))
    ket_coeffs2 = permutedims(ket_coeffs,perm)
    bra_coeffs2 = permutedims(bra_coeffs,perm)

    dim1 = size(ket_coeffs2)
    ket_coeffs2 = reshape(ket_coeffs2, dim1[1]*dim1[2]*dim1[3], prod(dim1[4:end]))

    dim2 = size(bra_coeffs2)
    bra_coeffs2 = reshape(bra_coeffs2, dim2[1]*dim2[2]*dim2[3], prod(dim2[4:end]))

    op = reshape(op, prod(size(op)[1:3]),prod(size(op)[4:6]))
    if state_sign == 1
        bra_coeffs2 .+= op' * ket_coeffs2
    elseif state_sign == -1
        bra_coeffs2 .-= op' * ket_coeffs2
    else
        error()
    end

    ket_coeffs2 = reshape(ket_coeffs2, dim1)
    bra_coeffs2 = reshape(bra_coeffs2, dim2)

    # now untranspose
    perm,_ = bubble_sort(perm)
    ket_coeffs2 = permutedims(ket_coeffs2,perm)
    bra_coeffs2 = permutedims(bra_coeffs2,perm)

    bra_coeffs .= bra_coeffs2
    return
#=}}}=#
end
"""
"""
function form_sigma_block!(term::ClusteredTerm4B,
                            cluster_ops::Vector{ClusterOps},
                            fock_bra::FockConfig, bra::TuckerConfig,
                            fock_ket::FockConfig, ket::TuckerConfig,
                            bra_coeffs, ket_coeffs)
#={{{=#
    #display(term)
    #println(bra, ket)

    c1 = term.clusters[1]
    c2 = term.clusters[2]
    c3 = term.clusters[3]
    c4 = term.clusters[4]
    length(fock_bra) == length(fock_ket) || throw(Exception)
    length(bra) == length(ket) || throw(Exception)
    n_clusters = length(bra)
    #
    # make sure inactive clusters are diagonal
    for ci in 1:n_clusters
        ci != c1.idx || continue
        ci != c2.idx || continue
        ci != c3.idx || continue
        ci != c4.idx || continue

        fock_bra[ci] == fock_ket[ci] || throw(Exception)
        bra[ci] == ket[ci] || return 0.0
    end

    #
    # make sure active clusters are correct transitions
    fock_bra[c1.idx] == fock_ket[c1.idx] .+ term.delta[1] || throw(Exception)
    fock_bra[c2.idx] == fock_ket[c2.idx] .+ term.delta[2] || throw(Exception)
    fock_bra[c3.idx] == fock_ket[c3.idx] .+ term.delta[3] || throw(Exception)
    fock_bra[c4.idx] == fock_ket[c4.idx] .+ term.delta[4] || throw(Exception)

    #
    # determine sign from rearranging clusters if odd number of operators
    state_sign = compute_terms_state_sign(term, fock_ket)

    #
    # op[IKMO,JLNP] = <I|p'|J> h(pqrs) <K|q|L> <M|r|N> <O|s|P>
    gamma1 = cluster_ops[c1.idx][term.ops[1]][(fock_bra[c1.idx],fock_ket[c1.idx])][:,bra[c1.idx],ket[c1.idx]]
    gamma2 = cluster_ops[c2.idx][term.ops[2]][(fock_bra[c2.idx],fock_ket[c2.idx])][:,bra[c2.idx],ket[c2.idx]]
    gamma3 = cluster_ops[c3.idx][term.ops[3]][(fock_bra[c3.idx],fock_ket[c3.idx])][:,bra[c3.idx],ket[c3.idx]]
    gamma4 = cluster_ops[c4.idx][term.ops[4]][(fock_bra[c4.idx],fock_ket[c4.idx])][:,bra[c4.idx],ket[c4.idx]]

    op = Array{Float64}[]
    @tensor begin
        op[J,L,N,P,I,K,M,O] := term.ints[p,q,r,s] * gamma1[p,I,J] * gamma2[q,K,L] * gamma3[r,M,N] * gamma4[s,O,P]
    end
    #@tensor begin
    #    op[q,r,I,J] := term.ints[p,q,r] * gamma1[p,I,J]
    #    op[r,I,J,K,L] := op[q,r,I,J] * gamma2[q,K,L]
    #    op[J,L,N,I,K,M] := op[r,I,J,K,L] * gamma2[r,M,N]
    #end

    # possibly cache some of these integrals
    # compress this
#    opsize = size(op)
#    op = reshape(op, prod(size(op)[1:4]), prod(size(op)[5:8]))
#    F = svd(op)
#    #display(F.S)
#    for si in 1:length(F.S)
#        if F.S[si] < 1e-3
#            F.S[si] = 0
#        end
#    end
#    op = F.U * Diagonal(F.S) * F.Vt
#    op = reshape(op,opsize)

    # now transpose state vectors and multiply, also, try without transposing to compare
    indices = collect(1:n_clusters+1)
    indices[c1.idx] = 0
    indices[c2.idx] = 0
    indices[c3.idx] = 0
    indices[c4.idx] = 0
    perm,_ = bubble_sort(indices)

    length(size(ket_coeffs)) == n_clusters + 1 || error(" tensors should be folded")

    n_roots = last(size(ket_coeffs))
    ket_coeffs2 = permutedims(ket_coeffs,perm)
    bra_coeffs2 = permutedims(bra_coeffs,perm)

    dim1 = size(ket_coeffs2)
    ket_coeffs2 = reshape(ket_coeffs2, dim1[1]*dim1[2]*dim1[3]*dim1[4], prod(dim1[5:end]))

    dim2 = size(bra_coeffs2)
    bra_coeffs2 = reshape(bra_coeffs2, dim2[1]*dim2[2]*dim2[3]*dim2[4], prod(dim2[5:end]))

    op = reshape(op, prod(size(op)[1:4]),prod(size(op)[5:8]))
    if state_sign == 1
        bra_coeffs2 .+= op' * ket_coeffs2
    elseif state_sign == -1
        bra_coeffs2 .-= op' * ket_coeffs2
    else
        error()
    end

    ket_coeffs2 = reshape(ket_coeffs2, dim1)
    bra_coeffs2 = reshape(bra_coeffs2, dim2)

    # now untranspose
    perm,_ = bubble_sort(perm)
    ket_coeffs2 = permutedims(ket_coeffs2,perm)
    bra_coeffs2 = permutedims(bra_coeffs2,perm)

    bra_coeffs .= bra_coeffs2
    return
#=}}}=#
end

"""
    build_sigma!(sigma_vector::TuckerState, ci_vector::TuckerState, cluster_ops, clustered_ham, nbody=4)
"""
function build_sigma!(sigma_vector::TuckerState, ci_vector::TuckerState, cluster_ops, clustered_ham; nbody=4)
    #={{{=#

    for (fock_bra, configs_bra) in sigma_vector
        for (fock_ket, configs_ket) in ci_vector
            fock_trans = fock_bra - fock_ket

            # check if transition is connected by H
            haskey(clustered_ham, fock_trans) == true || continue

            for (config_bra, coeff_bra) in configs_bra
                for (config_ket, coeff_ket) in configs_ket

                    for term in clustered_ham[fock_trans]

                        length(term.clusters) <= nbody || continue

                        FermiCG.form_sigma_block!(term, cluster_ops, fock_bra, config_bra,
                                                  fock_ket, config_ket,
                                                  coeff_bra, coeff_ket)

                    end
                end
            end
        end
    end
    return
    #=}}}=#
end

"""
    get_map(ci_vector, cluster_ops, clustered_ham)

Get LinearMap with takes a vector and returns action of H on that vector
"""
function get_map(ci_vector::TuckerState, cluster_ops, clustered_ham; shift = nothing)
    #={{{=#
    iters = 0

    dim = length(ci_vector)
    function mymatvec(v)
        iters += 1

        nr = 0
        if length(size(v)) == 1
            nr = 1
            v = reshape(v, length(v), nr)
        elseif length(size(v)) == 2
            nr = size(v)[2]
        else
            error(" is tensor not unfolded?")
        end

        set_vector!(ci_vector, v)

        fold!(ci_vector)
        sig = deepcopy(ci_vector)
        zero!(sig)
        build_sigma!(sig, ci_vector, cluster_ops, clustered_ham)

        unfold!(ci_vector)

        sig = get_vector(sig)

        if shift != nothing
            # this is how we do CEPA
            sig += shift * get_vector(ci_vector)
        end

        return sig
    end
    return LinearMap(mymatvec, dim, dim; issymmetric=true, ismutating=false, ishermitian=true)
end
#=}}}=#

function tucker_ci_solve!(ci_vector::TuckerState, cluster_ops, clustered_ham; tol=1e-5)
#={{{=#
    unfold!(ci_vector)
    Hmap = get_map(ci_vector, cluster_ops, clustered_ham)

    v0 = get_vector(ci_vector)
    nr = size(v0)[2]

    davidson = Davidson(Hmap,v0=v0,max_iter=80, max_ss_vecs=40, nroots=nr, tol=1e-5)
    #Adiag = StringCI.compute_fock_diagonal(problem,mf.mo_energy, e_mf)
    #FermiCG.solve(davidson)
    @printf(" Now iterate: \n")
    flush(stdout)
    #@time FermiCG.iteration(davidson, Adiag=Adiag, iprint=2)
    e,v = FermiCG.solve(davidson)
    set_vector!(ci_vector,v)
    return e,v
end
#=}}}=#

function compress_blocks(ts::TuckerState; thresh=1e-7)
    for (fock, tconfigs) in ts.data
        display(fock)
        for (tconfig, coeffs) in tconfigs
            display(tconfig)
            tucker_decompose(coeffs)
        end
    end
end

max_stars_repo_name: nmayhall-vt/FermiCG

max_issues_repo_name: nmayhall-vt/FermiCG

max_forks_repo_name: nmayhall-vt/FermiCG


text:
module Multimedia

export Display, display, pushdisplay, popdisplay, displayable, redisplay,
   MIME, @MIME, @MIME_str, writemime, reprmime, stringmime, istext,
   mimewritable, TextDisplay

###########################################################################
# We define a singleton type MIME{mime symbol} for each MIME type, so
# that Julia's dispatch and overloading mechanisms can be used to
# dispatch writemime and to add conversions for new types.

immutable MIME{mime} end

import Base: show, print, string, convert
MIME(s) = MIME{symbol(s)}()
show{mime}(io::IO, ::MIME{mime}) = print(io, "MIME type ", string(mime))
print{mime}(io::IO, ::MIME{mime}) = print(io, mime)

# needs to be a macro so that we can use ::@mime(s) in type declarations
macro MIME(s)
    Base.warn_once("@MIME(\"\") is deprecated, use MIME\"\" instead.")
    if isa(s,String)
        :(MIME{$(Expr(:quote, symbol(s)))})
    else
        :(MIME{symbol($s)})
    end
end

macro MIME_str(s)
    :(MIME{$(Expr(:quote, symbol(s)))})
end

###########################################################################
# For any type T one can define writemime(io, ::MIME"type", x::T) = ...
# in order to provide a way to export T as a given mime type.

mimewritable{mime}(::MIME{mime}, x) =
  method_exists(writemime, (IO, MIME{mime}, typeof(x)))

# it is convenient to accept strings instead of ::MIME
writemime(io, m::String, x) = writemime(io, MIME(m), x)
mimewritable(m::String, x) = mimewritable(MIME(m), x)

###########################################################################
# MIME types are assumed to be binary data except for a set of types known
# to be text data (possibly Unicode).  istext(m) returns whether
# m::MIME is text data, and reprmime(m, x) returns x written to either
# a string (for text m::MIME) or a Vector{Uint8} (for binary m::MIME),
# assuming the corresponding write_mime method exists.  stringmime
# is like reprmime except that it always returns a string, which in the
# case of binary data is Base64-encoded.
#
# Also, if reprmime is passed a String for a text type or Vector{Uint8} for
# a binary type, the argument is assumed to already be in the corresponding
# format and is returned unmodified.  This is useful so that raw data can be
# passed to display(m::MIME, x).

macro textmime(mime)
    quote
        mimeT = MIME{symbol($mime)}
        # avoid method ambiguities with the general definitions below:
        # (Q: should we treat Vector{Uint8} as a bytestring?)
        Base.Multimedia.reprmime(m::mimeT, x::Vector{Uint8}) = sprint(writemime, m, x)
        Base.Multimedia.stringmime(m::mimeT, x::Vector{Uint8}) = reprmime(m, x)

        Base.Multimedia.istext(::mimeT) = true
        if $(mime != "text/plain") # strings are shown escaped for text/plain
            Base.Multimedia.reprmime(m::mimeT, x::String) = x
        end
        Base.Multimedia.reprmime(m::mimeT, x) = sprint(writemime, m, x)
        Base.Multimedia.stringmime(m::mimeT, x) = reprmime(m, x)
    end
end

istext(::MIME) = false
function reprmime(m::MIME, x)
    s = IOBuffer()
    writemime(s, m, x)
    takebuf_array(s)
end
reprmime(m::MIME, x::Vector{Uint8}) = x
stringmime(m::MIME, x) = base64(writemime, m, x)
stringmime(m::MIME, x::Vector{Uint8}) = base64(write, x)

# it is convenient to accept strings instead of ::MIME
istext(m::String) = istext(MIME(m))
reprmime(m::String, x) = reprmime(MIME(m), x)
stringmime(m::String, x) = stringmime(MIME(m), x)

for mime in ["text/vnd.graphviz", "text/latex", "text/calendar", "text/n3", "text/richtext", "text/x-setext", "text/sgml", "text/tab-separated-values", "text/x-vcalendar", "text/x-vcard", "text/cmd", "text/css", "text/csv", "text/html", "text/javascript", "text/plain", "text/vcard", "text/xml", "application/atom+xml", "application/ecmascript", "application/json", "application/rdf+xml", "application/rss+xml", "application/xml-dtd", "application/postscript", "image/svg+xml", "application/x-latex", "application/xhtml+xml", "application/javascript", "application/xml", "model/x3d+xml", "model/x3d+vrml", "model/vrml"]
    @eval @textmime $mime
end

###########################################################################
# We have an abstract Display class that can be subclassed in order to
# define new rich-display output devices.  A typical subclass should
# overload display(d::Display, m::MIME, x) for supported MIME types m,
# (typically using reprmime or stringmime to get the MIME
# representation of x) and should also overload display(d::Display, x)
# to display x in whatever MIME type is preferred by the Display and
# is writable by x.  display(..., x) should throw a MethodError if x
# cannot be displayed.  The return value of display(...) is up to the
# Display type.

abstract Display

# it is convenient to accept strings instead of ::MIME
display(d::Display, mime::String, x) = display(d, MIME(mime), x)
display(mime::String, x) = display(MIME(mime), x)
displayable(d::Display, mime::String) = displayable(d, MIME(mime))
displayable(mime::String) = displayable(MIME(mime))

# simplest display, which only knows how to display text/plain
immutable TextDisplay <: Display
    io::IO
end
display(d::TextDisplay, M::MIME"text/plain", x) = writemime(d.io, M, x)
display(d::TextDisplay, x) = display(d, MIME"text/plain"(), x)

import Base: close, flush
flush(d::TextDisplay) = flush(d.io)
close(d::TextDisplay) = close(d.io)

###########################################################################
# We keep a stack of Displays, and calling display(x) uses the topmost
# Display that is capable of displaying x (doesn't throw an error)

const displays = Display[]
function pushdisplay(d::Display)
    global displays
    push!(displays, d)
end
popdisplay() = pop!(displays)
function popdisplay(d::Display)
    for i = length(displays):-1:1
        if d == displays[i]
            return splice!(displays, i)
        end
    end
    throw(KeyError(d))
end
function reinit_displays()
    empty!(displays)
    pushdisplay(TextDisplay(STDOUT))
end

function display(x)
    for i = length(displays):-1:1
        try
            return display(displays[i], x)
        catch e
            if !isa(e, MethodError)
                rethrow()
            end
        end
    end
    throw(MethodError(display, (x,)))
end

function display(m::MIME, x)
    for i = length(displays):-1:1
        try
            return display(displays[i], m, x)
        catch e
            if !isa(e, MethodError)
                rethrow()
            end
        end
    end
    throw(MethodError(display, (m, x)))
end

displayable{D<:Display,mime}(d::D, ::MIME{mime}) =
  method_exists(display, (D, MIME{mime}, Any))

function displayable(m::MIME)
    for d in displays
        if displayable(d, m)
            return true
        end
    end
    return false
end

###########################################################################
# The redisplay method can be overridden by a Display in order to
# update an existing display (instead of, for example, opening a new
# window), and is used by the IJulia interface to defer display
# until the next interactive prompt.  This is especially useful
# for Matlab/Pylab-like stateful plotting interfaces, where
# a plot is created and then modified many times (xlabel, title, etc.).

function redisplay(x)
    for i = length(displays):-1:1
        try
            return redisplay(displays[i], x)
        catch e
            if !isa(e, MethodError)
                rethrow()
            end
        end
    end
    throw(MethodError(redisplay, (x,)))
end

function redisplay(m::Union(MIME,String), x)
    for i = length(displays):-1:1
        try
            return redisplay(displays[i], m, x)
        catch e
            if !isa(e, MethodError)
                rethrow()
            end
        end
    end
    throw(MethodError(redisplay, (m, x)))
end

# default redisplay is simply to call display
redisplay(d::Display, x) = display(d, x)
redisplay(d::Display, m::Union(MIME,String), x) = display(d, m, x)

###########################################################################

end # module

max_stars_repo_name: alyst/julia

max_issues_repo_name: alyst/julia

max_forks_repo_name: alyst/julia


text:
using Gnuplot, Random
Random.seed!(145)
x, y = 1:2:20, 5*rand(10)
xerr = abs.(randn(10))
lcval = 1:10
@gp x y xerr lcval "w xerr t 'test' lc palette" palette(:berlin)
@gp :- xrange= (0,22) yrange = (0,6)
save(term="pngcairo size 600,400", output="./code/plt_errs19.png") # hide

max_stars_repo_name: lazarusA/gnuplot-examples

max_issues_repo_name: lazarusA/gnuplot-examples

max_forks_repo_name: lazarusA/gnuplot-examples


text:
export Contractor
export Bisection, Newton, Krawczyk

"""
    𝒩(f, f′, X, α)

Single-variable Newton operator.

The symbol for the operator is accessed with `\\scrN<tab>`.
"""
function 𝒩(f, f′, X::Interval{T}, α) where {T}
    m = Interval(mid(X, α))

    return m - (f(m) / f′(X))
end

"""
    𝒩(f, jacobian, X, α)

Multi-variable Newton operator.
"""
function 𝒩(f::Function, jacobian::Function, X::IntervalBox, α)  # multidimensional Newton operator
    m = IntervalBox(Interval.(mid(X, α)))
    J = jacobian(X)

    return IntervalBox(m .- (J \ f(m)))
end

"""
    𝒦(f, f′, X, α)

Single-variable Krawczyk operator.

The symbol for the operator is accessed with `\\scrK<tab>`.
"""
function 𝒦(f, f′, X::Interval{T}, α) where {T}
    m = Interval(mid(X, α))
    Y = 1 / f′(m)

    return m - Y*f(m) + (1 - Y*f′(X)) * (X - m)
end

"""
    𝒦(f, jacobian, X, α)

Multi-variable Krawczyk operator.
"""
function 𝒦(f, jacobian, X::IntervalBox{T}, α) where {T}
    m = mid(X, α)
    mm = IntervalBox(m)
    J = jacobian(X)
    Y = mid.(inv(jacobian(mm)))

    return m - Y*f(mm) + (I - Y*J) * (X.v - m)
end

"""
    Contractor{F}

Abstract type for contractors.
"""
abstract type Contractor{F} end

"""
    Bisection{F} <: Contractor{F}

Contractor type for the bisection method.
"""
struct Bisection{F} <: Contractor{F}
    f::F
end

function (contractor::Bisection)(r, tol)
    X = interval(r)
    image = (contractor.f)(X)

    if root_status(r) == :empty || !(contains_zero(image))
        return Root(X, :empty)
    end

    return Root(X, :unknown)
end

for (Method, 𝒪) in ((:Newton, 𝒩), (:Krawczyk, 𝒦))
    doc = """
        $Method{F, FP} <: Contractor{F}

    Contractor type for the $Method method.

    # Fields
      - `f::F`: function whose roots are searched
      - `f::FP`: derivative or jacobian of `f`

    -----

        (C::$Method)(X, tol, α=where_bisect)

    Contract an interval `X` using $Method operator and return the
    contracted interval together with its status.

    # Inputs
      - `R`: Root object containing the interval to contract.
      - `tol`: Precision to which unique solutions are refined.
      - `α`: Point of bisection of intervals.
    """

    @eval begin
        struct $Method{F, FP} <: Contractor{F}
            f::F
            f′::FP   # use \prime<TAB> for ′
        end

        function (C::$Method)(R, tol, α=where_bisect)
            op = x -> $𝒪(C.f, C.f′, x, α)
            rt = determine_region_status(op, C.f, R)
            return refine(op, rt, tol)
        end

        @doc $doc $Method
    end
end

"""
    safe_isempty(X)

Similar to `isempty` function for `IntervalBox`, but also works for `SVector`
of `Interval`.
"""
safe_isempty(X) = isempty(IntervalBox(X))

"""
    determine_region_status(contract, f, R)

Contraction operation for contractors using the first derivative of the
function.

Currently `Newton` and `Krawczyk` contractors use this.
"""
function determine_region_status(op, f, R)
    X = interval(R)
    former_status = root_status(R)

    imX = f(X)

    if former_status == :empty || !(contains_zero(imX))
        return Root(X, :empty)
    end

    safe_isempty(imX) && return Root(X, :empty)  # X is fully outside of the domain of f

    contracted_X = op(X)

    # Only happens if X is partially out of the domain of f
    safe_isempty(contracted_X) && return Root(X, :unknown)  # force bisection

    # given that have the Jacobian, can also do mean value form
    NX = contracted_X ∩ X

    isinf(X) && return Root(NX, :unknown)  # force bisection
    safe_isempty(NX) && return Root(X, :empty)

    if former_status == :unique || NX ⪽ X  # isinterior; know there's a unique root inside
        return Root(NX, :unique)
    end

    return Root(NX, :unknown)
end

"""
    refine(op, X::Region, tol)

Generic refine operation for Krawczyk and Newton.
This function assumes that it is already known that `X` contains a unique root.
"""
function refine(op, X::Region, tol)
    while diam(X) > tol  # avoid problem with tiny floating-point numbers if 0 is a root
        NX = op(X) ∩ X
        NX == X && break  # reached limit of precision
        X = NX
    end

    return X
end

"""
    refine(op, X::Tuple{Symbol, Region}, tol)

Wrap the refine method to leave unchanged intervals that are not guaranteed to
contain an unique solution.
"""
function refine(op, R::Root, tol)
    root_status(R) != :unique && return R
    return Root(refine(op, interval(R), tol), :unique)
end


text:
# Initialize 2D matrix w/rows corresponding to monomers
# and columns corresponding to ligands
# Indices are as follows:
# c[1,1] = M+, c[3,2] = M, c[2,3] = L, c[3,3] = ML
# Julia has 1-based indices therefore # of monomers = row index + 1 and
# # of ligands = column index + 1
function c0(imax, Ns, Mplus, L)

rows    = imax + 3
columns = Ns[imax] + 3

c0 = zeros(Float64, rows, columns)
c0[1,1] = Mplus
c0[2,3] = L

return c0

end


text:
const ALEInterface = Cvoid  # Using Cvoid inplace of Nothing because it is better suited to the context
const ALEPtr = Ptr{ALEInterface}
const ALEState = Cvoid
const ALEStatePtr = Ptr{ALEState}

const Logger = Dict(
    :info => 0,
    :warning => 1,
    :error => 2
)

"""
    ALE_new()

Creates a new ALE instance.
"""
ALE_new() = ccall((:ALE_new, libale_c), ALEPtr, ())

"""
    ALE_del(ale_instance::ALEPtr)

Destroys the passed ALE instance pointer. Acts as a destructor. Use with caution.
"""
ALE_del(ale::ALEPtr) = ccall((:ALE_del, libale_c), Cvoid, (ALEPtr,), ale)

"""
    getString(ale_instance::ALEPtr, key::String)

Returns the string value held by the ALE instance's field with the same name as `key`.
"""
function getString(ale::ALEPtr, key::String)
    res = ccall((:getString, libale_c), Ptr{Cchar}, (ALEPtr, Ptr{Cchar}),
        ale, key)
    bytestring(res)
end

"""
    getInt(ale_instance::ALEPtr, key::String)

Returns the integer value held by the ALE instance's field with the same name as `key`.
"""
getInt(ale::ALEPtr, key::String) = ccall((:getInt, libale_c), Cint,
    (ALEPtr, Ptr{Cchar}), ale, key)

"""
    getBool(ale_instance::ALEPtr, key::String)

Returns the boolean value held by the ALE instance's field with the same name as `key`.
"""
getBool(ale::ALEPtr, key::String) = ccall((:getBool, libale_c), Cint,
    (ALEPtr, Ptr{Cchar}), ale, key) > 0

"""
    getFloat(ale_instance::ALEPtr, key::String)

Returns the float value held by the ALE field with the same name as `key`.
"""
getFloat(ale::ALEPtr, key::String) = ccall((:getFloat, libale_c), Float32,
    (ALEPtr, Ptr{Cchar}), ale, key)

"""
    setString(ale_interface::ALEPtr, key::String, value::String)

Modifies the ALE instance field value with the same name as `key` to the passed String value.
"""
setString(ale::ALEPtr, key::String, value::String) =
    ccall((:setString, libale_c), Cvoid, (ALEPtr, Ptr{Cchar}, Ptr{Cchar}),
        ale, key, value)

"""
    setInt(ale_instance::ALEPtr, key::String, value::Int32)

Modifies the ALE instance field value with the same name as `key` to the passed integer value.
"""
setInt(ale::ALEPtr, key::String, value::Integer) = ccall((:setInt, libale_c),
    Cvoid, (ALEPtr, Ptr{Cchar}, Cint), ale, key, Cint(value))

"""
    setBool(ale_instance::ALEPtr, key::String, value::Bool)

Modifies the ALE instance field value with the same name as `key` to the passed boolean value.
"""
setBool(ale::ALEPtr, key::String, value::Bool) = ccall((:setBool, libale_c),
    Cvoid, (ALEPtr, Ptr{Cchar}, Cint), ale, key, value)

"""
    setFloat(ale_instance::ALEPtr, key::String, value::)

Modifies the ALE instance field with the same name as `key` to the passed floating point value.
"""
setFloat(ale::ALEPtr, key::String, value::Real) =
    ccall((:setFloat, libale_c), Cvoid, (ALEPtr, Ptr{Cchar}, Cfloat),
        ale, key, Cfloat(value))

"""
    loadROM(ale_instance::ALEPtr, rom_file::String)

Loads the binary of passed. `rom_file` can either be the absolute path to the binary,
or the name of the ROM that is present in the "deps/roms" directory. Access this list using
`getROMList()`.

# Examples
```julia-repl
julia> loadROM(ale, "pong")
Game console created:
  ROM file:  /Users/juliauser/.julia/artifacts/4af00c03a4bcddb3ce20c2d96c3d09527100767/ROMS/pong.bin
  Cart Name: Video Olympics (1978) (Atari)
  Cart MD5:  60e0ea3cbe0913d39803477945e9e5ec
  Display Format:  AUTO-DETECT ==> NTSC
  ROM Size:        2048
  Bankswitch Type: AUTO-DETECT ==> 2K

WARNING: Possibly unsupported ROM: mismatched MD5.
Cartridge_MD5: 60e0ea3cbe0913d39803477945e9e5ec
Cartridge_name: Video Olympics (1978) (Atari)

Running ROM file...
Random seed is 0

julia> loadROM(ale, "/Users/juilauser/Desktop/pewpew/roms/ms_pacman.bin")
Game console created:
  ROM file:  /Users/juilauser/Desktop/pewpew/roms/ms_pacman.bin
  Cart Name: Ms. Pac-Man (1982) (CCE)
  Cart MD5:  9469d18238345d87768e8965f9f4a6b2
  Display Format:  AUTO-DETECT ==> NTSC
  ROM Size:        8192
  Bankswitch Type: AUTO-DETECT ==> F8

WARNING: Possibly unsupported ROM: mismatched MD5.
Cartridge_MD5: 9469d18238345d87768e8965f9f4a6b2
Cartridge_name: Ms. Pac-Man (1982) (CCE)

Running ROM file...
Random seed is 0
```
"""
function loadROM(ale::ALEPtr, rom_file::String)
    if isfile(rom_file)
        ccall((:loadROM, libale_c), Nothing, (ALEPtr, Ptr{Cchar}), ale, rom_file)
    elseif isfile(joinpath(ROM_PATH, rom_file * ".bin"))
        loadROM(ale, joinpath(ROM_PATH, rom_file * ".bin"))
    else
        @error("ROM file $rom_file not found.")
    end
end

"""
    getROMList()

Returns an array of names of ROMs available. These ROMS can be loaded by simply passing
their name to the [`loadROM`](@ref) function, as opposed to passing their full path.

# Example
```julia-repl
julia> getROMList()
74-element Array{String,1}:
 "adventure"
 "air_raid"
 "alien"
 "amidar"
 "assault"
 "asterix"
 ⋮
 "venture"
 "video_pinball"
 "wizard_of_wor"
 "yars_revenge"
 "zaxxon"
```
"""
function getROMList()
    files = readdir(ROM_PATH)
    files = files[endswith.(files, ".bin")]
    file_names = [file_name[1:end-4] for file_name in files]
    file_names
end

"""
    act(ale_instance::ALEPtr, action::Integer)

Applies an action to the game and returns the reward. It is the user's
responsibility to check if the game has ended and reset when neccessary -
this function will keep pushing buttons on the game over screen.
"""
act(ale::ALEPtr, action::Integer)::Int =
    ccall((:act, libale_c), Cint, (ALEPtr, Cint), ale, Cint(action))

"""
    game_over(ale_instance::ALEPtr)

Returns a boolean value that tells whether the game is over. Any actions performed
from hereon until the next `reset_game()` have no effect.
"""
game_over(ale::ALEPtr) =
    ccall((:game_over, libale_c), Cint, (ALEPtr,), ale) > 0

"""
    reset_game(ale_instance::ALEPtr)

Resets the game and the reward counter.
"""
reset_game(ale::ALEPtr) = ccall((:reset_game, libale_c), Cvoid, (ALEPtr,), ale)

"""
    getLegalActionSet(ale_instance::ALEPtr)

Returns the vector of the legal set of actions needed the play the game. Not to be confused with `getMinimalActionSet`.

See also: [`getMinimalActionSet`](@ref)
"""
function getLegalActionSet(ale::ALEPtr)
    actions = Cint[]
    _getLegalActionSet!(ale, actions)
    Int.(actions)
end

function _getLegalActionSet!(ale::ALEPtr, actions::Vector{Cint})
    resize!(actions, getLegalActionSize(ale))
    ccall((:getLegalActionSet, libale_c), Cvoid, (ALEPtr, Ptr{Cint}),
        ale, actions)
end

"""
    getLegalActionSize(ale_instance::ALEPtr)

Returns the size of the legal action set. Not to be confused with `getMinimalActionSet`.

See also: [`getLegalActionSize`](@ref)
"""
getLegalActionSize(ale::ALEPtr)::Int =
    ccall((:getLegalActionSize, libale_c), Cint, (ALEPtr,), ale)

"""
    getMinimalActionSet(ale_instance::ALEPtr)

Returns the set of actions that actually affect the game. Not to be confused with `getLegalActionSet`.

See also: [`getLegalActionSet`](@ref)
"""
function getMinimalActionSet(ale::ALEPtr)
    actions = Cint[]
    _getMinimalActionSet!(ale, actions)
    Int.(actions)
end

function _getMinimalActionSet!(ale::ALEPtr, actions::Vector{Cint})
    resize!(actions, getMinimalActionSize(ale))
    ccall((:getMinimalActionSet, libale_c), Cvoid, (ALEPtr, Ptr{Cint}),
        ale, actions)
end

"""
    getMinimalActionSize(ale_instance::ALEPtr)

Returns the size of the minimal action set required to play the game. Not to be confused with `getLegalActionSize`.

See also: [`getLegalActionSize`](@ref)
"""
getMinimalActionSize(ale::ALEPtr)::Int =
    ccall((:getMinimalActionSize, libale_c), Cint, (ALEPtr,), ale)

"""
    getFrameNumber(ale_instance::ALEPtr)

Returns the frame number since the ALE instance was created.
"""
getFrameNumber(ale::ALEPtr) =
    ccall((:getFrameNumber, libale_c), Cint, (ALEPtr,), ale)

"""
    lives(ale_instance::ALEPtr)

Returns the number of lives remaining, by games that support the concept of lives.
"""
lives(ale::ALEPtr) = ccall((:lives, libale_c), Cint, (ALEPtr,), ale)

"""
    getEpisodeFrameNumber(ale_instance::ALEPtr)

Returns the frame number since the last reset.
"""
getEpisodeFrameNumber(ale::ALEPtr) =
    ccall((:getEpisodeFrameNumber, libale_c), Cint, (ALEPtr,), ale)

"""
    getRAMSize(ale_instance::ALEPtr)

Returns the size of the RAM.
"""
getRAMSize(ale::ALEPtr) = ccall((:getRAMSize, libale_c), Cint, (ALEPtr,), ale)
_getRAM(ale::ALEPtr, ram::Vector{Cuchar}) =
    ccall((:getRAM, libale_c), Cvoid, (ALEPtr, Ptr{Cuchar}), ale, ram)

"""
    getRAM(ale_instance::ALEPtr)

Returns the snapshot of the RAM as a vector.
"""
function getRAM(ale::ALEPtr)
    output = Array{UInt8, 1}(undef, getRAMSize(ale))
    _getRAM(ale, output)
    return output
end

# Screen function utilities
getScreenWidth(ale::ALEPtr)::Int =
    ccall((:getScreenWidth, libale_c), Cint, (ALEPtr,), ale)
getScreenHeight(ale::ALEPtr)::Int =
    ccall((:getScreenHeight, libale_c), Cint, (ALEPtr,), ale)
getScreen!(ale::ALEPtr, screen_data::Vector{Cuchar}) =
    ccall((:getScreen, libale_c), Cvoid, (ALEPtr, Ptr{Cuchar}), ale, screen_data)
getScreenRGB!(ale::ALEPtr, output_buffer::Vector{Cuchar}) =
    ccall((:getScreenRGB, libale_c), Cvoid, (ALEPtr, Ptr{Cuchar}),
        ale, output_buffer)
getScreenGrayscale!(ale::ALEPtr, output_buffer::Vector{Cuchar}) =
    ccall((:getScreenGrayscale, libale_c), Cvoid, (ALEPtr, Ptr{Cuchar}),
        ale, output_buffer)

"""
    getScreen(ale_instance::ALEPtr)

Returns the current game screen.
"""
function getScreen(ale::ALEPtr)
    w = getScreenWidth(ale)
    h = getScreenHeight(ale)
    screen_data = Array{Cuchar}(undef, w*h) # row-major order
    getScreen!(ale, screen_data)
    screen_data
end

"""
    getScreenRGB(ale_instance::ALEPtr)

Returns the RGB representation of the game screen as a (width * heigth * 3)-element Array{UInt, 1}.
"""
function getScreenRGB(ale::ALEPtr)
    w = getScreenWidth(ale)
    h = getScreenHeight(ale)
    screen_data = Array{Cuchar}(undef, w*h*3)
    getScreenRGB!(ale, screen_data)
    return screen_data
end

"""
    getScreenGrayscale(ale_instance::ALEPtr)

Returns the grayscale representation of the screen.
"""
function getScreenGrayscale(ale::ALEPtr)
    w = getScreenWidth(ale)
    h = getScreenHeight(ale)
    screen_data = Array{Cuchar}(undef, w*h)
    getScreenGrayscale!(ale, screen_data)
    return screen_data
end

"""
    saveState(ale_instance::ALEPtr)

Saves the state of the game
"""
saveState(ale::ALEPtr) = ccall((:saveState, libale_c), ALEState, (ALEPtr,), ale)

"""
    loadState(ale_instance::ALEPtr)

Loads the state of the game
"""
loadState(ale::ALEPtr) = ccall((:loadState, libale_c), ALEState, (ALEPtr,), ale)

"""
    cloneState(ale_instance::ALEPtr)

This makes a copy of the environment state. This copy does *not* include pseudorandomness,
making it suitable for planning purposes.

See also: [`cloneSystemState`](@ref)
"""
cloneState(ale::ALEPtr) = ccall((:cloneState, libale_c), ALEStatePtr,
    (ALEPtr,), ale)

"""
    restoreState(ale_instance::ALEPtr, state::ALEStatePtr)

Reverse operation of cloneState(). This doesn't restore pseudorandomness, so that
repeated calls to restoreState() in the stochastic controls settinig will not lead to the same outcomes.

See also: [`restoreSystemState`](@ref)
"""
restoreState(ale::ALEPtr, state_ref::ALEStatePtr) =
    ccall((:restoreState, libale_c), Cvoid, (ALEPtr, ALEStatePtr), ale, state_ref)

"""
    cloneSystemState(ale_instance::ALEPtr)

This makes a copy of the system & environment state, suitable for serialization. This includes
pseudorandomness and so is *not* suitable for planning purposes.
"""
cloneSystemState(ale::ALEPtr) = ccall((:cloneSystemState, libale_c),
    ALEStatePtr, (ALEPtr,), ale)

"""
    restoreSystemState(ale_instance::ALEPtr, state::ALEStatePtr)

Reverse operation of [`cloneSystemState`](@ref).
"""
restoreSystemState(ale::ALEPtr, state_ref::ALEStatePtr) =
    ccall((:restoreSystemState, libale_c), Cvoid, (ALEPtr, ALEStatePtr),
        ale, state_ref)
deleteState(state::ALEStatePtr) = ccall((:deleteState, libale_c), Cvoid,
    (ALEStatePtr,), state)

"""
    saveScreenPNG(ale_instance::ALEPtr, filename::String)

Save the current screen as a png file.

# Example
```julia-repl
julia> saveScreenPNG(ale, "pongscreenshot.png")
```
"""
saveScreenPNG(ale::ALEPtr, filename::String) = ccall((:saveScreenPNG, libale_c),
    Cvoid, (ALEPtr, Ptr{Cchar}), ale, filename)

function encodeState(state::ALEStatePtr)
    len = encodeStateLen(state)
    buf = Array{Cchar}(undef, len)
    encodeState!(state, buf)
end

function encodeState!(state::ALEStatePtr, buf::Array{Cchar})
    ccall((:encodeState, libale_c), Cvoid, (ALEStatePtr, Ptr{Cchar}, Cint),
        state, buf, length(buf))
    buf
end

encodeStateLen(state::ALEStatePtr) = ccall((:encodeStateLen, libale_c),
    Cint, (ALEStatePtr,), state)

decodeState(buf::Array{Cchar,1}) = ccall((:decodeState, libale_c),
    ALEStatePtr, (Ptr{Cchar}, Cint), buf, length(buf))

"""
    setLoggerMode!(mode::Symbol)

Sets the mode for the Logger for the ALE instance. Three modes(::Symbol) are available.
    :info    ==> logs all details and information
    :warning ==> logs warnings and errors
    :error   ==> logs errors only

# Example
```julia-repl
julia> setLoggerMode!(:info)
```
"""
function setLoggerMode!(mode::Symbol)
    @assert mode ∈ keys(Logger) "$mode is unavailable. Please select one of :info, :warning, :error"
    ccall((:setLoggerMode, libale_c), Cvoid, (Cint,), Logger[mode])
end


text:
# MadNLP.jl
# Created by Sungho Shin (sungho.shin@wisc.edu)

module MadNLPMa97

import ..MadNLPHSL:
    @kwdef, Logger, @debug, @warn, @error, libhsl,
    AbstractOptions, AbstractLinearSolver, set_options!, SparseMatrixCSC, SubVector, StrideOneVector,
    SymbolicException,FactorizationException,SolveException,InertiaException,
    introduce, factorize!, solve!, improve!, is_inertia, inertia
import ..MadNLPHSL: Mc68

const INPUT_MATRIX_TYPE = :csc

@enum(Ordering::Int32,AMD = 1, METIS = 3)
@enum(Scaling::Int32,SCALING_NONE = 0, MC64 = 1, MC77 = 2, MC30 = 4)

@kwdef mutable struct Options <: AbstractOptions
    ma97_num_threads::Int = 1
    ma97_print_level::Int = -1
    ma97_nemin::Int = 8
    ma97_small::Float64 = 1e-20
    ma97_order::Ordering = METIS
    ma97_scaling::Scaling = SCALING_NONE
    ma97_u::Float64 = 1e-8
    ma97_umax::Float64 = 1e-4
end

@kwdef mutable struct Control
    f_arrays::Cint = 0
    action::Cint = 0
    nemin::Cint = 0
    multiplier::Cdouble = 0.
    ordering::Cint = 0
    print_level::Cint = 0
    scaling::Cint = 0
    small::Cdouble = 0
    u::Cdouble = 0
    unit_diagnostics::Cint = 0
    unit_error::Cint = 0
    unit_warning::Cint = 0
    factor_min::Clong = 0
    solve_blas3::Cint = 0
    solve_min::Clong = 0
    solve_mf::Cint = 0
    consist_tol::Cdouble = 0
    ispare::Vector{Cint}
    rspare::Vector{Cdouble}
end

@kwdef mutable struct Info
    flag::Cint = 0
    flag68::Cint = 0
    flag77::Cint = 0
    matrix_dup::Cint = 0
    matrix_rank::Cint = 0
    matrix_outrange::Cint = 0
    matrix_missing_diag::Cint = 0
    maxdepth::Cint = 0
    maxfront::Cint = 0
    num_delay::Cint = 0
    num_factor::Clong = 0
    num_flops::Clong = 0
    num_neg::Cint = 0
    num_sup::Cint = 0
    num_two::Cint = 0
    ordering::Cint = 0
    stat::Cint = 0
    ispare::Vector{Cint}
    rspare::Vector{Cdouble}
end

mutable struct Solver <:AbstractLinearSolver
    n::Int32

    csc::SparseMatrixCSC{Float64,Int32}

    control::Control
    info::Info

    akeep::Vector{Ptr{Nothing}}
    fkeep::Vector{Ptr{Nothing}}

    opt::Options
    logger::Logger
end

ma97_default_control_d(control::Control) = ccall(
    (:ma97_default_control_d,libhsl),
    Nothing,
    (Ref{Control},),
    control)

ma97_analyse_d(check::Cint,n::Cint,ptr::StrideOneVector{Cint},row::StrideOneVector{Cint},
               val::Ptr{Nothing},akeep::Vector{Ptr{Nothing}},
               control::Control,info::Info,
               order::Ptr{Nothing}) = ccall(
                         (:ma97_analyse_d,libhsl),
                         Nothing,
                         (Cint,Cint,Ptr{Cint},Ptr{Cint},Ptr{Cdouble},
                          Ptr{Ptr{Nothing}},Ref{Control},Ref{Info},Ptr{Cint}),
                         check,n,ptr,row,val,akeep,control,info,order)
ma97_factor_d(matrix_type::Cint,ptr::Ptr{Nothing},row::Ptr{Nothing},
              val::StrideOneVector{Cdouble},akeep::Vector{Ptr{Nothing}},fkeep::Vector{Ptr{Nothing}},
              control::Control,info::Info,scale::Ptr{Nothing}) = ccall(
                  (:ma97_factor_d,libhsl),
                  Nothing,
                  (Cint,Ptr{Cint},Ptr{Cint},Ptr{Cdouble},Ptr{Ptr{Nothing}},
                   Ptr{Ptr{Nothing}},Ref{Control},Ref{Info},Ptr{Cdouble}),
                  matrix_type,ptr,row,val,akeep,fkeep,control,info,scale)
ma97_solve_d(job::Cint,nrhs::Cint,x::StrideOneVector{Cdouble},ldx::Cint,
             akeep::Vector{Ptr{Nothing}},fkeep::Vector{Ptr{Nothing}},
             control::Control,info::Info) = ccall(
                 (:ma97_solve_d,libhsl),
                 Nothing,
                 (Cint,Cint,Ptr{Cdouble},Cint,Ptr{Ptr{Nothing}},
                  Ptr{Ptr{Nothing}},Ref{Control},Ref{Info}),
                 job,nrhs,x,ldx,akeep,fkeep,control,info)
ma97_finalize_d(akeep::Vector{Ptr{Nothing}},fkeep::Vector{Ptr{Nothing}})=ccall(
    (:ma97_finalise_d,libhsl),
    Nothing,
    (Ptr{Ptr{Nothing}},Ptr{Ptr{Nothing}}),
    akeep,fkeep)
ma97_set_num_threads(n) = ccall((:omp_set_num_threads_,libhsl),
                                Cvoid,
                                (Ref{Cint},),
                                Cint(n))

function Solver(csc::SparseMatrixCSC{Float64,Int32};
                option_dict::Dict{Symbol,Any}=Dict{Symbol,Any}(),
                opt=Options(),logger=Logger(),
                kwargs...)

    set_options!(opt,option_dict,kwargs)

    ma97_set_num_threads(opt.ma97_num_threads)

    n = Int32(csc.n)

    info = Info(ispare=zeros(Int32,5),rspare=zeros(Float64,10))
    control=Control(ispare=zeros(Int32,5),rspare=zeros(Float64,10))
    ma97_default_control_d(control)

    control.print_level = opt.ma97_print_level
    control.f_arrays = 1
    control.nemin = opt.ma97_nemin
    control.ordering = Int32(opt.ma97_order)
    control.small = opt.ma97_small
    control.u = opt.ma97_u
    control.scaling = Int32(opt.ma97_scaling)

    akeep = [C_NULL]
    fkeep = [C_NULL]

    ma97_analyse_d(Int32(1),n,csc.colptr,csc.rowval,C_NULL,akeep,control,info,C_NULL)
    info.flag<0 && throw(SymbolicException())
    M = Solver(n,csc,control,info,akeep,fkeep,opt,logger)
    finalizer(finalize,M)
    return M
end
function factorize!(M::Solver)
    ma97_factor_d(Int32(4),C_NULL,C_NULL,M.csc.nzval,M.akeep,M.fkeep,M.control,M.info,C_NULL)
    M.info.flag<0 && throw(FactorizationException())
    return M
end
function solve!(M::Solver,rhs::StrideOneVector{Float64})
    ma97_solve_d(Int32(0),Int32(1),rhs,M.n,M.akeep,M.fkeep,M.control,M.info)
    M.info.flag<0 && throw(SolveException())
    return rhs
end
is_inertia(::Solver)=true
function inertia(M::Solver)
    return (M.info.matrix_rank-M.info.num_neg,M.n-M.info.matrix_rank,M.info.num_neg)
end
finalize(M::Solver) = ma97_finalize_d(M.akeep,M.fkeep)

function improve!(M::Solver)
    if M.control.u == M.opt.ma97_umax
        @debug(M.logger,"improve quality failed.")
        return false
    end
    M.control.u = min(M.opt.ma97_umax,M.control.u^.75)
    @debug(M.logger,"improved quality: pivtol = $(M.control.u)")
    return true
end
introduce(::Solver)="ma97"

end


text:
# ---
# title: 596. Classes More Than 5 Students
# id: problem596
# author: Tian Jun
# date: 2020-10-31
# difficulty: Easy
# categories:
# link: <https://leetcode.com/problems/classes-more-than-5-students/description/>
# hidden: true
# ---
#
# There is a table `courses` with columns: **student** and **class**
#
# Please list out all classes which have more than or equal to 5 students.
#
# For example, the table:
#
#
#
#     +---------+------------+
#     | student | class      |
#     +---------+------------+
#     | A       | Math       |
#     | B       | English    |
#     | C       | Math       |
#     | D       | Biology    |
#     | E       | Math       |
#     | F       | Computer   |
#     | G       | Math       |
#     | H       | Math       |
#     | I       | Math       |
#     +---------+------------+
#
#
# Should output:
#
#
#
#     +---------+
#     | class   |
#     +---------+
#     | Math    |
#     +---------+
#
#
#
#
# **Note:**
# The students should not be counted duplicate in each course.
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
module InputTextWithHints

Base.Experimental.@optlevel 1

using Redux
using CImGui
using ..InputTexts
import ..InputTexts: AbstractInputTextAction, get_string, get_size
import ..ReduxImGui: get_label

# actions
abstract type AbstractInputTextWithHintAction <: AbstractInputTextAction end

"""
    SetHint(label, hint)
Set hint to `hint`.
"""
struct SetHint <: AbstractInputTextAction
    label::String
    hint::String
end

# state
struct State <: AbstractImmutableState
    input_text::InputTexts.State
    hint::String
end
State(label::AbstractString, hint::AbstractString; varargs...) =
    State(InputTexts.State(label, varargs...), hint)

# reducers
input_text_with_hint(state::AbstractState, action::AbstractAction) = state
input_text_with_hint(state::Vector{<:AbstractState}, action::AbstractAction) = state
input_text_with_hint(state::Dict{String,<:AbstractState}, action::AbstractAction) = state

function input_text_with_hint(state::Dict{String,State}, action::AbstractInputTextAction)
    s = Dict{String,State}()
    for (k,v) in state
        s[k] = get_label(v) == action.label ? InputTexts.input_text(v, action) : v
    end
    return s
end

function input_text_with_hint(state::Dict{String,State}, action::AbstractInputTextWithHintAction)
    s = Dict{String,State}()
    for (k,v) in state
        s[k] = get_label(v) == action.label ? input_text_with_hint(v, action) : v
    end
    return s
end

input_text_with_hint(s::State, a::SetHint) = State(s.input_text, a.hint)

# helper
"""
    InputTextWithHint(store::AbstractStore, get_state=Redux.get_state) -> Bool
Return `true` when triggered.
"""
function InputTextWithHint(store::AbstractStore, get_state=Redux.get_state)
    s = get_state(store)
    is_triggered = CImGui.InputTextWithHint(get_label(s), s.hint, get_string(s), get_size(s))
    return is_triggered
end

get_label(s::State) = get_label(s.input_text)
get_size(s::State) = get_size(s.input_text)

"""
    get_string(s::State) -> String
Return the current string.
"""
get_string(s::State) = get_string(s.input_text)

end # module


text:
import DataStructures

function detect_anagrams(subject, candidates)
  letters = DataStructures.DefaultDict(0)
  results = []

  for letter in subject
    letters[lowercase(letter)] += 1
  end

  for candidate in candidates
    if lowercase(candidate) == lowercase(subject)
      continue
    end

    candidate_letters = DataStructures.DefaultDict(0)

    for letter in candidate
      candidate_letters[lowercase(letter)] += 1
    end

    if letters == candidate_letters
      append!(results, [ candidate ])
    end
  end

  return results
end


text:
toxla(x::XScalar) = Any[]
toxla(x::XArray) = Any[x]
toxla(x::Array{<:XScalar}) = Any[x]
toxla(x) = vcat(map(f -> toxla(getfield(x, f)), fieldnames(typeof(x)))...)

rebuild(T::AType{<:Union{XScalar,Array,Bool}}, xs) = T isa Const ? T.value : popfirst!(xs)

# TODO: this is a hack.
# Where types match we can construct the type directly.
# Where they don't (i.e. XArray) we should go via Functors.jl.
constructor(T::Type) = T.name.wrapper
constructor(T::Type{<:Tuple}) = tuple
constructor(T::Type{<:NamedTuple{K}}) where K = (xs...) -> NamedTuple{K}(xs)

function rebuild(T::Partial, xs)
  xs = map(T -> rebuild(T, xs), T.value)
  T = constructor(widen(T))
  return T(xs...)
end

function rebuild(T::Type{<:Tuple}, xs) # TODO generalise
  xs = rebuild.(T.parameters, (xs,))
  tuple(xs...)
end

rebuild(T::Mjolnir.Node, xs) = rebuild(widen(T), xs)
rebuild(T::Mjolnir.Shape, xs) = rebuild(widen(T), xs)
rebuild(T::Const, xs) = T.value

printstuff(x) = x
printstuff(x::Print) = nothing
function printstuff(x::Tuple{Print,Any})
  printstuff(x[1])
  println(x[1].data...)
  return printstuff(x[2])
end

function trace(Ts...)
  ir = Mjolnir.trace(Primitives(), Ts...)
  return ir |> broadcasts! |> prints!
end

struct XFunction
  func
  cache::Dict
end

xla(f) = XFunction(f, Dict())

function (f::XFunction)(args...)
  T = xtypeof.(args)
  if haskey(f.cache, T)
    (xla_f, out) = f.cache[T]
  else
    ir = trace(Const(f.func), T...)
    out = IRTools.returntype(blocks(ir)[end])
    deletearg!(ir, 1) # `f` is constant
    ir = convert_xla!(ir, ptuple(T...))
    xla_f, = f.cache[T] = XLA.compile(ir), out
  end
  return rebuild(out, xla_f(toxla(args)...)) |> printstuff
end

function Base.show(io::IO, f::XFunction)
  print(io, "xla(")
  show(io, f.func)
  print(io, ")")
end

invokeoriginal(f::XFunction, args...) = f.func(args...)

instead(::Operations, args, F::AType{<:XFunction}, xs...) =
  ([invokeoriginal, args...], (Const(invokeoriginal), F, xs...))

isxla() = false
@abstract Operations isxla() = Const(true)


text:
module PkgExtensions

export update

import Pkg, PrettyTables
using  Markdown, UUIDs

#
# # Constants
#
LAST_ACTIVE_PROJECT = Ref{Union{Nothing,String}}(nothing)

#
# # Functions
#

#
# ## Show deps
#
function show_deps(ctx, args::Vector, api_opts)
    _show_deps(;api_opts...)
end

disp(x::String)  = x
disp(x::Nothing) = "-"
disp(x::Bool)    = x ? "*" : " "
disp(x::UUID)    = string(x)[1:8]
disp(x)          = string(x)

hl_odd = PrettyTables.Highlighter(
    f = (data,i,j) -> (i % 2) == 0,
    crayon = PrettyTables.Crayon(foreground = :blue),
)

function _show_deps(;only_direct=true)
    deps   = Pkg.dependencies()
    if only_direct
        deps = Dict(uuid => info for (uuid, info) in deps if uuid in values(Pkg.Types.EnvCache().project.deps))
    end
    uuids  = keys(deps)
    fields = fieldnames(Pkg.Types.PackageInfo)
    table  = collect(uuids)
    for field in fields
        column = collect([getproperty(deps[uuid], field) for uuid in uuids])
        table  = hcat(table, column)
    end
    header = ["UUID", string.(fields)...]
    table = disp.(table)
    PrettyTables.pretty_table(table, header, PrettyTables.unicode;
                              crop=:horizontal, highlighters = hl_odd)
end

#
# ## git
#
function do_git(ctx, args, opts)
    cd(dirname(Base.active_project())) do
        try
            run(`git $args`)
        catch
        end
    end
end

#
# ## activate
#
function do_activate(ctx, args, opts)
    temp_env = Base.active_project()
    if haskey(opts, :temp) && opts[:temp]
        Pkg.activate(mktempdir())
    elseif isempty(args)
        Pkg.activate()
    elseif args[1] == "-"
        Pkg.activate(LAST_ACTIVE_PROJECT[])
    else
        Pkg.activate(args...; opts...)
    end
    LAST_ACTIVE_PROJECT[] = temp_env
end

#
# # REPL Interface
#
const CommandDeclaration   = Vector{Pair{Symbol,Any}}
const CompoundDeclarations = Vector{Pair{String,Vector{CommandDeclaration}}}
const OptionDeclaration    = Vector{Pair{Symbol,Any}}

#
# ## Declarations
#
compound_declarations = [
"package" => CommandDeclaration[
[   :name => "git",
    :short_name => "g",
    :handler => do_git,
    :arg_count => 0 => Inf,
    :description => "Run a git command in the active directory.",
    :help => md"""
    git <args>...

Run a `git` command in the active directory.
"""
],[
    :name => "activate",
    :short_name => "a",
    :handler => do_activate,
    :arg_count => 0 => 1,
    :option_spec => OptionDeclaration[
        [:name => "temp", :api => :temp => true],
    ],
    :description => "Activate.",
    :completions => Pkg.REPLMode.complete_activate,
    :help => md"""
    activate [-|<directory>]
    activate --temp

Make the directory the active project.
"""
]],
"show" => CommandDeclaration[
[   :name => "deps",
    :handler => show_deps,
    :description => "Show depdency information.",
    :option_spec => OptionDeclaration[
        [:name => "all", :api => :only_direct => false],
    ],
    :help => md"""
    show deps

Show dependency information.
""",
],
]]

#
# ## Utils
#
function update(declarations::CompoundDeclarations)
    compound_specs = Pkg.REPLMode.CompoundSpecs(declarations)
    for (super, specs) in compound_specs
        existing = get(Pkg.REPLMode.SPECS[], super, nothing)
        # merge specs, overshadow Pkg commands if conflict
        Pkg.REPLMode.SPECS[][super] = existing === nothing ? specs : merge(existing, specs)
    end
end

function __init__()
    update(compound_declarations)
end

end # module