text:
include("../src/Snappy.jl")

using Snappy
using Base.Test

@testset "round trip compression tests      " begin

    testfiles = [
        "alice29.txt",
        "asyoulik.txt",
        "html",
        "html_x_4",
        "kppkn.gtb",
        "lcet10.txt",
        "fireworks.jpeg",
        "geo.protodata",
        "paper-100k.pdf",
        "plrabn12.txt",
        "urls.10K",
        "random1.bin",
        "random2.bin",
        "random3.bin",
        "smallrandom1.bin",
    ]
    for file in testfiles
        raw = read("$(Base.source_dir())/testdata/$(file)")
        a = compress(raw)
        b = uncompress(a)
        @test hash(a) != hash(raw) # just make sure it isn't identity
        @test hash(b) != hash(a)
        @test hash(raw) == hash(b)
    end
end

@testset "random generated compression tests" begin

    wordsize = 1 << 4
    dictsize = 1 << 6
    maxwords = 1 << 16

    dictionary = [rand(UInt8, rand(1:wordsize)) for _ in 1:dictsize]
    for i in 1:50
        raw = vcat((dictionary[rand(1:dictsize)] for _ in 1:rand(1:maxwords))...)
        a = compress(raw)
        b = uncompress(a)
        @test hash(a) != hash(raw)
        @test hash(b) != hash(a)
        @test hash(raw) == hash(b)
    end
    # run some of the tests using the string wrapper
    for i in 1:50
        raw = vcat((dictionary[rand(1:dictsize)] for _ in 1:rand(1:maxwords))...)
        a = compress(String(raw))
        b = uncompress(a)
        @test hash(a) != hash(raw)
        @test hash(b) != hash(a)
        @test hash(raw) == hash(b)
    end

end

@testset "corrupted data tests              " begin

    src = "making sure we don't crash with corrupted input"
    dst = compress(src)

    @test length(dst) > 3

    dst[2] = ~dst[2]
    dst[4] = dst[3]

    @test_throws ErrorException uncompress(dst)

    # This is testing for a security bug - a buffer that decompresses to 100k
    # but we lie in the snappy header and only reserve 0 bytes of memory :)
    src = repeat("A", 100000)
    dst = compress(src)
    dst[1] = dst[2] = dst[3] = dst[4] = 0

    @test_throws ErrorException uncompress(dst)

    dst[1] = dst[2] = dst[3] = 0xff
    # This decodes to about 2 MB; much smaller, but should still fail.
    dst[4] = 0x00

    @test_throws ErrorException uncompress(dst)

    # try reading stuff in from a bad file.
    testfiles = [
        "baddata1.snappy",
        "baddata2.snappy",
        "baddata3.snappy",
    ]
    for file in testfiles
        raw = read("$(Base.source_dir())/testdata/$(file)")
        @test Snappy.length_uncompressed(raw)[1] < (1<<20)
        @test_throws ErrorException uncompress(raw)
    end

    # corrupted varint tests
    raw = [0xf0]
    @test_throws ErrorException Snappy.parse32(raw, 1)
    @test_throws ErrorException uncompress(raw)

    raw = [0x80, 0x80, 0x80, 0x80, 0x80, 0x0a]
    @test_throws ErrorException Snappy.parse32(raw, 1)
    @test_throws ErrorException uncompress(raw)

    raw = [0xfb, 0xff, 0xff, 0xff, 0x7f]
    @test_throws ErrorException Snappy.parse32(raw, 1)
    @test_throws ErrorException uncompress(raw)

    # check for an infinite loop caused by a copy with offset==0
    raw = [0x40, 0x12, 0x00, 0x00]
    #  0x40                 Length
    #  0x12, 0x00, 0x00     Copy with offset==0, length==5
    @test_throws ErrorException uncompress(raw)

    raw = [0x05, 0x12, 0x00, 0x00]
    #  0x05                 Length
    #  0x12, 0x00, 0x00     Copy with offset==0, length==5
    @test_throws ErrorException uncompress(raw)
end

@testset "simple tests                      " begin

    test_strings = map((e)->Vector{UInt8}(e), [
        "",
        "a",
        "ab",
        "abc",
        "aaaaaaa" * repeat("b", 16) * "aaaaa" * "abc",
        "aaaaaaa" * repeat("b", 256) * "aaaaa" * "abc",
        "aaaaaaa" * repeat("b", 2047) * "aaaaa" * "abc",
        "aaaaaaa" * repeat("b", 65536) * "aaaaa" * "abc",
        "abcaaaaaaa" * repeat("b", 65536) * "aaaaa" * "abc",
    ])

    for raw in test_strings
        a = compress(raw)
        b = uncompress(a)
        @test hash(a) != hash(raw)
        @test hash(b) != hash(a)
        @test hash(raw) == hash(b)
    end

    # verify max blowup (lots of four-byte copies)
    raw = reinterpret(UInt8, rand(UInt32, 20000))
    raw = vcat(raw, flipdim(raw, 1))
    a = compress(raw)
    b = uncompress(a)
    @test hash(a) != hash(raw)
    @test hash(b) != hash(a)
    @test hash(raw) == hash(b)

    # encode a range of varint values
    buf = zeros(UInt8, 5)
    for i in 0:30
        offset = Snappy.encode32!(buf, 1, UInt32(1 << i))
        val, offset2 = Snappy.parse32(buf, 1)
        @test val == UInt32(1 << i)
        @test offset == offset2
    end
end

@testset "find_match_length tests           " begin

    function test_find_match_length(a::String, b::String, limit::Int)
        a = Vector{UInt8}(a)
        b = Vector{UInt8}(b)
        c = vcat(a, b)
        return Snappy.find_match_length(c, start(a), endof(a)+1, endof(a)+limit)
    end

    # Hit s1_limit in 64-bit loop, hit s1_limit in single-character loop.
    @test 6 == test_find_match_length("012345", "012345", 6)
    @test 11 == test_find_match_length("01234567abc", "01234567abc", 11)

    # Hit s1_limit in 64-bit loop, find a non-match in single-character loop.
    @test 9 == test_find_match_length("01234567abc", "01234567axc", 9)

    # Same, but edge cases.
    @test 11 == test_find_match_length("01234567abc!", "01234567abc!", 11)
    @test 11 == test_find_match_length("01234567abc!", "01234567abc?", 11)

    # Find non-match at once in first loop.
    @test 0 == test_find_match_length("01234567xxxxxxxx", "?1234567xxxxxxxx", 16)
    @test 1 == test_find_match_length("01234567xxxxxxxx", "0?234567xxxxxxxx", 16)
    @test 4 == test_find_match_length("01234567xxxxxxxx", "01237654xxxxxxxx", 16)
    @test 7 == test_find_match_length("01234567xxxxxxxx", "0123456?xxxxxxxx", 16)

    # Find non-match in first loop after one block.
    @test 8 == test_find_match_length("abcdefgh01234567xxxxxxxx",
                                   "abcdefgh?1234567xxxxxxxx", 24)
    @test 9 == test_find_match_length("abcdefgh01234567xxxxxxxx",
                                   "abcdefgh0?234567xxxxxxxx", 24)
    @test 12 == test_find_match_length("abcdefgh01234567xxxxxxxx",
                                    "abcdefgh01237654xxxxxxxx", 24)
    @test 15 == test_find_match_length("abcdefgh01234567xxxxxxxx",
                                    "abcdefgh0123456?xxxxxxxx", 24)

    # 32-bit version:

    # Short matches.
    @test 0 == test_find_match_length("01234567", "?1234567", 8)
    @test 1 == test_find_match_length("01234567", "0?234567", 8)
    @test 2 == test_find_match_length("01234567", "01?34567", 8)
    @test 3 == test_find_match_length("01234567", "012?4567", 8)
    @test 4 == test_find_match_length("01234567", "0123?567", 8)
    @test 5 == test_find_match_length("01234567", "01234?67", 8)
    @test 6 == test_find_match_length("01234567", "012345?7", 8)
    @test 7 == test_find_match_length("01234567", "0123456?", 8)
    @test 7 == test_find_match_length("01234567", "0123456?", 7)
    @test 7 == test_find_match_length("01234567!", "0123456??", 7)

    # Hit s1_limit in 32-bit loop, hit s1_limit in single-character loop.
    @test 10 == test_find_match_length("xxxxxxabcd", "xxxxxxabcd", 10)
    @test 10 == test_find_match_length("xxxxxxabcd?", "xxxxxxabcd?", 10)

    # NOTE: this test is from the original C++, but I think it's invalid for julia.
    #   the strings it's comparing are only 12 characters - how can it have
    #   a match length of 13? I suspect the test passes in C++ because it
    #   reads 1 past the end of the string and ends up comparing the null
    #   terminators of each, and includes that as part of the "match length".
    @test_broken 13 == test_find_match_length("xxxxxxabcdef", "xxxxxxabcdef", 13)
    # repeating the test with explicitly including a terminator in the strings
    #   will test the same functionality as the above test is supposed to.
    @test 13 == test_find_match_length("xxxxxxabcdef\0", "xxxxxxabcdef\0", 13)

    # Same, but edge cases.
    @test 12 == test_find_match_length("xxxxxx0123abc!", "xxxxxx0123abc!", 12)
    @test 12 == test_find_match_length("xxxxxx0123abc!", "xxxxxx0123abc?", 12)

    # Hit s1_limit in 32-bit loop, find a non-match in single-character loop.
    @test 11 == test_find_match_length("xxxxxx0123abc", "xxxxxx0123axc", 13)

    # Find non-match at once in first loop.
    @test 6 == test_find_match_length("xxxxxx0123xxxxxxxx",
                                   "xxxxxx?123xxxxxxxx", 18)
    @test 7 ==test_find_match_length("xxxxxx0123xxxxxxxx",
                                   "xxxxxx0?23xxxxxxxx", 18)
    @test 8 == test_find_match_length("xxxxxx0123xxxxxxxx",
                                   "xxxxxx0132xxxxxxxx", 18)
    @test 9 == test_find_match_length("xxxxxx0123xxxxxxxx",
                                   "xxxxxx012?xxxxxxxx", 18)

    # Same, but edge cases.
    @test 6 == test_find_match_length("xxxxxx0123", "xxxxxx?123", 10)
    @test 7 == test_find_match_length("xxxxxx0123", "xxxxxx0?23", 10)
    @test 8 == test_find_match_length("xxxxxx0123", "xxxxxx0132", 10)
    @test 9 == test_find_match_length("xxxxxx0123", "xxxxxx012?", 10)

    # Find non-match in first loop after one block.
    @test 10 == test_find_match_length("xxxxxxabcd0123xx",
                                    "xxxxxxabcd?123xx", 16)
    @test 11 == test_find_match_length("xxxxxxabcd0123xx",
                                    "xxxxxxabcd0?23xx", 16)
    @test 12 == test_find_match_length("xxxxxxabcd0123xx",
                                    "xxxxxxabcd0132xx", 16)
    @test 13 == test_find_match_length("xxxxxxabcd0123xx",
                                    "xxxxxxabcd012?xx", 16)

    # Same, but edge cases.
    @test 10 == test_find_match_length("xxxxxxabcd0123", "xxxxxxabcd?123", 14)
    @test 11 == test_find_match_length("xxxxxxabcd0123", "xxxxxxabcd0?23", 14)
    @test 12 == test_find_match_length("xxxxxxabcd0123", "xxxxxxabcd0132", 14)
    @test 13 == test_find_match_length("xxxxxxabcd0123", "xxxxxxabcd012?", 14)

end


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder

name = "Ncurses"
version = v"6.1"

