text:
module RungeKutta

# Implements a Runge-Kutta solver for systems of ODEs.
# Currently only the fourth order fixed width version is implemented.

# Each of these functions defines the constants for the Runge-Kutta iteration.
k1{T<:AbstractFloat}(f::Array{Function,1}, x::Array{T,1},
    t::AbstractFloat, h::AbstractFloat) =
    map(fe -> fe(t,x), f);

k2{T<:AbstractFloat}(f::Array{Function,1}, x::Array{T,1}, t::AbstractFloat,
    h::AbstractFloat, k1val::Array{T,1}) =
    map(fe -> fe(t + h/2, x + (0.5*h).*k1val), f);

k3{T<:AbstractFloat}(f::Array{Function,1}, x::Array{T,1}, t::AbstractFloat,
    h::AbstractFloat, k2val::Array{T,1}) =
    map(fe -> fe(t + h/2, x + (0.5*h).*k2val), f);

k4{T<:AbstractFloat}(f::Array{Function,1}, x::Array{T,1}, t::AbstractFloat,
    h::AbstractFloat, k3val::Array{T,1}) =
   map(fe -> fe(t + h, x + h.*k3val), f);

# Computes the value of the next point in the Runge-Kutta iteration.
#
# Args:
#   f: The array of functions for computing x_n+1.
#   x: The current value of the points in the space.
#   t: The current value of time.
#   h: The step size.
#
# Returns: The next point in the solution.
function nextPoint{T<:AbstractFloat}(f::Array{Function,1}, x::Array{T,1},
    t::AbstractFloat, h::AbstractFloat)

    #Pre-calculate the constants.
    k1val = k1(f, x, t, h);
    k2val = k2(f, x, t, h, k1val);
    k3val = k3(f, x, t, h, k2val);
    k4val = k4(f, x, t, h, k3val);

   #This is the basic Runge-Kutta fixed width method.
   return t + h, x + (h/6).*(k1val + 2k2val + 2k3val + k4val);

end #Close nextPoint.

# Solves a system of ODEs using a fixed-width fourth-order Runge-Kutta method.
#
# Args:
#   f: An array of functions defining the system of first-order ODEs such that
#       f[i](t, x_n) = x[i]_{n+1}
#   x0: The initial point.
#   t0: The initial time.
#   h: The step size.
#   n: The number of iterations.
#
# Throws:
#   ArgumentError: If the length of f is not equal to the length of x0.
#   ArgumentError: If n is not greater than zero.
#
# Returns: The vector of times, and the vector of points produced by the solver.
function rk4f{T<:AbstractFloat}(f::Array{Function,1}, x0::Array{T,1},
    t0::AbstractFloat, h::AbstractFloat, n::Integer)

    #Validate that the lengths of f and x0 are the same.
    if length(f) != length(x0)
        throw(ArgumentError("There must be one function per element of x0."));
    end

    #Validate that n is greater than zero.
    if n <= 0
        throw(ArgumentError("Number of iterations must be positive."));
    end

    #Validate that the step size is greater than zero.
    if h <= 0.0
        throw(ArgumentError("Step size must be positive."));
    end

    x = zeros(length(x0), n+1);
    t = zeros(1, n+1);

    #The first points are the initial conditions.
    x[:,1] = x0;
    t[1] = t0;

    for ii in 2:(n+1)
        t[ii], x[:,ii] = nextPoint(f, x[:,ii-1], t[ii-1], h);
    end

    return t,x
end #Close rk4f.

export rk4f

end


text:
using TermWin
using Dates

TermWin.initsession()
#s = "a test"
v = newTwEntry( rootTwScreen, Date; width=25, posy=:center, posx=:center, showHelp=true, box=true )
v.title = "Date: "
v.data.tickSize = 1
activateTwObj( rootTwScreen )
ret = v.value
#unregisterTwObj( rootTwScreen, v )
TermWin.endsession()
println( "You entered ", string( ret ) )

ext: jl


text:
using BioMASS
using Test

@time begin
    @testset "BioMASS.jl" begin
        include("parameter_estimation.jl")
        include("bifurcation_analysis.jl")
    end
end


text:
using Test
using PowerSpectra
using Healpix
using IdentityRanges
using LinearAlgebra

##
@testset "spectralones" begin
    x = spectralones(5)
    @test all(parent(x) .== ones(5))

    x = spectralones(31:35)
    @test all(parent(x) .== ones(5))
    @test all(x[35] == 1.0)

    x = spectralones(0:4)
    @test all(parent(x) .== ones(5))
    @test all(x[0] == 1.0)
end

@testset "spectralzeros" begin
    x = spectralzeros(5)
    @test all(parent(x) .== zeros(5))

    x = spectralzeros(31:35)
    @test all(parent(x) .== zeros(5))
    @test all(x[35] == 0.0)

    x = spectralzeros(0:4)
    @test all(parent(x) .== zeros(5))
    @test all(x[0] == 0.0)

    @test all(x .== SpectralVector(zeros(5), -1))
end

##
@testset "UnitRange and AbstractRange slicing" begin
    A = spectralones(5,5)
    @test typeof(A[0:2,2:4]) == typeof(parent(A))
    @test all(A[0:2,2:4] .== ones(3,3))
    @test typeof(A[IdentityRange(0:2),IdentityRange(2:4)]) == typeof(A)
    @test all(parent(A[IdentityRange(0:2),IdentityRange(2:4)]) .== ones(3,3))
end

##
@testset "SpectralArray typing and index style" begin
    A = spectralones(35:39,5:9)
    @test Base.IndexStyle(typeof(A)) == Base.IndexStyle(typeof(ones(5,5)))
    @test PowerSpectra.parenttype(typeof(A)) == typeof(ones(5,5))
    @test PowerSpectra.parenttype(A) == typeof(ones(5,5))
end

##
@testset "SpectralArray constructor" begin
    A = spectralones(35:39,5:9)
    @test all(A.parent.offsets .== (34, 4))
    A = SpectralArray(ones(5,5), 35:39,5:9)
    @test all(A.parent.offsets .== (34, 4))
    A = SpectralArray(ones(5,5), 34, 4)
    @test all(A.parent.offsets .== (34, 4))
    A = SpectralArray(ones(5,5), (34, 4))
    @test all(A.parent.offsets .== (34, 4))

    @test first(axes(SpectralVector(ones(5),1),1)) == 2
end

##
@testset "SpectralArray linear algebra" begin
    A = spectralzeros(5:10,15:20)
    parent(A) .= 2Matrix(LinearAlgebra.I, 6, 6)
    Ainv = inv(A)
    @test all(parent(A) .== 4parent(Ainv))
    @test all(A.parent.offsets .== Ainv.parent.offsets)
end

max_stars_repo_name: xzackli/PlanckCov

max_issues_repo_name: xzackli/PlanckCov


text:
macro sum(i, loname, hiname, term)
    return quote
        lo = $loname
        hi = $hiname
        tmp = 0.0
        for i in lo:hi
            tmp += $term
        end
        return tmp
    end
end

i = 0
@sum(i, 1, 100, 1.0 / i)

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
# ------------------------------------------------------------------
# Licensed under the MIT License. See LICENSE in the project root.
# ------------------------------------------------------------------

# -----------
# GEOMETRIES
# -----------

"""
    boundingbox(geometry)

Axis-aligned bounding box of the `geometry`.
"""
boundingbox(geom::Geometry) = boundingbox(vertices(geom))

boundingbox(p::Primitive) = boundingbox(boundary(p))

boundingbox(b::Box) = b

function boundingbox(s::Sphere{Dim,T}) where {Dim,T}
  c = center(s)
  r = radius(s)
  r⃗ = Vec(ntuple(i->r, Dim))
  Box(c - r⃗, c + r⃗)
end

# --------
# DOMAINS
# --------

"""
    boundingbox(domain)

Axis-aligned bounding box of the `domain`.
"""
boundingbox(::Domain) = @error "fallback not implemented"

boundingbox(pset::PointSet) = boundingbox([pset[i] for i in 1:nelements(pset)])

boundingbox(grid::CartesianGrid) = Box(extrema(grid)...)

# ---------
# FALLBACK
# ---------

function boundingbox(points::AbstractVector{Point{Dim,T}}) where {Dim,T}
  xmin = MVector(ntuple(i->typemax(T), Dim))
  xmax = MVector(ntuple(i->typemin(T), Dim))
  for p in points
    x = coordinates(p)
    @. xmin = min(x, xmin)
    @. xmax = max(x, xmax)
  end
  Box(Point(xmin), Point(xmax))
end


