text:
# Abstract types

Signed <: Integer
Signed <: Number
Signed <: AbstractFloat

#+

isa(1, Int64) # equivalent to typeof(1) <: Int64
isa(1, Integer) # equivalent to typeof(1) <: Integer
isa(1, AbstractFloat) # equivalent to typeof(1) <: AbstractFloat

#+

isabstracttype(Real)
isabstracttype(Float64)

#+

isconcretetype(Real)
isconcretetype(Float64)

# ## Composite types

struct Rectangle
    bottomleft::Vector{Float64}
    width
    height
end

#+

r = Rectangle([1,2], 3, 4)
isa(r, Rectangle)

#+

methods(Rectangle)

#+

r.width
getproperty(r, :width)

#+

area(r::Rectangle) = r.width * r.height

function vertices(r::Rectangle)
    x, y = r.bottomleft
    w, h = r.width, r.height
    return [[x, y], [x + w, y], [x + w, y + h], [x, y + h]]
end

#+

area(r)
vertices(r)

#+

fieldnames(Rectangle)
fieldnames(typeof(r))

# ## Mutable composite types

r.bottomleft = [2;2]

#+

r.bottomleft[1] = 5
r.bottomleft

#+

area(r)
vertices(r)

#+

mutable struct MutableRectangle
    bottomleft::Vector{Float64}
    width
    height
end

#+

mr = MutableRectangle([1,2], 3, 4)
isa(mr, MutableRectangle)

#+

mr.width = 1.5
setproperty!(mr, :height, 2.5)
mr

#+

const AbstractRectangle = Union{Rectangle, MutableRectangle}

Rectangle <: AbstractRectangle
MutableRectangle <: AbstractRectangle

#+

perimeter(r::AbstractRectangle) = 2*(r.width + r.height)

perimeter(r)
perimeter(mr)

# ## Parametric types

abstract type AbstractPoint{T} end

struct Point{T <: Real} <: AbstractPoint{T}
    x::T
    y::T
end

#+

isconcretetype(Point{Float64})

#+

Point{Float64} <: Point <: AbstractPoint
Point{Int64} <: Point <: AbstractPoint

#+

Point{Float64} <: Point{Real}
Point{Float64} <: AbstractPoint{Float64}
Point{Float64} <: AbstractPoint{Real}

#+

coordinates(p::Point{Real}) = (p.x, p.y)

coordinates(Point(1,2))
coordinates(Point(1.0,2.0))

#+

coordinates(p::Point{<:Real}) = (p.x, p.y)

coordinates(Point(1,2))
coordinates(Point(1.0,2.0))

#+

Base.show(io::IO, p::AbstractPoint) = print(io, coordinates(p))

Point(4, 2)
Point(0.2, 1.3)

#+

Point(1, 2)
Point{Float32}(1, 2)

#+

Point(1, 2.0)

# ### Exercise:
# Define a structure that represents 3D-points. Do not forget to define it as a subtype of
# the AbstractPoint type. Then add a new method to the `coordinates` function.
#
# ---
# ### Solution:

# ---
#
# ## Constructors

Point(x::Real, y::Real) = Point(promote(x, y)...)
Point(1, 2.0)

#+

typeof(Point(1, 2.0))

#+

struct OrderedPair{T <: Real}
    x::T
    y::T

    function OrderedPair(x::Real, y::Real)
        x > y && error("the first argument must be less than or equal to the second one")
        xp, yp = promote(x, y)
        return new{typeof(xp)}(xp, yp)
    end
end

#+

OrderedPair(1,2)
OrderedPair(2,1)

# ### Exercise:
# Define a structure that represents ND-points and stores their coordinates as `Tuple`. Do
# not forget to define it as a subtype of the `AbstractPoint` type. Redefine the default
# inner constructor to create an instance of `PointND` from different types. Then add a new
# method to the `coordinates` function, and define function `dim` that returns the
# dimension of the point.
#
# **Hints:** use the `new` function in the definition of the new inner constructor.
#
# **Bonus:** Tuples with elements of the same type can be described by the special type
# `NTuple{N, T}`, where `N` is the number of elements and `T` their type.
#
# ---
# ### Solution:

# ---
#
# ## Default field values

Base.@kwdef struct MyType
    a::Int # required keyword
    b::Float64 = 2.3
    c::String = "hello"
end

methods(MyType)

#+

MyType(1, 2.3, "aaa")
MyType(; a = 3)
MyType(; a = 5, b = 4.5)

#+

(m::MyType)() = (m.a, m.b, m.c)

m = MyType(; a = 5, b = 4.5)
m()

#+

(m::MyType)(x::Real) = m.a*x + m.b
(m::MyType)(x::String) = "$(m.c), $(x)"

m(1)
m("world")

# ### Exercise:
# [Gaussian distribution](https://en.wikipedia.org/wiki/Normal_distribution) is uniquely
# represented by its mean $\mu$ and variance $\sigma^2>0$. Write a structure `Gauss`
# with the proper fields and an inner constructor that checks if the input parameters are
# correct. Initialization without arguments `Gauss()` should return the standardized normal
# distribution ($\mu = 0$ and $ \sigma = 1$).  Define a functor that computes the
# probability density function at a given point defined by
# $$
# f_{\mu, \sigma}(x) = \frac{1}{\sigma \sqrt{ 2\pi }} \exp\left\{ -\frac{1}{2} \left( \frac
# {x - \mu}{\sigma} \right) ^2 \right\},
# $$
# Verify that the probability density function is defined correctly, i.e., its integral
# equals 1.
#
# ---
# ### Solution:

# ---
#

using RecipesBase

@recipe function f(d::Gauss, x = (d.μ - 4d.σ):0.1:(d.μ + 4d.σ))
    seriestype  :=  :path
    label --> "Gauss(μ = $(d.μ), σ = $(d.σ))"
    xguide --> "x"
    yguide --> "f(x)"
    linewidth --> 2
    return x, d.(x)
end

#+

using Plots

plot(Gauss());
plot!(Gauss(4, 2); linewidth = 4, color = :red);
plot!(Gauss(-3, 2); label = "new label", linestyle = :dash)


text:
using Pkg
cd(@__DIR__)
Pkg.activate("."); Pkg.instantiate()

using Plots; gr()
using LaTeXStrings
using DataFrames
using CSV
using OrdinaryDiffEq
using LinearAlgebra
using DataInterpolations
using BenchmarkTools
using DiffEqFlux
using Optim
using Flux
using Zygote
using ApproxFun

# Model parameters
rc = 1.5e-2 # m
h = 0.5e-2 # m
m = 0.5 # kg
Jbc = 15e-3 # kg.m2
Js = 5e-3 # kg.m2
g = 9.81 # m/s2
K_th_th = 55e-3 # Nm/V
K_th_psi = 5e-3 # Nm/V
K_psi_th = 0.15 # Nm/V
K_psi_psi = 0.2 # Nm/V
D_th = 1e-2 # Nm/(rad/s)
D_psi = 8e-2 # Nm/(rad/s)

# Graphics
figpath = "../figs/"
LW1 = 2.5
LW2 = 1.5
LS1 = :solid
LS2 = :dot
LS3 = :dash
LC1 = :blue
LC2 = :red
LC3 = :green
LC4 = :black
LA1 = 1.0
LA2 = 0.7
LA3 = 0.3

# Parse data from table
df=CSV.read("../data/Lab-Helicopter_Experimental-data.csv")
ENV["COLUMNS"]=100 # change the number of IJulia columns from default 80 to 100
df[1:5,:]
u_th = df[:,1]
u_psi = df[:,2]
theta = df[:,3]
psi = df[:,4]
dt = df[1,5]
tm = 0:dt:(length(psi)-1)*dt
tspan = (0.0, (length(psi)-1)*dt)
Nd = length(tm);
timeseriesdata = [theta';psi']

A = fill(FourierBasis(2), 6)
nn1 = TensorLayer(A, 2)
N1 = size(nn1.p)[1]
α = randn(N1+6)/1000

function helicopter_uode(dx,x,Ps,P0,u,p,t)
# P0 = (rc=rc,h=h,m=m,Jbc=Jbc,Js=Js,g=g,K_th_th=K_th_th,K_th_psi=K_th_psi,K_psi_th=K_psi_th,K_psi_psi=K_psi_psi,D_th=D_th,D_psi=D_psi);
    th = x[1]
    psi = x[2]
    w_th = x[3]
    w_psi = x[4]

    U = [u.th(t),u.psi(t)]
    #
    M11 = Ps[4]*P0.Jbc + Ps[3]*P0.m*((Ps[1]*P0.rc)^2+(Ps[2]*P0.h)^2)
    M22 = Ps[4]*P0.Jbc*cos(th)^2 + Ps[3]*P0.m*(Ps[1]*P0.rc*cos(th)-Ps[2]*P0.h*sin(th))^2 + Ps[5]*P0.Js
    #
    dLdth = -w_psi^2*(((Ps[4]*P0.Jbc)+(Ps[3]*P0.m)*((Ps[1]*P0.rc)^2-(Ps[2]*P0.h)^2))*cos(th)*sin(th)-(Ps[1]*P0.rc)*(Ps[2]*P0.h)*(cos(th)^2-sin(th)^2)) - (Ps[3]*P0.m)*(Ps[6]*P0.g)*((Ps[1]*P0.rc)*cos(th) - (Ps[2]*P0.h)*sin(th))
    #
    dp_psi_dth = -2*(((Ps[4]*P0.Jbc)+(Ps[3]*P0.m)*((Ps[1]*P0.rc)^2-(Ps[2]*P0.h)^2))*cos(th)*sin(th) + (Ps[3]*P0.m)*(Ps[1]*P0.rc)*(Ps[2]*P0.h)*(cos(th)^2-sin(th)^2))*w_psi
    #
    F_th = Ps[7]*P0.K_th_th*U[1] - Ps[8]*P0.K_th_psi*U[2] - P0.D_th*(w_th) + nn1([p[1]*U[1], p[2]*U[2], p[3]*th, p[4]*psi, p[5]*w_th, p[6]*w_psi], p[7:end])[1]
    F_psi = Ps[9]*P0.K_psi_th*U[1] - Ps[10]*P0.K_psi_psi*U[2] - P0.D_psi*(w_psi) + nn1([p[1]*U[1], p[2]*U[2], p[3]*th, p[4]*psi, p[5]*w_th, p[6]*w_psi], p[7:end])[2]

    #
    dx[1] = w_th
    dx[2] = w_psi
    dx[3] = (dLdth + F_th)/M11
    dx[4] = (-dp_psi_dth*w_th + F_psi)/M22
end

P0 = (rc=rc,h=h,m=m,Jbc=Jbc,Js=Js,g=g,K_th_th=K_th_th,K_th_psi=K_th_psi,K_psi_th=K_psi_th,K_psi_psi=K_psi_psi,D_th=D_th,D_psi=D_psi)
Ps = fill(1.0,length(P0))
#
th0 = theta[1]
psi0 = psi[1]
w_th0 = (theta[2]-theta[1])/dt
w_psi0 = (psi[2]-psi[1])/dt
#
x0 = [th0,psi0,w_th0,w_psi0]
U_th = ZeroSpline(u_th,tm)
U_psi = ZeroSpline(u_psi,tm)
u_t = (th=U_th, psi=U_psi)
helicopter! = (dx,x,p,t) -> helicopter_uode(dx,x,Ps,P0,u_t,p,t)

