text:
using JuMP, EAGO

                     m = Model()

                     EAGO.register_eago_operators!(m)

                     @variable(m, -1 <= x[i=1:2] <= 1)

                     @variable(m, -20.345339284625858 <= q <= 16.479136897808104)

                     add_NL_constraint(m, :(gelu(0.24719801624822457 + 0.29438950671798914*gelu(0.7291821855933653 + 0.701331929085816*gelu(-0.7470183626892886 + -0.6899614176169906*$(x[1]) + -0.2208825424587011*$(x[2])) + 0.07780867500186073*gelu(-0.11062111710238876 + 0.5516588711870343*$(x[1]) + 0.9886940321128486*$(x[2])) + -0.19203457942980284*gelu(-0.28538228970678503 + -0.152956579710676*$(x[1]) + -0.25273256901725194*$(x[2])) + -0.7065670431246684*gelu(0.7652671472719401 + -0.8637339388616145*$(x[1]) + -0.2539724731598083*$(x[2]))) + -0.7293079760660532*gelu(0.05096633991424904 + 0.31501139336769945*gelu(-0.7470183626892886 + -0.6899614176169906*$(x[1]) + -0.2208825424587011*$(x[2])) + -0.3152903616610754*gelu(-0.11062111710238876 + 0.5516588711870343*$(x[1]) + 0.9886940321128486*$(x[2])) + -0.8990116916055966*gelu(-0.28538228970678503 + -0.152956579710676*$(x[1]) + -0.25273256901725194*$(x[2])) + 0.8373103942919853*gelu(0.7652671472719401 + -0.8637339388616145*$(x[1]) + -0.2539724731598083*$(x[2]))) + 0.8601808654748551*gelu(0.18255055897839467 + -0.7759686413687503*gelu(-0.7470183626892886 + -0.6899614176169906*$(x[1]) + -0.2208825424587011*$(x[2])) + 0.7813463540455929*gelu(-0.11062111710238876 + 0.5516588711870343*$(x[1]) + 0.9886940321128486*$(x[2])) + 0.215367127826533*gelu(-0.28538228970678503 + -0.152956579710676*$(x[1]) + -0.25273256901725194*$(x[2])) + 0.9673448256726829*gelu(0.7652671472719401 + -0.8637339388616145*$(x[1]) + -0.2539724731598083*$(x[2]))) + 0.9720568561360068*gelu(0.5199334546165302 + 0.4081309459110076*gelu(-0.7470183626892886 + -0.6899614176169906*$(x[1]) + -0.2208825424587011*$(x[2])) + 0.6995024782698374*gelu(-0.11062111710238876 + 0.5516588711870343*$(x[1]) + 0.9886940321128486*$(x[2])) + -0.3131098472642604*gelu(-0.28538228970678503 + -0.152956579710676*$(x[1]) + -0.25273256901725194*$(x[2])) + -0.4135863095985335*gelu(0.7652671472719401 + -0.8637339388616145*$(x[1]) + -0.2539724731598083*$(x[2])))) + gelu(-0.49402133773229506 + 0.2465550492028581*gelu(0.7291821855933653 + 0.701331929085816*gelu(-0.7470183626892886 + -0.6899614176169906*$(x[1]) + -0.2208825424587011*$(x[2])) + 0.07780867500186073*gelu(-0.11062111710238876 + 0.5516588711870343*$(x[1]) + 0.9886940321128486*$(x[2])) + -0.19203457942980284*gelu(-0.28538228970678503 + -0.152956579710676*$(x[1]) + -0.25273256901725194*$(x[2])) + -0.7065670431246684*gelu(0.7652671472719401 + -0.8637339388616145*$(x[1]) + -0.2539724731598083*$(x[2]))) + 0.17626257051164496*gelu(0.05096633991424904 + 0.31501139336769945*gelu(-0.7470183626892886 + -0.6899614176169906*$(x[1]) + -0.2208825424587011*$(x[2])) + -0.3152903616610754*gelu(-0.11062111710238876 + 0.5516588711870343*$(x[1]) + 0.9886940321128486*$(x[2])) + -0.8990116916055966*gelu(-0.28538228970678503 + -0.152956579710676*$(x[1]) + -0.25273256901725194*$(x[2])) + 0.8373103942919853*gelu(0.7652671472719401 + -0.8637339388616145*$(x[1]) + -0.2539724731598083*$(x[2]))) + 0.8688281780619751*gelu(0.18255055897839467 + -0.7759686413687503*gelu(-0.7470183626892886 + -0.6899614176169906*$(x[1]) + -0.2208825424587011*$(x[2])) + 0.7813463540455929*gelu(-0.11062111710238876 + 0.5516588711870343*$(x[1]) + 0.9886940321128486*$(x[2])) + 0.215367127826533*gelu(-0.28538228970678503 + -0.152956579710676*$(x[1]) + -0.25273256901725194*$(x[2])) + 0.9673448256726829*gelu(0.7652671472719401 + -0.8637339388616145*$(x[1]) + -0.2539724731598083*$(x[2]))) + 0.3904178182819744*gelu(0.5199334546165302 + 0.4081309459110076*gelu(-0.7470183626892886 + -0.6899614176169906*$(x[1]) + -0.2208825424587011*$(x[2])) + 0.6995024782698374*gelu(-0.11062111710238876 + 0.5516588711870343*$(x[1]) + 0.9886940321128486*$(x[2])) + -0.3131098472642604*gelu(-0.28538228970678503 + -0.152956579710676*$(x[1]) + -0.25273256901725194*$(x[2])) + -0.4135863095985335*gelu(0.7652671472719401 + -0.8637339388616145*$(x[1]) + -0.2539724731598083*$(x[2])))) + gelu(-0.6759660684871824 + 0.9311431633625067*gelu(0.7291821855933653 + 0.701331929085816*gelu(-0.7470183626892886 + -0.6899614176169906*$(x[1]) + -0.2208825424587011*$(x[2])) + 0.07780867500186073*gelu(-0.11062111710238876 + 0.5516588711870343*$(x[1]) + 0.9886940321128486*$(x[2])) + -0.19203457942980284*gelu(-0.28538228970678503 + -0.152956579710676*$(x[1]) + -0.25273256901725194*$(x[2])) + -0.7065670431246684*gelu(0.7652671472719401 + -0.8637339388616145*$(x[1]) + -0.2539724731598083*$(x[2]))) + -0.2015839225926066*gelu(0.05096633991424904 + 0.31501139336769945*gelu(-0.7470183626892886 + -0.6899614176169906*$(x[1]) + -0.2208825424587011*$(x[2])) + -0.3152903616610754*gelu(-0.11062111710238876 + 0.5516588711870343*$(x[1]) + 0.9886940321128486*$(x[2])) + -0.8990116916055966*gelu(-0.28538228970678503 + -0.152956579710676*$(x[1]) + -0.25273256901725194*$(x[2])) + 0.8373103942919853*gelu(0.7652671472719401 + -0.8637339388616145*$(x[1]) + -0.2539724731598083*$(x[2]))) + 0.01676614484829253*gelu(0.18255055897839467 + -0.7759686413687503*gelu(-0.7470183626892886 + -0.6899614176169906*$(x[1]) + -0.2208825424587011*$(x[2])) + 0.7813463540455929*gelu(-0.11062111710238876 + 0.5516588711870343*$(x[1]) + 0.9886940321128486*$(x[2])) + 0.215367127826533*gelu(-0.28538228970678503 + -0.152956579710676*$(x[1]) + -0.25273256901725194*$(x[2])) + 0.9673448256726829*gelu(0.7652671472719401 + -0.8637339388616145*$(x[1]) + -0.2539724731598083*$(x[2]))) + -0.07861065087286123*gelu(0.5199334546165302 + 0.4081309459110076*gelu(-0.7470183626892886 + -0.6899614176169906*$(x[1]) + -0.2208825424587011*$(x[2])) + 0.6995024782698374*gelu(-0.11062111710238876 + 0.5516588711870343*$(x[1]) + 0.9886940321128486*$(x[2])) + -0.3131098472642604*gelu(-0.28538228970678503 + -0.152956579710676*$(x[1]) + -0.25273256901725194*$(x[2])) + -0.4135863095985335*gelu(0.7652671472719401 + -0.8637339388616145*$(x[1]) + -0.2539724731598083*$(x[2])))) + gelu(-0.08101830092609674 + 0.5708385401734484*gelu(0.7291821855933653 + 0.701331929085816*gelu(-0.7470183626892886 + -0.6899614176169906*$(x[1]) + -0.2208825424587011*$(x[2])) + 0.07780867500186073*gelu(-0.11062111710238876 + 0.5516588711870343*$(x[1]) + 0.9886940321128486*$(x[2])) + -0.19203457942980284*gelu(-0.28538228970678503 + -0.152956579710676*$(x[1]) + -0.25273256901725194*$(x[2])) + -0.7065670431246684*gelu(0.7652671472719401 + -0.8637339388616145*$(x[1]) + -0.2539724731598083*$(x[2]))) + -0.9629021826429458*gelu(0.05096633991424904 + 0.31501139336769945*gelu(-0.7470183626892886 + -0.6899614176169906*$(x[1]) + -0.2208825424587011*$(x[2])) + -0.3152903616610754*gelu(-0.11062111710238876 + 0.5516588711870343*$(x[1]) + 0.9886940321128486*$(x[2])) + -0.8990116916055966*gelu(-0.28538228970678503 + -0.152956579710676*$(x[1]) + -0.25273256901725194*$(x[2])) + 0.8373103942919853*gelu(0.7652671472719401 + -0.8637339388616145*$(x[1]) + -0.2539724731598083*$(x[2]))) + 0.4780289946412233*gelu(0.18255055897839467 + -0.7759686413687503*gelu(-0.7470183626892886 + -0.6899614176169906*$(x[1]) + -0.2208825424587011*$(x[2])) + 0.7813463540455929*gelu(-0.11062111710238876 + 0.5516588711870343*$(x[1]) + 0.9886940321128486*$(x[2])) + 0.215367127826533*gelu(-0.28538228970678503 + -0.152956579710676*$(x[1]) + -0.25273256901725194*$(x[2])) + 0.9673448256726829*gelu(0.7652671472719401 + -0.8637339388616145*$(x[1]) + -0.2539724731598083*$(x[2]))) + -0.25224290978348307*gelu(0.5199334546165302 + 0.4081309459110076*gelu(-0.7470183626892886 + -0.6899614176169906*$(x[1]) + -0.2208825424587011*$(x[2])) + 0.6995024782698374*gelu(-0.11062111710238876 + 0.5516588711870343*$(x[1]) + 0.9886940321128486*$(x[2])) + -0.3131098472642604*gelu(-0.28538228970678503 + -0.152956579710676*$(x[1]) + -0.25273256901725194*$(x[2])) + -0.4135863095985335*gelu(0.7652671472719401 + -0.8637339388616145*$(x[1]) + -0.2539724731598083*$(x[2])))) - $q <= 0.0))

                     @objective(m, Min, q)

                     return m

max_stars_repo_name: PSORLab/RSActivationFunctions

max_issues_repo_name: PSORLab/RSActivationFunctions

max_forks_repo_name: PSORLab/RSActivationFunctions


text:
include("solver_preamble.jl")
import ProxSDP
factory = with_optimizer(ProxSDP.Optimizer, verbose=0)
config = MOI.Test.TestConfig(atol=1e-3, rtol=1e-3)
@testset "Linear" begin
    Tests.linear_test(factory, config, ["dsos_horn"])
end
@testset "SOC" begin
    Tests.soc_test(factory, config, ["sdsos_horn", "sdsos_bivariate_quadratic"])
end
@testset "SDP" begin
    Tests.sd_test(factory, config,
                  ["sos_horn",
                   "sos_bivariate_quadratic",
                   "quadratic_infeasible_lyapunov_switched_system",
                   "quartic_infeasible_lyapunov_switched_system"])
end

ext: jl


text:
## Julia 1.1.0

# Integer: (min, max)
println("\n\n---------------- Integer: (min, max) ----------------")
for T in [Int8,Int16,Int32,Int64,Int128,UInt8,UInt16,UInt32,UInt64,UInt128]
    println("$(lpad(T,7)): [$(typemin(T)),$(typemax(T))]")
end

# Integer: overflow
println("\n\n---------------- Integer: overflow ----------------")
max = typemax(Int64)
if (max + 1 == typemin(Int64))
    println("Int64 overflow detected ", max, " + 1 = ", max + 1)
end
println("BigInt(typemax(Int64) + 1) is $(BigInt(typemax(Int64)))")

# Float: zero
println("\n\n---------------- Float: zero ----------------")
println("0.0 == 0.0 is ", 0.0 == 0.0)
println("bitstring(0.0) is ", bitstring(0.0))
println("bitstring(-0.0) is ", bitstring(-0.0))
println("1/Inf is ", 1/Inf)

# Float: (min, max)
println("\n\n---------------- Float: (min, max) ----------------")
for T in [Float32, Float64]
    println("$(lpad(T,7)): [$(typemin(T)), $(typemax(T))], eps: $(eps(T))")
end

# Float: eps, nextfloat
println("\n\n---------------- Float: eps, nextfloat ----------------")
for V in [1e-10, 1e-5, 1.0, 2.0, 10.0, 100.0, 1000.0]
    println("eps($(V)) is $(eps(V)), nextfloat is $(nextfloat(V))")
end

# Operator: exponential
x = 3
println("\n\n---------------- Operator: exponential ----------------")
println("x = $(x), 2(x-1)^2 - 3(x-1) + 1 = $(2(x-1)^2 - 3(x-1) + 1)")
println("x = $(x), 2^2x = $(2^2x)")

max_stars_repo_name: fengyfei/sapphire

max_issues_repo_name: fengyfei/sapphire

max_forks_repo_name: fengyfei/sapphire


text:
"""
    Auto()

Use this in the place of a tangent/cotangent in [`test_frule`](@ref) or
[`test_rrule`](@ref) to have that tangent/cotangent generated automatically based on the
primal. Uses `FiniteDifferences.rand_tangent`.
"""
struct Auto end

