text:
function environment(
    f::Function;
    globalconfig::Vector = [],
    threadconfig::Vector = [],
)
    # Save the current settings
    #
    # CPLGetConfigOption() will return the value of the config option, be it
    #     either defined through environment variable, CPLSetConfigOption() or
    #     CPLSetThreadLocalConfigOption() (from the same thread).
    # CPLGetThreadLocalConfigOption() will return the value of the config
    #     option, but only if it has been set with
    #     CPLSetThreadLocalConfigOption()
    #
    # (ref https://github.com/mapbox/rasterio/pull/997#issuecomment-287117289)
    globalsettings = Dict(k => getconfigoption(k) for (k, v) in globalconfig)
    localsettings =
        Dict(k => getthreadconfigoption(k) for (k, v) in threadconfig)
    for (k, v) in threadconfig
        setthreadconfigoption(k, v)
    end
    for (k, v) in globalconfig
        setconfigoption(k, v)
    end

    return try
        f()
    finally
        # Restore previous settings
        for (k, v) in globalsettings
            if v == ""
                clearconfigoption(k)
            else
                setconfigoption(k, v)
            end
        end
        for (k, v) in localsettings
            if v == ""
                clearthreadconfigoption(k)
            else
                setthreadconfigoption(k, v)
            end
        end
    end
end

function executesql(f::Function, dataset::Dataset, args...)
    result = unsafe_executesql(dataset, args...)
    return try
        f(result)
    finally
        releaseresultset(dataset, result)
    end
end

function createfeature(f::Function, featuredefn::FeatureDefn)
    feature = unsafe_createfeature(featuredefn)
    reference(featuredefn)
    # the additional reference & dereference here is to deal with the case
    # where you start by (1) creating a featuredefn (0 references)
    # before (2) using it to create a feature here.
    # if we do not artificially increase the reference, then destroy(feature)
    # will release the featuredefn, when we're going to handle it ourselves
    # later. Therefore we dereference (rather than release) the featuredefn.
    return try
        f(feature)
    finally
        destroy(feature)
        dereference(featuredefn)
    end
end

"""
    addfielddefn!(layer::AbstractFeatureLayer, name, etype::OGRFieldType;
        <keyword arguments>)

Create a new field on a layer.

This function should not be called while there are feature objects in existence
that were obtained or created with the previous layer definition.

Not all drivers support this function. You can query a layer to check if it
supports it with the OLCCreateField capability. Some drivers may only support
this method while there are still no features in the layer. When it is
supported, the existing features of the backing file/database should be updated
accordingly.

Drivers may or may not support not-null constraints. If they support creating
fields with not-null constraints, this is generally before creating any feature
to the layer.

### Parameters
* `layer`:  the layer to write the field definition.
* `name`:   name of the field definition to write to disk.
* `etype`:  type of the field definition to write to disk.

### Keyword arguments
* `nwidth`:     the preferred formatting width. 0 (default) indicates undefined.
* `nprecision`: number of decimals for formatting. 0 (default) for undefined.
* `justify`:    the formatting justification ([OJUndefined], OJLeft or OJRight)
* `approx`:     If `true` (default `false`), the field may be created in a
                slightly different form depending on the limitations of the
                format driver.
"""
function addfielddefn!(
    layer::T,
    name::AbstractString,
    etype::OGRFieldType;
    nwidth::Integer = 0,
    nprecision::Integer = 0,
    justify::OGRJustification = OJUndefined,
    approx::Bool = false,
)::T where {T<:AbstractFeatureLayer}
    fielddefn = unsafe_createfielddefn(name, etype)
    setparams!(
        fielddefn,
        name,
        etype,
        nwidth = nwidth,
        nprecision = nprecision,
        justify = justify,
    )
    addfielddefn!(layer, fielddefn)
    destroy(fielddefn)
    return layer
end

"""
    writegeomdefn!(layer::AbstractFeatureLayer, name, etype::OGRwkbGeometryType,
        approx=false)

Write a new geometry field on a layer.

This function should not be called while there are feature objects in existence
that were obtained or created with the previous layer definition.

Not all drivers support this function. You can query a layer to check if it
supports it with the OLCCreateField capability. Some drivers may only support
this method while there are still no features in the layer. When it is
supported, the existing features of the backing file/database should be updated
accordingly.

Drivers may or may not support not-null constraints. If they support creating
fields with not-null constraints, this is generally before creating any feature
to the layer.

### Parameters
* `layer`:  the layer to write the field definition.
* `name`:   name of the field definition to write to disk.
* `etype`:  type of the geometry field defintion to write to disk.

### Keyword arguments
* `approx`: If `true` (default `false`), the geometry field may be created in a
            slightly different form depending on the limitations of the driver.
"""
function writegeomdefn!(
    layer::T,
    name::AbstractString,
    etype::OGRwkbGeometryType;
    approx::Bool = false,
)::T where {T<:AbstractFeatureLayer}
    geomdefn = unsafe_creategeomdefn(name, etype)
    addgeomdefn!(layer, geomdefn)
    destroy(geomdefn)
    return layer
end

function writegeomdefn(
    f::Function,
    layer::AbstractFeatureLayer,
    name::AbstractString,
    etype::OGRwkbGeometryType;
    approx::Bool = false,
)
    geomdefn = unsafe_creategeomdefn(name, etype)
    return try
        f(geomdefn)
        addgeomdefn!(layer, geomdefn)
    finally
        destroy(geomdefn)
    end
end

for gdalfunc in (
    :boundary,
    :buffer,
    :centroid,
    :clone,
    :convexhull,
    :create,
    :createcolortable,
    :createcoordtrans,
    :copy,
    :createfeaturedefn,
    :createfielddefn,
    :creategeom,
    :creategeomcollection,
    :creategeomfieldcollection,
    :creategeomdefn,
    :createlayer,
    :createlinearring,
    :createlinestring,
    :createmultilinestring,
    :createmultipoint,
    :createmultipolygon,
    :createmultipolygon_noholes,
    :createpoint,
    :createpolygon,
    :createRAT,
    :createstylemanager,
    :createstyletable,
    :createstyletool,
    :curvegeom,
    :delaunaytriangulation,
    :difference,
    :forceto,
    :fromGML,
    :fromJSON,
    :fromWKB,
    :fromWKT,
    :gdalbuildvrt,
    :gdaldem,
    :gdalgrid,
    :gdalnearblack,
    :gdalrasterize,
    :gdaltranslate,
    :gdalvectortranslate,
    :gdalwarp,
    :getband,
    :getcolortable,
    :getfeature,
    :getgeom,
    :getlayer,
    :getmaskband,
    :getoverview,
    :getpart,
    :getspatialref,
    :importCRS,
    :intersection,
    :importEPSG,
    :importEPSGA,
    :importESRI,
    :importPROJ4,
    :importWKT,
    :importXML,
    :importURL,
    :lineargeom,
    :newspatialref,
    :nextfeature,
    :pointalongline,
    :pointonsurface,
    :polygonfromedges,
    :polygonize,
    :read,
    :sampleoverview,
    :simplify,
    :simplifypreservetopology,
    :symdifference,
    :union,
    :update,
    :readraster,
)
    eval(quote
        function $(gdalfunc)(f::Function, args...; kwargs...)
            obj = $(Symbol("unsafe_$gdalfunc"))(args...; kwargs...)
            return try
                f(obj)
            finally
                destroy(obj)
            end
        end
    end)
end


text:
using Revise
using Lava
using VulkanCore
using LinearAlgebra
using GeometryTypes
using Lava: @class, @scope, @autodestroy

shaderFolder = String(@__DIR__) * "/../shaders/"

mutable struct CameraData
    view::Mat4f0
    proj::Mat4f0

    CameraData() = new(Mat4f0(1I), Mat4f0(1I))
    CameraData(view::Mat4f0, proj::Mat4f0) = new(view, proj)
end

@class CameraData

function getTransposeData(data::CameraData)
    return [transpose(data.view) transpose(data.proj)]
end

const u8vec4 = Vec{4, UInt8}
struct Vertex
    position::Vec3f0
    color::u8vec4
end

cubeVertices = [
    Vertex(Vec3f0(-1.0, -1.0, -1.0), u8vec4(0, 0, 0, 255)),
    Vertex(Vec3f0(-1.0, -1.0, 1.0), u8vec4(0, 0, 255, 255)),
    Vertex(Vec3f0(-1.0, 1.0, -1.0), u8vec4(0, 255, 0, 255)),
    Vertex(Vec3f0(-1.0, 1.0, 1.0), u8vec4(0, 255, 255, 255)),
    Vertex(Vec3f0(1.0, -1.0, -1.0), u8vec4(255, 0, 0, 255)),
    Vertex(Vec3f0(1.0, -1.0, 1.0), u8vec4(255, 0, 255, 255)),
    Vertex(Vec3f0(1.0, 1.0, -1.0), u8vec4(255, 255, 0, 255)),
    Vertex(Vec3f0(1.0, 1.0, 1.0), u8vec4(255, 255, 255, 255))
]

cubeIndices = Vector{UInt32}(  [0, 2, 3,
                                0, 3, 1,
                                0, 1, 5,
                                0, 5, 4,
                                2, 7, 3,
                                2, 6, 7,
                                1, 7, 5,
                                1, 3, 7,
                                0, 4, 6,
                                0, 6, 2,
                                4, 5, 7,
                                4, 7, 6])

function main()
    # Create a Vulkan instance, tell it we need glfw and the validation as
    # extension features
    @scope begin
    fs = Vector{features.IFeature}()
    glfw = features.create(features.GlfwOutput)
    push!(fs, features.create(features.Validation))
    push!(fs, glfw)
    instance = @autodestroy lava.create(lava.Instance, fs)
    queues = [lava.createGraphics(lava.QueueRequest, "graphics")]#::Vector{QueueRequest}

    # Create a logical device with a single graphics queue named "graphics" on
    # a discrete (non-integrated) GPU. For more complex apps you can request
    # additional queues (e.g. a separate transfer queue)
    device = @autodestroy instance.createDevice(queues,
                                    lava.NthOfTypeStrategy(vk.VK_PHYSICAL_DEVICE_TYPE_DISCRETE_GPU))
                                   # lava.NthOfTypeStrategy(vk.VK_PHYSICAL_DEVICE_TYPE_INTEGRATED_GPU

    # Create a pipeline layout that transfers the CameraData as push constants
    # and has no descriptors (textures or uniform buffers)
    plLayout = @autodestroy device.createPipelineLayout(CameraData)

    # Test descriptor set layout
    # descCreateInfo = lava.DescriptorSetLayoutCreateInfo()
    # lava.addUniformBuffer(descCreateInfo, vk.VK_SHADER_STAGE_VERTEX_BIT);
    # lava.addCombinedImageSampler(descCreateInfo, vk.VK_SHADER_STAGE_FRAGMENT_BIT);
    # descLayout = lava.createDescriptorSetLayout(device, descCreateInfo)
    # plLayout = lava.createPipelineLayout(device, CameraData, [descLayout])

    pass = @autodestroy device.createRenderPass(lava.createSimpleForward(lava.RenderPassCreateInfo, glfw.format()))

    vertShader = @autodestroy device.createShaderFromFile(shaderFolder * "cube_vert.spv")
    fragShader = @autodestroy device.createShaderFromFile(shaderFolder * "cube_frag.spv")
    stageVert = lava.defaults(lava.PipelineShaderStageCreateInfo, _module = vertShader)
    stageFrag = lava.defaults(lava.PipelineShaderStageCreateInfo, _module = fragShader)
    stages = [stageVert.handleRef()[], stageFrag.handleRef()[]]

    attributes = Vector{vk.VkVertexInputAttributeDescription}()
    bindings = Vector{vk.VkVertexInputBindingDescription}()
    lava.binding(attributes, bindings, UInt32(0), Vertex, :position)
    lava.binding(attributes, bindings, UInt32(0), Vertex, :color)
    vertexInputState = lava.PipelineVertexInputStateCreateInfo(attributes = attributes, bindings = bindings)

    ci = lava.defaults(lava.GraphicsPipelineCreateInfo,
        stages = stages,
        layout = plLayout,
        renderPass = lava.handleRef(pass)[],
        subpass = UInt32(0),
        depthTestEnable = VkExt.VK_TRUE,
        depthWriteEnable = VkExt.VK_TRUE,
        # Due to the flipped y axis in NDC, the triangle winding order is inverted, too
        frontFace = vk.VK_FRONT_FACE_CLOCKWISE,
        inputState = vertexInputState
    )
    pipeline = @autodestroy lava.GraphicsPipeline(ci)

    camera = Camera.GenericCamera()
    camera.setTarget(Vec3f0(0.0, 0.0, 0.0))
    fbos = @autodestroy Vector{lava.Framebuffer}()
    window = @autodestroy glfw.openWindow(UInt32(800), UInt32(600), true)
    imgList = @autodestroy []
    imgViewList = @autodestroy []
    window.buildSwapchainWith(views::Vector{lava.ImageView}->begin
        depthImageCreateInfo = lava.attachment2D(lava.getPhysicalDevice(device), window.mWidth, window.mHeight, vk.VK_FORMAT_D32_SFLOAT)
        depthImage = lava.createImage(depthImageCreateInfo, device)
        depthImage.realizeAttachment()
        depthView = depthImage.createView()
        lava.destroy!(imgList)
        lava.destroy!(imgViewList)
        push!(imgList, depthImage)
        push!(imgViewList, depthView)

        lava.destroy!(fbos)
        for view in views
            push!(fbos, lava.Framebuffer(pass, [depthView, view]))
        end

        camera.setAspectRatio(Float32(features.getWidth(window)) / Float32(features.getHeight(window)))
    end)

    #Upload of the cube mesh
    vab = @autodestroy lava.createBuffer(device, lava.arrayBuffer())
    vab.setDataVRAM(cubeVertices, Vertex)
    eab = @autodestroy lava.createBuffer(device, lava.indexBuffer())
    eab.setDataVRAM(cubeIndices, UInt32)

    camera.setPosition(Vec3f0(2.0, 2.0, 2.0))
    camera.setTarget(Vec3f0(0.0, 0.0, 0.0))
    camera.rotateAroundTarget_GlobalAxes(Float32(0.0), Float32(0.001), Float32(0.0))

    while true
        features.pollEvent()
        lava.graphicsQueue(device).catchUp(Int32(1))
        camera.rotateAroundTarget_GlobalAxes(Float32(0.0), Float32(0.001), Float32(0.0))
        #println(CameraData(Camera.getViewMatrix(camera), Camera.getProjectionMatrix(camera)))
        frame = window.startFrame()
        cmd = lava.graphicsQueue(device).beginCommandBuffer()

        # Vulkan doesn't take care of synchronization for us.
        # We need to tell it to wait with rendering until the image is ready
        # and to wait with presenting the image until the rendering is done.
        # The frame struct provides the necessary semaphores.
        cmd.wait(frame.imageReady())
        cmd.signal(frame.renderingComplete())

        forward = cmd.beginRenderpass(fbos[features.imageIndex(frame) + 1])

        sub = forward.startInlineSubpass()
        sub.bindPipeline(pipeline)
        sub.setViewports([vk.VkViewport(0, 0, features.getWidth(window), features.getHeight(window), 0.0, 1.0)])
        sub.bindVertexBuffers([vab])
        sub.bindIndexBuffer(eab)
        cameraData = CameraData(camera.getViewMatrix(), camera.getProjectionMatrix())
        sub.pushConstantBlock(cameraData.getTransposeData())

        sub.drawIndexed(UInt32(length(cubeIndices)))

        forward.endRenderPass()
        cmd.endCommandBuffer()
        frame.endFrame()
    end
    end