prob = ODEProblem(helicopter!,x0,tspan,p=nothing)
function uode_cost(p)
    #
    sol = solve(prob, p=p, AutoTsit5(TRBDF2(autodiff=false)), saveat=dt)
    #
    return norm(sol[1:2,1:end]-timeseriesdata[:,1:size(sol,2)])^2/Nd
    #
end

rad2deg(t,u) = (t,u*180/pi)

callback = function (p, l)
    @show p, l
    prob = ODEProblem(helicopter!,x0,tspan,p)
    sol = solve(prob, AutoTsit5(TRBDF2(autodiff=false)), saveat=dt)
    pl = plot(sol,vars=(rad2deg,0,2),lw=LW1,lc=LC1)
    plot!(pl,tm,psi*180/pi,lw=LW1,ls=LS2,lc=LC2)
    plot!(pl,xlim = tspan)
    plot!(pl)
    display(pl)
    return false
end
res = DiffEqFlux.sciml_train(uode_cost, α, ADAM(0.01), cb = callback, maxiters = 150)
p = res.minimizer

prob = ODEProblem(helicopter!,x0,tspan,p)
sol_final = solve(prob, AutoTsit5(TRBDF2(autodiff=false)), saveat=dt)

plot(sol_final,vars=(rad2deg,0,1),lw=LW1,lc=LC1,label=L"$\theta$")
plot!(tm,theta*180/pi,lw=LW1,ls=LS2,lc=LC2,label=L"\theta^\mathrm{d}")
plot!(xlim = tspan,xlabel=L"time $t$ [s]", ylabel=L"$\theta$ [$^\circ$]")
plot!(title="Pitch angle: model+TensorLayer(u) (blue) vs. data (red)",box=true)
figname="Helicopter_pitch-angle_model_fit_tensor_layer_u.svg"
savefig(figpath*figname)

plot(sol_final,vars=(rad2deg,0,2),lw=LW1,lc=LC1,label=L"$\theta$")
plot!(tm,psi*180/pi,lw=LW1,ls=LS2,lc=LC2,label=L"\theta^\mathrm{d}")
plot!(xlim = tspan,xlabel=L"time $t$ [s]", ylabel=L"$\theta$ [$^\circ$]")
plot!(title="Yaw angle: model+TensorLayer(u) (blue) vs. data (red)",box=true)
figname="Helicopter_yaw-angle_model_fit_tensor_layer_u.svg"
savefig(figpath*figname)

ext: jl


text:
#!/usr/bin/env julia
using Test, DynamicBoundspODEsDiscrete, McCormick, DynamicBoundsBase,
      DataStructures
using DynamicBoundspODEsDiscrete.StaticTaylorSeries
using DiffResults: JacobianResult

const DBB = DynamicBoundsBase
const DR = DynamicBoundspODEsDiscrete

@testset "Discretize and Relax" begin

    struct unit_test_name <: DR.AbstractStateContractorName
    end

    @test_throws ErrorException DR.state_contractor_k(unit_test_name())
    @test_throws ErrorException DR.state_contractor_γ(unit_test_name())
    @test_throws ErrorException DR.state_contractor_steps(unit_test_name())

    # test improvement condition for existence & uniqueness
    Yold = [Interval(1.0, 3.0); Interval(2.0, 4.0); Interval(1.0, 3.0)]
    Ynew = [Interval(1.5, 2.0); Interval(3.0, 3.5); Interval(0.5, 3.5)]
    nx_ic = 2
    @test DR.improvement_condition(Yold, Ynew, nx_ic)

    function J!(out, x, p, t)
        out[1, 1] = x
        nothing
    end

    function f!(dx, x, p, t)
        dx[1] = x[1]^2 + p[2]
        dx[2] = x[2] + p[1]^2
        nothing
    end
    np = 2
    nx = 2
    k = 3
    x = [0.1; 1.0]
    p = [0.2; 0.1]
    jtf! = DR.JacTaylorFunctor!(f!, nx, np, Val(k), Interval{Float64}(0.0), 0.0)

    xIntv = Interval{Float64}.(x)
    pIntv = Interval{Float64}.(p)
    yIntv = [xIntv; pIntv]
    DR.jacobian_taylor_coeffs!(jtf!, xIntv, pIntv, 0.0)

    jac = JacobianResult(jtf!.out, yIntv).derivs[1]
    tjac = zeros(Interval{Float64}, nx + np, nx * (k + 1))
    Jx = Matrix{Interval{Float64}}[
        zeros(Interval{Float64}, nx, nx) for i = 1:(k+1)
    ]
    Jp = Matrix{Interval{Float64}}[
        zeros(Interval{Float64}, nx, np) for i = 1:(k+1)
    ]

    DR.set_JxJp!(jtf!, xIntv, pIntv, 0.0)
    @test isapprox(jtf!.Jp[2][2, 1].lo, 0.4, atol = 1E-3)
    @test isapprox(jtf!.Jp[2][1, 2].lo, 1.0, atol = 1E-3)
    @test isapprox(jtf!.Jp[4][2, 1].lo, 0.0666666, atol = 1E-3)
    @test isapprox(jtf!.Jp[4][1, 2].lo, 0.079999, atol = 1E-3)
    @test isapprox(jtf!.Jx[2][1, 1].lo, 0.2, atol = 1E-3)
    @test isapprox(jtf!.Jx[2][2, 2].lo, 1.0, atol = 1E-3)
    @test isapprox(jtf!.Jx[4][1, 1].lo, 0.030666, atol = 1E-3)
    @test isapprox(jtf!.Jx[4][2, 2].lo, 0.1666, atol = 1E-3)

    # make/evaluate interval valued Taylor cofficient functor
    itf! = DR.TaylorFunctor!(
        f!,
        nx,
        np,
        Val(k),
        zero(Interval{Float64}),
        zero(Float64),
    )
    outIntv = Vector{Interval{Float64}}[zeros(Interval{Float64}, 2) for i = 1:4]
    itf!(outIntv, xIntv, pIntv, 0.0)
    @test isapprox(outIntv[1][1].lo, 0.10001, atol = 1E-3)
    @test isapprox(outIntv[2][2].lo, 1.0399999999999998, atol = 1E-3)
    @test isapprox(outIntv[3][1].lo, 0.011, atol = 1E-3)
    @test isapprox(outIntv[4][2].lo, 0.173334, atol = 1E-3)
    @test isapprox(outIntv[1][2].hi, 1.0, atol = 1E-3)
    @test isapprox(outIntv[2][1].hi, 0.1100000000000000, atol = 1E-3)
    @test isapprox(outIntv[3][2].hi, 0.52, atol = 1E-3)
    @test isapprox(outIntv[4][1].hi, 0.004766666666666669, atol = 1E-3)

    # make/evaluate real valued Taylor cofficient functor
    rtf! = DR.TaylorFunctor!(f!, nx, np, Val(k), zero(Float64), zero(Float64))
    out = Vector{Float64}[zeros(Float64, 2) for i = 1:4]
    rtf!(out, x, p, 0.0)
    @test isapprox(out[1][1], 0.10001, atol = 1E-3)
    @test isapprox(out[2][1], 0.11000000000000001, atol = 1E-3)
    @test isapprox(out[3][1], 0.011, atol = 1E-3)
    @test isapprox(out[4][1], 0.004766666666666668, atol = 1E-3)
    @test isapprox(out[1][2], 1.0, atol = 1E-3)
    @test isapprox(out[2][2], 1.04, atol = 1E-3)
    @test isapprox(out[3][2], 0.52, atol = 1E-3)
    @test isapprox(out[4][2], 0.17333333333333334, atol = 1E-3)

    # higher order existence tests
    hⱼ = 0.001
    hmin = 0.00001
    function euf!(out, x, p, t)
        out[1, 1] = -x[1]^2
        nothing
    end

    jtf_exist_unique! =
        DR.JacTaylorFunctor!(euf!, 1, 1, Val(k), Interval{Float64}(0.0), 0.0)
    xIntv_plus = xIntv .+ Interval(0, 1)
    DR.jacobian_taylor_coeffs!(jtf_exist_unique!, xIntv_plus, pIntv, 0.0)
    @test isapprox(
        jtf_exist_unique!.result.value[4].lo,
        -1.464100000000001,
        atol = 1E-5,
    )
    @test isapprox(
        jtf_exist_unique!.result.value[4].hi,
        -9.999999999999999e-5,
        atol = 1E-5,
    )
    @test isapprox(
        jtf_exist_unique!.result.derivs[1][3, 1].lo,
        0.03,
        atol = 1E-5,
    )
    @test isapprox(
        jtf_exist_unique!.result.derivs[1][3, 1].hi,
        3.6300000000000012,
        atol = 1E-5,
    )

    coeff_out = zeros(Interval{Float64}, 1, k)
    DR.set_JxJp!(jtf_exist_unique!, xIntv_plus, pIntv, 0.0)
    @test isapprox(jtf_exist_unique!.Jx[4][1, 1].lo, -5.32401, atol = 1E-5)
    @test isapprox(jtf_exist_unique!.Jx[4][1, 1].hi, -0.00399999, atol = 1E-5)
    @test jtf_exist_unique!.Jp[1][1, 1].lo ==
          jtf_exist_unique!.Jp[1][1, 1].hi ==
          0.0

    #u_result = DR.UniquenessResult(1,1)
    #DR.existence_uniqueness!(u_result, itf_exist_unique!, eufY, hⱼ, hmin,
    #                         coeff_out_exist_unique!, Jx, Jp)

    #@test u_result.step == 0.001
    #@test u_result.confirmed

    Y = [1.1 3.2; 4.0 -1.0]
    y = [1.1; 3.2]
    A = [1.1 3.2; 4.0 -1.0]
    B = [1.1 3.2; 4.0 -1.0]
    b = [1.1; -1.0]
    mul_split!(Y, A, B, 2)
    @test isapprox(Y[2,2], 13.8, atol=1E-5)
    mul_split!(y, A, b, 2)
    @test isapprox(Y[1,2], 0.3200000000000003, atol=1E-5)
end