"""
    PrimalAndTangent

A struct that represents a primal value paired with its tangent or cotangent.
For conciseness we refer to both tangent and cotangent as "tangent".
"""
struct PrimalAndTangent{P,D}
    primal::P
    tangent::D
end
primal(p::PrimalAndTangent) = p.primal
tangent(p::PrimalAndTangent) = p.tangent

"""
    primal ⊢ tangent

Infix shorthand method to construct a `PrimalAndTangent`.
Enter via `\\vdash` + tab on supporting editors.
"""
const ⊢ = PrimalAndTangent

"""
    auto_primal_and_tangent(primal; rng=Random.GLOBAL_RNG)
    auto_primal_and_tangent(::PrimalAndTangent; rng=Random.GLOBAL_RNG)

Convience constructor for `PrimalAndTangent` where the primal is provided

This function is idempotent. If you pass it a `PrimalAndTangent` it doesn't change it.
"""
auto_primal_and_tangent(primal; rng=Random.GLOBAL_RNG) = primal ⊢ rand_tangent(rng, primal)
auto_primal_and_tangent(both::PrimalAndTangent; kwargs...) = both


max_stars_repo_name: AnAverageHuman/competitive

max_issues_repo_name: AnAverageHuman/competitive

max_forks_repo_name: AnAverageHuman/competitive


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

@testset ExtendedTestSet "`Project`" begin
    tfm = Project(Translation(2, 2))
    bounds = Bounds((50, 50))
    @test projectionbounds(tfm, getprojection(tfm, bounds), bounds) == Bounds((3:52, 3:52))

    image = Image(rand(50, 50))
    @test_nowarn apply(tfm, image)
    timage = apply(tfm, image)
    @test timage.data[1] == image.data[1]
end


text:
__precompile__(true)
module Lifting

@assert VERSION >= v"0.6.0" "Lifting.jl requires Julia v0.6.0 or newer."

# lift types
lift(::Type{T}) where {T<:Nullable} = eltype(T)
lift(::Type{T}) where {T} = T

# lift values
function lift(x::T) where {T<:Nullable}
    @assert !isnull(x) "Cannot lift a null value."
    return get(x)
end

lift(x) = x

# unlift types
unlift(::Type{T}) where {T<:Nullable} = T
unlift(::Type{T}) where {T} = Nullable{T}

# unlift values
unlift(x::T) where {T<:Nullable} = x
unlift(x::T) where {T} = Nullable{T}(x)

# NullableArrays

using NullableArrays

lift(::Type{NullableArray{P,Q}}) where {P,Q} = Array{P,Q}
unlift(::Type{Array{P,Q}}) where {P,Q} = NullableArray{P,Q}

function lift(x::NullableVector{T}) where {T}
	LEN = length(x)
	out = Vector{T}(LEN)
	for i in 1:LEN
		@assert !isnull(x[i]) "Cannot lift a null value."
		out[i] = get(x[i])
	end
	return out
end

unlift(x::Vector{T}) where {T} = NullableArray(x)

export lift, unlift

end # module


text:
"""
Example of how to runs loops of OPFs and plots a series of graphs to facilitate the analisys
of the Transmission Constrained Residual Demand Derivative (TCRDD) approach using the
package MaximumGeneratorProfitTCRDD. The plots allow the user to visualize the actual
behaviour of the system under different ammounts of bids vs the system with approximations
for the calculation of the profit function. The TCRDD algorithm is taken from [1].

[1] L. Xu, R. Baldick and Y. Sutjandra, "Bidding Into Electricity Markets:
    A Transmission-Constrained Residual Demand Derivative Approach," in IEEE
    Transactions on Power Systems, vol. 26, no. 3, pp. 1380-1388, Aug. 2011,
    doi: 10.1109/TPWRS.2010.2083702.
"""
# Packages
using DataFrames
using Dates
using DataStructures
using Ipopt
using MaximumGeneratorProfitTCRDD
using PowerSystems
using PowerSimulations
using Test
using TimeSeries

# Build Test Case
include("data/example_case4/buildup_case4.jl")
include("data/example_case4/c_sys_case4.jl")
sys = c_sys_case4()

# Settings
# Number of OPF points to plot the profit function
div = 500
# Network representation
network = StandardPTDFModel
# Solver for OPF
solver = optimizer_with_attributes(Ipopt.Optimizer)
# Tolerance to identify a change in the TCRDD
tcrdd_tol = 1e-3
# Gap to force a recalculation of the TCRDD
gap = div #No forced recalculation
# Tolerance to identify binding constraints
dual_lines_tol = 1e-1
dual_gen_tol = 1e-1
# Flag to get an extra plot for Example Case 4
examplecase4 = true

# Run and plot the OPF Loops
plots_TCRDD(
    sys;
    div,
    network,
    solver,
    tcrdd_tol,
    gap,
    dual_lines_tol,
    dual_gen_tol,
    examplecase4
)

max_stars_repo_name: invenia/MaximumGeneratorProfitTCRDD

max_issues_repo_name: invenia/MaximumGeneratorProfitTCRDD

max_forks_repo_name: invenia/MaximumGeneratorProfitTCRDD


text:
using Mimi

include("rice2010.jl")
using .Rice2010

m = getrice()
run(m)

explore(m)


text:
# Note: This must be the first file executed in the tests.
#
# The following structure is used to test the problem reported at issue #507. We
# need to define it here because, since `recursivecopy!` and `copy_fields!` are
# generated functions, then the `ArrayInterface.isimmutable` method must be
# defined before the first `using DiffEqBase` in a Julia session.

struct Quaternion{T} <: AbstractVector{T}
    q0::T
    q1::T
    q2::T
    q3::T
end

using ArrayInterface
ArrayInterface.ismutable(::Type{<:Quaternion}) = false
Base.size(::Quaternion) = 4

# https://github.com/JuliaDiffEq/DifferentialEquations.jl/issues/525

using DEDataArrays, OrdinaryDiffEq, StochasticDiffEq
using StaticArrays, RecursiveArrayTools, LinearAlgebra, Test

mutable struct VectorType{T} <: DEDataVector{T}
    x::Vector{T}
    f::T
end

mutable struct MatrixType{T,S} <: DEDataMatrix{T}
    f::S
    x::Matrix{T}
end

mutable struct SimType2{T} <: DEDataVector{T}
  x::Vector{T}
  y::Vector{T}
  u::Vector{T}
end

A = [0.0; 1.0]
B = [1.0  2; 4  3]

a = VectorType{Float64}(copy(A), 1.0)
b = MatrixType{Float64,Float64}(2.0, copy(B))

# basic methods of AbstractArray interface
@test eltype(a) == Float64 && eltype(b) == Float64

# size
@test length(a) == 2 && length(b) == 4
@test size(a) == (2,) && size(b) == (2,2)

# iteration
@test first(a) == 0.0 && first(b) == 1
@test last(a) == 1.0 && last(b) == 3
for (i, (ai, Ai)) in enumerate(zip(a, A))
    @test ai == Ai
end
for (i, (bi, Bi)) in enumerate(zip(b, B))
    @test bi == Bi
end

# indexing
@test eachindex(a) == Base.LinearIndices(A) && eachindex(b) == vec(Base.LinearIndices(B))
@test a[2] == a[end] == 1.0
@test b[3] == b[1,2] == 2 && b[:, 1] == [1; 4]

a[1] = 3; b[2,1] = 1
@test a[1] == 3.0 && b[2] == 1

a[:] = A; b[:, 1] = B[1:2]
@test a.x == A && b.x == B