text:
using Plots
function grid_test(name; xlim=(-5,5),ylim=(-5,5), dims=[1,2], ϵ=-1e-1,order=5,geom=true,condition=100)

	hyperparameters = getParameters(X)
	parameters = copy(θ)

	function ∇cost(x::Number,y::Number)
		parameters[dims] .= x,y
		return ∇loss(F,parameters,X,hyperparameters)[2]
	end

	function cost(x::Number,y::Number)
		parameters[dims] .= x,y
		return loss(F,parameters,X,hyperparameters)
	end

	function finite_differences(x,y)
		parameters[dims] .= x,y
		return finite_difference_gradient(θ->loss(F,θ,X),parameters)
	end

	plot(size=(600,600), xlim=xlim, ylim=ylim, xlabel="parameters, θ")
	x,y = range(xlim...,length=50), range(ylim...,length=50)
	contour!( x, y, cost, alpha=0.5 )

	x,y = range(xlim...,length=12), range(ylim...,length=12)
	grid = collect(Iterators.product(x[2:end-1],y[2:end-1]))

	xGrid = vcat(map(x->x[1], grid)...)
	yGrid = vcat(map(x->x[2], grid)...)

	differences = map( finite_differences, xGrid, yGrid )
	quiver!( xGrid, yGrid, quiver=(ϵ*step(x)*getindex.(differences,dims[1]),ϵ*step(y)*getindex.(differences,dims[2])),
		color=:darkblue, lw=3 )

	gradients = map( ∇cost, xGrid, yGrid )
	quiver!( xGrid, yGrid, quiver=(ϵ*step(x)*getindex.(gradients,dims[1]),ϵ*step(y)*getindex.(gradients,dims[2])),
		color=:gold, lw=2 )

	errors = filter(x->~isnan(x), map( (x,y) -> acos(x'y/(norm(x)*norm(y)))/π, differences, gradients ))
	plot!([],[],color=:gold, lw=3, label="ForwardDiff", title="Mean Error $(100*round(sum(errors)/length(errors),digits=2))%")
	plot!([],[],color=:darkblue, lw=3, label="Central Differences")
	savefig(name)
end


text:
function build_regex_rules(ast::ASTNode, rulenameprefix="")
	assign_rulenames(ast, rulenameprefix)
    rules = Vector{RuleSource}()
    build_regex_rule(ast, rules)
    rules
end

function build_regex_rule(node::ASTNode, rules::Vector{RuleSource})

    if node.func in [:regex]
        build_regex_regex(node, rules)
    elseif node.func in [:or]
        build_regex_or(node, rules)
    elseif node.func in [:and]
        build_regex_and(node, rules)
    elseif node.func in [:quantifier]
        build_regex_quantifier(node, rules)
    elseif node.func in [:optional]
        build_regex_optional(node, rules)
    elseif node.func in [:bracket]
        build_regex_bracket(node, rules)
    elseif node.func in [:terminal]
        build_regex_terminal(node, rules)
    elseif node.func in [:range]
        # do nothing - handled by bracket
    else
	   error("Unexpected regex node with function $(node.func)")
    end

    for child in node.children
        build_regex_rule(child, rules)
    end

end

# "start" rule explicitly converts result to datatype
function build_regex_regex(node::ASTNode, rules::Vector{RuleSource})
    rule = build_rule_shortform_start(node)
    calledrulename = build_called_rulename(node.refs[1])
    push!(rule.source, "convert($(node.args[:datatype]), $(calledrulename)())")
    build_rule_shortform_end(rule, node)
    push!(rules, rule)
end

# apply rule choice point to children
function build_regex_or(node::ASTNode, rules::Vector{RuleSource})
    for child in node.children
        rule = build_rule_shortform_start(node)
        calledrulename = build_called_rulename(child)
        push!(rule.source, "$(calledrulename)()")
        build_rule_shortform_end(rule, node)
        push!(rules, rule)
    end
end

# concenate results of calls to the children
function build_regex_and(node::ASTNode, rules::Vector{RuleSource})
    rule = build_rule_shortform_start(node)
    push!(rule.source, "*(")
    for child in node.children
        calledrulename = build_called_rulename(child)
        push!(rule.source, "$(calledrulename)(),")
    end
    push!(rule.source, ")")
    build_rule_shortform_end(rule, node)
    push!(rules, rule)
end

# apply sequence choice point
function build_regex_quantifier(node::ASTNode, rules::Vector{RuleSource})
    rule = build_rule_shortform_start(node)
    calledrulename = build_called_rulename(node.children[1])
    if node.args[:upperbound] >= typemax(Int)
        push!(rule.source, "join(reps($(calledrulename)(), $(node.args[:lowerbound])))")
    else
        push!(rule.source, "join(reps($(calledrulename)(), $(node.args[:lowerbound]), $(node.args[:upperbound])))")
    end
    build_rule_shortform_end(rule, node)
    push!(rules, rule)
end

# apply choose(Bool)
function build_regex_optional(node::ASTNode, rules::Vector{RuleSource})
    rule = build_rule_shortform_start(node)
    calledrulename = build_called_rulename(node.children[1])
    push!(rule.source, "choose(Bool) ? $(calledrulename)() : \"\"")
    build_rule_shortform_end(rule, node)
    push!(rules, rule)
end

# return constant value
function build_regex_terminal(node::ASTNode, rules::Vector{RuleSource})
    rule = build_rule_shortform_start(node)
    push!(rule.source, "\"$(escape_string(node.args[:value]))\"")
    build_rule_shortform_end(rule, node)
    push!(rules, rule)
end

# use choose(Int) to choose uniformly across all possible ranges (which are the child nodes)
function build_regex_bracket(node::ASTNode, rules::Vector{RuleSource})

    rule = build_rule_start(node)

    cardinality = 0
    for child in node.children
        range = child.args[:value]
        cardinality += length(range)
    end

    if cardinality == 0
        push!(rule.source, "  \"\"")
    else
        push!(rule.source, "  idx = choose(Int, 0, $(cardinality-1))")

        condsource = Vector{AbstractString}()

        range = node.children[1].args[:value]
        unshift!(condsource, "  $(range[1]) + idx")
        offset = length(range)
        for i in 2:length(node.children)
            range = node.children[i].args[:value]
            unshift!(condsource, "  (idx >= $(offset)) ? $(range[1] - offset) + idx :")
            offset += length(range)
        end

        push!(rule.source, "  string(convert(Char, ")
        append!(rule.source, condsource)
        push!(rule.source, "  ))")

    end

    build_rule_end(rule, node)
    push!(rules, rule)

end


text:
# Floating point numbers - basic tests
@test isapprox(4.00000000000001, 4.0)
@test isapprox(5.0,4.999999999999993)
@test !isapprox(4.000000002, 4.00300002)

# Other tolerance levels
@test isapprox(4.32, 4.3; rtol=0.1, atol=0.01)
@test isapprox(1.001, 1.002; rtol=0.001, atol=0.0001)
@test !isapprox(4.5, 4.9; rtol=0.001, atol=0.001)

# Complex numbers
@test isapprox(1.0 + 1.0im, 1.0 + 1.00000000000001im)
@test isapprox(0.9999999999999 + 1.0im, 1.0 + 1.000000000000001im)
@test isapprox(0.9999 + 1.0im, 1.0 + 1.1im; rtol = 0.0001, atol=1.1)

# Complex <-> reals
@test isapprox(1.0 + 0im, 1.0000000000001)
@test isapprox(0.9999999999999, 1.0 + 0im)
@test !isapprox(1.0+1im, 1.000000000000001)

# Comparing NaNs
@test !isapprox(4.0,NaN)
@test !isapprox(NaN,4.0)
@test !isapprox(complex(2.3,NaN), complex(NaN,2.3))
@test !isapprox(NaN, NaN)
@test !isapprox(complex(NaN,NaN), complex(NaN,NaN))
@test !isapprox(complex(NaN,2.3), complex(NaN,2.3))
@test !isapprox(complex(2.3,NaN), complex(2.3,NaN))

# Tests for integers and rationals
@test isapprox(4,4)
@test isapprox(4,5; atol=1)
@test !isapprox(4,6; atol=1)
@test !isapprox(4,5)

@test isapprox(1//2+3im, 1//2+3im)
@test isapprox(1//3, 0.33333333333333333)
@test isapprox(1//3, 0.3333; rtol=0.0001, atol=0.0001)
@test isapprox(1+1im, 1im+1)

# Notably missing from this test suite at the moment
# * Tests for other magnitudes of numbers - very small, very big, and combinations of small and big
# * Tests for various odd combinations of types, e.g. isapprox(x::Integer, y::Rational)

max_stars_repo_name: jakebolewski/julia

max_issues_repo_name: jakebolewski/julia

max_forks_repo_name: jakebolewski/julia


text:
function load_past_months(bot_path; num=2)
    base_dir = joinpath(bot_path, "closed_orders/")
    months = []
    orders = Dict()
    for i in 0:1:num-1
        day = string(today() - Month(i))
        file = joinpath(joinpath(base_dir, split(day, "-")[1]), split(day, "-")[2]) * ".json"
        if isfile(file)
            tmp = JSON.parsefile(file)
            for (k,v) in tmp
                if !(k in ["profit", "trade_num"])
                    for (k1,v1) in v["trades"]
                        orders[k1] =  v1["orders"]
                    end
                end
            end
        end
    end
    return orders
end
export load_past_months


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

@doc raw"""CustomResourceDefinitionVersion describes a version for CRD.

    IoK8sApiextensionsApiserverPkgApisApiextensionsV1CustomResourceDefinitionVersion(;
        additionalPrinterColumns=nothing,
        name=nothing,
        schema=nothing,
        served=nothing,
        storage=nothing,
        subresources=nothing,
    )

    - additionalPrinterColumns::Vector{IoK8sApiextensionsApiserverPkgApisApiextensionsV1CustomResourceColumnDefinition} : additionalPrinterColumns specifies additional columns returned in Table output. See https://kubernetes.io/docs/reference/using-api/api-concepts/#receiving-resources-as-tables for details. If no columns are specified, a single column displaying the age of the custom resource is used.
    - name::String : name is the version name, e.g. “v1”, “v2beta1”, etc. The custom resources are served under this version at `/apis/<group>/<version>/...` if `served` is true.
    - schema::IoK8sApiextensionsApiserverPkgApisApiextensionsV1CustomResourceValidation : schema describes the schema used for validation, pruning, and defaulting of this version of the custom resource.
    - served::Bool : served is a flag enabling/disabling this version from being served via REST APIs
    - storage::Bool : storage indicates this version should be used when persisting custom resources to storage. There must be exactly one version with storage=true.
    - subresources::IoK8sApiextensionsApiserverPkgApisApiextensionsV1CustomResourceSubresources : subresources specify what subresources this version of the defined custom resource have.
"""
mutable struct IoK8sApiextensionsApiserverPkgApisApiextensionsV1CustomResourceDefinitionVersion <: SwaggerModel
    additionalPrinterColumns::Any # spec type: Union{ Nothing, Vector{IoK8sApiextensionsApiserverPkgApisApiextensionsV1CustomResourceColumnDefinition} } # spec name: additionalPrinterColumns
    name::Any # spec type: Union{ Nothing, String } # spec name: name
    schema::Any # spec type: Union{ Nothing, IoK8sApiextensionsApiserverPkgApisApiextensionsV1CustomResourceValidation } # spec name: schema
    served::Any # spec type: Union{ Nothing, Bool } # spec name: served
    storage::Any # spec type: Union{ Nothing, Bool } # spec name: storage
    subresources::Any # spec type: Union{ Nothing, IoK8sApiextensionsApiserverPkgApisApiextensionsV1CustomResourceSubresources } # spec name: subresources

    function IoK8sApiextensionsApiserverPkgApisApiextensionsV1CustomResourceDefinitionVersion(;additionalPrinterColumns=nothing, name=nothing, schema=nothing, served=nothing, storage=nothing, subresources=nothing)
        o = new()
        validate_property(IoK8sApiextensionsApiserverPkgApisApiextensionsV1CustomResourceDefinitionVersion, Symbol("additionalPrinterColumns"), additionalPrinterColumns)
        setfield!(o, Symbol("additionalPrinterColumns"), additionalPrinterColumns)
        validate_property(IoK8sApiextensionsApiserverPkgApisApiextensionsV1CustomResourceDefinitionVersion, Symbol("name"), name)
        setfield!(o, Symbol("name"), name)
        validate_property(IoK8sApiextensionsApiserverPkgApisApiextensionsV1CustomResourceDefinitionVersion, Symbol("schema"), schema)
        setfield!(o, Symbol("schema"), schema)
        validate_property(IoK8sApiextensionsApiserverPkgApisApiextensionsV1CustomResourceDefinitionVersion, Symbol("served"), served)
        setfield!(o, Symbol("served"), served)
        validate_property(IoK8sApiextensionsApiserverPkgApisApiextensionsV1CustomResourceDefinitionVersion, Symbol("storage"), storage)
        setfield!(o, Symbol("storage"), storage)
        validate_property(IoK8sApiextensionsApiserverPkgApisApiextensionsV1CustomResourceDefinitionVersion, Symbol("subresources"), subresources)
        setfield!(o, Symbol("subresources"), subresources)
        o
    end
end # type IoK8sApiextensionsApiserverPkgApisApiextensionsV1CustomResourceDefinitionVersion

const _property_map_IoK8sApiextensionsApiserverPkgApisApiextensionsV1CustomResourceDefinitionVersion = Dict{Symbol,Symbol}(Symbol("additionalPrinterColumns")=>Symbol("additionalPrinterColumns"), Symbol("name")=>Symbol("name"), Symbol("schema")=>Symbol("schema"), Symbol("served")=>Symbol("served"), Symbol("storage")=>Symbol("storage"), Symbol("subresources")=>Symbol("subresources"))
const _property_types_IoK8sApiextensionsApiserverPkgApisApiextensionsV1CustomResourceDefinitionVersion = Dict{Symbol,String}(Symbol("additionalPrinterColumns")=>"Vector{IoK8sApiextensionsApiserverPkgApisApiextensionsV1CustomResourceColumnDefinition}", Symbol("name")=>"String", Symbol("schema")=>"IoK8sApiextensionsApiserverPkgApisApiextensionsV1CustomResourceValidation", Symbol("served")=>"Bool", Symbol("storage")=>"Bool", Symbol("subresources")=>"IoK8sApiextensionsApiserverPkgApisApiextensionsV1CustomResourceSubresources")
Base.propertynames(::Type{ IoK8sApiextensionsApiserverPkgApisApiextensionsV1CustomResourceDefinitionVersion }) = collect(keys(_property_map_IoK8sApiextensionsApiserverPkgApisApiextensionsV1CustomResourceDefinitionVersion))
Swagger.property_type(::Type{ IoK8sApiextensionsApiserverPkgApisApiextensionsV1CustomResourceDefinitionVersion }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_IoK8sApiextensionsApiserverPkgApisApiextensionsV1CustomResourceDefinitionVersion[name]))}
Swagger.field_name(::Type{ IoK8sApiextensionsApiserverPkgApisApiextensionsV1CustomResourceDefinitionVersion }, property_name::Symbol) =  _property_map_IoK8sApiextensionsApiserverPkgApisApiextensionsV1CustomResourceDefinitionVersion[property_name]

function check_required(o::IoK8sApiextensionsApiserverPkgApisApiextensionsV1CustomResourceDefinitionVersion)
    (getproperty(o, Symbol("name")) === nothing) && (return false)
    (getproperty(o, Symbol("served")) === nothing) && (return false)
    (getproperty(o, Symbol("storage")) === nothing) && (return false)
    true
end

function validate_property(::Type{ IoK8sApiextensionsApiserverPkgApisApiextensionsV1CustomResourceDefinitionVersion }, name::Symbol, val)
end


text:
# ---
# title: 176. Second Highest Salary
# id: problem176
# author: Tian Jun
# date: 2020-10-31
# difficulty: Easy
# categories:
# link: <https://leetcode.com/problems/second-highest-salary/description/>
# hidden: true
# ---
#
# Write a SQL query to get the second highest salary from the `Employee` table.
#
#
#
#     +----+--------+
#     | Id | Salary |
#     +----+--------+
#     | 1  | 100    |
#     | 2  | 200    |
#     | 3  | 300    |
#     +----+--------+
#
#
# For example, given the above Employee table, the query should return `200` as
# the second highest salary. If there is no second highest salary, then the
# query should return `null`.
#
#
#
#     +---------------------+
#     | SecondHighestSalary |
#     +---------------------+
#     | 200                 |
#     +---------------------+
#
#
#
## @lc code=start
using LeetCode

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


text:
#   This file is part of Similitude.jl
#   It is licensed under the MIT license
#   Similitude Copyright (C) 2022 Michael Reed
#       _           _                         _
#      | |         | |                       | |
#   ___| |__   __ _| | ___ __ __ ___   ____ _| | __ _
#  / __| '_ \ / _` | |/ / '__/ _` \ \ / / _` | |/ _` |
# | (__| | | | (_| |   <| | | (_| |\ V / (_| | | (_| |
#  \___|_| |_|\__,_|_|\_\_|  \__,_| \_/ \__,_|_|\__,_|
#
#   https://github.com/chakravala
#   https://crucialflow.com

import UnitSystems: isconstant

include("$dir/constant.jl")

@pure factorize(x,a,b,c,d,e,f,g,h) = Similitude.factorize(x,param(a),param(b),param(c),param(d),param(e),param(f),param(g),param(h))

@pure Constant(N::Number) = 𝟏*N
@pure Constant(D::Dimension) = D
@pure Constant(::typeof(MathConstants.φ)) = φ
@pure Constant(::typeof(MathConstants.γ)) = Constant(valueat(35,vals))
@pure Constant(::typeof(ℯ)) = Constant(valueat(36,vals))
@pure Constant(::typeof(π)) = τ/𝟐
@pure Constant(::typeof(exp)) = Constant(ℯ)
@pure Constant(N::Float64) = Constant(factorize(N,τ,𝟐,𝟑,𝟓,𝟕,𝟏𝟏,𝟏𝟗,𝟒𝟑))
@pure Constant(N::Int) = Constant(factorize(N,τ,𝟐,𝟑,𝟓,𝟕,𝟏𝟏,𝟏𝟗,𝟒𝟑))
@pure Constant(N::AbelianGroup) = Constant{N}()

#printone(io::IO,::Val{vals}) = print(io, '𝟏')
Base.show(io::IO,x::Constant{N}) where N  = (showgroup(io,N,Natural,'𝟏'); print(io, " = ", constant(N)))

@pure constant(::Constant{N}) where N = constant(N)
constant(d::LogGroup{B},C=UnitSystems.Universe) where B = log(B,constant(value(d),C))
constant(d::ExpGroup{B},C=UnitSystems.Universe) where B = B^constant(value(d),C)

UnitSystems.unit(x::Constant,y=1) = x
@pure promoteint(v::Constant) = isone(v) ? 1 : v

Base.:^(::UnitSystems.Constant{A},::Constant{B}) where {A,B} = Constant{A^B}()
Base.:*(::UnitSystems.Constant{A},b::Constant) where A = Constant(A)*b
Base.:*(a::Constant,::UnitSystems.Constant{B}) where B = a*Constant(B)
Base.:/(::UnitSystems.Constant{A},b::Constant) where A = Constant(A)/b
Base.:/(a::Constant,::UnitSystems.Constant{B}) where B = a/Constant(B)

#Base.:*(a::Constant{N},b::Dimension{D}) where {N,D} = Dimension{N*D}()
#Base.:*(a::Dimension{D},b::Constant{N}) where {D,N} = Dimension{D*N}()
Base.:/(a::Constant,b::Dimension) = a*inv(b)
Base.:/(a::Dimension,b::Constant) = a*inv(b)

coefprod(a::Constant,b) = a*Constant(b)
coefprod(a,b::Constant) = Constant(a)*b
coefprod(a::Constant,b::Constant) = a*b

Base.:+(a::Constant,b::Similitude.Quantity{D,U}) where {D,U} = U(D)==𝟙 ? Similitude.Quantity{D,U}(a+b.v) : throw(error("$(U(D)) ≠ 𝟙 "))
Base.:+(a::Similitude.Quantity{D,U},b::Constant) where {D,U} = U(D)==𝟙 ? Similitude.Quantity{D,U}(a.v+b) : throw(error("$(U(D)) ≠ 𝟙 "))
Base.:-(a::Constant,b::Similitude.Quantity{D,U}) where {D,U} = U(D)==𝟙 ? Similitude.Quantity{D,U}(a-b.v) : throw(error("$(U(D)) ≠ 𝟙 "))
Base.:-(a::Similitude.Quantity{D,U},b::Constant) where {D,U} = U(D)==𝟙 ? Similitude.Quantity{D,U}(a.v-b) : throw(error("$(U(D)) ≠ 𝟙 "))
Base.:*(a::Constant,b::Similitude.Quantity{D,U}) where {D,U} = Similitude.Quantity{D,U}(a*b.v)
Base.:*(a::Similitude.Quantity{D,U},b::Constant) where {D,U} = Similitude.Quantity{D,U}(a.v*b)
Base.:/(a::Constant,b::Similitude.Quantity{D,U}) where {D,U} = Similitude.Quantity{inv(D),U}(a/b.v)
Base.:/(a::Similitude.Quantity{D,U},b::Constant) where {D,U} = Similitude.Quantity{D,U}(a.v/b)

phys(j,k=vals) = Constant(valueat(j,k))

for i ∈ 1:vals-11
    @eval begin
        export $(Symbol(basis[i]))
        const $(Symbol(basis[i])) = Constant(valueat($i,vals))
    end
end

export factorize

const golden = Constant(valueat(34,vals))
const eulergamma = Constant(valueat(35,vals))
const tau = Constant(valueat(37,vals))
const 𝟏 = Constant(valueat(0,vals))
const two = Constant(valueat(38,vals))
const three = Constant(valueat(39,vals))
const five = Constant(valueat(40,vals))
const seven = Constant(valueat(41,vals))
const eleven = Constant(valueat(42,vals))
const nineteen = Constant(valueat(43,vals))
const fourtythree = Constant(valueat(44,vals))
const 𝟐,𝟑,𝟓,𝟕,𝟏𝟏,𝟏𝟗,𝟒𝟑 = two,three,five,seven,eleven,nineteen,fourtythree
const zetta,zepto,yotta,yocto = (𝟐*𝟓)^21, (𝟐*𝟓)^-21, (𝟐*𝟓)^24, (𝟐*𝟓)^-24
const αinv,φ,τ = inv(α),golden,tau
const RK1990,KJ1990 = RK90,KJ90
const RK2014,KJ2014 = RK,KJ


text:
module BraiderTest
include("braider.jl")
end

module InternalsTest
include("internals.jl")
end

module SystemTest
include("system.jl")
end


text:
using Test
using TestSetExtensions
using LinearAlgebra
using Random
using Qaintessent

##==----------------------------------------------------------------------------------------------------------------------
@testset ExtendedTestSet "test helper functions to apply statevectors" begin
    @testset "test flipqubit!" begin
        N = 5
        n = rand(1:N)
        ψ = Statevector(N)
        ref = zeros(Int, 2^N)
        for i in 1:2^N
            bits = digits(i-1, base=2, pad=N)
            bits[n] = 1 - bits[n]
            ref[i] = sum(bits .* 2 .^[0:N-1...]) + 1
        end

        Qaintessent.flipqubit!(ψ, n)
        @test all(ψ.perm .== ref)
    end

    @testset "test 1 qubit orderqubit!" begin
        N = 5
        n = rand(1:N)
        ψ = Statevector(N)

        ref = zeros(Int, 2^N)
        for i in 1:2^N
            bits = digits(i-1, base=2, pad=N)
            temp = bits[n]
            bits[n:end-1] .= bits[n+1:end]
            bits[end] = temp
            ref[sum(bits .* 2 .^[0:N-1...]) + 1] = i
        end

        Qaintessent.orderqubit!(ψ, (n,))
        @test all(ψ.perm .== ref)
    end

    @testset "test n qubit orderqubit!" begin
        N = 5
        n1 = rand(1:N-1)
        n2 = rand(n1+1:N)
        ψ = Statevector(N)

        ref = zeros(Int, 2^N)
        for i in 1:2^N
            bits = digits(i-1, base=2, pad=N)
            temp1 = bits[n1]
            temp2 = bits[n2]
            bits[n1:n2-1]  .= bits[n1+1:n2]
            bits[n2-1:end-2] .= bits[n2+1:end]
            bits[end-1] = temp1
            bits[end] = temp2
            ref[sum(bits .* 2 .^[0:N-1...]) + 1] = i
        end

        Qaintessent.orderqubit!(ψ, (n1, n2))
        @test all(ψ.perm .== ref)
    end

    @testset "test swapqubit!" begin
        N = 5
        n1 = rand(1:N)
        n2 = rand(setdiff(1:N, n1))
        ψ = Statevector(N)

        ref = zeros(Int, 2^N)
        for i in 1:2^N
            bits = digits(i-1, base=2, pad=N)
            temp = bits[n1]
            bits[n1] = bits[n2]
            bits[n2] = temp
            ref[i] = sum(bits .* 2 .^[0:N-1...]) + 1
        end

        Qaintessent.swapqubit!(ψ, n1, n2)
        @test all(ψ.perm .== ref)
    end
end

"""Checks that tailored apply gives same result as general apply"""

@testset ExtendedTestSet "apply gates to state vector" begin
    N = 6
    θ = 0.7π
    ϕ = 0.4π
    n = randn(3); n /= norm(n)
    ψ = rand(ComplexF64, 2^N)

    @testset "apply basic gates to statevector" begin
        # single qubit gates
        for g in [X, Y, Z, HadamardGate(), SGate(), SdagGate(), TGate(), TdagGate(), RxGate(θ), RyGate(θ), RzGate(θ), RotationGate(θ, n), PhaseShiftGate(ϕ)]
            i = rand(1:N)
            cg = circuit_gate(i, g)
            cga = CircuitGate{1,AbstractGate}(cg.iwire, cg.gate) # generate same gate with type AbstractGate{1}
            s1 = Statevector(deepcopy(ψ))
            s2 = Statevector(deepcopy(ψ))
            apply!(s1, cg)
            apply!(s2, cga)
            @test s1.state ≈ s2.state
            @test s1.state ≈ sparse_matrix(cga, N) * ψ
        end
    end

    @testset "apply multiple basic gates to statevector" begin
        # single qubit gates
        g = [X, Y, Z, HadamardGate(), SGate(), TGate(), RxGate(θ), RyGate(θ), RzGate(θ), RotationGate(θ, n), PhaseShiftGate(ϕ)]
        i = randperm(N)[1:3]
        j = randperm(N)[1:3]
        cgs = circuit_gate.(j, g[i])
        cga = CircuitGate[]
        for cg in cgs
            push!(cga, CircuitGate{1,AbstractGate}(cg.iwire, cg.gate))
        end

        s1 = Statevector(deepcopy(ψ))
        s2 = Statevector(deepcopy(ψ))
        apply!(s1, cgs)
        apply!(s2, cga)
        @test s1.state ≈ s2.state
        @test s1.state ≈ sparse_matrix(cga, N) * ψ
    end

    @testset "apply moments to statevector" begin
        # single qubit gates
        for g in [X, Y, Z, HadamardGate(), SGate(), SdagGate(), TGate(), TdagGate(), RxGate(θ), RyGate(θ), RzGate(θ), RotationGate(θ, n), PhaseShiftGate(ϕ)]
            i = rand(1:N)
            cg = Moment([circuit_gate(i, g)])
            cga = Moment([CircuitGate{1,AbstractGate}(cg[1].iwire, cg[1].gate)]) # generate same gate with type AbstractGate{1}
            s1 = Statevector(deepcopy(ψ))
            s2 = Statevector(deepcopy(ψ))

            apply!(s1, cg)
            apply!(s2, cga)
            @test s1.state ≈ s2.state
            @test s1.state ≈ sparse_matrix(cga, N) * ψ
        end
    end

    @testset "apply multiple moments to statevector" begin
        # single qubit gates
        g = [X, Y, Z, HadamardGate(), SGate(), TGate(), RxGate(θ), RyGate(θ), RzGate(θ), RotationGate(θ, n), PhaseShiftGate(ϕ)]
        i = randperm(N)[1:3]
        j = randperm(N)[1:3]
        cgs = circuit_gate.(j, g[i])
        cga = CircuitGate[]
        for cg in cgs
            push!(cga, CircuitGate{1,AbstractGate}(cg.iwire, cg.gate))
        end
        s1 = Statevector(deepcopy(ψ))
        s2 = Statevector(deepcopy(ψ))

        apply!(s1, Moment(cgs))
        apply!(s2, Moment(cga))
        @test s1.state ≈ s2.state
        @test s1.state ≈ sparse_matrix(cga, N) * ψ
    end

    @testset "apply basic controlled gates to statevector" begin
        # control gate
        for g in [X, Y, Z, HadamardGate(), SGate(), TGate(), RxGate(θ), RyGate(θ), RzGate(θ), RotationGate(θ, n), PhaseShiftGate(ϕ)]
            i = rand(1:N)
            j = rand([1:i-1; i+1:N])
            cg = circuit_gate(i, g, j)
            cga = CircuitGate{2,ControlledGate{typeof(g)}}(cg.iwire, cg.gate)
            s1 = Statevector(deepcopy(ψ))
            s2 = Statevector(deepcopy(ψ))

            apply!(s1, cg)
            apply!(s2, cga)
            @test s1.state ≈ s2.state
            @test s1.state ≈ sparse_matrix(cga, N) * ψ
        end
    end

    @testset "apply swap gate to statevector" begin
        i = rand(1:N)
        j = rand([1:i-1; i+1:N])
        cg = circuit_gate((i, j), SwapGate())
        cga = CircuitGate{2,SwapGate}(cg.iwire, cg.gate)

        s1 = Statevector(deepcopy(ψ))
        s2 = Statevector(deepcopy(ψ))

        apply!(s1, cg)
        apply!(s2, cga)
        @test s1.state ≈ s2.state
        @test s1.state ≈ sparse_matrix(cga, N) * ψ
    end

    @testset "apply controlled swap gate to statevector" begin
        i = rand(1:N)
        j = rand([1:i-1; i+1:N])
        k = rand(setdiff(1:N, [i,j]))
        cg = circuit_gate((i, j), SwapGate(), k)
        cga = CircuitGate{3,ControlledGate{SwapGate}}(cg.iwire, cg.gate)

        s1 = Statevector(deepcopy(ψ))
        s2 = Statevector(deepcopy(ψ))

        apply!(s1, cg)
        apply!(s2, cga)
        @test s1.state ≈ s2.state
        @test s1.state ≈ sparse_matrix(cga, N) * ψ
    end

    @testset "apply 1-qubit MatrixGate to statevector" begin
        # MatrixGate: one qubit
        d = 2
        A = rand(ComplexF64, d, d)
        U, R = qr(A)
        U = Array(U);
        g = MatrixGate(U)
        i = rand(1:N)
        cg = circuit_gate(i, g)
        cga = CircuitGate{1,MatrixGate}(cg.iwire, cg.gate) # generate same gate with type AbstractGate{1}

        s1 = Statevector(deepcopy(ψ))
        s2 = Statevector(deepcopy(ψ))

        apply!(s1, cg)
        apply!(s2, cga)
        @test s1.state ≈ s2.state
        @test s1.state ≈ sparse_matrix(cga, N) * ψ

    end

    @testset "apply k-qubit MatrixGate to statevector" begin
        # MatrixGate: k qubits
        k = rand(1:N)
        A = rand(ComplexF64, 2^k, 2^k)
        U, R = qr(A)
        U = Array(U);
        g = MatrixGate(U)
        iwire = [rand(1:N)]
        for j in 1:k-1
            l = rand(1:N-j)
            i = setdiff([1:N...], iwire)[l]
            push!(iwire, i)
        end
        sort!(iwire)
        cg = circuit_gate((iwire...,), g)
        cga = CircuitGate{k,AbstractGate}((iwire...,), g)
        m = sparse_matrix(cga, N)

        s1 = Statevector(deepcopy(ψ))
        s2 = Statevector(deepcopy(ψ))

        apply!(s1, cg)
        apply!(s2, cga)
        @test s1.state ≈ s2.state
        @test s1.state ≈ sparse_matrix(cga, N) * ψ
    end

    @testset "apply controlled k-qubit MatrixGate to statevector" begin
        # MatrixGate: k qubits
        k = rand(1:N-2)
        A = rand(ComplexF64, 2^k, 2^k)
        U, R = qr(A)
        U = Array(U);
        g = MatrixGate(U)
        iwire = [rand(1:N)]
        for j in 1:k-1
            l = rand(1:N-j)
            i = setdiff([1:N...], iwire)[l]
            push!(iwire, i)
        end

        cntrl = (rand(setdiff([1:N...], iwire)),)

        sort!(iwire)
        cg = circuit_gate((iwire...,), g, cntrl)
        cga = CircuitGate{k+1,ControlledGate{typeof(g)}}(cg.iwire, cg.gate)
        m = sparse_matrix(cga, N)

        s1 = Statevector(deepcopy(ψ))
        s2 = Statevector(deepcopy(ψ))

        apply!(s1, cg)
        apply!(s2, cga)
        @test s1.state ≈ s2.state
        @test s1.state ≈ sparse_matrix(cga, N) * ψ
    end

end


text:
"""
    execute_function_safely(Func::Function, x::Array{T,1}; type_check::Bool = false, quiet_errors::Bool = true) where T<:Real
This function creates a function that executes the function for which a fixed point is sought. It is a helper function that is not exported.
### Takes
 * Func - The function input to fixed_point
 * x - The point at which to evaluate the function.
 * type_check - Should the type stability of the function be checked and a error reported in an input vector of Ints turns into a vector of floats (for instance)
### Returns
A FunctionEvaluationResult containing the following fields:
* Input_ - The input
* Output_ - The output of the FunctionExecutor. May be missing if function could not complete without error.
* Error_ - A enum of type FP_FunctionEvaluationError representing what error occured.
### Examples
Func(x) = sqrt(x)
execute_function_safely(Func, [-1.0,0.0,1.0])
execute_function_safely(Func,[Missing(),0.0,1.0])
execute_function_safely(Func,[7.0,0.0,1.0])
execute_function_safely(Func,[NaN,0.0,1.0])
execute_function_safely(Func,[Inf,0.0,1.0])
execute_function_safely(Func,-1.0)
execute_function_safely(Func,Missing())
execute_function_safely(Func,1.0)
execute_function_safely(Func,NaN)
execute_function_safely(Func,Inf)
"""
function execute_function_safely(Func::Function, x::Array{T,1}; type_check::Bool = false, quiet_errors::Bool = true) where T<:Real
    # Check input
    if sum(ismissing.(x)) > 0
        return FunctionEvaluationResult(x, missing, InputMissingsDetected)
    elseif sum(isnan.(x)) > 0
        return FunctionEvaluationResult(x, missing, InputNAsDetected)
    elseif sum(isinf.(x)) > 0
        return FunctionEvaluationResult(x, missing, InputInfsDetected)
    end
    # Run function.
    lenx = length(x)
    tf_result = Array{T,1}(undef,lenx)
    tf_full_result = missing
    if quiet_errors
        try
            tf_full_result =  Func(deepcopy(x))
        catch
            return FunctionEvaluationResult(x, missing, ErrorExecutingFunction)
        end
    else
        tf_full_result =  Func(deepcopy(x))
    end

    side_effect_to_report = missing
    # Now we if tf_full_result is an array then it is the normal case.
    if isa(tf_full_result,Array)
        tf_result = tf_full_result
    elseif isa(tf_full_result,Tuple)
        # In this case we assume that we have a tuple where the first output
        # is the fp array and the second is a NamedTuple
        tf_result = tf_full_result[1]
        side_effect_to_report = tf_full_result[2]
    else
        error("The Fixedpoint function can only return a vector or a tuple  of which the first entry is the vector for which a fixedpoint is sought and the second is a namedtuple (the contents of which are output for the user but are not used in fixed point acceleration).")
    end
    # Check Output and return.
    if sum(ismissing.(tf_result)) > 0
        return FunctionEvaluationResult(x, tf_result, OutputMissingsDetected      , side_effect_to_report)
    elseif sum(isnan.(tf_result)) > 0
        return FunctionEvaluationResult(x, tf_result, OutputNAsDetected           , side_effect_to_report)
    elseif sum(isinf.(tf_result)) > 0
        return FunctionEvaluationResult(x, tf_result, OutputInfsDetected          , side_effect_to_report)
    elseif (length(tf_result) != length(x))
        return FunctionEvaluationResult(x, tf_result, LengthOfOutputNotSameAsInput, side_effect_to_report)
    elseif type_check & (!isa(tf_result[1], T))
        return FunctionEvaluationResult(x, tf_result, FunctionIsNotTypeStable     , side_effect_to_report)
    else
        return FunctionEvaluationResult(x, tf_result, NoError                     , side_effect_to_report)
    end
end

"""
    fixed_point(func::Function, previous_FixedPointResults::FixedPointResults;
                    Algorithm::FixedPointAccelerationAlgorithm = Anderson,  ConvergenceMetric::Function  = supnorm(input::Array, output::Array) = maximum(abs.(output .- input)),
                    ConvergenceMetricThreshold::R = 1e-10, MaxIter::Integer = Integer(1000), MaxM::Integer = Integer(10), ExtrapolationPeriod::Integer = Integer(7), Dampening::Real = AbstractFloat(1.0),
                    PrintReports::Bool = false, ReportingSigFig::Integer = Integer(10), ReplaceInvalids::InvalidReplacement = NoAction, ConditionNumberThreshold::Real = 1e3, quiet_errors::Bool = true) where R<:Real
    fixed_point(func::Function, Inputs::Array{T, 1};
                    Algorithm::FixedPointAccelerationAlgorithm = Anderson,  ConvergenceMetric::Function  = supnorm(input::Array, output::Array) = maximum(abs.(output .- input)),
                    ConvergenceMetricThreshold::Real = 1e-10, MaxIter::Integer = Integer(1000), MaxM::Integer = Integer(10), ExtrapolationPeriod::Integer = Integer(7), Dampening::Real = AbstractFloat(1.0),
                    PrintReports::Bool = false, ReportingSigFig::Integer = Integer(10), ReplaceInvalids::InvalidReplacement = NoAction, ConditionNumberThreshold::Real = 1e3, quiet_errors::Bool = true) where T<:Real
    fixed_point(func::Function, Inputs::Real;
                    Algorithm::FixedPointAccelerationAlgorithm = Anderson,  ConvergenceMetric::Function  = supnorm(input::Array, output::Array) = maximum(abs.(output .- input)),
                    ConvergenceMetricThreshold::Real = 1e-10, MaxIter::Integer = Integer(1000), MaxM::Integer = Integer(10), ExtrapolationPeriod::Integer = Integer(7), Dampening::Real = AbstractFloat(1.0),
                    PrintReports::Bool = false, ReportingSigFig::Integer = Integer(10), ReplaceInvalids::InvalidReplacement = NoAction, ConditionNumberThreshold::Real = 1e3, quiet_errors::Bool = true)
    fixed_point(func::Function, Inputs::Array{T, 2}; Outputs::Array{T,2} = Array{T,2}(undef,size(Inputs)[1],0),
                    Algorithm::FixedPointAccelerationAlgorithm = Anderson,  ConvergenceMetric::Function  = supnorm(input::Array, output::Array) = maximum(abs.(output .- input)),
                    ConvergenceMetricThreshold::Real = 1e-10, MaxIter::Integer = Integer(1000), MaxM::Integer = Integer(10), ExtrapolationPeriod::Integer = Integer(7), Dampening::Real = AbstractFloat(1.0),
                    PrintReports::Bool = false, ReportingSigFig::Integer = Integer(10), ReplaceInvalids::InvalidReplacement = NoAction, ConditionNumberThreshold::Real = 1e3, quiet_errors::Bool = true) where T<:Real where R<:Real

A function for finding the fixed point of another function
### Takes
 *  func - This is the function for which a fixed point is sought. This function must take and return a vector of the same size dimension.
 *  Inputs - This can be either a 1D-vector of values that is an initial guess for a fixed point or it can be an N x A matrix of previous inputs for which corresponding outputs are available. In this case N is the dimensionality of the fixed point vector you are seeking (Hence each column is a matrix that is input to f) and A is the number of previous Inputs/Outputs that are being provided to the fixed point. Where a matrix is input, a corresponding outputs must be provided or the last column of the outputs matrix is taken as a startpoint guess and the rest of the inputs and output matrices are discarded.
 *  Outputs - This is a matrix of the Function values for each column of the input. It must be provided so that column k of the outputs matrix is equal to Function(Column k of inputs matrix).
 *  Algorithm - This is the fixed point Algorithm to be used. It can be Anderson, Simple, Aitken, Newton, MPE, RRE, VEA or SEA. See vignette and references to see explanations of these Algorithms.
 *  ConvergenceMetric This is a function that takes in a vector of inputs and a table of outputs and returns a scalar. This scalar should be low when convergence is close to being achieved. By default this is the maximum residual by absolute value (the sup norm in the space of residuals).
 *   ConvergenceMetricThreshold This is the threshold for terminating the algorithm. The algorithm will terminate when the scalar that ConvergenceMetric returns is less than ConvergenceMetricThreshold. This can be set to a negative number in which case the algorithm will run until MaxIter is hit or an error occurs (Note that an error is likely in trying to use any Algorithm other than "Simple" when a fixed point is already found).
 *  MaxIter - This is the maximum number of iterates that will be undertaken.
 *  MaxM - This is the maximum number of saved iterates that are used in the Anderson algorithm. It has no effect if another Algorithm is chosen. Note that the number of previous iterates that will actually be used is the minimum of MaxIter, the dimensionality of the f's vector and the number of inputs that have been tried to  previously (the width of the Outputs matrix at each given stage of the algorithm). If PrintReports = TRUE, the number of previous iterates actually used is reported as the algorithm is running.
 *  ExtrapolationPeriod - This is the number of simple iterates to perform before extrapolating. This is used for the MPE, RRE, VEA and SEA Algorithms and has no effect if another Algorithm is chosen Where an epsilon algorithm is used this should be a multiple of 2, ie (4,6,8,etc).
 *  Dampening - This is the dampening parameter. By default it is 1 which means no dampening takes place. It can also be less than 1 (indicating dampening) or more than 1 (indicating extrapolation).
 *  PrintReports - This is a boolean describing whether to print ongoing ConvergenceMetric values for each iterate.
 *  ReportingSigFig - This is the number of significant figures that will be used in printing the convergence values to the console (only if PrintReports is TRUE).
 *  ReplaceInvalids -Sometimes an acceleration algorithm proposed a vector with an invalid coordinate (NaN, Inf or missing). This parameter can be set to ReplaceInvalids (to replace invalid coordinates by the simple iterate values), ReplaceVector (to replace entire vector with a simple iterate) or NoAction (where an imminent error will occur).
 *  ConditionNumberThreshold - This is a threshold for what condition number is acceptable for solving the least squares problem for the Anderson Algorithm. If the condition number is larger than this threshold then fewer previous iterates will be used in solving the problem. This has no effect unless the "Anderson" Algorithm is used.
 *  quiet_errors - If true the function will return everything already calculated as soon as an error occurs. The callstack that lead to the error is not returned however. If false an error will be thrown with a callstack.
### Returns
 * A list containing the fixed_point, the Inputs and corresponding Outputs, and convergence values (which are computed under the "ConvergenceMetric").
   The list will also include a "Finish" statement describing why it has finished. This is often going to be due to either MaxIter or ConvergenceMetricThreshold being
   reached. It may also terminate due to an error in generating a new input guess or using the function with that guess. If this occurs the function will terminate early
   and the "Finish" statement will describe the issue. In this event there will also be additional objects returned in the list "NewInputVector" and possibly
   "NewOutputVector" that are useful in debugging the issue.
### Examples
 #' # For the simplest possible example we can seek the fixed point of the cos function with a scalar.
 #' Inputs = 0.3
 #' Func(x) = cos(x)
 #' A = fixed_point(Func, Inputs; Algorithm = Aitken, Dampening = 0.5)
 #' B = fixed_point(Func, Inputs; Algorithm = Anderson, Dampening = 1.0)
 #'
 #' # For this next one the ConvergenceMetricThreshold is negative so the algorithm
 #' # will keep running until MaxIter is met.
 #' C = fixed_point(Func, Inputs; Algorithm = Simple, MaxIter = 4, ConvergenceMetricThreshold = -1.0)
 #' # But we can continue solving for this fixed point but now switching to the Newton Algorithm.
 #' D = fixed_point(Func, C[:Inputs], C[:Outputs]; Algorithm = Newton)
 #'
 #' # We can also find a 4 dimensional fixed point vector of this function.
 #' Inputs = [0.3, 98, 0, pi]
 #' E = fixed_point(Func, Inputs; Algorithm = Anderson)
 #' F = fixed_point(Func, Inputs; Algorithm = Anderson, MaxM = 4, ReportingSigFig = 13)
"""
function fixed_point(func::Function, previous_FixedPointResults::FixedPointResults;
                    Algorithm::FixedPointAccelerationAlgorithm = Anderson,  ConvergenceMetric::Function  = supnorm(input, output) = maximum(abs.(output .- input)),
                    ConvergenceMetricThreshold::Real = 1e-10, MaxIter::Integer = Integer(1000), MaxM::Integer = Integer(10), ExtrapolationPeriod::Integer = Integer(7), Dampening::Real = AbstractFloat(1.0),
                    PrintReports::Bool = false, ReportingSigFig::Integer = Integer(10), ReplaceInvalids::InvalidReplacement = NoAction, ConditionNumberThreshold::Real = 1e3, quiet_errors::Bool = true)
    Inputs = previous_FixedPointResults.Inputs_
    Outputs = previous_FixedPointResults.Outputs_
    return fixed_point(func, Inputs; Outputs = Outputs, Algorithm = Algorithm, ConvergenceMetric = ConvergenceMetric, ConvergenceMetricThreshold = ConvergenceMetricThreshold,
                       MaxIter = MaxIter, MaxM = MaxM, ExtrapolationPeriod = ExtrapolationPeriod, Dampening = Dampening, PrintReports = PrintReports, ReportingSigFig = ReportingSigFig,
                       ReplaceInvalids = ReplaceInvalids, ConditionNumberThreshold = ConditionNumberThreshold, quiet_errors = quiet_errors)
end
function fixed_point(func::Function, Inputs::Array{T, 1};
                    Algorithm::FixedPointAccelerationAlgorithm = Anderson,  ConvergenceMetric::Function  = supnorm(input, output) = maximum(abs.(output .- input)),
                    ConvergenceMetricThreshold::Real = 1e-10, MaxIter::Integer = Integer(1000), MaxM::Integer = Integer(10), ExtrapolationPeriod::Integer = Integer(7), Dampening::Real = AbstractFloat(1.0),
                    PrintReports::Bool = false, ReportingSigFig::Integer = Integer(10), ReplaceInvalids::InvalidReplacement = NoAction, ConditionNumberThreshold::Real = 1e3, quiet_errors::Bool = true) where T<:Real
    Inputs2 = Array{Float64, 2}(undef,size(Inputs)[1],1)
    Inputs2[:,1] = Inputs
    return fixed_point(func, Inputs2; Algorithm = Algorithm, ConvergenceMetric = ConvergenceMetric, ConvergenceMetricThreshold = ConvergenceMetricThreshold,
                       MaxIter = MaxIter, MaxM = MaxM, ExtrapolationPeriod = ExtrapolationPeriod, Dampening = Dampening, PrintReports = PrintReports, ReportingSigFig = ReportingSigFig,
                       ReplaceInvalids = ReplaceInvalids, ConditionNumberThreshold = ConditionNumberThreshold, quiet_errors = quiet_errors)
end
function fixed_point(func::Function, Inputs::Real;
                    Algorithm::FixedPointAccelerationAlgorithm = Anderson,  ConvergenceMetric::Function  = supnorm(input, output) = maximum(abs.(output .- input)),
                    ConvergenceMetricThreshold::Real = 1e-10, MaxIter::Integer = Integer(1000), MaxM::Integer = Integer(10), ExtrapolationPeriod::Integer = Integer(7), Dampening::Real = AbstractFloat(1.0),
                    PrintReports::Bool = false, ReportingSigFig::Integer = Integer(10), ReplaceInvalids::InvalidReplacement = NoAction, ConditionNumberThreshold::Real = 1e3, quiet_errors::Bool = true)
    Inputs2 = Array{typeof(Inputs), 2}(undef,1,1)
    Inputs2[1,1] = Inputs
    return fixed_point(func, Inputs2; Algorithm = Algorithm, ConvergenceMetric = ConvergenceMetric, ConvergenceMetricThreshold = ConvergenceMetricThreshold,
                       MaxIter = MaxIter, MaxM = MaxM, ExtrapolationPeriod = ExtrapolationPeriod, Dampening = Dampening, PrintReports = PrintReports, ReportingSigFig = ReportingSigFig,
                       ReplaceInvalids = ReplaceInvalids, ConditionNumberThreshold = ConditionNumberThreshold, quiet_errors = quiet_errors)
end
function fixed_point(func::Function, Inputs::Array{T, 2}; Outputs::Array{T,2} = Array{T,2}(undef,size(Inputs)[1],0),
                    Algorithm::FixedPointAccelerationAlgorithm = Anderson,  ConvergenceMetric::Function  = supnorm(input, output) = maximum(abs.(output .- input)),
                    ConvergenceMetricThreshold::Real = 1e-10, MaxIter::Integer = Integer(1000), MaxM::Integer = Integer(10), ExtrapolationPeriod::Integer = Integer(7), Dampening::Real = AbstractFloat(1.0),
                    PrintReports::Bool = false, ReportingSigFig::Integer = Integer(10), ReplaceInvalids::InvalidReplacement = NoAction, ConditionNumberThreshold::Real = 1e3, quiet_errors::Bool = true) where T<:Real where R<:Real
    # This code first tests if the input point is a fixed point. Then if it is not a while loop runs to try to find a fixed point.
    if (ConditionNumberThreshold < 1) error("ConditionNumberThreshold must be at least 1.")  end
    SimpleStartIndex = Integer(size(Outputs)[2])
    if isempty(Outputs)
        if size(Inputs)[2] > 1
            @warn("If you do not give outputs to the function then you can only give one vector of inputs (in a 2d array) to the fixed_pointFunction. So for a function that takes an N dimensional array you should input a Array{Float64}(N,1) array.  As you have input an array of size Array{Float64}(N,k) with k > 1 we have discarded everything but the last column to turn it into a Array{Float64}(N,1) array.\n")
            Inputs = Inputs[:,size(Inputs)[2]]
        end
    else
        if size(Inputs) != size(Outputs)
            @warn("If you input a matrix of outputs as well as a matrix of inputs then inputs and outputs must be the same shape. As they differ in this case the last column of the inputs matrix has been
                  taken as the starting point and everything else discarded.")
            Inputs  = Inputs[:,size(Inputs)[2]]
            Outputs = Array{T,2}(undef,size(Inputs)[1],0)
            SimpleStartIndex = Integer(1)
        end
    end
    LengthOfArray = size(Inputs)[1]
    output_type = T
    # Do an initial run if no runs have been done:
    if isempty(Outputs)
        ExecutedFunction = execute_function_safely(func, Inputs[:,1]; quiet_errors = quiet_errors)
        if ExecutedFunction.Error_ != NoError
            return FixedPointResults(Inputs, Outputs, InvalidInputOrOutputOfIteration; FailedEvaluation_ = ExecutedFunction, Other_Output = ExecutedFunction.Other_Output_)
        end
        output_type = promote_type(typeof(Inputs[1]), typeof(ExecutedFunction.Output_[1]))
        converted_outputs = convert.(Ref(output_type), ExecutedFunction.Output_)
        Outputs = hcat(Outputs, converted_outputs)
        Inputs = convert.(Ref(output_type), Inputs)
    else
        # This ensures that MaxIter refers to max iter excluding any previous passed in results
        MaxIter = MaxIter + size(Outputs)[2]
        # This is to take into account previously passed in simple iterations (without jumps).
        SimpleStartIndex = SimpleStartIndex - (size(put_together_without_jumps(Inputs, Outputs))[2])
    end
    # First running through the last column of Inputs to test if we already have a fixed point.
    iter = Integer(size(Outputs)[2])
    ConvergenceVector = Array{output_type,1}(undef,iter)
    for i in 1:iter
        ConvergenceVector[i] = ConvergenceMetric(Inputs[:,i], Outputs[:,i])
    end
    if ConvergenceVector[iter] < ConvergenceMetricThreshold
        if (PrintReports)
            println("The last column of Inputs matrix is already a fixed point under input convergence metric and convergence threshold")
        end
        return FixedPointResults(Inputs, Outputs, ReachedConvergenceThreshold; ConvergenceVector_ = vec(ConvergenceVector), Other_Output = ExecutedFunction.Other_Output_)
    end
    # Printing a report for initial convergence
    Convergence = ConvergenceVector[iter]
    if (PrintReports)
        println("                                          Algorithm: ", lpad(Algorithm, 8)   , ". Iteration: ", lpad(iter, 5),". Convergence: ", lpad(round(Convergence, sigdigits=ReportingSigFig),ReportingSigFig+4), ". Time: ", now())
    end
    iter = iter + 1
    final_other_output = missing
    while (Convergence > ConvergenceMetricThreshold) & (iter <= MaxIter)
        # Generating new input and output.
        NewInputFunctionReturn = fixed_point_new_input(Inputs, Outputs, Algorithm; MaxM = MaxM,
                                     SimpleStartIndex = SimpleStartIndex,ExtrapolationPeriod = ExtrapolationPeriod,
                                     Dampening = Dampening, ConditionNumberThreshold = ConditionNumberThreshold,
                                     PrintReports = PrintReports, ReplaceInvalids = ReplaceInvalids)
        if PrintReports & (Algorithm != Anderson)
            print(lpad("",42))
        end
        ExecutedFunction = execute_function_safely(func, NewInputFunctionReturn; type_check = true, quiet_errors = quiet_errors)
        if ExecutedFunction.Error_ != NoError
            return FixedPointResults(Inputs, Outputs, InvalidInputOrOutputOfIteration; ConvergenceVector_  = vec(ConvergenceVector), FailedEvaluation_ = ExecutedFunction, Other_Output = ExecutedFunction.Other_Output_)
        end
        final_other_output = ExecutedFunction.Other_Output_
        Inputs  = hcat(Inputs, ExecutedFunction.Input_)
        Outputs = hcat(Outputs, convert(Array{output_type,1}, ExecutedFunction.Output_))
        # Checking and recording convergence
        Convergence = ConvergenceMetric(ExecutedFunction.Input_, ExecutedFunction.Output_)
        ConvergenceVector =  vcat(ConvergenceVector, Convergence)
        # Output of report and going to next iteration.
        if (PrintReports) println("Algorithm: ", lpad(Algorithm,8)   , ". Iteration: ", lpad(iter,5), ". Convergence: ", lpad(round(Convergence, sigdigits=ReportingSigFig),ReportingSigFig+4), ". Time: ", now()) end
        iter  = iter + 1
    end
    Finish = ReachedMaxIter
    if (Convergence < ConvergenceMetricThreshold) Finish = ReachedConvergenceThreshold end
    return FixedPointResults(Inputs, Outputs, Finish; ConvergenceVector_  = vec(ConvergenceVector),  Other_Output = final_other_output)
end

"""
    fixed_point_new_input(Inputs::AbstractArray{T,2}, Outputs::AbstractArray{T,2}, Algorithm::FixedPointAccelerationAlgorithm = Anderson;
                               MaxM::Integer = 10, SimpleStartIndex::Integer = 1, ExtrapolationPeriod::Integer = 1, Dampening::S = AbstractFloat(1),
                               ConditionNumberThreshold::R = AbstractFloat(1000), PrintReports::Bool = false, ReplaceInvalids::InvalidReplacement = NoAction) where R<:Real where S<:Real where T<:Real

This function takes the previous inputs and outputs from the fixed_point function and determines what vector to try next in seeking a fixed point.
### Takes
 *  Inputs - This is an N x A matrix of previous inputs for which corresponding outputs are available. In this case N is the dimensionality of the fixed point vector that is being sought (Hence each column is a matrix that is input to the "Function") and A is the number of previous Inputs/Outputs that are being provided to the fixed point.
* Outputs - This is a matrix of Function values for the each column of the "Inputs" matrix.
* Algorithm - This is the fixed point Algorithm to be used. It can be "Anderson", "Simple", "Aitken", "Newton", "MPE", "RRE", "VEA", "SEA".
* MaxM - This is the number of saved iterates that are used in the Anderson algorithm. This has no role if another Algorithm is used.
* SimpleStartIndex - This is the index for what column in the input/output matrices did the algorithm start doing simple iterates without jumps. This is used for all Algorithms except the simple and Anderson Algorithms where it has no effect.
* ExtrapolationPeriod - This is the number of simple iterates to perform before extrapolating. This is used for the MPE, RRE, VEA and SEA Algorithms and has no effect if another Algorithm is chosen.
* Dampening - This is the dampening parameter. By default it is 1 which means no dampening takes place. It can also be less than 1 (indicating dampening) or more than 1 (indicating extrapolation).
* ConditionNumberThreshold - This is what threshold should be chosen to drop previous iterates if the matrix is ill conditioned. Only used in Anderson acceleration.
* PrintReports - This is a boolean describing whether to print ongoing ConvergenceMetric values for each iterate.
### Returns
 * A nicely formatted string version of the input number for printing to the console.
### Examples
FPFunction = function(x){c(0.5*sqrt(abs(x[1] + x[2])), 1.5*x[1] + 0.5*x[2])}
A = fixed_point( Function = FPFunction, Inputs = [0.3,900], MaxIter = 6, Algorithm = Simple)
NewGuessAnderson = fixed_point_new_input(A[:Inputs], A[:Outputs], Algorithm = Anderson)
NewGuessVEA = fixed_point_new_input(A[:Inputs], A[:Outputs], Algorithm = VEA)
NewGuessMPE = fixed_point_new_input(A[:Inputs], A[:Outputs], Algorithm = MPE)
NewGuessAitken = fixed_point_new_input(A[:Inputs], A[:Outputs], Algorithm = Aitken)
"""
function fixed_point_new_input(Inputs::AbstractArray{T,2}, Outputs::AbstractArray{T,2}, Algorithm::FixedPointAccelerationAlgorithm = Anderson;
                               MaxM::Integer = 10, SimpleStartIndex::Integer = 1, ExtrapolationPeriod::Integer = 1, Dampening::S = AbstractFloat(1),
                               ConditionNumberThreshold::R = AbstractFloat(1000), PrintReports::Bool = false, ReplaceInvalids::InvalidReplacement = NoAction) where R<:Real where S<:Real where T<:Real
    CompletedIters = size(Outputs)[2]
    simple_iterate = Outputs[:,CompletedIters]
    proposed_input = repeat([NaN], size(simple_iterate)[1])
    if Algorithm == Simple
         proposed_input = simple_iterate
    elseif Algorithm == Anderson
        if CompletedIters < 2
            if (PrintReports) print("                           Used:",  lpad(0, 3)," lags. ") end
            proposed_input = simple_iterate
        else
            VectorLength   = size(Outputs)[1]
            M = min(MaxM-1,CompletedIters-1,VectorLength)

            recent_Outputs  = Outputs[:, (CompletedIters-M):CompletedIters]
            recent_Inputs   = Inputs[ :, (CompletedIters-M):CompletedIters]
            Resid           = recent_Outputs .-  recent_Inputs
            DeltaOutputs    = recent_Outputs[:,2:(M+1)] .- recent_Outputs[:,1:M]
            DeltaResids     = Resid[:,2:(M+1)]   .- Resid[:,1:M]
            LastResid       = Resid[:,M+1]
            LastOutput      = recent_Outputs[:,M+1]
            Coeffs          = repeat([NaN], size(DeltaOutputs)[2])
            ConditionNumber = NaN
            while sum(isnan.(Coeffs)) > 0
                if isempty(DeltaResids)
                    # This happens if there is convergence by constant increments and thus the
                    #  most recent DeltaResids is all zeros. So we end up dropping all DeltaResids and get an error here.
                    break
                end
                ConditionNumber = cond(DeltaResids)
                if ConditionNumber > ConditionNumberThreshold
                    M = M-1
                    DeltaOutputs= DeltaOutputs[:, 2:(M+1)]
                    DeltaResids = DeltaResids[ :, 2:(M+1)]
                    Coeffs      = repeat([NaN], size(DeltaOutputs)[2])
                    continue
                end
                Fit = fit(LinearModel,  hcat(DeltaResids), LastResid)
                Coeffs = Fit.pp.beta0
                if sum(isnan.(Coeffs)) > 0
                    M = M-1
                    if (M < 1.5)
                        # This happens occasionally in test cases where the iteration is very close to a fixed point.
                        if (PrintReports) print("                          Used:",  lpad(0, 3)," lags. ") end
                        break
                    end
                    DeltaOutputs = DeltaOutputs[:, 2:(M+1)]
                    DeltaResids  = DeltaResids[ :, 2:(M+1)]
                end
            end
            if isempty(Coeffs)
                if (PrintReports) print("Condition number is ", lpad("NaN", 5),". Used:",  lpad(0, 3)," lags. ") end
                proposed_input = repeat([NaN], VectorLength)
            else
                if (PrintReports) print("Condition number is ", lpad(round(ConditionNumber, sigdigits = 2), 5),". Used:",  lpad(M+1, 3)," lags. ") end
                proposed_input = LastOutput .- (Dampening .* vec(DeltaOutputs * Coeffs))
            end
        end
    elseif Algorithm == Aitken
        if ((CompletedIters + SimpleStartIndex) % 3) == 0
            # If we are in 3rd, 6th, 9th, 12th iterate from when we started Acceleration then we want to do a jumped Iterate,
            # First we extract the guess that started this run of 3 iterates (x), the Function applied to it (fx) and the function applied to that (ffx)
            x = Inputs[: ,(CompletedIters -1)]
            fx = Outputs[:,(CompletedIters -1)]
            ffx = Outputs[:,CompletedIters]
            # Now using the appropriate formula to make a new guess. Note that if a vector is input here it is used elementwise.
            proposed_input = x .- ((fx .- x).^2 ./ (ffx .- 2 .* fx .+ x))
        else
            # We just do a simple iterate. We do an attempt with the latest iterate.
            proposed_input = simple_iterate
        end
    elseif Algorithm == Newton
        if ((CompletedIters + SimpleStartIndex) % 2 == 1) & (CompletedIters > 1)
            # If we are in 3rd, 6th, 9th, 12th iterate from when we started Newton Acceleration then we want to do a Newton Iterate,
            # First we extract the guess that started this run of 3 iterates (x), the Function applied to it (fx) and the function applied to that (ffx)
            xk1  = Inputs[ :,(CompletedIters-1)]
            fxk1 = Outputs[:,(CompletedIters-1)]
            gxk1 = fxk1 .- xk1
            xk   = Inputs[ :,CompletedIters]
            fxk  = Outputs[:,CompletedIters]
            gxk  = fxk .- xk
            #ffx = Outputs[,CompletedIters ]
            # Now using the appropriate formula to make a new guess. Note that if a vector is input here it is used elementwise.
            derivative = (gxk.-gxk1)./(xk .- xk1)
            proposed_input   = xk .- (gxk./derivative)
        else
            # We just do a simple iterate.
            proposed_input = simple_iterate
        end
    elseif (Algorithm == MPE) | (Algorithm == RRE)
        SimpleIteratesMatrix = put_together_without_jumps(Inputs, Outputs)
        if (size(SimpleIteratesMatrix)[2] % ExtrapolationPeriod == 0)
            proposed_input = PolynomialExtrapolation(SimpleIteratesMatrix,Algorithm)
        else
            # We just do a simple iterate.
            proposed_input = simple_iterate
        end
    elseif (Algorithm == VEA) | (Algorithm == SEA)
        SimpleIteratesMatrix = put_together_without_jumps(Inputs, Outputs)
        if (size(SimpleIteratesMatrix)[2] % ExtrapolationPeriod) == 0
            proposed_input = EpsilonExtrapolation(SimpleIteratesMatrix, Algorithm)
        else
            # We just do a simple iterate.
            proposed_input = simple_iterate
        end
    end
    # Now the replacement strategies
    dodgy_entries = isnan.(proposed_input) .| ismissing.(proposed_input) .| isinf.(proposed_input)
    if sum(dodgy_entries) != 0
        if ReplaceInvalids == ReplaceElements
            proposed_input[dodgy_entries] = simple_iterate[dodgy_entries]
        elseif ReplaceInvalids == ReplaceVector
            proposed_input = simple_iterate
        end
    end
    return @. (Dampening .* proposed_input) .+ ((1-Dampening) .* simple_iterate)
end

"""
This function performs Minimal Polynomial extrapolation (MPE) or Reduced Rank Extrapolation (RRE) given a matrix of previous iterates of the function.
### Takes
 * Iterates - A matrix of inputs and outputs excluding jumps. Can be pieced together from Inputs and Outputs matrices of the fixed_point function using the put_together_without_jumps function
 * Algorithm The Algorithm for polynomial extrapolation. Should be either "MPE" for minimal polynomial extrapolation or "RRE" for reduced rank extrapolation.
### Returns
 * A vector containing the extrapolated vector.
"""
function PolynomialExtrapolation(Iterates::AbstractArray{R,2}, Algorithm::FixedPointAccelerationAlgorithm) where R<:Real
    if (Algorithm == MPE)
        TotalColumnsOfIterates = size(Iterates)[2]
        OldDifferences         = Iterates[:,2:(TotalColumnsOfIterates-1)] .- Iterates[:,1:(TotalColumnsOfIterates-2)]
        LastDifference         = Iterates[:,TotalColumnsOfIterates]       .- Iterates[:,(TotalColumnsOfIterates-1)]
        InverseOldDifferences  = pinv(OldDifferences)
        cVector                = -InverseOldDifferences * LastDifference
        cVector                = vcat(cVector,1)
        sumVec                 = sum(cVector)
        return (Iterates[:,2:TotalColumnsOfIterates] * cVector) ./ sumVec
    elseif (Algorithm == RRE)
        TotalColumnsOfIterates   = size(Iterates)[2]
        FirstColumn              = Iterates[:,1]
        Differences              = Iterates[:,2:(TotalColumnsOfIterates)]      - Iterates[:,1:(TotalColumnsOfIterates-1)]
        SecondDifferences        = Differences[:,2:(TotalColumnsOfIterates-1)] - Differences[:,1:(TotalColumnsOfIterates-2)]
        FirstDifference          = Differences[:,1]
        Differences              = Differences[:,1:(TotalColumnsOfIterates-2)]
        InverseSecondDifferences = pinv(SecondDifferences)
        return FirstColumn - ((Differences * InverseSecondDifferences) * FirstDifference)
    else
        error("Invalid Algorithm input. PolynomialExtrapolation function can only take Algorithm as MPE or RRE.")
    end
end

"""
This is a helper function for EpsilonExtrapolation
### Takes
 * Iterates - A matrix representing different iterates with one iterate per column. Can be pieced together from Inputs and Outputs matrices of the fixed_point function using the put_together_without_jumps function
 * Algorithm - Algorithm for epsilon extrapolation. Should be either "VEA" for the vector extrapolation algorithm or "SEA" for the scalar epsilon algorithm.
### Returns
 *  A vector with the extrapolated vector.
"""
function EpsilonExtrapolation(Iterates::AbstractArray{R,2}, Algorithm::FixedPointAccelerationAlgorithm) where R<:Real
    # The function cannot do anything to a one column input so will return input unchanged.
    if (size(Iterates)[2] == 1) return Iterates end
    if (size(Iterates)[2] % 2 == 0) Iterates = Iterates[:,2:size(Iterates)[2]] end
    if (!(Algorithm in [VEA, SEA])) error("Invalid Algorithm input. EpsilonExtrapolation function can only take Algorithm as VEA or SEA") end
    Mat = Iterates
    RowsOfMatrix    = size(Mat)[1]
    TotalColumnsOfMatrix = size(Mat)[2]
    PreviousMatrix = zeros(RowsOfMatrix,(TotalColumnsOfMatrix-1))
    for MatrixColumn in reverse(2:TotalColumnsOfMatrix)
        DiffMatrix = Mat[:,2:MatrixColumn] .- Mat[:,1:(MatrixColumn-1)]
        NewMatrix = PreviousMatrix + EpsilonExtrapolationVectorOfInverses(DiffMatrix, Algorithm)
        PreviousMatrix = Mat[:,2:(MatrixColumn-1)]
        Mat = NewMatrix
    end
    # The function can get NAs from the inversion (ie if differenceMatrix contains a zero
    # for SEA then there is division by zero). To avert this we try with 2 less columns.
    if (any(isnan.(Mat)) | any(ismissing.(Mat)))
        Mat = EpsilonExtrapolation(Iterates[:,3:size(Iterates)[2]],Algorithm)
    end
    return Mat[:,1]
end

"""
This is a helper function for EpsilonExtrapolation
### Takes
 * DifferenceMatrix - The matrix of the differences in elements to be inverted.
 * Algorithm - SEA or VEA.
### Returns
 * A vector of the result of inverting each (column) vector in a mmatrix.
"""
function EpsilonExtrapolationVectorOfInverses(DifferenceMatrix::AbstractArray{T,2}, Algorithm::FixedPointAccelerationAlgorithm) where T<:Real
    if (size(DifferenceMatrix)[1] == 1) | (Algorithm == SEA)
        return 1 ./ DifferenceMatrix
    else
        invs = transpose(pinv(DifferenceMatrix[:,1]))
        if size(DifferenceMatrix)[2] < 2 return invs end
        for i in 2:size(DifferenceMatrix)[2]
            invs = hcat(invs, transpose(pinv(DifferenceMatrix[:,i])))
        end
        return invs
    end
end

"""
This function takes the previous inputs and outputs and assembles a matrix with both excluding jumps.
### Takes
 * Inputs - This is an N x A matrix of previous inputs for which corresponding outputs are available. In this case N is the dimensionality of the fixed point vector that is being sought (and each column is a matrix that is input to the "Function") and A is the number of previous Inputs/Outputs that are being provided to the fixed point.
 * Outputs This is a matrix of "Function" values for each column of the "Inputs" matrix.
 * AgreementThreshold A parameter for determining when a column in Inputs and a column in Outputs match. They are deemed to match if the sum of the absolute values of the difference in the columns is less than AgreementThreshold.
### Returns
 * A matrix of inputs and outputs excluding jumps.
"""
function put_together_without_jumps(Inputs::AbstractArray{T,2}, Outputs::AbstractArray{T,2}, AgreementThreshold::Float64 = 1e-10) where T<:Real
  if (any(size(Inputs) != size(Outputs))) error("Inputs and Outputs matrices are not comformable.") end
  size_of_dims = size(Inputs)

  if (size_of_dims[2] == 1) return(hcat(Inputs, Outputs)) end
  Difference = (Inputs[:,2:(size_of_dims[2])] .- Outputs[:,1:(size_of_dims[2]-1)])
  Sum_Of_Differences = sum(Difference, dims = 1)[1,:]
  Agreements = abs.(Sum_Of_Differences) .< AgreementThreshold
  if (all(Agreements))
      return hcat(Inputs[:,1:(size_of_dims[2])], Outputs[:, size_of_dims[2]])
  else
      LocationsOfBreaks = findall(Agreements .== false)
      LastBreak = LocationsOfBreaks[length(LocationsOfBreaks)]
      return hcat(Inputs[:,(LastBreak+1):(size_of_dims[2])] , Outputs[:, size_of_dims[2]])
  end
end


text:
using Simplex
using CUDA
using Dates

"""
To get these mps files, you need to check ../netlib directory.
"""

fp = open("results.md", "w")
println(fp, "# Numerical Performance ($(now()))\n")

# for instance = ["afiro", "adlittle", "sc50a", "sc50b", "sc105", "sc205"]
for instance = ["sc50a", "sc50b", "sc105", "sc205"]
    println(fp, "## Instance: $instance\n")
    netlib = "../netlib/$instance.mps"
    c, xl, xu, bl, bu, A = Simplex.GLPK_read_mps(netlib)
    m, n = size(A); nnz = length(A.nzval);
    println(fp, "- nrows   : $m")
    println(fp, "- ncols   : $n")
    println(fp, "- nnz     : $nnz")
    println(fp, "- sparsity: $(nnz / (m*n))")
    for use_gpu in [false, true]
        if use_gpu
            println(fp, "\n### GPU\n")
        else
            println(fp, "\n### CPU\n")
        end
        println(fp, "```")
        lp = Simplex.StandardLpData(c, xl, xu, A, bl, bu)
        Simplex.run(lp, gpu = use_gpu, performance_io = fp)
        println(fp, "```")

        println("")
        lp = nothing
    end
end

close(fp)


text:
module Water

using ..Ahorn, Maple

const placements = Ahorn.PlacementDict(
    "Water" => Ahorn.EntityPlacement(
        Maple.Water,
        "rectangle"
    ),
)

Ahorn.minimumSize(entity::Maple.Water) = 8, 8
Ahorn.resizable(entity::Maple.Water) = true, true

Ahorn.selection(entity::Maple.Water) = Ahorn.getEntityRectangle(entity)

function Ahorn.render(ctx::Ahorn.Cairo.CairoContext, entity::Maple.Water, room::Maple.Room)
    x = Int(get(entity.data, "x", 0))
    y = Int(get(entity.data, "y", 0))

    width = Int(get(entity.data, "width", 32))
    height = Int(get(entity.data, "height", 32))

    Ahorn.drawRectangle(ctx, 0, 0, width, height, (0.0, 0.0, 1.0, 0.4), (0.0, 0.0, 1.0, 1.0))
end

end


text:
# Input/output routines

"""
    read(file; swap=true, terse=false) -> s::SACtr

Return the SAC trace as read from file `file` as a `SACtr` object.  If `swap` is false,
then auto-byteswapping is not performed and an error is returned if the file is not
of the assumed endianness.  Autoswapping is reported unless `terse` is `true`.
"""
function read(file; swap::Bool=true, terse::Bool=false, check_npts::Bool=true)
    data = open(file, "r") do f
        Base.read(f)
    end
    SACtr(data, file, swap=swap, terse=terse, check_npts=check_npts)
end

"""
    read_cut(file, b, e; swap=true, terse=false) -> s::SACtr
    read_cut(file, b_header, b_time, e_header, e_time; swap=true, terse=false) -> s::SACtr

Return the trace `s`, cut between beginning time `b` and end time `e` in seconds
relative to the zero time.

Optionally specify headers to which `b_time` and `e_time` are relative as
`b_header` and `e_header` respectively.  For example:

    s = read_cut(file, :a, -1, :f + 20)

If `b` is before the start of the trace, then the trace is read from the beginning;
similarly, if `e` is after the end of the trace, the trace is read up until the end.
In either case, a warning is issued unless `terse` is `true`.
"""
function read_cut(file, hb::Symbol, b::Real, he::Symbol, e::Real; swap=true, terse=false)
    b < e || throw(ArgumentError("Start cut must be before stop cut"))
    len = sac_byte_len
    clen = 2*sac_byte_len
    open(file, "r") do f
        header = Base.read(f, sac_header_len)
        # Determine endianness and act accordingly
        nvhdr = reinterpret(SACInt, header[sac_nvhdr_pos*len+1:(sac_nvhdr_pos+1)*len])[1]
        native = if nvhdr == sac_ver_num
            true
        elseif bswap(nvhdr) == sac_ver_num
            false
        else
            error("File does not appear to be SAC data (nvhdr = $nvhdr)")
        end
        native || swap || error("File $file is non-native-endian but `swap` is `false`")
        # Create new SACtr object, setting the trace to have one zero value
        append!(header, zeros(UInt8, 4))
        s = SACtr(header, check_npts=false)
        # Check args
        isundefined(s, hb) && error("Header $hb is undefined")
        isundefined(s, he) && error("Header $he is undefined")
        b < s.b && !terse &&
            @warn("Start cut $b is before trace start $(s.b); reading from trace start")
        e > s.e && !terse &&
            @warn("End cut $e is after trace end $(s.e); reading to trace end")
        # Now do cutting
        ib = round(Int, (b - getfield(s, hb))/s.delta) + 1
        ie = round(Int, (e - getfield(s, he))/s.delta) + 1
        npts = ie - ib + 1
        @assert ib >= 1 && ie <= s.npts
        s.b += (ib - 1)*s.delta
        s.npts = npts
        seek(f, sac_header_len + len*(ib - 1))
        s.t = reinterpret(SACFloat, Base.read(f, len*npts))
        native || (s.t .= bswap.(s.t))
        update_headers!(s)
    end
end
read_cut(file, b::Real, e::Real; kwargs...) = read_cut(file, :b, b, :b, e; kwargs...)

"""
    file_is_native_endian(file)

Return `true` if `file` is a native-endian (defined by a constant in the module to be
little-endian on this machine) SAC file, and `false` if not.

The heuristic is thus: native-endian files have bytes 305:308 which are
a representation of the value `6`.  `6` is the current magic SAC file version number,
hard-coded into the routine.
"""
function file_is_native_endian(file::String)
    nvhdr = try
        d = open(file, "r") do f
            seek(f, sac_nvhdr_pos*sac_byte_len)
            reinterpret(SACInt, Base.read(f, sac_byte_len))[1]
	    end
    catch err
        error("SAC.file_is_native_endian: Cannot open file '$file' for reading " *
            "or file is not the correct type (error $err)")
    end
    if nvhdr == sac_ver_num
        return true
    elseif bswap(nvhdr) == sac_ver_num
        return false
    else
        error("SAC.file_is_native_endian: File '$file' does not appear to be a " *
            "valid SAC file (nvhdr is $nvhdr)")
    end
end

"Write floats or integers either swapped or native-endian"
_write_swap(swap::Bool, F::IOStream, x) = Base.write(F, swap ? Base.bswap.(x) : x)
"Write a String as a SAC string of the correct length, padded with ' 's"
_write_string(F::IOStream, x::String, maxlen::Integer) =
    Base.write(F, x[1:min(length(x),maxlen)]*" "^(max(0,maxlen-length(x))))

@eval function write(s::SACtr, file; byteswap=sac_force_swap)
    open(file, "w") do f
        # Write header
        $([:(_write_swap(byteswap, f, s.$s)) for s in [sac_float_hdr; sac_int_hdr]]...)
        $([:(_write_swap(byteswap, f, SACInt(s.$s))) for s in sac_bool_hdr]...)
        # No byte-swapping needed for characters, but pad them to the correct length
        _write_string(f, s.kstnm, saccharlen)
        # Handle special case of double-length kevnm  header
        kevnm = isundefined(s.kevnm) ? "-12345  -12345  " : s.kevnm
        _write_string(f, kevnm, 2*saccharlen)
        $([:(_write_string(f, s.$s, saccharlen)) for s in sac_char_hdr[3:end]]...)
        # Trace
        _write_swap(byteswap, f, s.t)
    end
end
@doc """
    write(s::SACtr, file; byteswap)
    write(S::AbstractArray{SACtr}, files; byteswap)

Write a SAC trace `s` to `file`, or a set of traces `S` to a set of files `files`.
Set `byteswap` to `false` to force writing in native-endian format; set to `true`
to write bigendian files (MacSAC type).  The default is to write bigendian format.
""" write

function write(s::AbstractArray{SACtr}, file::Array{String}; args...)
    length(s) == length(file) || error("SAC.write: Arrays must be same length")
    for i = 1:length(s)
        write(s[i], file[i]; args...)
    end
    return
end

"""
    read_wild(pat, dir="."; echo=true) -> A, files

Read files matching globbing pattern `pat` from directory `dir`.
If `echo` is false, do not show which files are being read.

Returns an array of SACtr types `A`, and an array of file names `files`.
**NB:** If no files are matched, then empty `SACtr` and `String` arrays
are returned and a warning message printed.
"""
function read_wild(pat::String, dir::String="."; echo::Bool=true)
    files = Glob.glob(pat, dir)
    n = size(files, 1)
    if n == 0
        @info("SAC.read_wild: No files matching '$pat' in directory '$dir'")
        return SACtr[], String[]
    end
    A = Array{SACtr}(undef, n)
    for i = 1:n
        echo && @info("SAC.read: '$(files[i])'")
        A[i] = SAC.read(files[i]; terse=!echo)
    end
    return A, files
end

# List printing
function Base.show(io::IO, s::SACtr)
    out = "SAC.SACtr(delta=$(s[:delta])"
    for hdr in (:b, :npts, :kstnm, :gcarc, :az, :baz)
        if !isundefined(s[hdr])
            out *= ", " * string(hdr) * "=" * strip(string(s[hdr]))
        end
    end
    print(io, out * ")")
end

# Printing to the REPL, for instance
function Base.show(io::IO, ::MIME"text/plain", s::SACtr)
    hdr_string_len = maximum(length.(string.(sac_all_hdr)))
    print(io, "SAC.SACtr:")
    for hdr in sac_all_hdr
        !isundefined(s[hdr]) &&
            print(io, "\n", lpad(string(hdr), hdr_string_len, ' ') * ": ", s[hdr])
    end
end


text:
# implementation/stridedarray.jl
#
# High-level implementation of tensor operations for StridedArray from Julia
# Base Library. Checks dimensions and converts to StridedData before passing
# to low-level (recursive) function.

"""`add!(α, A, conjA, β, C, indCinA)`

Implements `C = β*C+α*permute(op(A))` where `A` is permuted according to `indCinA` and `op` is `conj` if `conjA=Val{:C}` or the identity map if `conjA=Val{:N}`. The indexable collection `indCinA` contains as nth entry the dimension of `A` associated with the nth dimension of `C`.
"""
function add!{CA}(α, A::StridedArray, ::Type{Val{CA}}, β, C::StridedArray, indCinA)
    for i = 1:ndims(C)
        size(A,indCinA[i]) == size(C,i) || throw(DimensionMismatch())
    end

    dims, stridesA, stridesC, minstrides = add_strides(size(C), _permute(_strides(A),indCinA), _strides(C))
    dataA = StridedData(A, stridesA, Val{CA})
    offsetA = 0
    dataC = StridedData(C, stridesC)
    offsetC = 0

    if α == 0
        β == 1 || _scale!(dataC,β,dims)
    elseif α == 1 && β == 0
        add_rec!(_one, dataA, _zero, dataC, dims, offsetA, offsetC, minstrides)
    elseif α == 1 && β == 1
        add_rec!(_one, dataA, _one, dataC, dims, offsetA, offsetC, minstrides)
    elseif β == 0
        add_rec!(α, dataA, _zero, dataC, dims, offsetA, offsetC, minstrides)
    elseif β == 1
        add_rec!(α, dataA, _one, dataC, dims, offsetA, offsetC, minstrides)
    else
        add_rec!(α, dataA, β, dataC, dims, offsetA, offsetC, minstrides)
    end
    return C
end

"""`trace!(α, A, conjA, β, C, indCinA, cindA1, cindA2)`

Implements `C = β*C+α*partialtrace(op(A))` where `A` is permuted and partially traced, according to `indCinA`, `cindA1` and `cindA2`, and `op` is `conj` if `conjA=Val{:C}` or the identity map if `conjA=Val{:N}`. The indexable collection `indCinA` contains as nth entry the dimension of `A` associated with the nth dimension of `C`. The partial trace is performed by contracting dimension `cindA1[i]` of `A` with dimension `cindA2[i]` of `A` for all `i in 1:length(cindA1)`.
"""
function trace!{CA}(α, A::StridedArray, ::Type{Val{CA}}, β, C::StridedArray, indCinA, cindA1, cindA2)
    NC = ndims(C)
    NA = ndims(A)

    for i = 1:NC
        size(A,indCinA[i]) == size(C,i) || throw(DimensionMismatch(""))
    end
    for i = 1:div(NA-NC,2)
        size(A,cindA1[i]) == size(A,cindA2[i]) || throw(DimensionMismatch(""))
    end

    pA = vcat(indCinA, cindA1, cindA2)
    dims, stridesA, stridesC, minstrides = trace_strides(_permute(size(A),pA), _permute(_strides(A),pA), _strides(C))
    dataA = StridedData(A, stridesA, Val{CA})
    offsetA = 0
    dataC = StridedData(C, stridesC)
    offsetC = 0

    if α == 0
        β == 1 || _scale!(dataC, β, dims)
    elseif α == 1 && β == 0
        trace_rec!(_one, dataA, _zero, dataC, dims, offsetA, offsetC, minstrides)
    elseif α == 1 && β == 1
        trace_rec!(_one, dataA, _one, dataC, dims, offsetA, offsetC, minstrides)
    elseif β == 0
        trace_rec!(α, dataA, _zero, dataC, dims, offsetA, offsetC, minstrides)
    elseif β == 1
        trace_rec!(α, dataA, _one, dataC, dims, offsetA, offsetC, minstrides)
    else
        trace_rec!(α, dataA, β, dataC, dims, offsetA, offsetC, minstrides)
    end
    return C
end

"""`contract!(α, A, conjA, B, conjB, β, C, oindA, cindA, oindB, cindB, indCinoAB, [method])`

Implements `C = β*C+α*contract(op(A),op(B))` where `A` and `B` are contracted according to `oindA`, `cindA`, `oindB`, `cindB` and `indCinoAB`. The operation `op` acts as `conj` if `conjA` or `conjB` equal `Val{:C}` or as the identity map if `conjA` (`conjB`) equal `Val{:N}`. The dimension `cindA[i]` of `A` is contracted with dimension `cindB[i]` of `B`. The `n`th dimension of C is associated with an uncontracted (open) dimension of `A` or `B` according to `indCinoAB[n] < NoA ? oindA[indCinoAB[n]] : oindB[indCinoAB[n]-NoA]` with `NoA=length(oindA)` the number of open dimensions of `A`.

The optional argument `method` specifies whether the contraction is performed using BLAS matrix multiplication by specifying `Val{:BLAS}` (default), or using a native algorithm by specifying `Val{:native}`. The native algorithm does not copy the data but is typically slower.
"""
function contract!{CA,CB,TC<:Base.LinAlg.BlasFloat}(α, A::StridedArray, ::Type{Val{CA}}, B::StridedArray, ::Type{Val{CB}}, β, C::StridedArray{TC}, oindA, cindA, oindB, cindB, indCinoAB, ::Type{Val{:BLAS}}=Val{:BLAS})
    NA = ndims(A)
    NB = ndims(B)
    NC = ndims(C)
    TA = eltype(A)
    TB = eltype(B)

    # dimension checking
    dimA = size(A)
    dimB = size(B)
    dimC = size(C)

    cdimsA = dimA[cindA]
    cdimsB = dimB[cindB]
    odimsA = dimA[oindA]
    odimsB = dimB[oindB]
    odimsAB = tuple(odimsA..., odimsB...)

    for i = 1:length(cdimsA)
        cdimsA[i] == cdimsB[i] || throw(DimensionMismatch())
    end
    cdims = cdimsA

    for i = 1:length(indCinoAB)
        dimC[i] == odimsAB[indCinoAB[i]] || throw(DimensionMismatch())
    end

    olengthA = prod(odimsA)
    olengthB = prod(odimsB)
    clength = prod(cdims)

    # permute A
    if CA == :C
        conjA = 'C'
        pA = vcat(cindA, oindA)
        if isa(A, Array{TC}) && pA == collect(1:NA)
            Amat = reshape(A, (clength, olengthA))
        else
            Apermuted = Array{TC}(tuple(cdims..., odimsA...))
            # tensorcopy!(A, 1:NA, Apermuted, pA)
            add!(1, A, Val{:N}, 0, Apermuted, pA)
            Amat = reshape(Apermuted, (clength, olengthA))
        end
    else
        conjA = 'N'
        pA = vcat(oindA, cindA)
        if isa(A, Array{TC}) && pA == collect(1:NA)
            Amat = reshape(A, (olengthA, clength))
        elseif isa(A, Array{TC}) && vcat(cindA, oindA) == collect(1:NA)
            conjA = 'T'
            Amat = reshape(A, (clength, olengthA))
        else
            Apermuted = Array{TC}(tuple(odimsA..., cdims...))
            # tensorcopy!(A, 1:NA, Apermuted, pA)
            add!(1, A, Val{:N}, 0, Apermuted, pA)
            Amat = reshape(Apermuted, (olengthA, clength))
        end
    end

    # permute B
    if CB == :C
        conjB = 'C'
        pB = vcat(oindB, cindB)
        if isa(B, Array{TC}) && pB == collect(1:NB)
            Bmat = reshape(B, (olengthB, clength))
        else
            Bpermuted = Array{TC}(tuple(odimsB..., cdims...))
            # tensorcopy!(B, 1:NB, Bpermuted, pB)
            add!(1, B, Val{:N}, 0, Bpermuted, pB)
            Bmat = reshape(Bpermuted, (olengthB, clength))
        end
    else
        conjB = 'N'
        pB = vcat(cindB, oindB)
        if  isa(B, Array{TC}) && pB == collect(1:NB)
            Bmat = reshape(B, (clength, olengthB))
        elseif isa(B, Array{TC}) && vcat(oindB, cindB) == collect(1:NB)
            conjB = 'T'
            Bmat = reshape(B, (olengthB, clength))
        else
            Bpermuted = Array{TC}(tuple(cdims..., odimsB...))
            # tensorcopy!(B, 1:NB, Bpermuted, pB)
            add!(1, B, Val{:N}, 0, Bpermuted, pB)
            Bmat = reshape(Bpermuted, (clength, olengthB))
        end
    end

    # calculate C
    if isa(C, Array) && indCinoAB == collect(1:NC)
        Cmat = reshape(C, (olengthA, olengthB))
        BLAS.gemm!(conjA, conjB, TC(α), Amat, Bmat, TC(β), Cmat)
    else
        Cmat = Array{TC}(olengthA, olengthB)
        BLAS.gemm!(conjA, conjB, TC(1), Amat, Bmat, TC(0), Cmat)
        # tensoradd!(α, reshape(Cmat, tuple(odimsA..., odimsB...)), pC, β, C, 1:NC)
        add!(α, reshape(Cmat, tuple(odimsA..., odimsB...)), Val{:N}, β, C, indCinoAB)
    end
    return C
end

function contract!{CA,CB}(α, A::StridedArray, ::Type{Val{CA}}, B::StridedArray, ::Type{Val{CB}}, β, C::StridedArray, oindA, cindA, oindB, cindB, indCinoAB, ::Type{Val{:native}}=Val{:native}())
    NA = ndims(A)
    NB = ndims(B)
    NC = ndims(C)

    # dimension checking
    dimA = size(A)
    dimB = size(B)
    dimC = size(C)

    cdimsA = dimA[cindA]
    cdimsB = dimB[cindB]
    odimsA = dimA[oindA]
    odimsB = dimB[oindB]
    odimsAB = tuple(odimsA..., odimsB...)

    # Perform contraction
    pA = vcat(oindA, cindA)
    pB = vcat(oindB, cindB)
    sA = _permute(_strides(A), pA)
    sB = _permute(_strides(B), pB)
    sC = _permute(_strides(C), invperm(indCinoAB))

    dimsA = _permute(size(A), pA)
    dimsB = _permute(size(B), pB)

    dims, stridesA, stridesB, stridesC, minstrides = contract_strides(dimsA, dimsB, sA, sB, sC)
    offsetA = offsetB = offsetC = 0
    dataA = StridedData(A, stridesA, Val{CA})
    dataB = StridedData(B, stridesB, Val{CB})
    dataC = StridedData(C, stridesC)

    # contract via recursive divide and conquer
    if α == 0
        β == 1 || _scale!(dataC, β, dims)
    elseif α == 1 && β == 0
        contract_rec!(_one, dataA, dataB, _zero, dataC, dims, offsetA, offsetB, offsetC, minstrides)
    elseif α == 1 && β == 1
        contract_rec!(_one, dataA, dataB, _one, dataC, dims, offsetA, offsetB, offsetC, minstrides)
    elseif β == 0
        contract_rec!(α, dataA, dataB, _zero, dataC, dims, offsetA, offsetB, offsetC, minstrides)
    elseif β == 1
        contract_rec!(α, dataA, dataB, _one, dataC, dims, offsetA, offsetB, offsetC, minstrides)
    else
        contract_rec!(α, dataA, dataB, β, dataC, dims, offsetA, offsetB, offsetC, minstrides)
    end
    return C
end


text:
param_data(p) =
  p isa Symbol ?
    (p, :Any, missing) :
    p isa Expr ?
      (p.head == :kw ?
         (param_data(p.args[1])[1:2]..., p.args[2]) :
         p.head == :(::) ?
           (p.args..., missing) :
           error("Unknown syntax $p")) :
      error("Unknown syntax $p")

def_data(expr) =
  let (name_params, body) =
        expr isa Expr ?
          (expr.head == :(=) ?
             expr.args :
             (expr.head == :call ?
                (expr, :(throw(UndefinedBackendException()))) :
                error("Unknown syntax $p"))) :
          error("Unknown syntax $p")
    name_params.args[1], name_params.args[2:end], body
  end

# Define for (just the) current backend
macro defcb(expr)
  name, params, body = def_data(expr)
  params_data = map(param_data, params)
  backend_name = Symbol("b_", name)
  esc(
    quote
      export $(name), $(backend_name)
      $(name)($(params...), backend::Backend=top_backend()) =
          $(backend_name)(backend, $(map(pd->pd[1], params_data)...))
      #$(backend_name)(backend::Backend, $(map(name_typ_init->Expr(:(::), name_typ_init[1], name_typ_init[2]), params_data)...)) =
      #    $(body)
    end)
end

# Define for (all the) current backends
macro defcbs(expr)
  name, params, body = def_data(expr)
  params_data = map(param_data, params)
  backend_name = Symbol("b_", name)
  esc(
    quote
      export $(name), $(backend_name)
      $(name)($(params...), backends::Backends=current_backends()) =
        for backend in backends
          $(backend_name)(backend, $(map(pd->pd[1], params_data)...))
        end
      #$(backend_name)(backend::Backend, $(map(name_typ_init->Expr(:(::), name_typ_init[1], name_typ_init[2]), params_data)...)) =
        #  $(body)
    end)
end

#=
@macroexpand @defcbs foo() =  x + 1
@macroexpand @defcbs foo(bar)
@macroexpand @defcbs foo(bar::Baz)
@macroexpand @defcbs foo(bar=quux)
@macroexpand @defcbs foo(bar::Baz=quux)
=#

# backend define

# backend call
macro bcall(backend, name_args)
  name = name_args.args[1]
  args = name_args.args[2:end]
  backend_name = esc(Symbol("b_$name"))
  quote
    $(backend_name)($(esc(backend)), $(esc.(args)...))
  end
end

# backends call
macro bscall(backends, name_args)
  name = name_args.args[1]
  args = name_args.args[2:end]
  backend_name = esc(Symbol("b_$name"))
  quote
    for backend in $(esc(backends))
      $(backend_name)(backend, $(esc.(args)...))
    end
  end
end

# current backend call
macro cbcall(name_args)
  name = name_args.args[1]
  args = name_args.args[2:end]
  backend_name = esc(Symbol("b_$name"))
  quote
	$(backend_name)(top_backend(), $(esc.(args)...))
  end
end

# current backends call
macro cbscall(name_args)
  name = name_args.args[1]
  args = name_args.args[2:end]
  backend_name = esc(Symbol("b_$name"))
  quote
    for backend in current_backends()
	    $(backend_name)(backend, $(esc.(args)...))
    end
  end
end

@defcbs delete_all_refs()
@defcbs unhighlight_all_refs()
#=
export current_layer
current_layer(backends::Backends=current_backends()) =
  [b_current_layer(b) for b in backends]
current_layer(layer, backends::Backends=current_backends()) =
  for (b, l) in zip(backends, layer)
	b_current_layer(b, l)
  end
=#
@defcbs set_layer_active(layer, status)
@defcbs switch_to_layer(layer)

@defcb disable_update()
@defcb enable_update()

@defcbs set_view(camera::Loc, target::Loc, lens::Real=50, aperture::Real=32)
@defcb get_view()
@defcbs b_zoom_extents()
@defcbs set_view_top()

@defcbs set_sun(altitude::Real, azimuth::Real)
@defcbs set_ground(level, material)

@defcbs zoom_extents()
#@defcb get_material(ref::Any)
#@defcbs create_material(name::String)
#@defcb current_material()
#@defcbs current_material(material)
@defcbs set_normal_sky()
@defcbs set_overcast_sky()

export b_render_pathname
b_render_pathname(b::Backend, name::String) =
  render_pathname(name)

export render_view
render_view(name::String="View", backend::Backend=top_backend()) =
  let path = b_render_pathname(backend, name)
    b_render_view(backend, prepare_for_saving_file(path))
    path
  end


text:
using OptimizationProblems

probs = filter(name -> name != :OptimizationProblems
               && name != :penalty2
               && name != :penalty3
               && name != :clplatea
               && name != :clplateb
               && name != :clplatec
               , names(OptimizationProblems))

mpb_probs_1D = (MathProgNLPModel(eval(p)(n_max),  name=string(p) )   for p in probs)