if !(VERSION < v"1.1" && testfile == "intervals.jl")
    using TaylorSeries

    function test_vs_Taylor1(x, y)
        flag = true
        for i = 0:2
            if x[i] !== y[i]
                flag = false
                break
            end
        end
        flag
    end

    @testset "Tests for STaylor1 expansions" begin

        @test STaylor1 <: DR.StaticTaylorSeries.AbstractSeries
        @test STaylor1{1,Float64} <: DR.StaticTaylorSeries.AbstractSeries{Float64}
        @test STaylor1([1.0, 2.0]) == STaylor1((1.0, 2.0))
        @test STaylor1(STaylor1((1.0, 2.0))) == STaylor1((1.0, 2.0))
        @test STaylor1(1.0, Val(2)) == STaylor1((1.0, 0.0, 0.0))

        @test +STaylor1([1.0, 2.0, 3.0]) == STaylor1([1.0, 2.0, 3.0])
        @test -STaylor1([1.0, 2.0, 3.0]) == -STaylor1([1.0, 2.0, 3.0])
        @test STaylor1([1.0, 2.0, 3.0]) + STaylor1([3.0, 2.0, 3.0]) ==
              STaylor1([4.0, 4.0, 6.0])
        @test STaylor1([1.0, 2.0, 3.0]) - STaylor1([3.0, 2.0, 4.0]) ==
              STaylor1([-2.0, 0.0, -1.0])
        @test STaylor1([1.0, 2.0, 3.0]) + 2.0 == STaylor1([3.0, 2.0, 3.0])
        @test STaylor1([1.0, 2.0, 3.0]) - 2.0 == STaylor1([-1.0, 2.0, 3.0])
        @test 2.0 + STaylor1([1.0, 2.0, 3.0]) == STaylor1([3.0, 2.0, 3.0])
        @test 2.0 - STaylor1([1.0, 2.0, 3.0]) == STaylor1([1.0, -2.0, -3.0])

        @test zero(STaylor1([1.0, 2.0, 3.0])) == STaylor1([0.0, 0.0, 0.0])
        @test one(STaylor1([1.0, 2.0, 3.0])) == STaylor1([1.0, 0.0, 0.0])

        @test_throws ArgumentError STaylor1([1.1, 2.1])/STaylor1([0.0, 2.1])

        @test isinf(STaylor1([Inf, 2.0, 3.0])) &&
              ~isinf(STaylor1([0.0, 0.0, 0.0]))
        @test isnan(STaylor1([NaN, 2.0, 3.0])) &&
              ~isnan(STaylor1([1.0, 0.0, 0.0]))
        @test iszero(STaylor1([0.0, 0.0, 0.0])) &&
              ~iszero(STaylor1([0.0, 1.0, 0.0]))

        @test length(STaylor1([0.0, 0.0, 0.0])) == 3
        @test size(STaylor1([0.0, 0.0, 0.0])) == 3
        @test firstindex(STaylor1([0.0, 0.0, 0.0])) == 0
        @test lastindex(STaylor1([0.0, 0.0, 0.0])) == 2

        st1 = STaylor1([1.0, 2.0, 3.0])
        @test st1(2.0) == 41.0
        @test st1() == 1.00
        st2 = typeof(st1)[st1; st1]
        @test st2(2.0)[1] == st2(2.0)[2] == 41.0
        @test st2()[1] == st2()[2] == 1.0
        @test StaticTaylorSeries.evaluate(st1, 2.0) == 41.0
        @test StaticTaylorSeries.evaluate(st1) == 1.00
        @test StaticTaylorSeries.evaluate(st2, 2.0)[1] ==
              StaticTaylorSeries.evaluate(st2, 2.0)[2] ==
              41.0
        @test StaticTaylorSeries.evaluate(st2)[1] ==
              StaticTaylorSeries.evaluate(st2)[2] ==
              1.0

        @test view(typeof(STaylor1([1.1, 2.1]))[STaylor1([1.1, 2.1]) STaylor1([1.1, 2.1]); STaylor1([1.1, 2.1]) STaylor1([1.1, 2.1])], :, 1)(0.0) == Float64[1.1; 1.1]
        @test view(typeof(STaylor1([1.1, 2.1]))[STaylor1([1.1, 2.1]) STaylor1([1.1, 2.1]); STaylor1([1.1, 2.1]) STaylor1([1.1, 2.1])], :, 1)() == Float64[1.1; 1.1]

        # check that STaylor1 and Taylor yeild same result
        t1 = STaylor1([1.1, 2.1, 3.1])
        t2 = Taylor1([1.1, 2.1, 3.1])
        for f in (exp, abs, log, sin, cos, sinh, cosh, mod2pi, sqrt, abs2,
                  deg2rad, rad2deg)
            @test test_vs_Taylor1(f(t1), f(t2))
        end

        @test DR.StaticTaylorSeries.get_order(t1) == 2
        @test axes(t1) isa Tuple{}
        @test iterate(t1)[1] == 1.10
        @test iterate(t1)[2] == 1
        @test eachindex(t1) == 0:2
        @test t1[:] == (1.10, 2.10, 3.10)
        @test t1[1:2] == (2.10, 3.10)
        @test t1[0:2:2] == (1.10, 3.10)
        @test rem(t1, 2) == t1
        @test mod(t1, 2) == t1

        @test STaylor1([1.1, 2.1, 3.1], Val(3),  Val(5)) == STaylor1([1.1, 2.1, 3.1, 0.0, 0.0, 0.0])

        t1_mod = mod(t1, 2.0)
        t2_mod = mod(t2, 2.0)
        @test isapprox(t1_mod[0], t2_mod[0], atol = 1E-10)
        @test isapprox(t1_mod[1], t2_mod[1], atol = 1E-10)
        @test isapprox(t1_mod[2], t2_mod[2], atol = 1E-10)

        t1_rem = rem(t1, 2.0)
        t2_rem = rem(t2, 2.0)
        @test isapprox(t1_rem[0], t2_rem[0], atol = 1E-10)
        @test isapprox(t1_rem[1], t2_rem[1], atol = 1E-10)
        @test isapprox(t1_rem[2], t2_rem[2], atol = 1E-10)

        t1a = STaylor1([2.1, 2.1, 3.1])
        t2a = Taylor1([2.1, 2.1, 3.1])

        for test_tup in (
            (/, t1, t1a, t2, t2a),
            (*, t1, t1a, t2, t2a),
            (/, t1, 1.3, t2, 1.3),
            (*, t1, 1.3, t2, 1.3),
            (+, t1, 1.3, t2, 1.3),
            (-, t1, 1.3, t2, 1.3),
            (*, 1.3, t1, 1.3, t2),
            (+, 1.3, t1, 1.3, t2),
            (-, 1.3, t1, 1.3, t2),
            (*, 1.3, t1, 1.3, t2),
            (^, t1, 0, t2, 0),
            (^, t1, 1, t2, 1),
            (^, t1, 2, t2, 2),
            (^, t1, 3, t2, 3),
            (^, t1, 4, t2, 4),
            (/, 1.3, t1, 1.3, t2),
            #(^, t1, -1, t2, -1),
            #(^, t1, -2, t2, -2),
            #(^, t1, -3, t2, -3),
            (^, t1, 0.6, t2, 0.6),
            (^, t1, 1 / 2, t2, 1 / 2),
        )
            temp1 = test_tup[1](test_tup[2], test_tup[3])
            temp2 = test_tup[1](test_tup[4], test_tup[5])
            check1 = isapprox(temp1[0], temp2[0], atol = 1E-10)
            check2 = isapprox(temp1[1], temp2[1], atol = 1E-10)
            check3 = isapprox(temp1[2], temp2[2], atol = 1E-10)
            @test check1
            @test check2
            @test check3
            if !check1 || !check2 || !check3
                println("$test_tup, $temp1, $temp2")
            end
        end

        @test isapprox(
            StaticTaylorSeries.square(t1)[0],
            (t2^2)[0],
            atol = 1E-10,
        )
        @test isapprox(
            StaticTaylorSeries.square(t1)[1],
            (t2^2)[1],
            atol = 1E-10,
        )
        @test isapprox(
            StaticTaylorSeries.square(t1)[2],
            (t2^2)[2],
            atol = 1E-10,
        )

        a = STaylor1([0.0, 1.2, 2.3, 4.5, 0.0])
        @test findfirst(a) == 1
        @test findlast(a) == 3

        eval_staylor = StaticTaylorSeries.evaluate(a, Interval(1.0, 2.0))
        @test isapprox(eval_staylor.lo, 7.99999, atol = 1E-4)
        @test isapprox(eval_staylor.hi, 47.599999999999994, atol = 1E-4)

        a = STaylor1([5.0, 1.2, 2.3, 4.5, 0.0])
        @test isapprox(deg2rad(a)[0], 0.087266, atol = 1E-5)
        @test isapprox(deg2rad(a)[2], 0.040142, atol = 1E-5)
        @test isapprox(rad2deg(a)[0], 286.4788975, atol = 1E-5)
        @test isapprox(rad2deg(a)[2], 131.7802928, atol = 1E-5)
        @test real(a) == STaylor1([5.0, 1.2, 2.3, 4.5, 0.0])
        @test imag(a) == STaylor1([0.0, 0.0, 0.0, 0.0, 0.0])
        @test adjoint(a) == STaylor1([5.0, 1.2, 2.3, 4.5, 0.0])
        @test conj(a) == STaylor1([5.0, 1.2, 2.3, 4.5, 0.0])
        @test a == abs(a)
        @test a == abs(-a)

        @test convert(
            STaylor1{3,Float64},
            STaylor1{3,Float64}((1.1, 2.2, 3.3)),
        ) == STaylor1{3,Float64}((1.1, 2.2, 3.3))
        @test convert(STaylor1{3,Float64}, 1) == STaylor1(1.0, Val(3))
        @test convert(STaylor1{3,Float64}, 1.2) == STaylor1(1.2, Val(3))

        #ta(a) = STaylor1(1.0, Val(15))
        @test_broken promote(1.0, STaylor1(1.0, Val(15)))[1] == STaylor1(1.0, Val(16))
        @test_broken promote(0, STaylor1(1.0, Val(15)))[1] == STaylor1(0.0, Val(16))
        @test_broken eltype(promote(STaylor1(1, Val(15)), 2)[2]) == Int
        @test_broken eltype(promote(STaylor1(1.0, Val(15)), 1.1)[2]) == Float64
        @test_broken eltype(promote(0, STaylor1(1.0, Val(15)))[1]) == Float64

        @test_broken promote_rule(typeof(STaylor1([1.1, 2.1])), typeof(STaylor1([1.1, 2.1]))) == STaylor1{2,Float64}
        @test_broken promote_rule(typeof(STaylor1([1.1, 2.1])), typeof(STaylor1([1, 2]))) == STaylor1{2,Float64}
        @test_broken promote_rule(typeof(STaylor1([1.1, 2.1])), typeof([1.1, 2.1])) == STaylor1{2,Float64}
        @test_broken promote_rule(typeof(STaylor1([1.1, 2.1])), typeof([1, 2])) == STaylor1{2,Float64}
        @test_broken promote_rule(typeof(STaylor1([1.1, 2.1])), typeof(1.1)) == STaylor1{2,Float64}
        @test_broken promote_rule(typeof(STaylor1([1.1, 2.1])), typeof(1)) == STaylor1{2,Float64} #TODO: FAILING

        @test convert(STaylor1{2,Float64}, [1; 2]) == STaylor1(Float64[1, 2])
        @test convert(STaylor1{2,Float64}, [1.1; 2.1]) == STaylor1([1.1, 2.1])
        @test convert(STaylor1{2,Rational{Int64}}, STaylor1(BigFloat[0.5, 0.75])) == STaylor1([0.5, 0.75])

        @test isapprox(StaticTaylorSeries.normalize_taylor(STaylor1([1.1, 2.1]), Interval(1.0, 2.0))[0], 13.7, atol = 1E-8)
        @test isapprox(StaticTaylorSeries._normalize(STaylor1([1.1, 2.1]), Interval(1.0, 2.0), Val(true))[0], 13.7, atol = 1E-8)
        @test isapprox(StaticTaylorSeries._normalize(STaylor1([1.1, 2.1]), Interval(1.0, 2.0), Val(false))[0], 13.7, atol = 1E-8)

        @test_nowarn Base.show(stdout, a)
        @test StaticTaylorSeries.coeffstring(a, 5) == "0.0t^4"
        @test StaticTaylorSeries.coeffstring(a, 2) == "1.2t"

        @test_throws ArgumentError abs(STaylor1([0.0, 2.1]))
    end