# simple broadcasts
@test [1//2] .* a == [0.0; 0.5]
@test b ./ 2 ==  [0.5  1.0; 2.0  1.5]
@test b .+ a  == [1.0  2.0; 5.0  4.0]
@test_broken a .+ b  == [1.0  2.0; 5.0  4.0] # Doesn't find the largest

# similar data arrays
a2 = similar(a); b2 = similar(b, (1,4)); b3 = similar(b, Float64, (1, 4))
@test typeof(a2.x) == Vector{Float64} && a2.f == 1.0
@test typeof(b2.x) == Matrix{Float64} && b2.f == 2.0
@test typeof(b3.x) == Matrix{Float64} && b3.f == 2.0

# copy all fields of data arrays
recursivecopy!(a2, a)
recursivecopy!(b2, b)
@test a2.x == A && vec(b2.x) == vec(B)
recursivecopy!(b3, b)

# copy all fields except of field x
a2.f = 3.0; b2.f = -1; b3.f == 0
DEDataArrays.copy_fields!(a, a2)
DEDataArrays.copy_fields!(b, b2)
@test a.f == 3.0 && b.f == -1.0
DEDataArrays.copy_fields!(b, b3)

# create data array with field x replaced by new array
a3 = DEDataArrays.copy_fields([1.0; 0.0], a)
@test a3 == VectorType{Float64}([1; 0], 3.0)

# broadcast assignments
a.f = 0.0
a .= a2 .+ a3
@test a == VectorType{Float64}([1.0; 1.0], 0.0)
@test a == a2 .+ a3
@test (a2 .+ a3) isa VectorType

old_b = copy(b)
b .= b .^ 2 .+ a3
@test b == MatrixType{Int,Float64}(-1.0, [2  5; 16  9])
@test b == old_b .^ 2 .+ a3
@test (b .^ 2 .+ a3) isa MatrixType

# Test ability to use MVectors in the solvers
mutable struct SimWorkspace{T} <: DEDataVector{T}
  x::MVector{2,T}
  a::T
end
s0   = SimWorkspace{Float64}(MVector{2,Float64}(1.0,4.0),1.)
similar(s0,Float64,size(s0))
s1   = SimWorkspace{Float64}(MVector{2,Float64}(2.0,1.0),1.)
s0 .+ s1 == SimWorkspace{Float64}(MVector{2,Float64}(3.0,5.0),1.)

mutable struct SimWorkspace2{T} <: DEDataVector{T}
  x::SVector{2,T}
  a::T
end
s0   = SimWorkspace2{Float64}(SVector{2,Float64}(1.0,4.0),1.)
s1   = SimWorkspace2{Float64}(SVector{2,Float64}(2.0,1.0),1.)
s0 .+ s1 == SimWorkspace2{Float64}(SVector{2,Float64}(3.0,5.0),1.)

# Test `recursivecopy!` in immutable structures derived from `AbstractArrays`.
# See issue #507.
mutable struct SimWorkspace3{T} <: DEDataVector{T}
    x::Vector{T}
    q::Quaternion{T}
end

a = SimWorkspace3([1.0,2.0,3.0], Quaternion(cosd(15), 0.0, 0.0, sind(15)))
b = SimWorkspace3([0.0,0.0,0.0], Quaternion(1.0, 0.0, 0.0, 0.0))

recursivecopy!(b,a)
@test b.x == a.x
@test b.q.q0 == a.q.q0
@test b.q.q1 == a.q.q1
@test b.q.q2 == a.q.q2
@test b.q.q3 == a.q.q3

a = SimWorkspace3([1.0,2.0,3.0], Quaternion(cosd(15), 0.0, 0.0, sind(15)))
b = SimWorkspace3([0.0,0.0,0.0], Quaternion(1.0, 0.0, 0.0, 0.0))

DEDataArrays.copy_fields!(b, a)
@test b.x == [0.0,0.0,0.0]
@test b.q.q0 == a.q.q0
@test b.q.q1 == a.q.q1
@test b.q.q2 == a.q.q2
@test b.q.q3 == a.q.q3

mutable struct SimType{T} <: DEDataVector{T}
    x::Array{T,1}
    f1::T
end

function f(u,p,t) # new out-of-place definition
    SimType([-0.5*u[1] + u.f1,
    -0.5*u[2]],u.f1)
end

function f!(du,u,p,t) # old in-place definition
    du[1] = -0.5*u[1] + u.f1
    du[2] = -0.5*u[2]
end

tstop1 = [5.]
tstop2 = [8.]

function condition(u,t,integrator)
  t in tstop1
end

function condition2(u,t,integrator)
  t in tstop2
end

function affect!(integrator)
  for c in full_cache(integrator)
    c.f1 = 1.5
  end
end

function affect2!(integrator)
  for c in full_cache(integrator)
    c.f1 = -1.5
  end
end

function affect!_oop(integrator)
  integrator.u.f1 = 1.5
end

function affect2!_oop(integrator)
  integrator.u.f1 = 1.5
end

save_positions = (true,true)
cb = DiscreteCallback(condition, affect!, save_positions=save_positions)
save_positions = (false,true)
cb2 = DiscreteCallback(condition2, affect2!, save_positions=save_positions)
cbs = CallbackSet(cb,cb2)

cb_oop = DiscreteCallback(condition, affect!_oop, save_positions=save_positions)
save_positions = (false,true)
cb2_oop = DiscreteCallback(condition2, affect2!_oop, save_positions=save_positions)
cbs_oop = CallbackSet(cb_oop,cb2_oop)

u0 = SimType([10.0;10.0], 0.0)

prob_inplace = ODEProblem(f!,u0,(0.0,10.0))
prob = ODEProblem(f,u0,(0.0,10.0))

tstop = [5.;8.]

sol = solve(prob_inplace,Tsit5(),callback = cbs, tstops=tstop)
sol = solve(prob,Tsit5(),callback = cbs_oop, tstops=tstop)

# https://github.com/JuliaDiffEq/DifferentialEquations.jl/issues/336

AA = SMatrix{2,2}([0 1;
                   0 0])

mutable struct MyStruct{T} <: DEDataVector{T}
    x::MVector{2,T}
    a::SVector{2,T}
end

function dyn(du,u,t,p)
    u.a = SVector{2}(0.0, 0.1)
    du .= AA*u.x + u.a

    return nothing
end

u0   = MyStruct(MVector{2}(0.,0.), SVector{2}(0.,0.))
prob = ODEProblem(dyn, u0, (0.,10.))

@test copy(u0) isa MyStruct
@test zero(u0) isa MyStruct
@test similar(u0) isa MyStruct
@test similar(u0,Float64) isa MyStruct
@test similar(u0,Float64,size(u0)) isa MyStruct

sol  = solve(prob, Tsit5())

# https://github.com/JuliaDiffEq/StochasticDiffEq.jl/issues/247

function f(du,u,p,t)
    du[1] = -0.5*u[1] + u.f1
    du[2] = -0.5*u[2]
end

tstop1 = [5.]
tstop2 = [8.]

function condition(u,t,integrator)
  t in tstop1
end

function condition2(u,t,integrator)
  t in tstop2
end

function affect!(integrator)
  for c in full_cache(integrator)
    c.f1 = 1.5
  end
end

function affect2!(integrator)
  for c in full_cache(integrator)
    c.f1 = -1.5
  end
end

save_positions = (true,true)

cb = DiscreteCallback(condition, affect!, save_positions=save_positions)

save_positions = (false,true)

cb2 = DiscreteCallback(condition2, affect2!, save_positions=save_positions)

cbs = CallbackSet(cb,cb2)

u0 = SimType([10.0;10.0], 0.0)
prob = ODEProblem(f,u0,(0.0,10.0))

tstop = [5.;8.]

# here the new part

function g(du,u,p,t)
  du[1] = 1.0
  du[2] = 1.2
end

dt = 1/2^4

prob2 = SDEProblem(f,g,u0,(0.0,10.0))

# this creates an error

sol = solve(prob2,callback = cbs,tstops=tstop,EM(),dt=dt,saveat=collect(8:0.1:10))

# https://github.com/JuliaDiffEq/DiffEqBase.jl/issues/327

struct SimulationState{T,U} <: DEDataArrays.DEDataArray{T, 1}
    x::Vector{T}
    u::U
end

#function Base.convert(::Type{SimulationState{T, U}}, a::AbstractArray{T}) where {T,U}
#    SimulationState(a, zero(eltype(a)))
#end

function open_loop(state, parameters, time)
    v = state[1] * -0.1 + state.u
    return SimulationState([v],state.u)
end

initial_conditions = SimulationState([10.0], 0.0)
time_span = (0.0, 20.0)
ode_prob = ODEProblem(open_loop, initial_conditions, time_span, nothing)

sol = solve(ode_prob, Tsit5(), reltol=1e-8, abstol=1e-8)

@testset "DEDataVector" begin
  f = function (du,u,p,t)
    du[1] = -0.5*u[1] + u.f1
    du[2] = -0.5*u[2]
  end

  tstop1 = [5.]
  tstop2 = [8.]
  tstop = [5.;8.]

  condition = function (u,t,integrator)
    t in tstop1
  end

  affect! = function (integrator)
    for c in OrdinaryDiffEq.full_cache(integrator)
      c.f1 = 1.5
    end
  end

  save_positions = (true,true)

  cb = DiscreteCallback(condition, affect!; save_positions=save_positions)

  condition2 = function (u,t,integrator)
    t in tstop2
  end

  affect2! = function (integrator)
    for c in OrdinaryDiffEq.full_cache(integrator)
      c.f1 = -1.5
    end
  end

  save_positions = (true,true)

  cb2 = DiscreteCallback(condition2, affect2!, save_positions=save_positions)

  cbs = CallbackSet(cb,cb2)

  u0 = SimType{Float64}([10;10], 0.0)
  prob = ODEProblem(f,u0,(0.0,10.0))
  sol = solve(prob,Tsit5(),callback = cbs, tstops=tstop)

  sol(1.5:0.5:2.5)

  @test [sol[i].f1 for i in eachindex(sol)] == [fill(0., 9);1.5*ones(5);-1.5*ones(4)]

  A = Matrix(Diagonal([0.3,0.6,0.9]))
  B = transpose([1 2 3])
  C = [1/3 1/3 1/3]

  function mysystem(t,x,dx,p,u)
      ucalc = u(x,p,t)
      x.u = ucalc
      x.y = C*x.x
      dx .= A*x.x + B*x.u
  end

  input = (x,p,t)->(1*one(t)≤t≤2*one(t) ? [one(t)] : [zero(t)])
  prob = DiscreteProblem((dx,x,p,t)->mysystem(t,x,dx,p,input), SimType2(fill(0., 3), fill(0., 1), fill(0., 1)), (0//1,4//1))
  sln = solve(prob, FunctionMap(scale_by_time=false), dt = 1//10)

  u1 = [sln[idx].u for idx in 1:length(sln)]
  u2 = [sln(t).u for t in range(0,stop=4,length=41)]
  @test any(x->x[1]>0, u1)
  @test any(x->x[1]>0, u2)

  sln = solve(prob, FunctionMap(scale_by_time=true), dt = 1//10)

  u1 = [sln[idx].u for idx in 1:length(sln)]
  u2 = [sln(t).u for t in range(0,stop=4,length=41)]
  @test any(x->x[1]>0, u1)
  @test any(x->x[1]>0, u2)

  sln = solve(prob, Euler(), dt = 1//10)

  @test u1 == [sln[idx].u for idx in 1:length(sln)] # Show that discrete is the same
  u1 = [sln[idx].u for idx in 1:length(sln)]
  u2 = [sln(t).u for t in range(0,stop=4,length=41)]
  @test any(x->x[1]>0, u1)
  @test any(x->x[1]>0, u2)

  sln = solve(prob, DP5(), dt = 1//10, adaptive=false)

  u1 = [sln[idx].u for idx in 1:length(sln)]
  u2 = [sln(t).u for t in range(0,stop=4,length=41)]
  @test any(x->x[1]>0, u1)
  @test any(x->x[1]>0, u2)
end

######################
# DEDataMatrix
mutable struct SimTypeg{T,T2,N} <: DEDataArray{T,N}
  x::Array{T,N} # two dimensional
  f1::T2
end

@testset "DEDataMatrix" begin
  tstop1 = [10.0]
  tstop2 = [300.]

  function mat_condition(u,t,integrator)
    t in tstop1
  end

  function mat_condition2(u,t,integrator)
    t in tstop2
  end

  function mat_affect!(integrator)
    for c in OrdinaryDiffEq.full_cache(integrator)
      c.f1 = +1.0
    end
  #  integrator.u[1,1] = 0.001
  end

  function mat_affect2!(integrator)
    for c in OrdinaryDiffEq.full_cache(integrator)
      c.f1 = 0.0
    end
  end

  save_positions = (true,true)
  cb = DiscreteCallback(mat_condition, mat_affect!, save_positions=save_positions)
  save_positions = (false,true)
  cb2 = DiscreteCallback(mat_condition2, mat_affect2!, save_positions=save_positions)
  cbs = CallbackSet(cb,cb2)

  function sigmoid(du,u,p,t)
    du[1,1] = 0.01*u[1,1]*(1-u[1,1]/20)
    du[1,2] = 0.01*u[1,2]*(1-u[1,2]/20)
    du[2,1] = 0.01*u[2,1]*(1-u[2,1]/20)
    du[2,2] = u.f1*du[1,1]
  end

  u0 = SimTypeg(fill(0.00001,2,2),0.0)
  tspan = (0.0,3000.0)
  prob = ODEProblem(sigmoid,u0,tspan)

  tstop =[tstop1;tstop2]
  sol = solve(prob,Tsit5(),callback = cbs, tstops=tstop)
  sol = solve(prob,Rodas4(),callback = cbs, tstops=tstop)
  sol = solve(prob,Kvaerno3(),callback = cbs, tstops=tstop)
  sol = solve(prob,Rodas4(autodiff=false),callback = cbs, tstops=tstop)
  sol = solve(prob,Kvaerno3(autodiff=false),callback = cbs, tstops=tstop)
end

mutable struct CtrlSimTypeScalar{T,T2} <: DEDataVector{T}
    x::Vector{T}
    ctrl_x::T2 # controller state
    ctrl_u::T2 # controller output
end

@testset "Interpolation Sides" begin
  function f(x,p,t)
      x.ctrl_u .- x
  end
  function f(dx,x,p,t)
      dx[1] = x.ctrl_u - x[1]
  end

  ctrl_f(e, x) = 0.85(e + 1.2x)

  x0 = CtrlSimTypeScalar([0.0], 0.0, 0.0)
  prob = ODEProblem{false}(f, x0, (0.0, 8.0))

  function ctrl_fun(int)
      e = 1 - int.u[1] # control error

      # pre-calculate values to avoid overhead in iteration over cache
      x_new = int.u.ctrl_x + e
      u_new = ctrl_f(e, x_new)

      # iterate over cache...
      if DiffEqBase.isinplace(int.sol.prob)
        for c in full_cache(int)
            c.ctrl_x = x_new
            c.ctrl_u = u_new
        end
      end
  end

  integrator = init(prob, Tsit5())

  # take 8 steps with a sampling time of 1s
  Ts = 1.0
  for i in 1:8
      ctrl_fun(integrator)
      DiffEqBase.step!(integrator,Ts,true)
  end

  sol = integrator.sol
  @test [u.ctrl_u for u in sol.u[2:end]] == [sol(t).ctrl_u for t in sol.t[2:end]]

  prob = ODEProblem{true}(f, x0, (0.0, 8.0))

  integrator = init(prob, Rosenbrock23())

  # take 8 steps with a sampling time of 1s
  Ts = 1.0
  for i in 1:8
      ctrl_fun(integrator)
      DiffEqBase.step!(integrator,Ts,true)
  end

  sol = integrator.sol
  @test [u.ctrl_u for u in sol.u[2:end]] == [sol(t).ctrl_u for t in sol.t[2:end]]
end


text:
export BoxConstrainedDifferentiableModel
using NLopt, SparseInverseProblems.Util
# A simple forward model with box constrained parameters.
# Assumes differentiablity of the forward operator.
#
# For concrete examples see examples/smi or examples/sysid.
#
abstract BoxConstrainedDifferentiableModel <: ForwardModel

#compute the measurement model
psi(model :: BoxConstrainedDifferentiableModel, theta :: Vector{Float64}) =
  error("psi not implemented for model $(typeof(model)).")

#compute the jacobian of the forward model
dpsi(model :: BoxConstrainedDifferentiableModel, theta :: Vector{Float64}) =
    error("dpsi not implemented for model $(typeof(model)).")

# Initial starting point for continuous optimization for the FW step.
# Should return a good guess for $\arg\min_\theta \langle \psi(theta), v \rangle.$
# Often computed using a grid.
getStartingPoint(model :: BoxConstrainedDifferentiableModel, v :: Vector{Float64}) =
  error("getStartingPoint not implemented for model $(typeof(model)).")

# Box constraints on the parameters.
# Returns a tuple of two vectors : lower bounds and upper bounds.
parameterBounds(model :: BoxConstrainedDifferentiableModel) =
  error("parameterBounds not implemented for model $(typeof(model)).")

#default implementation of phi in terms of psi --- you should probably overwride this!
phi(model :: BoxConstrainedDifferentiableModel, parameters :: Matrix{Float64}, weights :: Vector{Float64}) =
    sum([weights[i]*psi(model, vec(parameters[:,i])) for i = 1:size(parameters,2)])

# Computes $\nabla_{\theta_1,\ldots,\theta_k} \sum_i \langle w_i \psi(\theta_i), v \rangle.$
# Returns a matrix with the same shape as thetas (which is p by k).
# You should probably overwride this!
computeGradient(model :: BoxConstrainedDifferentiableModel, weights :: Vector{Float64},
  thetas :: Matrix{Float64}, v :: Vector{Float64}) =
  hcat([weights[i]*vec(dpsi(model, vec(thetas[:,i]))'*v) for i = 1:size(thetas,2)]...)

# Sets the parameters for the continuous optimizer.
# Can be overwridden.
function initializeOptimizer!(model :: BoxConstrainedDifferentiableModel, opt :: Opt)
  ftol_abs!(opt, 1e-6)
  xtol_rel!(opt, 0.0)
  maxeval!(opt, 200)
end

# Default implementation. Uses NLopt to do continuous optimization.
function lmo(model :: BoxConstrainedDifferentiableModel, v :: Vector{Float64})
  lb,ub = parameterBounds(model)
  initial_x = getStartingPoint(model, v)
  p = length(lb)

  function f_and_g!(point :: Vector{Float64}, gradient_storage :: Vector{Float64})
    output = phi(model,reshape(point,p,1), [1.0])
    ip = dot(output,v)
    s = sign(ip)
    gradient_storage[:] = -s*computeGradient(model, [1.0],reshape(point,length(point),1), v)
    return -s*ip
  end

  opt = Opt(:LD_MMA, p)
  initializeOptimizer!(model, opt)
  min_objective!(opt, f_and_g!)
  lower_bounds!(opt, lb)
  upper_bounds!(opt, ub)
  (optf,optx,ret) = optimize(opt, initial_x)
  return (optx,optf)
end

immutable SupportUpdateProblem
  nPoints :: Int64
  p :: Int64
  s :: BoxConstrainedDifferentiableModel
  y :: Vector{Float64}
  w :: Vector{Float64}
  lossFn :: Loss
end

# Default implementation. Uses NLopt to do continuous optimization.
function localDescent(s :: BoxConstrainedDifferentiableModel, lossFn :: Loss, thetas ::Matrix{Float64}, w :: Vector{Float64}, y :: Vector{Float64})
  lb,ub = parameterBounds(s)
  nPoints = size(thetas,2)
  p = size(thetas,1)
  su = SupportUpdateProblem(nPoints,p,s,y,w,lossFn)
  f_and_g!(x,g) = localDescent_f_and_g!(x,g,su)
  opt = Opt(NLopt.LD_MMA, length(thetas))
  initializeOptimizer!(s, opt)
  min_objective!(opt, f_and_g!)
  lower_bounds!(opt, vec(repmat(lb,1,nPoints)))
  upper_bounds!(opt, vec(repmat(ub,1,nPoints)))
  (optf,optx,ret) = optimize(opt, vec(thetas))
  return reshape(optx,p,nPoints)
end

function localDescent_f_and_g!(points :: Vector{Float64}, gradient_storage :: Vector{Float64}, s :: SupportUpdateProblem)
  points = reshape(points,s.p,s.nPoints)
  output = phi(s.s,points,s.w)
  residual = output - s.y
  l,v_star = loss(s.lossFn,residual)
  gradient_storage[:] = computeGradient(s.s, s.w, points, v_star)
  return l
end

# Only for squared loss
solveFiniteDimProblem(model :: BoxConstrainedDifferentiableModel, loss :: LSLoss, thetas :: Matrix{Float64}, y :: Vector{Float64}, tau :: Float64) = nnlasso(hcat([psi(model, vec(thetas[:,i])) for i = 1:size(thetas,2)]...),y,tau)


text:
using Test

import RandomizedPreconditioners
const RP = RandomizedPreconditioners

using LinearAlgebra
using Random

include("utils.jl")
include("sketch.jl")
include("preconditioner.jl")
include("rangefinder.jl")


text:
const my_dir = @__DIR__
const pluto_dir = joinpath(my_dir, "..", "..")

pushfirst!(LOAD_PATH, pluto_dir)

#
include(joinpath(my_dir, "PlutoRunner.jl"))
#

let
    i = findfirst(isequal(pluto_dir), LOAD_PATH)
    if i !== nothing
        deleteat!(LOAD_PATH, i)
    end
end


text:
module EquationUtilsModule

import ..CoreModule: CONST_TYPE, Node, copy_node, Options

# for symbolic constraints
using PyCall
const _sympy = PyNULL()

# Count the operators, constants, variables in an equation
function count_nodes(tree::Node)::Int
    if tree.degree == 0
        return 1
    elseif tree.degree == 1
        return 1 + count_nodes(tree.l)
    else
        return 1 + count_nodes(tree.l) + count_nodes(tree.r)
    end
end

# Count the max depth of a tree
function count_depth(tree::Node)::Int
    if tree.degree == 0
        return 1
    elseif tree.degree == 1
        return 1 + count_depth(tree.l)
    else
        return 1 + max(count_depth(tree.l), count_depth(tree.r))
    end
end

# Count the number of unary operators in the equation
function count_unary_operators(tree::Node)::Int
    if tree.degree == 0
        return 0
    elseif tree.degree == 1
        return 1 + count_unary_operators(tree.l)
    else
        return 0 + count_unary_operators(tree.l) + count_unary_operators(tree.r)
    end
end

# Count the number of binary operators in the equation
function count_binary_operators(tree::Node)::Int
    if tree.degree == 0
        return 0
    elseif tree.degree == 1
        return 0 + count_binary_operators(tree.l)
    else
        return 1 + count_binary_operators(tree.l) + count_binary_operators(tree.r)
    end
end

# Count the number of operators in the equation
function count_operators(tree::Node)::Int
    return count_unary_operators(tree) + count_binary_operators(tree)
end

# Count the number of constants in an equation
function count_constants(tree::Node)::Int
    if tree.degree == 0
        if tree.constant
            return 1
        else
            return 0
        end
    elseif tree.degree == 1
        return 0 + count_constants(tree.l)
    else
        return 0 + count_constants(tree.l) + count_constants(tree.r)
    end
end

"""
Compute the complexity of a tree.

By default, this is the number of nodes in a tree.
However, it could use the custom settings in options.complexity_mapping
if these are defined.
"""
function compute_complexity(tree::Node, options::Options)::Int
    if options.complexity_mapping.use
        return round(Int, _compute_complexity(tree, options))
    else
        return count_nodes(tree)
    end
end

function _compute_complexity(
    tree::Node, options::Options{A,B,dA,dB,C,complexity_type}
)::complexity_type where {A,B,dA,dB,C,complexity_type<:Real}
    if tree.degree == 0
        if tree.constant
            return options.complexity_mapping.constant_complexity
        else
            return options.complexity_mapping.variable_complexity
        end
    elseif tree.degree == 1
        return (
            options.complexity_mapping.unaop_complexities[tree.op] +
            _compute_complexity(tree.l, options)
        )
    else # tree.degree == 2
        return (
            options.complexity_mapping.binop_complexities[tree.op] +
            _compute_complexity(tree.l, options) +
            _compute_complexity(tree.r, options)
        )
    end
end

# Get all the constants from a tree
function get_constants(tree::Node)::AbstractVector{CONST_TYPE}
    if tree.degree == 0
        if tree.constant
            return [tree.val]
        else
            return CONST_TYPE[]
        end
    elseif tree.degree == 1
        return get_constants(tree.l)
    else
        both = [get_constants(tree.l), get_constants(tree.r)]
        return [constant for subtree in both for constant in subtree]
    end
end

# Set all the constants inside a tree
function set_constants(tree::Node, constants::AbstractVector{T}) where {T<:Real}
    if tree.degree == 0
        if tree.constant
            tree.val = convert(CONST_TYPE, constants[1])
        end
    elseif tree.degree == 1
        set_constants(tree.l, constants)
    else
        numberLeft = count_constants(tree.l)
        set_constants(tree.l, constants)
        set_constants(tree.r, constants[(numberLeft + 1):end])
    end
end

## Assign index to nodes of a tree
# This will mirror a Node struct, rather
# than adding a new attribute to Node.
mutable struct NodeIndex
    constant_index::Int  # Index of this constant (if a constant exists here)
    l::NodeIndex
    r::NodeIndex

    NodeIndex() = new()
end

function index_constants(tree::Node)::NodeIndex
    return index_constants(tree, 0)
end

function index_constants(tree::Node, left_index::Int)::NodeIndex
    index_tree = NodeIndex()
    index_constants(tree, index_tree, left_index)
    return index_tree
end

# Count how many constants to the left of this node, and put them in a tree
function index_constants(tree::Node, index_tree::NodeIndex, left_index::Int)
    if tree.degree == 0
        if tree.constant
            index_tree.constant_index = left_index + 1
        end
    elseif tree.degree == 1
        index_tree.constant_index = count_constants(tree.l)
        index_tree.l = NodeIndex()
        index_constants(tree.l, index_tree.l, left_index)
    else
        index_tree.l = NodeIndex()
        index_tree.r = NodeIndex()
        index_constants(tree.l, index_tree.l, left_index)
        index_tree.constant_index = count_constants(tree.l)
        left_index_here = left_index + index_tree.constant_index
        index_constants(tree.r, index_tree.r, left_index_here)
    end
end

# python function to check if expression is monontonically increasing
py"""
import sympy
def is_monotonic_increasing_test(expr, interval, var):

    # constant value never decreases
    if expr.is_constant():
        return True

    # get critical points as list
    turning_points = list(sympy.solveset(expr.diff(var), var, interval))
    turning_points.sort()
    # failed to find critical points
    # there could be 0 or infinite...
    if (turning_points == []):
        # fall back to simpler increasing function
        return bool(1 if (expr.limit(var, interval.end) - expr.limit(var, interval.start)) >= 0 else 0)
    increasing = 1
    # turn to false if interval from start of main interval to first critical point not increasing
    increasing = min(increasing, (1 if (expr.limit(var, turning_points[0]) - expr.limit(var, interval.start)) >= 0 else 0))
    # check intervals between all critical points
    for i in range(len(turning_points)-1):
        thisPoint = turning_points[i]
        nextPoint = turning_points[i+1]
        increasing = min(increasing, (1 if (expr.limit(var, nextPoint) - expr.limit(var, thisPoint)) >= 0 else 0))
        #increasing = min(increasing, sympy.is_increasing(expr, sympy.Interval(thisPoint, nextPoint, false, false), var))
    # check last interval
    increasing = min(increasing, (1 if (expr.limit(var, interval.end) - expr.limit(var, turning_points[-1])) >= 0 else 0))
    #increasing = min(increasing, sympy.is_increasing(expr, sympy.Interval(turning_points[-1], interval.end, false, false), var))
    return bool(increasing)
"""

function thermoConstraints(expr::PyObject, var::PyObject)

    # good way to check if lib needs to be loaded from PyCall/serialize.jl
    sympy = ispynull(_sympy) ? copy!(_sympy, pyimport_conda("sympy", "sympy")) : _sympy

    results = [true, true, true]

    # Axiom 1: the expr needs to pass through the origin
    try
        if sympy.limit(expr, var, 0, "+") != 0
            #println("constraint 1")
            results[1] = false
        end
    catch error
        #println(error)
        #println("SymPy cannot evaluate Axiom 1")
        results[1] = false
    end
        # Axiom 2: the expr needs to converge to Henry's Law at zero pressure
    try
        if (sympy.limit(sympy.diff(expr, var), var, 0) == sympy.oo
            || sympy.limit(sympy.diff(expr, var), var, 0) == -sympy.oo
            || sympy.limit(sympy.diff(expr, var), var, 0) == 0)
            #println("constraint 2")
            results[2] = false
        end
    catch error
        #println(error)
        #println("SymPy cannot evaluate Axiom 2")
        results[2] = false
    end

    # Axiom 3: the expr must be strictly increasing as pressure increases
    try
        # use custom function because sympy doesn't work as expected
        if not(is_monotonic_increasing_test(expr, sympy.Interval(0,sympy.oo), var))
            #println("constraint 3")
            results[3] = false
        end
    catch error
        #println("SymPy cannot evaluate Axiom 3")
        results[3] = false
    end

    return results
end

end


text:
function _prep_bound_err(ν, T, ν_indices, func, ν_bound, temperature_bound)
    signs = "-+"
    fmt(x) = isinf(x) ? "$(signs[1+(x>0)])∞" : "$x"

    if !contained(T, temperature_bound)
        DomainError(T, ("$func: invalid temperature. It should lie between "*
                        "$(temperature_bound.lower) & $(temperature_bound.upper)"))
    else
        bad_ν = (1 in ν_indices) ? ν[last(ν_indices) + 1] : ν[1]
        DomainError(bad_ν, ("$func: invalid frequency. It should lie between "*
                            "$(fmt(ν_bound.lower)) Hz & $(fmt(ν_bound.upper)) Hz"))
    end
end

"""
    bounds_checked_absorption(func; ν_bound, temp_bound)

Constructs a wrapped function that implements bounds checking and extrapolation

# Parameters
- `func`: a function that has a signature `f(ν::Real, T::Real, args...)::Real`, where `ν` is
  frequency (in Hz) and `T` is temperature (in K)
- `ν_bound::Interval`: Interval of frequencies (in Hz) over which `func` is valid.
- `temp_bound::Interval`: Interval of temperatures (in K) over which `func` is valid.

The resulting function will have this signature:
```
wrapped_func(ν::AbstractVector{<:Real}, T::Real, args...; kwargs...)
```
# Wrapped Function Parameters
- `ν::AbstractVector{<:Real}`: sorted vector of frequencies (in Hz)
- `T::Real`: temperature (in K)
- `args...`: function-specific arguments

For a description for the `kwargs...`, see [Continuum Absorption Kwargs](@ref).
"""
function bounds_checked_absorption(func; ν_bound = Interval(0, Inf), temp_bound = Interval(0, Inf))
    # if generating the inner function in this way involves too much overhead, we could:
    # - convert this outer function to a macro.
    # - perform the boundary checks inside of total_continuum_absorption

    # It's probably unnecessary to let T be a Vector since the outputs should be contiguous in ν to
    # fast interpolation. Broadcasting over T would be fine despite the fact that it entails
    # redundant bounds checking because real problems have only order ~100 T values
    function wrapped_func(ν::AbstractVector{<:Real}, T::Real, args...; error_oobounds::Bool = false,
                          out_α::Union{Nothing,AbstractVector} = nothing)
        @assert issorted(ν, rev = first(ν) > last(ν)) "ν should be sorted"

        α_type = promote_type(eltype(ν), typeof(T), typeof.(args)...)
        out_α = isnothing(out_α) ? zeros(α_type, length(ν)) : out_α
        @assert (eltype(out_α) == α_type) && (length(out_α) == length(ν))

        # find indices where out_α can be updated from arguments that are in bounds
        idx = (!contained(T, temp_bound)) ? (1:0) : contained_slice(ν, ν_bound)

        # handle any out-of-bounds arguments. If !error_oobounds, assume that there is 0 absorption
        # for these args. Otherwise, throw error
        if (idx != 1:length(ν)) && error_oobounds  # (it's okay for this to be slow)
            throw(_prep_bound_err(ν, T, idx, func, ν_bound, temp_bound))
        end

        view(out_α, idx) .+= func.(ν[idx], T, args...)
        out_α
    end

    wrapped_func
end


text:
push!(LOAD_PATH,"../src/")
using Documenter, MolHandler
makedocs(
    sitename="MolHandler.jl",
	  pages = [
        "Top" => "index.md",
        "Functions" => "functions.md",
        "Structs" => "structs.md"
    ]
)
deploydocs(
    repo = "github.com/yutakasi634/MolHandler.jl.git",
)


text:
## --- Images.jl

   using Colors: Color

   cmap = resize_colormap(viridis, 10)
   @test length(cmap) == 10
   @test isa(cmap, Array{<:Color,1})

   matrix = rand(10,10)
   # Specifiying limits
   img1 = imsc(matrix, viridis, 0, 1)
   @test isa(img1, Array{<:Color,2})
   img2 = imsci(matrix, viridis, 0, 1)
   @test isa(img2, AbstractArray{<:Color,2})
   @test all(img1 .== img2)

   # Auto-ranging
   img1 = imsc(matrix, viridis)
   @test isa(img1, Array{<:Color,2})
   img2 = imsci(matrix, viridis)
   @test isa(img2, AbstractArray{<:Color,2})
   @test all(img1 .== img2)

   # Other
   #@test display(colormaps) != NaN

## ---

max_stars_repo_name: brenhinkeller/StatGeochemBase

max_issues_repo_name: brenhinkeller/StatGeochemBase

max_forks_repo_name: brenhinkeller/StatGeochemBase


text:
__precompile__()

module Integrators

    using ..Equations
    using ..Solvers
    using ..Utils

    import ..Solvers.function_stages!

    include("utils/macro_utils.jl")

    export Tableau, TableauRK, TableauERK, TableauDIRK, TableauFIRK, TableauSIRK,
           TableauEPRK, TableauIPRK, TableauSARK, TableauSPARK, TableauGLM,
           showTableau, writeTableauToFile, readTableauERKFromFile

    include("integrators/tableaus.jl")

    export Timeseries, compute_timeseries!

    include("integrators/timeseries.jl")

    export Solution, SolutionODE, SolutionPODE, SolutionDAE, SolutionPDAE,
           reset!, set_initial_conditions!, compute_time_steps!,
           createHDF5, writeSolutionToHDF5

    include("integrators/solutions.jl")
    include("integrators/solutions_ode.jl")
    include("integrators/solutions_pode.jl")
    include("integrators/solutions_dae.jl")
    include("integrators/solutions_pdae.jl")

    export InitialGuess, evaluate, initialize!, update!

    export Integrator, IntegratorERK, IntegratorDIRK, IntegratorFIRK, IntegratorSIRK,
           IntegratorEPRK, IntegratorIPRK, IntegratorSARK, IntegratorSPARK,
           integrate, integrate!, function_stages!

    include("integrators/integrators.jl")
    include("integrators/integrators_erk.jl")
    include("integrators/integrators_eprk.jl")

    include("integrators/initial_guess.jl")

    include("integrators/integrators_dirk.jl")
    include("integrators/integrators_firk.jl")
    include("integrators/integrators_sirk.jl")
    include("integrators/integrators_iprk.jl")
    include("integrators/integrators_sark.jl")
    include("integrators/integrators_spark.jl")

end


text:
#=
    Same as f4.jl, but with tracing.
=#

mutable struct Tracer
    isredundant_iter::Vector{Int}
    pairset_size::Int
    basis_ntotal::Int
    ready::Bool

    function Tracer()
        new(Vector{Int}(undef, 0), 0, 0, false)
    end
end

#------------------------------------------------------------------------------

#=
    Input ivariants:
        - ring is set, and ring.ch == basis.ch, and ring ~ ht
        - divmasks in ht are set
        - basis is filled so that
            basis.ntotal = actual number of elements
            basis.ndone  = 0
            basis.nlead  = 0

    Output invariants:
        - basis.ndone == basis.ntotal == basis.nlead
        - basis.gens and basis.coeffs of size basis.ndone
        - basis elements are sorted increasingly wrt ordering on lead elements
        - divmasks in basis are filled and coincide to divmasks in hashtable

=#
function f4trace!(ring::PolyRing,
             basis::Basis{Coefftype},
             tracer::Tracer,
             pairset,
             ht,
             reduced,
             linalg,
             rng) where {Coefftype<:Coeff}

    # print("input: $(basis.ntotal) gens, $(ring.nvars) vars. ")

    @assert ring.ch == basis.ch
    @assert ring.ord == ht.ord && ring.nvars == ht.nvars && ring.explen == ht.explen
    # @error "hashtable divmasks"
    # println(ht.exponents[2:ht.load])
    # println("###########")
    # println(map(x->bitstring(x.divmask), ht.hashdata[2:ht.load]))
    # println("###########")
    @assert basis.ndone == 0

    # matrix storing coefficients in rows
    # wrt columns representing the current monomial basis
    matrix = initialize_matrix(ring, Coefftype)

    # initialize hash tables for update and symbolic preprocessing steps
    update_ht  = initialize_secondary_hash_table(ht)
    symbol_ht  = initialize_secondary_hash_table(ht)

    # a set to store critical pairs of polynomials to be reduced
    # if tracer.ready
    #     pairset_initial_size = tracer.pairset_size
    # else
    #     pairset_initial_size = 64
    # end
    # pairset = initialize_pairset(initial_size=pairset_initial_size)

    # @warn "ht initialily" ht.load ht.size

    # makes basis fields valid,
    # does not copy,
    # checks for redundancy of new elems
    plcm = Vector{Int}(undef, 0)
    if tracer.ready
        resize!(plcm, tracer.basis_ntotal + 1)
    end

    # @error "" length(pairset.pairs) length(plcm)

    pairset_size = update!(pairset, basis, ht, update_ht, plcm)
    if !tracer.ready
        tracer.pairset_size = pairset_size
    end

    # @warn "tracer info" tracer.ready tracer.pairset_size

    # @warn "ht update" ht.load ht.size

    d = 0
    # while there are pairs to be reduced
    while !isempty(pairset)
        d += 1
        @debug "F4 ITER $d"
        @debug "Available $(pairset.load) pairs"

        # TODO: learn, and select\discard S-polynomials
        if tracer.ready && tracer.isredundant_iter[d] == 0
            select_normal_and_discard!(pairset, basis, matrix, ht, symbol_ht)
            matrix    = initialize_matrix(ring, Coefftype)
            symbol_ht = initialize_secondary_hash_table(ht)
            continue
        end

        # selects pairs for reduction from pairset following normal strategy
        # (minimal lcm degrees are selected),
        # and puts these into the matrix rows
        select_normal!(pairset, basis, matrix, ht, symbol_ht)

        # @warn "ht select" ht.load ht.size

        symbolic_preprocessing!(basis, matrix, ht, symbol_ht)
        # symbolic_preprocessing_relaxed!(basis, matrix, ht, symbol_ht)
        @debug "Matrix of size $((matrix.nrows, matrix.ncols)), density TODO"

        # @warn "ht symbolic" ht.load ht.size

        # reduces polys and obtains new potential basis elements
        reduction!(basis, matrix, ht, symbol_ht, linalg, rng)
        @debug "Matrix reduced, density TODO"

        if !tracer.ready
            push!(tracer.isredundant_iter, 0)
            if matrix.npivots != 0
                tracer.isredundant_iter[end] = 1
            end
        end

        #=
        if tracer.ready && length(tracer.isredundant_iter) > d
            if tracer.isredundant_iter[d + 1] == 0
                matrix    = initialize_matrix(ring, Coefftype)
                symbol_ht = initialize_secondary_hash_table(ht)
                continue
            end
        end
        =#

        # update the current basis with polynomials produced from reduction,
        # does not copy,
        # checks for redundancy
        pairset_size = update!(pairset, basis, ht, update_ht, plcm)
        if !tracer.ready
            tracer.pairset_size = max(pairset_size, tracer.pairset_size)
        end
        # @warn "ht update" ht.load ht.size

        # TODO: is this okay hm ?
        # to be changed
        # TODO: clean hashtable
        matrix    = initialize_matrix(ring, Coefftype)
        symbol_ht = initialize_secondary_hash_table(ht)
        # clear symbolic hashtable
        # clear matrix

        if d > 10000
            @error "Something is probably wrong in f4. Please submit an issue."
            break
        end
    end

    tracer.ready = true
    tracer.basis_ntotal = basis.ntotal

    # remove redundant elements
    filter_redundant!(basis)

    if reduced
        reducegb_f4!(basis, matrix, ht, symbol_ht)
    end

    standardize_basis!(basis, ht, ht.ord)

    # assertion
    #=
    println("#####")
    println(basis.lead)
    println(map(x->bitstring(x.divmask), ht.hashdata[2:ht.load]))
    println("#####")
    =#
    #=
    for i in 1:basis.nlead
        #=
        println(basis.lead[i])
        println(basis.gens[i][1])
        println(ht.exponents[basis.gens[i][1]])
        println(UInt32(ht.hashdata[basis.gens[i][1]].divmask))
        =#
        @assert basis.lead[i] == ht.hashdata[basis.gens[i][1]].divmask
    end
    =#
    # println("HT: $(ht.load)/$(ht.size)")

    nothing
end


text:
import ZipFile

export NpzReader, NpzWriter, readnpz, writenpz

"""
    NpzReader(io::IO; own=false)
    NpzReader(fn::AbstractString)

Open the given .npz file for reading.

- `keys(reader)` returns the available keys.
- `readnpy(reader, key)` or `reader[key]` returns the array at the given key.
- `close(reader)` closes the underlying IO stream.
"""
struct NpzReader
    zipreader::ZipFile.Reader
    items::Dict{String,ZipFile.ReadableFile}
    function NpzReader(zr::ZipFile.Reader)
        items = Dict{String,ZipFile.ReadableFile}()
        for file in zr.files
            if endswith(file.name, ".npy")
                items[file.name[1:end-4]] = file
            else
                @debug "ignoring .npz entry" file.name
            end
        end
        new(zr, items)
    end
end

NpzReader(io::IO; own::Bool=false) = NpzReader(ZipFile.Reader(io, own))

NpzReader(fn::AbstractString) = NpzReader(ZipFile.Reader(fn))

Base.keys(r::NpzReader) = keys(r.items)

Base.close(r::NpzReader) = close(r.zipreader)

"""
    readnpy(r::NpzReader, k::AbstractString; transpose=false)

Return the array at key `k` from reader `r`.

Equivalent to `r[k; ...]`.
"""
function readnpy(r::NpzReader, k::AbstractString; kw...)
    f = r.items[k]
    seekstart(f)
    return readnpy(f; kw...)
end

"""
    readnpyheader(r::NpzReader, k::AbstractString)

Return the .npy header at key `k` from reader `r`.
"""
function readnpyheader(r::NpzReader, k::AbstractString; kw...)
    f = r.items[k]
    seekstart(f)
    return readnpyheader(f; kw...)
end

Base.getindex(r::NpzReader, k::AbstractString; kw...) = readnpy(r, k; kw...)

function Base.show(io::IO, ::MIME"text/plain", r::NpzReader)
    ks = sort(collect(keys(r)))
    print(io, "NpzReader with $(length(ks)) entries")
    if !isempty(ks)
        print(io, ":")
        for k in ks
            println(io)
            print(io, " ")
            show(io, k)
            print(io, " => ")
            h = readnpyheader(r, k)
            print(io, Base.dims2string(h.shape), " ")
            T = Array{h.descr.eltype, length(h.shape)}
            show(io, T)
        end
    end
end

"""
    readnpz(io_or_filename, [key_or_keys]; transpose=false)

Read arrays from a .npz file.

Returns a `Dict` of all arrays by default. The second argument can select a subset of keys
or a single key.
"""
function readnpz(src::Union{IO,AbstractString}, ks=nothing; kw...)
    readnpz(src) do r
        if ks isa AbstractString
            return readnpy(r, ks; kw...)
        else
            return Dict(k => readnpy(r, k; kw...) for k in (ks === nothing ? keys(r) : ks))
        end
    end
end

"""
    readnpz(f::Function, io_or_filename)

Open a .npz file and return `f(reader::NpzReader)`.
"""
function readnpz(f::Function, io::IO)
    return f(NpzReader(io))
end
function readnpz(f::Function, fn::AbstractString)
    r = NpzReader(fn)
    try
        return f(r)
    finally
        close(r)
    end
end

"""
    NpzWriter(io::IO; own=false)
    NpzWriter(fn::AbstractString)

Open the given .npz file for writing.

- `keys(writer)` returns the written keys.
- `writenpy(writer, key, array)` or `writer[key] = array` writes the given array.
- `close(writer)` closes the underlying IO stream.
"""
struct NpzWriter
    zipwriter::ZipFile.Writer
    items::Dict{String,ZipFile.WritableFile}
    function NpzWriter(zw::ZipFile.Writer)
        items = Dict{String,ZipFile.WritableFile}()
        new(zw, items)
    end
end

NpzWriter(io::IO; own::Bool=false) = NpzWriter(ZipFile.Writer(io, own))

NpzWriter(fn::AbstractString) = NpzWriter(ZipFile.Writer(fn))

"""
    writenpy(w::NpzWriter, k::AbstractString, x::AbstractArray; compress=false)

Write the array `x` to writer `w` at key `k`.

Equivalent to `w[k] = x`.
"""
function writenpy(w::NpzWriter, k::AbstractString, x::AbstractArray; compress::Bool=false, kw...)
    k = convert(String, k)
    method = compress ? ZipFile.Deflate : ZipFile.Store
    f = ZipFile.addfile(w.zipwriter, "$(k).npy", method=method)
    w.items[k] = f
    writenpy(f, x; kw...)
    close(f)
    return
end

function Base.setindex!(w::NpzWriter, x::AbstractArray, k::AbstractString; kw...)
    writenpy(w, k, x; kw...)
    return w
end

Base.keys(w::NpzWriter) = keys(w.items)

Base.length(w::NpzWriter) = length(keys(w))

Base.close(w::NpzWriter) = close(w.zipwriter)

function Base.show(io::IO, ::MIME"text/plain", w::NpzWriter)
    ks = sort(collect(keys(w)))
    print(io, "NpzWriter with $(length(ks)) items")
    if !isempty(ks)
        print(io, ":")
        for k in ks
            println(io)
            print(io, " ")
            show(io, k)
            print(io, " => ...")
        end
    end
end

"""
    writenpz(io_or_fn, arrays::Dict; compress=false)
    writenpz(io_or_fn, arrays::Pair...; compress=false)

Write the given arrays to a .npz file.

Any iterable of `Pair{<:AbstractString,<:AbstractArray}` can be written.
"""
function writenpz(src, arrays; kw...)
    writenpz(src) do w
        for (k, v) in arrays
            writenpy(w, k, v; kw...)
        end
    end
end
writenpz(src, arrays::Pair...; kw...) = writenpz(src, arrays; kw...)

"""
    writenpz(f::Function, io_or_filename)

Open a .npz file for writing and return `f(writer::NpzWriter)`.
"""
function writenpz(f::Function, io::IO)
    return f(NpzWriter(io))
end

function writenpz(f::Function, fn::AbstractString)
    w = NpzWriter(fn)
    try
        return f(w)
    finally
        close(w)
    end
end


text:
# Benchmark two solvers on a set of small problems and profile.
using Optimize

# Common setup
n = 5000
mpb_probs = filter(name -> name != :OptimizationProblems && name != :sbrybnd, names(OptimizationProblems))
ampl_prob_dir = "/home/local/USHERBROOKE/dusj1701/Documents/Recherche/TR-ARC/UnifiedImplementation/Julia/AMPLUnconstrained/"
ampl_probs = [Symbol(split(p, ".")[1]) for p in filter(x -> contains(x, ".nl"), readdir(ampl_prob_dir))]

# Example 1b: benchmark solvers on a set of problems
function compare_solvers(solvers,probs;title:: String = " ", kwargs...)
  bmark_args = Dict{Symbol, Any}(:skipif => model -> model.meta.ncon > 0, :max_f => 20000)#,:atol=>1e-10,:rtol=>1e-8)
  profile_args = Dict{Symbol, Any}(:title => title)
  stats, profiles = bmark_and_profile(solvers,
                                      #(MathProgNLPModel(eval(p)(n), name=string(p)) for p in mpb_probs),

                                      (MathProgNLPModel(eval(p)(n),
                                                        name=string(p)#,
                                                        #features = [:Grad]
                                                        )
                                       for p in probs),
                                      bmark_args=bmark_args, profile_args=profile_args)
    # Other available features :Jac, :Hess, :HessVec, :ExprGraph

  return stats, profiles
end


text:
"""
    IsColdEncoded(Y)

Returns a boolean true if the array Y is cold encoded, and false if not.
"""
IsColdEncoded(Y) = size( forceMatrix( Y ) )[2] == 1

struct ClassificationLabel
    ToHot::Dict
    ToCold::Dict
    LabelCount::Int
end

""""
    LabelEncoding( HotOrCold )

Determines if an Array, `Y`, is one hot encoded, or cold encoded by it's dimensions.
Returns a ClassificationLabel object/schema to convert between the formats.
"""
function LabelEncoding(HotOrCold)
    HotOrCold = forceMatrix(HotOrCold)
    if IsColdEncoded(HotOrCold)
        Lbls = unique(HotOrCold)
    else
        Lbls = (1:size(HotOrCold)[2])
    end
    return ClassificationLabel(Dict( Lbls .=> collect(1 : length(Lbls) ) ),
                                Dict( collect(1 : length(Lbls)) .=> Lbls ),
                                length(Lbls) )
end

"""
    ColdToHot(Y, Schema::ClassificationLabel)

Turns a cold encoded `Y` vector into a one hot encoded array.
"""
function ColdToHot(Y, Schema::ClassificationLabel)
    lenY = length( Y )
    Output = zeros( lenY, Schema.LabelCount )
    for y in 1 : lenY
        Output[y, Schema.ToHot[ Y[y] ] ] = 1
    end
    return Output
end

"""
    HotToCold(Y, Schema::ClassificationLabel)

Turns a one hot encoded `Y` array into a cold encoded vector.
"""
function HotToCold(Y, Schema::ClassificationLabel)
    Y = forceMatrix(Y)
    (lenY, Feats) = size( Y )
    @assert Feats == Schema.LabelCount
    Output = Array{Any,2}(undef, lenY, 1 )
    for y in 1 : lenY
        Output[ y ] =  Schema.ToCold[ findfirst( x -> x == 1, Y[ y , : ] ) ]
    end
    return Output
end

"""
    StatsFromTFPN(TP, TN, FP, FN)

Calculates many essential classification statistics based on the numbers of True Positive(`TP`), True Negative(`TN`),
 False Positive(`FP`), and False Negative(`FN`) examples.
"""
function StatsFromTFPN(TP, TN, FP, FN)
    ConfusionMatrix = reshape( [TP, FP, FN, TN], 2, 2 )
    Precision = TP / ( TP + FP )
    Recall = TP / ( TP + FN )
    Specificity = TN / ( TN + FP )
    Accuracy = ( TP + TN ) / ( TP + TN + FP + FN )
    FMeasure = 2.0 * ( ( Precision * Recall ) / ( Precision + Recall ) )
    FAR = FP / ( FP + TN )
    FNR = FN / ( FN + TP )
    Prevalence = TP + FN
    return Dict("ConfusionMatrix" => ConfusionMatrix,
                "TP" => TP, "FP" => FP, "TN" => TN, "FN" => FN,
                "Specificity" => Specificity,
                "Precision" => Precision,       "Recall" => Recall,
                "Accuracy" => Accuracy,         "FMeasure" => FMeasure,
                "FAR" => FAR,                   "FNR" => FNR,
                "Prevalence" => Prevalence )
end

"""
    MulticlassStats(Y, GT, schema; Microaverage = true)

Calculates many essential classification statistics based on predicted values `Y`, and ground truth values `GT`, using
the encoding `schema`. Returns a tuple whose first entry is a dictionary of averaged statistics, and whose second
entry is a dictionary of the form "Class" => Statistics Dictionary ...
"""
function MulticlassStats(Y, GT, schema; Microaverage = true)
    Y = forceMatrix(Y)
    GT = forceMatrix(GT)
    Y = IsColdEncoded(Y) ?  map(x -> Int(schema.ToHot[ x ]), Y) : HighestVote(Y)
    GT = IsColdEncoded(GT) ? map(x -> Int(schema.ToHot[ x ]), GT) : HighestVote(GT)
    @assert(size(Y)[1] == size(GT)[1])
    ConfusionMatrix = zeros( schema.LabelCount, schema.LabelCount )
    for y in 1 : size(Y)[1]
        ConfusionMatrix[ Y[ y ], GT[ y ] ] += 1
    end
    TP = zeros( schema.LabelCount ) ; TN = zeros( schema.LabelCount )
    FP = zeros( schema.LabelCount ) ; FN = zeros( schema.LabelCount )
    for c in 1 : schema.LabelCount
        TP[c] = ConfusionMatrix[c,c]
        FP[c] = sum(ConfusionMatrix[:,c]) - TP[c]
        FN[c] = sum(ConfusionMatrix[c,:]) - TP[c]
        TN[c] = sum(ConfusionMatrix) - TP[c] - FP[c] - FN[c]
    end

    #Compute Class-wise statistics.
    #Make a dictionary which contains for each class, a dictionary of statistics.
    ClasswiseStats = Dict()
    for class in 1 : schema.LabelCount
        ClasswiseStats[schema.ToCold[class]] = StatsFromTFPN(TP[class], TN[class], FP[class], FN[class])
    end
    GlobalStats = Dict()
    if Microaverage
         TP = StatsBase.mean( TP ); TN = StatsBase.mean( TN )
         FP = StatsBase.mean( FP ); FN = StatsBase.mean( FN )
         GlobalStats = Dict(:Global => StatsFromTFPN( TP, TN, FP, FN ) )
    else #Macro Average
        Precision = StatsBase.mean( TP ./ ( TP .+ FP ) )
        Recall = StatsBase.mean( TP ./ ( TP .+ FN ) )
        Specificity = StatsBase.mean( TN ./ ( TN .+ FP ) )
        Accuracy = StatsBase.mean( ( TP .+ TN ) ./ ( TP .+ TN .+ FP .+ FN ) )
        FMeasure = StatsBase.mean( 2.0 .* ( ( Precision .* Recall ) ./ ( Precision .+ Recall ) ) )
        FAR = StatsBase.mean( FP ./ ( FP .+ TN ) )
        FNR = StatsBase.mean( FN ./ ( FN .+ TP ) )
        GlobalStats = Dict(:Global => Dict( "ConfusionMatrix" => ConfusionMatrix,
                            "TP" => TP, "FP" => FP, "TN" => TN, "FN" => FN,
                            "Specificity" => Specificity,
                            "Precision" => Precision,       "Recall" => Recall,
                            "Accuracy" => Accuracy,         "FMeasure" => FMeasure,
                            "FAR" => FAR,                   "FNR" => FNR ) )
    end
    return (GlobalStats, ClasswiseStats)
end

"""
    StatsDictToDataFrame(DictOfStats, schema)

Converts a dictionary of statistics which is returned from `MulticlassStats` into a labelled dataframe.
This is an intermediate step for automated report generation.
"""
function StatsDictToDataFrame(ClasswiseStats; digits = 4,
                                StatsList = [   "FMeasure", "Accuracy", "Specificity",
                                                "Precision", "Recall", "FAR", "FNR" ])
    ClassName = repeat( [ "No Name" ], length(ClasswiseStats) );
    daf = DataFrame( :Statistics => StatsList)
    for (j, ( classnames, stats )) in enumerate(ClasswiseStats)
        if typeof(classnames) == Symbol
            classnames = string(classnames)
        end
        ClassName[ j ] = classnames
        ClassStats = zeros( length(StatsList) );
        for (i, stat) in enumerate( StatsList )
            ClassStats[ i ] = round( stats[ stat ]; sigdigits = digits )
        end
        daf[ Symbol( classnames ) ] = ClassStats
    end
    return daf
end

"""
    StatsToDataFrame(stats, schema, filepath, name)

Converts the 2-Tuple returned from `MulticlassStats()` (`stats`) to a CSV file with a specified `name`
in a specified `filepath` using the prescribed encoding `schema`.

The statistics associated with the global analysis will end in a file name  of "-global.csv"
and the local statistics for each class will end in a file named "-classwise.csv"
"""
function StatsToCSVs(Stats, filepath, name)
    globaldf = StatsDictToDataFrame(Stats[1])
    localdf = StatsDictToDataFrame(Stats[2])
    CSV.write(Base.joinpath(filepath, name * "-global.csv"), globaldf)
    CSV.write(Base.joinpath(filepath, name * "-classwise.csv"), localdf)
end

"""
    DataFrameToLaTeX( df, caption = "" )

Converts a DataFrame object to a LaTeX table (string).
"""
function DataFrameToLaTeX( df; caption = "")
    (Rows, Columns) = size(df)
    ColmFormat = reduce(*, ["c" for i in 1:Columns])
    ColmNames = join(string.(names(df))," & ")
    retstr = "\\begin{table}[h] \n" *
             "\t\\begin{tabular}{$ColmFormat} \n" *
             "\t\t $ColmNames \\\\ \\hline \n"

    TableInterior = [ "\t\t " * join(string.(values(df[row,:])), " & ") * " \\\\ \n" for row in 1:Rows]
    TableInterior = reduce(*, TableInterior)

    retstr *= TableInterior *
              "\t \\end{tabular} \n"
    retstr *= (length(caption) > 0) ? "\t \\caption{$caption} \n" : ""
    retstr *= "\\end{table}\n"
    return retstr
end

"""
    StatsToLaTeX(Stats, filepath = nothing, name = nothing,
                        digits = 3, maxcolumns = 6; Comment = "",
                        StatsList = [   "FMeasure", "Accuracy", "Specificity",
                                        "Precision", "Recall", "FAR", "FNR" ])

Converts a MulticlassStats object to a LaTeX table (string or saved file).
LaTeX tables contain rows of StatsList, and a maximum column number of maxcolumns.
Information is presented with a set number of decimals(digits).

"""
function StatsToLaTeX(Stats, filepath = nothing, name = nothing;
                        digits = 3, maxcolumns = 6, Comment = "",
                        StatsList = [   "FMeasure", "Accuracy", "Specificity",
                                        "Precision", "Recall", "FAR", "FNR"     ] )
    globaldf = StatsDictToDataFrame(Stats[1]; digits = digits, StatsList = StatsList)
    classes = [ k for k in keys( Stats[2] ) ]
    localdf = []
    if length(classes) > maxcolumns
        maxtbls = Int( round( length( classes ) / maxcolumns ) )
        tmpdf = StatsDictToDataFrame(Stats[2]; digits = digits, StatsList = StatsList)
        for rn in 1 : maxtbls
            colview = []
            if rn < maxtbls
                colview = classes[ ((maxcolumns * (rn - 1)) + 1) : (rn * maxcolumns)]
            else
                colview = classes[ ((maxcolumns * (rn - 1)) + 1) : end]
            end
            push!(localdf, tmpdf[:, Symbol.( vcat(["Statistics"], colview))] )
        end
    else
        localdf = StatsDictToDataFrame(Stats[2]; digits = digits, StatsList = StatsList)
    end
    TimeStamp = Dates.format(now(), "mm-dd-YYYY HH:MM")
    ReportStr = "\\documentclass[]{report}\n" *
                "% Report Generated from ChemometricsTools.jl ($TimeStamp)\n" *
                "% $Comment" *
                "\n\\begin{document}\n"
    ReportStr *= DataFrameToLaTeX( globaldf; caption = "Global Classification Statistics." )
    if length(classes) > maxcolumns
        for tdf in localdf
            ReportStr *= DataFrameToLaTeX( tdf; caption = "Classwise Classification Statistics." )
        end
    else
        ReportStr *= DataFrameToLaTeX( localdf; caption = "Classwise Classification Statistics." )
    end
    ReportStr *= "\n\\end{document}"
    #If no name is given -> Return the string
    if isa(filepath, Nothing) or isa(name, Nothing)
        return ReportStr
    else
        open( filepath * name * ".tex", "w" ) do f
            write( f, ReportStr )
        end
        return true
    end
end

"""
    Threshold(yhat; level = 0.5)

For a binary vector `yhat` this decides if the label is a 0 or a 1 based on it's value relative to a threshold `level`.
"""
Threshold(yhat; level = 0.5) = map( y -> (y >= level) ? 1 : 0, yhat)

"""
    MulticlassThreshold(yhat; level = 0.5)

Effectively does the same thing as Threshold() but per-row across columns.

*Warning this function can allow for no class assignments. HighestVote is preferred*
"""
function MulticlassThreshold(yhat; level = 0.5)
    newY = zeros(size(yhat))
    for obs in 1 : size(yhat)[1]
        (val, ind) = findmax( yhat[obs,:] )
        if val > level
            newY[ind] = val
        end
    end
    return newY
end

"""
    HighestVote(yhat)

Returns the column index for each row that has the highest value in one hot encoded `yhat`. Returns a one cold encoded vector.
"""
function HighestVote( yhat )
    return [ findmax( yhat[obs,:] )[2] for obs in 1 : size(yhat)[1]  ]
end

"""
    HighestVoteOneHot(yhat)

Turns the highest column-wise value to a 1 and the others to zeros per row in a one hot encoded `yhat`. Returns a one cold encoded vector.
"""
function HighestVoteOneHot( yhat )
    (Obs, Classes) = size( yhat )
    ret = zeros( ( Obs, Classes ) )
    for o in 1 : Obs
        ret[ o, argmax( yhat[ o, : ] ) ] = 1
    end
    return ret
end


text: struct S a; b; c; end

max_stars_repo_name: gabriel-doriath-dohler/pjulia

max_issues_repo_name: gabriel-doriath-dohler/pjulia

max_forks_repo_name: gabriel-doriath-dohler/pjulia


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

"""
    MetricBall(radii, [angles]; convention=TaitBryanExtr)
    MetricBall(radius, metric=Euclidean())

A metric ball is a neighborhood that can be expressed in terms
of a metric and a set of `radii`. The two main examples are the
Euclidean ball an the Mahalanobis (ellipsoid) ball.

When multiple `radii` are provided, they can be rotated by `angles`
according to a given rotation `convention`. Alternatively, a metric
from the [Distances.jl](https://github.com/JuliaStats/Distances.jl)
package can be specified together with a single `radius`.

## Examples

N-dimensional Euclidean ball with radius `1.0`:

```julia
julia> euclidean = MetricBall(1.0)
```

Axis-aligned 3D ellispoid with radii `[3.0, 2.0, 1.0]`:

```julia
julia> mahalanobis = MetricBall([3.0, 2.0, 1.0])
```

See also [`mahalanobis`](@ref).
"""
struct MetricBall{R,M} <: Neighborhood
  radii::R
  metric::M
end

function MetricBall(radii::AbstractVector{T}, angles=nothing;
                    convention=TaitBryanExtr) where {T<:Real}
  ndims  = length(radii)
  angles = isnothing(angles) ? (ndims == 2 ? [zero(T)] : zeros(T, 3)) : angles
  metric = mahalanobis(radii, angles, convention)
  MetricBall{typeof(radii),typeof(metric)}(radii, metric)
end

function MetricBall(radius::T, metric=Euclidean()) where {T<:Real}
  MetricBall{typeof(radius),typeof(metric)}(radius, metric)
end

"""
    radii(ball)

Return the radii of the metric `ball`.
"""
radii(ball::MetricBall) = ball.radii

"""
    metric(ball)

Return the metric of the metric `ball`.
"""
metric(ball::MetricBall) = ball.metric

"""
    range(ball)

Return the range of the metric `ball`, i.e.
the smallest value `r` such that `||v|| ≤ r`
for any `v ∈ ball`.
"""
Base.range(ball::MetricBall{<:Real}) = ball.radii
Base.range(::MetricBall{R,<:Mahalanobis}) where {R} = one(eltype(R))

function Base.show(io::IO, ball::MetricBall)
  r = length(ball.radii) > 1 ? Tuple(ball.radii) : ball.radii
  m = nameof(typeof(ball.metric))
  print(io, "MetricBall($r, $m)")
end


text:
# Example of robot kinematics constraints reported by John Gustafson
# http://www.johngustafson.net/presentations/Multicore2016-JLG.pdf

using ConstraintPropagation

d = Domain()

@add_constraint d s2*c5*s6 - s3*c5*s6 - s4*c5*s6 + c2*c6 + c3*c6 + c4*c6 == 0.4077
@add_constraint d c1*c2*s5 + c1*c3*s5 + c1*c4*s5 + s1*c5 == 1.9115
@add_constraint d s2*s5 + s3*s5 + s4*s5 == 1.9791
@add_constraint d c1*c2 + c1*c3 + c1*c4 + c1*c2 + c1*c3 + c1*c2 == 4.0616
@add_constraint d s1*c2 + s1*c3 + s1*c4 + s1*c2 + s1*c3 + s1*c3 == 1.7172
@add_constraint d s2 + s3 + s4 + s2 + s3 + s2 == 3.9701
@add_constraint d s1^2 + c1^2 == 1
@add_constraint d s2^2 + c2^2 == 1
@add_constraint d s3^2 + c3^2 == 1
@add_constraint d s4^2 + c4^2 == 1
@add_constraint d s5^2 + c5^2 == 1
@add_constraint d s6^2 + c6^2 == 1


text:
"""
    animate(vis::MechanismVisualizer,
            times::Vector{Float64},
            configurations::Vector{Vector{Float64}};
            fps::Float64=60, realtimerate::Float64=1.)

Animate the given mechanism passing through a time-coded series of
configurations by linearly interpolating the configuration vectors.
"""
function animate(vis::MechanismVisualizer,
                 times::Vector{Float64},
                 configurations::AbstractVector{<:AbstractVector{Float64}};
                 fps::Float64 = 60., realtimerate::Float64 = 1.)
    @assert fps > 0
    @assert 0 < realtimerate < Inf

    state = vis.state
    interpolated_configurations = interpolate((times,), configurations, Gridded(Linear()))
    t0, tf = first(times), last(times)
    framenum = 0
    walltime0 = time()
    @throttle framenum while true
        t = min(tf, t0 + (time() - walltime0) * realtimerate)
        q = interpolated_configurations(t)
        set_configuration!(state, q)
        rbd.normalize_configuration!(state)
        _render_state!(vis)
        framenum += 1
        t == tf && break
    end max_rate = fps
end

function MeshCat.Animation(mvis::MechanismVisualizer,
        times::AbstractVector{<:Real},
        configurations::AbstractVector{<:AbstractVector{<:Real}};
        fps::Integer=30)
    @assert axes(times) == axes(configurations)
    interpolated_configurations = interpolate((times,), configurations, Gridded(Linear()))
    animation = Animation()
    num_frames = floor(Int, (times[end] - first(times)) * fps)
    for frame in 0 : num_frames
        time = first(times) + frame / fps
        let mvis = mvis, interpolated_configurations = interpolated_configurations, time=time
            atframe(animation,  frame) do
                set_configuration!(mvis, interpolated_configurations(time))
            end
        end
    end
    return animation
end

MeshCat.setanimation!(mvis::MechanismVisualizer, args...; kw...) =
    setanimation!(visualizer(mvis), args...; kw...)

function MeshCat.setanimation!(mvis::MechanismVisualizer,
                      times::AbstractVector{<:Real},
                      configurations::AbstractVector{<:AbstractVector{<:Real}};
                      fps::Integer=30,
                      play::Bool=true,
                      repetitions::Integer=1)
    Base.depwarn("""
    `setanimation!(mvis, times, configurations; ..)` is deprecated. Instead, you can construct an `Animation` and then call `setanimation!` with the result.

    For example, if you previously did:

    ```
    setanimation!(mvis, times, configurations, fps=30, play=true)
    ```

    You should now do:

    ```
    animation = Animation(mvis, times, configurations, fps=25)
    setanimation!(mvis, animation, play=true)
    ```
    """, :setanimation_t_q)
    animation = Animation(mvis, times, configurations, fps=fps)
    setanimation!(mvis, animation, play=play, repetitions=repetitions)
end


text:
###### readtimearray ############

function readtimearray(fname::String)
    blob    = readcsv(fname)
    tstamps = Date[Date(i) for i in blob[2:end, 1]]
    vals    = insertNaN(blob[2:end, 2:end])
    cnames  = UTF8String[]
    for b in blob[1, 2:end]
        push!(cnames, b)
    end
    TimeArray(tstamps, vals, cnames)
end

function insertNaN{N}(aa::Array{Any,N})
    for i in 1:size(aa,1)
        for j in 1:size(aa,2)
            if !isa(aa[i,j], Real)
                aa[i,j] = NaN
            end
        end
    end
    convert(Array{Float64,N},aa)
end


text:
module BinaryProvider

using Libdl, Logging
using Pkg, Pkg.PlatformEngines, Pkg.BinaryPlatforms
import Pkg.PlatformEngines: package, download
import Pkg.BinaryPlatforms: Linux
export platform_key, platform_key_abi, platform_dlext, valid_dl_path,
       triplet, select_platform, platforms_match,
       Linux, MacOS, Windows, FreeBSD,
       parse_7z_list, parse_tar_list, verify,
       download_verify, unpack, package, download_verify_unpack,
       list_tarball_files, list_tarball_symlinks

# Some compatibility mapping
const choose_download = Pkg.BinaryPlatforms.select_platform
Linux(arch::Symbol, libc::Symbol) = Linux(arch; libc=libc)
function platform_key(machine::AbstractString = Sys.MACHINE)
    Base.depwarn("platform_key() is deprecated, use platform_key_abi() from now on", :binaryprovider_platform_key)
    platkey = platform_key_abi(machine)
    return typeof(platkey)(arch(platkey); libc=libc(platkey), call_abi=call_abi(platkey))
end

# Utilities for controlling verbosity
include("LoggingUtils.jl")
# Include our subprocess running functionality
include("OutputCollector.jl")
# Everything related to file/path management
include("Prefix.jl")
# Abstraction of "needing" a file, that would trigger an install
include("Products.jl")

function __init__()
    global global_prefix

    # Initialize our global_prefix
    global_prefix = Prefix(joinpath(@__DIR__, "..", "global_prefix"))

    # Find the right download/compression engines for this platform
    probe_platform_engines!()
end

# gotta go fast
include("precompile.jl")
_precompile_()

end # module


text:
# zt: should Linear Omega be default?
"""
LinearΩ: Stores data in single `AbstractVector{V}`

`lω[i] = lω.ωvec[lω.ids[i]]`

- No overhead for linearization / delinearization
- Good for optimization / inference methods which require linear interface
- If heterogenous elements are used, must convert (which can be slow)
"""
struct LinearΩ{I, SEG, V <: AbstractArray} <: ΩBase{I}
  ids::Dict{I, SEG} # zt: too precise, any associable will do. Map from Ω index I to segment of omvec
  ωvec::V
end

"`emp(::Type{T})` Empty container of type T"
function emp end

emp(a::Type{<:Array}) = a()

"$(SIGNATURES) Empty `LinearΩ`"
LinearΩ{I, SEG, V}() where {I, SEG, V} = LinearΩ{I, SEG, V}(Dict{I, SEG}(), emp(V))

linearize(lω::LinearΩ) = lω.ωvec
unlinearize(ωvec, lω::LinearΩ) = LinearΩ(lω.ids, ωvec)

# zt: fixme, dont use T, determine it by the input
flat(rv, ω::T) where T <: LinearΩ = floatvec -> rv(T(ω.ids, floatvec))

# flat2(rv, ω::LinearΩ) = LinearΩ(floatvec)

# zt: renamed this? Wanted to distinguish getting ith Omega value from
# projection. Maybe projection shouldnt overload getindex?
getdim(lω, i) = lω.ωvec[i]

# Update #

"Apply `kernel` to ith component"
function update(lω::LinearΩ, i::Int, kernel::Function) # zt: use name other than kernel
  # zt: `kernel` Might be a callable
  # zt: isn't there a more efficient version which just copies ωvec?
  lω_ = deepcopy(lω)
  lω_.ωvec[i] = kernel(lω_.ωvec[i])
  lω_
end

function update(lω::LinearΩ, i::Int, val)
  # This update != previous update, dont pun!
  lω_ = deepcopy(lω)
  lω_.ωvec[i] = val
  lω_
end

# randtype #

# Defualt behaviour: return type is T
randrtype(lω::LinearΩ{I, SEG, V}, ::Type{T}) where {T, I, SEG, V} = T

# Handle ForwardDiff: If `lω.ωvec` is a Vector of Duals then return type will be a DUal
randrtype(lω::LinearΩ{I, SEG, DT}, ::Type{T}) where {Q<:ForwardDiff.Dual, T, I, SEG, DT <: AbstractArray{Q}} = Q

# If Dims specified, return array of this type
randrtype(ω, T, ::Dims{N}) where N = Array{randrtype(ω, T), N}

# Specialization for when T is Type
randrtype(ω, ::Type{T}, ::Dims{N}) where {T, N} = Array{randrtype(ω, T), N}

# WHAT ABOUT Static Array
# sol? If V is a static vector then return a StaticArray
#    prob: reshaping it will return a reshape vector
#    prob: what if we want an actual staticarray instead
#  How to distingusih between staticarray rt and reshapedarray
# 1. specify in the model code
# 2. one issue is that Dims values aren't included in type information\: so there's no enough information in the model to specify the shape
# 3. It's might be possible with constant propagation that you get static types
# 4.
# WHAT ABOUT Float32
# I think that's handled fine from the type of T
# Is there ever a case when say we want values to be tracked (duals) and some not?
# Maybe if we have some discrete parameters and dont want them to continuized
# But what if we do want them to be.
#

nelem(lω) = length(lω.ωvec)
Base.values(lω::LinearΩ) = [lω.ωvec]
Base.isempty(lω::LinearΩ) = isempty(lω.ωvec)

# memrand #

# Flux-specific #

emp(a::Type{Tracker.TrackedArray{A, B, C}}) where {A, B, C} = C()

# zt: This is type piracy: use a different method name
Base.append!(ta::Tracker.TrackedArray, a::Array) =
  (append!(ta.data, a); append!(ta.grad, zero(a)))

Base.append!(ta::Tracker.TrackedArray, tb::Tracker.TrackedArray) =
  (append!(ta.data, tb.data); append!(ta.grad, tb.grad))

# using ZenUtils

function randrtype(lω::LinearΩ{I, SEG, V}, ::Type{T}, dims::Dims{N}) where {I, SEG, T, V <: Tracker.TrackedArray{T}, N}
  Tracker.TrackedArray{T, N, Array{T, N}}
end

# function randrtype(lω::LinearΩ{I, SEG, V}, ::Type{T}) where {I, SEG, V<: Tracker.TrackedArray, T <: AbstractFloat}
#   Tracker.TrackedReal{T}
# end

# coerce(::Type{T}, x) where T = convert(T, x)
# coerse(::Type{T}, x) where T

function memrand(lω::LinearΩ, id, X, dims::Dims; rng)
  RT = randrtype(lω,  X,  dims)
  seg = get(lω.ids, id, 0)
  if seg == 0
    res = rand(rng, X, dims)
    startidx = length(lω.ωvec) + 1
    lω.ids[id] = startidx:startidx + prod(dims) - 1
    res_::RT = res  # zt: issue si that might need to convert into correct format
    typeof(res_)
    typeof(lω.ωvec)
    # @grab res_
    p = lω.ωvec
    # @grab p
    append!(lω.ωvec, res_)
    res_
  else
    seg = lω.ids[id]
    subωvec = lω.ωvec[seg]
    res::RT = reshape(subωvec, dims) # reshaping a thing might produce a vector
    res
  end
end

memrand(lω::LinearΩ{I, SEG, V}, id, X; rng) where {I, SEG, V<: Tracker.TrackedArray} =
  first(memrand(lω, id, X, (1,); rng = rng))

randrtype(lω, x::Distributions.Distribution) = eltype(x)

function memrand(lω::LinearΩ, id, X; rng)
  RT = randrtype(lω, X)
  seg = get(lω.ids, id, 0)
  if seg == 0
    res::RT = rand(rng, X)
    res_::RT = res
    push!(lω.ωvec, res_)
    startidx = length(lω.ωvec)
    lω.ids[id] = startidx:startidx
    res_
  else
    # @show seg = lω.ids[id]
    # @show id
    subωvec = lω.ωvec[first(seg)]::RT
  end
end


text:
import SPTK: mfcc

function test_mfcc(len, order)
    srand(98765)
    dummy = rand(len)

    # By default, c0 is not contained
    cc = mfcc(dummy, order, czero=true, power=true)
    # @test length(cc) == order
    if any(isnan.(cc))
        @show cc
    end
    @test all(isfinite.(cc))
end

function test_mfcc_options()
    srand(98765)
    dummy = rand(512)

    try mfcc(dummy, 19, numfilterbunks=20) catch @test false end
    @test_throws ArgumentError mfcc(dummy, 20, numfilterbunks=20)

    # with c0
    println("-- test c0 is expected")
    cc = mfcc(dummy, 12; czero=true)
    @test length(cc) == 13

    # with c0 + power
    println("-- test c0 and power are expected")
    cc = mfcc(dummy, 12; czero=true, power=true)
    @test length(cc) == 14

    # with power
    println("-- test power is expected")
    cc = mfcc(dummy, 12; czero=false, power=true)
    @test length(cc) == 13
end

# TODO: Update binary dependency for windows
if !is_windows()
println("-- test_mfcc")
for len in [256, 512, 1024, 2048]
    for order in [12, 14, 16, 18]
        println(" where len = $len, order = $order")
        test_mfcc(len, order)
    end
end

test_mfcc_options()
end


text:
module JacobianTest

import Calculus

using Base.Test
using ForwardDiff

include(joinpath(dirname(@__FILE__), "utils.jl"))

##################
# hardcoded test #
##################

f! = (y, x) -> begin
    y[1] = x[1] * x[2]
    y[1] *= sin(x[3]^2)
    y[2] = y[1] + x[3]
    y[3] = y[1] / y[2]
    y[4] = x[3]
    return nothing
end
f = x -> (y = zeros(promote_type(eltype(x), Float64), 4); f!(y, x); return y)
y = zeros(4)
x = [1, 2, 3]
v = f(x)
j = [0.8242369704835132  0.4121184852417566  -10.933563142616123
     0.8242369704835132  0.4121184852417566  -9.933563142616123
     0.169076696546684   0.084538348273342   -2.299173530851733
     0.0                 0.0                 1.0]

for c in (Chunk{1}(), Chunk{2}(), Chunk{3}())
    # testing f(x)
    @test_approx_eq j ForwardDiff.jacobian(f, x, c)

    out = zeros(4, 3)
    ForwardDiff.jacobian!(out, f, x, c)
    @test_approx_eq out j

    # testing f!(y, x)
    @test_approx_eq j ForwardDiff.jacobian(f!, y, x, c)
    @test_approx_eq v y

    fill!(y, 0.0)
    out = zeros(4, 3)
    ForwardDiff.jacobian!(out, f!, y, x, c)
    @test_approx_eq out j
    @test_approx_eq y v

    out = JacobianResult(zeros(4), zeros(4, 3))
    ForwardDiff.jacobian!(out, f!, x, c)
    @test_approx_eq ForwardDiff.value(out) v
    @test_approx_eq ForwardDiff.jacobian(out) j
end

########################
# test vs. Calculus.jl #
########################

for (f!, f) in VECTOR_TO_VECTOR_FUNCS
    v = f(X)
    j = ForwardDiff.jacobian(f, X)
    @test_approx_eq_eps j Calculus.jacobian(f, X, :forward) FINITEDIFF_ERROR
    for c in CHUNK_SIZES
        for usecache in (true, false)
            chunk = Chunk{c}()

            # testing f(x)
            println("  ...testing $f with (chunk size = $c) and (usecache = $usecache)")
            out = ForwardDiff.jacobian(f, X, chunk; usecache = usecache)
            @test_approx_eq out j

            out = similar(Y, length(Y), length(X))
            ForwardDiff.jacobian!(out, f, X, chunk; usecache = usecache)
            @test_approx_eq out j

            # testing f!(y, x)
            println("  ...testing $(f!) with (chunk size = $c) and (usecache = $usecache)")
            y = zeros(Y)
            out = ForwardDiff.jacobian(f!, y, X, chunk; usecache = usecache)
            @test_approx_eq y v
            @test_approx_eq out j

            y = zeros(Y)
            out = similar(Y, length(Y), length(X))
            ForwardDiff.jacobian!(out, f!, y, X, chunk; usecache = usecache)
            @test_approx_eq y v
            @test_approx_eq out j

            out = JacobianResult(zeros(Y), similar(Y, length(Y), length(X)))
            ForwardDiff.jacobian!(out, f!, X, chunk; usecache = usecache)
            @test_approx_eq ForwardDiff.value(out) v
            @test_approx_eq ForwardDiff.jacobian(out) j
        end
    end
end

end # module


text:
@testset "GitLab.MergeRequest equality" begin
    a = GitLab.MergeRequest(; id=1)
    b = GitLab.MergeRequest(; id=1)

    @test a == b
    @test isequal(a, b)
end

@testset "get_pull_requests" begin
    @testset "GitHub" begin
        api = GitForge.GitHub.GitHubAPI()
        repo = GitHub.Repo(; name="Foo", owner=GitHub.User(; login="Bar"))

        apply(gh_unique_patch) do
            prs = CompatHelper.get_pull_requests(
                api, repo, "open"; per_page=10, page_limit=1
            )
            @test !isempty(prs)
        end
    end

    @testset "GitLab" begin
        api = GitForge.GitLab.GitLabAPI()
        repo = GitLab.Project(; id=1)

        apply(gl_unique_patch) do
            prs = CompatHelper.get_pull_requests(api, repo, "opened")
            @test !isempty(prs)
        end
    end
end

@testset "get_pr_titles" begin
    username = "foo"

    @testset "no forks" begin
        @testset "GitHub" begin
            origin_repo = GitHub.Repo(; id=1, fork=false)
            fork_repo = GitHub.Repo(; id=2, fork=true)

            pr_from_origin = GitHub.PullRequest(;
                head=GitHub.Head(; repo=origin_repo),
                user=GitHub.User(; login=username),
                title="PR A",
            )
            pr_from_fork = GitHub.PullRequest(;
                head=GitHub.Head(; repo=fork_repo),
                user=GitHub.User(; login=username),
                title="PR B",
            )

            apply(get_prs_patch([pr_from_origin, pr_from_fork])) do
                result = CompatHelper.get_pr_titles(
                    GitHub.GitHubAPI(), origin_repo, username
                )

                @test length(result) == 1
                @test pr_from_origin.title in result
                @test !(pr_from_fork.title in result)
            end
        end

        @testset "GitLab" begin
            origin_repo = GitLab.Project(; id=1, fork=false)
            fork_repo = GitLab.Project(; id=2, fork=true)

            pr_from_origin = GitLab.MergeRequest(;
                id=1, project_id=1, author=GitLab.User(; username=username), title="PR A"
            )
            pr_from_fork = GitLab.MergeRequest(;
                id=1, project_id=2, author=GitLab.User(; username=username), title="PR B"
            )

            apply(get_prs_patch([pr_from_origin, pr_from_fork])) do
                result = CompatHelper.get_pr_titles(
                    GitLab.GitLabAPI(), origin_repo, username
                )

                @test length(result) == 1
                @test pr_from_origin.title in result
                @test !(pr_from_fork.title in result)
            end
        end
    end

    @testset "only my prs" begin
        @testset "GitHub" begin
            origin_repo = GitHub.Repo(; id=1, fork=false)

            pr_from_me = GitHub.PullRequest(;
                head=GitHub.Head(; repo=origin_repo),
                user=GitHub.User(; login=username),
                title="PR A",
            )
            pr_from_other = GitHub.PullRequest(;
                head=GitHub.Head(; repo=origin_repo),
                user=GitHub.User(; login="bizbaz"),
                title="PR B",
            )

            apply(get_prs_patch([pr_from_me, pr_from_other])) do
                result = CompatHelper.get_pr_titles(
                    GitHub.GitHubAPI(), origin_repo, username
                )

                @test length(result) == 1
                @test pr_from_me.title in result
                @test !(pr_from_other.title in result)
            end
        end

        @testset "GitLab" begin
            origin_repo = GitLab.Project(; id=1, fork=false)

            pr_from_me = GitLab.MergeRequest(;
                id=1, project_id=1, author=GitLab.User(; username=username), title="PR A"
            )
            pr_from_other = GitLab.MergeRequest(;
                id=2, project_id=1, author=GitLab.User(; username="bizbaz"), title="PR B"
            )

            apply(get_prs_patch([pr_from_me, pr_from_other])) do
                result = CompatHelper.get_pr_titles(
                    GitLab.GitLabAPI(), origin_repo, username
                )

                @test length(result) == 1
                @test pr_from_me.title in result
                @test !(pr_from_other.title in result)
            end
        end
    end

    @testset "no forks and only my prs" begin
        @testset "GitHub" begin
            api = GitForge.GitHub.GitHubAPI()
            repo = GitHub.Repo(; id=1, fork=false)

            apply(gh_gpr_patch) do
                prs = CompatHelper.get_pr_titles(api, repo, "foobar")
                @test !isempty(prs)
                @test prs[1] == "title"
            end
        end

        @testset "GitLab" begin
            api = GitForge.GitLab.GitLabAPI()
            repo = GitLab.Project(; id=1)

            apply(gl_gpr_patch) do
                prs = CompatHelper.get_pr_titles(api, repo, "foobar")
                @test !isempty(prs)
                @test prs[1] == "title"
            end
        end
    end
end


text:
using PrintTools
using Test

@testset "PrintTools.jl" begin

    tfile = tempname()
    rm(tfile; force = true)
    token = string(rand(UInt))

    try
        # print functions
        for printfun in [
                print_inmw, println_inmw, print_ifmw, println_ifmw,
                tagprint_inmw, tagprintln_inmw, tagprint_ifmw, tagprintln_ifmw
            ]
            open(tfile, "w") do io
                printfun(io, token)
            end
            str = read(tfile, String)
            @test contains(str, token)
            rm(tfile; force = true)
        end

        # print err
        open(tfile, "w") do io
            try; error("ERROR: ", token)
                catch err; printerr(io, err)
            end
        end
        str = read(tfile, String)
        @test contains(str, token)
        rm(tfile; force = true)

        # sci
        println("\nsci_str")
        for f in [1.0, rand() * 1e3, rand() * 1e-12]
            println(sci_str(f; d = 1))
            println(sci_str(f; d = 3))
            println(sci_str(f; d = 5))
            @test true
        end
        println()

    finally
        rm(tfile; force = true)
    end
end;


text:
# ----------------------------------------------------------------------------------- #
# Copyright (c) 2020 Varnerlab
# Robert Frederick School of Chemical and Biomolecular Engineering
# Cornell University, Ithaca NY 14850

# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
#
# The above copyright notice and this permission notice shall be included in
# all copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
# THE SOFTWARE.
# ----------------------------------------------------------------------------------- #
function solve_static_fba_problem(problem::VLStaticFBAProblem, vargs...)::VLStaticSolution

    # TODO: check that the problem object is complete ...
    # TODO: error?

    # get stuff to solve the static fba problem -
    stoichiometric_matrix = problem.stoichiometric_matrix
    flux_bounds_array = problem.flux_bounds_array
    species_bounds_array = problem.species_bounds_array
    obj_coeff_vector = problem.obj_coeff_vector
    is_minimum_flag = problem.is_minimum_flag

    # solve the static fba problem -
    (objective_value, calculated_flux_array, dual_value_array, uptake_array, exit_flag, status_flag) = calculate_optimal_flux_distribution(stoichiometric_matrix,flux_bounds_array,
        species_bounds_array, obj_coeff_vector; min_flag=is_minimum_flag)

    # TODO: Should we throw an error if the calculation failed? -or- leave it up to the user to check for such stuff?
    # TODO: error -or- warning?

    # create a solution object and return to the user -
    return VLStaticSolution(objective_value, calculated_flux_array, exit_flag, status_flag)
end

function solve_dynamic_fba_problem(problem::VLDynamicFBAProblem, vargs...)::VLDynamicSolutionStep

    # TODO: solve dynamic flux balance analysis problem -
    # TODO: fill me in ...
    return VLDynamicSolutionStep()
end

function solve_simulation_problem(problem::VLAbstractProblem, vargs...)::VLAbstractSolution

    # switch on the type of problem -
    if typeof(problem) == VLStaticFBAProblem
        return solve_static_fba_problem(problem, vargs...)
    elseif (typeof(problem) == VLDynamicFBAProblem)
        return solve_dynamic_fba_problem(problem, vargs...)
    else
        throw(ArgumentError("Ooops! $(typeof(problem)) is not currently supported."))
    end
end

max_stars_repo_name: varnerlab/CFMG

max_issues_repo_name: varnerlab/CFMG

max_forks_repo_name: varnerlab/CFMG