end

ext: jl


text:
__precompile__()

module OCReract

export
    run_tesseract,
    run_and_get_output

include("tesseract.jl")

"""
OCReract is a simple Julia wrapper of the well-known OCR engine called Tesseract.

Here, a simple example of usage:

# Example
```julia-repl
julia> using Images;
julia> using OCReract;
julia> img_path = "/path/to/img.png";
julia> img = Images.load(img_path);
julia> res_text = run_and_get_output(img, psm=3, oem=1);
julia> println(res_text);
```

For more information, check the homepage in https://github.com/leferrad/OCReract.jl.
"""
OCReract

end # module


text:
"""
Several genetic operators chained together.
"""
struct OperatorPipeline <: GeneticOperator
    ops::Vector{GeneticOperator}

    function OperatorPipeline(operators...)
        ensure_nargs_match(operators) || throw("The number of ")
        new(operators)
    end
end

function apply(p::OperatorPipeline, parents::AbstractVector{<:AbstractVector{<:Real}})
    for i in 1:length(p.ops)
        parents = apply(p.ops[i], parents)
    end
    return parents
end

match(o1::GeneticOperator, o2::GeneticOperator) = numchildren(o1) == numparents(o2)
# Mutation ops can accept any number of parents
match(o1::GeneticOperator, o2::MutationOperator) = true

function ensure_nargs_match(operators)
    for i in 1:(length(operators)-1)
        match(operators[i], operators[i+1]) || return false
    end
    return true
end


text:
struct GrammarError <: Exception
    message::String
end

# These are concrete to allow messages to be attached
# to them.
struct ParseError <: Exception end
struct LexError <: Exception
    message::String
end

struct UnexpectedEOF <: Exception
    expected::Array
end

Base.showerror(io::IO, ue::UnexpectedEOF) = begin
    println("Unexpected end-of-input. Expected one of:"* join(expected,"\t\t*")*"\n")
end

#==
The UnexpectedInput type holds the common behaviour of
the UnexpectedCharacters and UnexpectedToken types
==#

abstract type UnexpectedInput <: Exception end

# TODO: are we accessing strings in a Unicode-aware way?
get_context(ui::UnexpectedInput,text;span=40) = begin
    pos = ui.pos_in_stream
    start = max(pos - span,1)
    _end = min(pos + span,length(text))
    before = rsplit(text[start:pos],"\n",limit=2)[end]
    after = split(text[pos+1:_end],"\n",limit=2)[1]
    return before*after*"\n"*repeat(" ",length(before)-1) * "^\n"
end

"""
    Given a parser instance and a dictionary mapping some label with
    some malformed syntax examples, it'll return the label for the
    example that bests matches the current error.

For Julia we need to make sure the correct parse_fn is dispatched for
this to make sense.
"""
match_examples(ui::UnexpectedInput,parse_call,examples::Dict,token_type_match_fallback=false,use_accepts=false) = begin
    match_examples(ui,parse_call,((a,b) for (a,b) in examples),
                   token_type_match_fallback=token_type_match_fallback,
                   use_accepts=use_accepts)
end

match_examples(ui::UnexpectedInput,parse_call,examples,token_type_match_fallback=false,use_accepts=false) = begin
    candidate = (nothing,false)
    for (i,(label,example)) in enumerate(examples)
        @assert !(example <: AbstractString)
        for (j,malformed) in enumerate(example)
            try
                parse_call(malformed)
            catch ut
                if ut isa UnexpectedInput
                    if ut.state == ui.state
                        if use_accepts && ut.accepts != ui.accepts
                            @debug "Different accepts with the same state[$(ui.state)]: $(ui.accepts) != $(ut.accepts) at example [$i][$j]"
                            continue
                        end
                        try
                            if ut.token == ui.token
                                @debug "Exact match at example [$i][$j]"
                                return label
                            end
                            if token_type_match_fallback
                                if (ut.token.type_ == ui.token.type_) && candidate[end]==nothing
                                    @debug "Token type fallback at example [$i][$j]"
                                    candidate = label,true
                                end
                            end
                        catch a
                            if !(a isa AttributeError)
                                rethrow(a)
                            end
                        end
                        if candidate[1] == nothing
                            @debug "Same state match at example [$i][$j]"
                            candidate = label,false
                        end
                    end
                else
                    rethrow(ut)
                end
            end
        end
    end
    return candidate[1]
end

"""
    UnexpectedCharacters

This exception is raised when the next characters in the input
stream cannot be matched with any acceptable tokens.
"""
struct UnexpectedCharacters <: UnexpectedInput
    seq
    pos_in_stream
    line
    column
    allowed
    considered_tokens
    state
    token_history
end

UnexpectedCharacters(seq,lex_pos,line,column;allowed=nothing,consider_tokens=nothing,state=nothing,token_history=nothing) = UnexpectedCharacters(seq,lex_pos,line,column,allowed,consider_tokens,state,token_history)

Base.showerror(io::IO,uc::UnexpectedCharacters) = begin
    s = uc.seq[uc.pos_in_stream]
    message = "No terminal defined for $s at line $(uc.line) col $(uc.column)"
    message = message * "\n\n" * get_context(uc,uc.seq)
    if uc.allowed != nothing
        message *= "\nExpecting: $(uc.allowed)"
    end
    if uc.token_history != nothing
        message *= "\nPrevious tokens: " * join(uc.token_history,", ")*"\n"
    end
    println(io,message)
end

"""
    UnexpectedToken

This exception is raised when a token is encountered in the input stream
that is not allowed by the grammar.
"""
struct UnexpectedToken <: UnexpectedInput
    token
    expected
    line
    column
    considered_rules
    state
    pos_in_stream
    puppet
    accepts
end

UnexpectedToken(token,expected;considered_rules=nothing,state=nothing,puppet=nothing) = begin
    line = token.line
    column = token.column
    pos_in_stream = token.pos_in_stream
    accepts = puppet != nothing && accepts(puppet)
    UnexpectedToken(token,expected,line,column,considered_rules,state,pos_in_stream,puppet,accepts)
end

Base.showerror(io::IO,ut::UnexpectedToken) = begin
    message = "Unexpected token $(ut.token) at line $(ut.line), column $(ut.column).\n"
    message *= "Expected one of: $(join(ut.expected,", "))"
    println(message)
end

struct VisitError <: Exception
    tree
    orig_exc
end

Base.showerror(io::IO,ve::VisitError) = begin
    println("Error trying to process rule \"$(ve.tree.data):\n\n$(ve.orig_exc)\" ")
end


text:
struct BS3ConstantCache{T,T2} <: OrdinaryDiffEqConstantCache
  a21::T
  a32::T
  a41::T
  a42::T
  a43::T
  c1::T2
  c2::T2
  btilde1::T
  btilde2::T
  btilde3::T
  btilde4::T
end

"""
constructBogakiShampine3()

Constructs the tableau object for the Bogakai-Shampine Order 2/3 method.
"""
function BS3ConstantCache(::Type{T},::Type{T2}) where {T<:CompiledFloats,T2<:CompiledFloats}
  a21 = T(0.5)
  a32 = T(0.75)
  a41 = T(0.2222222222222222)
  a42 = T(0.3333333333333333)
  a43 = T(0.4444444444444444)
  c1 = T2(0.5)
  c2 = T2(0.75)
  # b1 = T(0.2916666666666667)
  # b2 = T(0.25)
  # b3 = T(0.3333333333333333)
  # b4 = T(0.125)
  btilde1 = T(0.06944444444444445)
  btilde2 = T(-0.08333333333333333)
  btilde3 = T(-0.1111111111111111)
  btilde4 = T(0.125)
  BS3ConstantCache(a21,a32,a41,a42,a43,c1,c2,btilde1,btilde2,btilde3,btilde4)
end