end

@testset "Lohner's Method Interval Testset" begin

    ticks = 100.0
    steps = 100.0
    tend = steps / ticks

    x0(p) = [9.0]
    function f!(dx, x, p, t)
        dx[1] = p[1] - x[1]
        nothing
    end
    tspan = (0.0, tend)
    pL = [-1.0]
    pU = [1.0]
    prob = DBB.ODERelaxProb(f!, tspan, x0, pL, pU)

    integrator = DiscretizeRelax(
        prob,
        DR.LohnerContractor{7}(),
        h = 1 / ticks,
        repeat_limit = 1,
        skip_step2 = false,
        step_limit = steps,
        relax = false,
    )

    ratio = rand(1)
    pstar = pL .* ratio .+ pU .* (1.0 .- ratio)
    DBB.setall!(integrator, DBB.ParameterValue(), [0.0])
    DBB.relax!(integrator)

    lo_vec = getfield.(getindex.(integrator.storage[:], 1), :lo)
    hi_vec = getfield.(getindex.(integrator.storage[:], 1), :hi)

    @test isapprox(lo_vec[7], 8.417645335842485, atol = 1E-5)
    @test isapprox(hi_vec[7], 8.534116268673994, atol = 1E-5)
end

@testset "Lohner's Method Adaptive Interval Testset" begin

    ticks = 100.0
    steps = 100.0
    tend = steps / ticks

    x0(p) = [9.0]
    function f!(dx, x, p, t)
        dx[1] = p[1]  - x[1]*x[1]
        nothing
    end
    tspan = (0.0, tend)
    pL = [-1.0]
    pU = [1.0]
    prob = DBB.ODERelaxProb(f!, tspan, x0, pL, pU)

    integrator = DiscretizeRelax(
        prob,
        DR.LohnerContractor(7),
        repeat_limit = 1,
        skip_step2 = false,
        step_limit = steps,
        relax = false,
    )

    ratio = rand(1)
    pstar = pL .* ratio .+ pU .* (1.0 .- ratio)
    DBB.setall!(integrator, DBB.ParameterValue(), [0.0])
    DBB.relax!(integrator)

    lo_vec = getfield.(getindex.(integrator.storage[:], 1), :lo)
    hi_vec = getfield.(getindex.(integrator.storage[:], 1), :hi)

    @test isapprox(lo_vec[6], 5.884811279034315, atol = 1E-5)
    @test isapprox(hi_vec[6], 5.965419612738188, atol = 1E-5)

    support_set = DBB.get(integrator, DBB.SupportSet())
    outvec = zeros(length(support_set.s))

    DBB.getall!(outvec, integrator, DBB.Bound{Lower}())
    @test_broken isapprox(outvec[10], 0.9746606911231538, atol=1E-8)

    DBB.getall!(outvec, integrator, DBB.Bound{Upper}())
    @test_broken isapprox(outvec[10], 1.7009858838207106, atol=1E-8)

    DBB.getall!(outvec, integrator, DBB.Relaxation{Lower}())
    @test_broken isapprox(outvec[10], 0.9746606911231538, atol=1E-8)

    DBB.getall!(outvec, integrator, DBB.Relaxation{Upper}())
    @test_broken isapprox(outvec[10], 1.7009858838207106, atol=1E-8)
end

@testset "Lohner's Method MC Testset" begin

    ticks = 100.0
    steps = 100.0
    tend = steps / ticks

    x0(p) = [9.0]
    function f!(dx, x, p, t)
        dx[1] = p[1] - x[1]
        nothing
    end
    tspan = (0.0, tend)
    pL = [-0.3]
    pU = [0.3]
    prob = DynamicBoundsBase.ODERelaxProb(f!, tspan, x0, pL, pU)

    integrator = DiscretizeRelax(
        prob,
        DynamicBoundspODEsDiscrete.LohnerContractor{7}(),
        h = 1 / ticks,
        repeat_limit = 1,
        skip_step2 = false,
        step_limit = steps,
        relax = true,
    )

    ratio = rand(1)
    pstar = pL .* ratio .+ pU .* (1.0 .- ratio)
    DynamicBoundsBase.setall!(integrator, DynamicBoundsBase.ParameterValue(), [0.0])
    DynamicBoundsBase.relax!(integrator)

    lo_vec = getfield.(getfield.(getindex.(integrator.storage[:], 1), :Intv), :lo)
    hi_vec = getfield.(getfield.(getindex.(integrator.storage[:], 1), :Intv), :hi)

    @test isapprox(lo_vec[6], 8.546433647856642, atol = 1E-5)
    @test isapprox(hi_vec[6], 8.575695993156213, atol = 1E-5)

    support_set = DBB.get(integrator, DBB.SupportSet())
    out = Matrix{Float64}[]
    for i = 1:1
        push!(out, zeros(1,length(support_set.s)))
    end
    # DBB.getall!(out, integrator, DBB.Subgradient{Lower}())
    # @test isapprox(out[1][1,10], 0.08606881472877183, atol=1E-8)

    # DBB.getall!(out, integrator, DBB.Subgradient{Upper}())
    # @test isapprox(out[1][1,10], 0.08606881472877183, atol=1E-8)

    #out = zeros(1, length(support_set.s))
    #DBB.getall!(out, integrator, DBB.Bound{Lower}())
    #@test isapprox(out[1,10], 8.199560023022423, atol=1E-8)

    #DBB.getall!(out, integrator, DBB.Bound{Upper}())
    #@test isapprox(out[1,10], 8.251201311859687, atol=1E-8)

    #DBB.getall!(out, integrator, DBB.Relaxation{Lower}())
    #@test isapprox(out[1,10], 8.225380667441055, atol=1E-8) #

    #DBB.getall!(out, integrator, DBB.Relaxation{Upper}())
    #@test isapprox(out[1,10], 8.225380667441055, atol=1E-8) #

    #out = DBB.getall(integrator, DBB.Subgradient{Lower}())
    #@test isapprox(out[1][1,10], 0.08606881472877183, atol=1E-8)

    #out = DBB.getall(integrator, DBB.Subgradient{Upper}())
    #@test isapprox(out[1][1,10], 0.08606881472877183, atol=1E-8)

    #out = DBB.getall(integrator, DBB.Bound{Lower}())
    #@test isapprox(out[1,10], 8.199560023022423, atol=1E-8)

    #out = DBB.getall(integrator, DBB.Bound{Upper}())
    #@test isapprox(out[1,10], 8.251201311859687, atol=1E-8)

    #out = DBB.getall(integrator, DBB.Relaxation{Lower}())
    #@test isapprox(out[1,10], 8.225380667441055, atol=1E-8)

    #out = DBB.getall(integrator, DBB.Relaxation{Upper}())
    #@test isapprox(out[1,10], 8.225380667441055, atol=1E-8)

    #outvec = zeros(length(support_set.s))

    #DBB.getall!(outvec, integrator, DBB.Bound{Lower}())
    #@test isapprox(outvec[10], 8.199560023022423, atol=1E-8)

    #DBB.getall!(outvec, integrator, DBB.Bound{Upper}())
    #@test isapprox(outvec[10], 8.251201311859687, atol=1E-8)

    #DBB.getall!(outvec, integrator, DBB.Relaxation{Lower}())
    #@test isapprox(outvec[10], 8.225380667441055, atol=1E-8)

    #DBB.getall!(outvec, integrator, DBB.Relaxation{Upper}())
    #@test isapprox(outvec[10], 8.225380667441055, atol=1E-8)
end

@testset "Hermite-Obreshkoff Testset" begin

    ticks = 100.0
    steps = 100.0
    tend = steps / ticks

    x0(p) = [9.0]
    function f!(dx, x, p, t)
        dx[1] = p[1] - x[1]
        nothing
    end
    tspan = (0.0, tend)
    pL = [-1.0]
    pU = [1.0]
    prob = DBB.ODERelaxProb(f!, tspan, x0, pL, pU)

    integrator = DiscretizeRelax(prob, DR.HermiteObreschkoff(3, 3),
                                 h = 1/ticks, repeat_limit = 1,
                                 skip_step2 = false,
                                 step_limit = steps,
                                 relax = false)

    ratio = rand(1)
    pstar = pL .* ratio .+ pU .* (1.0 .- ratio)
    DBB.setall!(integrator, DBB.ParameterValue(), [0.0])
    DBB.relax!(integrator)

    lo_vec = getfield.(getindex.(integrator.storage[:], 1), :lo)
    hi_vec = getfield.(getindex.(integrator.storage[:], 1), :hi)

    @test isapprox(lo_vec[6], 8.510710620631412, atol = 1E-5)
    @test isapprox(hi_vec[6], 8.611419020357289, atol = 1E-5)
end

@testset "Wilhelm 2019 Integrator Testset" begin
    Y = Interval(0,5)
    X1 = Interval(1,2)
    X2 = Interval(-10,10)
    X3 = Interval(0,5)
    flag1 = DR.strict_x_in_y(X1,Y)
    flag2 = DR.strict_x_in_y(X2,Y)
    flag3 = DR.strict_x_in_y(X3,Y)
    @test flag1 == true
    @test flag2 == false
    @test flag3 == false

    A1 = Interval(0)
    A2 = Interval(0,3)
    A3 = Interval(-2,0)
    A4 = Interval(-3,2)
    ind1,B1,C1 = DR.extended_divide(A1)
    ind2,B2,C2 = DR.extended_divide(A2)
    ind3,B3,C3 = DR.extended_divide(A3)
    ind4,B4,C4 = DR.extended_divide(A4)
    @test ind1 == 0
    @test ind2 == 1
    @test ind3 == 2
    @test ind4 == 3
    @test B1 == Interval(-Inf,Inf)
    @test 0.33333 - 1E-4 <= B2.lo <= 0.33333 + 1E-4
    @test B2.hi == Inf
    @test B3 == Interval(-Inf,-0.5)
    @test B4.lo == -Inf
    @test -0.33333 - 1E-4 <= B4.hi <= -0.33333 + 1E-4
    @test C1 == Interval(-Inf,Inf)
    @test C2 == Interval(Inf,Inf)
    @test C3 == Interval(-Inf,-Inf)
    @test C4 == Interval(0.5,Inf)

    N =  Interval(-5,5)
    X = Interval(-5,5)
    Mii = Interval(-5,5)
    B = Interval(-5,5)
    rtol = 1E-4
    indx1,box11,box12 = DR.extended_process(N,X,Mii,B,rtol)
    Miib = Interval(0,5)
    S1b = Interval(1,5)
    S2b = Interval(1,5)
    Bb = Interval(1,5)
    indx2,box21,box22 = DR.extended_process(N,X,Mii,B,rtol)
    Miic = Interval(-5,0)
    S1c = Interval(1,5)
    S2c = Interval(1,5)
    Bc = Interval(1,5)
    indx3,box31,box32 = DR.extended_process(N,X,Mii,B,rtol)
    Miia = Interval(1,5)
    S1a = Interval(1,5)
    S2a = Interval(1,5)
    Ba = Interval(1,5)
    indx6,box61,box62 = DR.extended_process(N,X,Mii,B,rtol)
    Miid = Interval(0,0)
    S1d = Interval(1,5)
    S2d = Interval(1,5)
    Bd = Interval(1,5)
    indx8,box81,box82 = DR.extended_process(N,X,Mii,B,rtol)

    @test indx1 == 0
    @test box11 == Interval(-Inf,Inf)
    @test box12 == Interval(-5,5)

    @test indx2 == 0
    @test box21.hi > -Inf
    @test box22 == Interval(-5,5)

    @test indx3 == 0
    @test box31.lo < Inf
    @test box32 == Interval(-5,5)

    @test indx6 == 0
    @test box62.lo == -5.0
    @test box61.hi == Inf

    @test indx8 == 0
