text:
module RecurrenceAnalysis

using Distances, Statistics, LinearAlgebra, SparseArrays, DelayEmbeddings, StaticArrays

const DEFAULT_METRIC = Euclidean()

export RecurrenceMatrix, CrossRecurrenceMatrix, JointRecurrenceMatrix,
       AbstractRecurrenceMatrix, WithinRange, NeighborNumber, FAN

export embed,
       reconstruct,
       Dataset,
       distancematrix,
       textrecurrenceplot,
       recurrenceplot,
       recurrencerate,
       recurrencestructures,
       dl_average, dl_max, dl_entropy,
       vl_average, vl_max, vl_entropy,
       rt_average, rt_max, rt_entropy,
       determinism,
       divergence,
       trend,
       laminarity,
       trappingtime,
       meanrecurrencetime,
       nmprt,
       rqa,
       sorteddistances,
       skeletonize,
       @windowed,
       rna

include("matrices/distance_matrix.jl")
include("matrices/matrices.jl")
include("matrices/plot.jl")
include("matrices/skeletonization.jl")
include("rqa/histograms.jl")
include("rqa/rqa.jl")
include("rqa/radius.jl")
include("rna/graphs.jl")
include("rna/rna.jl")
include("rqa/windowed.jl")
include("deprecate.jl")

end


text:
module ChaosGameFractals

# Create chaosgame-like fractals
# Translated from https://github.com/python/performance/blob/2ebc46c24c5f816234d8026698470476a3c59d2e/performance/benchmarks/bm_chaos.py
# to Julia
# Original Copyright (C) 2005 Carl Friedrich Bolz

if VERSION >= v"0.7.0-DEV.3052"
    using Printf
end
if VERSION >= v"0.7.0-DEV.3406"
    using Random
end

const DEFAULT_THICKNESS = 0.25
const DEFAULT_WIDTH = 256
const DEFAULT_HEIGHT = 256
const DEFAULT_ITERATIONS = 500000 # 100x more iterations than Python
const DEFAULT_RNG_SEED = 1234

###########
# GVector #
###########

struct GVector{T <: Number}
    x::T
    y::T
    z::T
end

GVector(x=0.0, y=0.0, z=0.0) = GVector(x, y, z)

Base.abs(u::GVector) = sqrt(u.x^2 + u.y^2 + u.z^2)
dist(u::GVector, v::GVector) = sqrt((u.x - v.x)^2 +
                                    (u.y - v.y)^2 +
                                    (u.z - v.z)^2)
Base.:+(u::GVector, v::GVector) = GVector(u.x + v.x, u.y + v.y, u.z + v.z)
Base.:-(u::GVector, v::GVector) = u + (v * -1)
Base.:*(u::GVector, n::Number)= GVector(u.x * n, u.y * n, u.z * n)

function linear_combination(u::GVector, v::GVector, l1, l2=1-l1)
    return GVector(u.x * l1 + v.x * l2,
                   u.y * l1 + v.y * l2,
                   u.z * l1 + v.z * l2,)
end

Base.show(io::IO, u::GVector) = print(io, "<", u.x, ", ", u.y, ", ", u.z, ">")

function getknots(points::Vector{V}, degree::Int) where V <: GVector
    knots = [zeros(Int, degree); 1:(length(points) - degree)]
    append!(knots, fill(length(points) - degree, degree))
    return knots
end

##########
# Spline #
##########

struct Spline{T}
    knots::Vector{Int}
    points::Vector{GVector{T}}
    degree::Int
    d::Vector{GVector{T}}
end

function Spline(points, degree=3, knots = getknots(points, degree))
    if length(points) >  length(knots) - degree + 1
        throw(ArgumentError("too many control points"))
    elseif length(points) < length(knots) - degree + 1
        throw(ArgumentError("not enough control points"))
    elseif any(i -> knots[i] < knots[i-1], 2:length(knots))
        throw(ArgumentError("knots not strictly increasing"))
    end
    Spline(knots, points, degree, similar(points))
end

Base.length(s::Spline) = length(s.points)
getdomain(s::Spline) = s.knots[s.degree], s.knots[length(s.knots) - s.degree + 1]

function Base.getindex(s::Spline, u::Number)
    dom = getdomain(s)
    I = dom[2] - 1
    for ii in (s.degree - 1):(length(s.knots) - s.degree - 1)
        if u >= s.knots[ii + 1] && u < s.knots[ii + 2]
            I = ii
            break
        end
    end
    return I
end

function (s::Spline)(u::Number)
    dom = getdomain(s)
    if u < dom[1] || u > dom[2]
        throw(ArgumentError("Function value not in domain"))
    elseif u == dom[1]
        return s.points[1]
    elseif u == dom[2]
        return s.points[end]
    end
    I = getindex(s, u)
    @inbounds for ii in 0:s.degree
        s.d[ii+1] = s.points[I - s.degree + ii + 2]
    end
    U = s.knots
    @inbounds for ik in 1:s.degree
        for ii in (I - s.degree + ik + 1):(I + 1)
            ua = U[ii + s.degree - ik + 1]
            ub = U[ii]
            co1 = (ua - u) / (ua - ub)
            co2 = (u - ub) / (ua - ub)
            index = ii - I + s.degree - ik - 1
            s.d[index+1] = linear_combination(s.d[index+1], s.d[index + 2], co1, co2)
        end
    end
    return s.d[1]
end

#############
# ChaosGame #
#############

struct ChaosGame{T <: Number}
    splines::Vector{Spline{T}}
    thickness::T
    minx::T
    miny::T
    maxx::T
    maxy::T
    height::T
    width::T
    num_trafos::Vector{Int}
    num_total::Int
end

function ChaosGame(splines, thickness=0.1)
    minx = minimum(p.x for spl in splines for p in spl.points)
    miny = minimum(p.y for spl in splines for p in spl.points)
    maxx = maximum(p.x for spl in splines for p in spl.points)
    maxy = maximum(p.y for spl in splines for p in spl.points)
    height = maxy - miny
    width = maxx - minx
    num_trafos = Int[]
    maxlen = thickness * width / height
    for spl in splines
        len = 0
        curr = spl(0)
        for i in 1:999
            last = curr
            t = 1 / 999 * i
            curr = spl(t)
            len += dist(curr, last)
        end
        push!(num_trafos, max(1, trunc(Int, len / maxlen * 1.5)))
    end
    num_total = sum(num_trafos)
    return ChaosGame(splines, thickness, minx, miny, maxx, maxy, height, width, num_trafos, num_total)
end

function transform_point(point::GVector, game::ChaosGame, trafo = get_random_trafo(game))
    x = (point.x - game.minx) / game.width
    y = (point.y - game.miny) / game.height
    first_trafo = trafo[1] + 1
    start, stop = getdomain(game.splines[first_trafo])

    len = stop - start
    seg_len = len / game.num_trafos[first_trafo]
    t = start + seg_len * trafo[2] + seg_len * x
    basepoint = game.splines[first_trafo](t)
    if t + 1/50000 > stop
        neighbour = game.splines[first_trafo](t - 1/50000)
        derivative = neighbour - basepoint
    else
        neighbour = game.splines[first_trafo](t + 1/50000)
        derivative = basepoint - neighbour
    end

    if abs(derivative) != 0
        fact = 1 / abs(derivative) * (y - 0.5) * game.thickness
        basepoint = GVector(basepoint.x + derivative.y * fact,
                            basepoint.y - derivative.x * fact,
                            basepoint.z)
    end
    return truncate(basepoint, game)
end

function get_random_trafo(game::ChaosGame)
    r = rand(0:game.num_total)
    l = 0
    for i in 0:(length(game.num_trafos)-1)
        if r >= l && r < l + game.num_trafos[i+1]
            return i, rand(0:game.num_trafos[i+1]-1)
        end
        l += game.num_trafos[i+1]
    end
    return length(game.num_trafos) - 1, rand(0:game.num_trafos[end] - 1)
end

function Base.truncate(u::GVector, game::ChaosGame)
    GVector(clamp(u.x, game.minx, game.maxx),
            clamp(u.y, game.miny, game.maxy),
            u.z)
end

function create_image_chaos(game::ChaosGame, w = DEFAULT_HEIGHT,
                            h = DEFAULT_WIDTH, iterations = DEFAULT_ITERATIONS;
                            filename = nothing, seed=DEFAULT_RNG_SEED)

    if VERSION <= v"0.7.0-beta2.171"
        srand(seed)
    else
        Random.seed!(seed)
    end
    im = ones(Int, h, w)
    point = GVector((game.maxx + game.minx) / 2, (game.maxy + game.miny) / 2, 0.0)
    for _ in 1:iterations
        point = transform_point(point, game)

        x = (point.x - game.minx) / game.width * w
        y = (point.y - game.miny) / game.height * h
        xi = trunc(Int, x)
        yi = trunc(Int, y)
        xi == w && (xi -= 1)
        yi == h && (yi -= 1)
        im[h - yi, xi + 1] = 0
    end

    if filename != nothing
        write_ppm(im, filename)
    end
    return im
end

function write_ppm(im, filename::String)
    magic = "P6\n"
    maxval = 255
    h, w = size(im)

    open(filename, "w") do f
        write(f, magic)
        write(f, @sprintf "%i %i\n%i\n" w h maxval)
        for j in 1:w, i in 1:h
            c = round(Int, im[i, j] * 255)
            write(f, UInt8(c), UInt8(c), UInt8(c))
        end
    end
end

const splines = [
        Spline([
            GVector(1.597350, 3.304460, 0.000000),
            GVector(1.575810, 4.123260, 0.000000),
            GVector(1.313210, 5.288350, 0.000000),
            GVector(1.618900, 5.329910, 0.000000),
            GVector(2.889940, 5.502700, 0.000000),
            GVector(2.373060, 4.381830, 0.000000),
            GVector(1.662000, 4.360280, 0.000000)],
            3, [0, 0, 0, 1, 1, 1, 2, 2, 2]),
        Spline([
            GVector(2.804500, 4.017350, 0.000000),
            GVector(2.550500, 3.525230, 0.000000),
            GVector(1.979010, 2.620360, 0.000000),
            GVector(1.979010, 2.620360, 0.000000)],
            3, [0, 0, 0, 1, 1, 1]),
        Spline([
            GVector(2.001670, 4.011320, 0.000000),
            GVector(2.335040, 3.312830, 0.000000),
            GVector(2.366800, 3.233460, 0.000000),
            GVector(2.366800, 3.233460, 0.000000)],
            3, [0, 0, 0, 1, 1, 1])
    ]

const game = ChaosGame(splines)

perf_chaos() = create_image_chaos(game)

end

ext: jl


text:
include("spline_kernel.jl")

# Removing t = 0, such that Σ is invertible
t = Vector(0.1:0.1:1); p = 2;

# Creating generators U,V that result in a positive-definite matrix Σ
U, V = spline_kernel(t, p);
Σ    = spline_kernel_matrix(U, V);
chol = cholesky(Σ)

# Creating a symmetric extended generator representable semiseperable matrix
K = SymSemiseparable(U,V)
# Calculating its Cholesky factorization
L = SymSemiseparableChol(K)
# Creating a test vector
xt = ones(K.n);