# Collection of sources required to build Ncurses
sources = [
    "https://ftp.gnu.org/pub/gnu/ncurses/ncurses-$(version.major).$(version.minor).tar.gz" =>
    "aa057eeeb4a14d470101eff4597d5833dcef5965331be3528c08d99cebaa0d17",
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/ncurses-*/

# We need to run the native "tic" program
apk add ncurses

CONFIG_FLAGS=""
if [[ ${target} == x86_64-apple-darwin14 ]]; then
    CONFIG_FLAGS="${CONFIG_FLAGS} --disable-stripping"
elif [[ "${target}" == *-mingw* ]]; then
    CONFIG_FLAGS="--enable-sp-funcs --enable-term-driver"
    export CFLAGS="-lintl -liconv"
fi

./configure --prefix=${prefix} --build=${MACHTYPE} --host=${target} \
    --with-shared \
    --with-normal \
    --without-debug \
    --without-ada \
    --without-cxx-binding \
    --enable-widec \
    --enable-pc-files \
    --disable-rpath \
    --enable-colorfgbg \
    --enable-ext-colors \
    --enable-ext-mouse \
    --enable-warnings \
    --enable-assertions \
    --enable-database \
    ${CONFIG_FLAGS}
make -j${nproc}
make install

# Remove duplicates that don't work on case-insensitive filesystems
rm -f  ${prefix}/share/terminfo/2/2621a
rm -rf ${prefix}/share/terminfo/a
rm -rf ${prefix}/share/terminfo/e
rm -f  ${prefix}/share/terminfo/h/hp2621a
rm -f  ${prefix}/share/terminfo/h/hp70092a
rm -rf ${prefix}/share/terminfo/l
rm -rf ${prefix}/share/terminfo/m
rm -rf ${prefix}/share/terminfo/n
rm -rf ${prefix}/share/terminfo/p
rm -rf ${prefix}/share/terminfo/q
rm -rf ${prefix}/share/terminfo/X

# Install pc files and fool packages looking for non-wide-character ncurses
for lib in ncurses form panel menu; do
    if [[ "${target}" == *-mingw* ]]; then
        # O Windows, Windows, wherefore art thou Windows?
        abiver=6
    fi
    install -Dm644 "misc/${lib}w.pc" "${prefix}/lib/pkgconfig/${lib}w.pc"
    ln -s "${lib}w.pc" "${prefix}/lib/pkgconfig/${lib}.pc"
    ln -s "lib${lib}w${abiver}.${dlext}" "${libdir}/lib${lib}${abiver}.${dlext}"
done
ln -s ncursesw ${prefix}/include/ncurses
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = supported_platforms()

# The products that we will ensure are always built
products = Product[
    LibraryProduct(["libform", "libform6"], :libform),
    LibraryProduct(["libmenu", "libmenu6"], :libmenu),
    LibraryProduct(["libncurses", "libncurses6"], :libncurses),
    LibraryProduct(["libpanel", "libpanel6"], :libpanel),
]

# Dependencies that must be installed before this package can be built
dependencies = [
    "Gettext_jll",
]

# Build the tarballs.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies)

ext: jl

max_stars_repo_name: tbeason/Yggdrasil

max_issues_repo_name: tbeason/Yggdrasil

max_forks_repo_name: tbeason/Yggdrasil


text:
export Connector,getConnector

mutable struct Connector{T,TQ <: Union{Array{T,2},UniformScaling{Int}}, TK <: Union{Array{T,2},SparseMatrixCSC{T,Int}}} <: AbstractMeganetElement{T}
    K::TK
    b::T
    outTimes::Int
    Q::TQ # ???
end

nTheta(this::Connector) = 0
nFeatIn(this::Connector) = size(this.K,2)
nFeatOut(this::Connector) = size(this.K,1)
nDataOut(this::Connector) = ((this.Q==I) ? nFeatOut(this) : size(this.Q,1))
initTheta(this::Connector{T}) where {T <: Number} = zeros(T,0)

function getConnector(TYPE::Type, K; b = zero(TYPE),outTimes=0,Q=I)
	return Connector(K,b,outTimes,Q);
end

function apply(this::Connector{T},theta::Array{T},Y0::Array{T},tmp,doDerivative=true) where {T <: Number}
    nex = div(length(Y0),nFeatIn(this))
    Y0  = reshape(Y0,:,nex)

    if doDerivative
        if isempty(tmp)
            tmp = copy(Y0)
        else
            tmp .= Y0
        end
    end

    Y = this.K*Y0 .+ this.b
    Ydata::Array{T,2} = Array{T, 2}(0, 0) # Temporary fix until we know what type Q is
    if this.outTimes==1
        Ydata = this.Q*Y
    end

    return Ydata, Y, tmp
end

function Jmv(this::Connector{T},dtheta::Array{T},dY::Array{T},theta::Array{T},Y::Array{T},tmp=nothing) where {T <: Number}
    # ??? This doesn't seem to get used?
    nex = div(length(dY),nFeatIn(this))
    dY  = reshape(dY,:,nex)
    dY = this.K*dY
    dYdata::Array{T,2} = Array{T, 2}(0, 0) # Temporary fix until we know what type Q is
    if this.outTimes==1
        dYdata = this.Q*dY
    end

    return dYdata,dY
end

function JTmv(this::Connector{T},Wdata::Array{T},Win::Array{T},theta::Array{T},Y::Array{T},tmp=nothing) where {T <: Number}
    nex = div(length(Y),nFeatIn(this))
    if length(Win)==0
        W = zeros(T,1,1);
    else
        W = reshape(Win,:,nex);
    end

    if length(Wdata)>0
        Wdata = reshape(Wdata,:,nex);
        W = W .+ this.Q'*Wdata
    end

    dtheta = zeros(T,0);
    W   = this.K'*W;

    return dtheta,W
end


text:
using DiffEqProblemLibrary.SDEProblemLibrary: importsdeproblems; importsdeproblems()
using DiffEqProblemLibrary.SDEProblemLibrary: prob_sde_lorenz, oval2ModelExample

oval2ModelExample()

oval2ModelExample(useBigs=true,largeFluctuations=true)


text:
module CellFieldsTests

using Test
using Gridap.Arrays
using Gridap.Geometry
using Gridap.Fields
using Gridap.TensorValues
using FillArrays
using LinearAlgebra
using Gridap.CellData

import Gridap.Fields: ∇

domain = (0,4,0,4)
partition = (2,2)
model = CartesianDiscreteModel(domain,partition)

trian = get_triangulation(model)

qi = [Point(0.0,0.0),]
l = num_cells(trian)
q = Fill(qi,l)

fun(x) = x[1]
∇fun(x) = VectorValue(1,0)
∇(::typeof(fun)) = ∇fun

cf1 = CellField(fun,trian)
r = [[0.0], [2.0], [0.0], [2.0]]
g = Vector{VectorValue{2,Float64}}[[(1, 0)], [(1, 0)], [(1, 0)], [(1, 0)]]
test_cell_field(cf1,q,r,grad=g)

isa(cf1,Geometry.CellField)
@test RefStyle(cf1) == Val{true}()

indices = [3,2,3]
_cf1 = reindex(cf1,indices)
@test length(_cf1) == length(indices)
test_cell_field(_cf1,reindex(q,indices),reindex(r,indices),grad=reindex(g,indices))
isa(_cf1,Geometry.CellField)
@test RefStyle(_cf1) == Val{true}()

cf2 = cf1 + 4
r2 = map( (i) -> i.+4, r)
test_cell_field(cf2,q,r2)

cf3 = 2*cf2
r3 = map( (i) -> 2*i, r2)
test_cell_field(cf3,q,r3)

cf4 = -cf1
r4 = -r
g4 = -g
test_cell_field(cf4,q,r4,grad=g4)

u(x) = VectorValue(x[1],x[2])
∇u(x) = TensorValue(1,0,0,1)
∇(::typeof(u)) = ∇u

cfu = CellField(u,trian)
collect(evaluate(cfu,q))

grad_cfu = outer(∇,cfu)
collect(evaluate(grad_cfu,q))

t_grad_cfu = outer(cfu,∇)
collect(evaluate(t_grad_cfu,q))

div_cfu = ∇⋅cfu
collect(evaluate(div_cfu,q))

curl_cfu = cross(∇,cfu)
collect(evaluate(curl_cfu,q))

epsi_cfu = ε(cfu)
collect(evaluate(epsi_cfu,q))

btrian = BoundaryTriangulation(model)

s = CompressedArray([Point{1,Float64}[(0.25,),(0.75,)]],get_cell_type(btrian))

bcf1 = restrict(cf1,btrian)
@test isa(bcf1,CellField)

nvec = get_normal_vector(btrian)
@test isa(nvec,CellField)

z = 2*bcf1 + nvec
@test isa(z,CellField)

flux1 = nvec⋅∇u
collect(evaluate(flux1,s))

flux2 = ∇u⋅nvec
collect(evaluate(flux2,s))

strian = SkeletonTriangulation(model)

scf1 = restrict(cf1,strian)

@test scf1.left === scf1.inward
@test scf1.right === scf1.outward

@test isa(scf1,SkeletonCellField)

j = jump(scf1)

s = CompressedArray([Point{1,Float64}[(0.25,),(0.75,)]],get_cell_type(strian))
r = fill([0.0,0.0],length(j))
test_cell_field(j,s,r)

m = mean(scf1)
r = [[0.5, 1.5], [2.0, 2.0], [2.5, 3.5], [2.0, 2.0]]
test_cell_field(m,s,r)

j = jump(∇(scf1))
r = fill(VectorValue{2,Float64}[(0,0),(0,0)],length(j))
test_cell_field(j,s,r)

scfu = restrict(cfu,strian)

grad_cfu = jump(outer(∇,scfu))
collect(evaluate(grad_cfu,s))

t_grad_cfu = jump(outer(scfu,∇))
collect(evaluate(t_grad_cfu,s))

div_cfu = mean(∇⋅scfu)
collect(evaluate(div_cfu,s))

curl_cfu = jump(cross(∇,scfu))
collect(evaluate(curl_cfu,s))

epsi_cfu = mean(ε(scfu))
collect(evaluate(epsi_cfu,s))

#using Gridap.Visualization
#
#writevtk(btrian,"btrian",cellfields=["bcf1"=>bcf1,"nvec"=>nvec,"z"=>z])

#using Gridap.Visualization
#
#writevtk(strian,"strian",cellfields=["j"=>j])

end # module


text:
export GaussianMeanVariance, GaussianMeanPrecision, GaussianWeighteMeanPrecision
export MvGaussianMeanCovariance, MvGaussianMeanPrecision, MvGaussianWeightedMeanPrecision
export UnivariateNormalDistributionsFamily, MultivariateNormalDistributionsFamily, NormalDistributionsFamily
export UnivariateGaussianDistributionsFamily, MultivariateGaussianDistributionsFamily, GaussianDistributionsFamily

const GaussianMeanVariance            = NormalMeanVariance
const GaussianMeanPrecision           = NormalMeanPrecision
const GaussianWeighteMeanPrecision    = NormalWeightedMeanPrecision
const MvGaussianMeanCovariance        = MvNormalMeanCovariance
const MvGaussianMeanPrecision         = MvNormalMeanPrecision
const MvGaussianWeightedMeanPrecision = MvNormalWeightedMeanPrecision

const UnivariateNormalDistributionsFamily{T}   = Union{NormalMeanPrecision{T}, NormalMeanVariance{T}, NormalWeightedMeanPrecision{T}}
const MultivariateNormalDistributionsFamily{T} = Union{MvNormalMeanPrecision{T}, MvNormalMeanCovariance{T}, MvNormalWeightedMeanPrecision{T}}
const NormalDistributionsFamily{T}             = Union{UnivariateNormalDistributionsFamily{T}, MultivariateNormalDistributionsFamily{T}}

const UnivariateGaussianDistributionsFamily   = UnivariateNormalDistributionsFamily
const MultivariateGaussianDistributionsFamily = MultivariateNormalDistributionsFamily
const GaussianDistributionsFamily             = NormalDistributionsFamily

import Base: prod, convert

# Variate forms promotion

promote_variate_type(::Type{ Univariate },   ::Type{ F }) where { F <: UnivariateNormalDistributionsFamily }   = F
promote_variate_type(::Type{ Multivariate }, ::Type{ F }) where { F <: MultivariateNormalDistributionsFamily } = F

promote_variate_type(::Type{ Univariate }, ::Type{ <: MvNormalMeanCovariance })        = NormalMeanVariance
promote_variate_type(::Type{ Univariate }, ::Type{ <: MvNormalMeanPrecision })         = NormalMeanPrecision
promote_variate_type(::Type{ Univariate }, ::Type{ <: MvNormalWeightedMeanPrecision }) = NormalWeightedMeanPrecision

promote_variate_type(::Type{ Multivariate }, ::Type{ <: NormalMeanVariance })          = MvNormalMeanCovariance
promote_variate_type(::Type{ Multivariate }, ::Type{ <: NormalMeanPrecision })         = MvNormalMeanPrecision
promote_variate_type(::Type{ Multivariate }, ::Type{ <: NormalWeightedMeanPrecision }) = MvNormalWeightedMeanPrecision

# Conversion to mean - variance parametrisation

function Base.convert(::Type{ NormalMeanVariance{T} }, dist::UnivariateNormalDistributionsFamily) where { T <: Real }
    return NormalMeanVariance(convert(T, mean(dist)), convert(T, var(dist)))
end

function Base.convert(::Type{ MvNormalMeanCovariance{T} }, dist::MultivariateNormalDistributionsFamily) where { T <: Real }
    return convert(MvNormalMeanCovariance{T, AbstractArray{T, 1}}, dist)
end

function Base.convert(::Type{ MvNormalMeanCovariance{T, M} }, dist::MultivariateNormalDistributionsFamily) where { T <: Real, M <: AbstractArray{T} }
    return convert(MvNormalMeanCovariance{T, AbstractArray{T, 1}, AbstractArray{T, 2}}, dist)
end

function Base.convert(::Type{ MvNormalMeanCovariance{T, M, P} }, dist::MultivariateNormalDistributionsFamily) where { T <: Real, M <: AbstractArray{T}, P <: AbstractArray{T} }
    return MvNormalMeanCovariance(convert(M, mean(dist)), convert(P, cov(dist)))
end

function Base.convert(::Type{ NormalMeanVariance }, dist::UnivariateNormalDistributionsFamily{T}) where { T <: Real }
    return convert(NormalMeanVariance{T}, dist)
end

function Base.convert(::Type{ MvNormalMeanCovariance }, dist::MultivariateNormalDistributionsFamily{T}) where { T <: Real }
    return convert(MvNormalMeanCovariance{T}, dist)
end

# Conversion to mean - precision parametrisation

function Base.convert(::Type{ NormalMeanPrecision{T} }, dist::UnivariateNormalDistributionsFamily) where { T <: Real }
    return NormalMeanPrecision(convert(T, mean(dist)), convert(T, precision(dist)))
end

function Base.convert(::Type{ MvNormalMeanPrecision{T} }, dist::MultivariateNormalDistributionsFamily) where { T <: Real }
    return convert(MvNormalMeanPrecision{T, AbstractArray{T, 1}}, dist)
end

function Base.convert(::Type{ MvNormalMeanPrecision{T, M} }, dist::MultivariateNormalDistributionsFamily) where { T <: Real, M <: AbstractArray{T} }
    return convert(MvNormalMeanPrecision{T, AbstractArray{T, 1}, AbstractArray{T, 2}}, dist)
end

function Base.convert(::Type{ MvNormalMeanPrecision{T, M, P} }, dist::MultivariateNormalDistributionsFamily) where { T <: Real, M <: AbstractArray{T}, P <: AbstractArray{T} }
    return MvNormalMeanPrecision(convert(M, mean(dist)), convert(P, precision(dist)))
end

function Base.convert(::Type{ NormalMeanPrecision }, dist::UnivariateNormalDistributionsFamily{T}) where { T <: Real }
    return convert(NormalMeanPrecision{T}, dist)
end

function Base.convert(::Type{ MvNormalMeanPrecision }, dist::MultivariateNormalDistributionsFamily{T}) where { T <: Real }
    return convert(MvNormalMeanPrecision{T}, dist)
end

# Conversion to weighted mean - precision parametrisation

function Base.convert(::Type{ NormalWeightedMeanPrecision{T} }, dist::UnivariateNormalDistributionsFamily) where { T <: Real }
    return NormalWeightedMeanPrecision(convert(T, weightedmean(dist)), convert(T, precision(dist)))
end

function Base.convert(::Type{ MvNormalWeightedMeanPrecision{T} }, dist::MultivariateNormalDistributionsFamily) where { T <: Real }
    return convert(MvNormalWeightedMeanPrecision{T, AbstractArray{T, 1}}, dist)
end

function Base.convert(::Type{ MvNormalWeightedMeanPrecision{T, M} }, dist::MultivariateNormalDistributionsFamily) where { T <: Real, M <: AbstractArray{T} }
    return convert(MvNormalWeightedMeanPrecision{T, AbstractArray{T, 1}, AbstractArray{T, 2}}, dist)
end

function Base.convert(::Type{ MvNormalWeightedMeanPrecision{T, M, P} }, dist::MultivariateNormalDistributionsFamily) where { T <: Real, M <: AbstractArray{T}, P <: AbstractArray{T} }
    return MvNormalWeightedMeanPrecision(convert(M, weightedmean(dist)), convert(P, precision(dist)))
end

function Base.convert(::Type{ NormalWeightedMeanPrecision }, dist::UnivariateNormalDistributionsFamily{T}) where { T <: Real }
    return convert(NormalWeightedMeanPrecision{T}, dist)
end

function Base.convert(::Type{ MvNormalWeightedMeanPrecision }, dist::MultivariateNormalDistributionsFamily{T}) where { T <: Real }
    return convert(MvNormalWeightedMeanPrecision{T}, dist)
end

# Basic prod fallbacks to weighted mean precision and converts first argument back

prod_analytical_rule(::Type{ <: UnivariateNormalDistributionsFamily }, ::Type{ <: UnivariateNormalDistributionsFamily }) = ProdAnalyticalRuleAvailable()

function Base.prod(::ProdAnalytical, left::L, right::R) where { L <: UnivariateNormalDistributionsFamily, R <: UnivariateNormalDistributionsFamily }
    wleft  = convert(NormalWeightedMeanPrecision, left)
    wright = convert(NormalWeightedMeanPrecision, right)
    return prod(ProdAnalytical(), wleft, wright)
end

prod_analytical_rule(::Type{ <: MultivariateNormalDistributionsFamily }, ::Type{ <: MultivariateNormalDistributionsFamily }) = ProdAnalyticalRuleAvailable()

function Base.prod(::ProdAnalytical, left::L, right::R) where { L <: MultivariateNormalDistributionsFamily, R <: MultivariateNormalDistributionsFamily }
    wleft  = convert(MvNormalWeightedMeanPrecision, left)
    wright = convert(MvNormalWeightedMeanPrecision, right)
    return prod(ProdAnalytical(), wleft, wright)
end


text:
## Integrator Dispatches

# Can get rid of an allocation here with a function
# get_tmp_arr(integrator.cache) which gives a pointer to some
# cache array which can be modified.

@inline function _searchsortedfirst(v::AbstractVector, x, lo::Integer, forward::Bool)
    u = oftype(lo, 1)
    lo = lo - u
    hi = length(v) + u
    @inbounds while lo < hi - u
        m = (lo + hi) >>> 1
        @inbounds if (forward && v[m] < x) || (!forward && v[m] > x)
            lo = m
        else
            hi = m
        end
    end
    return hi
end

@inline function _searchsortedlast(v::AbstractVector, x, lo::Integer, forward::Bool)
  u = oftype(lo, 1)
  lo = lo - u
  hi = length(v) + u
  @inbounds while lo < hi - u
      m = (lo + hi) >>> 1
      @inbounds if (forward && v[m] > x) || (!forward && v[m] < x)
          hi = m
      else
          lo = m
      end
  end
  return lo
end

@inline function _ode_addsteps!(integrator,f=integrator.f,always_calc_begin = false,allow_calc_end = true,force_calc_end = false)
  if !(typeof(integrator.cache) <: CompositeCache)
    DiffEqBase.addsteps!(integrator.k,integrator.tprev,integrator.uprev,integrator.u,
                  integrator.dt,f,integrator.p,integrator.cache,
                  always_calc_begin,allow_calc_end,force_calc_end)
  else
    DiffEqBase.addsteps!(integrator.k,integrator.tprev,integrator.uprev,integrator.u,
                  integrator.dt,f,integrator.p,
                  @inbounds(integrator.cache.caches[integrator.cache.current]),
                  always_calc_begin,allow_calc_end,force_calc_end)
  end
end
@inline DiffEqBase.addsteps!(integrator::ODEIntegrator,args...) = _ode_addsteps!(integrator,args...)

@inline function ode_interpolant(Θ,integrator::DiffEqBase.DEIntegrator,idxs,deriv)
  DiffEqBase.addsteps!(integrator)
  if !(typeof(integrator.cache) <: CompositeCache)
    val = ode_interpolant(Θ,integrator.dt,integrator.uprev,integrator.u,integrator.k,integrator.cache,idxs,deriv)
  else
    val = composite_ode_interpolant(Θ,integrator,integrator.cache.caches,integrator.cache.current,idxs,deriv)
  end
  val
end

@generated function composite_ode_interpolant(Θ, integrator, caches::T, current, idxs, deriv) where {T <: Tuple}
  expr = Expr(:block)
  for i in 1:length(T.types)
    push!(expr.args, quote
      if $i == current
        return ode_interpolant(Θ,integrator.dt,integrator.uprev,integrator.u,integrator.k,caches[$i],idxs,deriv)
      end
    end)
  end
  push!(expr.args, quote
    throw("Cache $current is not available. There are only $(length(caches)) caches.")
  end)
  return expr
end

@inline function ode_interpolant!(val,Θ,integrator::DiffEqBase.DEIntegrator,idxs,deriv)
  DiffEqBase.addsteps!(integrator)
  if !(typeof(integrator.cache) <: CompositeCache)
    ode_interpolant!(val,Θ,integrator.dt,integrator.uprev,integrator.u,integrator.k,integrator.cache,idxs,deriv)
  else
    ode_interpolant!(val,Θ,integrator.dt,integrator.uprev,integrator.u,integrator.k,integrator.cache.caches[integrator.cache.current],idxs,deriv)
  end
end

@generated function composite_ode_interpolant!(val, Θ, integrator, caches::T, current, idxs, deriv) where {T <: Tuple}
  expr = Expr(:block)
  for i in 1:length(T.types)
    push!(expr.args, quote
      if $i == current
        return ode_interpolant!(val,Θ,integrator.dt,integrator.uprev,integrator.u,integrator.k,caches[$i],idxs,deriv)
      end
    end)
  end
  push!(expr.args, quote
    throw("Cache $current is not available. There are only $(length(caches)) caches.")
  end)
  return expr
end

@inline function current_interpolant(t::Number,integrator::DiffEqBase.DEIntegrator,idxs,deriv)
  Θ = (t-integrator.tprev)/integrator.dt
  ode_interpolant(Θ,integrator,idxs,deriv)
end

@inline function current_interpolant(t,integrator::DiffEqBase.DEIntegrator,idxs,deriv)
  Θ = (t.-integrator.tprev)./integrator.dt
  [ode_interpolant(ϕ,integrator,idxs,deriv) for ϕ in Θ]
end

@inline function current_interpolant!(val,t::Number,integrator::DiffEqBase.DEIntegrator,idxs,deriv)
  Θ = (t-integrator.tprev)/integrator.dt
  ode_interpolant!(val,Θ,integrator,idxs,deriv)
end

@inline function current_interpolant!(val,t,integrator::DiffEqBase.DEIntegrator,idxs,deriv)
  Θ = (t.-integrator.tprev)./integrator.dt
  [ode_interpolant!(val,ϕ,integrator,idxs,deriv) for ϕ in Θ]
end

@inline function current_extrapolant(t::Number,integrator::DiffEqBase.DEIntegrator,idxs=nothing,deriv=Val{0})
  Θ = (t-integrator.tprev)/(integrator.t-integrator.tprev)
  ode_extrapolant(Θ,integrator,idxs,deriv)
end

@inline function current_extrapolant!(val,t::Number,integrator::DiffEqBase.DEIntegrator,idxs=nothing,deriv=Val{0})
  Θ = (t-integrator.tprev)/(integrator.t-integrator.tprev)
  ode_extrapolant!(val,Θ,integrator,idxs,deriv)
end

@inline function current_extrapolant(t::AbstractArray,integrator::DiffEqBase.DEIntegrator,idxs=nothing,deriv=Val{0})
  Θ = (t.-integrator.tprev)./(integrator.t-integrator.tprev)
  [ode_extrapolant(ϕ,integrator,idxs,deriv) for ϕ in Θ]
end

@inline function current_extrapolant!(val,t,integrator::DiffEqBase.DEIntegrator,idxs=nothing,deriv=Val{0})
  Θ = (t.-integrator.tprev)./(integrator.t-integrator.tprev)
  [ode_extrapolant!(val,ϕ,integrator,idxs,deriv) for ϕ in Θ]
end

@inline function ode_extrapolant!(val,Θ,integrator::DiffEqBase.DEIntegrator,idxs,deriv)
  DiffEqBase.addsteps!(integrator)
  if !(typeof(integrator.cache) <: CompositeCache)
    ode_interpolant!(val,Θ,integrator.t-integrator.tprev,integrator.uprev2,integrator.uprev,integrator.k,integrator.cache,idxs,deriv)
  else
    composite_ode_extrapolant!(val,Θ,integrator,integrator.cache.caches,integrator.cache.current,idxs,deriv)
  end
end

@generated function composite_ode_extrapolant!(val, Θ, integrator, caches::T, current, idxs, deriv) where {T <: Tuple}
  expr = Expr(:block)
  for i in 1:length(T.types)
    push!(expr.args, quote
      if $i == current
        return ode_interpolant!(val,Θ,integrator.t-integrator.tprev,integrator.uprev2,integrator.uprev,integrator.k,caches[$i],idxs,deriv)
      end
    end)
  end
  push!(expr.args, quote
    throw("Cache $current is not available. There are only $(length(caches)) caches.")
  end)
  return expr
end

@inline function ode_extrapolant(Θ,integrator::DiffEqBase.DEIntegrator,idxs,deriv)
  DiffEqBase.addsteps!(integrator)
  if !(typeof(integrator.cache) <: CompositeCache)
    ode_interpolant(Θ,integrator.t-integrator.tprev,integrator.uprev2,integrator.uprev,integrator.k,integrator.cache,idxs,deriv)
  else
    composite_ode_extrapolant(Θ,integrator,integrator.cache.caches,integrator.cache.current,idxs,deriv)
  end
end

@generated function composite_ode_extrapolant(Θ, integrator, caches::T, current, idxs, deriv) where {T <: Tuple}
  expr = Expr(:block)
  for i in 1:length(T.types)
    push!(expr.args, quote
      if $i == current
        return ode_interpolant(Θ,integrator.t-integrator.tprev,integrator.uprev2,integrator.uprev,integrator.k,caches[$i],idxs,deriv)
      end
    end)
  end
  push!(expr.args, quote
    throw("Cache $current is not available. There are only $(length(caches)) caches.")
  end)
  return expr
end

"""
ode_interpolation(tvals,ts,timeseries,ks)

Get the value at tvals where the solution is known at the
times ts (sorted), with values timeseries and derivatives ks
"""
function ode_interpolation(tvals,id,idxs,deriv,p,continuity::Symbol=:left)
  @unpack ts,timeseries,ks,f,cache = id
  @inbounds tdir = sign(ts[end]-ts[1])
  idx = sortperm(tvals,rev=tdir<0)

  # start the search thinking it's ts[1]-ts[2]
  i₋ = 1
  i₊ = 2
  vals = map(idx) do j
    t = tvals[j]

    if continuity === :left
      # we have i₋ = i₊ = 1 if t = ts[1], i₊ = i₋ + 1 = lastindex(ts) if t > ts[end],
      # and otherwise i₋ and i₊ satisfy ts[i₋] < t ≤ ts[i₊]
      i₊ = min(lastindex(ts), _searchsortedfirst(ts,t,i₊,tdir > 0))
      i₋ = i₊ > 1 ? i₊ - 1 : i₊
    else
      # we have i₋ = i₊ - 1 = 1 if t < ts[1], i₊ = i₋ = lastindex(ts) if t = ts[end],
      # and otherwise i₋ and i₊ satisfy ts[i₋] ≤ t < ts[i₊]
      i₋ = max(1, _searchsortedlast(ts,t,i₋,tdir > 0))
      i₊ = i₋ < lastindex(ts) ? i₋ + 1 : i₋
    end

    dt = ts[i₊] - ts[i₋]
    Θ = iszero(dt) ? oneunit(t) / oneunit(dt) : (t-ts[i₋]) / dt

    if typeof(cache) <: (FunctionMapCache) || typeof(cache) <: FunctionMapConstantCache
      return ode_interpolant(Θ,dt,timeseries[i₋],timeseries[i₊],0,cache,idxs,deriv)
    elseif !id.dense
      return linear_interpolant(Θ,dt,timeseries[i₋],timeseries[i₊],idxs,deriv)
    elseif typeof(cache) <: CompositeCache
      DiffEqBase.addsteps!(ks[i₊],ts[i₋],timeseries[i₋],timeseries[i₊],dt,f,p,cache.caches[id.alg_choice[i₊]]) # update the kcurrent
      return ode_interpolant(Θ,dt,timeseries[i₋],timeseries[i₊],ks[i₊],cache.caches[id.alg_choice[i₊]],idxs,deriv)
    else
      DiffEqBase.addsteps!(ks[i₊],ts[i₋],timeseries[i₋],timeseries[i₊],dt,f,p,cache) # update the kcurrent
      return ode_interpolant(Θ,dt,timeseries[i₋],timeseries[i₊],ks[i₊],cache,idxs,deriv)
    end
  end
  invpermute!(vals, idx)
  DiffEqArray(vals, tvals)
end

"""
ode_interpolation(tvals,ts,timeseries,ks)

Get the value at tvals where the solution is known at the
times ts (sorted), with values timeseries and derivatives ks
"""
function ode_interpolation!(vals,tvals,id,idxs,deriv,p,continuity::Symbol=:left)
  @unpack ts,timeseries,ks,f,cache = id
  @inbounds tdir = sign(ts[end]-ts[1])
  idx = sortperm(tvals,rev=tdir<0)

  # start the search thinking it's in ts[1]-ts[2]
  i₋ = 1
  i₊ = 2
  @inbounds for j in idx
    t = tvals[j]

    if continuity === :left
      # we have i₋ = i₊ = 1 if t = ts[1], i₊ = i₋ + 1 = lastindex(ts) if t > ts[end],
      # and otherwise i₋ and i₊ satisfy ts[i₋] < t ≤ ts[i₊]
      i₊ = min(lastindex(ts), _searchsortedfirst(ts,t,i₊,tdir > 0))
      i₋ = i₊ > 1 ? i₊ - 1 : i₊
    else
      # we have i₋ = i₊ - 1 = 1 if t < ts[1], i₊ = i₋ = lastindex(ts) if t = ts[end],
      # and otherwise i₋ and i₊ satisfy ts[i₋] ≤ t < ts[i₊]
      i₋ = max(1, _searchsortedlast(ts,t,i₋,tdir > 0))
      i₊ = i₋ < lastindex(ts) ? i₋ + 1 : i₋
    end

    dt = ts[i₊] - ts[i₋]
    Θ = iszero(dt) ? oneunit(t) / oneunit(dt) : (t-ts[i₋]) / dt

    if typeof(cache) <: (FunctionMapCache) || typeof(cache) <: FunctionMapConstantCache
      if eltype(timeseries) <: AbstractArray
        ode_interpolant!(vals[j],Θ,dt,timeseries[i₋],timeseries[i₊],0,cache,idxs,deriv)
      else
        vals[j] = ode_interpolant(Θ,dt,timeseries[i₋],timeseries[i₊],0,cache,idxs,deriv)
      end
    elseif !id.dense
      if eltype(timeseries) <: AbstractArray
        linear_interpolant!(vals[j],Θ,dt,timeseries[i₋],timeseries[i₊],idxs,deriv)
      else
        vals[j] = linear_interpolant(Θ,dt,timeseries[i₋],timeseries[i₊],idxs,deriv)
      end
    elseif typeof(cache) <: CompositeCache
      DiffEqBase.addsteps!(ks[i₊],ts[i₋],timeseries[i₋],timeseries[i₊],dt,f,p,cache.caches[id.alg_choice[i₊]]) # update the kcurrent
      if eltype(timeseries) <: AbstractArray
        ode_interpolant!(vals[j],Θ,dt,timeseries[i₋],timeseries[i₊],ks[i₊],cache.caches[id.alg_choice[i₊]],idxs,deriv)
      else
        vals[j] = ode_interpolant(Θ,dt,timeseries[i₋],timeseries[i₊],ks[i₊],cache.caches[id.alg_choice[i₊]],idxs,deriv)
      end
    else
      DiffEqBase.addsteps!(ks[i₊],ts[i₋],timeseries[i₋],timeseries[i₊],dt,f,p,cache) # update the kcurrent
      if eltype(vals[j]) <: AbstractArray
        ode_interpolant!(vals[j],Θ,dt,timeseries[i₋],timeseries[i₊],ks[i₊],cache,idxs,deriv)
      else
        vals[j] = ode_interpolant(Θ,dt,timeseries[i₋],timeseries[i₊],ks[i₊],cache,idxs,deriv)
      end
    end
  end

  vals
end

"""
ode_interpolation(tval::Number,ts,timeseries,ks)

Get the value at tval where the solution is known at the
times ts (sorted), with values timeseries and derivatives ks
"""
function ode_interpolation(tval::Number,id,idxs,deriv,p,continuity::Symbol=:left)
  @unpack ts,timeseries,ks,f,cache = id
  @inbounds tdir = sign(ts[end]-ts[1])

  if continuity === :left
    # we have i₋ = i₊ = 1 if tval = ts[1], i₊ = i₋ + 1 = lastindex(ts) if tval > ts[end],
    # and otherwise i₋ and i₊ satisfy ts[i₋] < tval ≤ ts[i₊]
    i₊ = min(lastindex(ts), _searchsortedfirst(ts,tval,2,tdir > 0))
    i₋ = i₊ > 1 ? i₊ - 1 : i₊
  else
    # we have i₋ = i₊ - 1 = 1 if tval < ts[1], i₊ = i₋ = lastindex(ts) if tval = ts[end],
    # and otherwise i₋ and i₊ satisfy ts[i₋] ≤ tval < ts[i₊]
    i₋ = max(1, _searchsortedlast(ts,tval,1,tdir > 0))
    i₊ = i₋ < lastindex(ts) ? i₋ + 1 : i₋
  end

  @inbounds begin
    dt = ts[i₊] - ts[i₋]
    Θ = iszero(dt) ? oneunit(tval) / oneunit(dt) : (tval-ts[i₋]) / dt

    if typeof(cache) <: (FunctionMapCache) || typeof(cache) <: FunctionMapConstantCache
      val = ode_interpolant(Θ,dt,timeseries[i₋],timeseries[i₊],0,cache,idxs,deriv)
    elseif !id.dense
      val = linear_interpolant(Θ,dt,timeseries[i₋],timeseries[i₊],idxs,deriv)
    elseif typeof(cache) <: CompositeCache
      DiffEqBase.addsteps!(ks[i₊],ts[i₋],timeseries[i₋],timeseries[i₊],dt,f,p,cache.caches[id.alg_choice[i₊]]) # update the kcurrent
      val = ode_interpolant(Θ,dt,timeseries[i₋],timeseries[i₊],ks[i₊],cache.caches[id.alg_choice[i₊]],idxs,deriv)
    else
      DiffEqBase.addsteps!(ks[i₊],ts[i₋],timeseries[i₋],timeseries[i₊],dt,f,p,cache) # update the kcurrent
      val = ode_interpolant(Θ,dt,timeseries[i₋],timeseries[i₊],ks[i₊],cache,idxs,deriv)
    end
  end

  val
end

"""
ode_interpolation!(out,tval::Number,ts,timeseries,ks)

Get the value at tval where the solution is known at the
times ts (sorted), with values timeseries and derivatives ks
"""
function ode_interpolation!(out,tval::Number,id,idxs,deriv,p,continuity::Symbol=:left)
  @unpack ts,timeseries,ks,f,cache = id
  @inbounds tdir = sign(ts[end]-ts[1])

  if continuity === :left
    # we have i₋ = i₊ = 1 if tval = ts[1], i₊ = i₋ + 1 = lastindex(ts) if tval > ts[end],
    # and otherwise i₋ and i₊ satisfy ts[i₋] < tval ≤ ts[i₊]
    i₊ = min(lastindex(ts), _searchsortedfirst(ts,tval,2,tdir > 0))
    i₋ = i₊ > 1 ? i₊ - 1 : i₊
  else
    # we have i₋ = i₊ - 1 = 1 if tval < ts[1], i₊ = i₋ = lastindex(ts) if tval = ts[end],
    # and otherwise i₋ and i₊ satisfy ts[i₋] ≤ tval < ts[i₊]
    i₋ = max(1, _searchsortedlast(ts,tval,1,tdir > 0))
    i₊ = i₋ < lastindex(ts) ? i₋ + 1 : i₋
  end

  @inbounds begin
    dt = ts[i₊] - ts[i₋]
    Θ = iszero(dt) ? oneunit(tval) / oneunit(dt) : (tval-ts[i₋]) / dt

    if typeof(cache) <: (FunctionMapCache) || typeof(cache) <: FunctionMapConstantCache
      ode_interpolant!(out,Θ,dt,timeseries[i₋],timeseries[i₊],0,cache,idxs,deriv)
    elseif !id.dense
      linear_interpolant!(out,Θ,dt,timeseries[i₋],timeseries[i₊],idxs,deriv)
    elseif typeof(cache) <: CompositeCache
      DiffEqBase.addsteps!(ks[i₊],ts[i₋],timeseries[i₋],timeseries[i₊],dt,f,p,cache.caches[id.alg_choice[i₊]]) # update the kcurrent
      ode_interpolant!(out,Θ,dt,timeseries[i₋],timeseries[i₊],ks[i₊],cache.caches[id.alg_choice[i₊]],idxs,deriv)
    else
      DiffEqBase.addsteps!(ks[i₊],ts[i₋],timeseries[i₋],timeseries[i₊],dt,f,p,cache) # update the kcurrent
      ode_interpolant!(out,Θ,dt,timeseries[i₋],timeseries[i₊],ks[i₊],cache,idxs,deriv)
    end
  end

  out
end

"""
By default, Hermite interpolant so update the derivative at the two ends
"""
function DiffEqBase.addsteps!(k,t,uprev,u,dt,f,p,cache,always_calc_begin = false,allow_calc_end = true,force_calc_end = false)
  if length(k)<2 || always_calc_begin
    if typeof(cache) <: OrdinaryDiffEqMutableCache
      rtmp = similar(u, eltype(eltype(k)))
      f(rtmp,uprev,p,t)
      copyat_or_push!(k,1,rtmp)
      f(rtmp,u,p,t+dt)
      copyat_or_push!(k,2,rtmp)
    else
      copyat_or_push!(k,1,f(uprev,p,t))
      copyat_or_push!(k,2,f(u,p,t+dt))
    end
  end
  nothing
end

"""
ode_interpolant and ode_interpolant! dispatch
"""
function ode_interpolant(Θ,dt,y₀,y₁,k,cache,idxs,T::Type{Val{TI}}) where TI
  _ode_interpolant(Θ,dt,y₀,y₁,k,cache,idxs,T)
end

function ode_interpolant(Θ,dt,y₀,y₁,k,cache::OrdinaryDiffEqMutableCache,idxs,T::Type{Val{TI}}) where TI
  if typeof(idxs) <: Number || typeof(y₀) <: Union{Number,SArray}
    # typeof(y₀) can be these if saveidxs gives a single value
    _ode_interpolant(Θ,dt,y₀,y₁,k,cache,idxs,T)
  elseif typeof(idxs) <: Nothing
    out = oneunit(Θ) .* y₁
    _ode_interpolant!(out,Θ,dt,y₀,y₁,k,cache,idxs,T)
  else
    out = oneunit(Θ) .* y₁[idxs]
    _ode_interpolant!(out,Θ,dt,y₀,y₁,k,cache,idxs,T)
  end
end

function ode_interpolant!(out,Θ,dt,y₀,y₁,k,cache,idxs,T::Type{Val{TI}}) where TI
  _ode_interpolant!(out,Θ,dt,y₀,y₁,k,cache,idxs,T)
end

##################### Hermite Interpolants

# If no dispatch found, assume Hermite
function _ode_interpolant(Θ,dt,y₀,y₁,k,cache,idxs,T::Type{Val{TI}}) where TI
  hermite_interpolant(Θ,dt,y₀,y₁,k,Val{typeof(cache) <: OrdinaryDiffEqMutableCache},idxs,T)
end

function _ode_interpolant!(out,Θ,dt,y₀,y₁,k,cache,idxs,T::Type{Val{TI}}) where TI
  hermite_interpolant!(out,Θ,dt,y₀,y₁,k,idxs,T)
end

"""
Hairer Norsett Wanner Solving Ordinary Differential Euations I - Nonstiff Problems Page 190

Herimte Interpolation, chosen if no other dispatch for ode_interpolant
"""
@muladd function hermite_interpolant(Θ,dt,y₀,y₁,k,::Type{Val{false}},idxs::Nothing,T::Type{Val{0}}) # Default interpolant is Hermite
  #@.. (1-Θ)*y₀+Θ*y₁+Θ*(Θ-1)*((1-2Θ)*(y₁-y₀)+(Θ-1)*dt*k[1] + Θ*dt*k[2])
  @inbounds (1-Θ)*y₀+Θ*y₁+Θ*(Θ-1)*((1-2Θ)*(y₁-y₀)+(Θ-1)*dt*k[1] + Θ*dt*k[2])
end

@muladd function hermite_interpolant(Θ,dt,y₀,y₁,k,::Type{Val{true}},idxs::Nothing,T::Type{Val{0}}) # Default interpolant is Hermite
  #@.. (1-Θ)*y₀+Θ*y₁+Θ*(Θ-1)*((1-2Θ)*(y₁-y₀)+(Θ-1)*dt*k[1] + Θ*dt*k[2])
  @inbounds @.. (1-Θ)*y₀+Θ*y₁+Θ*(Θ-1)*((1-2Θ)*(y₁-y₀)+(Θ-1)*dt*k[1] + Θ*dt*k[2])
end

@muladd function hermite_interpolant(Θ,dt,y₀,y₁,k,cache,idxs,T::Type{Val{0}}) # Default interpolant is Hermite
  # return @.. (1-Θ)*y₀[idxs]+Θ*y₁[idxs]+Θ*(Θ-1)*((1-2Θ)*(y₁[idxs]-y₀[idxs])+(Θ-1)*dt*k[1][idxs] + Θ*dt*k[2][idxs])
  return (1-Θ)*y₀[idxs]+Θ*y₁[idxs]+Θ*(Θ-1)*((1-2Θ)*(y₁[idxs]-y₀[idxs])+(Θ-1)*dt*k[1][idxs] + Θ*dt*k[2][idxs])
end

@muladd function hermite_interpolant!(out,Θ,dt,y₀,y₁,k,idxs::Nothing,T::Type{Val{0}}) # Default interpolant is Hermite
  @inbounds @.. out = (1-Θ)*y₀+Θ*y₁+Θ*(Θ-1)*((1-2Θ)*(y₁-y₀)+(Θ-1)*dt*k[1] + Θ*dt*k[2])
  #@inbounds for i in eachindex(out)
  #  out[i] = (1-Θ)*y₀[i]+Θ*y₁[i]+Θ*(Θ-1)*((1-2Θ)*(y₁[i]-y₀[i])+(Θ-1)*dt*k[1][i] + Θ*dt*k[2][i])
  #end
  #out
end

@muladd function hermite_interpolant!(out,Θ,dt,y₀,y₁,k,idxs,T::Type{Val{0}}) # Default interpolant is Hermite
  @views @.. out = (1-Θ)*y₀[idxs]+Θ*y₁[idxs]+Θ*(Θ-1)*((1-2Θ)*(y₁[idxs]-y₀[idxs])+(Θ-1)*dt*k[1][idxs] + Θ*dt*k[2][idxs])
  #@inbounds for (j,i) in enumerate(idxs)
  #  out[j] = (1-Θ)*y₀[i]+Θ*y₁[i]+Θ*(Θ-1)*((1-2Θ)*(y₁[i]-y₀[i])+(Θ-1)*dt*k[1][i] + Θ*dt*k[2][i])
  #end
  #out
end

"""
Herimte Interpolation, chosen if no other dispatch for ode_interpolant
"""
@muladd function hermite_interpolant(Θ,dt,y₀,y₁,k,::Type{Val{false}},idxs::Nothing,T::Type{Val{1}}) # Default interpolant is Hermite
  #@.. k[1] + Θ*(-4*dt*k[1] - 2*dt*k[2] - 6*y₀ + Θ*(3*dt*k[1] + 3*dt*k[2] + 6*y₀ - 6*y₁) + 6*y₁)/dt
  @inbounds k[1] + Θ*(-4*dt*k[1] - 2*dt*k[2] - 6*y₀ + Θ*(3*dt*k[1] + 3*dt*k[2] + 6*y₀ - 6*y₁) + 6*y₁)/dt
end

@muladd function hermite_interpolant(Θ,dt,y₀,y₁,k,::Type{Val{true}},idxs::Nothing,T::Type{Val{1}}) # Default interpolant is Hermite
  @inbounds @.. k[1] + Θ*(-4*dt*k[1] - 2*dt*k[2] - 6*y₀ + Θ*(3*dt*k[1] + 3*dt*k[2] + 6*y₀ - 6*y₁) + 6*y₁)/dt
end

@muladd function hermite_interpolant(Θ,dt,y₀,y₁,k,cache,idxs,T::Type{Val{1}}) # Default interpolant is Hermite
  # return @.. k[1][idxs] + Θ*(-4*dt*k[1][idxs] - 2*dt*k[2][idxs] - 6*y₀[idxs] + Θ*(3*dt*k[1][idxs] + 3*dt*k[2][idxs] + 6*y₀[idxs] - 6*y₁[idxs]) + 6*y₁[idxs])/dt
  return k[1][idxs] + Θ*(-4*dt*k[1][idxs] - 2*dt*k[2][idxs] - 6*y₀[idxs] + Θ*(3*dt*k[1][idxs] + 3*dt*k[2][idxs] + 6*y₀[idxs] - 6*y₁[idxs]) + 6*y₁[idxs])/dt
end

@muladd function hermite_interpolant!(out,Θ,dt,y₀,y₁,k,idxs::Nothing,T::Type{Val{1}}) # Default interpolant is Hermite
  @inbounds @.. out = k[1] + Θ*(-4*dt*k[1] - 2*dt*k[2] - 6*y₀ + Θ*(3*dt*k[1] + 3*dt*k[2] + 6*y₀ - 6*y₁) + 6*y₁)/dt
  #@inbounds for i in eachindex(out)
  #  out[i] = k[1][i] + Θ*(-4*dt*k[1][i] - 2*dt*k[2][i] - 6*y₀[i] + Θ*(3*dt*k[1][i] + 3*dt*k[2][i] + 6*y₀[i] - 6*y₁[i]) + 6*y₁[i])/dt
  #end
  #out
end

@muladd function hermite_interpolant!(out,Θ,dt,y₀,y₁,k,idxs,T::Type{Val{1}}) # Default interpolant is Hermite
  @views @.. out = k[1][idxs] + Θ*(-4*dt*k[1][idxs] - 2*dt*k[2][idxs] - 6*y₀[idxs] + Θ*(3*dt*k[1][idxs] + 3*dt*k[2][idxs] + 6*y₀[idxs] - 6*y₁[idxs]) + 6*y₁[idxs])/dt
  #@inbounds for (j,i) in enumerate(idxs)
  #  out[j] = k[1][i] + Θ*(-4*dt*k[1][i] - 2*dt*k[2][i] - 6*y₀[i] + Θ*(3*dt*k[1][i] + 3*dt*k[2][i] + 6*y₀[i] - 6*y₁[i]) + 6*y₁[i])/dt
  #end
  #out
end

"""
Herimte Interpolation, chosen if no other dispatch for ode_interpolant
"""
@muladd function hermite_interpolant(Θ,dt,y₀,y₁,k,::Type{Val{false}},idxs::Nothing,T::Type{Val{2}}) # Default interpolant is Hermite
  #@.. (-4*dt*k[1] - 2*dt*k[2] - 6*y₀ + Θ*(6*dt*k[1] + 6*dt*k[2] + 12*y₀ - 12*y₁) + 6*y₁)/(dt*dt)
  @inbounds (-4*dt*k[1] - 2*dt*k[2] - 6*y₀ + Θ*(6*dt*k[1] + 6*dt*k[2] + 12*y₀ - 12*y₁) + 6*y₁)/(dt*dt)
end

@muladd function hermite_interpolant(Θ,dt,y₀,y₁,k,::Type{Val{true}},idxs::Nothing,T::Type{Val{2}}) # Default interpolant is Hermite
  #@.. (-4*dt*k[1] - 2*dt*k[2] - 6*y₀ + Θ*(6*dt*k[1] + 6*dt*k[2] + 12*y₀ - 12*y₁) + 6*y₁)/(dt*dt)
  @inbounds @.. (-4*dt*k[1] - 2*dt*k[2] - 6*y₀ + Θ*(6*dt*k[1] + 6*dt*k[2] + 12*y₀ - 12*y₁) + 6*y₁)/(dt*dt)
end

@muladd function hermite_interpolant(Θ,dt,y₀,y₁,k,cache,idxs,T::Type{Val{2}}) # Default interpolant is Hermite
  #out = similar(y₀,axes(idxs))
  #@views @.. out = (-4*dt*k[1][idxs] - 2*dt*k[2][idxs] - 6*y₀[idxs] + Θ*(6*dt*k[1][idxs] + 6*dt*k[2][idxs] + 12*y₀[idxs] - 12*y₁[idxs]) + 6*y₁[idxs])/(dt*dt)
  @views out = (-4*dt*k[1][idxs] - 2*dt*k[2][idxs] - 6*y₀[idxs] + Θ*(6*dt*k[1][idxs] + 6*dt*k[2][idxs] + 12*y₀[idxs] - 12*y₁[idxs]) + 6*y₁[idxs])/(dt*dt)
  out
end

@muladd function hermite_interpolant!(out,Θ,dt,y₀,y₁,k,idxs::Nothing,T::Type{Val{2}}) # Default interpolant is Hermite
  @inbounds @.. out = (-4*dt*k[1] - 2*dt*k[2] - 6*y₀ + Θ*(6*dt*k[1] + 6*dt*k[2] + 12*y₀ - 12*y₁) + 6*y₁)/(dt*dt)
  #@inbounds for i in eachindex(out)
  #  out[i] = (-4*dt*k[1][i] - 2*dt*k[2][i] - 6*y₀[i] + Θ*(6*dt*k[1][i] + 6*dt*k[2][i] + 12*y₀[i] - 12*y₁[i]) + 6*y₁[i])/(dt*dt)
  #end
  #out
end

@muladd function hermite_interpolant!(out,Θ,dt,y₀,y₁,k,idxs,T::Type{Val{2}}) # Default interpolant is Hermite
  @views @.. out = (-4*dt*k[1][idxs] - 2*dt*k[2][idxs] - 6*y₀[idxs] + Θ*(6*dt*k[1][idxs] + 6*dt*k[2][idxs] + 12*y₀[idxs] - 12*y₁[idxs]) + 6*y₁[idxs])/(dt*dt)
  #@inbounds for (j,i) in enumerate(idxs)
  #  out[j] = (-4*dt*k[1][i] - 2*dt*k[2][i] - 6*y₀[i] + Θ*(6*dt*k[1][i] + 6*dt*k[2][i] + 12*y₀[i] - 12*y₁[i]) + 6*y₁[i])/(dt*dt)
  #end
  #out
end

"""
Herimte Interpolation, chosen if no other dispatch for ode_interpolant
"""
@muladd function hermite_interpolant(Θ,dt,y₀,y₁,k,::Type{Val{false}},idxs::Nothing,T::Type{Val{3}}) # Default interpolant is Hermite
  #@.. (6*dt*k[1] + 6*dt*k[2] + 12*y₀ - 12*y₁)/(dt*dt*dt)
  @inbounds (6*dt*k[1] + 6*dt*k[2] + 12*y₀ - 12*y₁)/(dt*dt*dt)
end

@muladd function hermite_interpolant(Θ,dt,y₀,y₁,k,::Type{Val{true}},idxs::Nothing,T::Type{Val{3}}) # Default interpolant is Hermite
  #@.. (6*dt*k[1] + 6*dt*k[2] + 12*y₀ - 12*y₁)/(dt*dt*dt)
  @inbounds @.. (6*dt*k[1] + 6*dt*k[2] + 12*y₀ - 12*y₁)/(dt*dt*dt)
end

@muladd function hermite_interpolant(Θ,dt,y₀,y₁,k,cache,idxs,T::Type{Val{3}}) # Default interpolant is Hermite
  #out = similar(y₀,axes(idxs))
  #@views @.. out = (6*dt*k[1][idxs] + 6*dt*k[2][idxs] + 12*y₀[idxs] - 12*y₁[idxs])/(dt*dt*dt)
  @views out = (6*dt*k[1][idxs] + 6*dt*k[2][idxs] + 12*y₀[idxs] - 12*y₁[idxs])/(dt*dt*dt)
  out
end

@muladd function hermite_interpolant!(out,Θ,dt,y₀,y₁,k,idxs::Nothing,T::Type{Val{3}}) # Default interpolant is Hermite
  @inbounds @.. out = (6*dt*k[1] + 6*dt*k[2] + 12*y₀ - 12*y₁)/(dt*dt*dt)
  #for i in eachindex(out)
  #  out[i] = (6*dt*k[1][i] + 6*dt*k[2][i] + 12*y₀[i] - 12*y₁[i])/(dt*dt*dt)
  #end
  #out
end

@muladd function hermite_interpolant!(out,Θ,dt,y₀,y₁,k,idxs,T::Type{Val{3}}) # Default interpolant is Hermite
  @views @.. out = (6*dt*k[1][idxs] + 6*dt*k[2][idxs] + 12*y₀[idxs] - 12*y₁[idxs])/(dt*dt*dt)
  #for (j,i) in enumerate(idxs)
  #  out[j] = (6*dt*k[1][i] + 6*dt*k[2][i] + 12*y₀[i] - 12*y₁[i])/(dt*dt*dt)
  #end
  #out
end

######################## Linear Interpolants

@muladd @inline function linear_interpolant(Θ,dt,y₀,y₁,idxs::Nothing,T::Type{Val{0}})
  Θm1 = (1-Θ)
  @.. Θm1*y₀ + Θ*y₁
end

@muladd @inline function linear_interpolant(Θ,dt,y₀,y₁,idxs,T::Type{Val{0}})
  Θm1 = (1-Θ)
  @.. Θm1*y₀[idxs] + Θ*y₁[idxs]
end

@muladd @inline function linear_interpolant!(out,Θ,dt,y₀,y₁,idxs::Nothing,T::Type{Val{0}})
  Θm1 = (1-Θ)
  @.. out = Θm1*y₀ + Θ*y₁
  out
end

@muladd @inline function linear_interpolant!(out,Θ,dt,y₀,y₁,idxs,T::Type{Val{0}})
  Θm1 = (1-Θ)
  @views @.. out = Θm1*y₀[idxs] + Θ*y₁[idxs]
  out
end

"""
Linear Interpolation
"""
@inline function linear_interpolant(Θ,dt,y₀,y₁,idxs::Nothing,T::Type{Val{1}})
  (y₁ - y₀)/dt
end

@inline function linear_interpolant(Θ,dt,y₀,y₁,idxs,T::Type{Val{1}})
  @.. (y₁[idxs] - y₀[idxs])/dt
end

@inline function linear_interpolant!(out,Θ,dt,y₀,y₁,idxs::Nothing,T::Type{Val{1}})
  @.. out = (y₁ - y₀)/dt
  out
end

@inline function linear_interpolant!(out,Θ,dt,y₀,y₁,idxs,T::Type{Val{1}})
  @views @.. out = (y₁[idxs] - y₀[idxs])/dt
  out
end


text:
module ConstInferred
    export @constinferred, @constinferred_broken

    const ConstantValue = Union{Number,Char,QuoteNode}

    using InteractiveUtils: gen_call_with_extracted_types
    using Test
    using Test: Returned, Threw

    @static if isdefined(Base, :typesplit)
        const typesplit = Base.typesplit
    else
        Base.@pure function typesplit(@nospecialize(a), @nospecialize(b))
            if a <: b
                return Base.Bottom
            end
            if isa(a, Union)
                return Union{typesplit(a.a, b),
                             typesplit(a.b, b)}
            end
            return a
        end
    end

    _enabled = Ref(true)
    enable() = (_enabled[] = true; return nothing)
    disable() = (_enabled[] = false; return nothing)

    _materialize_broadcasted(f, args...) =
        Broadcast.materialize(Broadcast.broadcasted(f, args...))

    """
        @constinferred [AllowedType] f(x)
    Tests that the call expression `f(x)` returns a value of the same type inferred by the
    compiler. It is useful to check for type stability. Similar to `Test.@inferred`, but differs in two important ways.

    Firstly, `@constinferred` tries to test type stability under constant propagation by
    testing type stability of a new function, where all arguments or keyword arguments of
    the original function `f` that have a constant value (in the call expression) of type
    `Union{Number,Char,Symbol}` are hard coded.

    If you want to test for constant propagation in a variable which is not hard-coded in the call expression, you can interpolate it into the expression.

    Secondly, @constinferred returns the value if type inference succeeds, like `@inferred`,
    but used the `Test.@test` mechanism and shows up as an actual test error when type
    inference fails.
    ```
    """
    macro constinferred(ex)
        _constinferred(ex, __module__, __source__, Test.do_test)
    end
    macro constinferred(allow, ex)
        _constinferred(ex, __module__, __source__, Test.do_test, allow)
    end

    macro constinferred_broken(ex)
        _constinferred(ex, __module__, __source__, Test.do_broken_test)
    end

    macro constinferred_broken(allow, ex)
        _constinferred(ex, __module__, __source__, Test.do_broken_test, allow)
    end

    function _constinferred(ex, mod, src, test_f, allow = :(Union{}))
        if Meta.isexpr(ex, :ref)
            ex = Expr(:call, :getindex, ex.args...)
        end
        Meta.isexpr(ex, :call)|| error("@constinferred requires a call expression")
        farg = ex.args[1]
        if isa(farg, Symbol) && first(string(farg)) == '.'
            farg = Symbol(string(farg)[2:end])
            ex = Expr(:call, GlobalRef(Test, :_materialize_broadcasted),
                farg, ex.args[2:end]...)
        end
        pre = quote
                $(esc(allow)) isa Type || throw(ArgumentError("@constinferred requires a type as second argument"))
            end
        if length(ex.args) > 1 && Meta.isexpr(ex.args[2], :parameters)
            kwargs = ex.args[2].args
            args = ex.args[3:end]
        elseif length(ex.args) > 1 && Meta.isexpr(ex.args[2], :kw)
            kwargs = ex.args[2:end]
            args = Any[]
        else
            kwargs = Any[]
            args = ex.args[2:end]
        end
        newf = gensym()
        rightargs = Any[]
        rightkwargs = Any[]
        leftargs = Any[]
        callargs = Any[]
        quoteargs = Any[]
        for x in args
            if x isa ConstantValue
                push!(rightargs, x)
            elseif Meta.isexpr(x, :$)
                s = gensym()
                push!(rightargs, s)
                push!(quoteargs, Expr(:(=), s, x))
            else
                s = gensym()
                push!(leftargs, s)
                if Meta.isexpr(x, :...)
                    push!(rightargs, Expr(:..., s))
                    push!(callargs, Expr(:tuple, esc(x)))
                else
                    push!(rightargs, s)
                    push!(callargs, esc(x))
                end
            end
        end
        for x in kwargs
            xkey = x.args[1]
            xval = x.args[2]
            if xval isa ConstantValue
                push!(rightkwargs, x)
            elseif Meta.isexpr(xval, :$)
                s = gensym()
                push!(rightkwargs, Expr(:kw, xkey, s))
                push!(quoteargs, Expr(:(=), s, xval))
            else
                s = gensym()
                push!(rightkwargs, Expr(:kw, xkey, s))
                push!(leftargs, s)
                push!(callargs, esc(xval))
            end
        end
        f = Expr(:$, ex.args[1])
        fundefhead = Expr(:call, newf, leftargs...)
        fundefbody = Expr(:block, quoteargs..., isempty(kwargs) ?
                            Expr(:call, f, rightargs...) :
                            Expr(:call, f, Expr(:parameters, rightkwargs...), rightargs...))
        fundefex = esc(Expr(:quote, Expr(:(=), fundefhead, fundefbody)))

        inftypes = gensym()
        rettype = gensym()
        result = esc(gensym())
        newfcall = Expr(:., mod, QuoteNode(newf))
        main = quote
            callargs = ($(callargs...),)
            $result = $newfcall(callargs...)
            $(esc(inftypes)) = Base.return_types($newfcall, Base.typesof(callargs...))
            $(esc(rettype)) = $result isa Type ? Type{$result} : typeof($result)
        end
        orig_ex = Expr(:inert, Expr(:macrocall, Symbol("@constinferred"), nothing, ex))
        post = quote
            if $(_enabled[])
                if length($(esc(inftypes))) > 1
                    testresult = Threw(ArgumentError("more than one inferred type"), Base.catch_stack(), $(QuoteNode(src)))
                else
                    v = $(esc(rettype)) <: $(esc(allow)) || $(esc(rettype)) == typesplit($(esc(inftypes))[1], $(esc(allow)))
                    testresult = Returned(v, Expr(:call, :!=, $(esc(rettype)), $(esc(inftypes))[1]), $(QuoteNode(src)))
                end

                $test_f(testresult, $orig_ex)
            end
            $result
        end
        finalex = Base.remove_linenums!(quote
            $pre
            Core.eval($mod, $fundefex)
            let
                $main
                $post
            end
        end)
        finalex
    end
end


text:
using AmplNLWriter, JuMP, Ipopt
using Base.Test

include("nl_convert.jl")
include("nl_linearity.jl")
include("nl_write.jl")
include("sol_file_parser.jl")

# needed for the scoping of `solver` in the examples
solver = JuMP.UnsetSolver()
solvers = Any[]
push!(solvers, AmplNLSolver(Ipopt.amplexe, ["print_level=0"]))

examples_path = joinpath(dirname(dirname(@__FILE__)), "examples")

for solver in solvers
    solvername = getsolvername(solver)
    @testset "[examples] test solver $solvername" begin
        for example in [
                "jump_nltrig.jl", "jump_nlexpr.jl", "jump_pruning.jl",
                "jump_minlp.jl", "jump_nonlinearbinary.jl", "jump_no_obj.jl",
                "jump_const_obj.jl", "jump_maxmin.jl"
            ]
            include(joinpath(examples_path, example))
        end
    end
end

include(Pkg.dir("JuMP","test","solvers.jl"))
include(Pkg.dir("JuMP","test","nonlinear.jl"))


text:
using Pluto
Pluto.run("0.0.0.0", 1234;
	configuration=Pluto.ServerConfiguration(launch_browser=false),
	security=Pluto.ServerSecurity(false),
	)


text:
module HPM

using ..LIKWID: LibLikwid, access_initialized

"""
Sets the mode how the MSR and PCI registers should be accessed. Available options:
  * `0` or `LibLikwid.ACCESSMODE_DIRECT`: direct access (propably root priviledges required)
  * `1` or `LibLikwid.ACCESSMODE_DAEMON`: accesses through the access daemon
Must be called **before** [`HPM.init`](@ref).
"""
function mode(mode::Union{Integer, LibLikwid.AccessMode})
    if mode == LibLikwid.ACCESSMODE_DIRECT || mode == LibLikwid.ACCESSMODE_DAEMON
        LibLikwid.HPMmode(mode)
        return true
    end
    return false
end

"""
Initialize the access module internals to either the MSR/PCI files or the access daemon
"""
function init()
    err = LibLikwid.HPMinit()
    if err == 0
        access_initialized[] = true
        return true
    end
    return false
end

"""
Add the given CPU to the access module. This opens the commnunication to either the MSR/PCI files or the access daemon.
"""
function add_thread(cpuid)
    if !access_initialized[]
        return -1
    end
    ret = LibLikwid.HPMaddThread(cpuid)
    return Int(ret)
end

"""
Close the connections to the MSR/PCI files or the access daemon.
"""
function finalize()
    LibLikwid.HPMfinalize()
    access_initialized[] = false
    return nothing
end

end # module


text:
"""
`Block` is used for get index of a block matrix.  For example,
```julia
A[Block(1),Block(2)]
```
retrieves the 1 x 2 block
"""
struct Block
    K::Int
end

convert(::Type{T}, B::Block) where T<:Integer = convert(T, B.K)::T

for OP in (:(one), :(zero), :(+), :(-))
    @eval $OP(B::Block) = Block($OP(B.K))
end

for OP in (:(==), :(!=), :(isless))
    @eval $OP(A::Block,B::Block) = $OP(A.K,B.K)
end

for OP in (:(div), :(rem))
    @eval begin
        $OP(A::Block,B::Block) = Block($OP(A.K,B.K))
        $OP(A::Integer,B::Block) = Block($OP(A,B.K))
        $OP(A::Block,B::Integer) = Block($OP(A.K,B))
    end
end


text:
# This file was generated, do not modify it.

using MLJ
X, y = @load_iris
@load DecisionTreeClassifier

dtc = DecisionTreeClassifier()
r   = range(dtc, :max_depth, lower=1, upper=5)

tm = TunedModel(model=dtc, ranges=[r, ], measure=cross_entropy)

m = machine(tm, X, y)
fit!(m)

fitted_params(m).best_model.max_depth

X = (x1=rand(100), x2=rand(100), x3=rand(100))
y = 2X.x1 - X.x2 + 0.05 * randn(100);

dtr = @load DecisionTreeRegressor
forest = EnsembleModel(atom=dtr)

r1 = range(forest, :(atom.n_subfeatures), lower=1, upper=3)
r2 = range(forest, :bagging_fraction, lower=0.4, upper=1.0)
tm = TunedModel(model=forest, tuning=Grid(resolution=12),
                resampling=CV(nfolds=6), ranges=[r1, r2],
                measure=rms)
m = machine(tm, X, y)
fit!(m);

r = report(m)
r.best_measurement

max_stars_repo_name: tlienart/MLJTutorials

max_issues_repo_name: tlienart/MLJTutorials

max_forks_repo_name: tlienart/MLJTutorials


text:
"""
    ObjectiveFunction{VT}

Represents an objective function optimized by algorithms.
Subtypes of `ObjectiveFunction` must implement at least
* `compute_value(::ObjectiveFunction, x)` for primal value evaluation
* `compute_gradient(::ObjectiveFunction, x)` for gradient evaluation.
and optionally `compute_value_gradient(::ObjectiveFunction, x)` returning the (primal, gradient) pair.
"""
abstract type ObjectiveFunction end

function compute_value end

function compute_gradient end

"""
    compute_value_gradient(f::ObjectiveFunction, x; kwargs)

Computes in one call the pair `(function_value, function_grad)` evaluated at `x`.
By default, calls `compute_value` and `compute_gradient` with keyword `kwargs`
passed to both.
"""
compute_value_gradient(f::ObjectiveFunction, x; kwargs...) =
    (compute_value(f, x; kwargs...), compute_gradient(f, x; kwargs...))

struct SimpleFunctionObjective{F,G} <: ObjectiveFunction
    f::F
    grad::G
end

compute_value(f::SimpleFunctionObjective, x) = f.f(x)
compute_gradient(f::SimpleFunctionObjective, x) = f.grad(x)

"""
    StochasticObjective{F, G, XT}(f::F, grad::G, xs::XT)

Represents an objective function evaluated with stochastic gradient.
`f(θ, x)` evaluates the loss for data point `x` and parameter `θ`.
`grad(θ, x)` evaluates the loss gradient with respect to data point `x` at parameter `θ`.
`xs` must be an indexable iterable (`Vector{Vector{Float64}}` for instance).
Functions using a `StochasticObjective` have optional keyword arguments `rng`, `batch_size`
and `full_evaluation` controlling whether the function should be evaluated over all data points.
"""
struct StochasticObjective{F,G,XT} <: ObjectiveFunction
    f::F
    grad::G
    xs::XT
end

function compute_value(
    f::StochasticObjective,
    θ;
    batch_size::Integer=length(f.xs),
    rng=Random.GLOBAL_RNG,
    full_evaluation=false,
)
    (batch_size, rand_indices) = if full_evaluation
        (length(f.xs), eachindex(f.xs))
    else
        (batch_size, rand(rng, eachindex(f.xs), batch_size))
    end
    return sum(f.f(θ, f.xs[idx]) for idx in rand_indices) / batch_size
end

function compute_gradient(
    f::StochasticObjective,
    θ;
    batch_size::Integer=length(f.xs) ÷ 10 + 1,
    rng=Random.GLOBAL_RNG,
    full_evaluation=false,
)
    (batch_size, rand_indices) = if full_evaluation
        (length(f.xs), eachindex(f.xs))
    else
        (batch_size, rand(rng, eachindex(f.xs), batch_size))
    end
    return sum(f.grad(θ, f.xs[idx]) for idx in rand_indices) / batch_size
end

function compute_value_gradient(
    f::StochasticObjective,
    θ;
    batch_size::Integer=length(f.xs) ÷ 10 + 1,
    rng=Random.GLOBAL_RNG,
    full_evaluation=false,
)
    (batch_size, rand_indices) = if full_evaluation
        (length(f.xs), eachindex(f.xs))
    else
        (batch_size, rand(rng, eachindex(f.xs), batch_size))
    end
    # map operation, for each index, computes value and gradient
    function map_op(idx)
        @inbounds x = f.xs[idx]
        return (f.f(θ, x) / batch_size, f.grad(θ, x) / batch_size)
    end
    # reduce: take partial value and gradient, adds value and gradient wrt new point
    function reduce_op(left_tup, right_tup)
        (f_val, g_val) = left_tup
        (f_new, g_new) = right_tup
        return (f_val + f_new, g_val + g_new)
    end

    return mapfoldr(map_op, reduce_op, rand_indices)
end


text:
### A Pluto.jl notebook ###
# v0.12.17

using Markdown
using InteractiveUtils

# This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error).
macro bind(def, element)
    quote
        local el = $(esc(element))
        global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : missing
        el
    end
end

# ╔═╡ fff0ff7a-4a94-11eb-1821-fbc6796d2e2d
begin
	using FSM
	using PlutoUI
	using Plots
end

# ╔═╡ 7e7db788-4a96-11eb-3a41-2115cc7d3b58
md"""
# Albedo models
"""

# ╔═╡ 7595e612-4ab0-11eb-3a50-eb114ed9507d
md"""
## Model setup
"""

# ╔═╡ 24c749a8-4a95-11eb-1fc8-f34a2d08bdbe
ebm = EBM{Float64}()

# ╔═╡ d66ab284-4a96-11eb-261d-b73fc2de22c2
cn = Constants{Float64}()

# ╔═╡ fa0dd0ec-4ab2-11eb-34ce-d16073a73f39
md"""
## Albedo function
"""

# ╔═╡ 4e6b1406-4a95-11eb-0827-432daf4b223e
function compute_albedo(ebm::EBM, cn::Constants, Sf)

	if ebm.am == 0  # Diagnosed snow albedo
		ebm.albs = ebm.asmn + (ebm.asmx - ebm.asmn) * (ebm.Tsurf - cn.Tm) / ebm.Talb
	elseif ebm.am == 1  # Prognostic snow albedo
		tau = 3600 * ebm.tcld
		if (ebm.Tsurf >= cn.Tm)
			tau = 3600 * ebm.tmlt
		end
		rt = 1 / tau + Sf / ebm.Salb
		alim = (ebm.asmn / tau + Sf * ebm.asmx / ebm.Salb) / rt
		ebm.albs = alim + (ebm.albs - alim) * exp(-rt * ebm.dt)
	end

	if (ebm.albs < min(ebm.asmx, ebm.asmn))
		ebm.albs = min(ebm.asmx, ebm.asmn)
	end

	if (ebm.albs > max(ebm.asmx, ebm.asmn))
		ebm.albs = max(ebm.asmx, ebm.asmn)
	end

end

# ╔═╡ 0b59e04c-4aaf-11eb-33ff-f9e05476a28f
begin
	Tsurf = cn.Tm-10:0.01:cn.Tm+10
	Sf = 0

	albedo = []
	for T in Tsurf
		eb = EBM{Float64}(Tsurf = T, am = 0)
		compute_albedo(eb, cn, Sf)
		push!(albedo, eb.albs)
	end
end

# ╔═╡ f866f932-4ab0-11eb-0c58-0dd947ff1341
md"""
## Diagnostic albedo model
"""

# ╔═╡ be13ab46-4aaf-11eb-1f96-752bcbef74f6
plot(Tsurf, albedo, xlabel = "Temperature (K)", ylabel = "Albedo (-)", title = "Albedo model = 0")

# ╔═╡ 1427ffea-4ab1-11eb-1862-254bcabf97ac
md"""
## Prognostic albedo model
"""

# ╔═╡ 798de214-4ab3-11eb-0d8e-f7fe85ae089e
md"""
Add snow at hour = 100
"""

# ╔═╡ 5512e630-4ab3-11eb-2347-8f4afbd5f228
@bind Sf_one Slider(0:200, show_value=true)

# ╔═╡ 20190330-4ab1-11eb-33b4-e9ccd85a469a
begin
	eb_cold = EBM{Float64}(Tsurf = 271.15, am = 1)
	eb_warm = EBM{Float64}(Tsurf = 273.15, am = 1)
	albs_cold = [eb_cold.albs]
	albs_warm = [eb_warm.albs]
	t_prog = [0]
	for t in 1:24*10

		if t == 100
			Sf = Sf_one/3600
		else
			Sf = 0
		end

		compute_albedo(eb_cold, cn, Sf)
		compute_albedo(eb_warm, cn, Sf)
		push!(albs_cold, eb_cold.albs)
		push!(albs_warm, eb_warm.albs)
		push!(t_prog, t)
	end
end

# ╔═╡ 2016de88-4ab2-11eb-3ca1-a75c6752bc71
begin
	plot(t_prog, albs_cold)
	plot!(t_prog, albs_warm)
	plot!(xlabel = "Temperature (K)", ylabel = "Albedo (-)", title = "Albedo model = 1")
	plot!(ylim = (0, 1))
end

# ╔═╡ Cell order:
# ╟─7e7db788-4a96-11eb-3a41-2115cc7d3b58
# ╟─7595e612-4ab0-11eb-3a50-eb114ed9507d
# ╟─24c749a8-4a95-11eb-1fc8-f34a2d08bdbe
# ╟─d66ab284-4a96-11eb-261d-b73fc2de22c2
# ╟─fff0ff7a-4a94-11eb-1821-fbc6796d2e2d
# ╟─fa0dd0ec-4ab2-11eb-34ce-d16073a73f39
# ╠═4e6b1406-4a95-11eb-0827-432daf4b223e
# ╟─0b59e04c-4aaf-11eb-33ff-f9e05476a28f
# ╟─f866f932-4ab0-11eb-0c58-0dd947ff1341
# ╟─be13ab46-4aaf-11eb-1f96-752bcbef74f6
# ╟─1427ffea-4ab1-11eb-1862-254bcabf97ac
# ╟─20190330-4ab1-11eb-33b4-e9ccd85a469a
# ╟─798de214-4ab3-11eb-0d8e-f7fe85ae089e
# ╟─5512e630-4ab3-11eb-2347-8f4afbd5f228
# ╟─2016de88-4ab2-11eb-3ca1-a75c6752bc71


text:
# Apply default faceting look to a grid of `AxisEntries`

function facet_wrap!(fig, aes::AbstractMatrix{AxisEntries})
    scale = get(aes[1].scales, :layout, nothing)
    isnothing(scale) && return
    linkaxes!(aes...)
    for ae in aes
        ax = ae.axis
        entries = ae.entries
        vs = Iterators.filter(
            !isnothing,
            (get(entry.primary, :layout, nothing) for entry in entries)
        )
        it = iterate(vs)
        if isnothing(it)
            delete!(ax)
        else
            v, _ = it
            ax.title[] = string(v)
        end
    end
    return
end

# F. Greimel implementation from AlgebraOfGraphics

# TODO: add configuration options here (esp. to determine when / how to link)
function facet_grid!(fig, aes::AbstractMatrix{AxisEntries})
    M, N = size(aes)
    row_scale, col_scale = map(sym -> get(aes[1].scales, sym, nothing), (:row, :col))
    all(isnothing, (row_scale, col_scale)) && return
    hideinnerdecorations!(aes)
    linkaxes!(aes...)

    nonempty_aes = filter(ae -> !isempty(ae.entries), aes)

    ax = first(nonempty_aes).axis
    titlegap = ax.titlegap
    titlecolor = ax.titlecolor
    titlefont = ax.titlefont
    titlesize = ax.titlesize
    facetlabelattributes = (
        color=titlecolor,
        font=titlefont,
        textsize=titlesize,
    )

    consistent_xlabels = all(ae -> ae.axis.xlabel[] == ax.xlabel[], nonempty_aes)
    consistent_ylabels = all(ae -> ae.axis.ylabel[] == ax.ylabel[], nonempty_aes)

    if !isnothing(row_scale) && consistent_ylabels
        for ae in aes
            ae.axis.ylabelvisible[] = false
        end
        row_dict = Dict(zip(plotvalues(row_scale), datavalues(row_scale)))
        facetlabelpadding = lift(titlegap) do gap
            return (gap, 0f0, 0f0, 0f0)
        end
        for m in 1:M
            Label(fig[m, N, Right()], string(row_dict[m]);
                rotation=-π/2, padding=facetlabelpadding,
                facetlabelattributes...)
        end
        protrusion = lift(
            (xs...) -> maximum(x -> x.left, xs),
            (MakieLayout.protrusionsobservable(ae.axis) for ae in aes[:, 1])...
        )
        # TODO: here and below, set in such a way that one can change padding after the fact?
        ylabelpadding = lift(protrusion, ax.ylabelpadding) do val, p
            return (0f0, val + p, 0f0, 0f0)
        end
        ylabelcolor = ax.ylabelcolor
        ylabelfont = ax.ylabelfont
        ylabelsize = ax.ylabelsize
        ylabelattributes = (
            color=ylabelcolor,
            font=ylabelfont,
            textsize=ylabelsize,
        )
        Label(fig[:, 1, Left()], ax.ylabel;
            rotation=π/2, padding=ylabelpadding, ylabelattributes...)
    end
    if !isnothing(col_scale) && consistent_xlabels
        for ae in aes
            ae.axis.xlabelvisible[] = false
        end
        col_dict = Dict(zip(plotvalues(col_scale), datavalues(col_scale)))
        labelpadding = lift(titlegap) do gap
            return (0f0, 0f0, gap, 0f0)
        end
        for n in 1:N
            Label(fig[1, n, Top()], string(col_dict[n]);
                padding=labelpadding, facetlabelattributes...)
        end
        protrusion = lift(
            (xs...) -> maximum(x -> x.bottom, xs),
            (MakieLayout.protrusionsobservable(ae.axis) for ae in aes[M, :])...
        )
        xlabelpadding = lift(protrusion, ax.xlabelpadding) do val, p
            return (0f0, 0f0, 0f0, val + p)
        end
        xlabelcolor = ax.xlabelcolor
        xlabelfont = ax.xlabelfont
        xlabelsize = ax.xlabelsize
        xlabelattributes = (
            color=xlabelcolor,
            font=xlabelfont,
            textsize=xlabelsize,
        )
        Label(fig[M, :, Bottom()], ax.xlabel;
            padding=xlabelpadding, xlabelattributes...)
    end
    return
end

function facet!(fig, aes::AbstractMatrix{AxisEntries})
    facet_wrap!(fig, aes)
    facet_grid!(fig, aes)
    return
end

function facet!(fg::FigureGrid)
    facet!(fg.figure, fg.grid)
    return fg
end

## Layout helpers

isaxis2d(::Axis) = true
isaxis2d(::Any) = false
isaxis2d(ae::AxisEntries) = isaxis2d(ae.axis)

for sym in [:hidexdecorations!, :hideydecorations!, :hidedecorations!]
    @eval function $sym(ae::AxisEntries; kwargs...)
        axis = ae.axis
        isaxis2d(axis) && $sym(axis; kwargs...)
    end
end

for sym in [:linkxaxes!, :linkyaxes!, :linkaxes!]
    @eval function $sym(ae::AxisEntries, aes::AxisEntries...)
        axs = filter(isaxis2d, map(ae->ae.axis, (ae, aes...)))
        isempty(axs) || $sym(axs...)
    end
end

function hideinnerdecorations!(aes::Matrix{AxisEntries})
    options = (label=true, ticks=true, minorticks=true, grid=false, minorgrid=false)
    foreach(ae -> hidexdecorations!(ae; options...), aes[1:end-1, :])
    foreach(ae -> hideydecorations!(ae; options...), aes[:, 2:end])
end

function deleteemptyaxes!(aes::Matrix{AxisEntries})
    for ae in aes
        if isempty(ae.entries)
            delete!(ae.axis)
        end
    end
end

function resizetocontent!(fig::Figure)
    figsize = size(fig.scene)
    sz = map((Col(), Row()), figsize) do dir, currentsize
        inferredsize = determinedirsize(fig.layout, dir)
        return ceil(Int, something(inferredsize, currentsize))
    end
    sz == figsize || resize!(fig.scene, sz)
    return fig
end

function resizetocontent!(fg::FigureGrid)
    resizetocontent!(fg.figure)
    return fg
end


text:
#NNTest.jl

#NN.jl testbed
include("NN.jl")
using NN, PyPlot

x = linspace(-0.5,0.5,100)#(rand(1,10000)-0.5)
f = (x)->x
#f = (x)->exp(x)
#f = (x)->x.*x.*x
#f = (x)->sin(x)
#f = (x)->x.*x+1
#f = (x)->sin(exp(x.*x.*x-x+12)).^2
y = f(x)

alpha_final = 0.1
nn = buildNN([1;2;1],0.00001,minimum(y),maximum(y),neuron_type="relu",reg_type="L2",verbose=false,init_type="rand",output_neuron_type="none")
nn,losses, losses_true,alpha_final = train(x,y,nn,alpha=alpha_final,num_epochs=3,minibatch_size=1,update_type="rmsprop",half_rate=1000)
x_ = linspace(-0.5,0.5,100)
y_ = [predict(nn,[x_[i]]) for i = 1:length(x_)]
#plotting stuff
#plot(x_,f(x_),x_,y_)


text:
using MultipleViewGeometry, Test
using MultipleViewGeometry.ModuleTypes

ℳ = [Point2D(rand(1,2)) for i = 1:10000]

moments(FundamentalMatrix(), (ℳ,ℳ))

@benchmark moments(HomographyMatrix(), (ℳ,ℳ))

# TODO

FundamentalMatrix()
HomographyMatrix()


text:
#=
Modia main module.

* Developers: Hilding Elmqvist, Mogram AB, Toivo Henningsson, Lund and Martin Otter, DLR
* Copyright (c) 2016-2018: Hilding Elmqvist, Toivo Henningsson, Martin Otter
* License: MIT (expat)
=#

#import ModiaMath

"""
Modia - Dynamic Modeling and Simulation in Julia

To define a model:

```julia
  @model FirstOrder begin
     x = Variable(start=1)   # start means x(0)
     T = Parameter(0.5)      # Time constant
     u = 2.0                 # Same as Parameter(2.0)
  @equations begin
     T*der(x) + x = u        # der() means time derivative
     end
  end;
```

To simulate a model:
```julia
  result = simulate(FirstOrder, 2);
  @show result["x"][end];
  ModiaMath.plot(result, "x")
```

To run examples:
```julia
  include("$(Modia.ModiaDir)/examples/runexamples.jl")
```

To run tests:
```julia
  include("$(Modia.ModiaDir)/test/runtests.jl")
```
For more information, see (https://github.com/ModiaSim/Modia.jl/blob/master/README.md)
"""
module Modia

const Version = "0.2.1"
const Date = "2018-10-07"

#println(" \n\nWelcome to Modia - Dynamic MODeling and Simulation in julIA")
print(" \n\nWelcome to ")
print("Mod")
@static if VERSION < v"0.7.0-DEV.2005"
    print_with_color(:red, "ia", bold=true)
    print(" - ")
    print_with_color(:light_black, "Dynamic ")
    print("Mod")
    print_with_color(:light_black, "eling and Simulation with Jul")
    print_with_color(:red, "ia", bold=true)
else
    printstyled("ia", bold=true, color=:red)
    print(" - ")
    printstyled("Dynamic ", color=:light_black)
    print("Mod")
    printstyled("eling and Simulation with Jul", color=:light_black)
    printstyled("ia", bold=true, color=:red)
end

println()
println("Version $Version ($Date)")
println("Type \"?Modia\" for help.\n\n")

const ModiaDir = dirname(Base.@__DIR__)

export ModiaDir

export @model, simulateModel, simulate, checkSimulation, simulateMultiModeModel
export Variable, Float, Boolean, Integ, Str, Parameter, Var, Par, undefined
export Variability, constant, parameter, discrete, continuous
export Property, general, symmetric, orthogonal, rotationGroup3D
export skew, skewCoords
export allInstances
export @component, addComponent!

# Import packages that are used in examples and tests
# (in order that there are no requirements on the environment
#  in which the examples and tests are executed).
import DataStructures
import ModiaMath
import StaticArrays
import Unitful

@static if VERSION >= v"0.7.0-DEV.2005"
    import LinearAlgebra
    import Markdown
    import SparseArrays
    import Test
end

# Make string macro U_str available by "using Modia"
include("language/Unitful_U_str.jl")
using  .Unitful_U_str
export @U_str

# Include all sub-modules of Modia
include("language/ModiaLogging.jl")
include("language/Instantiation.jl")
include("language/Execution.jl")

include("symbolic/DAEquations/Synchronous.jl") # Before models/Electric, etc
include("symbolic/symbolic.jl")

using .Instantiation

include("models/ModiaBase.jl")
include("models/models.jl")  # Before symbolic because MultiBody is used in BasicStructuralTransform

end


text:
using Test

using Stonks: DataClientError
using Stonks.APIClients: APIClient, APIResource, YahooClient
using Stonks.APIClients:
  get_type_param,
  get_supported_types,
  get_resource,
  get_resource_from_clients,
  build_clients_from_env
using Stonks.Models:
  AbstractStonksRecord, AssetInfo, AssetPrice, ExchangeRate, IncomeStatement, BalanceSheet, CashflowStatement, Earnings

@testset "Client building utilities" begin
  test_client = YahooClient("abc")
  function unset_env_vars()
    for (k, _) in ENV
      (endswith(k, "_TOKEN") || endswith(k, "_APIKEY")) && delete!(ENV, k)
    end
  end
  set_env_vars() = foreach(k -> ENV[k] = "secret", ["ALPHAVANTAGE_TOKEN", "YAHOO_TOKEN"])
  unset_env_vars()

  @testset "Get type parameter of API resource" begin
    T = get_type_param(test_client.resources["price"])
    @test T === AssetPrice
  end

  @testset "Get supported types" begin
    types = get_supported_types(test_client)
    @test isempty(
      setdiff(types, [ExchangeRate, AssetPrice, AssetInfo, IncomeStatement, BalanceSheet, CashflowStatement, Earnings])
    )
  end

  @testset "Get resource from client of a certain type" begin
    resource = get_resource(test_client, AssetInfo)
    @test get_type_param(resource) === AssetInfo
  end

  @testset "Build clients from environmental vars, no env vars set" begin
    unset_env_vars()
    @test isa(build_clients_from_env(), DataClientError)
  end

  @testset "Build clients from environmental vars, 2 env vars set" begin
    set_env_vars()
    clients = build_clients_from_env()
    @test length(clients) == 2
    @test isa(clients, Vector{APIClient})
    unset_env_vars()
  end

  @testset "Get one resource from multiple clients" begin
    set_env_vars()
    clients = build_clients_from_env()
    resource = get_resource_from_clients(clients, AssetPrice)
    @test isa(resource, APIResource{AssetPrice})
    unset_env_vars()
  end

  @testset "Get resource with no configured clients" begin
    set_env_vars()
    clients = build_clients_from_env()
    struct MyCustomType <: AbstractStonksRecord
      value::String
    end
    resource = get_resource_from_clients(clients, MyCustomType)
    @test isa(resource, DataClientError)
  end
end


text:
function _equalarrays(a1::AbstractArray{T1}, a2::AbstractArray{T2}) where T1 where T2
    #=
    Compare two arrays for equality.
    =#

    size(a1) == size(a2) ? all(a1 .== a2) : false

end

# exact numeric type
ExactNumber = Union{Complex{Rational{I1}}, Rational{I2}, I3, Complex{I4}} where I1 <: Integer where
I2 <: Integer where I3 <: Integer where I4 <: Integer

function simplify(n::ExactNumber)
    #=
    Represent an ExactNumber in the simplest way possible.
    =#

    newn = n
    ir = imag(n) == 0

    if ir
        try
            newn = convert(Int64, n)
        catch err
            nothing
        end
    else
        try
            newn = convert(Complex{Int64}, n)
        catch err
            nothing
        end
    end

    newn

end

function _compresssum(v::AbstractVector{T}) where T

    # if elements can be added together, do so
    cv = Vector{T}()
    for e1 in v

        # try to add this element to one of the previous ones
        matched = false
        for (i, e2) in enumerate(cv)

            if addable(e1, e2)  # we found a place to add to
                cv[i] = e1 + e2
                matched = true
                break
            end

        end

        if !matched # we didn't find a spot to add to, so make a new entry
            push!(cv, e1)
        end

    end

    cv

end

function _cancelpairs!(v::AbstractVector{T}, u::AbstractVector{T}) where T
    #=
    Delete pairs of equal elements in the two vectors.
    =#

    lv = length(v)
    lu = length(u)
    delv = falses(lv)
    delu = falses(lu)

    for (i, x) in enumerate(v), (j, y) in enumerate(u)

        if x==y && !delv[i] && !delu[j] # delete if not already
            delv[i] = true
            delu[j] = true
        end

    end

    deleteat!(v, delv)
    deleteat!(u, delu)

    nothing

end

max_stars_repo_name: aormorningstar/SymbolicMajoranas

max_issues_repo_name: aormorningstar/SymbolicMajoranas

max_forks_repo_name: aormorningstar/SymbolicMajoranas


text:
function checkout_branch(dir::AbstractString,
    branch::AbstractString;
    git_command::AbstractString = "git")
    Base.run(Cmd(`$(git_command) checkout $(branch)`; dir=dir))
end

clone_repo(repo::GitHub.Repo) = clone_repo(repo_url(repo))

function clone_repo(url::AbstractString)
    parent_dir = mktempdir(; cleanup=true)
    repo_dir = joinpath(parent_dir, "REPO")
    my_retry(() -> _clone_repo_into_dir(url, repo_dir))
    @info("Clone was successful")
    return repo_dir
end

function _clone_repo_into_dir(url::AbstractString, repo_dir)
    @info("Attempting to clone...")
    rm(repo_dir; force = true, recursive = true)
    mkpath(repo_dir)
    LibGit2.clone(url, repo_dir)
    return repo_dir
end

"""
    load_files_from_url_and_tree_hash(f, destination::String, url::String, tree_hash::String) -> Bool

Attempts to clone a git repo from `url` into a temporary directory, runs `f(dir)` on that directory,
then extract the files and folders from a given `tree_hash`, placing them in `destination`.

Returns a boolean indicating if the cloning succeeded.
"""
function load_files_from_url_and_tree_hash(f, destination::String, url::String, tree_hash::String)
    pkg_clone_dir = mktempdir()
    clone_success = try
        _clone_repo_into_dir(url, pkg_clone_dir)
        true
    catch e
        @error "Cloning $url failed" e
        false
    end
    # if cloning failed, bail now
    !clone_success && return clone_success

    f(pkg_clone_dir)
    Tar.extract(Cmd(`git archive $tree_hash`; dir=pkg_clone_dir), destination)
    return clone_success
end

"""
    parse_registry_pkg_info(registry_path, pkg, version=nothing) -> @NamedTuple{uuid::String, repo::String, subdir::String, tree_hash::Union{Nothing, String}}

Searches the registry located at `registry_path` for a package with name `pkg`. Upon finding it,
it parses the associated `Package.toml` file and returns the UUID and repository URI, and `subdir`.

If `version` is supplied, then the associated `tree_hash` will be returned. Otherwise, `tree_hash` will be `nothing`.
"""
function parse_registry_pkg_info(registry_path, pkg, version=nothing)
    # We know the name of this package but not its uuid. Look it up in
    # the registry that includes the current PR.
    packages = TOML.parsefile(joinpath(registry_path, "Registry.toml"))["packages"]
    filter!(packages) do (key, value)
        value["name"] == pkg
    end
    # For Julia >= 1.4 this can be simplified with the `only` function.
    always_assert(length(packages) == 1)
    uuid = convert(String, first(keys(packages)))
    # Also need to find out the package repository.
    package = TOML.parsefile(joinpath(registry_path, packages[uuid]["path"], "Package.toml"))
    repo = convert(String, package["repo"])
    subdir = convert(String, get(package, "subdir", ""))
    if version === nothing
        tree_hash = nothing
    else
        versions = TOML.parsefile(joinpath(registry_path, packages[uuid]["path"], "Versions.toml"))
        tree_hash = convert(String, versions[string(version)]["git-tree-sha1"])
    end
    return (; uuid=uuid, repo=repo, subdir=subdir, tree_hash=tree_hash)
end

function _comment_disclaimer()
    result = string("\n\n",
                    "Note that the guidelines are only required for the pull request ",
                    "to be merged automatically. However, it is **strongly recommended** ",
                    "to follow them, since otherwise the pull request needs to be ",
                    "manually reviewed and merged by a human.",
                    "\n\n",
                    "After you have fixed the AutoMerge issues, simple retrigger Registrator, ",
                    "which will automatically update this pull request. ",
                    "You do not need to change the version number in your `Project.toml` file ",
                    "(unless of course the AutoMerge issue is that you skipped a version number, ",
                    "in which case you should change the version number).",
                    "",
                    "\n\n",
                    "If you do not want to fix the AutoMerge issues, please post a comment ",
                    "explaining why you would like this pull request to be manually merged.",
                    "",
                    )
    return result
end

function _comment_noblock()
    result = string("\n\n---\n",
                    "If you want to prevent this pull request from ",
                    "being auto-merged, simply leave a comment. ",
                    "If you want to post a comment without blocking ",
                    "auto-merging, you must include the text ",
                    "`[noblock]` in your comment.")
    return result
end

function comment_text_pass(::NewVersion,
                           suggest_onepointzero::Bool,
                           version::VersionNumber)
    result = string("Your `new version` pull request met all of the ",
                    "guidelines for auto-merging and is scheduled to ",
                    "be merged in the next round.",
                    _comment_noblock(),
                    _onepointzero_suggestion(suggest_onepointzero, version),
                    "\n<!-- [noblock] -->")
    return result
end

function comment_text_pass(::NewPackage,
                           suggest_onepointzero::Bool,
                           version::VersionNumber,
                           is_jll::Bool)
    if is_jll
      result = string("Your `new _jll package` pull request met all of the ",
                      "guidelines for auto-merging and is scheduled to ",
                      "be merged in the next round.",
                      "\n\n",
                      _comment_noblock(),
                      _onepointzero_suggestion(suggest_onepointzero, version),
                      "\n<!-- [noblock] -->")
    else
      result = string("Your `new package` pull request met all of the ",
                      "guidelines for auto-merging and is scheduled to ",
                      "be merged when the mandatory waiting period (3 days) has elapsed.",
                      "\n\n",
                      "Since you are registering a new package, ",
                      "please make sure that you have read the ",
                      "package naming guidelines: ",
                      "https://julialang.github.io/Pkg.jl/dev/creating-packages/#Package-naming-guidelines-1",
                      "\n\n",
                      _comment_noblock(),
                      _onepointzero_suggestion(suggest_onepointzero, version),
                      "\n<!-- [noblock] -->")
    end
    return result
end

function comment_text_fail(::NewPackage,
                           reasons::Vector{String},
                           suggest_onepointzero::Bool,
                           version::VersionNumber)
    reasons_formatted = join(string.("- ", reasons), "\n")
    result = string("Your `new package` pull request does not meet ",
                    "the guidelines for auto-merging. ",
                    "Please make sure that you have read the ",
                    "[General registry README](https://github.com/JuliaRegistries/General/blob/master/README.md). ",
                    "The following guidelines were not met:\n\n",
                    reasons_formatted,
                    _comment_disclaimer(),
                    "\n\n",
                    "Since you are registering a new package, ",
                    "please make sure that you have also read the ",
                    "package naming guidelines: ",
                    "https://julialang.github.io/Pkg.jl/dev/creating-packages/#Package-naming-guidelines-1",
                    "\n\n",
                    _comment_noblock(),
                    _onepointzero_suggestion(suggest_onepointzero, version),
                    "\n<!-- [noblock] -->")
    return result
end

function comment_text_fail(::NewVersion,
                           reasons::Vector{String},
                           suggest_onepointzero::Bool,
                           version::VersionNumber)
    reasons_formatted = join(string.("- ", reasons), "\n")
    result = string("Your `new version` pull request does not meet ",
                    "the guidelines for auto-merging. ",
                    "Please make sure that you have read the ",
                    "[General registry README](https://github.com/JuliaRegistries/General/blob/master/README.md). ",
                    "The following guidelines were not met:\n\n",
                    reasons_formatted,
                    _comment_disclaimer(),
                    _comment_noblock(),
                    _onepointzero_suggestion(suggest_onepointzero, version),
                    "\n<!-- [noblock] -->")
    return result
end

function comment_text_merge_now()
    result = string("The mandatory waiting period has elapsed.\n\n",
                    "Your pull request is ready to merge.\n\n",
                    "I will now merge this pull request.",
                    "\n<!-- [noblock] -->")
    return result
end

is_julia_stdlib(name) = name in julia_stdlib_list()

function julia_stdlib_list()
    stdlib_list = readdir(Pkg.Types.stdlib_dir())
    # Before Julia v1.6 Artifacts.jl isn't a standard library, but
    # we want to include it because JLL packages depend on the empty
    # placeholder https://github.com/JuliaPackaging/Artifacts.jl
    # in older versions for compatibility.
    if VERSION < v"1.6.0"
        push!(stdlib_list, "Artifacts")
    end
    return stdlib_list
end

function now_utc()
    utc = TimeZones.tz"UTC"
    return Dates.now(utc)
end

function _onepointzero_suggestion(suggest_onepointzero::Bool,
                                  version::VersionNumber)
    if suggest_onepointzero && version < v"1.0.0"
        result = string("\n\n---\n",
                        "On a separate note, I see that you are registering ",
                        "a release with a version number of the form ",
                        "`v0.X.Y`.\n\n",
                        "Does your package have a stable public API? ",
                        "If so, then it's time for you to register version ",
                        "`v1.0.0` of your package. ",
                        "(This is not a requirement. ",
                        "It's just a recommendation.)\n\n",
                        "If your package does not yet have a stable public ",
                        "API, then of course you are not yet ready to ",
                        "release version `v1.0.0`.")
        return result
    else
        return ""
    end
end

function time_is_already_in_utc(dt::Dates.DateTime)
    utc = TimeZones.tz"UTC"
    return TimeZones.ZonedDateTime(dt, utc; from_utc = true)
end

"""
    get_all_non_jll_package_names(registry_dir::AbstractString) -> Vector{String}

Given a path to the directory holding a registry, returns the names of all the non-JLL packages
defined in that registry, along with the names of Julia's standard libraries.
"""
function get_all_non_jll_package_names(registry_dir::AbstractString)
    packages = [x["name"] for x in values(TOML.parsefile(joinpath(registry_dir, "Registry.toml"))["packages"])]
    sort!(packages)
    append!(packages, values(RegistryTools.stdlibs()))
    filter!(x -> !endswith(x, "_jll"), packages)
    unique!(packages)
    return packages
end


text:
import Base: +, -

const genericwrappers = (
    :LowerTriangular,
    :UpperTriangular,
    :Hermitian,
    :Symmetric
)

const unittriangularwrappers = (
    (:UnitUpperTriangular, :UpperTriangular),
    (:UnitLowerTriangular, :LowerTriangular)
)

function kernel_generic(ctx, B, J, min_size)
    lin_idx = linear_index(ctx)
    lin_idx > min_size && return nothing
    @inbounds diag_idx = diagind(B)[lin_idx]
    @inbounds B[diag_idx] += J
    return nothing
end

function kernel_unittriangular(ctx, B, J, diagonal_val, min_size)
    lin_idx = linear_index(ctx)
    lin_idx > min_size && return nothing
    @inbounds diag_idx = diagind(B)[lin_idx]
    @inbounds B[diag_idx] = diagonal_val + J
    return nothing
end

for (t1, t2) in unittriangularwrappers
    @eval begin
        function (+)(A::$t1{T, <:AbstractGPUMatrix}, J::UniformScaling) where T
            B = similar(parent(A), typeof(oneunit(T) + J))
            copyto!(B, parent(A))
            min_size = minimum(size(B))
            gpu_call(kernel_unittriangular, B, J, one(eltype(B)), min_size; elements=min_size)
            return $t2(B)
        end

        function (-)(J::UniformScaling, A::$t1{T, <:AbstractGPUMatrix}) where T
            B = similar(parent(A), typeof(J - oneunit(T)))
            B .= .- parent(A)
            min_size = minimum(size(B))
            gpu_call(kernel_unittriangular, B, J, -one(eltype(B)), min_size; elements=min_size)
            return $t2(B)
        end
    end
end

for t in genericwrappers
    @eval begin
        function (+)(A::$t{T, <:AbstractGPUMatrix}, J::UniformScaling) where T
            B = similar(parent(A), typeof(oneunit(T) + J))
            copyto!(B, parent(A))
            min_size = minimum(size(B))
            gpu_call(kernel_generic, B, J, min_size; elements=min_size)
            return $t(B)
        end

        function (-)(J::UniformScaling, A::$t{T, <:AbstractGPUMatrix}) where T
            B = similar(parent(A), typeof(J - oneunit(T)))
            B .= .- parent(A)
            min_size = minimum(size(B))
            gpu_call(kernel_generic, B, J, min_size; elements=min_size)
            return $t(B)
        end
    end
end

# Breaking Hermiticity when adding a complex value to the diagonal
function (+)(A::Hermitian{T,<:AbstractGPUMatrix}, J::UniformScaling{<:Complex}) where T
    B = similar(parent(A), typeof(oneunit(T) + J))
    copyto!(B, parent(A))
    min_size = minimum(size(B))
    gpu_call(kernel_generic, B, J, min_size; elements=min_size)
    return B
end

function (-)(J::UniformScaling{<:Complex}, A::Hermitian{T,<:AbstractGPUMatrix}) where T
    B = similar(parent(A), typeof(J - oneunit(T)))
    B .= .-parent(A)
    min_size = minimum(size(B))
    gpu_call(kernel_generic, B, J, min_size; elements=min_size)
    return B
end

# Finally the generic matrix version
function (+)(A::AbstractGPUMatrix{T}, J::UniformScaling) where T
    B = similar(A, typeof(oneunit(T) + J))
    copyto!(B, A)
    min_size = minimum(size(B))
    gpu_call(kernel_generic, B, J, min_size; elements=min_size)
    return B
end

function (-)(J::UniformScaling, A::AbstractGPUMatrix{T}) where T
    B = similar(A, typeof(J - oneunit(T)))
    B .= .-A
    min_size = minimum(size(B))
    gpu_call(kernel_generic, B, J, min_size; elements=min_size)
    return B
end