end

@testset "Discretize and Relax - Access Functions" begin

    use_relax = false
    lohners_type = 1
    prob_num = 1
    ticks = 100.0
    steps = 100.0
    tend = steps / ticks

    x0(p) = [1.2; 1.1]
    function f!(dx, x, p, t)
        dx[1] = p[1] * x[1] * (one(typeof(p[1])) - x[2])
        dx[2] = p[1] * x[2] * (x[1] - one(typeof(p[1])))
        nothing
    end
    tspan = (0.0, tend)
    pL = [2.95]
    pU = [3.05]

    prob = DBB.ODERelaxProb(f!, tspan, x0, pL, pU)

    integrator = DiscretizeRelax(
        prob,
        DynamicBoundspODEsDiscrete.LohnerContractor{7}(),
        h = 1 / ticks,
        repeat_limit = 1,
        skip_step2 = false,
        step_limit = steps,
        relax = use_relax,
    )

    @test DBB.supports(integrator, DBB.IntegratorName())
    @test !DBB.supports(integrator, DBB.Gradient())
    @test DBB.supports(integrator, DBB.Subgradient())
    @test DBB.supports(integrator, DBB.Bound())
    @test DBB.supports(integrator, DBB.Relaxation())
    @test DBB.supports(integrator, DBB.IsNumeric())
    @test DBB.supports(integrator, DBB.IsSolutionSet())
    @test DBB.supports(integrator, DBB.TerminationStatus())
    @test DBB.supports(integrator, DBB.Value())
    @test DBB.supports(integrator, DBB.ParameterValue())
    @test DBB.supports(integrator, DBB.SupportSet())

    @test DBB.get(integrator, DBB.IntegratorName()) == "Discretize & Relax Integrator"
    @test !DBB.get(integrator, DBB.IsNumeric())
    @test DBB.get(integrator, DBB.IsSolutionSet())
    @test DBB.get(integrator, DBB.TerminationStatus()) == RELAXATION_NOT_CALLED

    DBB.set!(integrator, DBB.SupportSet(Float64[i for i in range(0.0, tend, length = 200)]))

    ratio = rand(1)
    pstar = pL .* ratio .+ pU .* (1.0 .- ratio)
    DBB.setall!(integrator, DBB.ParameterValue(), [0.0])
    DBB.relax!(integrator)

    DBB.setall!(integrator, DBB.ParameterBound{Lower}(), [2.99])
    DBB.setall!(integrator, DBB.ParameterBound{Upper}(), [3.01])

    support_set = DBB.get(integrator, DBB.SupportSet())
    #@test support_set.s[3] == 0.02

    #=
    out = Matrix{Float64}[]
    for i in 1:1
        push!(out, zeros(1,length(support_set.s)))
    end
    DBB.getall!(out, integrator, DBB.Subgradient{Lower}())
    @test out[1][1,10] == 0.0

    DBB.getall!(out, integrator, DBB.Subgradient{Upper}())
    @test out[1][1,10] == 0.0

    out = zeros(1,length(support_set.s))
    DBB.getall!(out, integrator, DBB.Bound{Lower}())
    @test isapprox(out[1,10], 1.1507186500504751, atol=1E-8)

    DBB.getall!(out, integrator, DBB.Bound{Upper}())
    @test isapprox(out[1,10], 1.1534467709985823, atol=1E-8)

    DBB.getall!(out, integrator, DBB.Relaxation{Lower}())
    @test isapprox(out[1,10], 1.1507186500504751, atol=1E-8)

    DBB.getall!(out, integrator, DBB.Relaxation{Upper}())
    @test isapprox(out[1,10], 1.1534467709985823, atol=1E-8)

    out = DBB.getall(integrator, DBB.Subgradient{Lower}())
    @test out[1][1,10] == 0.0

    out = DBB.getall(integrator, DBB.Subgradient{Upper}())
    @test out[1][1,10] == 0.0

    out = DBB.getall(integrator, DBB.Bound{Lower}())
    @test isapprox(out[1,10], 1.1507186500504751, atol=1E-8)

    out = DBB.getall(integrator, DBB.Bound{Upper}())
    @test isapprox(out[1,10], 1.1534467709985823, atol=1E-8)

    out = DBB.getall(integrator, DBB.Relaxation{Lower}())
    @test isapprox(out[1,10], 1.1507186500504751, atol=1E-8)

    out = DBB.getall(integrator, DBB.Relaxation{Upper}())
    @test isapprox(out[1,10], 1.1534467709985823, atol=1E-8)
    =#
end


text:
using Kinetic, ProgressMeter, Plots, LinearAlgebra, JLD2

###
# dataset
###

X = Float32.([[1e-4, 0.0, 1e-4]; zeros(3); 1.0])
Y = Float32.([1.0, 0.0])

function init_field!(ks, ctr, face)
    for i in eachindex(ctr)
        prim = [2.0 * rand(), 0.0, 1 / rand()]

        ctr[i].prim .= prim
        ctr[i].w .= prim_conserve(prim, ks.gas.γ)
        ctr[i].f .= maxwellian(ks.vs.u, prim)
    end
    for i in eachindex(face)
        face[i].fw .= 0.0
        face[i].ff .= 0.0
    end
end

function regime_data(ks, w, prim, sw, f)
    Mu, Mxi, _, _1 = gauss_moments(prim, ks.gas.K)
    a = pdf_slope(prim, sw, ks.gas.K)
    swt = -prim[1] .* moments_conserve_slope(a, Mu, Mxi, 1)
    A = pdf_slope(prim, swt, ks.gas.K)
    tau = vhs_collision_time(prim, ks.gas.μᵣ, ks.gas.ω)
    fr = chapman_enskog(ks.vs.u, prim, a, A, tau)
    L = norm((f .- fr) ./ prim[1])

    x = [w; sw; tau]
    y = ifelse(L <= 0.005, [1.0, 0.0], [0.0, 1.0])
    return x, y
end

###
# initialize kinetic solver
###

cd(@__DIR__)
D = Dict{Symbol,Any}()
begin
    D[:matter] = "gas"
    D[:case] = "sod"
    D[:space] = "1d1f1v"
    D[:flux] = "kfvs"
    D[:collision] = "bgk"
    D[:nSpecies] = 1
    D[:interpOrder] = 2
    D[:limiter] = "vanleer"
    D[:boundary] = "period"
    D[:cfl] = 0.5
    D[:maxTime] = 5.0

    D[:x0] = 0.0
    D[:x1] = 1.0
    D[:nx] = 100
    D[:pMeshType] = "uniform"
    D[:nxg] = 1

    D[:umin] = -5.0
    D[:umax] = 5.0
    D[:nu] = 100
    D[:vMeshType] = "rectangle"
    D[:nug] = 0

    D[:knudsen] = 0.05
    D[:mach] = 0.0
    D[:prandtl] = 1.0
    D[:inK] = 0.0
    D[:omega] = 0.81
    D[:alphaRef] = 1.0
    D[:omegaRef] = 0.5
end

###
# prepare training set
###

@showprogress for loop = 1:0.2:10
    D[:knudsen] = exp(-Float64(loop))
    ks = SolverSet(D)
    ctr, face = init_fvm(ks, ks.ps, :dynamic_array; structarray = true)
    init_field!(ks, ctr, face)

    t = 0.0
    dt = timestep(ks, ctr, t)
    nt = Int(ks.set.maxTime ÷ dt) + 1
    res = zeros(3)
    for iter = 1:2000#nt
        if iter%13 == 0
            for i = 1:ks.ps.nx
                sw = (ctr[i+1].w .- ctr[i-1].w) / ks.ps.dx[i] / 2.0
                x, y = regime_data(ks, ctr[i].w, ctr[i].prim, sw, ctr[i].f)
                X = hcat(X, x)
                Y = hcat(Y, y)
            end
        end

        reconstruct!(ks, ctr)
        evolve!(ks, ctr, face, dt; mode = Symbol(ks.set.flux), bc = Symbol(ks.set.boundary))
        KitBase.update!(ks, ctr, face, dt, res; coll = Symbol(ks.set.collision), bc = Symbol(ks.set.boundary))

        t += dt
        #if t > ks.set.maxTime || maximum(res) < 5.e-7
        #    break
        #end
    end
end

plot_line(ks, ctr)

function split_dataset(X, Y, ratio = 9::Integer)
    idx2 = rand(size(X, 2) ÷ ratio) * size(X, 2) .|> floor .|> Int
    sort!(idx2)
    unique!(idx2)
    idx1 = setdiff(collect(1:size(X, 2)), idx2)

    x_train = [X[:, j] for j in idx1]
    x_test = [X[:, j] for j in idx2]
    y_train = [Y[:, j] for j in idx1]
    y_test = [Y[:, j] for j in idx2]

    X_train = zeros(eltype(X), size(X, 1), length(x_train))
    Y_train = zeros(eltype(X), size(Y, 1), length(x_train))
    X_test = zeros(eltype(X), size(X, 1), length(x_test))
    Y_test = zeros(eltype(X), size(Y, 1), length(x_test))
    for j in axes(X_train, 2)
        X_train[:, j] .= x_train[j]
        Y_train[:, j] .= y_train[j]
    end
    for j in axes(X_test, 2)
        X_test[:, j] .= x_test[j]
        Y_test[:, j] .= y_test[j]
    end

    return X_train, Y_train, X_test, Y_test
end

X1, Y1, X2, Y2 = split_dataset(X, Y)

@save "data.jld2" X1 Y1 X2 Y2