# Testing multiplication
@test isapprox(L*xt, chol.L*xt, atol=1e-6)
@test isapprox(L'*xt, chol.U*xt, atol=1e-6)

# Testing inverses (Using Cholesky factorizations)
@test isapprox(L'\(L\xt), chol.U\(chol.L\xt), atol=1e-6)
@test isapprox(inv(L), chol.U\(chol.L\Diagonal(ones(K.n))), atol=1e-6)
@test isapprox(inv(L,xt),  chol.U\(chol.L\xt), atol=1e-6)

# Testing logdet
@test isapprox(logdet(L), logdet(chol.L), atol=1e-10)
@test isapprox(logdet(L'), logdet(chol.U), atol=1e-10)

# Testing the Cholesky as a lower triangular matrix
C = SymSemiseparableChol(U,V);
@test isapprox(C*xt, tril(U*V')*xt)

# Testing going from Cholesky back to the matrix K
Kl = SymSemiseparable(L)
@test isapprox(K*x, Kl*x)


text:
module Aeris

end # module


text:
using EconometricsLH, CategoricalArrays, DataFrames, GLM, StatsBase
using Test

function make_regr_table(startValue :: Float64,  nameV :: Vector{Symbol})
    n = length(nameV);
    coeffV = collect(range(startValue, startValue .+ 3.0, length = n));
    seV = collect(range(startValue .* 0.1, startValue .* 0.4, length = n));
    rt = RegressionTable(nameV, coeffV, seV);
    return rt
end

function make_regr_info(j)
    ri = RegressorInfo(Symbol("r$j"),  0.5 + j,  0.1 + j);
    validate_regressor(ri);
    return ri
end

function regressor_info_test()
    @testset "RegressorInfo" begin
        n = 3;
        rV = [make_regr_info(j)  for j = 1 : n];
        coeffV = [regr_coeff(rV[j])  for j = 1 : n];
        seV = [regr_se(rV[j])  for j = 1 : n];

        rSum = reduce_regr_infos(rV, sum);
        @test regr_coeff(rSum) ≈ sum(coeffV)
        @test regr_se(rSum) ≈ sum(seV)

        r1 = make_regr_info(1);
        r2 = make_regr_info(1);
        @test isapprox(r1, r2)
        scale_regressor(r2, 3.0);
        @test !isapprox(r1, r2)
        @test isapprox(regr_coeff(r1) * 3.0, regr_coeff(r2))
        println(r1);
    end
end

function regression_table_test()
    @testset "RegressionTable" begin
        rt = RegressionTable();
        @test n_regressors(rt) == 0

        nameV = [:c, :b, :a];
        rt = make_regr_table(1.0, nameV);
        coeff2 = get_coefficient(rt, nameV[2])
        se2 = get_std_error(rt, nameV[2])
        @test all(get_coeff_se(rt, nameV[2]) .≈ (coeff2, se2))

        coeff2V, se2V = get_coeff_se_multiple(rt, nameV[1:2]);
        @test coeff2V[2] ≈ coeff2
        @test se2V[2] ≈ se2

        # The order is indeterminate
        coeff2V, se2V = get_coeff_se_multiple(rt, nameV);
        name3V, coeff3V, se3V = get_all_coeff_se(rt);
        @test all(sort(coeff3V) .≈ sort(coeff2V))
        @test all(sort(se3V) .≈ sort(se2V))

        @test sort(get_names(rt)) == sort(nameV)

        # Rename
        rt3 = make_regr_table(3.0, [:a, :b, :c]);
        rt4 = make_regr_table(3.0, [:a, :b, :c]);
        rename_regressor(rt4, :a, :aNew);
        @test isapprox(get_regressor(rt3, :a), get_regressor(rt4, :aNew))

        # Set missing regressors
        rt5 = make_regr_table(3.0, [:a, :b, :c]);
        set_missing_regressors!(rt3, [:a, :c]);
        @test have_same_regressors([rt3, rt5])
        set_missing_regressors!(rt3, [:a, :d]);
        @test sort(get_names(rt3)) == [:a, :b, :c, :d]

        println(rt)

        # Drop regressor
        drop_regressor!(rt5, :a);
        @test !has_regressor(rt5, :a)
    end
end

function regr_table_mult_test()
    @testset "Multiple RegressionTables" begin
        rt1 = make_regr_table(2.0, [:a, :c, :b]);
        rt2 = make_regr_table(2.0, [:b, :c, :a]);
        @test have_same_regressors([rt1, rt2])
        @test !isapprox(rt1, rt2)

        rt3 = make_regr_table(2.1, [:a, :b, :c, :d]);
        @test !have_same_regressors([rt1, rt2, rt3])

        # Reduce
        rtOut = reduce_regr_tables([rt1, rt2], mean);
        for name in [:a, :b, :c]
            coeffOut, seOut = get_coeff_se(rtOut, name);
            @test coeffOut ≈ (get_coefficient(rt1, name) + get_coefficient(rt2, name)) / 2.0;
            @test seOut ≈ (get_std_error(rt1, name) + get_std_error(rt2, name)) / 2.0
        end

        rtOut2 = mean([rt1, rt2]);
        @test isapprox(rtOut, rtOut2)
        for name in [:a, :b, :c]
            @test get_coefficient(rtOut, name) ≈ get_coefficient(rtOut2, name)
            @test get_std_error(rtOut, name) ≈ get_std_error(rtOut2, name)
        end

        rtOut3 = reduce_regr_tables([rt1, rt2], StatsBase.std);
        rtOut4 = StatsBase.std([rt1, rt2]);
        @test isapprox(rtOut3, rtOut4)
    end
end

function from_lin_model_test()
    @testset "From linear model" begin
        n = 16;
        zV = categorical(round.(Int, collect(range(1, 4, length = n))));
        data = DataFrame(X = collect(range(1.0, 3.0, length = n)) .^ 2,
            Y = collect(range(0.1, 3.5, length = n)),
            Z = zV);
        ols = lm(@formula(Y ~ X + Z), data);
        coeffV = coef(ols);

        rt = RegressionTable(ols);
        @test get_coefficient(rt, :X) ≈ coeffV[2]
        @test get_coefficient(rt, :constant) ≈ coeffV[1]
        @test get_coefficient(rt, :Z2) ≈ coeffV[3]
    end
end

@testset "EconometricsLH" begin
    regressor_info_test();
    regression_table_test();
    regr_table_mult_test();
    from_lin_model_test();
    include("regression_test_test.jl");
end

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


text:
# Define an abstract type to represent color difference metrics
abstract type DifferenceMetric end

abstract type EuclideanDifferenceMetric{T<:Color3} <: DifferenceMetric end

# TODO?: make the DifferenMetrics parametric, to preserve type-stability

struct DE_2000 <: DifferenceMetric
    kl::Float64
    kc::Float64
    kh::Float64
    DE_2000(kl=1, kc=1, kh=1) = new(kl, kc, kh)
end
"""
    DE_2000(kl=1, kc=1, kh=1)

Construct a metric of the CIE Delta E 2000 recommendation, with weighting
parameters `kl`, `kc` and `kh` as provided for in the recommendation. When not
provided, these parameters default to `1`.
"""
DE_2000()

struct DE_94 <: DifferenceMetric
    kl::Float64
    kc::Float64
    kh::Float64
    k1::Float64
    k2::Float64
    DE_94(kl=1, kc=1, kh=1, k1=0.045, k2=0.015) = new(kl, kc, kh, k1, k2)
end
"""
    DE_94(kl=1, kc=1, kh=1, k1=0.045, k2=0.015)

Construct a metric of CIE Delta E 94 recommendation (1994), with weighting
parameters `kl`, `kc`, `kh`, `k1`, and `k2` as provided for in the
recommendation. The `kl`, `k1`, and `k2` depend on the application:

|    |Graphic Arts|Textiles|
|:--:|:-----------|:-------|
|`kl`|`1`         |`2`     |
|`k1`|`0.045`     |`0.048` |
|`k2`|`0.015`     |`0.014` |

and the default values are for graphic arts. The `kc` and `kh` default to `1`.

The `DE_94` is more perceptually uniform than the [`DE_AB`](@ref), but has some
non-uniformities resolved by the [`DE_2000`](@ref).

!!! note
    The `DE_94` is a quasimetric, i.e. violates symmetry. Therefore,
    `colordiff(a, b, metric=DE_94())` may not equal to
    `colordiff(b, a, metric=DE_94())`.
    The first argument of `colordiff` is taken as the reference (standard)
    color.
"""
DE_94()

struct DE_JPC79 <: DifferenceMetric

end
"""
    DE_JPC79()

Construct a metric using McDonald's "JP Coates Thread Company" color difference
formula.
"""
DE_JPC79()

struct DE_CMC <: DifferenceMetric
    kl::Float64
    kc::Float64
    DE_CMC(kl=1, kc=1) = new(kl, kc)
end
"""
    DE_CMC(kl=1, kc=1)

Construct a metric using the CMC equation (CMC l:c), with weighting parameters
`kl` and `kc`. When not provided, these parameters default to `1`.

!!! note
    The `DE_CMC` is a quasimetric, i.e. violates symmetry. Therefore,
    `colordiff(a, b, metric=DE_CMC())` may not equal to
    `colordiff(b, a, metric=DE_CMC())`.
    The first argument of `colordiff` is taken as the reference (standard)
    color.
"""
DE_CMC()

struct DE_BFD <: DifferenceMetric
    wp::XYZ{Float64}
    kl::Float64
    kc::Float64
    DE_BFD(wp::XYZ, kl=1, kc=1) = new(wp, kl, kc)
end
"""
    DE_BFD([wp,] kl=1, kc=1)

Construct a metric using the BFD equation, with weighting parameters `kl` and
`kc`. Additionally, a whitepoint `wp` can be specified, because the BFD equation
must convert between `XYZ` and `Lab` during the computation. When not provided,
`kl` and `kc` default to `1`, and `wp` defaults to CIE D65 (`Colors.WP_D65`).
"""
DE_BFD(kl=1, kc=1) = DE_BFD(WP_DEFAULT, kl, kc)

struct DE_AB <: EuclideanDifferenceMetric{Lab}

end
"""
    DE_AB()

Construct a metric of the original CIE Delta E equation (ΔE*ab), or Euclidean
color difference equation in the `Lab` (CIELAB) colorspace.
"""
DE_AB()

struct DE_DIN99 <: EuclideanDifferenceMetric{DIN99}

end
"""
    DE_DIN99()

Construct a metric using Euclidean color difference equation applied in the
`DIN99` colorspace.
"""
DE_DIN99()

struct DE_DIN99d <: EuclideanDifferenceMetric{DIN99d}

end
"""
    DE_DIN99d()

Construct a metric using Euclidean color difference equation applied in the
`DIN99d` colorspace.
"""
DE_DIN99d()

struct DE_DIN99o <: EuclideanDifferenceMetric{DIN99o}

end
"""
    DE_DIN99o()

Construct a metric using Euclidean color difference equation applied in the
`DIN99o` colorspace.
"""
DE_DIN99o()

# Color difference metrics
# ------------------------

# Evaluate the CIEDE2000 color difference formula, implemented according to:
#   Klaus Witt, CIE Color Difference Metrics, Colorimetry: Understanding the CIE
#   System. 2007
#
# Args:
#   a, b: Any two colors.
#
# Returns:
#   The CIEDE2000 color difference metric evaluated between a and b.
#

# Delta E 2000
function _colordiff(ai::Color, bi::Color, m::DE_2000)
    twentyfive7 = pow7(25)

    # Ensure that the input values are in L*a*b* space
    a_Lab = convert(Lab, ai)
    b_Lab = convert(Lab, bi)

    # Calculate some necessary factors from the L*a*b* values
    mc = (chroma(a_Lab) + chroma(b_Lab))/2
    g = (1 - sqrt(pow7(mc) / (pow7(mc) + twentyfive7))) / 2

    a_Lab_r = Lab(a_Lab.l, a_Lab.a * (1 + g), a_Lab.b)
    b_Lab_r = Lab(b_Lab.l, b_Lab.a * (1 + g), b_Lab.b)

    # Convert to L*C*h, where the remainder of the calculations are performed
    a = convert(LCHab, a_Lab_r)
    b = convert(LCHab, b_Lab_r)

    # Calculate the delta values for each channel
    dl, dc, dh = b.l - a.l, b.c - a.c, delta_h(b_Lab_r, a_Lab_r)

    # Calculate mean L*, C* and hue values
    ml, mc, mh = (a.l + b.l) / 2, (a.c + b.c) / 2, mean_hue(a, b)

    # lightness weight
    mls = (ml - 50)^2
    sl = 1.0 + 0.015 * mls / sqrt(20 + mls)

    # chroma weight
    sc = 1 + 0.045mc

    # hue weight
    t = 1 - 0.17 * cosd(mh - 30) +
            0.24 * cosd(2mh) +
            0.32 * cosd(3mh + 6) -
            0.20 * cosd(4mh - 63)
    sh = 1 + 0.015 * mc * t

    # rotation term
    dtheta = 60 * exp(-((mh - 275)/25)^2)
    cr = 2 * sqrt(pow7(mc) / (pow7(mc) + twentyfive7))
    tr = -sind(dtheta) * cr

    # Final calculation
    sqrt((dl/(m.kl*sl))^2 + (dc/(m.kc*sc))^2 + (dh/(m.kh*sh))^2 +
         tr * (dc/(m.kc*sc)) * (dh/(m.kh*sh)))
end

# Delta E94
function _colordiff(ai::Color, bi::Color, m::DE_94)
    a = convert(LCHab, ai)
    b = convert(LCHab, bi)

    # Calculate the delta values for each channel
    dl, dc, dh = b.l - a.l, b.c - a.c, delta_h(b, a)

    # Lightness, hue, chroma correction terms
    sl = 1
    sc = 1 + m.k1 * a.c
    sh = 1 + m.k2 * a.c

    sqrt((dl/(m.kl*sl))^2 + (dc/(m.kc*sc))^2 + (dh/(m.kh*sh))^2)
end

# Metric form of jpc79 color difference equation (mostly obsolete)
function _colordiff(ai::Color, bi::Color, m::DE_JPC79)
    # Convert directly into LCh
    a = convert(LCHab, ai)
    b = convert(LCHab, bi)

    # Calculate deltas in each direction
    dl, dc, dh = b.l - a.l, b.c - a.c, delta_h(b, a)

    # Calculate mean lightness, chroma and hue
    ml, mc, mh = (a.l + b.l) / 2, (a.c + b.c) / 2, mean_hue(a, b)

    # L* adjustment term
    sl = 0.08195*ml/(1+0.01765*ml)

    # C* adjustment term
    sc = 0.638+(0.0638*mc/(1+0.0131*mc))

    # H* adjustment term
    if (mc < 0.38)
        sh = sc
    elseif (mh >= 164 && mh <= 345)
        sh = sc*(0.56+abs(0.2*cosd(mh+168)))
    else
        sh = sc*(0.38+abs(0.4*cosd(mh+35)))
    end

    # Calculate the final difference
    sqrt((dl/sl)^2 + (dc/sc)^2 + (dh/sh)^2)
end

# Metric form of the cmc color difference
function _colordiff(ai::Color, bi::Color, m::DE_CMC)
    # Convert directly into LCh
    a = convert(LCHab, ai)
    b = convert(LCHab, bi)

    # Calculate deltas in each direction
    dl, dc, dh = b.l - a.l, b.c - a.c, delta_h(b, a)

    # L* adjustment term
    if (a.l <= 16)
        sl = 0.511
    else
        sl = 0.040975 * a.l / (1 + 0.01765 * a.l)
    end

    # C* adjustment term
    sc = 0.0638 * a.c / (1 + 0.0131 * a.c) + 0.638

    f = sqrt(a.c^4 / (a.c^4 + 1900))

    if 164 <= a.h <= 345
        t = 0.56 + abs(0.2 * cosd(a.h + 168))
    else
        t = 0.36 + abs(0.4 * cosd(a.h + 35))
    end

    # H* adjustment term
    sh = sc*(t*f+1-f)

    # Calculate the final difference
    sqrt((dl/(m.kl*sl))^2 + (dc/(m.kc*sc))^2 + (dh/sh)^2)
end

# The BFD color difference equation
function _colordiff(ai::Color, bi::Color, m::DE_BFD)
    # Currently, support for the `wp` argument of `convert` is limited.
    function to_xyz(c::Color, wp)
        c isa XYZ && return c
        (c isa xyY || c isa LMS) && return convert(XYZ, c)
        (c isa Lab || c isa Luv) && return convert(XYZ, c, wp)
        c isa LCHuv && return convert(XYZ, convert(Luv, c), wp)
        convert(XYZ, convert(Lab, c), wp)
    end

    # We have to start back in XYZ because BFD uses a different L equation
    a_XYZ = to_xyz(ai, m.wp)
    b_XYZ = to_xyz(bi, m.wp)

    la = 54.6*log10(a_XYZ.y+1.5)-9.6
    lb = 54.6*log10(b_XYZ.y+1.5)-9.6

    # Convert into LCh with the proper white point
    a_Lab = convert(Lab, a_XYZ, m.wp)
    b_Lab = convert(Lab, b_XYZ, m.wp)

    # Substitute in the different L values into the L*C*h values
    a = LCHab(la, chroma(a_Lab), hue(a_Lab))
    b = LCHab(lb, chroma(b_Lab), hue(b_Lab))

    # Calculate deltas in each direction
    dl, dc, dh = b.l - a.l, b.c - a.c, delta_h(b, a)

    # Find the mean value of the inputs to use as the "standard"
    ml, mc, mh = (a.l + b.l) / 2, (a.c + b.c) / 2, mean_hue(a, b)

    # Correction terms for a variety of nonlinearities in CIELAB.
    g = sqrt(mc^4/(mc^4 + 14000))

    t = 0.627 + 0.055 * cosd( mh - 245) -
                0.040 * cosd(2mh - 136) +
                0.070 * cosd(3mh -  32) +
                0.049 * cosd(4mh + 114) -
                0.015 * cosd(5mh + 103)

    rc = sqrt(mc^6/(mc^6 + 7e7))

    rh = -0.260 * cosd( mh - 308) -
          0.379 * cosd(2mh - 160) -
          0.636 * cosd(3mh - 254) +
          0.226 * cosd(4mh + 140) -
          0.194 * cosd(5mh + 280)

    dcc = 0.035*mc/(1+0.00365*mc) + 0.521
    dhh = dcc*(g*t+1-g)
    rt = rc*rh

    # Final calculation
    sqrt((dl/m.kl)^2 + (dc/(m.kc*dcc))^2 + (dh/dhh)^2 + rt*((dc*dh)/(dcc*dhh)))
end

function _colordiff(ai::Color, bi::Color,
                    m::EuclideanDifferenceMetric{T}) where {T <: Color3}
    a, b = convert(T, ai), convert(T, bi)

    d1, d2, d3 = comp1(a) - comp1(b), comp2(a) - comp2(b), comp3(a) - comp3(b)

    sqrt(d1^2 + d2^2 + d3^2)
end

# Default to Delta E 2000
"""
    colordiff(a, b; metric=DE_2000())

Compute an approximate measure of the perceptual difference between colors `a`
and `b`. Optionally, a `metric` may be supplied, chosen among [`DE_2000`](@ref)
(the default), [`DE_94`](@ref), [`DE_JPC79`](@ref), [`DE_CMC`](@ref),
[`DE_BFD`](@ref), [`DE_AB`](@ref), [`DE_DIN99`](@ref), [`DE_DIN99d`](@ref) and
[`DE_DIN99o`](@ref).

The return value is a non-negative number in a type depending on the colors and
metric.

!!! note
    The supported metrics measure the difference within `Lab` or its variant
    colorspaces. When the input colors are not in the colorspace internally used
    by the metric, the colors (e.g. in `RGB`) are converted with the default
    whitepoint CIE D65 (`Colors.WP_D65`). If you want to use another whitepoint,
    convert the colors into the colorspace used by metric (e.g. `Lab` for
    [`DE_2000`](@ref)) in advance.
"""
colordiff(ai::Union{Number, Color},
          bi::Union{Number, Color};
          metric::DifferenceMetric=DE_2000()) = _colordiff(ai, bi, metric)
@deprecate colordiff(ai::Color, bi::Color, metric::DifferenceMetric) colordiff(ai, bi; metric=metric)

function colordiff(ai::Colorant, bi::Colorant; metric::DifferenceMetric=DE_2000())
    alpha(ai) == 1 && alpha(bi) == 1 && return _colordiff(color(ai), color(bi), metric)
    throw(ArgumentError("""
        cannot evaluate the difference in transparent colors.
          Their appearance depends on the backdrop."""))
end

_colordiff(ai::AbstractGray, bi::Number, metric::DifferenceMetric) = _colordiff(ai, Gray(bi), metric)
_colordiff(ai::Number, bi::AbstractGray, metric::DifferenceMetric) = _colordiff(Gray(ai), bi, metric)
_colordiff(ai::Number, bi::Number, metric::DifferenceMetric) = _colordiff(Gray(ai), Gray(bi), metric)


text:
using Test, PolynomialBases
using LinearAlgebra

ufunc(x) = sinpi(x)
uprim(x) = π*cospi(x)

function tolerance(p, T=Float64)
    if p <= 2
        5.
    elseif p <= 6
        0.5
    elseif p <= 8
        0.04
    elseif p <= 10
        0.002
    elseif p <= 12
        if T==Float32
            5.e-5
        else
            3.e-5
        end
    elseif p <= 14
        if T==Float32
            3.e-5
        else
            6.e-7
        end
    elseif p <= 16
        if T==Float32
            4.e-5
        else
            7.e-9
        end
    else
        if T==Float32
            6.e-5
        else
            7.e-11
        end
    end
end

# regression test
for basis_type in subtypes(PolynomialBases.NodalBasis{PolynomialBases.Line})
    println(devnull, "  ", basis_type(1))
    basis_type <: ClosedNewtonCotes && continue
    for p in 5:20, T in (Float32, Float64)
        basis = basis_type(p, T)
        u = ufunc.(basis.nodes)
        @test basis.D * u ≈ uprim.(basis.nodes) atol=tolerance(p, T)
    end
end

# compare direct derivative evaluation and matrix
for basis_type in subtypes(PolynomialBases.NodalBasis{PolynomialBases.Line})
    for p in 0:15
        basis = try
            basis_type(p)
        catch m
            isa(m, DomainError) && continue
            throw(m)
        end
        u = ufunc.(basis.nodes)
        xplot = range(-1, stop=1, length=100)
        u1 = derivative_at(xplot, u, basis)
        u2 = interpolate(xplot, basis.D*u, basis)
        @test u1 ≈ u2
    end

    for p in 0:15
        basis = try
            basis_type(p, Float32)
        catch m
            isa(m, DomainError) && continue
            throw(m)
        end
        u = ufunc.(basis.nodes)
        xplot = range(-1, stop=1, length=100)
        u1 = derivative_at(xplot, u, basis)
        u2 = interpolate(xplot, basis.D*u, basis)
        @test norm(u1 - u2, Inf) < 5.f-3
    end
end


text:
@testset "Debug Options" begin
    # helper to get debug as string
    io = IOBuffer()
    M = ManifoldsBase.DefaultManifold(2)
    x = [4.,2.]
    o = GradientDescentOptions(x, StopAfterIteration(20), ConstantStepsize(1.))
    f = y -> distance(M,y,x).^2
    ∇f = y -> -2*log(M,y,x)
    p = GradientProblem(M,f,∇f)
    a1 = DebugDivider("|",x -> print(io,x))
    @test Manopt.dispatch_options_decorator(DebugOptions(o,a1)) === Val{true}()
    # constructors
    @test DebugOptions(o,a1).debugDictionary[:All] == a1
    @test DebugOptions(o,[a1]).debugDictionary[:All].group[1] == a1
    @test DebugOptions(o,Dict(:A => a1)).debugDictionary[:A] == a1
    @test DebugOptions(o, ["|"]).debugDictionary[:All].group[1].divider == a1.divider
    # single AbstractOptionsActions
    # DebugDivider
    a1(p,o,0); s =
    @test String(take!(io)) == "|"
    #DebugGroup
    DebugGroup([a1,a1])(p,o,0)
    @test String(take!(io)) == "||"
    # Debug Every
    DebugEvery(a1,10,false)(p,o,9)
    @test String(take!(io)) == ""
    DebugEvery(a1,10,true)(p,o,10)
    @test String(take!(io)) == "|"
    @test DebugEvery(a1,10,true)(p,o,-1) == nothing
    # Debug Cost
    @test DebugCost("A").prefix == "A"
    DebugCost(false,x -> print(io,x))(p,o,0)
    @test String(take!(io)) == "F(x): 0.0"
    DebugCost(false,x -> print(io,x))(p,o,-1)
    @test String(take!(io)) == ""
    # entry
    DebugEntry(:x,"x:",x->print(io,x))(p,o,0)
    @test String(take!(io)) == "x: $x"
    DebugEntry(:x,"x:",x->print(io,x))(p,o,-1)
    @test String(take!(io)) == ""
    # Change
    a2 = DebugChange(StoreOptionsAction( (:x,) ), "Last: ", x -> print(io,x))
    a2(p,o,0) # init
    o.x = [3.,2.]
    a2(p,o,1)
    @test String(take!(io)) == "Last: 1.0"
    # Iterate
    DebugIterate(x -> print(io,x))(p,o,0)
    @test String(take!(io)) == "x:$(o.x)"
    DebugIterate(x -> print(io,x))(p,o,1)
    @test String(take!(io)) == "x:$(o.x)"
    # Iteration
    DebugIteration(x -> print(io,x))(p,o,0)
    @test String(take!(io)) == "Initial"
    DebugIteration(x -> print(io,x))(p,o,23)
    @test String(take!(io)) == "# 23"
    # DEbugEntryChange - reset
    o.x = x
    a3 = DebugEntryChange(:x,(p,o,x,y)-> distance(p.M,x,y), StoreOptionsAction( (:x,) ), "Last: ", x -> print(io,x))
    a4 = DebugEntryChange(x,:x,(p,o,x,y)-> distance(p.M,x,y), StoreOptionsAction( (:x,) ), "Last: ", x -> print(io,x))
    a3(p,o,0) # init
    @test String(take!(io)) == ""
    a4(p,o,0) # init
    @test String(take!(io)) == ""
    #change
    o.x = [3.,2.]
    a3(p,o,1)
    @test String(take!(io)) == "Last: 1.0"
    a4(p,o,1)
    @test String(take!(io)) == "Last: 1.0"
    # StoppingCriterion
    DebugStoppingCriterion(x -> print(io,x))(p, o, 1)
    @test String(take!(io)) == ""
    o.stop(p,o,20)
    DebugStoppingCriterion(x -> print(io,x))(p, o, 20)
    @test String(take!(io)) == ""
    o.stop(p,o,21)
    DebugStoppingCriterion(x -> print(io,x))(p, o, 21)
    @test String(take!(io)) == "The algorithm reached its maximal number of iterations (20).\n"
    #DebugFactory
    df = DebugFactory([:Stop,"|"])
    @test isa( df[:Stop], DebugStoppingCriterion)
    @test isa(df[:All], DebugGroup)
    @test isa( df[:All].group[1], DebugDivider )
    @test length(df[:All].group) == 1
    df = DebugFactory([:Stop,"|",20])
    @test isa(df[:All], DebugEvery)
    @test all(  isa.( DebugFactory([:Change, :Iteration, :Iterate, :Cost, :x])[:All].group,
        [DebugChange, DebugIteration, DebugIterate, DebugCost, DebugEntry]) )
    @test DebugActionFactory(a3) == a3
end


text:
#Wrappers for functions directly in the rospy namespace
export init_node, is_shutdown, spin,
       get_param, has_param, set_param, delete_param,
       logdebug, loginfo, logwarn, logerr, logfatal

"""
    init_node(name; args...)

Initialize this node, registering it with the ROS master. All arguments are passed on directly to
the rospy init_node function.
"""
init_node(name::AbstractString; args...) =
    __rospy__[:init_node](ascii(name); args...)

"""
    is_shutdown()

Return the shutdown status of the node.
"""
is_shutdown()          = __rospy__[:is_shutdown]()

get_published_topics() = __rospy__[:get_published_topics]()
get_ros_root()         = __rospy__[:get_ros_root]()

"""
    spin()

Block execution and process callbacks/service calls until the node is shut down.
"""
function spin()
    #Have to make sure both Julia tasks and python threads can wake up so
    #can't just call rospy's spin
    while ! is_shutdown()
        rossleep(Duration(0.001))
    end
end

#Parameter server API
"""
    get_param(param_name, default=nothing)

Request the value of a parameter from the parameter server, with optional default value. If no
default is given, throws a `KeyError` if the parameter cannot be found.
"""
function get_param(param_name::AbstractString, def=nothing)
    try
        if def == nothing
            __rospy__[:get_param](ascii(param_name))
        else
            __rospy__[:get_param](ascii(param_name), def)
        end
    catch ex
        throw(KeyError(pycall(pybuiltin("str"), PyAny, ex.val)[2:end-1]))
    end
end

"""
    set_param(param_name, val)

Set the value of a parameter on the parameter server.
"""
set_param(param_name::AbstractString, val) =
    __rospy__[:set_param](ascii(param_name), val)

"""
    has_param(param_name)

Return a boolean specifying if a parameter exists on the parameter server.
"""
has_param(param_name::AbstractString) =
    __rospy__[:has_param](ascii(param_name))

"""
    delete_param(param_name)

Delete a parameter from the parameter server. Throws a `KeyError` if no such parameter exists.
"""
function delete_param(param_name::AbstractString)
    try
        __rospy__[:delete_param](ascii(param_name))
    catch ex
        throw(KeyError(pycall(pybuiltin("str"), PyAny, ex.val)[2:end-1]))
    end
end

#Doesn't work for some reason
#rospy_search_param(param_name::AbstractString) =
#    __rospy__[:rospy_search_param](ascii(param_name))
get_param_names() = __rospy__[:get_param_names]()

#Logging API
logdebug(msg, args...) = __rospy__[:logdebug](msg, args...)
loginfo(msg, args...)  = __rospy__[:loginfo](msg, args...)
logwarn(msg, args...)  = __rospy__[:logwarn](msg, args...)
logerr(msg, args...)   = __rospy__[:logerr](msg, args...)
logfatal(msg, args...) = __rospy__[:logfatal](msg, args...)

"""
    logdebug, loginfo, logwarn, logerr, logfatal

Call the rospy logging system at the corresponding message level, passing a message and other
arguments directly.
"""
logdebug, loginfo, logwarn, logerr, logfatal

#Node information
get_name()             = __rospy__[:get_name]()
get_namespace()        = __rospy__[:get_namespace]()
get_node_uri()         = __rospy__[:get_node_uri]()
get_caller_id()        = __rospy__[:get_caller_id]()


text:
struct DifferentBC
    left
    right
end

function calcLeftFlux(bc::DifferentBC,
                      state::HydroSnapshot,
                      edge_velocities::Array{Float64,1},
                      eos,
                      cached)::Conserved
    calcLeftFlux(bc.left,
                 state,
                 edge_velocities,
                 eos,
                 cached)
end

function calcRightFlux(bc::DifferentBC,
                       state::HydroSnapshot,
                       edge_velocities::Array{Float64,1},
                       eos,
                       cached)::Conserved
    calcRightFlux(bc.right,
                  state,
                  edge_velocities,
                  eos,
                  cached)
end

max_stars_repo_name: bolverk/Sirion

max_issues_repo_name: bolverk/Sirion

max_forks_repo_name: bolverk/Sirion


text:
"""
    function get_Pg_thermal_constprice(args..., kwargs...)

Identify generators which cost function is a constant price, locate their buses and returns
a matrix with their Name, Power output and Bus number (sorted by Bus). Skips Slack
generators, Generators which Active Power is binding and Generators with non zero slope.
Note: For a quadratic cost function of the form: αPg²+ βPg + γ, a constant price generator
will have the parameter α == 0.

# Arguments
- `sys::System`:                            Power system in p.u. (from PowerSystems.jl)
- `res::PowerSimulations.
    OperationsProblemResults`:              Results of the solved OPF for the system
                                             (from PowerSimulations.jl)

# Keywords
- `dual_gen_tol::Float64 = 1e-1`:           Tolerance to identify any binding generators
"""
function get_Pg_thermal_constprice(
    sys::System,
    res::PowerSimulations.OperationsProblemResults;
    dual_gen_tol::Float64 = 1e-1
    )
    # Define elements
    gen_loc = 0
    gen_constprice = 0 #number of generators with constant price (counter starter)
    gens_thermal_constprice_Pg = Array{Float64}(undef,0)
    gens_thermal_constprice_names = Array{String}(undef,0)
    gens_thermal_constprice_busnumbers = Array{Int64}(undef,0)

    # Get all thermal Generators
    gens_thermal = get_components(ThermalStandard, sys)
    all_PGenThermal = get_variables(res)[:P__ThermalStandard] #Optimised PGen

    # Get dual variables from the OPF results
    dualPgub = get_duals(res)[:P_ub__ThermalStandard__RangeConstraint]
    dualPglb = get_duals(res)[:P_lb__ThermalStandard__RangeConstraint]

    # Identify Generators and save their information
    for gen_thermal in gens_thermal
        gen_loc = gen_loc +1
        if get_bustype(gen_thermal.bus) ≠ BusTypes.REF  #skips slack
            bindub = !isapprox(dualPgub[1,gen_loc], 0.0; atol = dual_gen_tol) #is binding?
            bindlb = !isapprox(dualPglb[1,gen_loc], 0.0; atol = dual_gen_tol) #is binding?
            if !bindub && !bindlb #skips binding Pg gens Dual = 0
                α, β = get_cost(get_variable(get_operation_cost(gen_thermal)))
                γ = get_fixed(get_operation_cost(gen_thermal))
                if α == 0 && (β ≠ 0 || γ ≠ 0)#skips non zero slope
                    gen_constprice = gen_constprice + 1
                    gen_thermal_constprice_name = gen_thermal.name
                    resize!(gens_thermal_constprice_Pg,gen_constprice)
                    resize!(gens_thermal_constprice_names,gen_constprice)
                    resize!(gens_thermal_constprice_busnumbers,gen_constprice)
                    gens_thermal_constprice_Pg[gen_constprice] =
                        all_PGenThermal[1, gen_thermal_constprice_name]
                    gens_thermal_constprice_names[gen_constprice] =
                        gen_thermal_constprice_name
                    gens_thermal_constprice_busnumbers[gen_constprice] =
                        gen_thermal.bus.number
                end
            end
        end
    end

    # Matrix of generators with their information sorted
    gens_thermal_constprice_Matrix = Array{Any}(undef,(gen_constprice, 3))
    gens_thermal_constprice_Matrix[:, 1] = gens_thermal_constprice_names
    gens_thermal_constprice_Matrix[:, 2] = gens_thermal_constprice_Pg
    gens_thermal_constprice_Matrix[:, 3] = gens_thermal_constprice_busnumbers
    # Sort by bus to match with PTDF structure
    gens_thermal_constprice_Matrix =
        gens_thermal_constprice_Matrix[sortperm(gens_thermal_constprice_Matrix[:, 3]), :]

    return gens_thermal_constprice_Matrix
end

max_stars_repo_name: invenia/MaximumGeneratorProfitTCRDD

max_issues_repo_name: invenia/MaximumGeneratorProfitTCRDD

max_forks_repo_name: invenia/MaximumGeneratorProfitTCRDD


text:
import FractionalCalculus.FracIntAlg

"""
Riemann-Liouville sense fractional integral algorithms

Note this two algorithms belong to direct compute, precise are ensured, but maybe cause more memory allocation and take more compilation time.
"""
abstract type RLInt <: FracIntAlg end

"""
# Riemann Liouville sense fractional integral

    fracint(f::Function, α, start_point, end_point, h, RL_Direct())

Riemann_Liouville fractional integral using complex step differentiation.

```math
(J^αf)(t)=\\frac{1}{\\Gamma(α)} \\int_0^t(t-τ)^{α-1}f(τ)dτ
```

By using [QuadGK](https://github.com/JuliaMath/QuadGK.jl) calculate the integration and obtain the value.

### Example:

```julia-repl
julia> fracint(x->x^5, 0.5, 0, 2.5, 1e-8, RL_Direct())
```

Returns a tuple (1.1639316474512205, 1.0183453796725215e-8), which contains the value of this derivative is 1.1639316474512205, and the error estimate is 1.0183453796725215e-8
"""
struct RL_Direct <: RLInt end

"""
# Riemann Liouville sense fractional integral with first diff known.

    fracint(f, fd, α, start_point, end_point, RL_Direct_First_Diff_Known())

### Example

```julia-repl
julia> fracint(x->x^5, x->5*x^4, 0.5, 0, 2.5, RL_Direct_First_Diff_Known())
```

With first derivative known, we can use the Riemann Liouville sense to obtain the fractional integral more effcient.
"""
struct RL_Direct_First_Diff_Known <: RLInt end

"""
@article{LI20113352,
title = {Numerical approaches to fractional calculus and fractional ordinary differential equation},
author = {Changpin Li and An Chen and Junjie Ye},
}
"""

"""
# Riemann Liouville sense fractional integral using piecewise interpolation.

    fracint(f, α, end_point, h, RL_Piecewise())

### Example

```julia-repl
julia> fracint(x->x^5, 0.5, 2.5, 0.0001, RL_Piecewise())
```

By deploying Piecewise interpolation to approximate the original function, with small step size, this method is fast and take little memory allocation.

Using piecewise linear interpolation:

```math
    y_n(t)=\\frac{t-t_{i+1}}{t_i-t_{i+1}}y(t_i)+\\frac{t-t_i}{t_{i+1}-t_i}y(t_{i+1})
```

Constitute the original function with the interpolation and implement the summation.
"""
struct RL_Piecewise <: RLInt end

"""
@book{oldham_spanier_1984,
title={The fractional calculus: Theory and applications of differentiation and integration to arbitrary order},
author={Oldham, Keith B. and Spanier, Jerome},
year={1984}}
"""

"""
# Riemann Liouville sense fractional integral approximation.

    fracint(f, α, end_point, h, RLInt_Approx())

### Example

```julia-repl
julia> fracint(x->x^5, 0.5, 2.5, 0.0001, RLInt_Approx())
```

Using the **Staircase approximation** to approximate the original function and implement the summation.
"""
struct RLInt_Approx <: RLInt end

"""
# Riemann Liouville sense fractional integral using **Linear interpolation**.

    fracint(f, α, end_point, h, RL_LinearInterp())

### Example

```julia-repl
julia> fracint(x->x^5, 0.5, 2.5, 0.0001, RL_LinearInterp())
```

**RL_LinearInterp** is more complex compared with *RLInt_Approx* but more precise.

Deploying the [Linear Interpolation](https://en.wikipedia.org/wiki/Linear_interpolation) between ``f_{j+1}`` and ``f_j``, **RL_LinearInterp** method is more precise than **RLInt_Approx**.
"""
struct RL_LinearInterp <: RLInt end

"""
@article{2009,
title={Matrix approach to discrete fractional calculus II: Partial fractional differential equations},
DOI={10.1016/j.jcp.2009.01.014},
author={Podlubny, Igor and Chechkin, Aleksei and Skovranek, Tomas and Chen, YangQuan and Vinagre Jara, Blas M.},
}

"""

"""
# Riemann Liouville sense integral using Triangular Strip Matrix to discrete.

    fracint(f, α, end_point, h, RLInt_Matrix())

Using Triangular Strip Matrix to approximate fractional integral.

### Example

```julia-repl
julia> fracint(x->x^5, 0.5, 2.5, 0.0001, RLInt_Matrix())
```

!!! info
    Triangular Strip Matrix method returns the derivative in the interval ``[0, T]`` in ```Vector```

!!! tip
    With the advancing Triangular Strip Matrix method, you can not only compute fractional integrals, integer order, higher order integral is also supported!!
Try to set α as an integer, arbitrary integer of course! I promise you would enjoy it😏

Using [Triangular Strip Matrix](https://en.wikipedia.org/wiki/Triangle_strip) to discrete the integral.
"""
struct RLInt_Matrix <: RLInt end

"""
# Riemann Liouville sense fractional integral using fractional Simpson's formula to approximate.

### Example

```julia-repl
julia> fracint(x->x, 0.5, 1, 0.001, RLInt_Simpson)
```

@inproceedings{Li2015NumericalMF,
  title={Numerical Methods for Fractional Calculus},
  author={Changpin Li and Fanhai Zeng},
  year={2015}
}

Using fractional Simpson's formula to discrete fractional integral.
"""
struct RLInt_Simpson <: RLInt end

"""
# Riemann Liouville sense fractional integral using fractional trapezoidal formula to approximate.

### Example

```julia-repl
julia> fracint(x->x, 0.5, 1, 0.001, RLInt_Trapezoidal())
```

《Numerical methods for fractional calculus》.
Using Trapezoidal method to approximate fractional integral
"""
struct RLInt_Trapezoidal <: RLInt end

"""
# Riemann Liouville sense fractional integral using fractional rectangular formula to approximate.

### Example

```julia-repl
julia> fracint(x->x, 0.5, 1, 0.001, RLInt_Rectangular())
```
"""
struct RLInt_Rectangular <: RLInt end

"""
# Riemann Liouville sense fractional integral using cubic spline interpolation to approximate.

### Example

```julia-repl
julia> fracint(x->x, 0.5, 1, 0.0001, RLInt_Cubic_Spline_Interp())
```

Error estimate of this method is ``\\mathcal{O(h^4)}``, it is determined by the error of the cubic spline interpolation.

!!! warning "Set h as 0.001 or bigger"
    For some reason, in the **RLInt_Cubic_Spline_Interp** method, set **h** as 0.001 would get better result.

> Numerical Methods for Fractional Calculus Page 34
"""
struct RLInt_Cubic_Spline_Interp <: RLInt end

####################################
###     Type defination done     ###
####################################

"""
Check if the format of nargins is correct
"""
function checks(f, α, start_point, end_point)
    α % 1 != 0 ? nothing : error("Decimal number only!")
    if isa(end_point, Number)
        start_point < end_point ? nothing : DomainError("Domain error! Start point must smaller than end point")
    elseif isa(end_point, AbstractArray)
        start_point < minimum(end_point) ? nothing : DomainError("Vector domain error! Start point must smaller than end point")
    else
        ErrorException("Please input correct point you wan tto compute")
    end

    #The fractional integral of number is relating with the end_point.
    if isa(f, Number)
        f == 0 ? 0 : 2*f*sqrt(end_point/pi)
    end
end

function fracint(f::Union{Function, Number}, α, start_point, end_point::Real, h, ::RL_Direct)
    #checks(f, α, start_point, end_point)

    temp1 = f(start_point) .* (end_point-start_point) .^α
    #Use Complex differentiation to obtain the differentiation
    g(τ) = imag(f(τ .+ 1*im*h) ./ h) .* (end_point-τ) .^α
    temp2 = quadgk(g, start_point, end_point)
    result = (temp1 .+ temp2) ./gamma(α+1)
    return result
end

function fracint(f::Function, fd::Function, α, start_point, end_point, ::RL_Direct_First_Diff_Known)
    #checks(f, α, start_point, end_point)

    #The fractional integral of number is relating with the end_point.
    if typeof(f) <: Number
        if f == 0
            return 0
        else
            return 2*f*sqrt(end_point/pi)
        end
    end

    temp1 = f(start_point) .* (end_point-start_point) .^α
    g(τ) = fd(τ) .* (end_point-τ) .^α
    temp2 = quadgk(g, start_point, end_point)
    result = (temp1 .+ temp2) ./ gamma(α+1)
    return result
end

function fracint(f::Union{Function, Number}, α::Float64, start_point, end_point::AbstractArray, ::RL_Direct_First_Diff_Known)::Vector
    ResultArray = Float64[]

    for (_, value) in enumerate(end_point)
        append!(ResultArray, fracint(f, α, start_point, value, RL_Direct_First_Diff_Known())[1])
    end
    return ResultArray
end

function fracint(f::Union{Function, Number}, α::Float64, end_point, h, ::RL_Piecewise)::Float64
    #checks(f, α, 0, end_point)

    #Init
    n = Int64(floor(end_point/h))
    result = zero(Float64)

    for i ∈ 0:n
        result += W(i, n, α)*f(i*h)
    end

    result1 = result*h^α/gamma(α+2)
    return result1
end
function W(i, n, α)
    if i == 0
        return n^α*(α+1-n) + (n-1)^(α+1)
    elseif i == n
        return 1
    else
        return (n-i-1)^(α+1) + (n-i+1)^(α+1) - 2*(n-i)^(α+1)
    end
end

function fracint(f::Union{Number, Function}, α::Float64, end_point::AbstractArray, h, ::RL_Piecewise)::Vector
    result = map(x->fracint(f, α, x, h, RL_Piecewise()), end_point)
    return result
end

function fracint(f::Union{Function, Number}, α::Float64, end_point, h, ::RLInt_Approx)::Float64
    #checks(f, α, 0, end_point)

    α = -α
    n = Int64(floor(end_point/h))
    result = zero(Float64)

    for i ∈ 0:n-1
        result += (f(end_point - i*h) + f(end_point - (i+1)*h))/2*((i+1)^(-α) - i^(-α))
    end

    result1 = result*end_point^(-α)*n^α/gamma(1-α)
    return result1
end

function fracint(f::Union{Function, Number}, α::Float64, end_point::AbstractArray, h, ::RLInt_Approx)::Vector
    ResultArray = Float64[]

    for (_, value) in enumerate(end_point)
        append!(ResultArray, fracint(f, α, value, h, RLInt_Approx())[1])
    end
    return ResultArray
end

function fracint(f::Union{Number, Function}, α::Float64, end_point::AbstractArray, h, ::RLInt_Approx)::Vector
    result = map(x->fracint(f, α, x, h, RLInt_Approx()), end_point)
    return result
end

function fracint(f::Union{Function, Number}, α::Float64, end_point::Number, h, ::RL_LinearInterp)::Float64

    α = -α
    n = Int64(floor(end_point/h))
    result = zero(Float64)

    for i ∈ 0:n-1
        result += ((i+1)*f(end_point-i*h)-i*f(end_point-(i+1)*h))/(-α)*((i+1)^(-α) - i^(-α))+(f(end_point - (i+1)*h) - f(end_point-i*h))/(1-α)*((i+1)^(1-α) - i^(1-α))
    end

    result1 = result*end_point^(-α)*n^α/gamma(-α)
    return result1
end

function fracint(f::Union{Number, Function}, α::Float64, end_point::AbstractArray, h, ::RL_LinearInterp)::Vector
    result = map(x->fracint(f, α, x, h, RL_LinearInterp()), end_point)
    return result
end

function fracint(f, α::Number, end_point, h::Float64, ::RLInt_Matrix)
    N = Int64(floor(end_point/h+1))
    tspan = collect(0:h:end_point)
    return J(N, α, h)*f.(tspan)
end

function omega(n, p)
    omega = zeros(n+1)

    omega[1]=1
    @fastmath @inbounds @simd for i ∈ 1:n
        omega[i+1]=(1-(p+1)/i)*omega[i]
    end

    return omega

end
function J(N, p, h::Float64)
    result=zeros(N, N)
    temp=omega(N, -p)

    @fastmath @inbounds @simd for i ∈ 1:N
        result[i, 1:i]=reverse(temp[1:i])
    end

    return h^p*result
end

#=
RLInt_Simpson Algorithm
=#
function fracint(f, α, point, h, ::RLInt_Simpson)
    N = Int64(floor(point/h))
    temp1 = 0
    temp2 = 0

    @fastmath @inbounds @simd for i ∈ 0:N
        temp1 += cₖₙ(i, N, α)*f(i*h)
    end

    @fastmath @inbounds @simd for k ∈ 0:N-1
        temp2 += ĉₖₙ(k, N, α)*f((k+0.5)*h)
    end

    return h^α/gamma(α+3)*temp1 + 4*h^α/gamma(α+3)*temp2
end

function cₖₙ(k, n, α)
    if k == 0
        return 4*((n+1)^(2+α)-n^(2+α))-(α+2)*(3*(n+1)^(1+α)+n^(1+α))+(α+2)*(α+1)*(n+1)^α
    elseif 1 ≤ k ≤ n-1
        return -(α+2)*((n+1-k)^(1+α)+6*(n-k)^(1+α)+(n-k-1)^(1+α))+4*((n+1-k)^(2+α)-(n-1-k)^(2+α))
    elseif k == n
        return 2-α
    end
end
function ĉₖₙ(k, n, α)
    return ((α+2)*((n+1-k)^(1+α)+(n-k)^(1+α))-2*((n+1-k)^(2+α)-(n-k)^(2+α)))
end

function fracint(f::Union{Number, Function}, α::Float64, end_point::AbstractArray, h, ::RLInt_Simpson)::Vector
    result = map(x->fracint(f, α, x, h, RLInt_Simpson()), end_point)
    return result
end

#=
RLInt_Trapezoidal Algorithm
=#
function fracint(f, α, point, h, ::RLInt_Trapezoidal)
    N = Int64(floor(point/h))
    result = zero(Float64)

    @fastmath @inbounds @simd for i ∈ 0:N
        result += aₖₙ(i, N, α)*f(i*h)
    end

    return h^α/gamma(α+2)*result
end

function aₖₙ(k, n, α)
    if k == 0
        return (n-1)^(α+1) - (n-1-α)*n^α
    elseif 1 ≤ k ≤ n-1
        return (n-k+1)^(α+1)+(n-1-k)^(α+1)-2*(n-k)^(α+1)
    else
        return 1
    end
end

function fracint(f, α, point, h, ::RLInt_Rectangular)
    N=Int64(floor(point/h))
    result = zero(Float64)

    @fastmath @inbounds @simd for i ∈ 0:N-1
        result += bₖ(N-i-1, α)*f(i*h)
    end

    return h^α/gamma(α+1)*result
end

function bₖ(k, α)
    return (k+1)^α-k^α
end

function fracint(f::Union{Number, Function}, α::Float64, end_point::AbstractArray, h, ::RLInt_Rectangular)::Vector
    result = map(x->fracint(f, α, x, h, RLInt_Rectangular()), end_point)
    return result
end

#=
RLInt_Cubic_Spline_Interp Algorithm, when h is 0.01 behave best
=#
function fracint(f, α, point, h, ::RLInt_Cubic_Spline_Interp)
    N = Int64(floor(point/h))
    result = zero(Float64)

    @fastmath @inbounds @simd for j ∈ 0:N
        result += eⱼₙ(j, N, α)*f(j*h) + h*êⱼₙ(j, N, α)*first_order(f, j*h, h)
    end

    return h^α/gamma(α+4)*result
end

function eⱼₙ(j, n, α)
    if j == 0
        return -6*(n-1)^(2+α)*(1+2*n+α) + n^α*(12*n^3-6*(3+α)*n^2 + (1+α)*(2+α)*(3+α))
    elseif j == n
        return 6*(1+α)
    else
        return 6*(4*(n-j)^(3+α) + (n-j-1)^(2+α)*(2*j-2*n-1-α) + (1+n-j)^(2+α)*(2*j-2*n+1+α))
    end
end
function êⱼₙ(j, n, α)
    if j == 0
        return -2*(n-1)^(2+α)*(3*n+α) + n^(1+α)*(6*n^2-4*(3+α)*n + (2+α)*(3+α))
    elseif j == n
        return -2*α
    else
        return 2*(n-j-1)^(α+2)*(3*j-3*n-α)-2*(n-j+1)^(α+2)*(3*j-3*n+α) - (n-j)^(α+2)*(24+8α)
    end
end

function fracint(f::Union{Number, Function}, α::Float64, end_point::AbstractArray, h, ::RLInt_Cubic_Spline_Interp)::Vector
    result = map(x->fracint(f, α, x, h, RLInt_Cubic_Spline_Interp()), end_point)
    return result
end
# Deploy Complex Step Differentiation to compute the first order derivative.
function first_order(f, point, h)
    return imag(f(point + im*h))/h
end


text:
module TestBioSequences

using Test
using Random
using StableRNGs
using LinearAlgebra: normalize
import BioSymbols
using BioSequences
using StatsBase
using YAML
using BioGenerics

const codons = [
    "AAA", "AAC", "AAG", "AAU",
    "ACA", "ACC", "ACG", "ACU",
    "AGA", "AGC", "AGG", "AGU",
    "AUA", "AUC", "AUG", "AUU",
    "CAA", "CAC", "CAG", "CAU",
    "CCA", "CCC", "CCG", "CCU",
    "CGA", "CGC", "CGG", "CGU",
    "CUA", "CUC", "CUG", "CUU",
    "GAA", "GAC", "GAG", "GAU",
    "GCA", "GCC", "GCG", "GCU",
    "GGA", "GGC", "GGG", "GGU",
    "GUA", "GUC", "GUG", "GUU",
    "UAA", "UAC", "UAG", "UAU",
    "UCA", "UCC", "UCG", "UCU",
    "UGA", "UGC", "UGG", "UGU",
    "UUA", "UUC", "UUG", "UUU",
    # translatable ambiguities in the standard code
    "CUN", "CCN", "CGN", "ACN",
    "GUN", "GCN", "GGN", "UCN"
]

function random_translatable_rna(n)
    probs = fill(1.0 / length(codons), length(codons))
    cumprobs = cumsum(probs)
    r = rand()
    x = Vector{String}(undef, n)
    for i in 1:n
        x[i] = codons[searchsorted(cumprobs, rand()).start]
    end

    return string(x...)
end

function get_bio_fmt_specimens()
    path = joinpath(dirname(@__FILE__), "BioFmtSpecimens")
    if !isdir(path)
        run(`git clone --depth 1 https://github.com/BioJulia/BioFmtSpecimens.git $(path)`)
    end
end

# The generation of random test cases...

function random_array(n::Integer, elements, probs)
    cumprobs = cumsum(probs)
    x = Vector{eltype(elements)}(undef, n)
    for i in 1:n
        x[i] = elements[searchsorted(cumprobs, rand()).start]
    end
    return x
end

# Return a random DNA/RNA sequence of the given length.
function random_seq(n::Integer, nts, probs)
    cumprobs = cumsum(probs)
    x = Vector{Char}(undef, n)
    for i in 1:n
        x[i] = nts[searchsorted(cumprobs, rand()).start]
    end
    return String(x)
end

function random_seq(::Type{A}, n::Integer) where {A<:Alphabet}
    # TODO: Resolve the use of symbols(A()).
    nts = symbols(A())
    probs = Vector{Float64}(undef, length(nts))
    fill!(probs, 1 / length(nts))
    return LongSequence{A}(random_seq(n, nts, probs))
end

function random_dna(n, probs=[0.24, 0.24, 0.24, 0.24, 0.04])
    return random_seq(n, ['A', 'C', 'G', 'T', 'N'], probs)
end

function random_rna(n, probs=[0.24, 0.24, 0.24, 0.24, 0.04])
    return random_seq(n, ['A', 'C', 'G', 'U', 'N'], probs)
end

function random_aa(len)
    return random_seq(len,
        ['A', 'R', 'N', 'D', 'C', 'Q', 'E', 'G', 'H', 'I',
         'L', 'K', 'M', 'F', 'P', 'S', 'T', 'W', 'Y', 'V', 'X' ],
        push!(fill(0.049, 20), 0.02))
end

function intempdir(fn::Function, parent=tempdir())
    dirname = mktempdir(parent)
    try
        cd(fn, dirname)
    finally
        rm(dirname, recursive=true)
    end
end

function random_dna_kmer(len)
    return random_dna(len, [0.25, 0.25, 0.25, 0.25])
end

function random_rna_kmer(len)
    return random_rna(len, [0.25, 0.25, 0.25, 0.25])
end

function random_dna_kmer_nucleotides(len)
    return random_array(len, [DNA_A, DNA_C, DNA_G, DNA_T],
                        [0.25, 0.25, 0.25, 0.25])
end

function random_rna_kmer_nucleotides(len)
    return random_array(len, [RNA_A, RNA_C, RNA_G, RNA_U],
                        [0.25, 0.25, 0.25, 0.25])
end

function dna_complement(seq::AbstractString)
    seqc = Vector{Char}(undef, length(seq))
    complementer = Dict(zip("-ACGTSWYRKMDVHBN", "-TGCASWRYMKHBDVN"))
    for (i, c) in enumerate(seq)
        seqc[i] = complementer[c]
    end
    return String(seqc)
end

function rna_complement(seq::AbstractString)
    seqc = Vector{Char}(undef, length(seq))
    complementer = Dict(zip("-ACGUSWYRKMDVHBN", "-UGCASWRYMKHBDVN"))
    for (i, c) in enumerate(seq)
        seqc[i] = complementer[c]
    end
    return String(seqc)
end

function random_interval(minstart, maxstop)
    start = rand(minstart:maxstop)
    return start:rand(start:maxstop)
end

include("symbols.jl")

@testset "BioSequences" begin
    a = dna"A-CG-G"; b = rna"A-CG-G"; c = aa"AK-MV-";
    @test ungap(a) == dna"ACGG"
    @test ungap(b) == rna"ACGG"
    @test ungap(c) == aa"AKMV"
    @test ungap!(a) === a && a == dna"ACGG"
    @test ungap!(b) === b && b == rna"ACGG"
    @test ungap!(c) === c && c == aa"AKMV"
end

@testset "LongSequences" begin
    include("longsequences/conversion.jl")
    include("longsequences/basics.jl")
    include("longsequences/seqview.jl")
    include("longsequences/hashing.jl")
    include("longsequences/iteration.jl")
    include("longsequences/subseq.jl")
    include("longsequences/mutability.jl")
    include("longsequences/print.jl")
    include("longsequences/transformations.jl")
    include("longsequences/mutability.jl")
    include("longsequences/predicates.jl")
    include("longsequences/find.jl")
    include("longsequences/randseq.jl")
    include("longsequences/shuffle.jl")
end

@testset "ReferenceSequences" begin
    include("refseq/conversion.jl")
    include("refseq/basics.jl")
    include("refseq/long.jl")
    include("refseq/print.jl")
    include("refseq/random.jl")
end

include("composition.jl")

@testset "Mers" begin
    include("mers/conversion.jl")
    include("mers/comparisons.jl")
    include("mers/length.jl")
    include("mers/arithmetic.jl")
    include("mers/access.jl")
    include("mers/random.jl")
    include("mers/find.jl")
    include("mers/print.jl")
    include("mers/transformations.jl")
    include("mers/mismatches.jl")
    include("mers/debruijn_neighbors.jl")
    include("mers/shuffle.jl")
end

@testset "Iterators" begin
    include("iterators/condition.jl")
    include("iterators/eachmer.jl")
end

@testset "Search" begin
    include("search/exact.jl")
    include("search/approximate.jl")
    include("search/regex.jl")
    include("search/pwm.jl")
end

include("counting.jl")
include("translation.jl")
include("demultiplexer.jl")

end


text:
### CONVENTIONS ###
# quaternion: q = [rho q4] where rho is unit vector
# MRP: stereographic projection for MRP uses domain S^3 \ {(0,0,0,-1)}

# References
# A Survey of Attitude Representations; Shuster 1993
# Spacecraft Dynamics and Control; de Ruiter, Damaren, & Forbes 2013
# Time-Optimal Reorientation of a Spacecraft Using an Inverse Dynamics Optimization Method; Boyarko, Romano and Yakimenko 2011
# A General Construction Scheme for Unit Quaternion Curves with Simple High Order Derivatives; Kim et al. 1995
# Modified Rodrigues Parameters: An Efficient Representation of Orientation in 3D VIsion and Graphics; Terzakis et al. 2018

function skew(v::Vector{T}) where T
  # generates 3x3 skew symmetric matrix for vector
  return [0     -v[3]   v[2];
          v[3]   0      -v[1];
         -v[2]   v[1]   0]
end

function dcm2angle(dcm::Matrix{T}) where T
  # 1.3.2.2 on p21 in Spacecraft Dynamics and Control by Anton de Ruiter
  # Eq. 102 in Shuster
  return acos(0.5*(trace(dcm) - 1))
end

function dcm2vec(dcm::Matrix{T}) where T
  # 1.3.2.2 on p21 in Spacecraft Dynamics and Control by Anton de Ruiter
  # Eq. 103a in Shuster
  phi = dcm2angle(dcm)

  if abs(phi^2) < 0.00001
    # vector undefined if phi == 0, but physically meaningless
    return zeros(T,3)
  elseif abs(phi) == pi
    # three possibilities if phi == pi, pick one
    # Eq. 104 in Shuster
    # Eq. 1.29 in SDAC_deRuiter on p21

    a1 = sqrt(0.5*(dcm[1,1] + 1))
    a2 = sqrt(0.5*(dcm[2,2] + 1))
    a3 = sqrt(0.5*(dcm[3,3] + 1))

    # a1 > 0
    a = [abs(a1); sign(dcm[1,2])*abs(a2); sign(dcm[1,3])*abs(a3)]
    # a2 > 0
    [sign(dcm[1,2])*abs(a1); abs(a2); sign(dcm[2,3])*abs(a3)]
    # a3 > 0
    [sign(dcm[1,3])*abs(a1); sign(dcm[2,3])*abs(a2); abs(a3)]
    return a./norm(a)
  end

  a = zeros(T,3)
  a[1] = 0.5*(dcm[2,3] - dcm[3,2])/sin(phi)
  a[2] = 0.5*(dcm[3,1] - dcm[1,3])/sin(phi)
  a[3] = 0.5*(dcm[1,2] - dcm[2,1])/sin(phi)
  return a
end

function dcm2quat(dcm::Matrix{T}, canonicalize::Bool=false) where T
  # canonicalize: true corresponds to positive scalar component q for dual representation
  # Eq. 163-168 in Shuster

  phi,axis = dcm2angle(dcm),dcm2vec(dcm)

  q = [sin(0.5*phi).*axis; cos(0.5*phi)]
  q = canonicalize ? sign(q[4])*q./norm(q) : q./norm(q)
  return q
end

function dcm2mrp(dcm::Matrix{T}, canonicalize::Bool=false) where T
  return quat2mrp(dcm2quat(dcm, canonicalize))
end

function vec2quat(axis::Vector{T}, phi::T, canonicalize::Bool=false) where T
  q = [sin(0.5*phi)*axis/norm(axis); cos(0.5*phi)]
  q = canonicalize ? sign(q[4])*q./norm(q) : q./norm(q)
  return q
end

function quat2angle(q::Vector{T}) where T
  return 2*atan2(norm(q[1:3]), q[4])
end

function quat2vec(q::Vector{T}, canonicalize::Bool=false) where T
  # canonicalize: true corresponds to positive scalar component q for dual representation
  if norm(q[1:3])^2 < 0.00001
    return zeros(T,3)
  end

  a = canonicalize ? sign(q[4])*q[1:3]/norm(q[1:3]) : q[1:3]/norm(q[1:3])
  return a
end

function quat2dcm(q::Vector{T}, attitude_convention::String="orientation") where T
  # Eq. 97 in Shuster (note sin term sign flipped based on skew matrix convention)
  # Eq. 1.26 on p20 in Spacecraft Dynamics and Control by Anton de Ruiter
  phi = quat2angle(q)
  n = quat2vec(q)
  dcm = cos(phi)*eye(3) + (1-cos(phi))*n*n' - sin(phi)*skew(n)

  if attitude_convention == "rotation"
    dcm = dcm'
  elseif attitude_convention != "orientation"
    @warn "Must specify if attitude is rotation or orientation! Defaulting to orientation."
  end
  return dcm
end

function vec2dcm(e::Vector{T}, phi::T) where T
  q = vec2quat(e, phi)
  return quat2dcm(q)
end

function quat_derivative(q::Vector{T}, w::Vector{T}) where T
  # Eq. 305 in Shuster
  omega_skew = [-skew(w) w; -w' 0]
  return 0.5*omega_skew*q
end

function quat_multiply(q1::Vector{T}, q2::Vector{T}) where T
  # Eq. 169-171 in Shuster
  # q1 * q2 <==> R(q1) * R(q2)

  q_product = zeros(T,4)
  q_product[1:3] = q1[4]*q2[1:3] + q2[4]*q1[1:3] - cross(Vector(q1[1:3]), Vector(q2[1:3]))
  q_product[4] = q1[4]*q2[4] - dot(q1[1:3], q2[1:3])
  return q_product
end

function quat_inv(q::Vector{T}) where T
  # Eq. 177 in Shuster
  qinv = copy(q)
  qinv[1:3]*=-1
  return qinv
end

function quat_error(q1::Vector{T}, q2::Vector{T}) where T
  # corresponds to  q2 = quat_multiply(dq,q1)
  # quat_multiply(q2, quat_inv(q1))

  # Spacecraft Dynamics and Control - Sidi
  #prod_matrix = [q1[4] q1[3] -q1[2] q1[1]; -q1[3] q1[4] q1[1] q1[2]; q1[2] -q1[1] q1[4] q1[3]; -q1[1] -q1[2] -q1[3] q1[4]]
  #dq = prod_matrix * [-q2[1:3]; q2[4]]

  # Space Vehicle Dynamics and Control - Wie
  prod_matrix = [q1[4] q1[3] -q1[2] -q1[1]; -q1[3] q1[4] q1[1] -q1[2]; q1[2] -q1[1] q1[4] -q1[3]; q1[1] q1[2] q1[3] q1[4]]
  dq = prod_matrix * q2
  return dq
end

function quat_rotate(q::Vector{T}, v::Vector{T}, attitude_convention::String="orientation") where T
  # Eq. 183 in Shuster
  # v: 3 element vector to be rotated

  qr = q
  if attitude_convention == "rotation"
    qr = quat_inv(q)
  elseif attitude_convention != "orientation"
    @warn "Must specify if quaternion represents rotation or orientation! Defaulting to orientation."
  end

  v_bar = [v; 0]
  return quat_multiply(qr, quat_multiply(v_bar, quat_inv(qr)))[1:3]
end

function quat_interp(q0::Vector{T}, q1::Vector{T}, t::T) where T
  # p371 in "A General Construction Scheme for Unit Quaternion Curves with Simple High Order Derivatives"
  # NOTE: q1 * q2 <==> R(q1) * R(q2)

  if (t>1 || t<0)
    print("Supplied time t must lie in [0,1]!\n")
    return t<0 ? q0 : q1
  end

  qerr = quat_multiply(quat_inv(q0),q1)
  return quat_multiply(q0, quat_exp(t*quat_log(qerr)))
end

function quat_exp(v::Vector{T}) where T
  # gives quaternion q corresponding to orientation obtained
  # from an initial orientation [0 0 0 1] by rotating of an angle norm(v)
  # around fixed direction v/norm(v)
  # Eq. 9 in "Time-Optimal Reorientation of a Spacecraft Using an Inverse Dynamics Optimization Method"
  phi = norm(v)
  if phi^2 < 0.00001
    return Vector{T}([0;0;0;1])
  end
  return Vector{T}([sin(0.5*phi)*v/phi;cos(0.5*phi)])
end

function quat_log(q::Vector{T}) where T
  # gives vector having direction of the Euler's axis and magnitude equal
  # to Euler's angle of the orientation from [0 0 0 1] to q
  # Eq. 9 in "Time-Optimal Reorientation of a Spacecraft Using an Inverse Dynamics Optimization Method"

  if norm(q[1:3])^2 < 0.00001
    return zeros(T,3)
  end
  return quat2angle(q)*quat2vec(q)
end

function quat2mrp(q::Vector{T}) where T
  # Eq. 249 in Shuster
  # Eq. 22 in Terzakis et al.
  if abs(q[4]+1) < 0.000001
    @warn "Stereographic projection undefined for (0,0,0,-1)"
    return -q[1:3]/(1-q[4])
  end
  return q[1:3]/(1+q[4])
end

function mrp2quat(p::Vector{T}) where T
  # Eq. 25 & 26 in Terzakis et al.
  s = (1-norm(p)^2)/(1+norm(p)^2)
  v = 2*p/(1+norm(p)^2)
  return [v;s]
end

function mrp2dcm(p::Vector{T}, attitude_convention::String="orientation") where T
  # Eq. 255b in Shuster
  # NOTE (6/26): switched +4 term to -4 term because that's what I
  # found everywhere else and unit tests pass with this only

  dcm = eye(3) + (8*skew(p)^2 - 4*(1 - norm(p)^2)*skew(p))/(1+norm(p)^2)^2

  if attitude_convention == "rotation"
    dcm = dcm'
  elseif attitude_convention != "orientation"
    @warn "Must specify if attitude is rotation or orientation! Defaulting to orientation."
  end
  return dcm
end

function mrp_derivative(p::Vector{T}, w::Vector{T}) where T
  # Eq. 338 in Shuster (in terms of the body-referenced angular velocity)
  return 0.25*((1-norm(p)^2)*w - 2*cross(w,p) + 2*dot(w,p)*p)
end

function mrp_multiply(p1::Vector{T}, p2::Vector{T}) where T
  # Eq. 257 in Shuster
  # p1 * p2 <==> R(p1) * R(p2)
  # TODO(acauligi): check if denominator is close to 0. for edge case

  return ((1-norm(p2)^2)*p1 + (1-norm(p1)^2)*p2 - 2*cross(p1,p2)) /
    (1+norm(p1)^2*norm(p2)^2 - 2*dot(p1,p2))
end

# function Aq(J::Matrix{T},X::Vector{T})
#   Jxx,Jyy,Jzz = diag(J)
#   qx,qy,qz,qw = X[1:4]
#   wx,wy,wz = X[5:7]
#   Jq = 0.5*[0 wz -wy wx qw -qz qy;
#           -wz 0 wx wy qz qw -qx;
#            wy -wx 0 wz -qy qx qw;
#           -wx -wy -wz 0 -qx -qy -qz]
#   Jw = [0 0 0 0 0                     (Jyy-Jzz)*wz/Jxx   (Jyy-Jzz)*wy/Jxx;
#         0 0 0 0 -(Jxx-Jzz)*wz/Jyy  0                     -(Jxx-Jzz)*wx/Jyy;
#         0 0 0 0 (Jxx-Jyy)*wy/Jzz   (Jxx-Jyy)*wx/Jzz                     0]
#   return [Jq;Jw]
# end

# function Apb(J::Matrix{T},X::Vector{T})
#   Jxx,Jyy,Jzz = diag(J)
#   px,py,pz = X[1:3]
#   wx,wy,wz = X[4:6]
#   Jpb = [(px*wx)/2+(py*wy)/2+(pz*wz)/2 wz/2+(px*wy)/2-(py*wx)/2 (px*wz)/2-wy/2-(pz*wx)/2 px^2/4-py^2/4-pz^2/4+1/4 (px*py)/2-pz/2 py/2+(px*pz)/2;
#         (py*wx)/2-(px*wy)/2-wz/2 (px*wx)/2+(py*wy)/2+(pz*wz)/2 wx/2+(py*wz)/2-(pz*wy)/2 pz/2+(px*py)/2 -px^2/4+py^2/4-pz^2/4+1/4 (py*pz)/2-px/2;
#         wy/2-(px*wz)/2+(pz*wx)/2 (pz*wy)/2-(py*wz)/2-wx/2 (px*wx)/2+(py*wy)/2+(pz*wz)/2 (px*pz)/2-py/2 px/2+(py*pz)/2 -px^2/4-py^2/4+pz^2/4+1/4]
#   Jw = [0 0 0 0 (Jyy*wz-Jzz*wz)/Jxx (Jyy*wy-Jzz*wy)/Jxx;
#         0 0 0 -(Jxx*wz-Jzz*wz)/Jyy 0 -(Jxx*wx-Jzz*wx)/Jyy;
#         0 0 0 (Jxx*wy-Jyy*wy)/Jzz (Jxx*wx-Jyy*wx)/Jzz 0]
#   return [Jpb;Jw]
# end

# function Api(J::Matrix{T},X::Vector{T})
#   Jxx,Jyy,Jzz = diag(J)
#   px,py,pz = X[1:3]
#   wx,wy,wz = X[4:6]
#   Jpi = [(px*wx)/2+(py*wy)/2+(pz*wz)/2 (px*wy)/2-wz/2-(py*wx)/2 wy/2+(px*wz)/2-(pz*wx)/2 px^2/4-py^2/4-pz^2/4+1/4 pz/2+(px*py)/2 (px*pz)/2-py/2;
#           wz/2-(px*wy)/2+(py*wx)/2 (px*wx)/2+(py*wy)/2+(pz*wz)/2 (py*wz)/2-wx/2-(pz*wy)/2 (px*py)/2-pz/2 -px^2/4+py^2/4-pz^2/4+1/4 px/2+(py*pz)/2;
#           (pz*wx)/2-(px*wz)/2-wy/2 wx/2-(py*wz)/2+(pz*wy)/2 (px*wx)/2+(py*wy)/2+(pz*wz)/2 py/2+(px*pz)/2 (py*pz)/2-px/2 -px^2/4-py^2/4+pz^2/4+1/4]
#   Jw = [0 0 0 0                     (Jyy*wz-Jzz*wz)/Jxx   (Jyy*wy-Jzz*wy)/Jxx;
#         0 0 0 -(Jxx*wz-Jzz*wz)/Jyy  0                     -(Jxx*wx-Jzz*wx)/Jyy;
#         0 0 0 (Jxx*wy-Jyy*wy)/Jzz   (Jxx*wx-Jyy*wx)/Jzz   0]
#   return [Jpi;Jw]
# end

# function B(J::Matrix{T},X::Vector{T})
#   Jxx,Jyy,Jzz = diag(J)
#   Bw = length(X) == 6 ? zeros(T,6,3) : zeros(T,7,3)
#   Bw[end-2:end,:] = diagm([1/Jxx,1/Jyy,1/Jzz])
#   return Bw
# end

# function mrp_jacobian(p,w)
#   # w in inertial frame
#   Ji = zeros(eltype(p), (3,6))
#   p1,p2,p3 = p
#   wi1,wi2,wi3 = w
#
#   # Jb[1,:] = [(p1*wb1)/2+(p2*wb2)/2+(p3*wb3)/2, wb3/2+(p1*wb2)/2-(p2*wb1)/2, (p1*wb3)/2-wb2/2-(p3*wb1)/2, p1^2/4-p2^2/4-p3^2/4+1/4, (p1*p2)/2-p3/2, p2/2 + (p1*p3)/2]
#   # Jb[2,:] = [(p2*wb1)/2-(p1*wb2)/2-wb3/2, (p1*wb1)/2+(p2*wb2)/2+(p3*wb3)/2, wb1/2+(p2*wb3)/2-(p3*wb2)/2, p3/2+(p1*p2)/2, -p1^2/4+p2^2/4-p3^2/4+1/4, (p2*p3)/2 - p1/2]
#   # Jb[3,:] = [wb2/2-(p1*wb3)/2+(p3*wb1)/2, (p3*wb2)/2-(p2*wb3)/2-wb1/2, (p1*wb1)/2+(p2*wb2)/2+(p3*wb3)/2, (p1*p3)/2-p2/2, p1/2+(p2*p3)/2, -p1^2/4-p2^2/4+p3^2/4+1/4]
#
#   Ji[1,:] = [(p1*wi1)/2+(p2*wi2)/2+(p3*wi3)/2, (p1*wi2)/2-wi3/2-(p2*wi1)/2, wi2/2+(p1*wi3)/2-(p3*wi1)/2, p1^2/4-p2^2/4-p3^2/4+1/4, p3/2+(p1*p2)/2, (p1*p3)/2-p2/2]
#   Ji[2,:] = [wi3/2-(p1*wi2)/2+(p2*wi1)/2, (p1*wi1)/2+(p2*wi2)/2+(p3*wi3)/2, (p2*wi3)/2-wi1/2-(p3*wi2)/2, (p1*p2)/2-p3/2, -p1^2/4+p2^2/4-p3^2/4+1/4, p1/2+(p2*p3)/2]
#   Ji[3,:] = [(p3*wi1)/2-(p1*wi3)/2-wi2/2, wi1/2-(p2*wi3)/2+(p3*wi2)/2, (p1*wi1)/2+(p2*wi2)/2+(p3*wi3)/2, p2/2+(p1*p3)/2, (p2*p3)/2-p1/2, -p1^2/4-p2^2/4+p3^2/4+1/4]
#   return Ji
# end

# function quat2rpy(q::Vector{T})
#   sinr = 2*(q[4]*q[1] + q[2]*q[3])
#   cosr = 1 - 2*(q[1]*q[1] + q[2]*q[2])
#   roll = atan2(sinr, cosr)
#
#   sinp = 2*(q[4]*q[2] - q[3]*q[1])
#   if abs(sinp) >= 1
#     pitch = sign(sinp)*pi/2
#   else
#     pitch = asin(sinp)
#   end
#
#   siny = 2*(q[4]*q[3] - q[1]*q[2])
#   cosy = 1 - 2*(q[2]*q[2] + q[3]*q[3])
#   yaw = atan2(siny, cosy)
#
#   return [roll; pitch; yaw]
# end

# function plot_attitude_sphere(q_matrix::Array, plot_live::Bool=false, meshing::Int=20)
#   if size(q_matrix)[1] == 4
#     q_matrix = q_matrix'
#   end
#
#   n_samples = size(q_matrix)[1]
#
#   vector_matrix = zeros(n_samples, 3)
#   [vector_matrix[idx,:] = quat2angle(q_matrix[idx,:]) * quat2vec(q_matrix[idx,:], false) for idx in range(1, n_samples)]
#
#
#   # sphere of radius pi
#   phi = linspace(0, 2*pi, meshing)
#   theta = linspace(0, pi, meshing)
#   x = pi*cos(phi)*sin(theta)'
#   y = pi*sin(phi)*sin(theta)'
#   z = pi*ones(meshing)*cos(theta)'
#
#   figure()
#   plot_wireframe(x,y,z, linestyle="-.")
#
#   for idx in 1:n_samples
#     scatter3D(vector_matrix[idx, 1], vector_matrix[idx, 2], vector_matrix[idx, 3], color="red")
#     plot_live && pause(1/n_samples)
#   end
# end

# function plot_vector_rotation(Qs::Matrix{T}, v::Matrix{T}=eye(T,3)) where T
#   # Qs: 4xN matrix of quaternions
#   # v: 3xN vectors to be plotted; defaults to body frame
#
#   n_quat = size(Qs,2)
#   n_vec = size(v,2)
#   slew_matrix = zeros(T,3,n_quat,n_vec)
#
#   for i in 1:n_quat
#     q = Qs[:,i]
#     for j in 1:n_vec
#       slew_matrix[:,i,j] = quat_rotate(q,v[:,j],"rotation")
#     end
#   end
#
#   # unit sphere
#   meshing = 25
#   phi = linspace(0,2*pi,meshing)
#   theta = linspace(0,pi,meshing)
#   x = cos.(phi)*sin.(theta)'
#   y = sin.(phi)*sin.(theta)'
#   z = ones(meshing)*cos.(theta)'
#
#   cs = ["red","green","blue"]
#
#   fig = PyPlot.figure()
#   for i in 1:n_quat
#     clf()
#     PyPlot.plot_wireframe(x,y,z, linestyle="-.", color="blue")
#     for j in 1:n_vec
#       v = slew_matrix[:,i,j]
#       PyPlot.plot3D([0; v[1]], [0; v[2]], [0, v[3]], color=cs[j%3+1])
#     end
#     pause(0.0000001)
#   end
#
#   return fig, slew_matrix
# end

# function plot_q_slew(q_matrix::Array, v::Array=[-1])
#   # q_matrix: Nx4 matrix
#   # v: 3D vector to be rotated by DCM
#
#   # 4xN to Nx4
#   if size(q_matrix) == (4,4)
#     # check case of 4x4 matrix and transpose if quaternions are column wise
#     if size(q_matrix)[1] == 4 && (norm(q_matrix[1,:]) - 1 > norm(q_matrix[:,1]) - 1)
#       # NOTE: cannot identify if matrix is identity matrix i.e. eye(4)
#       q_matrix = q_matrix'
#     end
#   elseif size(q_matrix)[1] == 4
#     q_matrix = q_matrix'
#   end
#
#   numel = size(q_matrix)[1]
#
#   # basis vectors to slew
#   v1 = [1;0;0]
#   v2 = [0;1;0]
#   v3 = [0;0;1]
#
#   # generate slew matrix for basis vectors
#   slew_matrix = zeros(numel,3,3)
#   for row_idx in 1:numel
#     slew_matrix[row_idx,:, 1] =  quat_rotate(q_matrix[row_idx,:], v1, "rotation")
#     slew_matrix[row_idx,:, 2] =  quat_rotate(q_matrix[row_idx,:], v2, "rotation")
#     slew_matrix[row_idx,:, 3] =  quat_rotate(q_matrix[row_idx,:], v3, "rotation")
#   end
#
#   vector_matrix = zeros(numel,3)
#   if v != [-1]
#     v = v/norm(v)
#     for row_idx in 1:numel
#       vector_matrix[row_idx,:] =  quat_rotate(q_matrix[row_idx,:], v, "rotation")
#     end
#   end
#
#   # unit sphere
#   mag = norm(v)
#   meshing = 50
#   phi = linspace(0, 2*pi, meshing)
#   theta = linspace(0, pi, meshing)
#   x = mag*cos(phi)*sin(theta)'
#   y = mag*sin(phi)*sin(theta)'
#   z = mag*ones(meshing)*cos(theta)'
#
#   fig1 = figure()
#   for idx in 1:size(slew_matrix)[1]
#     clf()
#     plot_wireframe(x,y,z, linestyle="-.", color="white")
#
#     v_here = slew_matrix[idx,:,1]
#     plot3D([0; v_here[1]], [0; v_here[2]], [0, v_here[3]], color="red")
#
#     v_here = slew_matrix[idx,:,2]
#     plot3D([0; v_here[1]], [0; v_here[2]], [0, v_here[3]], color="blue")
#
#     v_here = slew_matrix[idx,:,3]
#     plot3D([0; v_here[1]], [0; v_here[2]], [0, v_here[3]], color="green")
#     pause(0.00001)
#
#
#     if v != [-1]
#       plot3D(vector_matrix[1:idx,1], vector_matrix[1:idx,2], vector_matrix[1:idx,3], color="black")
#     end
#   end
# end


text:
using CSV
using Plots
using DataFrames
using PlutoUI

function read_df(path)
    d = DataFrame(
                  N = Int[],
                  Z = Int[],
                  A = Int[],
                  binding_energy = Float32[],
                  binding_energy_per_nucleon = Float32[],
                  )
    row = 0
    start = 40
    a = -1
    for line in eachline(open(path))
        row += 1
        if row < start
            continue
        end
        n = parse(Int, line[7:9])
        z = parse(Int, line[12:14])
        # a can be empty (It is only there if it differs from the last a)
        # It looks like 1 - - - 2 - - - 3 ... (- mark lines)
        a_ = tryparse(Int, line[17:19])
        a = isnothing(a_) ? a : a_
        #a = tryparse(Int, line[17:19])
        # sometimes # is used instead of .
        # for this purpose we consider them equal
        b_e = parse(Float32, replace(line[49:60], "#" => "."))
        push!(d, [n, z, a, b_e, b_e/a])
    end
    d
end

function max_binding_energy(df, col)
    gdf = groupby(df, col)
    combine(gdf, :binding_energy_per_nucleon => maximum)
end

df = read_df("data")
max_bindings = max_binding_energy(df, "A")
A = max_bindings.A
E = max_bindings.binding_energy_per_nucleon_maximum
highlights = [1,4,8,12,16,56]
text_y = E[highlights]
y_ = 7800
upper = text_y .> y_
under = text_y .< y_
text_y[upper] = text_y[upper] .- 600
text_y[under] = text_y[under] .+ 500
p = plot(
         A,
         E,
         xscale=:log,
         label="",
         xlabel="A = Z + N",
         ylabel="Binding Energy / A [kEV]",
         )
scatter!(
         A[highlights],
         E[highlights],
         label="",
         markersize=6,
         markershape=:star4,
         )
annotate!(
          highlights,
          text_y,
          # these have to be filled manually because I didnt bother reading them
          ["¹H", "⁴HE", "⁸Be", "¹²C  ", "  ¹⁶O", "⁵⁶FE"],
          )

savefig(p, "binding_energy.png")
savefig(p, "binding_energy.svg")


text:
module DeferredFutures

export DeferredFuture

using AutoHashEquals

@auto_hash_equals immutable DeferredFuture <: Base.AbstractRemoteRef
    outer::Future
end

DeferredFuture() = DeferredFuture(Future())
DeferredFuture(pid::Integer) = DeferredFuture(Future(pid))

function Base.put!(df::DeferredFuture, val)
    rc = RemoteChannel()
    put!(rc, val)
    put!(df.outer, rc)
    return df
end

function Base.isready(df::DeferredFuture)
    isready(df.outer) && isready(fetch(df.outer))
end

function Base.fetch(df::DeferredFuture)
    fetch(fetch(df.outer))
end

function Base.wait(df::DeferredFuture)
    wait(df.outer)
    wait(fetch(df.outer))
    return df
end

end # module


text:
using ProgrammableAPI
using Test


text:
using MorePolynomials
using Intervals
using Test

@testset "Global.jl" begin
    x = LGRPoly([1.0,2,3])
    y = GlobalPoly(x,0,1)

    @test_throws ErrorException push!(y,x,1,2)
    @test_throws ErrorException push!(y,x,0,1)
    @test_throws ErrorException push!(y,x,-1,0)
    @test_nowarn push!(y,x,2,3)

    x = LGRPoly([1.0,2,1])
    y = GlobalPoly(x,0,1)
    @test_nowarn update!(y,[1.0,1,1.0])
    @test_nowarn push!(y,x,2,3)
    @test_nowarn push!(y,x,1,2)
    @test_nowarn push!(y,x,3,4)
    @test_nowarn push!(y,x,-1,0)
    @test_throws DomainError y(8)
    x1 = LagrangePoly([-1,0,1],[0.0,1,0])
    x2 = LagrangePoly([-1,0,1],[0.0,2,0])
    x3 = LagrangePoly([-1,0,1],[0.0,3,0])
    y1 = GlobalPoly(x1,0,1)
    push!(y1,x2,1,2)
    push!(y1,x3,2,3)
    @test_throws DimensionMismatch update!(y1,[1.0,1,1])
    @test_nowarn update!(y1,[0.0,1,2,3,4,5,6])

    x1 = LGRPoly([0.0,1,0])
    x2 = LGRPoly([0.0,2,0])
    x3 = LGRPoly([0.0,3,0])
    y1 = GlobalPoly(x1,0,1)
    push!(y1,x2,2,4)
    push!(y1,x3,4,5)
    @test_nowarn update!(y1,[0.0,1,2,4,5,6,7,8])

    # generate psudorandom values for this
    checkinfsTestFunc(a,b,c,d) = MorePolynomials.checkinfs(Interval(c,d),Interval(a,b)) # a b local -> c d global
    @test checkinfsTestFunc(-Inf,Inf,-Inf,Inf) == true
    @test_throws DomainError checkinfsTestFunc(-Inf,-2,-Inf,Inf)
    @test checkinfsTestFunc(-2,2,-Inf,Inf) == true
    @test checkinfsTestFunc(-Inf,Inf,-Inf,2) == true
    @test checkinfsTestFunc(-Inf,Inf,-2,Inf) == true
    @test checkinfsTestFunc(-Inf,Inf,-2,2) == true
    @test checkinfsTestFunc(-Inf,2,-Inf,1) == true
    @test checkinfsTestFunc(-Inf,2,-Inf,3) == true
    @test_throws DomainError checkinfsTestFunc(-Inf,2,2,Inf)
    @test_throws DomainError checkinfsTestFunc(-2, Inf,-Inf,2)
    @test checkinfsTestFunc(2,Inf,3,Inf) == true
    @test checkinfsTestFunc(-2,2,3,Inf) == true
    @test checkinfsTestFunc(3,4,3,Inf) == true
    @test checkinfsTestFunc(3,4,-Inf,2) == true
    @test checkinfsTestFunc(3,Inf,3,4) == true
    @test_throws DomainError checkinfsTestFunc(3,Inf,2,4)
    @test_throws DomainError checkinfsTestFunc(3,Inf,1,2)
    @test checkinfsTestFunc(3,4,5,6) == true

    bignumber = 999999999999999999999999999999999
    global2localTestFunc(a,b,c,d,x) = MorePolynomials.global2local(Interval(a,b),Interval(c,d),x)# a b local -> c d global
    @test global2localTestFunc(-Inf,Inf,-Inf,Inf,4) == 4
    @test global2localTestFunc(-2,2,-Inf,Inf,Inf) == 2
    @test global2localTestFunc(-2,2,-Inf,Inf,-Inf) == -2
    @test global2localTestFunc(-2,2,-Inf,Inf,0) == 0
    @test global2localTestFunc(-Inf,Inf,-Inf,2,1) == 1
    @test global2localTestFunc(-Inf,Inf,-2,Inf,1) == 1
    @test global2localTestFunc(-Inf,Inf,-2,2,1) == 1
    @test global2localTestFunc(2,Inf,3,Inf,3) == 2
    @test global2localTestFunc(2,Inf,3,Inf,Inf) == Inf
    @test global2localTestFunc(-Inf,2,-Inf,3,3) == 2
    @test global2localTestFunc(-Inf,2,-Inf,3,-Inf) == -Inf
    @test global2localTestFunc(3,4,-Inf,2,-bignumber) ≈ 3
    @test global2localTestFunc(3,4,2,Inf,2) == 3
    @test global2localTestFunc(3,4,2,Inf,bignumber) ≈ 4
    @test global2localTestFunc(3,4,-Inf,2,2) == 4
    @test global2localTestFunc(3,Inf,3,4,3) == 3
    @test global2localTestFunc(3,Inf,3,4,4) == 4
    @test global2localTestFunc(3,4,5,6,5) == 3
    @test global2localTestFunc(3,4,5,6,6) == 4
    @test global2localTestFunc(3,4,5,6,5.5) == 3.5

end


text:
using MPSGE
using Test
using MPSGE.JuMP.Containers

@testset "MPSGE" begin

    @testset "TWOBYTWO (functional version)" begin
        m = Model()

        inputcoeff = add!(m, Parameter(:inputcoeff, value=2.))
        endow = add!(m, Parameter(:endow, value=2.))
        elascoeff = add!(m, Parameter(:elascoeff, value=2.))
        outputmult = add!(m, Parameter(:outputmult, value=2.))

        X = add!(m, Sector(:X))
        Y = add!(m, Sector(:Y))
        U = add!(m, Sector(:U))

        PX = add!(m, Commodity(:PX))
        PY = add!(m, Commodity(:PY))
        PU = add!(m, Commodity(:PU))
        PL = add!(m, Commodity(:PL))
        PK = add!(m, Commodity(:PK))

        RA = add!(m, Consumer(:RA, benchmark=150.))

        add!(m, Production(X, 1, [Output(PX, 100)], [Input(PL, :(25 * $inputcoeff)), Input(PK, 50)]))
        add!(m, Production(Y, :(0.5 * $elascoeff), [Output(PY, 50)], [Input(PL, 20), Input(PK, 30)]))
        add!(m, Production(U, 1, [Output(PU, :(75 * $outputmult))], [Input(PX, 100), Input(PY, 50)]))

        add!(m, DemandFunction(RA, [Demand(PU,150)], [Endowment(PL, :(35 * $endow)), Endowment(PK, 80)]))

        avm = algebraic_version(m)
        @test typeof(avm) == MPSGE.AlgebraicWrapper

        solve!(m)

        @test value(m, :X) ≈ 1.
        @test MPSGE.Complementarity.result_value(m._jump_model[:Y]) ≈ 1.
        @test MPSGE.Complementarity.result_value(m._jump_model[:U]) ≈ 1.
        @test MPSGE.Complementarity.result_value(m._jump_model[:RA]) ≈ 150.

        @test MPSGE.Complementarity.result_value(m._jump_model[:PX]) ≈ 1.
        @test MPSGE.Complementarity.result_value(m._jump_model[:PY]) ≈ 1.
        @test MPSGE.Complementarity.result_value(m._jump_model[:PU]) ≈ 1.
        @test MPSGE.Complementarity.result_value(m._jump_model[:PL]) ≈ 1.
        @test MPSGE.Complementarity.result_value(m._jump_model[:PK]) ≈ 1.

        @test MPSGE.Complementarity.result_value(m._jump_model[Symbol("PL†X")]) ≈ 50.
        @test MPSGE.Complementarity.result_value(m._jump_model[Symbol("PL†Y")]) ≈ 20.
        @test MPSGE.Complementarity.result_value(m._jump_model[Symbol("PK†X")]) ≈ 50.
        @test MPSGE.Complementarity.result_value(m._jump_model[Symbol("PK†Y")]) ≈ 30.
        @test MPSGE.Complementarity.result_value(m._jump_model[Symbol("PX†U")]) ≈ 100.
        @test MPSGE.Complementarity.result_value(m._jump_model[Symbol("PY†U")]) ≈ 50.

        avm2 = algebraic_version(m)
        @test typeof(avm2) == MPSGE.AlgebraicWrapper

        # For now just run these functions, we might add tests for the results
        # at a later point
        repr(MIME("text/plain"), m)
        repr(MIME("text/plain"), avm2)
        repr(MIME("text/latex"), avm2)

        set_fixed!(PX, true)
        set_value(endow, 2.2)
        solve!(m)

        @test value(m, :X) ≈ 1.04880885
        @test MPSGE.Complementarity.result_value(m._jump_model[:Y]) ≈ 1.03886012
        @test MPSGE.Complementarity.result_value(m._jump_model[:U]) ≈ 1.04548206
        @test MPSGE.Complementarity.result_value(m._jump_model[:RA]) ≈ 157.321327225523
        @test MPSGE.Complementarity.result_value(m._jump_model[:PX]) ≈ 1.0000000000
        @test MPSGE.Complementarity.result_value(m._jump_model[:PY]) ≈ 1.00957658
        @test MPSGE.Complementarity.result_value(m._jump_model[:PU]) ≈ 1.00318206
        @test MPSGE.Complementarity.result_value(m._jump_model[:PL]) ≈ 0.95346259
        @test MPSGE.Complementarity.result_value(m._jump_model[:PK]) ≈ 1.04880885
        @test MPSGE.Complementarity.result_value(m._jump_model[Symbol("PL†X")]) ≈ 52.4404424085075
        @test MPSGE.Complementarity.result_value(m._jump_model[Symbol("PL†Y")]) ≈ 21.1770570584356
        @test MPSGE.Complementarity.result_value(m._jump_model[Symbol("PK†X")]) ≈ 47.6731294622795
        @test MPSGE.Complementarity.result_value(m._jump_model[Symbol("PK†Y")]) ≈ 28.877805079685
        @test MPSGE.Complementarity.result_value(m._jump_model[Symbol("PX†U")]) ≈ 100.318205802571
        @test MPSGE.Complementarity.result_value(m._jump_model[Symbol("PY†U")]) ≈ 49.6833066029729
    end

    @testset "TWOBYTWO (macro version)" begin
        m = Model()

        inputcoeff = @parameter(m, inputcoeff, 2.)
        endow = @parameter(m, endow, 2.)
        elascoeff = @parameter(m, elascoeff, 2.)
        outputmult = @parameter(m, outputmult, 2.)

        @sector(m, X)
        @sector(m, Y)
        @sector(m, U)

        @commodity(m, PX)
        @commodity(m, PY)
        @commodity(m, PU)
        @commodity(m, PL)
        @commodity(m, PK)

        @consumer(m, RA, benchmark=150.)

        @production(m, X, 1, [Output(PX, 100)], [Input(PL, :(25 * $inputcoeff)), Input(PK, 50)])
        @production(m, Y, :(0.5 * $elascoeff), [Output(PY, 50)], [Input(PL, 20), Input(PK, 30)])
        @production(m, U, 1, [Output(PU, :(75 * $outputmult))], [Input(PX, 100), Input(PY, 50)])

        @demand(m, RA, [Demand(PU, 150)], [Endowment(PL, :(35 * $endow)), Endowment(PK, 80)])

        avm = algebraic_version(m)
        @test typeof(avm) == MPSGE.AlgebraicWrapper

        solve!(m)

        @test value(m, :X) ≈ 1.
        @test MPSGE.Complementarity.result_value(m._jump_model[:Y]) ≈ 1.
        @test MPSGE.Complementarity.result_value(m._jump_model[:U]) ≈ 1.
        @test MPSGE.Complementarity.result_value(m._jump_model[:RA]) ≈ 150.

        @test MPSGE.Complementarity.result_value(m._jump_model[:PX]) ≈ 1.
        @test MPSGE.Complementarity.result_value(m._jump_model[:PY]) ≈ 1.
        @test MPSGE.Complementarity.result_value(m._jump_model[:PU]) ≈ 1.
        @test MPSGE.Complementarity.result_value(m._jump_model[:PL]) ≈ 1.
        @test MPSGE.Complementarity.result_value(m._jump_model[:PK]) ≈ 1.

        @test MPSGE.Complementarity.result_value(m._jump_model[Symbol("PL†X")]) ≈ 50.
        @test MPSGE.Complementarity.result_value(m._jump_model[Symbol("PL†Y")]) ≈ 20.
        @test MPSGE.Complementarity.result_value(m._jump_model[Symbol("PK†X")]) ≈ 50.
        @test MPSGE.Complementarity.result_value(m._jump_model[Symbol("PK†Y")]) ≈ 30.
        @test MPSGE.Complementarity.result_value(m._jump_model[Symbol("PX†U")]) ≈ 100.
        @test MPSGE.Complementarity.result_value(m._jump_model[Symbol("PY†U")]) ≈ 50.

        avm2 = algebraic_version(m)
        @test typeof(avm2) == MPSGE.AlgebraicWrapper

        # For now just run these functions, we might add tests for the results
        # at a later point
        repr(MIME("text/plain"), m)
        repr(MIME("text/plain"), avm2)
        repr(MIME("text/latex"), avm2)

        set_fixed!(PX, true)
        set_value(endow, 2.2)
        solve!(m)

        @test value(m, :X) ≈ 1.04880885
        @test MPSGE.Complementarity.result_value(m._jump_model[:Y]) ≈ 1.03886012
        @test MPSGE.Complementarity.result_value(m._jump_model[:U]) ≈ 1.04548206
        @test MPSGE.Complementarity.result_value(m._jump_model[:RA]) ≈ 157.321327225523
        @test MPSGE.Complementarity.result_value(m._jump_model[:PX]) ≈ 1.0000000000
        @test MPSGE.Complementarity.result_value(m._jump_model[:PY]) ≈ 1.00957658
        @test MPSGE.Complementarity.result_value(m._jump_model[:PU]) ≈ 1.00318206
        @test MPSGE.Complementarity.result_value(m._jump_model[:PL]) ≈ 0.95346259
        @test MPSGE.Complementarity.result_value(m._jump_model[:PK]) ≈ 1.04880885
        @test MPSGE.Complementarity.result_value(m._jump_model[Symbol("PL†X")]) ≈ 52.4404424085075
        @test MPSGE.Complementarity.result_value(m._jump_model[Symbol("PL†Y")]) ≈ 21.1770570584356
        @test MPSGE.Complementarity.result_value(m._jump_model[Symbol("PK†X")]) ≈ 47.6731294622795
        @test MPSGE.Complementarity.result_value(m._jump_model[Symbol("PK†Y")]) ≈ 28.877805079685
        @test MPSGE.Complementarity.result_value(m._jump_model[Symbol("PX†U")]) ≈ 100.318205802571
        @test MPSGE.Complementarity.result_value(m._jump_model[Symbol("PY†U")]) ≈ 49.6833066029729
    end

    @testset "TWOBYTWO (indexed version)" begin
        m = Model()
        goods = [:x, :y]
        factors = [:l, :k]
        factor = DenseAxisArray(Float64[50 50; 20 30], goods, factors)
        supply = DenseAxisArray(Float64[100, 50], goods)
        @parameter(m, endow, 1.0)
        Y = add!(m, Sector(:Y, indices=(goods,)))
        U = add!(m, Sector(:U))
        PC = add!(m, Commodity(:PC, indices=(goods,)))
        PU = add!(m, Commodity(:PU))
        PF = add!(m, Commodity(:PF, indices=(factors,)))
        RA = add!(m, Consumer(:RA, benchmark=150.))

        for i in goods
            @production(m, Y[i], 1, [Output(PC[i], supply[i])], [Input(PF[:l], factor[i,:l]), Input(PF[:k], factor[i,:k])])
        end
        @production(m, U, 1, [Output(PU, 150)], [Input(PC[:x], 100), Input(PC[:y], 50)])
        @demand(m, RA, [Demand(PU, 150)], [Endowment(PF[:l], :(70 * $endow)), Endowment(PF[:k], 80.)])

        solve!(m, cumulative_iteration_limit=0)

        @test MPSGE.Complementarity.result_value(m._jump_model[:Y][:x]) ≈ 1.
        @test MPSGE.Complementarity.result_value(m._jump_model[:Y][:y]) ≈ 1.
        @test MPSGE.Complementarity.result_value(m._jump_model[:U]) ≈ 1.
        @test MPSGE.Complementarity.result_value(m._jump_model[:RA]) ≈ 150.
        @test MPSGE.Complementarity.result_value(m._jump_model[:PC][:x]) ≈ 1.
        @test MPSGE.Complementarity.result_value(m._jump_model[:PC][:y]) ≈ 1.
        @test MPSGE.Complementarity.result_value(m._jump_model[:PU]) ≈ 1.
        @test MPSGE.Complementarity.result_value(m._jump_model[:PF][:l]) ≈ 1.
        @test MPSGE.Complementarity.result_value(m._jump_model[:PF][:k]) ≈ 1.
        @test MPSGE.Complementarity.result_value(m._jump_model[Symbol("PF[l]†Y[x]")]) ≈ 50.
        @test MPSGE.Complementarity.result_value(m._jump_model[Symbol("PF[l]†Y[y]")]) ≈ 20.
        @test MPSGE.Complementarity.result_value(m._jump_model[Symbol("PF[k]†Y[x]")]) ≈ 50.
        @test MPSGE.Complementarity.result_value(m._jump_model[Symbol("PF[k]†Y[y]")]) ≈ 30.
        @test MPSGE.Complementarity.result_value(m._jump_model[Symbol("PC[x]†U")]) ≈ 100.
        @test MPSGE.Complementarity.result_value(m._jump_model[Symbol("PC[y]†U")]) ≈ 50.

        set_fixed!(PC[:x], true)
        set_value(endow, 1.1)

        solve!(m)
        @test MPSGE.Complementarity.result_value(m._jump_model[:Y][:x]) ≈ 1.04880885
        @test MPSGE.Complementarity.result_value(m._jump_model[:Y][:y]) ≈ 1.03886012
        @test MPSGE.Complementarity.result_value(m._jump_model[:U]) ≈ 1.04548206
        @test MPSGE.Complementarity.result_value(m._jump_model[:RA]) ≈ 157.321327225523
        @test MPSGE.Complementarity.result_value(m._jump_model[:PC][:x]) ≈ 1.0000000000
        @test MPSGE.Complementarity.result_value(m._jump_model[:PC][:y]) ≈ 1.00957658
        @test MPSGE.Complementarity.result_value(m._jump_model[:PU]) ≈ 1.00318206
        @test MPSGE.Complementarity.result_value(m._jump_model[:PF][:l]) ≈ 0.95346259
        @test MPSGE.Complementarity.result_value(m._jump_model[:PF][:k]) ≈ 1.04880885
        @test MPSGE.Complementarity.result_value(m._jump_model[Symbol("PF[l]†Y[x]")]) ≈ 52.4404424085075
        @test MPSGE.Complementarity.result_value(m._jump_model[Symbol("PF[l]†Y[y]")]) ≈ 21.1770570584356
        @test MPSGE.Complementarity.result_value(m._jump_model[Symbol("PF[k]†Y[x]")]) ≈ 47.6731294622795
        @test MPSGE.Complementarity.result_value(m._jump_model[Symbol("PF[k]†Y[y]")]) ≈ 28.877805079685
        @test MPSGE.Complementarity.result_value(m._jump_model[Symbol("PC[x]†U")]) ≈ 100.318205802571
        @test MPSGE.Complementarity.result_value(m._jump_model[Symbol("PC[y]†U")]) ≈ 49.6833066029729

    end
    @testset "JPMGE (Joint Production Intermediate Demand)" begin

        m = Model()
        goods = [:g1, :g2]
        factors = [:l, :k]
        sectors = [:s1, :s2]
        make0 = DenseAxisArray(Float64[6 2; 2 10], goods, sectors)
        use0 = DenseAxisArray(Float64[4 2; 2 6], goods, sectors)
        fd0 = DenseAxisArray(Float64[1 3; 1 1], factors, sectors)
        c0 = DenseAxisArray(Float64[2, 4], goods)
        e0 = DenseAxisArray(Float64[sum(fd0[f,:]) for f in factors], factors)
        @parameter(m, endow, 1.0)
        X = add!(m, Sector(:X, indices=(sectors,)))
        P = add!(m, Commodity(:P, indices=(goods,)))
        PF = add!(m, Commodity(:PF, indices=(factors,)))
        Y = add!(m, Consumer(:Y, benchmark=sum(fd0)))
        for j in sectors
            @production(m, X[j], 1, [Output(P[i], make0[i,j]) for i in goods], [[Input(P[i], use0[i,j]) for i in goods]; [Input(PF[f], fd0[f,j]) for f in factors]])
        end
        @demand(m, Y, [Demand(P[i], c0[i]) for i in goods], [Endowment(PF[:k], e0[:k]), Endowment(PF[:l], :($(endow) * $(e0[:l]))) ])

        avm = algebraic_version(m)
        @test typeof(avm) == MPSGE.AlgebraicWrapper

        solve!(m)

        @test MPSGE.Complementarity.result_value(m._jump_model[:X][:s1]) ≈ 1.
        @test MPSGE.Complementarity.result_value(m._jump_model[:X][:s2]) ≈ 1.
        @test MPSGE.Complementarity.result_value(m._jump_model[:P][:g1]) ≈ 1.
        @test MPSGE.Complementarity.result_value(m._jump_model[:P][:g2]) ≈ 1.
        @test MPSGE.Complementarity.result_value(m._jump_model[:PF][:l]) ≈ 1.
        @test MPSGE.Complementarity.result_value(m._jump_model[:PF][:k]) ≈ 1.
        @test MPSGE.Complementarity.result_value(m._jump_model[:Y]) ≈ 6.
        @test MPSGE.Complementarity.result_value(m._jump_model[Symbol("P[g1]‡X[s1]")]) ≈ 6.
        @test MPSGE.Complementarity.result_value(m._jump_model[Symbol("P[g2]‡X[s1]")]) ≈ 2.
        @test MPSGE.Complementarity.result_value(m._jump_model[Symbol("P[g1]‡X[s2]")]) ≈ 2.
        @test MPSGE.Complementarity.result_value(m._jump_model[Symbol("P[g2]‡X[s2]")]) ≈10.
        @test MPSGE.Complementarity.result_value(m._jump_model[Symbol("P[g1]†X[s1]")]) ≈ 4.
        @test MPSGE.Complementarity.result_value(m._jump_model[Symbol("P[g2]†X[s1]")]) ≈ 2.
        @test MPSGE.Complementarity.result_value(m._jump_model[Symbol("PF[l]†X[s1]")]) ≈ 1.
        @test MPSGE.Complementarity.result_value(m._jump_model[Symbol("PF[k]†X[s1]")]) ≈ 1.
        @test MPSGE.Complementarity.result_value(m._jump_model[Symbol("P[g1]†X[s2]")]) ≈ 2.
        @test MPSGE.Complementarity.result_value(m._jump_model[Symbol("P[g2]†X[s2]")]) ≈ 6.
        @test MPSGE.Complementarity.result_value(m._jump_model[Symbol("PF[l]†X[s2]")]) ≈ 3.
        @test MPSGE.Complementarity.result_value(m._jump_model[Symbol("PF[k]†X[s2]")]) ≈ 1.
        @test MPSGE.Complementarity.result_value(m._jump_model[Symbol("P[g1]ρY")]) ≈ 2.
        @test MPSGE.Complementarity.result_value(m._jump_model[Symbol("P[g2]ρY")]) ≈ 4.
        solve!(m)

        #Counter-factual
        set_value(endow, 1.1)
        set_value(Y, 6.4)
        set_fixed!(Y, true)

        solve!(m)
        @test MPSGE.Complementarity.result_value(m._jump_model[:X][:s1]) ≈ 0.99692562
        @test MPSGE.Complementarity.result_value(m._jump_model[:X][:s2]) ≈ 1.09975731
        @test MPSGE.Complementarity.result_value(m._jump_model[:P][:g1]) ≈ 1.01306317
        @test MPSGE.Complementarity.result_value(m._jump_model[:P][:g2]) ≈ 0.99467892
        @test MPSGE.Complementarity.result_value(m._jump_model[:PF][:l]) ≈ 0.97665932
        @test MPSGE.Complementarity.result_value(m._jump_model[:PF][:k]) ≈ 1.0513495
        @test MPSGE.Complementarity.result_value(m._jump_model[:Y]) ≈ 6.4
        @test MPSGE.Complementarity.result_value(m._jump_model[Symbol("P[g1]‡X[s1]")]) ≈ 6.0271570795595
        @test MPSGE.Complementarity.result_value(m._jump_model[Symbol("P[g2]‡X[s1]")]) ≈ 1.97259368637181
        @test MPSGE.Complementarity.result_value(m._jump_model[Symbol("P[g1]‡X[s2]")]) ≈ 2.03066186698742
        @test MPSGE.Complementarity.result_value(m._jump_model[Symbol("P[g2]‡X[s2]")]) ≈ 9.96905522562002 #note - digits after 9.969 added from MPSGE.jl results bc GAMS not showing
        @test MPSGE.Complementarity.result_value(m._jump_model[Symbol("P[g1]†X[s1]")]) ≈ 3.98197684247759
        @test MPSGE.Complementarity.result_value(m._jump_model[Symbol("P[g2]†X[s1]")]) ≈ 2.02778707866119
        @test MPSGE.Complementarity.result_value(m._jump_model[Symbol("PF[l]†X[s1]")]) ≈ 1.03260012599366
        @test MPSGE.Complementarity.result_value(m._jump_model[Symbol("PF[k]†X[s1]")]) ≈ 0.959241925360056
        @test MPSGE.Complementarity.result_value(m._jump_model[Symbol("P[g1]†X[s2]")]) ≈ 1.96980110341978
        @test MPSGE.Complementarity.result_value(m._jump_model[Symbol("P[g2]†X[s2]")]) ≈ 6.018624497730061 #note - digits after 6.0186 added from MPSGE.jl results bc GAMS not showing
        @test MPSGE.Complementarity.result_value(m._jump_model[Symbol("PF[l]†X[s2]")]) ≈ 3.06483480444652
        @test MPSGE.Complementarity.result_value(m._jump_model[Symbol("PF[k]†X[s2]")]) ≈ 0.94903406
        @test MPSGE.Complementarity.result_value(m._jump_model[Symbol("P[g1]ρY")]) ≈ 2.1058245831306
        @test MPSGE.Complementarity.result_value(m._jump_model[Symbol("P[g2]ρY")]) ≈ 4.28949135539776

    end
end