"""
constructBogakiShampine3()

Constructs the tableau object for the Bogakai-Shampine Order 2/3 method.
"""
function BS3ConstantCache(T::Type,T2::Type)
  a21 = T(1//2)
  a32 = T(3//4)
  a41 = T(2//9)
  a42 = T(1//3)
  a43 = T(4//9)
  c1 = T2(1//2)
  c2 = T2(3//4)
  # b1 = T(7//24)
  # b2 = T(1//4)
  # b3 = T(1//3)
  # b4 = T(1//8)
  btilde1 = T(5//72)
  btilde2 = T(-1//12)
  btilde3 = T(-1//9)
  btilde4 = T(1//8)
  BS3ConstantCache(a21,a32,a41,a42,a43,c1,c2,btilde1,btilde2,btilde3,btilde4)
end

struct OwrenZen3ConstantCache{T,T2} <: OrdinaryDiffEqConstantCache
  a21::T
  a31::T
  a32::T
  a41::T
  a42::T
  a43::T
  c1::T2
  c2::T2
  btilde1::T
  btilde2::T
  btilde3::T
  r13::T
  r12::T
  r23::T
  r22::T
  r33::T
  r32::T
end

function OwrenZen3ConstantCache(::Type{T},::Type{T2}) where {T<:CompiledFloats,T2<:CompiledFloats}
  a21 = T(0.5217391304347826)
  a31 = T(-0.18133333333333335)
  a32 = T(0.9813333333333333)
  a41 = T(0.2152777777777778)
  a42 = T(0.4592013888888889)
  a43 = T(0.3255208333333333)
  c1 = T2(0.5217391304347826)
  c2 = T2(0.8)
  # b1 = T(0.041666666666666664)
  # b2 = T(0.9583333333333334)
  btilde1 = T(-0.1736111111111111)
  btilde2 = T(0.4991319444444444)
  btilde3 = T(-0.3255208333333333)
  r13 = T(0.5694444444444444)
  r12 = T(-1.3541666666666667)
  r23 = T(-0.9184027777777778)
  r22 = T(1.3776041666666667)
  r33 = T(-0.6510416666666666)
  r32 = T(0.9765625)
  OwrenZen3ConstantCache(a21,a31,a32,a41,a42,a43,c1,c2,btilde1,btilde2,btilde3,
                         r13,r12,r23,r22,r33,r32)
end

function OwrenZen3ConstantCache(::Type{T},::Type{T2}) where {T,T2}
  a21 = T(12//23)
  a31 = T(-68//375)
  a32 = T(368//375)
  a41 = T(31//144)
  a42 = T(529//1152)
  a43 = T(125//384)
  c1 = T2(12//23)
  c2 = T2(4//5)
  # b1 = T(1//24)
  # b2 = T(23//24)
  btilde1 = T(-25//144)
  btilde2 = T(575//1152)
  btilde3 = T(-125//384)
  r13 = T(41//72)
  r12 = T(-65//48)
  r23 = T(-529//576)
  r22 = T(529//384)
  r33 = T(-125//192)
  r32 = T(125//128)
  OwrenZen3ConstantCache(a21,a31,a32,a41,a42,a43,c1,c2,btilde1,btilde2,btilde3,
                         r13,r12,r23,r22,r33,r32)
end

struct OwrenZen4ConstantCache{T,T2} <: OrdinaryDiffEqConstantCache
  a21::T
  a31::T
  a32::T
  a41::T
  a42::T
  a43::T
  a51::T
  a52::T
  a53::T
  a54::T
  a61::T
  a63::T
  a64::T
  a65::T
  c1::T2
  c2::T2
  c3::T2
  c4::T2
  btilde1::T
  btilde3::T
  btilde4::T
  btilde5::T
  r14::T
  r13::T
  r12::T
  r34::T
  r33::T
  r32::T
  r44::T
  r43::T
  r42::T
  r54::T
  r53::T
  r52::T
  r64::T
  r63::T
  r62::T
end

function OwrenZen4ConstantCache(::Type{T},::Type{T2}) where {T<:CompiledFloats,T2<:CompiledFloats}
    a21 = T(0.16666666666666666)
    a31 = T(0.03214024835646457)
    a32 = T(0.26515704894083275)
    a41 = T(0.6895990230002036)
    a42 = T(-1.6993690209647874)
    a43 = T(1.6568288214939955)
    a51 = T(-0.09002509947964493)
    a52 = T(0.6817777777777778)
    a53 = T(-0.2402791551882461)
    a54 = T(0.5151931435567799)
    a61 = T(0.08990252172070354)
    a63 = T(0.4360623278236915)
    a64 = T(0.1842858372687918)
    a65 = T(0.2897493131868132)
    c1 = T2(0.16666666666666666)
    c2 = T2(0.2972972972972973)
    c3 = T2(0.6470588235294118)
    c4 = T2(0.8666666666666667)
    # b1 = T(0.27823691460055094)
    # b3 = T(-0.09428374655647383)
    # b4 = T(0.8160468319559229)
    btilde1 = T(0.18833439287984743)
    btilde3 = T(-0.5303460743801653)
    btilde4 = T(0.6317609946871311)
    btilde5 = T(-0.2897493131868132)
    r14 = T(-1.0513495872621479)
    r13 = T(2.922894131082889)
    r12 = T(-2.7816420221000375)
    r34 = T(2.4266369235379472)
    r33 = T(-5.725398502723277)
    r32 = T(3.7348239070090217)
    r44 = T(-0.7705135914137909)
    r43 = T(1.1724555082899983)
    r42 = T(-0.21765607960741548)
    r54 = T(-2.8643157295948325)
    r53 = T(5.149132832816039)
    r52 = T(-1.9950677900343932)
    r64 = T(2.2595419847328246)
    r63 = T(-3.519083969465649)
    r62 = T(1.2595419847328244)
    OwrenZen4ConstantCache(a21,a31,a32,a41,a42,a43,a51,a52,a53,a54,
                           a61,a63,a64,a65,c1,c2,c3,c4,btilde1,btilde3,btilde4,btilde5,
                           r14,r13,r12,r34,r33,r32,r44,r43,r42,
                           r54,r53,r52,r64,r63,r62)
end

function OwrenZen4ConstantCache(::Type{T},::Type{T2}) where {T,T2}
  a21 = T(1//6)
  a31 = T(44//1369)
  a32 = T(363//1369)
  a41 = T(3388//4913)
  a42 = T(-8349//4913)
  a43 = T(8140//4913)
  a51 = T(-36764//408375)
  a52 = T(767//1125)
  a53 = T(-32708//136125)
  a54 = T(210392//408375)
  a61 = T(1697//18876)
  a63 = T(50653//116160)
  a64 = T(299693//1626240)
  a65 = T(3375//11648)
  c1 = T2(1//6)
  c2 = T2(11//37)
  c3 = T2(11//17)
  c4 = T2(13//15)
  # b1 = T(101//363)
  # b3 = T(-1369//14520)
  # b4 = T(11849//14520)
  btilde1 = T(1185//6292)
  btilde3 = T(-4107//7744)
  btilde4 = T(68493//108416)
  btilde5 = T(-3375//11648)
  r14 = T(-866577//824252)
  r13 = T(1806901//618189)
  r12 = T(-104217//37466)
  r34 = T(12308679//5072320)
  r33 = T(-2178079//380424)
  r32 = T(861101//230560)
  r44 = T(-7816583//10144640)
  r43 = T(6244423//5325936)
  r42 = T(-63869//293440)
  r54 = T(-624375//217984)
  r53 = T(982125//190736)
  r52 = T(-1522125//762944)
  r64 = T(296//131)
  r63 = T(-461//131)
  r62 = T(165//131)
  OwrenZen4ConstantCache(a21,a31,a32,a41,a42,a43,a51,a52,a53,a54,
                         a61,a63,a64,a65,c1,c2,c3,c4,btilde1,btilde3,btilde4,btilde5,
                         r14,r13,r12,r34,r33,r32,r44,r43,r42,
                         r54,r53,r52,r64,r63,r62)
end

struct OwrenZen5ConstantCache{T,T2} <: OrdinaryDiffEqConstantCache
  a21::T
  a31::T
  a32::T
  a41::T
  a42::T
  a51::T
  a52::T
  a53::T
  a54::T
  a61::T
  a62::T
  a63::T
  a64::T
  a65::T
  a71::T
  a72::T
  a73::T
  a74::T
  a75::T
  a76::T
  a81::T
  a83::T
  a84::T
  a85::T
  a86::T
  a87::T
  c1::T2
  c2::T2
  c3::T2
  c4::T2
  c5::T2
  c6::T2
  btilde1::T
  btilde3::T
  btilde4::T
  btilde5::T
  btilde6::T
  btilde7::T
  r15::T
  r14::T
  r13::T
  r12::T
  r35::T
  r34::T
  r33::T
  r32::T
  r45::T
  r44::T
  r43::T
  r42::T
  r55::T
  r54::T
  r53::T
  r52::T
  r65::T
  r64::T
  r63::T
  r62::T
  r75::T
  r74::T
  r73::T
  r72::T
  r85::T
  r84::T
  r83::T
  r82::T
end

function OwrenZen5ConstantCache(::Type{T},::Type{T2}) where {T<:CompiledFloats,T2<:CompiledFloats}
    a21 = T(0.16666666666666666)
    a31 = T(0.0625)
    a32 = T(0.1875)
    a41 = T(0.25)
    a42 = T(-0.75)
    a51 = T(-0.75)
    a52 = T(3.75)
    a53 = T(-3.0)
    a54 = T(0.5)
    a61 = T(0.26895043731778423)
    a62 = T(-0.7084548104956269)
    a63 = T(0.8658892128279884)
    a64 = T(0.15462307371928363)
    a65 = T(0.06184922948771345)
    a71 = T(-0.02947695035460993)
    a72 = T(0.18500664893617022)
    a73 = T(0.4802345261121857)
    a74 = T(-0.5337849069148937)
    a75 = T(-0.013090093085106383)
    a76 = T(0.7861107753062541)
    a81 = T(0.08783068783068783)
    a83 = T(0.3006060606060606)
    a84 = T(0.22777777777777777)
    a85 = T(0.027777777777777776)
    a86 = T(0.06218596218596219)
    a87 = T(0.2938217338217338)
    c1 = T2(0.16666666666666666)
    c2 = T2(0.25)
    c3 = T2(0.5)
    c4 = T2(0.5)
    c5 = T2(0.6428571428571429)
    c6 = T2(0.875)
    # b1 = T(-0.1111111111111111)
    # b3 = T(1.2121212121212122)
    # b4 = T(-1.75)
    # b5 = T(-0.08333333333333333)
    # b6 = T(1.7323232323232323)
    btilde1 = T(-0.19894179894179895)
    btilde3 = T(0.9115151515151515)
    btilde4 = T(-1.9777777777777779)
    btilde5 = T(-0.1111111111111111)
    btilde6 = T(1.67013727013727)
    btilde7 = T(-0.2938217338217338)
    r15 = T(1.892063492063492)
    r14 = T(-6.067155067155067)
    r13 = T(7.282458282458283)
    r12 = T(-4.0195360195360195)
    r35 = T(-7.214545454545455)
    r34 = T(20.676923076923078)
    r33 = T(-20.31142191142191)
    r32 = T(7.14965034965035)
    r45 = T(7.866666666666666)
    r44 = T(-18.78205128205128)
    r43 = T(13.508547008547009)
    r42 = T(-2.3653846153846154)
    r55 = T(2)
    r54 = T(-5.294871794871795)
    r53 = T(4.534188034188034)
    r52 = T(-1.2115384615384615)
    r65 = T(-1.4924630924630924)
    r64 = T(1.5785667324128863)
    r63 = T(1.1958838881915805)
    r62 = T(-1.219801565955412)
    r75 = T(-7.051721611721612)
    r74 = T(16.273203719357564)
    r73 = T(-11.978886071193763)
    r72 = T(3.0512256973795435)
    r85 = T(4)
    r84 = T(-8.384615384615385)
    r83 = T(5.769230769230769)
    r82 = T(-1.3846153846153846)
    OwrenZen5ConstantCache(a21,a31,a32,a41,a42,a51,a52,a53,
                           a54,a61,a62,a63,a64,a65,a71,a72,a73,
                           a74,a75,a76,a81,a83,a84,a85,a86,a87,
                           c1,c2,c3,c4,c5,c6,btilde1,btilde3,btilde4,btilde5,btilde6,btilde7,
                           r15,r14,r13,r12,r35,r34,r33,r32,r45,r44,
                           r43,r42,r55,r54,r53,r52,r65,r64,r63,
                           r62,r75,r74,r73,r72,r85,r84,r83,r82)
end

function OwrenZen5ConstantCache(::Type{T},::Type{T2}) where {T,T2}
  a21 = T(1//6)
  a31 = T(1//16)
  a32 = T(3//16)
  a41 = T(1//4)
  a42 = T(-3//4)
  a43 = T(1)
  a51 = T(-3//4)
  a52 = T(15//4)
  a53 = T(-3)
  a54 = T(1//2)
  a61 = T(369//1372)
  a62 = T(-243//343)
  a63 = T(297//343)
  a64 = T(1485//9604)
  a65 = T(297//4802)
  a71 = T(-133//4512)
  a72 = T(1113//6016)
  a73 = T(7945//16544)
  a74 = T(-12845//24064)
  a75 = T(-315//24064)
  a76 = T(156065//198528)
  a81 = T(83//945)
  a83 = T(248//825)
  a84 = T(41//180)
  a85 = T(1//36)
  a86 = T(2401//38610)
  a87 = T(6016//20475)
  c1 = T2(1//6)
  c2 = T2(1//4)
  c3 = T2(1//2)
  c4 = T2(1//2)
  c5 = T2(9//14)
  c6 = T2(7//8)
  # b1 = T(-1//9)
  # b3 = T(40//33)
  # b4 = T(-7//4)
  # b5 = T(-1//12)
  # b6 = T(343//198)
  btilde1 = T(-188//945)
  btilde3 = T(752//825)
  btilde4 = T(-89//45)
  btilde5 = T(-1//9)
  btilde6 = T(32242//19305)
  btilde7 = T(-6016//20475)
  r15 = T(596//315)
  r14 = T(-4969//819)
  r13 = T(17893//2457)
  r12 = T(-3292//819)
  r35 = T(-1984//275)
  r34 = T(1344//65)
  r33 = T(-43568//2145)
  r32 = T(5112//715)
  r45 = T(118//15)
  r44 = T(-1465//78)
  r43 = T(3161//234)
  r42 = T(-123//52)
  r55 = T(2)
  r54 = T(-413//78)
  r53 = T(1061//234)
  r52 = T(-63//52)
  r65 = T(-9604//6435)
  r64 = T(2401//1521)
  r63 = T(60025//50193)
  r62 = T(-40817//33462)
  r75 = T(-48128//6825)
  r74 = T(96256//5915)
  r73 = T(-637696//53235)
  r72 = T(18048//5915)
  r85 = T(4)
  r84 = T(-109//13)
  r83 = T(75//13)
  r82 = T(-18//13)
  OwrenZen5ConstantCache(a21,a31,a32,a41,a42,a51,a52,a53,
                         a54,a61,a62,a63,a64,a65,a71,a72,a73,
                         a74,a75,a76,a81,a83,a84,a85,a86,a87,
                         c1,c2,c3,c4,c5,c6,btilde1,btilde3,btilde4,btilde5,btilde6,btilde7,
                         r15,r14,r13,r12,r35,r34,r33,r32,r45,r44,
                         r43,r42,r55,r54,r53,r52,r65,r64,r63,
                         r62,r75,r74,r73,r72,r85,r84,r83,r82)
end

struct Tsit5ConstantCache{T,T2} <: OrdinaryDiffEqConstantCache
  c1::T2
  c2::T2
  c3::T2
  c4::T2
  c5::T2
  c6::T2
  a21::T
  a31::T
  a32::T
  a41::T
  a42::T
  a43::T
  a51::T
  a52::T
  a53::T
  a54::T
  a61::T
  a62::T
  a63::T
  a64::T
  a65::T
  a71::T
  a72::T
  a73::T
  a74::T
  a75::T
  a76::T
  btilde1::T
  btilde2::T
  btilde3::T
  btilde4::T
  btilde5::T
  btilde6::T
  btilde7::T
  r11::T
  r12::T
  r13::T
  r14::T
  r22::T
  r23::T
  r24::T
  r32::T
  r33::T
  r34::T
  r42::T
  r43::T
  r44::T
  r52::T
  r53::T
  r54::T
  r62::T
  r63::T
  r64::T
  r72::T
  r73::T
  r74::T
end

function Tsit5ConstantCache(::Type{T},::Type{T2}) where {T<:CompiledFloats,T2<:CompiledFloats}
  c1 =       T2(0.161)
  c2 =       T2(0.327)
  c3 =       T2(0.9)
  c4 =       T2(0.9800255409045097)
  c5 =       T2(1)
  c6 =       T2(1)
  a21=       T(0.161)
  a31=       T(-0.008480655492356989)
  a32=       T(0.335480655492357)
  a41=       T(2.8971530571054935)
  a42=       T(-6.359448489975075)
  a43=       T(4.3622954328695815)
  a51=       T(5.325864828439257)
  a52=       T(-11.748883564062828)
  a53=       T(7.4955393428898365)
  a54=       T(-0.09249506636175525)
  a61=       T(5.86145544294642)
  a62=       T(-12.92096931784711)
  a63=       T(8.159367898576159)
  a64=       T(-0.071584973281401)
  a65=       T(-0.028269050394068383)
  a71=       T(0.09646076681806523)
  a72=       T(0.01)
  a73=       T(0.4798896504144996)
  a74=       T(1.379008574103742)
  a75=       T(-3.290069515436081)
  a76=       T(2.324710524099774)
  # b1 =       T(0.09468075576583945)
  # b2 =       T(0.009183565540343254)
  # b3 =       T(0.4877705284247616)
  # b4 =       T(1.234297566930479)
  # b5 =       T(-2.7077123499835256)
  # b6 =       T(1.866628418170587)
  # b7 =       T(0.015151515151515152)
  btilde1 =  T(-0.00178001105222577714)
  btilde2 =  T(-0.0008164344596567469)
  btilde3 =  T(0.007880878010261995)
  btilde4 =  T(-0.1447110071732629)
  btilde5 =  T(0.5823571654525552)
  btilde6 =  T(-0.45808210592918697)
  btilde7 =  T(0.015151515151515152)

  r11,r12,r13,r14,r22,r23,r24,r32,r33,r34,r42,r43,r44,r52,r53,r54,r62,r63,r64,r72,r73,r74 = Tsit5Interp(T)
  Tsit5ConstantCache(c1,c2,c3,c4,c5,c6,a21,a31,a32,a41,a42,a43,a51,a52,a53,a54,a61,a62,a63,a64,a65,a71,a72,a73,a74,a75,a76,btilde1,btilde2,btilde3,btilde4,btilde5,btilde6,btilde7,r11,r12,r13,r14,r22,r23,r24,r32,r33,r34,r42,r43,r44,r52,r53,r54,r62,r63,r64,r72,r73,r74)
end

function Tsit5ConstantCache(T::Type,T2::Type)
  c1 =       T2(161//1000)
  c2 =       T2(327//1000)
  c3 =       T2(9//10)
  c4 =       T2(parse(BigFloat,".9800255409045096857298102862870245954942137979563024768854764293221195950761080302604"))
  c5 =       T2(1)
  c6 =       T2(1)
  a21=       T(161//1000)
  a31=       T(parse(BigFloat,"-.8480655492356988544426874250230774675121177393430391537369234245294192976164141156943e-2"))
  a32=       T(parse(BigFloat,".3354806554923569885444268742502307746751211773934303915373692342452941929761641411569"))
  a41=       T(parse(BigFloat,"2.897153057105493432130432594192938764924887287701866490314866693455023795137503079289"))
  a42=       T(parse(BigFloat,"-6.359448489975074843148159912383825625952700647415626703305928850207288721235210244366"))
  a43=       T(parse(BigFloat,"4.362295432869581411017727318190886861027813359713760212991062156752264926097707165077"))
  a51=       T(parse(BigFloat,"5.325864828439256604428877920840511317836476253097040101202360397727981648835607691791"))
  a52=       T(parse(BigFloat,"-11.74888356406282787774717033978577296188744178259862899288666928009020615663593781589"))
  a53=       T(parse(BigFloat,"7.495539342889836208304604784564358155658679161518186721010132816213648793440552049753"))
  a54=       T(parse(BigFloat,"-.9249506636175524925650207933207191611349983406029535244034750452930469056411389539635e-1"))
  a61=       T(parse(BigFloat,"5.861455442946420028659251486982647890394337666164814434818157239052507339770711679748"))
  a62=       T(parse(BigFloat,"-12.92096931784710929170611868178335939541780751955743459166312250439928519268343184452"))
  a63=       T(parse(BigFloat,"8.159367898576158643180400794539253485181918321135053305748355423955009222648673734986"))
  a64=       T(parse(BigFloat,"-.7158497328140099722453054252582973869127213147363544882721139659546372402303777878835e-1"))
  a65=       T(parse(BigFloat,"-.2826905039406838290900305721271224146717633626879770007617876201276764571291579142206e-1"))
  a71=       T(parse(BigFloat,".9646076681806522951816731316512876333711995238157997181903319145764851595234062815396e-1"))
  a72=       T(1//100)
  a73=       T(parse(BigFloat,".4798896504144995747752495322905965199130404621990332488332634944254542060153074523509"))
  a74=       T(parse(BigFloat,"1.379008574103741893192274821856872770756462643091360525934940067397245698027561293331"))
  a75=       T(parse(BigFloat,"-3.290069515436080679901047585711363850115683290894936158531296799594813811049925401677"))
  a76=       T(parse(BigFloat,"2.324710524099773982415355918398765796109060233222962411944060046314465391054716027841"))
  # b1 =       T(parse(BigFloat,".9468075576583945807478876255758922856117527357724631226139574065785592789071067303271e-1"))
  # b2 =       T(parse(BigFloat,".9183565540343253096776363936645313759813746240984095238905939532922955247253608687270e-2"))
  # b3 =       T(parse(BigFloat,".4877705284247615707855642599631228241516691959761363774365216240304071651579571959813"))
  # b4 =       T(parse(BigFloat,"1.234297566930478985655109673884237654035539930748192848315425833500484878378061439761"))
  # b5 =       T(parse(BigFloat,"-2.707712349983525454881109975059321670689605166938197378763992255714444407154902012702"))
  # b6 =       T(parse(BigFloat,"1.866628418170587035753719399566211498666255505244122593996591602841258328965767580089"))
  # b7 =       T(1//66)
  btilde1 =  T(parse(BigFloat,"-1.780011052225771443378550607539534775944678804333659557637450799792588061629796e-03"))
  btilde2 =  T(parse(BigFloat,"-8.164344596567469032236360633546862401862537590159047610940604670770447527463931e-04"))
  btilde3 =  T(parse(BigFloat,"7.880878010261996010314727672526304238628733777103128603258129604952959142646516e-03"))
  btilde4 =  T(parse(BigFloat,"-1.44711007173262907537165147972635116720922712343167677619514233896760819649515e-01"))
  btilde5 =  T(parse(BigFloat,"5.823571654525552250199376106520421794260781239567387797673045438803694038950012e-01"))
  btilde6 =  T(parse(BigFloat,"-4.580821059291869466616365188325542974428047279788398179474684434732070620889539e-01"))
  btilde7 =  T(1//66)

  r11,r12,r13,r14,r22,r23,r24,r32,r33,r34,r42,r43,r44,r52,r53,r54,r62,r63,r64,r72,r73,r74 = Tsit5Interp(T)
  Tsit5ConstantCache(c1,c2,c3,c4,c5,c6,a21,a31,a32,a41,a42,a43,a51,a52,a53,a54,a61,a62,a63,a64,a65,a71,a72,a73,a74,a75,a76,btilde1,btilde2,btilde3,btilde4,btilde5,btilde6,btilde7,r11,r12,r13,r14,r22,r23,r24,r32,r33,r34,r42,r43,r44,r52,r53,r54,r62,r63,r64,r72,r73,r74)
end

"""
Coefficients for the polynomial
bᵢΘ = ri1*Θ + ri2*Θ^2 + ri3*Θ^3 + ...

These are the coefficients of the expanded form of the polynomials from

Runge–Kutta pairs of order 5(4) satisfying only the first column
simplifying assumption

Ch. Tsitouras
"""
function Tsit5Interp(::Type{T}) where {T<:CompiledFloats}
  r11 = T(1.0)
  r12 = T(-2.763706197274826)
  r13 = T(2.9132554618219126)
  r14 = T(-1.0530884977290216)

  r22 = T(0.13169999999999998)
  r23 = T(-0.2234)
  r24 = T(0.1017)

  r32 = T(3.9302962368947516)
  r33 = T(-5.941033872131505)
  r34 = T(2.490627285651253)

  r42 = T(-12.411077166933676)
  r43 = T(30.33818863028232)
  r44 = T(-16.548102889244902)

  r52 = T(37.50931341651104)
  r53 = T(-88.1789048947664)
  r54 = T(47.37952196281928)

  r62 = T(-27.896526289197286)
  r63 = T(65.09189467479366)
  r64 = T(-34.87065786149661)

  r72 = T(1.5)
  r73 = T(-4)
  r74 = T(2.5)

  return r11,r12,r13,r14,r22,r23,r24,r32,r33,r34,r42,r43,r44,r52,r53,r54,r62,r63,r64,r72,r73,r74
end

"""
Coefficients for the polynomial
bᵢΘ = ri1*Θ + ri2*Θ^2 + ri3*Θ^3 + ...

These are the coefficients of the expanded form of the polynomials from

Runge–Kutta pairs of order 5(4) satisfying only the first column
simplifying assumption

Ch. Tsitouras
"""
function Tsit5Interp(T::Type)
  r11 = T(parse(BigFloat,"0.999999999999999974283372471559910888475488471328"))
  r12 = T(parse(BigFloat,"-2.763706197274825911336735930481400260916070804192"))
  r13 = T(parse(BigFloat,"2.91325546182191274375068099306808"))
  r14 = T(-1.0530884977290216)

  r22 = T(parse(BigFloat,"0.13169999999999999727"))
  r23 = T(parse(BigFloat,"-0.22339999999999999818"))
  r24 = T(0.1017)

  r32 = T(parse(BigFloat,"3.93029623689475152850687446709813398"))
  r33 = T(parse(BigFloat,"-5.94103387213150473470249202589458001"))
  r34 = T(parse(BigFloat,"2.490627285651252793"))

  r42 = T(parse(BigFloat,"-12.411077166933676983734381540685453484102414134010752"))
  r43 = T(parse(BigFloat,"30.3381886302823215981729903691836576"))
  r44 = T(parse(BigFloat,"-16.54810288924490272"))

  r52 = T(parse(BigFloat,"37.50931341651103919496903965334519631242339792120440212"))
  r53 = T(parse(BigFloat,"-88.1789048947664011014276693541209817"))
  r54 = T(parse(BigFloat,"47.37952196281928122"))

  r62 = T(parse(BigFloat,"-27.896526289197287805948263144598643896"))
  r63 = T(parse(BigFloat,"65.09189467479367152629021928716553658"))
  r64 = T(parse(BigFloat,"-34.87065786149660974"))

  r72 = T(1.5)
  r73 = T(-4)
  r74 = T(2.5)

  return r11,r12,r13,r14,r22,r23,r24,r32,r33,r34,r42,r43,r44,r52,r53,r54,r62,r63,r64,r72,r73,r74
end

struct BS5ConstantCache{T,T2} <: OrdinaryDiffEqConstantCache
  c1::T2
  c2::T2
  c3::T2
  c4::T2
  c5::T2
  a21::T
  a31::T
  a32::T
  a41::T
  a42::T
  a43::T
  a51::T
  a52::T
  a53::T
  a54::T
  a61::T
  a62::T
  a63::T
  a64::T
  a65::T
  a71::T
  a72::T
  a73::T
  a74::T
  a75::T
  a76::T
  a81::T
  a83::T
  a84::T
  a85::T
  a86::T
  a87::T
  bhat1::T
  bhat3::T
  bhat4::T
  bhat5::T
  bhat6::T
  btilde1::T
  btilde3::T
  btilde4::T
  btilde5::T
  btilde6::T
  btilde7::T
  btilde8::T
  c6::T2
  c7::T2
  c8::T2
  a91::T
  a92::T
  a93::T
  a94::T
  a95::T
  a96::T
  a97::T
  a98::T
  a101::T
  a102::T
  a103::T
  a104::T
  a105::T
  a106::T
  a107::T
  a108::T
  a109::T
  a111::T
  a112::T
  a113::T
  a114::T
  a115::T
  a116::T
  a117::T
  a118::T
  a119::T
  a1110::T
  r016::T
  r015::T
  r014::T
  r013::T
  r012::T
  r036::T
  r035::T
  r034::T
  r033::T
  r032::T
  r046::T
  r045::T
  r044::T
  r043::T
  r042::T
  r056::T
  r055::T
  r054::T
  r053::T
  r052::T
  r066::T
  r065::T
  r064::T
  r063::T
  r062::T
  r076::T
  r075::T
  r074::T
  r073::T
  r072::T
  r086::T
  r085::T
  r084::T
  r083::T
  r082::T
  r096::T
  r095::T
  r094::T
  r093::T
  r106::T
  r105::T
  r104::T
  r103::T
  r102::T
  r116::T
  r115::T
  r114::T
  r113::T
  r112::T
end

"""

An Efficient Runge-Kutta (4,5) Pair by P.Bogacki and L.F.Shampine
 Computers and Mathematics with Applications, Vol. 32, No. 6, 1996, pages 15 to 28
"""
function BS5ConstantCache(::Type{T},::Type{T2}) where {T<:CompiledFloats,T2<:CompiledFloats}
  c1     =T2(0.16666666666666666)
  c2     =T2(0.2222222222222222)
  c3     =T2(0.42857142857142855)
  c4     =T2(0.6666666666666666)
  c5     =T2(0.75)

  a21    =T(0.16666666666666666)
  a31    =T(0.07407407407407407)
  a32    =T(0.14814814814814814)
  a41    =T(0.13338192419825073)
  a42    =T(-0.47230320699708456)
  a43    =T(0.7674927113702624)
  a51    =T(0.22895622895622897)
  a52    =T(-0.36363636363636365)
  a53    =T(0.2937062937062937)
  a54    =T(0.5076405076405076)
  a61    =T(0.026500355113636364)
  a62    =T(0.23011363636363635)
  a63    =T(0.10772747760052448)
  a64    =T(0.1602190777972028)
  a65    =T(0.225439453125)
  a71    =T(0.18159821692916506)
  a72    =T(-0.38707982536247293)
  a73    =T(0.41288268560074054)
  a74    =T(0.6409913191253497)
  a75    =T(-1.0121439383514517)
  a76    =T(1.1637515420586695)
  a81    =T(0.07279265873015874)
  a83    =T(0.28662437773692473)
  a84    =T(0.19513621794871794)
  a85    =T(0.008638392857142857)
  a86    =T(0.3595655806182122)
  a87    =T(0.07724277210884353)

  bhat1  =T(-0.00234375)
  bhat3  =T(0.0103760754048583)
  bhat4  =T(-0.016490384615384615)
  bhat5  =T(0.018984375)
  bhat6  =T(-0.010526315789473684)
#  btilde1=T(0.07084476451760402)
#  btilde2=T(0)
#  btilde3=T(0.2956730769230769)
#  btilde4=T(0.17965747482208388)
#  btilde5=T(0.029861111111111113)
#  btilde6=T(0.3462886755067825)
#  btilde7=T(0.07176240133870539)
  btilde1=T(-0.0019478942125547064)
  btilde3=T(0.009048699186152193)
  btilde4=T(-0.015478743126634073)
  btilde5=T(0.021222718253968254)
  btilde6=T(-0.013276905111429694)
  btilde7=T(-0.005480370770138146)
  btilde8=T(0.005912495780636172)
  c6,c7,c8,a91,a92,a93,a94,a95,a96,a97,a98,a101,a102,a103,a104,a105,a106,a107,a108,a109,a111,a112,a113,a114,a115,a116,a117,a118,a119,a1110 = BS5Interp(T,T2)
  r016,r015,r014,r013,r012,r036,r035,r034,r033,r032,r046,r045,r044,r043,r042,r056,r055,r054,r053,r052,r066,r065,r064,r063,r062,r076,r075,r074,r073,r072,r086,r085,r084,r083,r082,r096,r095,r094,r093,r106,r105,r104,r103,r102,r116,r115,r114,r113,r112 = BS5Interp_polyweights(T)
  BS5ConstantCache(c1,c2,c3,c4,c5,a21,a31,a32,a41,a42,a43,a51,a52,a53,a54,a61,a62,a63,a64,a65,a71,a72,a73,a74,a75,a76,a81,a83,a84,a85,a86,a87,bhat1,bhat3,bhat4,bhat5,bhat6,btilde1,btilde3,btilde4,btilde5,btilde6,btilde7,btilde8,c6,c7,c8,a91,a92,a93,a94,a95,a96,a97,a98,a101,a102,a103,a104,a105,a106,a107,a108,a109,a111,a112,a113,a114,a115,a116,a117,a118,a119,a1110,r016,r015,r014,r013,r012,r036,r035,r034,r033,r032,r046,r045,r044,r043,r042,r056,r055,r054,r053,r052,r066,r065,r064,r063,r062,r076,r075,r074,r073,r072,r086,r085,r084,r083,r082,r096,r095,r094,r093,r106,r105,r104,r103,r102,r116,r115,r114,r113,r112)
end

"""

An Efficient Runge-Kutta (4,5) Pair by P.Bogacki and L.F.Shampine
 Computers and Mathematics with Applications, Vol. 32, No. 6, 1996, pages 15 to 28
"""
function BS5ConstantCache(T::Type,T2::Type)
  c1     =T2(1//6)
  c2     =T2(2//9)
  c3     =T2(3//7)
  c4     =T2(2//3)
  c5     =T2(3//4)

  a21    =T(1//6)
  a31    =T(2//27)
  a32    =T(4//27)
  a41    =T(183//1372)
  a42    =T(-162//343)
  a43    =T(1053//1372)
  a51    =T(68//297)
  a52    =T(-4//11)
  a53    =T(42//143)
  a54    =T(1960//3861)
  a61    =T(597//22528)
  a62    =T(81//352)
  a63    =T(63099//585728)
  a64    =T(58653//366080)
  a65    =T(4617//20480)
  a71    =T(174197//959244)
  a72    =T(-30942//79937)
  a73    =T(8152137//19744439)
  a74    =T(666106//1039181)
  a75    =T(-29421//29068)
  a76    =T(482048//414219)
  a81    =T(587//8064)
  a83    =T(4440339//15491840)
  a84    =T(24353//124800)
  a85    =T(387//44800)
  a86    =T(2152//5985)
  a87    =T(7267//94080)

  bhat1  =T(-3//1280)
  bhat3  =T(6561//632320)
  bhat4  =T(-343//20800)
  bhat5  =T(243//12800)
  bhat6  =T(-1//95)
  #  btilde1=T(2479//34992)
  #  btilde2=T(0)
  #  btilde3=T(123//416)
  #  btilde4=T(612941//3411720)
  #  btilde5=T(43//1440)
  #  btilde6=T(2272//6561)
  #  btilde7=T(79937//1113912)
  btilde1=T(-3817//1959552)
  btilde3=T(140181//15491840)
  btilde4=T(-4224731//272937600)
  btilde5=T(8557//403200)
  btilde6=T(-57928//4363065)
  btilde7=T(-23930231//4366535040)
  btilde8=T(3293//556956)
  c6,c7,c8,a91,a92,a93,a94,a95,a96,a97,a98,a101,a102,a103,a104,a105,a106,a107,a108,a109,a111,a112,a113,a114,a115,a116,a117,a118,a119,a1110 = BS5Interp(T,T2)
  r016,r015,r014,r013,r012,r036,r035,r034,r033,r032,r046,r045,r044,r043,r042,r056,r055,r054,r053,r052,r066,r065,r064,r063,r062,r076,r075,r074,r073,r072,r086,r085,r084,r083,r082,r096,r095,r094,r093,r106,r105,r104,r103,r102,r116,r115,r114,r113,r112 = BS5Interp_polyweights(T)
  BS5ConstantCache(c1,c2,c3,c4,c5,a21,a31,a32,a41,a42,a43,a51,a52,a53,a54,a61,a62,a63,a64,a65,a71,a72,a73,a74,a75,a76,a81,a83,a84,a85,a86,a87,bhat1,bhat3,bhat4,bhat5,bhat6,btilde1,btilde3,btilde4,btilde5,btilde6,btilde7,btilde8,c6,c7,c8,a91,a92,a93,a94,a95,a96,a97,a98,a101,a102,a103,a104,a105,a106,a107,a108,a109,a111,a112,a113,a114,a115,a116,a117,a118,a119,a1110,r016,r015,r014,r013,r012,r036,r035,r034,r033,r032,r046,r045,r044,r043,r042,r056,r055,r054,r053,r052,r066,r065,r064,r063,r062,r076,r075,r074,r073,r072,r086,r085,r084,r083,r082,r096,r095,r094,r093,r106,r105,r104,r103,r102,r116,r115,r114,r113,r112)
end

"""

An Efficient Runge-Kutta (4,5) Pair by P.Bogacki and L.F.Shampine
 Computers and Mathematics with Applications, Vol. 32, No. 6, 1996, pages 15 to 28

Used in the lazy construction of the dense output

k9, k10, k11 are not computed until called in the dense routine
"""
function BS5Interp(::Type{T},::Type{T2}) where {T<:CompiledFloats,T2<:CompiledFloats}

  c6    = T2(0.5)
  c7    = T2(0.8333333333333334)
  c8    = T2(0.1111111111111111)
  a91   = T(0.07405598958333333)
  a92   = T(0)
  a93   = T(0.28964485093442743)
  a94   = T(0.12839214966168092)
  a95   = T(-0.003779296875)
  a96   = T(0.014230019493177388)
  a97   = T(-0.03379371279761905)
  a98   = T(0.03125)
  a101  = T(-0.06358724036162344)
  a102  = T(0.5742461924818869)
  a103  = T(-0.06365063007249953)
  a104  = T(0.043159777438314964)
  a105  = T(0.8370112883898733)
  a106  = T(-0.34045447246719235)
  a107  = T(0.04926503818334922)
  a108  = T(-0.006882677669165967)
  a109  = T(-0.19577394258960973)
  a111  = T(0.0636090772400987)
  a112  = T(0.01057854182854183)
  a113  = T(0.06600100945670531)
  a114  = T(0.02048391555358402)
  a115  = T(0.003682270330219549)
  a116  = T(0.155258632271002)
  a117  = T(-0.08509702513818027)
  a118  = T(0.1)
  a119  = T(-0.1)
  a1110 = T(-0.12340531043086005)

  return c6,c7,c8,a91,a92,a93,a94,a95,a96,a97,a98,a101,a102,a103,a104,a105,a106,a107,a108,a109,a111,a112,a113,a114,a115,a116,a117,a118,a119,a1110
end

"""

An Efficient Runge-Kutta (4,5) Pair by P.Bogacki and L.F.Shampine
 Computers and Mathematics with Applications, Vol. 32, No. 6, 1996, pages 15 to 28

Used in the lazy construction of the dense output

k9, k10, k11 are not computed until called in the dense routine
"""
function BS5Interp(T::Type,T2::Type)

  c6    = T2(1//2)
  c7    = T2(5//6)
  c8    = T2(1//9)
  a91   = T(455//6144)
  a92   = T(0)
  a93   = T(10256301//35409920)
  a94   = T(2307361//17971200)
  a95   = T(-387//102400)
  a96   = T(73//5130)
  a97   = T(-7267//215040)
  a98   = T(1//32)
  a101  = T(-837888343715//13176988637184)
  a102  = T(30409415//52955362)
  a103  = T(-48321525963//759168069632)
  a104  = T(8530738453321//197654829557760)
  a105  = T(1361640523001//1626788720640)
  a106  = T(-13143060689//38604458898)
  a107  = T(18700221969//379584034816)
  a108  = T(-5831595//847285792)
  a109  = T(-5183640//26477681)
  a111  = T(98719073263//1551965184000)
  a112  = T(1307//123552)
  a113  = T(4632066559387//70181753241600)
  a114  = T(7828594302389//382182512025600)
  a115  = T(40763687//11070259200)
  a116  = T(34872732407//224610586200)
  a117  = T(-2561897//30105600)
  a118  = T(1//10)
  a119  = T(-1//10)
  a1110 = T(-1403317093//11371610250)

  return c6,c7,c8,a91,a92,a93,a94,a95,a96,a97,a98,a101,a102,a103,a104,a105,a106,a107,a108,a109,a111,a112,a113,a114,a115,a116,a117,a118,a119,a1110
end

"""
Coefficients for the polynomial
bᵢΘ = ri1*Θ + ri2*Θ^2 + ri3*Θ^3 + ...

These coefficients are taken from RKSuite

Note that RKSuite has an error: r081 should be 0
and r011 should be 1. This is pretty easy to spot
since the first order interpolation is linear from y₀.
"""
function BS5Interp_polyweights(::Type{T}) where T<:CompiledFloats
  r016   = T(-11.547607240534195)
  r015   = T(36.89579791207878)
  r014   = T(-45.470343197183475)
  r013   = T(27.298136302807084)
  r012   = T(-8.103191118438032)

  r036   = T(-27.245925284262768)
  r035   = T(74.85879078710211)
  r034   = T(-72.2759709558427)
  r033   = T(28.38602193195626)
  r032   = T(-3.4362921012159933)

  r046   = T(-14.307769126529058)
  r045   = T(38.240038148817945)
  r044   = T(-35.469854845413806)
  r043   = T(13.060399314592578)
  r042   = T(-1.3276772735189397)

  r056   = T(-0.7296779223862557)
  r055   = T(1.9817123385873385)
  r054   = T(-1.964240395021645)
  r053   = T(0.8847786781120115)
  r052   = T(-0.16393430643430643)

  r066   = T(19.121088881250603)
  r065   = T(-65.9928405785889)
  r064   = T(77.4690381021765)
  r063   = T(-34.163041154825144)
  r062   = T(3.9253203306051474)

  r076   = T(-14.676126700680273)
  r075   = T(42.17455357142857)
  r074   = T(-42.4062818877551)
  r073   = T(16.83892431972789)
  r072   = T(-1.853826530612245)

  r086   = T(19.453125)
  r085   = T(-54.359375)
  r084   = T(53.671875)
  r083   = T(-21.078125)
  r082   = T(2.3125)

  r096   = T(18.333333333333332)
  r095   = T(-39)
  r094   = T(23)
  r093   = T(-2.3333333333333335)

  r106   = T(-23.400440940191388)
  r105   = T(70.20132282057416)
  r104   = T(-73.55422182095978)
  r103   = T(30.106238940962648)
  r102   = T(-3.3528990003856314)

  r116   = T(35)
  r115   = T(-105)
  r114   = T(117)
  r113   = T(-59)
  r112   = T(12)

  return r016,r015,r014,r013,r012,r036,r035,r034,r033,r032,r046,r045,r044,r043,r042,r056,r055,r054,r053,r052,r066,r065,r064,r063,r062,r076,r075,r074,r073,r072,r086,r085,r084,r083,r082,r096,r095,r094,r093,r106,r105,r104,r103,r102,r116,r115,r114,r113,r112
end

"""
Coefficients for the polynomial
bᵢΘ = ri1*Θ + ri2*Θ^2 + ri3*Θ^3 + ...

These coefficients are taken from RKSuite

Note that RKSuite has an error: r081 should be 0
and r011 should be 1. This is pretty easy to spot
since the first order interpolation is linear from y₀.
"""
function BS5Interp_polyweights(T::Type)
  r016   = T(-12134338393//1050809760)
  r015   = T(12923488183//350269920)
  r014   = T(-2722545893//59875200)
  r013   = T(35856435071//1313512200)
  r012   = T(-3547880131//437837400)

  r036   = T(-33197340367//1218433216)
  r035   = T(65150312289//870309440)
  r034   = T(-27096444225//374902528)
  r033   = T(4323308999//152304152)
  r032   = T(-1046723109//304608304)

  r046   = T(-284800997201//19905339168)
  r045   = T(6343174409579//165877826400)
  r044   = T(-201150852119//5671036800)
  r043   = T(3249645975331//248816739600)
  r042   = T(-55058055073//41469456600)

  r056   = T(-540919//741312)
  r055   = T(85695583//43243200)
  r054   = T(-2903933//1478400)
  r053   = T(3586937//4054050)
  r052   = T(-1772261//10810800)

  r066   = T(7157998304//374350977)
  r065   = T(-41174140576//623918295)
  r064   = T(413114104//5332635)
  r063   = T(-9134977024//267393555)
  r062   = T(2449079168//623918295)

  r076   = T(-138073//9408)
  r075   = T(94471//2240)
  r074   = T(-1329861//31360)
  r073   = T(792103//47040)
  r072   = T(-7267//3920)

  r086   = T(1245//64)
  r085   = T(-3479//64)
  r084   = T(3435//64)
  r083   = T(-1349//64)
  r082   = T(37//16)

  r096   = T(55//3)
  r095   = T(-39)
  r094   = T(23)
  r093   = T(-7//3)

  r106   = T(-1774004627//75810735)
  r105   = T(1774004627//25270245)
  r104   = T(-26477681//359975)
  r103   = T(11411880511//379053675)
  r102   = T(-423642896//126351225)

  r116   = T(35)
  r115   = T(-105)
  r114   = T(117)
  r113   = T(-59)
  r112   = T(12)

  return r016,r015,r014,r013,r012,r036,r035,r034,r033,r032,r046,r045,r044,r043,r042,r056,r055,r054,r053,r052,r066,r065,r064,r063,r062,r076,r075,r074,r073,r072,r086,r085,r084,r083,r082,r096,r095,r094,r093,r106,r105,r104,r103,r102,r116,r115,r114,r113,r112
end

struct DP5ConstantCache{T,T2} <: OrdinaryDiffEqConstantCache
  a21::T
  a31::T
  a32::T
  a41::T
  a42::T
  a43::T
  a51::T
  a52::T
  a53::T
  a54::T
  a61::T
  a62::T
  a63::T
  a64::T
  a65::T
  a71::T
  a73::T
  a74::T
  a75::T
  a76::T
  btilde1::T
  btilde3::T
  btilde4::T
  btilde5::T
  btilde6::T
  btilde7::T
  c1::T2
  c2::T2
  c3::T2
  c4::T2
  c5::T2
  c6::T2
  d1::T
  d3::T
  d4::T
  d5::T
  d6::T
  d7::T
end

function DP5ConstantCache(::Type{T},::Type{T2}) where {T<:CompiledFloats,T2<:CompiledFloats}
  a21 = T(0.2)
  a31 = T(0.075)
  a32 = T(0.225)
  a41 = T(0.9777777777777777)
  a42 = T(-3.7333333333333334)
  a43 = T(3.5555555555555554)
  a51 = T(2.9525986892242035)
  a52 = T(-11.595793324188385)
  a53 = T(9.822892851699436)
  a54 = T(-0.2908093278463649)
  a61 = T(2.8462752525252526)
  a62 = T(-10.757575757575758)
  a63 = T(8.906422717743473)
  a64 = T(0.2784090909090909)
  a65 = T(-0.2735313036020583)
  a71 = T(0.09114583333333333)
  a73 = T(0.44923629829290207)
  a74 = T(0.6510416666666666)
  a75 = T(-0.322376179245283)
  a76 = T(0.13095238095238096)
  # b1  = T(0.08991319444444444)
  # b3  = T(0.4534890685834082)
  # b4  = T(0.6140625)
  # b5  = T(-0.2715123820754717)
  # b6  = T(0.08904761904761904)
  # b7  = T(0.025)
  btilde1 = T(-0.0012326388888888888)
  btilde3 = T(0.0042527702905061394)
  btilde4 = T(-0.03697916666666667)
  btilde5 = T(0.05086379716981132)
  btilde6 = T(-0.0419047619047619)
  btilde7 = T(0.025)
  c1 = T2(0.2)
  c2 = T2(0.3)
  c3 = T2(0.8)
  c4 = T2(0.8888888888888888)
  c5 = T2(1)
  c6 = T2(1)
  d1,d3,d4,d5,d6,d7 = DP5_dense_ds(T)
  DP5ConstantCache(a21,a31,a32,a41,a42,a43,a51,a52,a53,a54,a61,a62,a63,a64,a65,a71,a73,a74,a75,a76,btilde1,btilde3,btilde4,btilde5,btilde6,btilde7,c1,c2,c3,c4,c5,c6,d1,d3,d4,d5,d6,d7)
end

function DP5_dense_ds(::Type{T}) where T<:CompiledFloats
  d1  = T(-1.1270175653862835)
  d3  = T(2.675424484351598)
  d4  = T(-5.685526961588504)
  d5  = T(3.5219323679207912)
  d6  = T(-1.7672812570757455)
  d7  = T(2.382468931778144)
  return d1,d3,d4,d5,d6,d7
end

function DP5ConstantCache(T::Type,T2::Type)
  a21 = T(1//5)
  a31 = T(3//40)
  a32 = T(9//40)
  a41 = T(44/45)
  a42 = T(-56//15)
  a43 = T(32//9)
  a51 = T(19372//6561)
  a52 = T(-25360//2187)
  a53 = T(64448//6561)
  a54 = T(-212//729)
  a61 = T(9017//3168)
  a62 = T(-355//33)
  a63 = T(46732//5247)
  a64 = T(49//176)
  a65 = T(-5103//18656)
  a71 = T(35//384)
  a73 = T(500//1113)
  a74 = T(125//192)
  a75 = T(-2187//6784)
  a76 = T(11//84)
  #  b1  = T(5179//57600)
  #  b3  = T(7571//16695)
  #  b4  = T(393//640)
  #  b5  = T(-92097//339200)
  #  b6  = T(187//2100)
  #  b7  = T(1//40)
  btilde1 = T(-71//57600)
  btilde3 = T(71//16695)
  btilde4 = T(-71//1920)
  btilde5 = T(17253//339200)
  btilde6 = T(-22//525)
  btilde7 = T(1//40)
  c1 = T2(1//5)
  c2 = T2(3//10)
  c3 = T2(4//5)
  c4 = T2(8//9)
  c5 = T2(1)
  c6 = T2(1)
  d1,d3,d4,d5,d6,d7 = DP5_dense_ds(T)
  DP5ConstantCache(a21,a31,a32,a41,a42,a43,a51,a52,a53,a54,a61,a62,a63,a64,a65,a71,a73,a74,a75,a76,btilde1,btilde3,btilde4,btilde5,btilde6,btilde7,c1,c2,c3,c4,c5,c6,d1,d3,d4,d5,d6,d7)
end

function DP5_dense_ds(T)
  d1  = T(-12715105075//11282082432)
  d3  = T(87487479700//32700410799)
  d4  = T(-10690763975//1880347072)
  d5  = T(701980252875//199316789632)
  d6  = T(-1453857185//822651844)
  d7  = T(69997945//29380423)
  return d1,d3,d4,d5,d6,d7
end

#=
function DP5_dense_bs(T)
  b1  = T(5179//57600)
  b3  = T(7571//16695)
  b4  = T(393//640)
  b5  = T(-92097//339200)
  b6  = T(187//2100)
  b7  = T(1//40)
  return b1,b3,b4,b5,b6,b7
end
=#

"""
An Optimized Runge-Kutta method for the solution of Orbital Problems
by Z.A. Anastassi and T.E. Simos
Journal of Computational and Applied Mathematics, Volume 175, Issue 1, 1 March 2005, Pages 1 to 9.
"""
struct Anas5ConstantCache{T,T2} <: OrdinaryDiffEqConstantCache
  a21::T
  a31::T
  a32::T
  a41::T
  a42::T
  a43::T
  a51::T
  a52::T
  a53::T
  a54::T
  a61::T
  a62::T
  a63::T
  a64::T
  a65::T
  c2::T2
  c3::T2
  c4::T2
  c5::T2
  c6::T2
  b1::T
  b3::T
  b4::T
  b5::T
  b6::T
end

function Anas5ConstantCache(::Type{T},::Type{T2}) where {T<:CompiledFloats,T2<:CompiledFloats}
    a21 = T(0.1)
    a31 = T(-0.2222222222222222)
    a32 = T(0.5555555555555556)
    a41 = T(3.111111111111111)
    a42 = T(-4.444444444444444)
    a43 = T(2.0)
    a51 = T(-1.409625)
    a52 = T(2.1375)
    a53 = T(-0.2295)
    a54 = T(0.401625)
    a61 = T(-4.0)
    a62 = T(5.0)
    a63 = T(0.0)
    a64 = T(0.0)
    a65 = T(0.0)
    c2 = T2(0.1)
    c3 = T2(0.3333333333333333)
    c4 = T2(0.6666666666666667)
    c5 = T2(0.9)
    c6 = T2(1)
    b1 = T(0.1064814814814815)
    b3 = T(0.4632352941176471)
    b4 = T(0.1607142857142857)
    b5 = T(0.3112356053532524)
    b6 = T(-.04166666666666667)
    Anas5ConstantCache(a21,a31,a32,a41,a42,a43,a51,a52,a53,a54,a61,a62,a63,a64,a65,c2,c3,c4,c5,c6,b1,b3,b4,b5,b6)
end

function Anas5ConstantCache(::Type{T},::Type{T2}) where {T,T2}
    a21 = T(1//10)
    a31 = T(-2//9)
    a32 = T(5//9)
    a41 = T(28//9)
    a42 = T(-40//9)
    a43 = T(2//1)
    a51 = T(-11277//8000)
    a52 = T(171//80)
    a53 = T(-459//2000)
    a54 = T(3213//8000)
    a61 = T(-4//1)
    a62 = T(5//1)
    a63 = T(0//1)
    a64 = T(0//1)
    a65 = T(0//1)
    c2 = T2(1//10)
    c3 = T2(1//3)
    c4 = T2(2//3)
    c5 = T2(9//10)
    c6 = T2(1//1)
    b1 = T(23//216)
    b3 = T(63//136)
    b4 = T(9//56)
    b5 = T(1000//3213)
    b6 = T(-1//24)
    Anas5ConstantCache(a21,a31,a32,a41,a42,a43,a51,a52,a53,a54,a61,a62,a63,a64,a65,c2,c3,c4,c5,c6,b1,b3,b4,b5,b6)
end


text:
# # Internals
#
# AlgebraOfGraphics is based on *contexts*, which can be extended. A context, together with
# a named tuple (used for `x`, `y` axes or attributes in the plot) forms a `Style`.
# Each `Style` can actually encode many traces, encoded as a list of `NamedTuple => Style`
# pairs, accessible using `pairs(st::Style)`.
#
# ## Contexts
#
# The `DataContext` is invoked with `data(df)`, where `df` respects the Tables.jl interface.

using RDatasets: dataset
using AlgebraOfGraphics
mpg = dataset("ggplot2", "mpg")
t = data(mpg)
st = style(:Cyl, color = :Year => categorical)

pairs(t * st)

# In the `dims(i)` context, all variables in `style` are sliced along the `i`-th dimension.
# The resulting arrays of arrays broadcasted to a common shape, and
# each "inner array" corresponds to a separate trace. The syntax `dims` exists to allow setting
# discrete attributes variables that only vary with one of the dimensions.

d = dims(1) * style(rand(5, 3, 2), rand(5, 3), color = dims(2))
pairs(d)

# ## Combining styles and contexts together
#
# All outputs of `style`, `data`, and `dims` inherit can be combined using `+`
# (adding a new layer), or `*` (merge information in existing layer).

using AbstractPlotting, CairoMakie
mpg1 = copy(mpg)
mpg1.Displ = mpg.Displ .* 0.1
ts = data(mpg) * spec(markersize = 5) + data(mpg1) * spec(markersize=10)
sl = ts * style(:Hwy, :Displ, color = :Cyl => categorical)

# The result can then be plotted using the `draw` function:

sl * spec(Scatter) |> draw
AbstractPlotting.save("combine.svg", AbstractPlotting.current_scene()); nothing #hide

# ![](combine.svg)
#
# ## Implementing a new context
#
# To implement a new context, one can overload:
#
# - `AlgebraOfGraphics._merge(c::MyContext, s1::Style, s2::Style)`, to allow applying `style` to `MyContext`,
#
# - `AlgebraOfGraphics._pairs(c::MyContext, s::Style)`, which must iterates `group::NamedTuple => style::Style` pairs.
#
# See example implementation in the [context file](https://github.com/JuliaPlots/AlgebraOfGraphics.jl/blob/master/src/context.jl).


text:
function __init__mnist()
    DEPNAME = "MNIST"
    TRAINIMAGES = "train-images-idx3-ubyte.gz"
    TRAINLABELS = "train-labels-idx1-ubyte.gz"
    TESTIMAGES  = "t10k-images-idx3-ubyte.gz"
    TESTLABELS  = "t10k-labels-idx1-ubyte.gz"

    register(DataDep(
            DEPNAME,
            """
            Dataset: THE MNIST DATABASE of handwritten digits
            Authors: Yann LeCun, Corinna Cortes, Christopher J.C. Burges
            Website: http://yann.lecun.com/exdb/mnist/

            [LeCun et al., 1998a]
                Y. LeCun, L. Bottou, Y. Bengio, and P. Haffner.
                "Gradient-based learning applied to document recognition."
                Proceedings of the IEEE, 86(11):2278-2324, November 1998

            The files are available for download at the offical
            website linked above. Note that using the data
            responsibly and respecting copyright remains your
            responsibility. The authors of MNIST aren't really
            explicit about any terms of use, so please read the
            website to make sure you want to download the
            dataset.
            """,
            "https://ossci-datasets.s3.amazonaws.com/mnist/" .* [TRAINIMAGES, TRAINLABELS, TESTIMAGES, TESTLABELS],
            "0bb1d5775d852fc5bb32c76ca15a7eb4e9a3b1514a2493f7edfcf49b639d7975",
            # post_fetch_method = DataDeps.unpack  # TODO should we unzip instead of reading directly from the zip?
        ))
end

"""
    MNIST(; Tx=Float32, split=:train, dir=nothing)
    MNIST([Tx, split])

The MNIST database of handwritten digits.

- Authors: Yann LeCun, Corinna Cortes, Christopher J.C. Burges
- Website: http://yann.lecun.com/exdb/mnist/

MNIST is a classic image-classification dataset that is often
used in small-scale machine learning experiments. It contains
70,000 images of handwritten digits. Each observation is a 28x28
pixel gray-scale image that depicts a handwritten version of 1 of
the 10 possible digits (0-9).

# Arguments

$ARGUMENTS_SUPERVISED_ARRAY
- `split`: selects the data partition. Can take the values `:train:` or `:test`.

# Fields

$FIELDS_SUPERVISED_ARRAY
- `split`.

# Methods

$METHODS_SUPERVISED_ARRAY
- [`convert2image`](@ref) converts features to `Gray` images.

# Examples

The images are loaded as a multi-dimensional array of eltype `Tx`.
If `Tx <: Integer`, then all values will be within `0` and `255`,
otherwise the values are scaled to be between `0` and `1`.
`MNIST().features` is a 3D array (i.e. a `Array{Tx,3}`), in
WHN format (width, height, num_images). Labels are stored as
a vector of integers in `MNIST().targets`.

```julia-repl
julia> using MLDatasets: MNIST

julia> dataset = MNIST()
MNIST:
  metadata    =>    Dict{String, Any} with 3 entries
  split       =>    :train
  features    =>    28×28×60000 Array{Float32, 3}
  targets     =>    60000-element Vector{Int64}

julia> dataset[1:5].targets
5-element Vector{Int64}:
7
2
1
0
4

julia> X, y = dataset[];

julia> dataset = MNIST(UInt8, :test)
MNIST:
  metadata    =>    Dict{String, Any} with 3 entries
  split       =>    :test
  features    =>    28×28×10000 Array{UInt8, 3}
  targets     =>    10000-element Vector{Int64}
```
"""
struct MNIST <: SupervisedDataset
    metadata::Dict{String, Any}
    split::Symbol
    features::Array{<:Any,3}
    targets::Vector{Int}
end

MNIST(; split=:train, Tx=Float32, dir=nothing) = MNIST(Tx, split; dir)

function MNIST(Tx::Type, split::Symbol=:train; dir=nothing)
    @assert split in [:train, :test]
    if split === :train
        IMAGESPATH = "train-images-idx3-ubyte.gz"
        LABELSPATH = "train-labels-idx1-ubyte.gz"
    else
        IMAGESPATH  = "t10k-images-idx3-ubyte.gz"
        LABELSPATH  = "t10k-labels-idx1-ubyte.gz"
    end

    features_path = datafile("MNIST", IMAGESPATH, dir)
    features = bytes_to_type(Tx, MNISTReader.readimages(features_path))

    targets_path = datafile("MNIST", LABELSPATH, dir)
    targets = Vector{Int}(MNISTReader.readlabels(targets_path))
    # targets = reshape(targets, 1, :)

    metadata = Dict{String,Any}()
    metadata["n_observations"] = size(features)[end]
    metadata["features_path"] = features_path
    metadata["targets_path"] = targets_path

    return MNIST(metadata, split, features, targets)
end

convert2image(::Type{<:MNIST}, x::AbstractArray{<:Integer}) =
    convert2image(MNIST, reinterpret(N0f8, convert(Array{UInt8}, x)))

function convert2image(::Type{<:MNIST}, x::AbstractArray{T,N}) where {T,N}
    @assert N == 2 || N == 3
    x = permutedims(x, (2, 1, 3:N...))
    return ImageCore.colorview(Gray, x)
end

# DEPRECATED INTERFACE, REMOVE IN v0.7 (or 0.6.x)
# This is has an hack to deprecate the old datasets-are-modules interface
# in favor of a datasets-are-types interface.
# Unfortunately overloading getproperty for a a type doesn't play well for
# parameterized types (e.g. typeof(dataset) hangs).
# After we remove this deprecation path, we can parameterize the type like this:
# MNIST{Tx, A <: AbstractArray{Tx, 3}}
function Base.getproperty(::Type{MNIST}, s::Symbol)
    if s === :traintensor
        @warn "MNIST.traintensor() is deprecated, use `MNIST(split=:train).features` instead." maxlog=2
        traintensor(T::Type=N0f8; kws...) = traintensor(T, :; kws...)
        traintensor(i; kws...) = traintensor(N0f8, i; kws...)
        function traintensor(T::Type, i; dir=nothing)
            MNIST(; split=:train, Tx=T, dir)[i][1]
        end
        return traintensor
    elseif s === :testtensor
        @warn "MNIST.testtensor() is deprecated, use `MNIST(split=:test).features` instead."  maxlog=2
        testtensor(T::Type=N0f8; kws...) = testtensor(T, :; kws...)
        testtensor(i; kws...) = testtensor(N0f8, i; kws...)
        function testtensor(T::Type, i; dir=nothing)
            MNIST(; split=:test, Tx=T, dir)[i][1]
        end
        return testtensor
    elseif s === :trainlabels
        @warn "MNIST.trainlabels() is deprecated, use `MNIST(split=:train).targets` instead."  maxlog=2
        trainlabels(; kws...) = trainlabels(:; kws...)
        function trainlabels(i; dir=nothing)
            MNIST(; split=:train, dir)[i][2]
        end
        return trainlabels
    elseif s === :testlabels
        @warn "MNIST.testlabels() is deprecated, use `MNIST(split=:test).targets` instead." maxlog=2
        testlabels(; kws...) = testlabels(:; kws...)
        function testlabels(i; dir=nothing)
            MNIST(; split=:test, dir)[i][2]
        end
        return testlabels
    elseif s === :traindata
        @warn "MNIST.traindata() is deprecated, use `MNIST(split=:train)[]` instead." maxlog=2
        traindata(T::Type=N0f8; kws...) = traindata(T, :; kws...)
        traindata(i; kws...) = traindata(N0f8, i; kws...)
        function traindata(T::Type, i; dir=nothing)
            MNIST(; split=:train, Tx=T, dir)[i]
        end
        return traindata
    elseif s === :testdata
        @warn "MNIST.testdata() is deprecated, use `MNIST(split=:test)[]` instead."  maxlog=2
        testdata(T::Type=N0f8; kws...) = testdata(T, :; kws...)
        testdata(i; kws...) = testdata(N0f8, i; kws...)
        function testdata(T::Type, i; dir=nothing)
            MNIST(; split=:test, Tx=T, dir)[i]
        end
        return testdata
    elseif s === :convert2image
        @warn "MNIST.convert2image(x) is deprecated, use `convert2image(MNIST, x)` instead"
        return x -> convert2image(MNIST, x)
    else
        return getfield(MNIST, s)
    end
end


text:
# This file is a part of TextSearch.jl

export vectorize_corpus, vectorize, bow

"""
    vectorize(voc::Vocabulary, tokenlist::AbstractVector, buff=BOW())
    vectorize(voc::Vocabulary, textconfig::TextConfig, text::AbstractString, buff=BOW())

Creates a bag of words from the given text (a string or a list of strings).
If bow is given then updates the bag with the text.
When `config` is given, the text is parsed according to it.
"""
function vectorize(voc::Vocabulary, tokenlist::AbstractVector, bow=BOW())
    z = zero(UInt32)
    for token in tokenlist
        tokenID = get(voc, token, z)
        if z != tokenID
            bow[tokenID] = get(bow, tokenID, zero(Int32)) + one(Int32)
        end
    end

    bow
end

function vectorize(voc::Vocabulary, textconfig::TextConfig, text::AbstractString, buff)
    tokens = tokenize(textconfig, text, buff)
    v = vectorize(voc, tokens, buff.bow)
    v
end

function vectorize(voc::Vocabulary, textconfig::TextConfig, text::AbstractString)
    buff = take!(CACHES)
    empty!(buff)
    try
        copy(vectorize(voc, textconfig, text, buff))
    finally
        put!(CACHES, buff)
    end
end

"""
    vectorize(voc::Vocabulary, textconfig::TextConfig, messages::AbstractVector)

Computes a bag of words from messages
"""
function vectorize(voc::Vocabulary, textconfig::TextConfig, messages::AbstractVector, buff)
    empty!(buff.bow)
    for text in messages
        empty!(buff.normtext); empty!(buff.tokens); empty!(buff.unigrams);
        tokens = tokenize(textconfig, text)
        vectorize(voc, tokens, buff.bow)
    end

    buff.bow
end

function vectorize(voc::Vocabulary, textconfig::TextConfig, messages::AbstractVector)
    buff = take!(CACHES)
    try
        copy(vectorize(voc, textconfig, messages, buff))
    finally
        put!(CACHES, buff)
    end
end

"""
    vectorize_corpus(textconfig::TextConfig, textconfig::TextConfig, corpus::AbstractVector; minbatch=0)

Computes a list of bag of words from a corpus
"""
function vectorize_corpus(voc::Vocabulary, textconfig::TextConfig, corpus::AbstractVector; minbatch=0)
    n = length(corpus)
    X = Vector{BOW}(undef, n)
    minbatch = getminbatch(minbatch, n)

    #@batch minbatch=minbatch per=thread
    Threads.@threads for i in 1:n
        X[i] = vectorize(voc, textconfig, corpus[i])
    end

    X
end

max_stars_repo_name: sadit/TextModel

max_issues_repo_name: sadit/TextModel

max_forks_repo_name: sadit/TextModel


text:
function valid_pid(id_number)
   is_v = true
   #has to be 9 digit number to be valid
   n = match(r"[0-9]{9}",id_number)

   if (n == nothing) || (length(id_number)!=9)
      is_v = false
   end

   return is_v
end

function valid_byr(year)
   is_v = true

   #has to be 4 digit 1920<number<2002
   n = match(r"[0-9]{4}",year)
   if (n == nothing) || (length(year)!=4)
      return is_v = false
   end

   y = parse(Int,n.match)
   if (y<1920) || (y>2002)
      is_v = false
   end

   return is_v
end

function valid_iyr(year)
   is_v = true

   #has to be 4 digit 2010<number<2020
   n = match(r"[0-9]{4}",year)
   if (n == nothing) || (length(year)!=4)
      return is_v = false
   end

   y = parse(Int,n.match)
   if (y<2010) || (y>2020)
      is_v = false
   end

   return is_v
end

function valid_eyr(year)
   is_v = true

   #has to be 4 digit 2020<number<2030
   n = match(r"[0-9]{4}",year)
   if (n == nothing) || (length(year)!=4)
      return is_v = false
   end

   y = parse(Int,n.match)
   if (y<2020) || (y>2030)
      is_v = false
   end

   return is_v
end

function valid_ecl(color)
   is_v = true

   #has to be exactly one of:
   #amb blu brn gry grn hzl oth
   n = eachmatch(r"(amb|blu|brn|gry|grn|hzl|oth)",color)
   if (size(collect(n))[1] != 1)
      return is_v = false
   end

   return is_v
end

function valid_hcl(color)
   is_v = true

   #6 characters 0-9 or a-f
   n = match(r"[0-9a-f]{6}",color)
   if (n == nothing) || (length(color) != 6)
      return is_v = false
   end

   return is_v
end

function valid_hgt(height)
   is_v = true

   #number followed by cm  or in
   #150<height_cm<193, 59<height_in<76
   n = match(r"\d+(?=cm)",height)
   if (n == nothing)
      n = match(r"\d+(?=in)",height)
      if (n == nothing)
         is_v = false
      else
         h = parse(Int,n.match)
         if (h<59) || (h>76)
            is_v = false
         end
      end
   else
      h = parse(Int,n.match)
      if (h<150) || (h>193)
         is_v = false
      end
   end

   return is_v
end

function valid(string_to_validate)
   is_valid = true

   req_words  = ["byr","iyr","eyr","hgt", "hcl", "ecl", "pid"]

   for word in req_words
       if ! occursin(word, string_to_validate)
          is_valid = false
          return is_valid
       end
   end

   m = match(r"(?<=pid:)\w+",string_to_validate)
   if (m == nothing) || ! (valid_pid(m.match))
      is_valid = false
      return is_valid
   end

   m = match(r"(?<=byr:)\w+",string_to_validate)
   if (m == nothing) || ! (valid_byr(m.match))
      is_valid = false
      return is_valid
   end

   m = match(r"(?<=iyr:)\w+",string_to_validate)
   if (m == nothing) || ! (valid_iyr(m.match))
      is_valid = false
      return is_valid
   end

   m = match(r"(?<=eyr:)\w+",string_to_validate)
   if (m == nothing) || ! (valid_eyr(m.match))
      is_valid = false
      return is_valid
   end

   m = match(r"(?<=ecl:)\w+",string_to_validate)
   if (m == nothing) || ! (valid_ecl(m.match))
      is_valid = false
      return is_valid
   end

   m = match(r"(?<=hcl:#)\w+",string_to_validate)
   if (m == nothing) || ! (valid_hcl(m.match))
      is_valid = false
      return is_valid
   end

   m = match(r"(?<=hgt:)\w+",string_to_validate)
   if (m == nothing) || ! (valid_hgt(m.match))
      is_valid = false
      return is_valid
   end

   return is_valid
end

#could have done the valid one with Regex
#m = eachmatch(r"(?<=(eyr|pid):)\w+",string_to_study)
#collect(m)[1] would give "pid" and collect(m)[1].match the value after :
#size(collect(m)[1]) would return size of iterator

f = open("input_day4","r")

temp_s = ""
valid_p = 0

#note that I modified my inputs to have an empty line at the end
for lines in readlines(f)
   #if not empty line concatenate de string with space
   if lines != ""
      global temp_s = temp_s*" "*lines
   else
      #println(temp_s)
      if valid(temp_s)
         global valid_p += 1
      end
      global temp_s = ""

   end
end

println(valid_p)


text:
Basic_Materials = 1
Conglomerates = 2
Consumer_Goods = 3
Financial = 4
Healthcare = 5
Industrial_Goods = 6
Services = 7
Technology = 8
Utilities = 9


text:
module ForeignCallbacks

import Base: Libc

struct Node{T}
    next::Ptr{Node{T}} # needs to be first field
    data::T
end
Node(data::T) where T = Node{T}(C_NULL, data)

setnext!(node::Ptr{Node{T}}, next::Ptr{Node{T}}) where {T} =
    unsafe_store!(convert(Ptr{Ptr{Node{T}}}, node), next)

function calloc(::Type{T}) where T
    ptr = Libc.malloc(sizeof(T))
    ccall(:memset, Ptr{Cvoid}, (Ptr{Cvoid}, Cint, Csize_t), ptr, 0, sizeof(T))
    return convert(Ptr{T}, ptr)
end

"""
    SingleConsumerStack{T}

A variant of Treiber stack that assumes single-consumer for simple (hopefully)
correct implementation.

Safety notes: `push!` and `unsafe_push!` on `SingleConsumerStack{T}` can be
called from multiple tasks/threads. `unsafe_push! can be called from foreign
threads without the Julia runtime.  Only one Julia task is allowed to call
`popall!` at any point.  This simplifies the implementation and avoids the use
of "counted pointer" (and hence 128 bit CAS) that would typically be required
for general Treiber stack with manual memory management.
"""
mutable struct SingleConsumerStack{T}
    @atomic top::Ptr{Node{T}}
    function SingleConsumerStack{T}() where T
        @assert Base.datatype_pointerfree(Node{T})
        new{T}(C_NULL)
    end
end

function Base.push!(stack::SingleConsumerStack{T}, data::T) where T
    node = Node(data)
    p_node = convert(Ptr{Node{T}}, Libc.malloc(sizeof(Node{T})))
    Base.unsafe_store!(p_node, node)

    top = @atomic :monotonic stack.top
    while true
        setnext!(p_node, top)
        top, ok = @atomicreplace :acquire_release :monotonic stack.top top => p_node
        ok && return nothing
    end
end

# Manual implementation of `push!` since `unsafe_pointer_to_objref(ptr)::T` creates a gcframe
function unsafe_push!(ptr::Ptr{Cvoid}, data::T) where T
    node = Node(data)
    p_node = convert(Ptr{Node{T}}, Libc.malloc(sizeof(Node{T})))
    Base.unsafe_store!(p_node, node)

    p_top = Ptr{Ptr{Node{T}}}(ptr)
    top = Core.Intrinsics.atomic_pointerref(p_top, :monotonic)
    while true
        setnext!(p_node, top)
        top, ok = Core.Intrinsics.atomic_pointerreplace(
            p_top,
            top,
            p_node,
            :acquire_release,
            :monotonic,
        )
        ok && return nothing
    end
end

popall!(stack::SingleConsumerStack{T}) where T = moveto!(T[], stack)

function moveto!(results::AbstractVector{T}, stack::SingleConsumerStack{T}) where T
    p_node = @atomic :monotonic stack.top
    while true
        p_node, ok = @atomicreplace :acquire_release :monotonic stack.top p_node => C_NULL
        ok && break
    end

    # Copy the node data into `results` vector
    while p_node != C_NULL
        node = unsafe_load(p_node)
        Libc.free(p_node)
        push!(results, node.data)
        p_node = node.next
    end

    return results
end

mutable struct ForeignCallback{T}
    queue::SingleConsumerStack{T}
    cond::Base.AsyncCondition
    task::Task

    function ForeignCallback{T}(callback; fifo::Bool = true) where T
        stack = SingleConsumerStack{T}()
        cond = Base.AsyncCondition()
        mayreverse = fifo ? Iterators.reverse : identity
        task = Threads.@spawn begin
            local results = T[]
            while isopen(cond)
                wait(cond)
                moveto!(results, stack)
                for data in mayreverse(results)
                    Base.errormonitor(Threads.@spawn callback(data))
                end
                empty!(results)
            end
        end
        this = new{T}(stack, cond, task)
        finalizer(this) do this
            close(this.cond)
            # TODO: free queue we are leaking at least one node here
        end
    end
end

struct ForeignToken
    handle::Ptr{Cvoid}
    queue::Ptr{Cvoid}
end

"""
    Create a foreign token for a callback.

Note that the passed `ForeignCallback` must be GC preserved while this token is being used by
the foreign code. The user can pass the token around and use it to invoke the callback with the data
by calling `notify!`.
"""
ForeignToken(fc::ForeignCallback) = ForeignToken(fc.cond.handle, Base.pointer_from_objref(fc.queue))

function notify!(token::ForeignToken, data::T) where T
    unsafe_push!(token.queue, data)
    ccall(:uv_async_send, Cvoid, (Ptr{Cvoid},), token.handle)
    return
end

end # module


text:
const PREFIX = joinpath(@__DIR__, "usr")
const DOWNLOADS = joinpath(PREFIX, "downloads")
const DEPS_FILE = joinpath(@__DIR__, "deps.jl")
const SRC_DIR = joinpath(PREFIX, "src")
const LIB_DIR = joinpath(PREFIX, "lib")

# upstream odpi
const ODPI_VERSION_NUMBER = v"4.2.1"
const ODPI_SOURCE_URL = "https://github.com/oracle/odpi/archive/v$(ODPI_VERSION_NUMBER).tar.gz"

# patched odpi: see issue #21
#const ODPI_VERSION_NUMBER = "3.1.2-dev-numeric"
#const ODPI_SOURCE_URL = "https://github.com/felipenoris/odpi/archive/v$(ODPI_VERSION_NUMBER).tar.gz"

const ODPI_SOURCE_LOCAL_FILEPATH = joinpath(DOWNLOADS, "odpi_source.tar.gz")

@static if Sys.islinux()
    const SHARED_LIB = joinpath(PREFIX, "lib", "libdpi.so.$(ODPI_VERSION_NUMBER)")
elseif Sys.isapple()
    const SHARED_LIB = joinpath(PREFIX, "lib", "libdpi.$(ODPI_VERSION_NUMBER).dylib")
else
    error("Target system not supported.")
end

function mkdir_if_not_exists(dir; verbose::Bool=false)
    if !isdir(dir)
        verbose && println("Creating $dir.")
        mkdir(dir)
    end
end

function rmdir_if_exists(dir; verbose::Bool=false)
    if isdir(dir)
        verbose && println("Removing directory $dir.")
        rm(dir, recursive=true)
    end
end

function download_source_files(; verbose::Bool=false)

    function untar_source_files(; verbose::Bool=false)
        # tar -xf $ODPI_SOURCE_LOCAL_FILEPATH -C $SRC_DIR --strip-components=1
        mkdir_if_not_exists(SRC_DIR, verbose=verbose)
        cmd_array = ["tar", "-xf", ODPI_SOURCE_LOCAL_FILEPATH, "-C", SRC_DIR, "--strip-components=1"]
        actual_cmd = Cmd(cmd_array)
        verbose && println(actual_cmd)
        run(Cmd(actual_cmd, dir=PREFIX))
    end

    mkdir_if_not_exists(PREFIX, verbose=verbose)
    mkdir_if_not_exists(DOWNLOADS, verbose=verbose)

    if !isfile(ODPI_SOURCE_LOCAL_FILEPATH)
        verbose && println("Downloading $ODPI_SOURCE_URL...")
        download(ODPI_SOURCE_URL, ODPI_SOURCE_LOCAL_FILEPATH)
    end

    untar_source_files(verbose=verbose)
end

function clean_src_files(;verbose::Bool=false)
    rmdir_if_exists(SRC_DIR, verbose=verbose)
    rmdir_if_exists(DOWNLOADS, verbose=verbose)
end

function clean_all(;verbose::Bool=false)
    rmdir_if_exists(PREFIX, verbose=verbose)
    if isfile(DEPS_FILE)
        verbose && println("Removing file $DEPS_FILE.")
        rm(DEPS_FILE)
    end
end

function patch(original_file, patch_file, output_file; verbose::Bool=false)
    verbose && println("Applying patch $patch_file to $original_file.")
    @assert isfile(original_file) && isfile(patch_file)

    function copy_content(io_in, io_out)
        while !eof(io_in)
            write(io_out, read(io_in, UInt8))
        end
    end

    open(output_file, "w") do io_out
        open(original_file, "r") do io_in
            copy_content(io_in, io_out)
        end

        open(patch_file, "r") do io_in
            copy_content(io_in, io_out)
        end
    end

    verbose && println("$output_file was generated.")
end

function build_shared_library(; verbose::Bool=false)
    mkdir_if_not_exists(LIB_DIR, verbose=verbose)

    # apply patch
    original_file = joinpath(SRC_DIR, "embed", "dpi.c")
    patch_file = joinpath(@__DIR__, "dpi_patch.c")
    patched_file = joinpath(SRC_DIR, "embed", "dpi_patched.c")
    patch(original_file, patch_file, patched_file, verbose=verbose)

    if Sys.islinux()
        #=
        cc -c -fPIC -I ../include -ldl -o dpi.o dpi_patched.c
        cc -shared -fPIC -Wl,-soname,libdpi.so.3 -o ../../lib/libdpi.so.3.1.0 dpi.o -lc
        =#

        build_script = [
            ["cc", "-c", "-fPIC", "-I", joinpath(SRC_DIR, "include"), "-ldl", "-o", joinpath(SRC_DIR, "embed", "dpi.o"), patched_file],
            ["cc", "-shared", "-fPIC", "-Wl,-soname,libdpi.so.3", "-o", SHARED_LIB, joinpath(SRC_DIR, "embed", "dpi.o"), "-lc"]
        ]

    elseif Sys.isapple()

        # cc -dynamiclib -I ../include -o ../../lib/libdpi.3.1.0.dylib dpi_patched.c

        build_script = [
            ["cc", "-dynamiclib", "-I", joinpath(SRC_DIR, "include"), "-o", SHARED_LIB, patched_file]
        ]

    else
        error("Target system not supported.")
    end

    for cmd_array in build_script
        actual_cmd = Cmd(cmd_array)
        verbose && println(actual_cmd)
        run(Cmd(actual_cmd, dir=PREFIX))
    end

    @assert isfile(SHARED_LIB) "Failed building libdpi shared library."
end

function write_deps_file(;verbose::Bool=false)
    @assert isfile(SHARED_LIB) "Couldn't find shared library $SHARED_LIB."
    lib_file = basename(SHARED_LIB)

    deps_file_content = """
# This file is generated by build.jl and should be called
# from `check_deps()` from within your module's `__init__()` method
import Libdl

const libdpi = joinpath(@__DIR__, "usr", "lib", "$lib_file")
function check_deps()
    global libdpi
    if !isfile(libdpi)
        error("\$(libdpi) does not exist, Please re-run Pkg.build(\\"Oracle\\"), and restart Julia.")
    end

    if Libdl.dlopen_e(libdpi) in (C_NULL, nothing)
        error("\$(libdpi) cannot be opened, Please re-run Pkg.build(\\"Oracle\\"), and restart Julia.")
    end
end
"""

    open(DEPS_FILE, "w") do f
        write(f, deps_file_content)
    end
    verbose && println("Created deps file $DEPS_FILE.")
end

function main(;verbose::Bool=false)
    clean_all(verbose=verbose)
    download_source_files(verbose=verbose)
    build_shared_library(verbose=verbose)
    write_deps_file(verbose=verbose)
    clean_src_files(verbose=verbose)
end

main()


text:
module Torch

using ..Ahorn, Maple

const placements = Ahorn.PlacementDict(
    "Torch" => Ahorn.EntityPlacement(
        Maple.Torch,
        "point",
        Dict{String, Any}(
            "startLit" => false
        )
    ),
    "Torch (Lit)" => Ahorn.EntityPlacement(
        Maple.Torch,
        "point",
        Dict{String, Any}(
            "startLit" => true
        )
    ),
)

function torchSprite(entity::Maple.Torch)
    lit = get(entity.data, "startLit", false)

    return lit ? "objects/temple/litTorch03.png" : "objects/temple/torch00.png"
end

function Ahorn.selection(entity::Maple.Torch)
    x, y = Ahorn.position(entity)
    sprite = torchSprite(entity)

    return Ahorn.getSpriteRectangle(sprite, x, y)
end

function Ahorn.render(ctx::Ahorn.Cairo.CairoContext, entity::Maple.Torch, room::Maple.Room)
    sprite = torchSprite(entity)
    Ahorn.drawSprite(ctx, sprite, 0, 0)
end

end


text:
###############################################
#                                             #
# Main file for starting solving functions    #
# and generating figures                      #
#                                             #
###############################################

# Load packages
using NLopt
using FiniteDifferences
using GridInterpolations
using Plots
using JLD2

# Include main functions
include("nopolicy.jl")
include("withpolicy.jl")
include("figures_paper.jl")

# Include helper functions
include("helperfunctions.jl")

# Solve model with standard parameters
@time noint    = nopolicy();                # no intervention at all
@time notest   = withpolicy("notest");         # no testing
@time untarget = withpolicy("untargettest");   # untargeted testing
@time target   = withpolicy("targettest");     # targeted testing
@time isolate  = withpolicy("isolate");        # isolation

# Save important objects
#save("src/data/calib1_nointervention.jld2", noint);
#save("src/data/calib1_notest_small.jld2", notest);
#save("src/data/calib1_untarget_small.jld2", untarget);
#save("src/data/calib1_target_small.jld2", target);
#save("src/data/calib1_isolate_small.jld2", isolate);

# Load objects
#noint = load("src/data/calib1_nointervention.jld2");
#notest = load("src/data/calib1_notest_small.jld2");
#untarget = load("src/data/calib1_untarget_small.jld2");
#target = load("src/data/calib1_target_small.jld2");
#isolate = load("src/data/calib1_isolate_small.jld2");

# Figures (using freshly-generated data)
createFig5(noint, notest, untarget)
#savefig("src/figures/Fig5_dim20")
createFig6(notest, untarget, target)
#savefig("src/figures/Fig6_dim20")
createFig7(untarget, target, isolate)
#savefig("src/figures/Fig7_dim20")

# Figures (using stored data)
#createFig5()
#createFig6()
#createFig7()


text:
# This file is a part of AstroLib.jl. License is MIT "Expat".
# Copyright (C) 2016 Mosè Giordano.

"""
    juldate(date::DateTime) -> reduced_julia_days

### Purpose ###

Convert from calendar to Reduced Julian Days.

### Explanation ###

Julian Day Number is a count of days elapsed since Greenwich mean noon on 1
January 4713 B.C.  Julian Days are the number of Julian days followed by the
fraction of the day elapsed since the preceding noon.

This function takes the given `date` and returns the number of Julian calendar
days since epoch `1858-11-16T12:00:00` (Reduced Julian Days = Julian Days -
2400000).

### Argument ###

* `date`: date in Julian Calendar, UTC standard.  Each element can be given in `DateTime`
  type or anything that can be converted to that type.

### Output ###

The number of Reduced Julian Days is returned.

### Example ###

Get number of Reduced Julian Days at 2016-03-20T15:24:00.

```jldoctest
julia> using AstroLib, Dates

julia> juldate(DateTime(2016, 03, 20, 15, 24))
57468.14166666667

julia> juldate(2016, 03, 20, 15, 24)
57468.14166666667

julia> juldate("2016-03-20T15:24")
57468.14166666667
```

### Notes ###

Julian Calendar is assumed, thus before `1582-10-15T00:00:00` this function is
*not* the inverse of `daycnv`.  For the conversion proleptic Gregorian date to
number of Julian days, use `jdcnv`, which is the inverse of `daycnv`.

Code of this function is based on IDL Astronomy User's Library.
"""
function juldate(dt::DateTime)
    # Before 1582-10-15 Dates.datetime2julian uses proleptic Gregorian Calendar, instead
    # AstroLib's juldate uses Julian Calendar.  In addition, after that day using
    # Dates.datetime2julian(date)-2.4e6 would give results not precise enough.  For all
    # these reasons, we use here the same algorithm as AstroLib's juldate.
    year, month, day = Dates.yearmonthday(dt)
    hours, minutes, seconds, milliseconds =
        Dates.hour(dt), Dates.minute(dt),
        Dates.second(dt), Dates.millisecond(dt)

    if year < 0
        year += 1
    elseif year == 0
        error("There is no year zero in Julian Calendar")
    end
    # Do not include leap year in January and February.
    if month < 3
        year -= 1
        month += 12
    end
    jd = fld(year, 4) + 365*(year - 1860) + floor(306001 * (month + 1) / 10000) +
        (day + hours/24 + minutes/1440 + seconds/86400 + milliseconds/86400000) - 1055//10

    # Adjust for Gregorian Calendar, started on 1582-10-15 (= RJD -100830.5).
    if jd > -100830.5
        a = div(year, 100)
        jd += 2 - a + fld(a, 4)
    end
    return jd
end

juldate(dt...) = juldate(DateTime(dt...))


text:
import Pkg
Pkg.add("DifferentialEquations")
using DifferentialEquations, Plots
function indep_growth(u, i)
    #define r
    if i == 1
        r = 1.5
        k = 5000
    else
        r = 1.8
        k = 2000
    end
    # return product of r and abundance of focal u
    return (r*u[i]*(1-(u[i]/k)))
end

function dep_growth(u, i, a_pass)
    # determine if a12 or a21 should be used
    if i == 1 # use a12
        #a = 0.01
        a = a_pass
        k = 5000
        return u[i]*(a*u[2])*(1-(u[i]/k))
    else
        a = 0.01
        k = 2000
        return u[i]*(a*u[1])*(1-(u[i]/k))
    end
end

function dep_death(u, i, a_pass)
    if i == 1
        death = dep_growth(u, 2, a_pass)
    else
        death = dep_growth(u, 1, a_pass)
    end
    return death
end

function indep_death(u, i)
    if i == 1
        m = 1.4
    else
        m = 1.3
    end

    return m*u[i]
end

function eqns(du, u, p, t)
    a_pass = p
    for i in 1:1
        du[i] = indep_growth(u, i) +
        dep_growth(u, i, a_pass) -
        dep_death(u, i, a_pass) -
        indep_death(u, i)
    end
    for i in 2:2
        du[i] = indep_growth(u, i) +
        dep_growth(u, i, a_pass) -
        dep_death(u, i, a_pass) -
        indep_death(u, i)
    end
end

function extinction_threshold(out,u,t,integrator)
    for i in 1:2
        out[i] = 10 - u[i]
    end
end

function extinction_affect!(integrator, event_idx)
    # loop again through all species
    for i in 1:2
        if event_idx == i
            integrator.u[i] = 0
        end
    end
end

# function eqns_test!(du, u, p, t)
#     for i in 1:2
#         du[i] = indep_growth(u, i) +
#         dep_growth(u, i) -
#         dep_death(u, i) -
#         dep_death(u, i)
#     end
# end

# set constant stop values
const stop1 = [5.0, 60.0]
const stop2 = [8.0, 75.0]

function inc_mort_cond(u,t,integrator)
  t in stop1
end

# decrease mortality (stop pulse)
function dec_mort_cond(u,t,integrator)
  t in stop2
end

# apply an affect when condition is reached
function affect_inc!(integrator)
  integrator.p += 0.001
end

function affect_dec!(integrator)
  integrator.p -= 0.001
end
# set vector callback
extinction_callback =
    VectorContinuousCallback(extinction_threshold,
                             extinction_affect!,
                             2,
                             save_positions = ntuple(x->true,2))

save_positions = (true, true)
cb_inc = DiscreteCallback(inc_mort_cond,
                            affect_inc!,
                            save_positions=save_positions)
save_positions = (true, true)
cb_dec = DiscreteCallback(dec_mort_cond,
                            affect_inc!,
                            save_positions=save_positions)
steady_callback = TerminateSteadyState(1e-8,
                            1e-6,
                            DiffEqCallbacks.allDerivPass)

callback_set = CallbackSet(cb_inc, cb_dec, extinction_callback, steady_callback)

const tstop = [5.;8.;60.;75.;]
save_values = [range(5.0, 20.0, step = 0.1);
                range(20.0, 59.9, step = 5);
                range(60, 100, step = 1)]
save_values = [range(0,100, step = 10)]

network_id = parse(Int64, ARGS[1])
p = 0.01
u0 = [100.0; 100.0]
tspan = (0.0, 1000.0)
prob = ODEProblem(eqns,
                    u0,
                    tspan,
                    p,
                    tstops=tstop)
@time sol = solve(prob,
                Tsit5(),
                reltol = 1e-08,
                callback = callback_set)
object = hcat(sol.t, sol.u)
writedlm("./out/out$network_id.csv", species_list, ",")
plot(sol)

max_stars_repo_name: colebrookson/trait-based-rewiring

max_issues_repo_name: colebrookson/trait-based-rewiring

max_forks_repo_name: colebrookson/trait-based-rewiring


text:
using UniformScalingMatrixFunctions
using Test

@testset "UniformScalingMatrixFunctions.jl" begin
    @test exp(0I) === 1.0I
    @test exp(-im*π/2*I) ≈ -im*I
    @test exp(-im*π/4*I) ≈ (1-im)/√(2)*I

    @test sin(0I) === 0.0I # This assumes iszero(sin(0)), always true?
    @test sin(π/2*I) ≈ I
    @test cos(0I) === 1.0I # This assumes isone(cos(0)), always true?
    @test cos(π/2*I) ≈ 0I atol=1e-14
    @test tan(π/4*I) ≈ I

    @test asin(I) ≈ π/2*I
    @test acos(I) ≈ 0*I atol=1e-14
    @test atan(I) ≈ π/4*I

    # Test hyperbolics + inverses

    @test log(I) === 0.0I
    @test log(exp(I)) === 1.0I

    @test √(4I) === 2.0I
    @test √(4I) == 2I
    @test √(4.0I) === 2.0I
end


text:
using EPOCHInput
using Test
using Unitful

@testset "EPOCHInput.jl" begin
    folders = [
        "gauss",
        "Nuter2018",
        "Wang-Arefiev2020"
    ]
    results = Dict(
        "gauss" => Dict(
            :nx => 10,
            :a0 => 70,
            :laser => Dict(
                :lambda => 800u"nm",
            )
        ),
        "Nuter2018" => Dict(
            :nx => 640,
            :a0 => 1,
            :laser => Dict(
                :lambda => 800u"nm",
            )
        ),
        "Wang-Arefiev2020" => Dict(
            :a0 => 85,
            :laser => Dict(
                :lambda => 1u"μm",
            )
        )
    )
    @testset "$folder" for folder in folders
        p = parse_input("$folder/input.deck")
        @testset "$k" for (k,v) in results[folder]
            if v isa Dict
                for (_k, _v) in pairs(v)
                    if unit(_v) == NoUnits
                        @test getindex(p[k], _k) == _v
                    else
                        @test getindex(p[k], _k) ≈ _v
                    end
                end
            else
                @test p[k] == v
            end
        end
    end
end