###
# test
###
#=
i = 50
sw = (ctr[i+1].w .- ctr[i-1].w) / ks.ps.dx[i] / 2.0
tau = vhs_collision_time(ctr[i].prim, ks.gas.μᵣ, ks.gas.ω)
x, y = regime_data(ks, ctr[i].w, ctr[i].prim, sw, ctr[i].f)

Mu, Mxi, _, _1 = gauss_moments(ctr[i].prim, ks.gas.K)
a = pdf_slope(ctr[i].prim, sw, ks.gas.K)
sw = -ctr[i].prim[1] .* moments_conserve_slope(a, Mu, Mxi, 1)
A = pdf_slope(ctr[i].prim, sw, ks.gas.K)
fr = chapman_enskog(ks.vs.u, ctr[i].prim, a, A, tau)

plot(ks.vs.u, ctr[i].f)
plot!(ks.vs.u, fr, line=:dash)
=#

max_stars_repo_name: CSMMLab/Flowmachine

max_issues_repo_name: CSMMLab/Flowmachine

max_forks_repo_name: CSMMLab/Flowmachine


text:
function binary_string_to_decimal(binstring)
    total = 0
    num_bits = length(binstring)
    for i in 1:num_bits
        if binstring[i] == '1'
            total += 2^(num_bits-i)
        end
    end
    total
end

# given an array of binary strings, count whether there are
# more '1's or '0's (or equal) at a given position in the string
function most_common_bit_per_position(input_array, position)
    # return '1' if more ones, or '0' if more zeros, or 'e' if equal
    zero_count = 0
    one_count = 0
    for row in input_array
        if row[position] == '1'
            one_count += 1
        else
            zero_count += 1
        end
    end
    if one_count > zero_count
        return '1'
    elseif zero_count > one_count
        return '0'
    else
        return 'e'
    end
end

# see if we keep a given row, based on conditions for oxygen or co2 filters
# by comparing a bit with the most common bit in that position.
function filter_row(bit, most_common, mode)
    # mode can be 'oxygen' or 'CO2'
    if mode == "oxygen"
        if (bit == most_common) ||
            ((most_common == 'e') && (bit == '1'))
            return true
        else
            return false
        end
    else # CO2
        if most_common == 'e'
            if bit == '0'
                return true
            end
        elseif (bit != most_common)
            return true
        end
    end
    return false
end

# recursive function to filter the list of binary digits down
function filter_by_most_common_bit(input_array, position, mode)
    # mode can be 'oxygen' or 'CO2'
    # position is where we are in the binary string,
    # e.g. for position=4 we are looking at the 4th bit.
    println("Checking position ",position, " for ",mode, ", have ",size(input_array), " rows")

    # loop over input_array, and put selected rows into filtered_array
    filtered_array = []
    for row in input_array
        most_common =  most_common_bit_per_position(input_array, position)
        if filter_row(row[position], most_common, mode)
            push!(filtered_array, row)
        end
    end
    # check if we can break out of the recursion loop
    if (size(filtered_array)[1] == 1)
        return filtered_array
    end
    # if not, move to next position
    filter_by_most_common_bit(filtered_array, position+1, mode)
end

# read the input data into an array (strips newlines for us)
rows = readlines("input.txt")

# get the oxygen and CO2 ratings
# the recursive 'filter_by_most_common_bit' function
# should return a one-element array
oxygen_binary = filter_by_most_common_bit(rows, 1, "oxygen")
co2_binary = filter_by_most_common_bit(rows, 1, "CO2")

if size(oxygen_binary)[1] != 1
    println("Something went wrong filtering for oxygen! ",size(oxygen_binary))
end
if size(co2_binary)[1] != 1
    println("Something went wrong filtering for CO2! ",size(co2_binary))
end

println("oxygen ", oxygen_binary)
println("CO2 ", co2_binary)

# convert to decimal
oxygen_dec = binary_string_to_decimal(oxygen_binary[1])
co2_dec = binary_string_to_decimal(co2_binary[1])
println("oxygen (decimal) ",oxygen_dec)
println("CO2 (decimal) ",co2_dec)

# multiply
println("total life support ",oxygen_dec * co2_dec)


text:
# FCA.jl

"""
Module providing the free componant analysis (FCA) algorithm.
"""
module FCA

# Package requried
using Optim # for manifold optimization
using LinearAlgebra
using GenericLinearAlgebra
using Statistics
using Combinatorics

import StatsBase: sample # for sample in embed.jl

# export functions
export  freecf, # main function

        # loss function
        κ₄,
        neg_abs_sum_free_kurt,
        free_ent,
        sum_free_ent,

        # gradient function
        grad_neg_abs_sum_free_kurt,
        grad_sum_free_ent,

        # freely center and whiten
        mat_center,
        free_whiten,

        # embed
        mat_embed,

        # icf
        icf,
        kurt,
        neg_abs_sum_kurt,
        grad_neg_abs_sum_kurt,
        negent,
        ent_sum,
        grad_ent_sum,

        # errpd
        err_pd

# include sources
include("loss.jl")
include("gradient.jl")
include("fcf.jl")
include("free_whiten.jl")
include("embed.jl")
include("icf.jl")
include("orthopt.jl")
include("errpd.jl")

end


text:
# These definitions are necessary to ensure that the right functions are called for
# different input types. In particular, the besselj for non-integer orders is not
# available for DoubleFloat, so we use the one from ArbNumerics.

# define operations on ArbFloat
HypergeometricFunctions.logabsgamma(x::ArbNumerics.ArbFloat) =
	ArbNumerics.lgamma(abs(x)), sign(gamma(x))

lgamma(x::ArbNumerics.ArbFloat) = ArbNumerics.lgamma(x)
lgamma(x::DoubleFloats.DoubleFloat) = DoubleFloats.lgamma(x)
lgamma(x::Number) = SpecialFunctions.lgamma(x)

J(ν::ArbNumerics.ArbFloat, x) = ArbNumerics.besselj(ν, x)
J(ν, x::ArbNumerics.ArbFloat) = ArbNumerics.besselj(ν, x)
J(ν::ArbNumerics.ArbFloat, x::ArbNumerics.ArbFloat) = ArbNumerics.besselj(ν, x)
J(ν::Float64, x::Float64) = SpecialFunctions.besselj(ν, x)

# doublefloat is missing half-integer bessel
J(ν::DoubleFloats.DoubleFloat, x::DoubleFloats.DoubleFloat) =
    DoubleFloats.Double64(
        (ArbNumerics.besselj(ArbNumerics.ArbFloat(ν), ArbNumerics.ArbFloat(x))))
J(ν::DoubleFloats.DoubleFloat, x::Complex{DoubleFloats.DoubleFloat}) =
    DoubleFloats.Double64(
        (ArbNumerics.besselj(ArbNumerics.ArbFloat(ν), ArbNumerics.ArbFloat(x))))


text:
# Manchester Code

# Author : Irreq

struct Manchester
    main::String
    bit_map::Dict
end

options = Manchester(
"normal",

# bit map for functions
Dict(
    "normal"=>Dict(
    "0"=>"01",
    "1"=>"10"
    ),

    "invert"=>Dict(
    "0"=>"10",
    "1"=>"01"
    ),

    "differential"=>Dict(
    "0"=>"01"

    ),
    "decode"=>Dict(
    "01"=>"0",
    "10"=>"1"
    )
)
)

struct Encoded
    data::String
    method::String
end

struct Decoded
    data::String
    method::String
end

function _parser(data, n)

     return data[1:n], data[n+1:end]
end

function _findall(pattern,string::AbstractString)
    toReturn = UnitRange{Int64}[]
    s = 1
    while true

        range = findnext(pattern,string,s)
        if range == nothing
             break
        else
            push!(toReturn, range)
            s = first(range)+1
        end
    end
    return toReturn
end

function encode(data, encoding=options.main; args=nothing)

    method = options.bit_map[encoding]

    encoded = []

    while length(data) > 0

        bit, data = _parser(data, 1)

        value = method[bit]

        push!(encoded, value)

    end

    encoded = join(encoded)

    result = Encoded(encoded, encoding)

    return result
end

function decode(data, encoding=options.main; args=nothing)

    errors = Dict([i=>_findall(i, data) for (i) = ["000","111"]])

    method = options.bit_map[encoding]

    bits = []

    while length(data) > 0

        bit, data = _parser(data, 2)

        push!(bits, bit)
    end

    decoded = []

    parity = []

    for (i) = bits
        if length(i)%2!=0
            push!(parity, nothing)
            break
        end
        if i[1] == i[2]
            push!(parity, false)

        else
            push!(parity, true)

        end

    end

    for (i) = bits

        bit = ""
        if (encoding) == "invert"
            if (i) == "01"
                bit = "1"
            end
            if (i) == "10"
                bit = "0"
            end

        elseif (encoding) == "normal"
            if (i) == "01"
                bit = "0"
            end
            if (i) == "10"
                bit = "1"
            end
        else
            bit = ""
        end

        push!(decoded, bit)
    end

    decoded = join(decoded)

    result = Encoded(decoded, encoding)
end

max_stars_repo_name: Irreq/manchester-coding

max_issues_repo_name: Irreq/manchester-coding

max_forks_repo_name: Irreq/manchester-coding


text:
using Distributed
addprocs(16)
@everywhere include("../src/get_sens_lorenz.jl")
si = LinRange(27.5,31.5,10)
s = zeros(3,size(si)[1])
s[2,:] .= si
s[1,:] .= 10.
s[3,:] .= 8.0/3
get_sens(s)

max_stars_repo_name: nishaChandramoorthy/linearResponse

max_issues_repo_name: nishaChandramoorthy/linearResponse

max_forks_repo_name: nishaChandramoorthy/linearResponse


text:
"""
    ImageConfig

    ImageConfig(init; kw...)

Common configuration component for all [`ImageOutput`](@ref).

# Keywords

- `init` output init object, used to generate other arguments automatically.
- `minval`: Minimum value in the grid(s) to normalise for conversion to an RGB pixel.
    Number or `Tuple` for multiple grids.
- `maxval`: Maximum value in the grid(s) to normalise for conversion to an RGB pixel.
    Number or `Tuple` for multiple grids.
- `font`: `String` name of font to search for. A default will be guessed.
- `text`: `TextConfig()` or `nothing` for no text. Default is `TextConfig(; font=font)`.
- `scheme`: ColorSchemes.jl scheme, or `Greyscale()`. ObjectScheme() by default.
- `imagegen`: [`ImageGenerator`](@ref) like [`Image`](@ref) or [`Layout`](@ref) Will
    be detected automatically
"""
struct ImageConfig{P,Min,Max,IB,TC}
    imagegen::P
    minval::Min
    maxval::Max
    imagebuffer::IB
    textconfig::TC
end
function ImageConfig(init;
    font=autofont(), text=TextConfig(; font=font), textconfig=text,
    scheme=ObjectScheme(), imagegen=autoimagegen(init, scheme),
    minval=nothing, maxval=nothing, kw...
)
    imagebuffer = _allocimage(imagegen, init)
    ImageConfig(imagegen, minval, maxval, imagebuffer, textconfig)
end

imagegen(ic::ImageConfig) = ic.imagegen
minval(ic::ImageConfig) = ic.minval
maxval(ic::ImageConfig) = ic.maxval
imagebuffer(ic::ImageConfig) = ic.imagebuffer
textconfig(ic::ImageConfig) = ic.textconfig

"""
    ImageOutput <: GraphicOutput

Abstract supertype for Graphic outputs that display the simulation frames as RGB images.

`ImageOutput`s must have [`Extent`](@ref), [`GraphicConfig`](@ref)
and [`ImageConfig`](@ref) components, and define a [`showimage`](@ref) method.

See [`GifOutput`](@ref) for an example.

Although the majority of the code is maintained here to enable sharing
and reuse, most `ImageOutput`s are not provided in DynamicGrids.jl to avoid
heavy dependencies on graphics libraries. See
[DynamicGridsGtk.jl](https://github.com/cesaraustralia/DynamicGridsGtk.jl)
and [DynamicGridsInteract.jl](https://github.com/cesaraustralia/DynamicGridsInteract.jl)
for implementations.
"""
abstract type ImageOutput{T,F} <: GraphicOutput{T,F} end

# Generic `ImageOutput` constructor that construct an `ImageOutput` from another `Output`.
function (::Type{F})(o::T;
    frames=frames(o), extent=extent(o), graphicconfig=graphicconfig(o),
    imageconfig=imageconfig(o), textconfig=textconfig(o), kw...
) where F <: ImageOutput where T <: Output
    F(;
        frames=frames, running=false, extent=extent, graphicconfig=graphicconfig,
        imageconfig=imageconfig, textconfig=textconfig, kw...
    )
end

# Generic `ImageOutput` constructor. Converts an init `AbstractArray` or `NamedTuple`
# to a vector of `AbstractArray`s, uses `kw` to constructs required
# [`Extent`](@ref), [`GraphicConfig`](@ref) and [`ImageConfig`](@ref) objects unless
# they are specifically passed in using `extent`, `graphicconfig`, `imageconfig`.

# All other keyword arguments are passed to these constructors.
# Unused or mis-spelled keyword arguments are ignored.
function (::Type{T})(init::Union{NamedTuple,AbstractMatrix};
    extent=nothing, graphicconfig=nothing, imageconfig=nothing, kw...
) where T <: ImageOutput
    extent = extent isa Nothing ? Extent(; init=init, kw...) : extent
    graphicconfig = graphicconfig isa Nothing ? GraphicConfig(; kw...) : extent
    imageconfig = imageconfig isa Nothing ? ImageConfig(init; kw...) : imageconfig
    T(;
        frames=[deepcopy(init)], running=false, extent=extent,
        graphicconfig=graphicconfig, imageconfig=imageconfig, kw...
    )
end

imageconfig(o::Output) = ImageConfig(init(o))
imageconfig(o::ImageOutput) = o.imageconfig

imagegen(o::Output) = imagegen(imageconfig(o))
minval(o::Output) = minval(imageconfig(o))
maxval(o::Output) = maxval(imageconfig(o))
imagebuffer(o::Output) = imagebuffer(imageconfig(o))
textconfig(o::Output) = textconfig(imageconfig(o))

showframe(o::ImageOutput, data) = showimage(grid_to_image!(o, data), o, data)
showimage(image, o, data) = showimage(image, o)

# Headless image output
mutable struct NoDisplayImageOutput{T,F<:AbstractVector{T},E,GC,IC} <: ImageOutput{T,F}
    frames::F
    running::Bool
    extent::E
    graphicconfig::GC
    imageconfig::IC
end
function NoDisplayImageOutput(;
    frames, running, extent, graphicconfig, imageconfig, kw...
)
    NoDisplayImageOutput(frames, running, extent, graphicconfig, imageconfig)
end


text:
# cstructs

function ArbFloat{P}(x::Mag) where {P}
    z = ArbFloat{P}()
    ccall(@libarb(arf_set_mag), Cvoid, (Ref{ArbFloat}, Ref{Mag},), z, x)
    return z
end

ArbReal{P}(x::Mag) where {P} = ArbReal{P}(ArbFloat{P}(x))

function ArbComplex{P}(rea::Mag) where {P}
    a = ArbFloat{P}(rea)
    return ArbComplex{P}(a)
end

function ArbComplex{P}(rea::ArbFloat{P}, ima::ArbFloat{P}) where {P}
    realpart = ArbReal{P}(rea)
    imagpart = ArbReal{P}(ima)
    result = ArbComplex{P}(realpart, imagpart)
    return result
end
convert(::Type{ArbComplex{P}}, rea::ArbFloat{P}, ima::ArbFloat{P}) where {P} = ArbComplex{P}(rea, ima)

function ArbComplex{P}(rea::ArbReal{P}, ima::ArbReal{P}) where {P}
    z = ArbComplex{P}()
    ccall(@libarb(acb_set_arb_arb), Cvoid,
          (Ref{ArbComplex}, Ref{ArbReal}, Ref{ArbReal}), z, rea, ima)
    return z
end
convert(::Type{ArbComplex{P}}, rea::ArbReal{P}, ima::ArbReal{P}) where {P} = ArbComplex{P}(rea, ima)

function ArbComplex{P}(rea::Mag, ima::Mag) where {P}
    realpart = ArbReal{P}(rea)
    imagpart = ArbReal{P}(ima)
    result = ArbComplex{P}(realpart, imagpart)
    return result
end
convert(::Type{ArbComplex{P}}, rea::Mag, ima::Mag) where {P} = ArbComplex{P}(rea, ima)

function Mag(x::ArbFloat{P}) where {P}
    z = Mag()
    ccall(@libarb(arf_get_mag), Cvoid, (Ref{Mag}, Ref{ArbFloat}), z, x)
    return z
end

function Mag(x::ArbReal{P}) where {P}
    z = Mag()
    ccall(@libarb(arb_get_mag), Cvoid, (Ref{Mag}, Ref{ArbReal}), z, x)
    return z
end

function Mags(x::ArbComplex{P}) where {P}
    rea = Mag(real(x))
    ima = Mag(imag(x))
    return (rea, ima)
end

convert(::Type{Tuple{Mag,Mag}}, x::ArbComplex{P}) where {P} = Mag(real(x)), Mag(imag(y))
convert(::Type{Tuple{ArbFloat{P}, ArbFloat{P}}}, x::ArbComplex{P}) where {P} = (ArbFloat{P}(real(x)), ArbFloat{P}(imag(x)))

function Mag(x::ArbFloat{P}, ::Type{UpperBound}) where {P}
    z = ArbFloat{P}()
    ccall(@libarb(arf_get_ubound_mag), Cvoid,
          (Ref{ArbFloat}, Ref{Mag}), z, x)
    return z
end

function Mag(x::ArbFloat{P}, ::Type{LowerBound}) where {P}
    z = ArbFloat{P}()
    ccall(@libarb(arf_get_lbound_mag), Cvoid,
          (Ref{ArbFloat}, Ref{Mag}), z, x)
    return z
end

ArbFloat{P}(x::ArbReal{P}) where {P} = ArbFloat{P}(midpoint(x), UpperBound)

function ArbFloat{P}(x::ArbReal{P}, ::Type{UpperBound}) where {P}
    z = ArbFloat{P}()
    ccall(@libarb(arb_get_ubound_arf), Cvoid,
          (Ref{ArbFloat}, Ref{ArbReal}, Clong), z, x, P)
    return z
end
convert(::Type{ArbFloat{P}}, x::ArbReal{P}, ::Type{UpperBound}) where {P} = ArbFloat{P}(x)

function ArbFloat{P}(x::ArbReal{P}, ::Type{LowerBound}) where {P}
    z = ArbFloat{P}()
    ccall(@libarb(arb_get_lbound_arf), Cvoid,
          (Ref{ArbFloat}, Ref{ArbReal}, Clong), z, x, P)
    return z
end
convert(::Type{ArbFloat{P}}, x::ArbReal{P}, ::Type{LowerBound}) where {P} = ArbFloat{P}(x)

ArbFloat{P}(x::ArbComplex{P}) where {P} = ArbFloat{P}(x, UpperBound)

function ArbFloat{P}(x::ArbComplex{P}, ::Type{UpperBound}, ) where {P}
    z = ArbFloat{P}()
    ccall(@libarb(arb_get_ubound_arf), Cvoid,
          (Ref{ArbFloat}, Ref{ArbComplex}, Clong), z, x, P)
    return z
end

function ArbFloat{P}(x::ArbComplex{P}, ::Type{LowerBound}) where {P}
    z = ArbFloat{P}()
    ccall(@libarb(arb_get_lbound_arf), Cvoid,
          (Ref{ArbFloat}, Ref{ArbComplex}, Clong), z, x, P)
    return z
end

ArbFloat{P}(x::ArbComplex{P}, ::Type{Radius}) where {P} = ArbFloat{P}(x, UpperBound, Radius)
convert(::Type{ArbFloat{P}}, x::ArbComplex{P}, ::Type{Radius}) where {P} = ArbFloat{P}(x, Radius)

function ArbFloat{P}(x::ArbComplex{P}, ::Type{UpperBound}, ::Type{Radius}) where {P}
    z = ArbFloat{P}()
    ccall(@libarb(arb_get_rad_ubound_arf), Cvoid,
          (Ref{ArbFloat}, Ref{ArbComplex}, Clong), z, x, P)
    return z
end

function ArbFloat{P}(x::ArbComplex{P}, ::Type{LowerBound}, ::Type{Radius}) where {P}
    z = ArbFloat{P}()
    ccall(@libarb(arb_get_rad_lbound_arf), Cvoid,
          (Ref{ArbFloat}, Ref{ArbComplex}, Clong), z, x, P)
    return z
end

Mag(x::ArbComplex{P}) where {P} = Mag(x, UpperBound)

function Mag(x::ArbComplex{P}, ::Type{UpperBound}) where {P}
    z = Mag()
    ccall(@libarb(acb_get_mag), Cvoid,
          (Ref{Mag}, Ref{ArbComplex}, Clong), z, x, P)
    return z
end

function Mag(x::ArbComplex{P}, ::Type{LowerBound}) where {P}
    z = Mag()
    ccall(@libarb(acb_get_mag_lower), Cvoid,
          (Ref{Mag}, Ref{ArbReal}, Clong), z, x, P)
    return z
end


text:
using LinearAlgebra

# Postcompose HuberLoss

f = HuberLoss(1.0, 1.0)
g = Postcompose(f, 2.5)
h = HuberLoss(1.0, 2.5)

x = randn(10)
x = 0.5*x/norm(x)

gx = call_test(g, x)
hx = call_test(h, x)

@test abs(gx-hx)/(1+abs(gx)) <= 1e-12

grad_gx, gx1 = gradient_test(g, x)
grad_hx, hx1 = gradient_test(h, x)

@test abs(gx1-gx)/(1+abs(gx)) <= 1e-12
@test abs(hx1-hx)/(1+abs(hx)) <= 1e-12
@test norm(grad_gx-grad_hx, Inf)/(1+norm(grad_gx, Inf)) <= 1e-12

yg, gyg = prox_test(g, x, 1.3)
yh, hyh = prox_test(h, x, 1.3)

@test abs(gyg-hyh)/(1+abs(gyg)) <= 1e-12
@test norm(yg-yh, Inf)/(1+norm(yg, Inf)) <= 1e-12

x = randn(10)
x = 1.2*x/norm(x)

gx = call_test(g, x)
hx = call_test(h, x)

@test abs(gx-hx)/(1+abs(gx)) <= 1e-12

grad_gx, gx1 = gradient_test(g, x)
grad_hx, hx1 = gradient_test(h, x)

@test abs(gx1-gx)/(1+abs(gx)) <= 1e-12
@test abs(hx1-hx)/(1+abs(hx)) <= 1e-12
@test norm(grad_gx-grad_hx, Inf)/(1+norm(grad_gx, Inf)) <= 1e-12

yg, gyg = prox_test(g, x, 1.3)
yh, hyh = prox_test(h, x, 1.3)

@test abs(gyg-hyh)/(1+abs(gyg)) <= 1e-12
@test norm(yg-yh, Inf)/(1+norm(yg, Inf)) <= 1e-12


text:
module HJBSolverTests
using FactCheck
using HJBSolver
using Base.Test

include("mertontest.jl")
include("boundarytest.jl")
include("test2d.jl")

FactCheck.exitstatus()
end

max_stars_repo_name: anriseth/HJBSolver

max_issues_repo_name: anriseth/HJBSolver

max_forks_repo_name: anriseth/HJBSolver


text:
struct MatrixLogNormal{T<:Real,M<:AbstractMatrix,C<:AbstractPDMat} <:
    Distribution{Matrixvariate,Continuous}
    normal::MatrixNormal{T, M, C}

    function MatrixLogNormal(normal::MatrixNormal{T, M, C}) where {T, M, C}
        new{T, M, C}(normal)
    end
end

function MatrixLogNormal(M, U, V)
    MatrixLogNormal(MatrixNormal(M, U, V))
end

rand(MLN::MatrixLogNormal) = exp.(rand(MLN.normal))

function location(::Type{D}, s::Symbol,
                  M::AbstractMatrix,
                  U::AbstractMatrix, V::AbstractMatrix) where {D<:MatrixLogNormal}
    loc = similar(M)
    @inbounds for jdx in size(loc, 2)
        M[:, jdx] .= location(MvLogNormal, s,
                              M[:, jdx],
                              diagm(V[jdx, jdx] .* diag(U)))
    end
    loc
end

function location(::Type{D}, s::Symbol,
                  M::AbstractMatrix,
                  U::AbstractPDMat, V::AbstractPDMat) where {D<:MatrixLogNormal}
    location(D, s, M, U.mat, V.mat)
end


text:
# Ennui and Endurance (Panda Magazine Jan. 2017)
@testset "ennui and endurance" begin
    cluster = Collective.best_cluster(corpus, Collective.wordlist(IOBuffer("""
        pule
        ramata
        eta
        mien
        put
        fees
        tiy
        buy
        dale
        oro
        jay
        organa
        pia
        freeing
        car
        pried
        """)))
    @test sort(cluster.words) == ["car", "dale", "eta", "fees", "freeing", "mien", "organa", "pried"]
    @test cluster.feature.description == "has a transaddition with letter 'z'"
end


text:
export sumdropdims

"""
    sumdropdims(arr; dims)

Alias for `dropdims(sum(arr, dims=dims), dims=dims)`
"""
function sumdropdims(arr::AbstractArray; dims)
    dropdims(sum(arr, dims=dims), dims=dims)
end

"""
    sumdropdims(f, arr; dims)

Alias for `dropdims(sum(f, arr, dims=dims), dims=dims)`
"""
function sumdropdims(f, arr::AbstractArray; dims)
    dropdims(sum(f, arr, dims=dims), dims=dims)
end


text:
immutable ExponentialIntercept{T} <: Model
    β₀::T
    a::T
    δ::T
end

function Base.show(io::IO, ::Type{ExponentialIntercept})
    return print(io, "Intercept Exponential")
end

function Base.show(io::IO, m::ExponentialIntercept)
    @printf(
        io,
        "Intercept Exponential (β₀ = %.3f, a = %.3f, δ = %.3f)",
        m.β₀,
        m.a,
        m.δ,
    )
end

function initial_parameters(::Type{ExponentialIntercept})
    return Float64[0.0, log(1.0), logit(0.5)]
end

canonical_parameters(m::ExponentialIntercept) = (m.β₀, m.a, m.δ)

function ExponentialIntercept(parameters::Vector)
    β₀ = parameters[1]
    a = exp(parameters[2])
    δ = invlogit(parameters[3])
    return ExponentialIntercept(β₀, a, δ)
end

@inline function predict(
    m::ExponentialIntercept,
    x2::Real,
    t2::Real,
    x1::Real,
    t1::Real,
)
    β₀, a, δ = m.β₀, m.a, m.δ
    u1 = x1 * δ^t1
    u2 = x2 * δ^t2
    z = β₀ + a * (u2 - u1)
    p = invlogit(z)
    return p
end

@inline function gradient_component(
    m::ExponentialIntercept,
    x2::Real,
    t2::Real,
    x1::Real,
    t1::Real,
)
    β₀, a, δ = m.β₀, m.a, m.δ
    u2 = x2 * δ^t2
    u1 = x1 * δ^t1
    gr = (
        1.0,
        (u2 - u1) * a,
        a * ((x2 * t2 * δ^t2 - x1 * t1 * δ^t1) / δ) * δ * (1 - δ),
    )
    return gr
end


text:
using Thrift
import Thrift.process, Thrift.meta

include("gen-jl/floatops/floatops_constants.jl");
include("gen-jl/floatops/floatops_types.jl");
include("gen-jl/floatops/FloatCalc.jl");

include("gen-jl/arithmetic/arithmetic_constants.jl");
include("gen-jl/arithmetic/arithmetic_types.jl");
include("gen-jl/arithmetic/Calc.jl");

function calcclnt(niter::Int)
    clnt_transport = TSocket(19999)
    #proto = TCompactProtocol(clnt_transport)
    proto = TBinaryProtocol(clnt_transport)

    clnt = CalcClient(proto)

    println("opening connection. $(clnt_transport)")
    open(clnt_transport)
    println("opened connection. $(clnt_transport)")
    nops = length(OPS)
    for idx in 1:niter
        try
            op = string(OPS[round(Int, (nops-1)*rand())+1])
            p1 = round(Int32, 100*rand())
            p2 = round(Int32, 100*rand())
            res = calculate(clnt, op, p1, p2)
            println("$op ( $p1, $p2 ) = $res")
        catch ex
            !isa(ex, InvalidOperation) && rethrow()
            println(ex.oper)
        end
    end
    close(clnt_transport)
    println("closed connection. $(clnt_transport)")
end

function floatcalcclnt(niter::Int)
    clnt_transport = TSocket(19999)
    #proto = TCompactProtocol(clnt_transport)
    proto = TBinaryProtocol(clnt_transport)

    clnt = CalcClient(proto)

    println("opening connection. $(clnt_transport)")
    open(clnt_transport)
    println("opened connection. $(clnt_transport)")
    nops = length(FLOAT_OPS)
    for idx in 1:niter
        try
            m = (idx > 3) ? 50 : 8
            op = string(FLOAT_OPS[round(Int, (nops-1)*rand())+1])
            p1 = (2^m) * randn()
            p2 = (2^m) * randn()
            res = float_calculate(clnt, op, p1, p2)
            println("$op ( $p1, $p2 ) = $res")
        catch ex
            !isa(ex, InvalidFloatOperation) && rethrow()
            println(ex.oper)
        end
    end
    close(clnt_transport)
    println("closed connection. $(clnt_transport)")
end

for idx in 1:5
    calcclnt(10)
    floatcalcclnt(10)
end


text:
import Documenter

Documenter.makedocs(
    modules = [IntervalOptimisation],
    format = :html,
    sitename = "IntervalOptimisation.jl",
    root = joinpath(dirname(dirname(@__FILE__)), "docs"),
    pages = Any["Home" => "index.md"],
    strict = true,
    linkcheck = true,
    checkdocs = :exports,
    authors = "David Sanders"
)

Documenter.deploydocs(
    repo = "github.com/dpsanders/IntervalOptimisation.jl.git",
    target = "build",
    deps = nothing,
    make = nothing
)


text:
"""
```
t = find_threshold(Otsu(),  histogram, edges)
```

Under the assumption that the histogram is bimodal the threshold is
set so that the resultant inter-class variance is maximal.

# Output

Returns a real number `t` that specifies the threshold.

# Arguments

The function arguments are described in more detail below.

##  `histogram`

An `AbstractArray` storing the frequency distribution.

##  `edges`

An `AbstractRange` specifying how the intervals for the frequency distribution
are divided.

# Example

Compute the threshold for the "camerman" image in the `TestImages` package.

```julia
using TestImages, ImageContrastAdjustment, HistogramThresholding

img = testimage("cameraman")
edges, counts = build_histogram(img,256)
#=
  The `counts` array stores at index 0 the frequencies that were below the
  first bin edge. Since we are seeking a threshold over the interval
  partitioned by `edges` we need to discard the first bin in `counts`
  so that the dimensions of `edges` and `counts` match.
=#
t = find_threshold(Otsu(),counts[1:end], edges)
```

# Reference

1. Nobuyuki Otsu (1979). "A threshold selection method from gray-level histograms". *IEEE Trans. Sys., Man., Cyber.* 9 (1): 62–66. [doi:10.1109/TSMC.1979.4310076](http://dx.doi.org/doi:10.1109/TSMC.1979.4310076)
"""
function find_threshold(algorithm::Otsu,  histogram::AbstractArray, edges::AbstractRange)
  N = sum(histogram)
  pdf = histogram / N
  histogram_indices = first(axes(pdf))
  first_bin = first(histogram_indices)
  last_bin = last(histogram_indices)
  zeroth_cummulative_moment = cumsum(pdf)
  first_cummulative_moment = cumsum((first_bin:last_bin) .* pdf)
  μ_T = first_cummulative_moment[end]
  σ²_T = sum( ((first_bin:last_bin) .- μ_T).^2  .* pdf )
  maxval = zero(eltype(first(pdf)))

  # Equation (6) for determining the probability of the first class.
  function ω(k)
    let zeroth_cummulative_moment = zeroth_cummulative_moment
      return zeroth_cummulative_moment[k]
    end
  end

  # Equation (7) for determining the mean of the first class.
  function μ(k)
   let first_cummulative_moment = first_cummulative_moment
     return first_cummulative_moment[k]
   end
  end

  # Equation (18) for determining the between-cass variance.
  function σ²(k)
    let μ_T = μ_T
      return (μ_T*ω(k) - μ(k))^2 / (ω(k) * (1-ω(k)))
    end
  end

  t = first_bin
  for k = first_bin:last_bin - 1
    val = σ²(k)
    if val  > maxval
       maxval = val
       t = k
    end
  end
  return edges[t]
end