text:
const CVec = Ptr{Cvoid}

struct PetscVec
    ptr::Ref{CVec}

    PetscVec() = new(Ref{Ptr{Cvoid}}())
end

# allows us to pass PetscVec objects directly into CVec ccall signatures
Base.cconvert(::Type{CVec}, vec::PetscVec) = vec.ptr[]

"""
    VecCreate(comm::MPI.Comm, vec::PetscVec)

Wrapper to VecCreate
"""
function VecCreate(comm::MPI.Comm, vec::PetscVec)
    error = ccall((:VecCreate, libpetsc), PetscErrorCode, (MPI.MPI_Comm, Ptr{CVec}), comm, vec.ptr)
    @assert iszero(error)
end

function VecCreate(comm::MPI.Comm = MPI.COMM_WORLD)
    vec = PetscVec()
    VecCreate(comm, vec)
    return vec
end

"""
    VecSetValue(vec::PetscVec, i::PetscInt, v::PetscScalar, mode::InsertMode = INSERT_VALUES)

Wrapper to `VecSetValue`. Indexing starts at 0 (as in PETSc).

# Implementation
For an unknow reason, calling PETSc.VecSetValue leads to an "undefined symbol: VecSetValue" error.
So this wrapper directly call VecSetValues (anyway, this is what is done in PETSc...)
"""
function VecSetValue(vec::PetscVec, i::PetscInt, v::PetscScalar, mode::InsertMode = INSERT_VALUES)
    VecSetValues(vec, PetscIntOne, [i], [v], mode)
end

function VecSetValue(vec::PetscVec, i, v, mode::InsertMode = INSERT_VALUES)
    VecSetValue(vec, PetscInt(i), PetscScalar(v), mode)
end

"""
    VecSetValues(vec::PetscVec, nI::PetscInt, I::Vector{PetscInt}, V::Array{PetscScalar}, mode::InsertMode = INSERT_VALUES)

Wrapper to `VecSetValues`. Indexing starts at 0 (as in PETSc)
"""
function VecSetValues(vec::PetscVec, nI::PetscInt, I::Vector{PetscInt}, V::Array{PetscScalar}, mode::InsertMode = INSERT_VALUES)
    error = ccall((:VecSetValues, libpetsc), PetscErrorCode,
        (CVec, PetscInt, Ptr{PetscInt}, Ptr{PetscScalar}, InsertMode),
        vec, nI, I, V, mode)
    @assert iszero(error)
end

function VecSetValues(vec::PetscVec, I::Vector{PetscInt}, V::Array{PetscScalar}, mode::InsertMode = INSERT_VALUES)
    VecSetValues(vec, PetscInt(length(I)), I, V, mode)
end

function VecSetValues(vec::PetscVec, I, V, mode::InsertMode = INSERT_VALUES)
    VecSetValues(vec, PetscInt.(I), PetscScalar.(V), mode)
end

"""
    VecSetSizes(vec::PetscVec, nrows_loc, nrows_glo)

Wrapper to VecSetSizes
"""
function VecSetSizes(vec::PetscVec, nrows_loc, nrows_glo)
    nr_loc = PetscInt(nrows_loc)
    nr_glo = PetscInt(nrows_glo)
    error = ccall((:VecSetSizes, libpetsc),
                PetscErrorCode,
                (CVec, PetscInt, PetscInt),
                vec, nr_loc, nr_glo
            )
    @assert iszero(error)
end

"""
    VecSetFromOptions(vec::PetscVec)

Wrapper to VecSetFromOptions
"""
function VecSetFromOptions(vec::PetscVec)
    error = ccall((:VecSetFromOptions, libpetsc), PetscErrorCode, (CVec,), vec)
    @assert iszero(error)
end

"""
    VecSetUp(vec::PetscVec)

Wrapper to VecSetUp
"""
function VecSetUp(vec::PetscVec)
    error = ccall((:VecSetUp, libpetsc), PetscErrorCode, (CVec,), vec)
    @assert iszero(error)
end

"""
    VecGetOwnershipRange(vec::PetscVec)

Wrapper to `VecGetOwnershipRange`

The result `(rstart, rend)` is a Tuple indicating the rows handled by the local processor.

# Warning
`PETSc` indexing starts at zero (so `rstart` may be zero) and `rend-1` is the last row
handled by the local processor.
"""
function VecGetOwnershipRange(vec::PetscVec)
    rstart = Ref{PetscInt}(0)
    rend = Ref{PetscInt}(0)

    error = ccall((:VecGetOwnershipRange, libpetsc), PetscErrorCode, (CVec, Ref{PetscInt}, Ref{PetscInt}), vec, rstart, rend)
    @assert iszero(error)

    return rstart[], rend[]
end

"""
    VecGetSize(vec::PetscVec)

Wrapper for VecGetSize
"""
function VecGetSize(vec::PetscVec)
    n = Ref{PetscInt}()

    error = ccall((:VecGetSize, libpetsc), PetscErrorCode, (CVec, Ref{PetscInt}), vec, n)
    @assert iszero(error)

    return n[]
end

"""
    VecGetLocalSize(vec::PetscVec)

Wrapper for VecGetLocalSize
"""
function VecGetLocalSize(vec::PetscVec)
    n = Ref{PetscInt}()

    error = ccall((:VecGetLocalSize, libpetsc), PetscErrorCode, (CVec, Ref{PetscInt}), vec, n)
    @assert iszero(error)

    return n[]
end

"""
    VecAssemblyBegin(vec::PetscVec)

Wrapper to VecAssemblyBegin
"""
function VecAssemblyBegin(vec::PetscVec)
    error = ccall((:VecAssemblyBegin, libpetsc), PetscErrorCode, (CVec,), vec)
    @assert iszero(error)
end

"""
    VecAssemblyEnd(vec::PetscVec)

Wrapper to VecAssemblyEnd
"""
function VecAssemblyEnd(vec::PetscVec)
    error = ccall((:VecAssemblyEnd, libpetsc), PetscErrorCode, (CVec,), vec)
    @assert iszero(error)
end

"""
    VecDuplicate(vec::PetscVec)

Wrapper for VecDuplicate, except that it returns the new vector instead of taking it as an input.
"""
function VecDuplicate(vec::PetscVec)
    x = PetscVec()
    error = ccall((:VecDuplicate, libpetsc), PetscErrorCode, (CVec, Ptr{CVec}), vec, x.ptr)
    @assert iszero(error)

    return x
end

"""
    VecGetArray(vec::PetscVec, own = false)

Wrapper for VecGetArray.

# Warning
I am not confortable at all with memory management, both on the C side and on the Julia side. Use
this at you own risk.

According to Julia documentation, `own` optionally specifies whether Julia should take ownership
of the memory, calling free on the pointer when the array is no longer referenced."

"""
function VecGetArray(vec::PetscVec, own = false)
    # Get array pointer
    array_ref = Ref{Ptr{PetscScalar}}()
    error = ccall((:VecGetArray, libpetsc), PetscErrorCode, (CVec, Ref{Ptr{PetscScalar}}), vec, array_ref)
    @assert iszero(error)

    # Get array size
    rstart, rend = VecGetOwnershipRange(vec)
    n = rend - rstart # this is not `rend - rstart + 1` because of VecGetOwnershipRange convention

    array = unsafe_wrap(Array, array_ref[], n; own = own)

    return array, array_ref
end

"""
    VecRestoreArray(vec::PetscVec, array_ref)

Wrapper for VecRestoreArray. `array_ref` is obtained from `VecGetArray`
"""
function VecRestoreArray(vec::PetscVec, array_ref)
    error = ccall((:VecRestoreArray, libpetsc), PetscErrorCode, (CVec, Ref{Ptr{PetscScalar}}), vec, array_ref)
    @assert iszero(error)
end

"""
    VecView(vec::PetscVec, viewer::PetscViewer)

Wrapper to `VecView`. Currently leads to a fatal error when using with `viewer = PetscViewerStdWorld()`.
Use `VecView(vec)` instead.
"""
function VecView(vec::PetscVec, viewer::PetscViewer)
    error = ccall( (:VecView, libpetsc), PetscErrorCode, (CVec, PetscViewer), vec, viewer);
    @assert iszero(error)
end

"""
    VecView(vec::PetscVec)

Wrapper to `VecView` with the default PETSc viewer. This a fallback for the currently buggy
`VecView(vec::PetscVec, viewer::PetscViewer)`
"""
function VecView(vec::PetscVec)
    error = ccall( (:VecView, libpetsc), PetscErrorCode, (CVec, Ptr{Cvoid}), vec, C_NULL)
    @assert iszero(error)
end

"""
    VecDestroy(vec::PetscVec)

Wrapper to VecDestroy
"""
function VecDestroy(vec::PetscVec)
    error = ccall( (:VecDestroy, libpetsc),
                    PetscErrorCode,
                    (Ptr{CVec},),
                    vec.ptr)
    @assert iszero(error)
end


text:
using Base.Filesystem
using Test
using ModelicaScriptingTools: withOMC, testmodel, loadModel, installAndLoad
using OMJulia: sendExpression

moroot = dirname(@__DIR__)
outdir = joinpath(moroot, "out")
if !ispath(outdir)
    mkdir(outdir)
end
refdir = joinpath(moroot, "regRefData")

# relative tolerance for unit tests
# This is chosen a little larger than the tolerance used for most simulations
# (1e-6) to allow for some error propagation.
# Models that do not have a controlled timing of AP events (namely experiments
# including the N cell model) may still fail tests after small and harmless
# changes, because the error accumulates over time.
rrtol = 1e-5
withOMC(outdir, moroot) do omc
    sendExpression(omc, "setCommandLineOptions(\"-d=newInst,nfAPI\")")
    installAndLoad(omc, "Modelica"; version="3.2.3")
    @testset "Simulate examples" begin
        @testset "SodiumChannelSteady" begin
            testmodel(omc, "InaMo.Examples.ComponentTests.SodiumChannelSteady"; refdir=refdir, regRelTol=rrtol)
        end
        @testset "SodiumChannelIV" begin
            testmodel(omc, "InaMo.Examples.ComponentTests.SodiumChannelIV"; refdir=refdir, regRelTol=rrtol)
        end
        @testset "InwardRectifierLin" begin
            testmodel(omc, "InaMo.Examples.ComponentTests.InwardRectifierLin"; refdir=refdir, regRelTol=rrtol)
        end
        @testset "GHKFlux" begin
            testmodel(omc, "InaMo.Examples.ComponentTests.GHKFlux"; override=Dict(
                "stopTime"=>100, "numberOfIntervals"=>1000
            ), refdir=refdir, regRelTol=rrtol)
        end
        @testset "LTypeCalciumSteady" begin
            testmodel(omc, "InaMo.Examples.ComponentTests.LTypeCalciumSteady"; refdir=refdir, regRelTol=rrtol)
        end
        @testset "LTypeCalciumIV" begin
            testmodel(omc, "InaMo.Examples.ComponentTests.LTypeCalciumIV"; refdir=refdir, regRelTol=rrtol)
        end
        @testset "LTypeCalciumIVN" begin
            testmodel(omc, "InaMo.Examples.ComponentTests.LTypeCalciumIVN"; refdir=refdir, regRelTol=rrtol)
        end
        @testset "LTypeCalciumStep" begin
            testmodel(omc, "InaMo.Examples.ComponentTests.LTypeCalciumStep"; refdir=refdir, regRelTol=rrtol)
        end
        @testset "TransientOutwardSteady" begin
            testmodel(omc, "InaMo.Examples.ComponentTests.TransientOutwardSteady"; refdir=refdir, regRelTol=rrtol)
        end
        @testset "TransientOutwardIV" begin
            testmodel(omc, "InaMo.Examples.ComponentTests.TransientOutwardIV"; refdir=refdir, regRelTol=rrtol)
        end
        @testset "RapidDelayedRectifierSteady" begin
            testmodel(omc, "InaMo.Examples.ComponentTests.RapidDelayedRectifierSteady"; refdir=refdir, regRelTol=rrtol)
        end
        @testset "RapidDelayedRectifierIV" begin
            testmodel(omc, "InaMo.Examples.ComponentTests.RapidDelayedRectifierIV"; refdir=refdir, regRelTol=rrtol)
        end
        @testset "HyperpolarizationActivatedSteady" begin
            testmodel(omc, "InaMo.Examples.ComponentTests.HyperpolarizationActivatedSteady"; refdir=refdir, regRelTol=rrtol)
        end
        @testset "HyperpolarizationActivatedIV" begin
            testmodel(omc, "InaMo.Examples.ComponentTests.HyperpolarizationActivatedIV"; refdir=refdir, regRelTol=rrtol)
        end
        @testset "SustainedInwardSteady" begin
            testmodel(omc, "InaMo.Examples.ComponentTests.SustainedInwardSteady"; refdir=refdir, regRelTol=rrtol)
        end
        @testset "SustainedInwardIV" begin
            testmodel(omc, "InaMo.Examples.ComponentTests.SustainedInwardIV"; refdir=refdir, regRelTol=rrtol)
        end
        @testset "SustainedInwardIVKurata" begin
            testmodel(omc, "InaMo.Examples.ComponentTests.SustainedInwardIVKurata"; refdir=refdir, regRelTol=rrtol)
        end
        @testset "SodiumPotassiumPumpLin" begin
            testmodel(omc, "InaMo.Examples.ComponentTests.SodiumPotassiumPumpLin"; refdir=refdir, regRelTol=rrtol)
        end
        @testset "SodiumCalciumExchangerRampInada" begin
            testmodel(omc, "InaMo.Examples.ComponentTests.SodiumCalciumExchangerRampInada"; refdir=refdir, regRelTol=rrtol)
        end
        @testset "SodiumCalciumExchangerLinMatsuoka" begin
            testmodel(omc, "InaMo.Examples.ComponentTests.SodiumCalciumExchangerLinMatsuoka"; refdir=refdir, regRelTol=rrtol)
        end
        @testset "SodiumCalciumExchangerLinKurata" begin
            testmodel(omc, "InaMo.Examples.ComponentTests.SodiumCalciumExchangerLinKurata"; refdir=refdir, regRelTol=rrtol)
        end
        @testset "AllCells" begin
            # we use a different tolerance here, because errors in N cell model accumulate over time
            testmodel(omc, "InaMo.Examples.FullCell.AllCells"; refdir=refdir, regRelTol=1e-3)
        end
        @testset "AllCellsC" begin
            # we use a different tolerance here, because errors in N cell model accumulate over time
            testmodel(omc, "InaMo.Examples.FullCell.AllCellsC"; refdir=refdir, regRelTol=1e-3)
        end
        @testset "FullCellSpon" begin
            testmodel(omc, "InaMo.Examples.FullCell.FullCellSpon"; override=Dict(
                "stopTime"=>0.5, "numberOfIntervals"=>5000, "variableFilter"=>raw"cell\.(naca|cal)\.i",
                "tolerance"=>1e-6
            ), refdir=refdir, regRelTol=rrtol)
        end
        @testset "CaHandlingApprox" begin
            testmodel(omc, "InaMo.Examples.ComponentTests.CaHandlingApprox"; refdir=refdir, regRelTol=rrtol)
        end
        @testset "CaBuffer" begin
            testmodel(omc, "InaMo.Examples.ComponentTests.CaBuffer"; refdir=refdir, regRelTol=rrtol)
        end
        @testset "CaBuffer2" begin
            testmodel(omc, "InaMo.Examples.ComponentTests.CaBuffer2"; refdir=refdir, regRelTol=rrtol)
        end
        @testset "CaDiffusion" begin
            testmodel(omc, "InaMo.Examples.ComponentTests.CaDiffusion"; refdir=refdir, regRelTol=rrtol)
        end
        @testset "CaSERCA" begin
            testmodel(omc, "InaMo.Examples.ComponentTests.CaSERCA"; refdir=refdir, regRelTol=rrtol)
        end
        @testset "CaRyanodineReceptor" begin
            testmodel(omc, "InaMo.Examples.ComponentTests.CaRyanodineReceptor"; refdir=refdir, regRelTol=rrtol)
        end
        @testset "AtrialCellSmokeTest" begin
            # only test if model can be loaded, simulation not working yet
            @test isnothing(loadModel(omc, "InaMo.Examples.FullCell.AllCells"))
            # testmodel(omc, "InaMo.Examples.FullCell.AllCells"; refdir=refdir, regRelTol=rrtol)
        end
        #=
        @testset "SteadyStates" begin
            # only test if model can be simulated
            # TODO this currently does not work, because internal functions cannot be referenced from outside a component
            testmodel(omc, "InaMo.Examples.ComponentTests.SteadyStates"; override=Dict(
                "stopTime"=>0.01
            ), refdir=nothing, regRelTol=rrtol)
        end
        =#
    end
end

ext: jl

max_stars_repo_name: CSchoel/inamo

max_issues_repo_name: CSchoel/inamo

max_forks_repo_name: CSchoel/inamo


text:
#normalized root mean square error
function nrms(res, w)
    @assert length(res) == length(w)
    return sqrt( chi2(res, w)/length(res) )
end

# chi square
function chi2(res, w)
    @assert length(res) == length(w)
    return sum(w .* (res.^2))
end

# B and C matrices from Milani and Gronchi (2010)
function BHC(res, w, npar)
    nobs = length(res)
    B_mat = Matrix{TaylorN{Float64}}(undef, nobs, npar)
    H_mat = Array{TaylorN{Float64}}(undef, nobs, npar, npar)
    C_mat = Array{TaylorN{Float64}}(undef, npar, npar)
    for i in 1:nobs
        B_mat[i,:] .= TaylorSeries.gradient(res[i])
    end
    for i in 1:nobs
        for j in 1:npar
                H_mat[i,j,:] .= TaylorSeries.gradient(B_mat[i,j])
        end
    end
    sqrtw_B = sqrt.(w) .* B_mat
    C_mat .= (sqrtw_B') * sqrtw_B
    return B_mat, H_mat, C_mat
end

# transpose(ξ)*H matrix, with weights
function ξTH(w, res, H_mat, npar)
    ξTHv = Array{Float64}(undef, npar, npar)
    for j in 1:npar
        for i in 1:npar
           ξTHv[i,j] = (w .* res)' * (H_mat[:,i,j])
        end
    end
    return ξTHv
end

# differential corrections subroutine for least-squares fitting
function diffcorr(res, w, x0, niters=5)
    @assert length(res) == length(w)
    nobs = length(res)
    npar = length(x0)
    B_mat, H_mat, C_mat = BHC(res, w, npar)
    D_mat = B_mat' * (w .* res)
    # ξTH_mat = ξTH(w, res, H_mat, npar)
    x_new = x0
    for i in 1:niters
        D = D_mat(x_new)
        C = C_mat(x_new) #.+ ξTH_mat(x_new)
        Δx = - inv(C)*D
        x_new = x_new .+ Δx
        @show sqrt(((Δx')*(C*Δx))/npar)
    end
    C = C_mat(x_new)
    Γ = inv(C)
    return x_new, Γ
end

# Newton method subroutine for least-squares fitting
function newtonls(res, w, x0, niters=5)
    @assert length(res) == length(w)
    nobs = length(res)
    npar = length(x0)
    Q = chi2(res, w)/nobs
    x_new = x0
    for i in 1:niters
        dQ = TaylorSeries.gradient(Q)(x_new)
        d2Q = TaylorSeries.hessian(Q, x_new)
        Δx = - inv(d2Q)*dQ
        x_new = x_new + Δx
        C = d2Q/(2/nobs) # C = d2Q/(2/m)
        @show sqrt(((Δx')*(C*Δx))/npar)
    end
    C = TaylorSeries.hessian(Q, x_new)/(2/nobs) # C = d2Q/(2/m)
    Γ = inv(C)
    return x_new, Γ
end

function newtonls_Q(Q, nobs, x0, niters=5)
    npar = length(x0)
    x_new = x0
    for i in 1:niters
        dQ = TaylorSeries.gradient(Q)(x_new)
        d2Q = TaylorSeries.hessian(Q, x_new)
        Δx = - inv(d2Q)*dQ
        x_new = x_new + Δx
        C = d2Q/(2/nobs) # C = d2Q/(2/m)
        @show sqrt(((Δx')*(C*Δx))/npar)
    end
    C = TaylorSeries.hessian(Q, x_new)/(2/nobs) # C = d2Q/(2/m)
    Γ = inv(C)
    return x_new, Γ
end

# specialized version of newtonls on 6 variables for parametrized orbit determination wrt A2 nongrav coefficient
function newtonls_6v(res, w, x0, niters=5)
    @assert length(res) == length(w)
    nobs = length(res)
    npar = 6 # length(x0)
    Q = chi2(res, w)/nobs
    x_new = x0
    for i in 1:niters
        dQ = TaylorSeries.gradient(Q)(x_new)[1:6]
        d2Q = TaylorSeries.hessian(Q, x_new)[1:6,1:6]
        Δx = - inv(d2Q)*dQ
        x_new[1:6] = x_new[1:6] + Δx
        C = d2Q/(2/nobs) # C = d2Q/(2/m)
        @show sqrt(((Δx')*(C*Δx))/npar)
    end
    C = TaylorSeries.hessian(Q, x_new)/(2/nobs) # C = d2Q/(2/m)
    Γ = inv(C[1:6,1:6])
    return x_new, Γ
end

function newtonls_A2(res, w, x0, niters=5)
    @assert length(res) == length(w)
    nobs = length(res)
    npar = length(x0)
    Q = chi2(res, w)/nobs
    x_new = x0
    for i in 1:niters
        dQ = TaylorSeries.gradient(Q)(x_new)
        d2Q = TaylorSeries.hessian(Q, x_new)
        Δx = - inv(d2Q)*dQ
        x_new[7] = x_new[7] + Δx[7]
        C = d2Q/(2/nobs) # C = d2Q/(2/m)
        @show sqrt(((Δx')*(C*Δx))/npar)
    end
    C = TaylorSeries.hessian(Q, x_new)/(2/nobs) # C = d2Q/(2/m)
    Γ = inv(C)
    return x_new, Γ, C
end


text:
module SimpleSymbolic

struct S
    x::Any
end

Base.show(io::IO, s::S) = print(io, s.x)

function Base.:+(s1::S, s2::S)
    if s1.x == 0
        return s2
    elseif s2.x == 0
        return s1
    else
        if isa(s1.x, Number) && isa(s2.x, Number)
            return S(s1.x + s2.x)
        else
            return S(string(s1.x) * " + " * string(s2.x))
        end
    end
end

function Base.:*(s1::S, s2::S)
    if s1.x == 1
        return s2
    elseif s2.x == 1
        return s1
    elseif s1.x == 0 || s2.x == 0
        return S(0)
    else
        if isa(s1.x, Number) && isa(s2.x, Number)
            return S(s1.x*s2.x)
        else
            if isa(s1.x, String) && isa(s2.x, String)
                if s1.x[1] == '-' && s2.x[1] == '-'
                    return S(string(s1.x[2:end]) * "*" * string(s2.x[2:end]))
                end
            end
            return S(string(s1.x) * "*" * string(s2.x))
        end
    end
end

using StaticArrays

mx1 = @SMatrix [S(1)  S(0)        S(0);
                S(0)  S("cosθ₁")  S("-sinθ₁");
                S(0)  S("sinθ₁")  S("cosθ₁")]

my1 = @SMatrix [S("cosθ₁")   S(0)  S("sinθ₁");
                S(0)         S(1)  S(0);
                S("-sinθ₁")  S(0)  S("cosθ₁")]

mz1 = @SMatrix [S("cosθ₁")  S("-sinθ₁")  S(0);
                S("sinθ₁")  S("cosθ₁")   S(0);
                S(0)        S(0)         S(1)]

mx2 = @SMatrix [S(1)  S(0)        S(0);
                S(0)  S("cosθ₂")  S("-sinθ₂");
                S(0)  S("sinθ₂")  S("cosθ₂")]

my2 = @SMatrix [S("cosθ₂")   S(0)  S("sinθ₂");
                S(0)         S(1)  S(0);
                S("-sinθ₂")  S(0)  S("cosθ₂")]

mz2 = @SMatrix [S("cosθ₂")  S("-sinθ₂")  S(0);
                S("sinθ₂")  S("cosθ₂")   S(0);
                S(0)        S(0)         S(1)]

mx3 = @SMatrix [S(1)  S(0)        S(0);
                S(0)  S("cosθ₃")  S("-sinθ₃");
                S(0)  S("sinθ₃")  S("cosθ₃")]

my3 = @SMatrix [S("cosθ₃")   S(0)  S("sinθ₃");
                S(0)         S(1)  S(0);
                S("-sinθ₃")  S(0)  S("cosθ₃")]

mz3 = @SMatrix [S("cosθ₃")  S("-sinθ₃")  S(0);
                S("sinθ₃")  S("cosθ₃")   S(0);
                S(0)        S(0)         S(1)]

v = @SVector [S("v[1]"), S("v[2]"), S("v[3]")]

myx = my1 * mx2
mxy = mx1 * my2
mxz = mx1 * mz2
mzx = mz1 * mx2
mzy = mz1 * my2
myz = my1 * mz2

myxy = my1 * mx2 * my3
myxz = my1 * mx2 * mz3
mxyx = mx1 * my2 * mx3
mxyz = mx1 * my2 * mz3
mxzx = mx1 * mz2 * mx3
mxzy = mx1 * mz2 * my3
mzxz = mz1 * mx2 * mz3
mzxy = mz1 * mx2 * my3
mzyz = mz1 * my2 * mz3
mzyx = mz1 * my2 * mx3
myzy = my1 * mz2 * my3
myzx = my1 * mz2 * mx3

export S, mx1, my1, mz1, mx2, my2, mz2, mx3, my3, mz3, v, mxy, myx, mxz, mzx, mzy, myz,
    myxy, myxz, mxyx, mxyz, mxzx, mxzy, mzxz, mzxy, mzyz, mzyx, myzy, myzx

end # module


text:
using SafeTestsets: @safetestset

@time @safetestset "rfam" begin include("rfam.jl") end


text:
export ×, ⊂

# pure operators
set_in = LogicalOperation((a, b)->set_includes(a, b), :∈, 2, false, false)
Base.:(∈)(x, S::MathematicalSet) = LogicalExpression(AbstractExpression[x, S], set_in)
Base.:(∉)(x, S::MathematicalSet) = ¬(x ∈ S)
# subset operator
⊂ = LogicalOperation((a, b)->false, :⊂, 2)

# ×
function ×(A::MathematicalSet, B::MathematicalSet)
    @unique_symbols a b
    IntensionalSet(orderedpair(a, b), (a ∈ A, b ∈ B))
end
function ×(x, A::MathematicalSet)
    @unique_symbols a
    IntensionalSet(orderedpair(x, a), a ∈ A)
end
function ×(A::MathematicalSet, x)
    @unique_symbols a
    IntensionalSet(orderedpair(a, x), a ∈ A)
end

# ∪
function Base.:(∪)(A::MathematicalSet, B::MathematicalSet)
    @unique_symbols a
    IntensionalSet(a, (a ∈ A) ∨ (a ∈ B))
end

# ∩
function Base.:(∩)(A::MathematicalSet, B::MathematicalSet)
    @unique_symbols a
    IntensionalSet(a, (a ∈ A) ∧ (a ∈ B))
end

# \ (set difference)
function Base.:(\)(A::MathematicalSet, B::MathematicalSet)
    @unique_symbols a
    IntensionalSet(a, (a ∈ A) ∧ (a ∉ B))
end


text:
using GR

# Simple particle simulation

const N = 300                   # number of particles
M = 0.05 * ones(Float64, N)     # masses
const S = 0.04                  # size of particles

function step(dt, p, v)
  # update positions
  p[:,:] += dt * v[:,:]

  # find pairs of particles undergoing a collision
  for i in 1:N
    for j in i+1:N
      dx = p[i,1] - p[j,1]
      dy = p[i,2] - p[j,2]
      d = sqrt(dx*dx + dy*dy)

      if d < 2*S
        # relative location & velocity vectors
        r_rel = p[i] - p[j]
        v_rel = v[i] - v[j]

        # momentum vector of the center of mass
        v_cm = (M[i] * v[i] + M[j] * v[j]) / (M[i] + M[j])

        # collisions of spheres reflect v_rel over r_rel
        rr_rel = dot(r_rel, r_rel)
        vr_rel = dot(v_rel, r_rel)
        v_rel = 2 * r_rel * vr_rel / rr_rel - v_rel

        # assign new velocities
        v[i] = v_cm + v_rel * M[j] / (M[i] + M[j])
        v[j] = v_cm - v_rel * M[i] / (M[i] + M[j])
      end
    end
  end

  # check for crossing boundary
  for i in 1:N
    if p[i,1] < -2 + S
      p[i,1] = -2 + S
      v[i,1] *= -1
    elseif p[i,1] > 2 - S
      p[i,1] = 2 - S
      v[i,1] *= -1
    end
    if p[i,2] < -2 + S
      p[i,2] = -2 + S
      v[i,2] *= -1
    elseif p[i,2] > 2 - S
      p[i,2] = 2 - S
      v[i,2] *= -1
    end
  end

  return p, v
end

setwindow(-2, 2, -2, 2)
setviewport(0, 1, 0, 1)
setmarkertype(GR.MARKERTYPE_SOLID_CIRCLE)
setmarkersize(1.0)

n = 0
t = 0.0

p = (rand(N,2) - 0.5) * (4-2*S)     # initial positions
v = rand(N,2) - 0.5                 # initial velocities

const dt = 1.0 / 30

while t < 3

  start = tic()
  p, v = step(dt, p, v)

  clearws()
  setmarkercolorind(983)
  polymarker(p[:,1], p[:,2])

  if n > 0
    text(0.01, 0.95, @sprintf("%10s: %4d fps", "Julia", round(n / t)))
  end
  updatews()

  n += 1
  t += 1.0 / (toq() * 1000000)

end


text:
function newbound_rounding(U,V,bmatchval)
  U_sortperm = sortcolsperm(U,true)
  V_sortperm = sortcolsperm(V,true)
  nU = size(U,1)
  nV = size(V,1)
  r = size(U,2)
  @assert r==size(V,2)

  d = min(size(U_sortperm,1),size(V_sortperm,1))
  U_weights = sort(U,1,rev=true)
  V_weights = sort(V,1,rev=true)

#   U_weights = U_weights[1:d,:]
#   V_weights = V_weights[1:d,:]
#
#   U_sortperm = U_sortperm[1:d,:]
#   V_sortperm = V_sortperm[1:d,:]

  P = spzeros(nU,nV)
  allrecoveries = zeros(size(U,2))
  for i = 1:size(U,2)
    ui = U_weights[:,i]
    vi = V_weights[:,i]
    lastid_ui = findfirst(ui.<0)
    lastid_vi = findfirst(vi.<0)

    if lastid_ui == 0 && lastid_vi == 0
      lastidpos = d
      lneg = -1
    elseif lastid_vi == 0
      lastidpos = min(d,lastid_ui-1)
      lneg = -1
    elseif lastid_ui == 0
      lastidpos = min(d,lastid_vi-1)
      lneg = -1
    else
      lastidpos = min(lastid_ui,lastid_vi)-1
      lneg = min(nU-lastid_ui,nV-lastid_vi)
    end

    ei1 = U_sortperm[1:lastidpos,i]
    ej1 = V_sortperm[1:lastidpos,i]
    ei2 = U_sortperm[nU-lneg:nU,i]
    ej2 = V_sortperm[nV-lneg:nV,i]
    ei = vcat(ei1,ei2)
    ej = vcat(ej1,ej2)
    # allrecoveries[i] = evaluate_erdosreyni_experiment(A,B,ei,ej)
    P = P + sparse(ei,ej,1,nU,nV) #+ generate_b_match_overlapping(ei,ej,2,nU,nV)
    # with bmatching:
    # bmatchval = 20
    for bm = 1:bmatchval
    if !isempty(ej)
      popfirst!(ej)
      P = P + sparse(ei[1:end-bm],ej,1,nU,nV)
    end
    end

  end
  @show nnz(P)/prod(size(P))
  return P
end

function newbound_rounding_lowrank_evaluation_relaxed(U,V,bmatchval)
  # U = Float32.(U)
  # V = Float32.(V)
  U_sortperm = sortcolsperm(U,true)
  V_sortperm = sortcolsperm(V,true)
  nU = size(U,1)
  nV = size(V,1)
  r = size(U,2)
  @assert r == size(V,2)

  d = min(size(U_sortperm,1),size(V_sortperm,1))
  U_weights = sort(U,dims=1,rev=true)
  V_weights = sort(V,dims=1,rev=true)

#   U_weights = U_weights[1:d,:]
#   V_weights = V_weights[1:d,:]
#
#   U_sortperm = U_sortperm[1:d,:]
#   V_sortperm = V_sortperm[1:d,:]

#   P = spzeros(nU,nV)
  U1 = []
  V1 = []
  allrecoveries = zeros(size(U,2))
  for i = 1:size(U,2)
    ui = U_weights[:,i]
    vi = V_weights[:,i]
    lastid_ui = findfirst(ui.<0)
    lastid_vi = findfirst(vi.<0)

    lastid_ui = (lastid_ui === nothing) ? 0 : lastid_ui
    lastid_vi = (lastid_vi === nothing) ? 0 : lastid_vi

    if lastid_ui == 0 && lastid_vi == 0
      lastidpos = d
      lneg = -1
    elseif lastid_vi == 0
      lastidpos = min(d,lastid_ui-1)
      lneg = -1
    elseif lastid_ui == 0
      lastidpos = min(d,lastid_vi-1)
      lneg = -1
    else
      lastidpos = min(lastid_ui,lastid_vi)-1
      lneg = min(nU-lastid_ui,nV-lastid_vi)
    end

    ei1 = U_sortperm[1:lastidpos,i]
    ej1 = V_sortperm[1:lastidpos,i]
    ei2 = U_sortperm[nU-lneg:nU,i]
    ej2 = V_sortperm[nV-lneg:nV,i]
    ei = vcat(ei1,ei2)
    ej = vcat(ej1,ej2)
    # allrecoveries[i] = evaluate_erdosreyni_experiment(A,B,ei,ej)
#     P = P + sparse(ei,ej,1,nU,nV)

    # P = P + sparse(ei,ej,1,nU,nV) #+ generate_b_match_overlapping(ei,ej,2,nU,nV)
    # with bmatching:
    append!(U1,ei)
    append!(V1,ej)
    # bmatchval = 6
    #println("bmatchval is $bmatchval")
    for bm = 1:bmatchval
      if !isempty(ej)
        popfirst!(ej)
        # P = P + sparse(ei[1:end-bm],ej,1,nU,nV)
        append!(U1,ei[1:end-bm])
        append!(V1,ej)
      end
    end
  end

  all_matches = [U1 V1]
  unique_matches = unique(all_matches,dims=1)
  U1unique = unique_matches[:,1]
  V1unique = unique_matches[:,2]
  uo = U[U1unique,:]
  vo = V[V1unique,:]
  weights = vec(sum(uo.*vo,dims=2))
  X = sparse(U1unique,V1unique,weights,nU,nV)

  return X
end

max_stars_repo_name: chuckcol/LambdaTAME

max_issues_repo_name: chuckcol/LambdaTAME

max_forks_repo_name: chuckcol/LambdaTAME


text:
struct GenericVendor <: AbstractVendor end
struct GenericProduct <: AbstractProduct end

mutable struct CSVImporter{T <: AbstractSchema} <: AbstractImporter
    isenabled::Bool
    schema::T
end

function isenabled(importer::CSVImporter)
    return importer.isenabled
end

function enable!(importer::CSVImporter)
    importer.isenabled = true
end

function disable!(importer::CSVImporter)
    importer.isenabled  = false
end

function (load::CSVImporter{<:AbstractSchema})(path::Path)
    @show "Inside generic CSV importer..."
    disable!(load)
    load_dataframe(path)
end

function load_dataframe(path::Path)
    path₁ = joinpath(get_directory(path), get_filename(path))
    data = nothing
    if is_readable_file(path₁)
        try
            data = CSV.File(path₁) |> DataFrame
        catch e
            println(e)
            CImGui.OpenPopup("Has the file been corrupted?")
        end
    else
        CImGui.OpenPopup("Do you have permission to read the file?")
    end
    handle_import_error_messages()
    data
end

function load_dataframe(path::Path, header::Vector{String})
    path₁ = joinpath(get_directory(path), get_filename(path))
    data = nothing
    if is_readable_file(path₁)
        try
            data = CSV.File(path₁; header = header) |> DataFrame
        catch e
            println(e)
            CImGui.OpenPopup("Has the file been corrupted?")
        end
    else
        CImGui.OpenPopup("Do you have permission to read the file?")
    end
    handle_import_error_messages()
    data
end

function handle_import_error_messages()
    if CImGui.BeginPopupModal("Has the file been corrupted?", C_NULL, CImGui.ImGuiWindowFlags_AlwaysAutoResize)
        CImGui.Text("Unable to open the specified file.\nPlease verify that: \n   (1) the file has not been corrupted; \n   (2) you have permission to access the file.\n\n")
        CImGui.Separator()
        CImGui.Button("OK", (120, 0)) && CImGui.CloseCurrentPopup()
        CImGui.SetItemDefaultFocus()
        CImGui.EndPopup()
    end

    if CImGui.BeginPopupModal("Do you have permission to read the file?", C_NULL, CImGui.ImGuiWindowFlags_AlwaysAutoResize)
        CImGui.Text("Unable to access the specified file.\nPlease verify that: \n   (1) the file exists; \n   (2) you have permission to read the file.\n\n")
        CImGui.Separator()
        CImGui.Button("OK", (120, 0)) && CImGui.CloseCurrentPopup()
        CImGui.SetItemDefaultFocus()
        CImGui.EndPopup()
    end

    if CImGui.BeginPopupModal("Have you opened the appropriate file?", C_NULL, CImGui.ImGuiWindowFlags_AlwaysAutoResize)
    CImGui.Text("The file format does not to conform to the expected schema.\nPlease verify that the data in the file follows the expected convention.\n\n")
    CImGui.Separator()
    CImGui.Button("OK", (120, 0)) && CImGui.CloseCurrentPopup()
    CImGui.SetItemDefaultFocus()
    CImGui.EndPopup()
end
end

# Import Labelled Intervals
function (load::CSVImporter{<:CSVSchema{<: GenericVendor, <: GenericProduct, <: IntervalLabels}})(path::Path)
    disable!(load)
    df = load_dataframe(path)
    _, ncol = size(df)
    if ncol != 6
        CImGui.OpenPopup("Have you opened the appropriate file?")
        return nothing
    else
        dict = Dict{String, LabelledInterval}()
        for row in eachrow(df)
            label = row[:label]
            start = Float64(row[:start])
            stop = Float64(row[:stop])
            i₀ = row[:first]
            s = row[:step]
            i₁ = row[:last]
            interval = i₀:s:i₁
            dict[label] = LabelledInterval(label, NestedInterval(start = start, stop = stop, interval = interval))
        end
        return LabelledIntervals("Conditions", dict)
    end
end

struct ImportContext{T₁ <: AbstractDialogControl,   T₂ <: AbstractDialogModel,  T₃ <: AbstractDisplayProperties, T₄ <: Union{AbstractImporter}} <: AbstractContext
    control::T₁
    model::T₂
    display_properties::T₃
    action::T₄
end

function (context::ImportContext)()
        control = context.control
        model = context.model
        display_properties = context.display_properties
        action = context.action
        isenabled(control) ? control(model, display_properties, action) : nothing
        data = isenabled(action) ? action(get_path(model, ConfirmedStatus())) : nothing
end

function isrunning(context::ImportContext)
        isenabled(context.control) || isenabled(context.action)
end

function enable!(context::ImportContext)
    enable!(context.control)
end


text:
module ACEds

# Write your package code here.

#include("./equivRot3d.jl")
include("./rpi/rpi.jl")
#include("./test.jl")
end


text:
#! format: off
@testset "Common" begin
    xs =[3,5,1,6,2,9,8,nothing, nothing,7,4]
    t = TreeNode(3,
            TreeNode(5,
                TreeNode(6),
                TreeNode(2,
                    TreeNode(7),
                    TreeNode(4)
                )
            ),
            TreeNode(1,
                TreeNode(9),
                TreeNode(8)
            )
        )
    @test t == convert(TreeNode{Int}, xs)
    @test t == TreeNode{Int}(xs)
    @test TreeNode{Int}([6,3,5,nothing,2,0,nothing,nothing,1]) ==
    TreeNode(6, TreeNode(3, nothing, TreeNode(2, nothing, TreeNode(1))), TreeNode(5, TreeNode(0)))
end


text:
using Baseball
using Test

@testset "Baseball.jl" begin
end


text:
@with_kw struct cylinder_inputs{T <: AbstractFloat}
    μsp::T = 10.0
    μa::T = 0.1
    n_ext::T = 1.0  #surrounding index of refraction
    n_med::T = 1.4 # layers index of refraction

    #source, detector
    lz::T = 5.0 # length of cylinder
    ρ::T = 1.0 #
    a::T = 5.0 # radius of cylinder
    z::T = 0.0
    pw::T = 0.1

    ω::T = 0.0 #
end

# look at get_afac for some reason gives type instabillity
function diffusionparams(μsp, n_med, n_ext)
    ## Diffusion parameters
    D = 1/3μsp
    ν = 29.9792345/n_med
    A = 1.0 #get_afac(n_ext/n_med) # need to calculate reflection between layers and surrounding medium
    zb = 2*A*D
    z0 = 1/(μsp)

    return D, ν, A, zb, z0
end

function _greencylin(α, D, z, z0, zb, lz)
    g = (exp(-α*abs(z - z0)) - exp(-α*(z + z0 + 2*zb)))
    g -= exp(α*(z0 + z - 2*lz - 2*zb))*(1 - exp(-2*α*(z0 + zb)))*(1 - exp(-2*α*(zb + z)))/(1 - exp(-2*α*(lz + 2*zb)))
    return g / (2*D*α)
end

function _greencylin_CW(sn, μa, D, z, z0, zb, lz)
    α = sqrt(μa/D + sn^2)
    return _greencylin(α, D, z, z0, zb, lz)
end

function _greencylin_FD(sn, μa, D, z, z0, zb, lz, ν, ω)
    α = sqrt(μa/D + sn^2 + im*ω/(D*ν))
    return _greencylin(α, D, z, z0, zb, lz)
end

function _greencylin_Laplace(sn, μa, D, z, z0, zb, lz, ν, s)
    α = sqrt(μa/D + sn^2 + s/(D*ν))
    return _greencylin(α, D, z, z0, zb, lz)
end

function fluence_DA_cylinder_CW(data::cylinder_inputs, besselroots)

    D, ν, A, zb, z0 = diffusionparams(data.μsp, data.n_med, data.n_ext)

    ϕ = zero(eltype(data.ρ))
    for bessel in besselroots
        ϕ += _greencylin_CW(bessel/(data.a + zb), data.μa, D, data.z, z0, zb, data.lz)*besselj0(bessel/(data.a + zb)*data.ρ)/(besselj1(bessel))^2
    end

    return ϕ / (π*(data.a + zb)^2)
end

function fluence_DA_cylinder_FD(data::cylinder_inputs, besselroots)

    D, ν, A, zb, z0 = diffusionparams(data.μsp, data.n_med, data.n_ext)

    ϕ = zero(Complex{eltype(data.ρ)})

    for n in eachindex(besselroots)
        ϕ += _greencylin_FD(besselroots[n]/(data.a + zb), data.μa, D, data.z, z0, zb, data.lz, ν, data.ω)*besselj0(besselroots[n]/(data.a + zb)*data.ρ)/(besselj1(besselroots[n]))^2
    end

    return ϕ / (π*(data.a + zb)^2)
end

### Compute Time Domain Signal

# calulcate fluence in laplace space: s -> im*ω
function _fluence_DA_cylinder_Laplace(s, data::cylinder_inputs, besselroots)

    D, ν, A, zb, z0 = diffusionparams(data.μsp, data.n_med, data.n_ext)

    ϕ = zero(eltype(s))

    for bessel in besselroots
        ϕ += _greencylin_Laplace(bessel/(data.a + zb), data.μa, D, data.z, z0, zb, data.lz, ν, s)*besselj0(bessel/(data.a + zb)*data.ρ)/(besselj1(bessel))^2
    end

    return ϕ/(π*(data.a + zb)^2)
end

# use adaptive hyperbola contour for each time point Laplace Transform
# much more accurate but slower. N = 8 is usually sufficient
# need larger N for very short times (< 0.1 ns) and very long times (> 3ns)
function fluence_DA_cylinder_TD(t, N, data::cylinder_inputs, besselroots)
    Rt = zeros(Float64, length(t))
    Rt = LT_hyperbola(s -> _fluence_DA_cylinder_Laplace(s, data, besselroots), N, t)

    return Rt
end

# use fixed contour for every time point
# signifanctly faster (10-100x) than adaptive contour but less accurate
# utilizing this is usually ok but further testing is needed
# N = 28 is usuall pretty accurate
# Need about 600 Bessel Zeros
# the 59.9 is dividing it by some constant to match the fluence in semi-inf. Unsure why this is needed (not needed in the steady-state: match exactly)
function fluence_DA_cylinder_TD(t, N, data::cylinder_inputs, besselroots)
    Rt = zeros(Float64, length(t))
    Rt = LT_hyper_fixed(s -> _fluence_DA_cylinder_Laplace(s, data, besselroots), N, t)

    return Rt
end

### for a finite flat beam incident onto center of cylinder top
# have data structure contain pw (radius of beam)

function fluence_DA_cylinder_CW_beam(data::cylinder_inputs1, besselroots)

    D, ν, A, zb, z0 = diffusionparams(data.μsp, data.n_med, data.n_ext)

    ϕ = zero(eltype(data.ρ))
    sn = zero(eltype(data.ρ))

    for n in eachindex(besselroots)
        sn = besselroots[n]/(data.a + zb)
        ϕ += _greencylin_CW(sn, data.μa, D, data.z, z0, zb, data.lz)*besselj0(sn*data.ρ)*besselj1(sn*data.pw)/(sn*(besselj1(sn*(data.a + zb))^2))
    end

    return 2*ϕ/(data.pw*π*(data.a + zb)^2)
end

function _fluence_DA_cylinder_Laplace_beam(s, data::cylinder_inputs1, besselroots)

    D, ν, A, zb, z0 = diffusionparams(data.μsp, data.n_med, data.n_ext)
    ϕ = zero(eltype(data.ρ))

    sn = zero(eltype(data.ρ))
    for n in eachindex(besselroots)
        sn = besselroots[n]/(data.a + zb)
        ϕ += _greencylin_Laplace(sn, data.μa, D, data.z, z0, zb, data.lz, ν, s)*besselj0(sn*data.ρ)*besselj1(sn*data.pw)/(sn*(besselj1(sn*(data.a + zb))^2))
    end

    return 2*ϕ/(data.pw*π*(data.a + zb)^2)
end

function fluence_DA_cylinder_TD_beam(t, N, data::cylinder_inputs1, besselroots)
    Rt = zeros(Float64, length(t))
    Rt = LT_hyper_fixed(s -> _fluence_DA_cylinder_Laplace_beam(s, data, besselroots), N, t)

    return Rt
end

### arbitrary position

@with_kw struct cylinder_inputs_arb3
    μsp::Float64 = 10.0
    μa::Float64 = 0.1
    n_ext::Float64 = 1.0  #surrounding index of refraction
    n_med::Float64 = 1.4 # layers index of refraction

    #source, detector
    lz::Float64 = 5.0 # length of cylinder
    ρ::Float64 = 5.0 #
    ρ0::Float64 = 5.0
    a::Float64 = 5.0 # radius of cylinder
    z::Float64 = 2.0
    z0::Float64 = 3.0

    θ = 0.0

    ω::Float64 = 0.0 #
end

function fluence_DA_cylinder_CW(data::cylinder_inputs_arb2, besselj_roots)

    D, ν, A, zb, z0 = diffusionparams(data.μsp, data.n_med, data.n_ext)

    ϕ = zero(eltype(data.ρ))
    ϕ1 = zero(eltype(data.ρ))

    for m = 0:100
        ϕ1 = 0.0
        for bessel in besselj_roots[:, abs(m) + 1] # 1 column is 0 order, bessel roots symmetric J(-nu, x) = (-1)^2 J(nu, x)
            ϕ1 += _greencylin_CW(bessel/(data.a + zb), data.μa, D, data.z, z0, zb, data.lz)*besselj(m, bessel/(data.a + zb)*data.ρ)*besselj(m, bessel/(data.a + zb)*data.ρ0)/(besselj(m + 1, bessel))^2

        end
        ϕ1 *= ((2 - (m == 0 ? 1 : 0))*cos(m*data.θ))

        ϕ += ϕ1
    end

    return ϕ / (π*(data.a + zb)^2)
end

function fluence_DA_cylinder_CW(data::cylinder_inputs_arb3, besselj_roots)

    D, ν, A, zb, z0 = diffusionparams(data.μsp, data.n_med, data.n_ext)

    ρ0 = data.ρ0 -  1/data.μsp[1]
    z0 = data.z0

    ϕ = 0.0
    ϕ1 = 0.0

    for m = 0:1200
        ϕ1 = 0.0
        for bessel in besselj_roots[:, abs(m) + 1] # 1 column is 0 order, bessel roots symmetric J(-nu, x) = (-1)^2 J(nu, x)
            ϕ1 += _greencylin_CW(bessel/(data.a + zb), data.μa, D, data.z, z0, zb, data.lz)*besselj(m, bessel/(data.a + zb)*data.ρ)*besselj(m, bessel/(data.a + zb)*ρ0)/(besselj(m + 1, bessel))^2

        end
        #println(ϕ1)
        ϕ1 *= ((2 - (m == 0 ? 1 : 0))*cos(m*data.θ))

        ϕ += ϕ1
        #println(ϕ)
    end

    return ϕ / (π*(data.a + zb)^2)
end

## work in process: need to test best convergence
## can't break for order m when theta is non-zero. You get oscillating term
##appears to be a bug when z = z0. Slow convergence or something off? Seems when theta = 0 everything is ok...
function fluence_DA_cylinder_CW(data::cylinder_inputs_arb3, besselj_roots)

    D, ν, A, zb, z0 = diffusionparams(data.μsp, data.n_med, data.n_ext)

    ρ0 = data.ρ0 -  1/data.μsp[1]
    z0 = data.z0

    ϕ = 0.0
    ϕ1 = 0.0

    root_old = 0.0
    order_old = 0.0
    count = 0

    TOL = 1e-12

    for m = -1600:1600#0:8000
        ϕ1 = 0.0
        for bessel in besselj_roots[:, abs(m) + 1] # 1 column is 0 order, bessel roots symmetric J(-nu, x) = (-1)^2 J(nu, x)
            root_old = ϕ1
            ϕ1 += _greencylin_CW(bessel/(data.a + zb), data.μa, D, data.z, z0, zb, data.lz)*besselj(m, bessel/(data.a + zb)*data.ρ)*besselj(m, bessel/(data.a + zb)*ρ0)/(besselj(m + 1, bessel))^2

            if (abs(ϕ1 - root_old))/ϕ1 < TOL
                #println(bessel)
                break
            end
        end

        ϕ1 *= #=((2 - (m == 0 ? 1 : 0))*=#cos(m*data.θ)
        order_old = ϕ1

        ϕ += ϕ1
    end

    return ϕ / (π*(data.a + zb)^2)
end

function _fluence_DA_cylinder_Laplace(s, data::cylinder_inputs_arb3, besselj_roots)

    D, ν, A, zb, z0 = diffusionparams(data.μsp, data.n_med, data.n_ext)

    ρ0 = data.ρ0 -  1/data.μsp[1]
    z0 = data.z0

    ϕ = 0.0
    ϕ1 = 0.0

    for m = 0:600
        ϕ1 = 0.0
        for bessel in besselj_roots[:, abs(m) + 1] # 1 column is 0 order, bessel roots symmetric J(-nu, x) = (-1)^2 J(nu, x)
            ϕ1 += _greencylin_Laplace(bessel/(data.a + zb), data.μa, D, data.z, z0, zb, data.lz, ν, s)*besselj(m, bessel/(data.a + zb)*data.ρ)*besselj(m, bessel/(data.a + zb)*ρ0)/(besselj(m + 1, bessel))^2

        end
        ϕ1 *= ((2 - (m == 0 ? 1 : 0))*cos(m*data.θ))

        ϕ += ϕ1
    end

    return ϕ / (π*(data.a + zb)^2)
end
function fluence_DA_cylinder_TD(t, N, data::cylinder_inputs_arb3, besselroots)
    Rt = zeros(Float64, length(t))
    Rt = LT_hyper_fixed(s -> _fluence_DA_cylinder_Laplace(s, data, besselroots), N, t)

    return Rt
end
#=
Calculate Bessel Roots

using GSL

s = 1:1200 # number of roots to calculate
m = 0:1300 # order of besselj
besselj_roots = [sf_bessel_zero_Jnu_e(m, s).val for s in s, m in m]

=#

max_issues_repo_name: heltonmc/LightPropagation


text:
using Test
using CLIMA
using StaticArrays
using LinearAlgebra
using KernelAbstractions
using CLIMA.MPIStateArrays: device

include("ode_tests_common.jl")

CLIMA.init()
const ArrayType = CLIMA.array_type()

@testset "ODE Solvers" begin
    @testset "Convergence/extensive" begin
        @testset "1-rate ODE" begin
            function rhs!(dQ, Q, ::Nothing, time; increment)
                if increment
                    dQ .+= Q * cos(time)
                else
                    dQ .= Q * cos(time)
                end
            end
            exactsolution(q0, time) = q0 * exp(sin(time))

            @testset "Explicit methods" begin
                finaltime = 20.0
                dts = [2.0^(-k) for k in 0:7]
                errors = similar(dts)
                q0 =
                    ArrayType === Array ? [1.0] : range(-1.0, 1.0, length = 303)
                for (method, expected_order) in explicit_methods
                    for (n, dt) in enumerate(dts)
                        Q = ArrayType(q0)
                        solver = method(rhs!, Q; dt = dt, t0 = 0.0)
                        solve!(Q, solver; timeend = finaltime)
                        Q = Array(Q)
                        errors[n] =
                            maximum(@. abs(Q - exactsolution(q0, finaltime)))
                    end
                    rates = log2.(errors[1:(end - 1)] ./ errors[2:end])
                    @test isapprox(rates[end], expected_order; atol = 0.17)
                end
            end
        end

        @testset "Two-rate ODE with a linear stiff part" begin
            c = 100.0
            function rhs_full!(dQ, Q, ::Nothing, time; increment)
                if increment
                    dQ .+= im * c * Q .+ exp(im * time)
                else
                    dQ .= im * c * Q .+ exp(im * time)
                end
            end

            function rhs_nonlinear!(dQ, Q, ::Nothing, time; increment)
                if increment
                    dQ .+= exp(im * time)
                else
                    dQ .= exp(im * time)
                end
            end
            rhs_slow! = rhs_nonlinear!

            function rhs_linear!(dQ, Q, ::Nothing, time; increment)
                if increment
                    dQ .+= im * c * Q
                else
                    dQ .= im * c * Q
                end
            end
            rhs_fast! = rhs_linear!

            function exactsolution(q0, time)
                q0 * exp(im * c * time) +
                (exp(im * time) - exp(im * c * time)) / (im * (1 - c))
            end

            @testset "IMEX methods" begin
                finaltime = pi / 2
                dts = [2.0^(-k) for k in 2:13]
                errors = similar(dts)

                q0 = ArrayType <: Array ? [1.0] : range(-1.0, 1.0, length = 303)
                for (method, expected_order) in imex_methods
                    for split_explicit_implicit in (false, true)
                        for variant in (LowStorageVariant(), NaiveVariant())
                            for (n, dt) in enumerate(dts)
                                Q = ArrayType{ComplexF64}(q0)
                                rhs! =
                                    split_explicit_implicit ? rhs_nonlinear! :
                                    rhs_full!
                                solver = method(
                                    rhs!,
                                    rhs_linear!,
                                    LinearBackwardEulerSolver(
                                        DivideLinearSolver(),
                                        isadjustable = true,
                                    ),
                                    Q;
                                    dt = dt,
                                    t0 = 0.0,
                                    split_explicit_implicit = split_explicit_implicit,
                                    variant = variant,
                                )
                                solve!(Q, solver; timeend = finaltime)
                                Q = Array(Q)
                                errors[n] = maximum(@. abs(
                                    Q - exactsolution(q0, finaltime),
                                ))
                            end

                            rates = log2.(errors[1:(end - 1)] ./ errors[2:end])
                            @test errors[1] < 2.0
                            @test isapprox(
                                rates[end],
                                expected_order;
                                atol = 0.1,
                            )
                        end
                    end
                end
            end

            @testset "MRRK methods (no substeps)" begin
                finaltime = pi / 2
                dts = [2.0^(-k) for k in 2:11]
                errors = similar(dts)
                for (slow_method, slow_expected_order) in slow_mrrk_methods
                    for (fast_method, fast_expected_order) in fast_mrrk_methods
                        q0 = ArrayType === Array ? [1.0] :
                            range(-1.0, 1.0, length = 303)
                        for (n, dt) in enumerate(dts)
                            Q = ArrayType{ComplexF64}(q0)
                            solver = MultirateRungeKutta(
                                (
                                    slow_method(rhs_slow!, Q),
                                    fast_method(rhs_fast!, Q),
                                );
                                dt = dt,
                                t0 = 0.0,
                            )
                            solve!(Q, solver; timeend = finaltime)
                            Q = Array(Q)
                            errors[n] = maximum(@. abs(
                                Q - exactsolution(q0, finaltime),
                            ))
                        end

                        rates = log2.(errors[1:(end - 1)] ./ errors[2:end])
                        min_order =
                            min(slow_expected_order, fast_expected_order)
                        max_order =
                            max(slow_expected_order, fast_expected_order)
                        @test (
                            isapprox(rates[end], min_order; atol = 0.1) ||
                            isapprox(rates[end], max_order; atol = 0.1) ||
                            min_order <= rates[end] <= max_order
                        )
                    end
                end
            end

            @testset "MRRK methods (with substeps)" begin
                finaltime = pi / 2
                dts = [2.0^(-k) for k in 2:15]
                errors = similar(dts)
                for (slow_method, slow_expected_order) in slow_mrrk_methods
                    for (fast_method, fast_expected_order) in fast_mrrk_methods
                        q0 = ArrayType === Array ? [1.0] :
                            range(-1.0, 1.0, length = 303)
                        for (n, fast_dt) in enumerate(dts)
                            slow_dt = c * fast_dt
                            Q = ArrayType{ComplexF64}(q0)
                            solver = MultirateRungeKutta((
                                slow_method(rhs_slow!, Q; dt = slow_dt),
                                fast_method(rhs_fast!, Q; dt = fast_dt),
                            ))
                            solve!(Q, solver; timeend = finaltime)
                            Q = Array(Q)
                            errors[n] = maximum(@. abs(
                                Q - exactsolution(q0, finaltime),
                            ))
                        end

                        rates = log2.(errors[1:(end - 1)] ./ errors[2:end])
                        min_order =
                            min(slow_expected_order, fast_expected_order)
                        max_order =
                            max(slow_expected_order, fast_expected_order)
                        @test (
                            isapprox(rates[end], min_order; atol = 0.1) ||
                            isapprox(rates[end], max_order; atol = 0.1) ||
                            min_order <= rates[end] <= max_order
                        )
                    end
                end
            end

            @testset "MIS methods (with substeps)" begin
                finaltime = pi / 2
                dts = [2.0^(-k) for k in 2:11]
                errors = similar(dts)
                for (mis_method, mis_expected_order) in mis_methods
                    for fast_method in (LSRK54CarpenterKennedy,)
                        q0 = ArrayType === Array ? [1.0] :
                            range(-1.0, 1.0, length = 303)
                        for (n, dt) in enumerate(dts)
                            Q = ArrayType{ComplexF64}(q0)
                            solver = mis_method(
                                rhs_slow!,
                                rhs_fast!,
                                fast_method,
                                4,
                                Q;
                                dt = dt,
                                t0 = 0.0,
                            )
                            solve!(Q, solver; timeend = finaltime)
                            Q = Array(Q)
                            errors[n] = maximum(@. abs(
                                Q - exactsolution(q0, finaltime),
                            ))
                        end
                        rates = log2.(errors[1:(end - 1)] ./ errors[2:end])
                        @test isapprox(
                            rates[end],
                            mis_expected_order;
                            atol = 0.1,
                        )
                    end
                end
            end
        end

        #=
        Test problem (4.2) from RobertsSarsharSandu2018arxiv
        @article{RobertsSarsharSandu2018arxiv,
          title={Coupled Multirate Infinitesimal GARK Schemes for Stiff Systems with
                 Multiple Time Scales},
          author={Roberts, Steven and Sarshar, Arash and Sandu, Adrian},
          journal={arXiv preprint arXiv:1812.00808},
          year={2019}
        }

        Note: The actual rates are all over the place with this test and passing largely
              depends on final dt size
        =#
        @testset "2-rate ODE from RobertsSarsharSandu2018arxiv" begin
            ω = 100
            λf = -10
            λs = -1
            ξ = 1 // 10
            α = 1
            ηfs = ((1 - ξ) / α) * (λf - λs)
            ηsf = -ξ * α * (λf - λs)
            Ω = @SMatrix [
                λf ηfs
                ηsf λs
            ]

            function rhs_fast!(dQ, Q, param, t; increment)
                @inbounds begin
                    increment || (dQ .= 0)
                    yf = Q[1]
                    ys = Q[2]
                    gf = (-3 + yf^2 - cos(ω * t)) / 2yf
                    gs = (-2 + ys^2 - cos(t)) / 2ys
                    dQ[1] += Ω[1, 1] * gf + Ω[1, 2] * gs - ω * sin(ω * t) / 2yf
                end
            end

            function rhs_slow!(dQ, Q, param, t; increment)
                @inbounds begin
                    increment || (dQ .= 0)
                    yf = Q[1]
                    ys = Q[2]
                    gf = (-3 + yf^2 - cos(ω * t)) / 2yf
                    gs = (-2 + ys^2 - cos(t)) / 2ys
                    dQ[2] += Ω[2, 1] * gf + Ω[2, 2] * gs - sin(t) / 2ys
                end
            end

            exactsolution(t) = [sqrt(3 + cos(ω * t)); sqrt(2 + cos(t))]

            @testset "MRRK methods (no substeps)" begin
                finaltime = 5π / 2
                dts = [2.0^(-k) for k in 2:8]
                error = similar(dts)
                for (slow_method, slow_expected_order) in slow_mrrk_methods
                    for (fast_method, fast_expected_order) in fast_mrrk_methods
                        for (n, dt) in enumerate(dts)
                            Q = exactsolution(0)
                            solver = MultirateRungeKutta(
                                (
                                    slow_method(rhs_slow!, Q),
                                    fast_method(rhs_fast!, Q),
                                );
                                dt = dt,
                                t0 = 0.0,
                            )
                            solve!(Q, solver; timeend = finaltime)
                            error[n] = norm(Q - exactsolution(finaltime))
                        end

                        rate = log2.(error[1:(end - 1)] ./ error[2:end])
                        min_order =
                            min(slow_expected_order, fast_expected_order)
                        max_order =
                            max(slow_expected_order, fast_expected_order)
                        @test (
                            isapprox(rate[end], min_order; atol = 0.3) ||
                            isapprox(rate[end], max_order; atol = 0.3) ||
                            min_order <= rate[end] <= max_order
                        )
                    end
                end
            end

            @testset "MRRK methods (with substeps)" begin
                finaltime = 5π / 2
                dts = [2.0^(-k) for k in 2:9]
                error = similar(dts)
                for (slow_method, slow_expected_order) in slow_mrrk_methods
                    for (fast_method, fast_expected_order) in fast_mrrk_methods
                        for (n, fast_dt) in enumerate(dts)
                            Q = exactsolution(0)
                            slow_dt = ω * fast_dt
                            solver = MultirateRungeKutta((
                                slow_method(rhs_slow!, Q; dt = slow_dt),
                                fast_method(rhs_fast!, Q; dt = fast_dt),
                            ))
                            solve!(Q, solver; timeend = finaltime)
                            error[n] = norm(Q - exactsolution(finaltime))
                        end

                        rate = log2.(error[1:(end - 1)] ./ error[2:end])
                        min_order =
                            min(slow_expected_order, fast_expected_order)
                        max_order =
                            max(slow_expected_order, fast_expected_order)
                        @test (
                            isapprox(rate[end], min_order; atol = 0.3) ||
                            isapprox(rate[end], max_order; atol = 0.3) ||
                            min_order <= rate[end] <= max_order
                        )
                    end
                end
            end

            @testset "MRRK methods with IMEX fast solver" begin
                function rhs_zero!(dQ, Q, param, t; increment)
                    if !increment
                        dQ .= 0
                    end
                end

                finaltime = 5π / 2
                dts = [2.0^(-k) for k in 2:9]
                error = similar(dts)
                for (slow_method, slow_expected_order) in slow_mrrk_methods
                    for (fast_method, fast_expected_order) in imex_methods
                        for (n, fast_dt) in enumerate(dts)
                            Q = exactsolution(0)
                            slow_dt = ω * fast_dt
                            solver = MultirateRungeKutta((
                                slow_method(rhs_slow!, Q; dt = slow_dt),
                                fast_method(
                                    rhs_fast!,
                                    rhs_zero!,
                                    LinearBackwardEulerSolver(
                                        DivideLinearSolver(),
                                        isadjustable = true,
                                    ),
                                    Q;
                                    dt = fast_dt,
                                    split_explicit_implicit = false,
                                ),
                            ))
                            solve!(Q, solver; timeend = finaltime)
                            error[n] = norm(Q - exactsolution(finaltime))
                        end

                        rate = log2.(error[1:(end - 1)] ./ error[2:end])
                        min_order =
                            min(slow_expected_order, fast_expected_order)
                        max_order =
                            max(slow_expected_order, fast_expected_order)
                        atol = fast_method == ARK2GiraldoKellyConstantinescu ?
                            0.5 : 0.37
                        @test (
                            isapprox(rate[end], min_order; atol = atol) ||
                            isapprox(rate[end], max_order; atol = atol) ||
                            min_order <= rate[end] <= max_order
                        )
                    end
                end
            end

            @testset "MIS methods (with substeps)" begin
                finaltime = 5π / 2
                dts = [2.0^(-k) for k in 2:10]
                error = similar(dts)
                for (mis_method, mis_expected_order) in mis_methods
                    for fast_method in (LSRK54CarpenterKennedy,)
                        for (n, dt) in enumerate(dts)
                            Q = exactsolution(0)
                            solver = mis_method(
                                rhs_slow!,
                                rhs_fast!,
                                fast_method,
                                4,
                                Q;
                                dt = dt,
                                t0 = 0.0,
                            )
                            solve!(Q, solver; timeend = finaltime)
                            error[n] = norm(Q - exactsolution(finaltime))
                        end

                        rate = log2.(error[1:(end - 1)] ./ error[2:end])
                        @test isapprox(
                            rate[end],
                            mis_expected_order;
                            atol = 0.1,
                        )
                    end
                end
            end
        end

        # Simple 3-rate problem based on test of RobertsSarsharSandu2018arxiv
        #
        # NOTE: Since we have no theory to say this ODE solver is accurate, the rates
        #      suggest that things are really only 2nd order.
        #
        # TODO: This is not great, but no theory to say we should be accurate!
        @testset "3-rate ODE" begin
            ω1, ω2, ω3 = 10000, 100, 1
            λ1, λ2, λ3 = -100, -10, -1
            β1, β2, β3 = 2, 3, 4

            ξ12 = λ2 / λ1
            ξ13 = λ3 / λ1
            ξ23 = λ3 / λ2

            α12, α13, α23 = 1, 1, 1

            η12 = ((1 - ξ12) / α12) * (λ1 - λ2)
            η13 = ((1 - ξ13) / α13) * (λ1 - λ3)
            η23 = ((1 - ξ23) / α23) * (λ2 - λ3)

            η21 = ξ12 * α12 * (λ2 - λ1)
            η31 = ξ13 * α13 * (λ3 - λ1)
            η32 = ξ23 * α23 * (λ3 - λ2)

            Ω = @SMatrix [
                λ1 η12 η13
                η21 λ2 η23
                η31 η32 λ3
            ]

            function rhs1!(dQ, Q, param, t; increment)
                @inbounds begin
                    increment || (dQ .= 0)
                    y1, y2, y3 = Q[1], Q[2], Q[3]
                    g = @SVector [
                        (-β1 + y1^2 - cos(ω1 * t)) / 2y1,
                        (-β2 + y2^2 - cos(ω2 * t)) / 2y2,
                        (-β3 + y3^2 - cos(ω3 * t)) / 2y3,
                    ]
                    dQ[1] += Ω[1, :]' * g - ω1 * sin(ω1 * t) / 2y1
                end
            end
            function rhs2!(dQ, Q, param, t; increment)
                @inbounds begin
                    increment || (dQ .= 0)
                    y1, y2, y3 = Q[1], Q[2], Q[3]
                    g = @SVector [
                        (-β1 + y1^2 - cos(ω1 * t)) / 2y1,
                        (-β2 + y2^2 - cos(ω2 * t)) / 2y2,
                        (-β3 + y3^2 - cos(ω3 * t)) / 2y3,
                    ]
                    dQ[2] += Ω[2, :]' * g - ω2 * sin(ω2 * t) / 2y2
                end
            end
            function rhs3!(dQ, Q, param, t; increment)
                @inbounds begin
                    increment || (dQ .= 0)
                    y1, y2, y3 = Q[1], Q[2], Q[3]
                    g = @SVector [
                        (-β1 + y1^2 - cos(ω1 * t)) / 2y1,
                        (-β2 + y2^2 - cos(ω2 * t)) / 2y2,
                        (-β3 + y3^2 - cos(ω3 * t)) / 2y3,
                    ]
                    dQ[3] += Ω[3, :]' * g - ω3 * sin(ω3 * t) / 2y3
                end
            end
            function rhs12!(dQ, Q, param, t; increment)
                rhs1!(dQ, Q, param, t; increment = increment)
                rhs2!(dQ, Q, param, t; increment = true)
            end

            exactsolution(t) = [
                sqrt(β1 + cos(ω1 * t)),
                sqrt(β2 + cos(ω2 * t)),
                sqrt(β3 + cos(ω3 * t)),
            ]

            @testset "MRRK methods (no substeps)" begin
                finaltime = π / 2
                dts = [2.0^(-k) for k in 2:10]
                error = similar(dts)
                for (rate3_method, rate3_order) in slow_mrrk_methods
                    for (rate2_method, rate2_order) in slow_mrrk_methods
                        for (rate1_method, rate1_order) in fast_mrrk_methods
                            for (n, dt) in enumerate(dts)
                                Q = exactsolution(0)
                                solver = MultirateRungeKutta(
                                    (
                                        rate3_method(rhs3!, Q),
                                        rate2_method(rhs2!, Q),
                                        rate1_method(rhs1!, Q),
                                    );
                                    dt = dt,
                                    t0 = 0.0,
                                )
                                solve!(Q, solver; timeend = finaltime)
                                error[n] = norm(Q - exactsolution(finaltime))
                            end

                            rate = log2.(error[1:(end - 1)] ./ error[2:end])
                            min_order =
                                min(rate3_order, rate2_order, rate1_order)
                            max_order =
                                max(rate3_order, rate2_order, rate1_order)
                            @test 2 <= rate[end]
                        end
                    end
                end
            end

            @testset "MRRK methods (with substeps)" begin
                finaltime = π / 2
                dts = [2.0^(-k) for k in 10:17]
                error = similar(dts)
                for (rate3_method, rate3_order) in slow_mrrk_methods
                    for (rate2_method, rate2_order) in slow_mrrk_methods
                        for (rate1_method, rate1_order) in fast_mrrk_methods
                            for (n, dt1) in enumerate(dts)
                                Q = exactsolution(0)
                                dt2 = (ω1 / ω2) * dt1
                                dt3 = (ω2 / ω3) * dt2
                                solver = MultirateRungeKutta((
                                    rate3_method(rhs3!, Q; dt = dt3),
                                    rate2_method(rhs2!, Q; dt = dt2),
                                    rate1_method(rhs1!, Q; dt = dt1),
                                ))
                                solve!(Q, solver; timeend = finaltime)
                                error[n] = norm(Q - exactsolution(finaltime))
                            end

                            rate = log2.(error[1:(end - 1)] ./ error[2:end])
                            min_order =
                                min(rate3_order, rate2_order, rate1_order)
                            max_order =
                                max(rate3_order, rate2_order, rate1_order)

                            @test 2 <= rate[end]
                        end
                    end
                end
            end
        end

        #=
        Test problem (8.2) from Sandu (2019) for MRI-GARK Schemes
            @article{Sandu2019,
                title={A class of multirate infinitesimal gark methods},
                author={Sandu, Adrian},
                journal={SIAM Journal on Numerical Analysis},
                volume={57},
                number={5},
                pages={2300--2327},
                year={2019},
                publisher={SIAM},
                doi={10.1137/18M1205492}
            }
        =#
        @testset "2-rate problem" begin
            ω = 20
            λf = -10
            λs = -1
            ξ = 1 // 10
            α = 1
            ηfs = ((1 - ξ) / α) * (λf - λs)
            ηsf = -ξ * α * (λf - λs)
            Ω = @SMatrix [
                λf ηfs
                ηsf λs
            ]

            function rhs_fast!(dQ, Q, param, t; increment)
                @kernel function knl!(dQ, Q, t, increment)
                    @inbounds begin
                        increment || (dQ .= 0)
                        yf = Q[1]
                        ys = Q[2]
                        gf = (-3 + yf^2 - cos(ω * t)) / 2yf
                        gs = (-2 + ys^2 - cos(t)) / 2ys
                        dQ[1] +=
                            Ω[1, 1] * gf + Ω[1, 2] * gs - ω * sin(ω * t) / 2yf
                    end
                end
                event = Event(device(Q))
                event = knl!(device(Q), 1)(
                    dQ,
                    Q,
                    t,
                    increment;
                    ndrange = 1,
                    dependencies = (event,),
                )
                wait(device(Q), event)
            end

            function rhs_slow!(dQ, Q, param, t; increment)
                @kernel function knl!(dQ, Q, t, increment)
                    @inbounds begin
                        increment || (dQ .= 0)
                        yf = Q[1]
                        ys = Q[2]
                        gf = (-3 + yf^2 - cos(ω * t)) / 2yf
                        gs = (-2 + ys^2 - cos(t)) / 2ys
                        dQ[2] += Ω[2, 1] * gf + Ω[2, 2] * gs - sin(t) / 2ys
                    end
                end
                event = Event(device(Q))
                event = knl!(device(Q), 1)(
                    dQ,
                    Q,
                    t,
                    increment;
                    ndrange = 1,
                    dependencies = (event,),
                )
                wait(device(Q), event)
            end

            struct ODETestConvNonLinBE <: AbstractBackwardEulerSolver end
            ODESolvers.Δt_is_adjustable(::ODETestConvNonLinBE) = true
            function (::ODETestConvNonLinBE)(Q, Qhat, α, p, t)
                @kernel function knl!(Q, Qhat, α, p, t)
                    @inbounds begin
                        # Slow RHS has zero tendency of yf so just copy Qhat
                        Q[1] = yf = Qhat[1]

                        gf = (-3 + yf^2 - cos(ω * t)) / 2yf

                        # solves: Q = Qhat[2] + α * rhs_slow(Q, t)
                        # (simplifies to a quadratic equation)
                        a = 2 - α * Ω[2, 2]
                        b = -2 * (Qhat[2] + α * Ω[2, 1] * gf)
                        c = α * (Ω[2, 2] * (2 + cos(t)) + sin(t))
                        Q[2] = (-b + sqrt(b^2 - 4 * a * c)) / (2a)
                    end
                end
                event = Event(device(Q))
                event = knl!(device(Q), 1)(
                    Q,
                    Qhat,
                    α,
                    p,
                    t;
                    ndrange = 1,
                    dependencies = (event,),
                )
                wait(device(Q), event)
            end

            exactsolution(t) =
                ArrayType([sqrt(3 + cos(ω * t)); sqrt(2 + cos(t))])

            finaltime = 1
            dts = [2.0^(-k) for k in 2:7]
            error = similar(dts)
            @testset "Explicit" begin
                for (mri_method, mri_expected_order) in mrigark_erk_methods
                    for (fast_method, fast_expected_order) in
                        fast_mrigark_methods
                        for (n, slow_dt) in enumerate(dts)
                            Q = exactsolution(0)
                            fast_dt = slow_dt / ω
                            fastsolver = fast_method(rhs_fast!, Q; dt = fast_dt)
                            solver = mri_method(
                                rhs_slow!,
                                fastsolver,
                                Q,
                                dt = slow_dt,
                            )
                            solve!(Q, solver; timeend = finaltime)
                            error[n] = norm(Q - exactsolution(finaltime))
                        end

                        rate = log2.(error[1:(end - 1)] ./ error[2:end])
                        order = mri_expected_order
                        @test isapprox(
                            rate[end],
                            mri_expected_order;
                            atol = 0.3,
                        )
                    end
                end
            end

            @testset "Implicit" begin
                for (mri_method, mri_expected_order) in mrigark_irk_methods
                    for (fast_method, fast_expected_order) in
                        fast_mrigark_methods
                        for (n, slow_dt) in enumerate(dts)
                            Q = exactsolution(0)
                            fast_dt = slow_dt / ω
                            fastsolver = fast_method(rhs_fast!, Q; dt = fast_dt)
                            solver = mri_method(
                                rhs_slow!,
                                ODETestConvNonLinBE(),
                                fastsolver,
                                Q,
                                dt = slow_dt,
                            )
                            solve!(Q, solver; timeend = finaltime)
                            error[n] = norm(Q - exactsolution(finaltime))
                        end

                        rate = log2.(error[1:(end - 1)] ./ error[2:end])
                        order = mri_expected_order
                        @test isapprox(
                            rate[end],
                            mri_expected_order;
                            atol = 0.3,
                        )
                    end
                end
            end

            @testset "IMEX" begin
                for (method, expected_order) in imex_methods
                    for (n, dt) in enumerate(dts)
                        Q = exactsolution(0)
                        solver = method(
                            rhs_fast!,
                            rhs_slow!,
                            ODETestConvNonLinBE(),
                            Q;
                            dt = dt,
                            t0 = 0.0,
                            split_explicit_implicit = true,
                            variant = NaiveVariant(),
                        )
                        solve!(Q, solver; timeend = finaltime)
                        error[n] = norm(Q - exactsolution(finaltime))
                    end
                    rate = log2.(error[1:(end - 1)] ./ error[2:end])
                    order = expected_order
                    @test isapprox(rate[end], expected_order; atol = 0.3)
                end
            end
        end

        # 3-rate modification of the above test problem
        @testset "3-rate problem" begin
            ω1, ω2, ω3 = 20, 5, 1
            λ1, λ2, λ3 = -20, -5, -1
            β1, β2, β3 = 2, 2, 2

            ξ12 = λ2 / (λ1 + λ2)
            ξ13 = λ3 / (λ1 + λ3)
            ξ23 = λ3 / (λ2 + λ3)

            α12, α13, α23 = 1, 1, 1

            η12 = ((1 - ξ12) / α12) * (λ1 - λ2)
            η13 = 0#((1 - ξ13) / α13) * (λ1 - λ3)
            η23 = ((1 - ξ23) / α23) * (λ2 - λ3)

            η21 = ξ12 * α12 * (λ2 - λ1)
            η31 = 0#ξ13 * α13 * (λ3 - λ1)
            η32 = ξ23 * α23 * (λ3 - λ2)

            Ω = @SMatrix [
                λ1 η12 η13
                η21 λ2 η23
                η31 η32 λ3
            ]

            function rhs1!(dQ, Q, param, t; increment)
                @kernel function knl!(dQ, Q, t, increment)
                    @inbounds begin
                        increment || (dQ .= 0)
                        y1, y2, y3 = Q[1], Q[2], Q[3]
                        g = @SVector [
                            (-β1 + y1^2 - cos(ω1 * t)) / 2y1,
                            (-β2 + y2^2 - cos(ω2 * t)) / 2y2,
                            (-β3 + y3^2 - cos(ω3 * t)) / 2y3,
                        ]
                        dQ[1] += Ω[1, :]' * g - ω1 * sin(ω1 * t) / 2y1
                    end
                end
                event = Event(device(Q))
                event = knl!(device(Q), 1)(
                    dQ,
                    Q,
                    t,
                    increment;
                    ndrange = 1,
                    dependencies = (event,),
                )
                wait(device(Q), event)
            end
            function rhs2!(dQ, Q, param, t; increment)
                @kernel function knl!(dQ, Q, t, increment)
                    @inbounds begin
                        increment || (dQ .= 0)
                        y1, y2, y3 = Q[1], Q[2], Q[3]
                        g = @SVector [
                            (-β1 + y1^2 - cos(ω1 * t)) / 2y1,
                            (-β2 + y2^2 - cos(ω2 * t)) / 2y2,
                            (-β3 + y3^2 - cos(ω3 * t)) / 2y3,
                        ]
                        dQ[2] += Ω[2, :]' * g - ω2 * sin(ω2 * t) / 2y2
                    end
                end
                event = Event(device(Q))
                event = knl!(device(Q), 1)(
                    dQ,
                    Q,
                    t,
                    increment;
                    ndrange = 1,
                    dependencies = (event,),
                )
                wait(device(Q), event)
            end
            function rhs3!(dQ, Q, param, t; increment)
                @kernel function knl!(dQ, Q, t, increment)
                    @inbounds begin
                        increment || (dQ .= 0)
                        y1, y2, y3 = Q[1], Q[2], Q[3]
                        g = @SVector [
                            (-β1 + y1^2 - cos(ω1 * t)) / 2y1,
                            (-β2 + y2^2 - cos(ω2 * t)) / 2y2,
                            (-β3 + y3^2 - cos(ω3 * t)) / 2y3,
                        ]
                        dQ[3] += Ω[3, :]' * g - ω3 * sin(ω3 * t) / 2y3
                    end
                end
                event = Event(device(Q))
                event = knl!(device(Q), 1)(
                    dQ,
                    Q,
                    t,
                    increment;
                    ndrange = 1,
                    dependencies = (event,),
                )
                wait(device(Q), event)
            end
            struct ODETestConvNonLinBE3Rate <: AbstractBackwardEulerSolver end
            ODESolvers.Δt_is_adjustable(::ODETestConvNonLinBE3Rate) = true
            function (::ODETestConvNonLinBE3Rate)(Q, Qhat, α, p, t)
                @kernel function knl!(Q, Qhat, α, p, t)
                    @inbounds begin
                        # Slower RHS has zero tendency of yf so just copy Qhat
                        Q[1] = y1 = Qhat[1]
                        Q[2] = y2 = Qhat[2]

                        g = @SVector [
                            (-β1 + y1^2 - cos(ω1 * t)) / 2y1,
                            (-β2 + y2^2 - cos(ω2 * t)) / 2y2,
                        ]

                        # solves: Q = Qhat + α * rhs_slow(Q, t)
                        # (simplifies to a quadratic equation)
                        a = 2 - α * Ω[3, 3]
                        b =
                            -2 *
                            (Qhat[3] + α * Ω[3, 1] * g[1] + α * Ω[3, 2] * g[2])
                        c = α * (Ω[3, 3] * (β3 + cos(t)) + sin(t))
                        Q[3] = (-b + sqrt(b^2 - 4 * a * c)) / (2a)
                    end
                end
                event = Event(device(Q))
                event = knl!(device(Q), 1)(
                    Q,
                    Qhat,
                    α,
                    p,
                    t;
                    ndrange = 1,
                    dependencies = (event,),
                )
                wait(device(Q), event)
            end

            exactsolution(t) = ArrayType([
                sqrt(β1 + cos(ω1 * t)),
                sqrt(β2 + cos(ω2 * t)),
                sqrt(β3 + cos(ω3 * t)),
            ])

            finaltime = 1
            dts = [2.0^(-k) for k in 1:7]
            error = similar(dts)
            @testset "Explicit" begin
                for (slow_method, slow_order) in mrigark_erk_methods
                    for (mid_method, mid_order) in mrigark_erk_methods
                        for (fast_method, fast_order) in fast_mrigark_methods
                            for (n, slow_dt) in enumerate(dts)
                                Q = exactsolution(0)
                                mid_dt = slow_dt / ω2
                                fast_dt = slow_dt / ω1
                                fastsolver = fast_method(rhs1!, Q; dt = fast_dt)
                                midsolver = mid_method(
                                    rhs2!,
                                    fastsolver,
                                    Q,
                                    dt = mid_dt,
                                )
                                slowsolver = slow_method(
                                    rhs3!,
                                    midsolver,
                                    Q,
                                    dt = slow_dt,
                                )
                                solve!(Q, slowsolver; timeend = finaltime)
                                error[n] = norm(Q - exactsolution(finaltime))
                            end

                            rate = log2.(error[1:(end - 1)] ./ error[2:end])
                            expected_order =
                                min(slow_order, mid_order, fast_order)
                            @test isapprox(
                                rate[end],
                                expected_order;
                                atol = 0.3,
                            )
                        end
                    end
                end
            end
            @testset "Implicit-Explicit" begin
                for (slow_method, slow_order) in mrigark_irk_methods
                    for (mid_method, mid_order) in mrigark_erk_methods
                        for (fast_method, fast_order) in fast_mrigark_methods
                            for (n, slow_dt) in enumerate(dts)
                                Q = exactsolution(0)
                                mid_dt = slow_dt / ω2
                                fast_dt = slow_dt / ω1
                                fastsolver = fast_method(rhs1!, Q; dt = fast_dt)
                                midsolver = mid_method(
                                    rhs2!,
                                    fastsolver,
                                    Q,
                                    dt = mid_dt,
                                )
                                slowsolver = slow_method(
                                    rhs3!,
                                    ODETestConvNonLinBE3Rate(),
                                    midsolver,
                                    Q,
                                    dt = slow_dt,
                                )
                                solve!(Q, slowsolver; timeend = finaltime)
                                error[n] = norm(Q - exactsolution(finaltime))
                            end

                            rate = log2.(error[1:(end - 1)] ./ error[2:end])
                            expected_order =
                                min(slow_order, mid_order, fast_order)
                            @test isapprox(
                                rate[end],
                                expected_order;
                                atol = 0.4,
                            )
                        end
                    end
                end
            end
        end
    end
end


text:
using IntelVectorMath
using Distributions, Statistics, BenchmarkTools # for benchmark
using Plots # for plotting
using JLD2, FileIO # to save file

cd(dirname(@__FILE__))
include(joinpath(dirname(dirname(@__FILE__)), "test", "common.jl"))

################################################################
complex = !isempty(ARGS) && ARGS[1] == "complex"
complex = false

# First generate some random data and test functions in Base on it
const NVALS = 10_000
base_unary = complex ? base_unary_complex : base_unary_real
base_binary = complex ? base_binary_complex : base_binary_real
types = complex ? (Complex64, Complex128) : (Float32, Float64)

# arrays of inputs are stored in a Tuple. So later for calling use inp... to get the content of the Tuple
input = Dict( t =>
[
 [(randindomain(t, NVALS, domain),) for (_, _, domain) in base_unary];
 [(randindomain(t, NVALS, domain1), randindomain(t, NVALS, domain2)) for (_, _, domain1, domain2) in base_binary];
 # (randindomain(t, NVALS, (0, 100)), randindomain(t, 1, (-1, 20))[1])
]
    for t in types)

fns = [[x[1:2] for x in base_unary_real];
       [x[1:2] for x in base_binary_real]]
################################################################

"""
    bench(fns, input)

benchmark function for IntelVectorMath.jl. Calls both Base and IntelVectorMath functions and stores the benchmarks in two nested Dict. First layer specifies type, and second layer specifies the function name. The result is a Tuple, 1st element being benchmark for Base/SpecialFunctions and 2nd element being for IntelVectorMath.

# Examples
```julia
fns = [(:Base, :acos); (:Base, :atan); (:SpecialFunctions, :ref)] # array of tuples
input = Dict( Float64 => [(rand(1000)); (rand(1000), rand(1000)); (rand(1000))]) # Dict of array of tuples
times = bench(fns, input)

times[Float64][:acos][1] # Base.acos benchmark for Float64
times[Float64][:acos][2] # IntelVectorMath.acos benchmark for Float64
```
"""
function bench(fns, input)
    Dict(t => begin
        Dict( fn[2] => begin
            base_fn = eval(:($(fn[1]).$(fn[2])))
            vml_fn = eval(:(IntelVectorMath.$(fn[2])))
            println("benchmarking $vml_fn for type $t")
            timesBase = @benchmark $base_fn.($inp...)
            timesVML = @benchmark $vml_fn($inp...)

            [timesBase, timesVML]
        end for (fn, inp) in zip(fns, input[t]) )
    end for t in types)
end
################################################################

# do benchmark
benches = bench(fns, input)

@save "benchmarkData.jld" benches
# @save "benchmarkData-complex.jld" benches

# benches = load("benchmarkData.jld2", "benches")

# something is wrong with these
deleteat!(fns, [18, 31])
delete!(benches[Float32], :lgamma)
delete!(benches[Float64], :lgamma)
delete!(benches[Float32], :log10)
delete!(benches[Float64], :log10)

################################################################

"""
benchmark analysis function
"""
function ratioci(y, x, alpha = 0.05)
    tq² = abs2(quantile(TDist(length(x) + length(y) - 2), alpha))
    μx = mean(x)
    σx² = varm(x, μx)
    μy = mean(y)
    σy² = varm(y, μy)
    a = sqrt((μx * μy)^2 - (μx^2 - tq² * σx²) * (μy^2 - tq² * σy²))
    b = μx^2 - tq² * σx²
    return (((μx * μy) - a) / b, ((μx * μy) + a) / b)
end

################################################################

"""
Does analysis of the benchmark data and plots them as bars.
"""
function plotBench()

# Print ratio
    colors = [:blue, :red]
    for itype = 1:length(types)

        # creating arrays of times from benchmarks
        benchVals = collect(values(benches[types[itype]]))
        builtint = [x[1].times for x in benchVals]
        vmlt = [x[2].times for x in benchVals]

        # calculating mean of run times
        μ = vec(map(mean, builtint) ./ map(mean, vmlt))

        # error bar disabled because benchmark tools takes care of errors

        # # calculating error bars
        # ci = zeros(Float64, 2, length(fns))
        # for ifn = 1:length(builtint)
        #     lower, upper = ratioci(builtint[ifn], vmlt[ifn])
        #     ci[1, ifn] = μ[ifn] - lower
        #     ci[2, ifn] = upper - μ[ifn]
        # end

        # adding bar
        plt = bar!(
            0.2+(0.4*itype):length(fns),
            μ,
            # yerror = ci,
            fillcolor = colors[itype],
            labels = [string(x) for x in types],
            dpi = 600
        )
    end
    fname = [string(fn[2]) for fn in fns]
    # if !complex
    #     fname[end-1] = "A.^B"
    #     fname[end] = "A.^b"
    # end
    xlims!(0, length(fns) + 1)
    xticks!(1:length(fns)+1, fname, rotation = 70, fontsize = 10)
    title!("IntelVectorMath Performance for array of size $NVALS")
    ylabel!("Relative Speed (IntelVectorMath/Base)")
    hline!([1], line=(4, :dash, 0.6, [:green]), labels = 1)
    savefig("performance$(complex ? "_complex" : "").png")

end

################################################################

# do plotting
plotBench()


text:
using Revise
using SpineOpt

run_spineopt(ARGS...)


text:
# Shaded line plots
function PlotShaded(label_id, x::AbstractArray{<:Real}, y::AbstractArray{<:Real}, count,
                    y_ref::Real, args...)
    (x, y) = promote(x, y)
    return PlotShaded(label_id, x, y, count, Float64(y_ref), args...)
end

function PlotShaded(label_id, x::AbstractArray{<:Real}, y1::AbstractArray{<:Real},
                    y2::AbstractArray{<:Real}, args...)
    return PlotShaded(label_id, promote(x, y1, y2)..., args...)
end

function PlotShaded(x::AbstractArray{T}, y1::AbstractArray{T}, y2::AbstractArray{T};
                    label_id::String="",
                    count::Integer=min(length(x), length(y2), length(y2)),
                    offset::Integer=0, stride::Integer=1) where {T<:ImPlotData}
    return PlotShaded(label_id, x, y1, y2, count, offset, stride * sizeof(T))
end

function PlotShaded(x::AbstractArray{T1}, y1::AbstractArray{T2}, y2::AbstractArray{T3};
                    kwargs...) where {T1<:Real,T2<:Real,T3<:Real}
    return PlotShaded(promote(x, y1, y2)...; kwargs...)
end

function PlotShaded(x::AbstractArray{T}, y::AbstractArray{T}, y_ref::T;
                    count::Integer=min(length(x), length(y)), offset::Integer=0,
                    stride::Integer=1, label_id::String="") where {T<:ImPlotData}
    return PlotShaded(label_id, x, y, count, y_ref, offset, stride * sizeof(T))
end

function PlotShaded(x::AbstractArray{T1}, y::AbstractArray{T2}, y_ref::T3;
                    kwargs...) where {T1<:Real,T2<:Real,T3<:Real}
    (x, y) = promote(x, y)
    y_ref = eltype(x)(y_ref)
    return PlotShaded(x, y, y_ref; kwargs...)
end

function PlotShaded(x::UnitRange{<:Integer}, y1::AbstractArray{T}, y2::AbstractArray{T};
                    label_id::String="") where {T<:ImPlotData}
    count = length(x)
    offset = x.start >= 1 ? x.start - 1 : throw("Offset out of bounds")
    x = collect(T, x)
    return PlotShaded(label_id, x, y1, y2, count, offset, sizeof(T))
end

function PlotShaded(x::UnitRange{<:Integer}, y::AbstractArray{T}, y_ref::T;
                    label_id::String="") where {T<:ImPlotData}
    count = length(x)
    offset = x.start >= 1 ? x.start - 1 : throw("Offset out of bounds")
    x = collect(T, x)
    return PlotShaded(label_id, x, y, count, y_ref, offset, sizeof(T))
end

function PlotShaded(x::StepRange, y1::AbstractArray{T}, y2::AbstractArray{T};
                    label_id::String="") where {T<:ImPlotData}
    x.stop < 1 && throw("Range out of bounds")
    count = length(x)
    offset = x.start >= 1 ? x.start - 1 : throw("Range out of bounds")
    stride = x.step >= 1 ? x.step : throw("Stride must be greater than zero.")
    x = collect(T, x)
    return PlotShaded(label_id, x, y1, y2, count, offset, stride * sizeof(T))
end

function PlotShaded(x::StepRange, y::AbstractArray{T}, y_ref::T;
                    label_id::String="") where {T<:ImPlotData}
    x.stop < 1 && throw("Range out of bounds")
    count = length(x)
    offset = x.start >= 1 ? x.start - 1 : throw("Range out of bounds")
    stride = x.step >= 1 ? x.step : throw("Stride must be greater than zero.")

    x = collect(T, x)
    return PlotShaded(label_id, x, y, count, y_ref, offset, stride * sizeof(T))
end

function PlotShaded(x::OrdinalRange, y1::AbstractArray{T1}, y2::AbstractArray{T2};
                    kwargs...) where {T1<:Real,T2<:Real}
    return PlotShaded(x, promote(y1, y2)...; kwargs...)
end

function PlotShaded(x::OrdinalRange, y::AbstractArray{T1}, y_ref::T2;
                    kwargs...) where {T1<:Real,T2<:Real}
    y_ref = eltype(y)(y_ref)
    return PlotShaded(x, y, y_ref; kwargs...)
end

# xfield, yfield should be propertynames of eltype(structvec)
function PlotShaded(structvec::Vector{T}, xfield::Symbol, yfield::Symbol, y_ref::Real;
                    count::Integer=length(structvec), offset::Integer=0, stride::Integer=1,
                    label_id::String="") where {T}
    Tx = fieldtype(T, xfield)
    Ty = fieldtype(T, yfield)
    x_offset = fieldoffset(T, Base.fieldindex(T, xfield))
    y_offset = fieldoffset(T, Base.fieldindex(T, yfield))
    x_ptr = Ptr{Tx}((pointer(structvec, 1) + x_offset))
    y_ptr = Ptr{Ty}((pointer(structvec, 1) + y_offset))

    if !T.mutable
        # this is somewhat illegal and is used only to pass a pointer through AbstractArray argument into ccall
        x = unsafe_wrap(Vector{Tx}, x_ptr, size(structvec); own=false)
        y = unsafe_wrap(Vector{Ty}, y_ptr, size(structvec); own=false)
        stride = stride * sizeof(T)
    else # two new vectors every 1/60 second...
        x = Vector{Tx}(undef, length(structvec))
        y = Vector{Ty}(undef, length(structvec))
        for (i, val) in enumerate(structvec)
            x[i] = getproperty(val, xfield)
            y[i] = getproperty(val, yfield)
        end
        if Tx !== Ty
            x, y = promote(x, y)
        end
        stride = stride * sizeof(eltype(x))
    end

    return PlotShaded(label_id, x, y, count, y_ref, offset, stride)
end


text:
#! /usr/bin/env julia

using MeshSlicer
using Lint
using Base.Test
using ImmutableArrays

# Make pwd MeshSlicer/test
data_path = Pkg.dir("MeshSlicer")*"/test/data/"

# getmesh/getface test
println("Testing PolygonMesh...")

# get the first face from cube.stl
testFace = Face(Vector3[Vector3(10.0,0.0,10.0),
                      Vector3(10.0,0.0,0.0),
                      Vector3(10.0,10.0,10.0)],
                      Vector3(1.0,0.0,0.0))

testBounds = Bounds3{Float64}(10.0,10.0,10.0,0,0,0)

resultmesh = PolygonMesh(data_path*"cube.stl")

@test resultmesh.faces[end].vertices == testFace.vertices
@test resultmesh.faces[end].normal == testFace.normal

println("Testing Bounds...")
@test resultmesh.bounds.xmax == testBounds.xmax
@test resultmesh.bounds.ymax == testBounds.ymax
@test resultmesh.bounds.zmax == testBounds.zmax
@test resultmesh.bounds.xmin == testBounds.xmin
@test resultmesh.bounds.ymin == testBounds.ymin
@test resultmesh.bounds.zmin == testBounds.zmin
@test resultmesh.bounds == testBounds

# Test bounds on translated cube
# translate([1,2,3])cube([4,5,6]);
testBounds2 = Bounds3(4.0,6.0,9.0,1.0,2.0,4.0)
updateFace = Face(Vector3[Vector3(5.0,3.0,4.0),
                      Vector3(2.0,7.0,4.0),
                      Vector3(2.0,3.0,3.0)],
                      Vector3(-1.0,0.0,0.0)) # contains correct bounds for cube
update!(testBounds2, updateFace)
resultmesh2 = PolygonMesh(data_path*"translated_cube.stl")

@test resultmesh2.bounds.xmax == testBounds2.xmax
@test resultmesh2.bounds.ymax == testBounds2.ymax
@test resultmesh2.bounds.zmax == testBounds2.zmax
@test resultmesh2.bounds.xmin == testBounds2.xmin
@test resultmesh2.bounds.ymin == testBounds2.ymin
@test resultmesh2.bounds.zmin == testBounds2.zmin
@test resultmesh2.bounds == testBounds2

# Test patio detection
println("Tesing patio detection...")

@test resultmesh2.patios == [3.0, 9.0]
@test resultmesh.patios == [0.0, 10.0]

# Test LineSegment
testpoint1 = Vector3[Vector3(0.0,0.0,0.0),Vector3(0.0,1.0,0.0),Vector3(0.0,1.0,1.0)]
normal = Vector3(1.0, 1.0, 1.0)
l1 = LineSegment(testpoint1[3],testpoint1[1],testpoint1[2], 0.5, normal)
l2 = LineSegment(testpoint1[3],testpoint1[2],testpoint1[1], 0.5, normal)
l3 = LineSegment(testpoint1[3],testpoint1[2],testpoint1[1], 0.0, normal)

testpoint2 = Vector3[Vector3(0.0,0.0,0.0),Vector3(1.0,1.0,1.0),Vector3(0.0,1.0,0.0)]
l4 = LineSegment(testpoint2[2],testpoint2[1],testpoint2[3], 0.5, normal)
l5 = LineSegment(testpoint2[2],testpoint2[3],testpoint2[1], 0.5, normal)

l6 = LineSegment(testFace, 5.0)
l7 = LineSegment(testpoint2[2],testpoint2[3],testpoint2[1], 1.0, normal)

testpoint3 = Vector3[Vector3(0.0,0.0,0.0),Vector3(0.0,1.0,1.0),Vector3(0.0,0.0,2.0)]
testface3 = Face(testpoint3, normal)
l8 = LineSegment(testface3, 0.9)
l9 = LineSegment(Face(Vector3[testpoint3[3], testpoint3[1], testpoint3[2]], normal), 1.5)

@test l1.start == Vector2(0.0, 0.5)
@test l1.finish == Vector2(0.0, 1.0)
@test l2.start == Vector2(0.0, 1.0)
@test l2.finish == Vector2(0.0, 0.5)
@test l3.start == Vector2(0.0, 1.0)
@test l3.finish == Vector2(0.0, 0.0)
@test l4.start == Vector2(0.5, 0.5)
@test l4.finish == Vector2(0.5, 1.0)
@test l5.start == Vector2(0.5, 1.0)
@test l6.start == Vector2(10.0, 0.0)
@test l6.finish == Vector2(10.0, 5.0)
@test l6.normal == Vector3(1.0,0.0,0.0)
@test l7.start == Vector2(1.0, 1.0)
@test l7.finish == Vector2(1.0, 1.0)
@test l8.start == Vector2(0.0, 0.0)
@test l8.finish == Vector2(0.0, 0.9)
@test l9.start == Vector2(0.0, 0.0)
@test l9.finish == Vector2(0.0, 0.5)

# test binary STL
println("Testing Binary STL import...")
binarymesh = PolygonMesh(data_path*"cube_binary.stl")
@test binarymesh.bounds == Bounds3(10.0,10.0,10.0,0.0,0.0,0.0)

println("Testing slicing...")
slice1 = MeshSlice(binarymesh, [0.0,5.0,10.0])
for slice in slice1
    for poly in slice.polygons
        n = length(poly.segments)
        for i = 1:n-1
            @test poly.segments[i].finish == poly.segments[i+1].start
        end
    end
end

println("Tesing Polygon construction...")
ep = 0.1
l1 = LineSegment(Vector2(0.0+ep,0.0), Vector2(0.0,1.0+ep), Vector3(0.0,0.0,0.0))
l2 = LineSegment(Vector2(0.0+ep,1.0), Vector2(0.5-ep,1.0), Vector3(0.0,0.0,0.0))
l3 = LineSegment(Vector2(0.5+ep,1.0), Vector2(0.0,0.0+ep), Vector3(0.0,0.0,0.0))
@test length(Polygon([l1,l2,l3], eps=ep*2, autoeps=false)) == 1
@test length(Polygon([l1,l2,l3], eps=ep/2, autoeps=false)) == 0

# test mesh rotation
println("Testing Mesh Rotation...")
r1 = rotate(1, 1, 0, 0, 0, 0, 1, 0, 0, pi/2)
r2 = rotate(1, 1, 0, 0, 0, 0, 0, 1, 0, pi/2)
r3 = rotate(1, 1, 0, 0, 0, 0, 0, 0, 1, pi/2)
rr1 = [1.0, 0, 1.0]
rr2 = [0.0, 1.0, -1.0]
rr3 = [-1.0, 1.0, 0.0]
for i = 1:3
    @test isapprox(r1[i] ,rr1[i])
    @test isapprox(r2[i] ,rr2[i])
    @test isapprox(r3[i] ,rr3[i])
end

rotate!(binarymesh, 45.0, [1.0,0.0,0.0], [0.0,0.0,0.0])
@test binarymesh.bounds == Bounds3(10.0,10.0,13.762255133518483,0.0,-8.509035245341185,0.0)

rotate!(binarymesh, 45.0, [0.0,1.0,0.0], [0.0,0.0,0.0])
@test binarymesh.bounds == Bounds3(16.96357128682594,10.0,13.762255133518483,0.0,-8.509035245341185,-8.509035245341185)

# Issue 21
mesh = PolygonMesh(data_path*"cube.stl")
a = MeshSlice(mesh, [eps()])
@test length(a) == 1
@test length(a[1].polygons) == 1
@test length(a[1].polygons[1].segments) == 4

# run lint
lintpkg("MeshSlicer")


text:
module MLUtilities
using StatsBase
export logsumexp, checkgrad, checkgrad_highdim
function logsumexp(x...)
  m = maximum([(x[i])[1] for i = 1:length(x)])
  return m + log(+([exp(a-m) for a in x]...))
end

function checkgrad(x::Array{Float64},func::Function,grad::Function; eps=1e-6)
  x = copy(x)
  ndim = length(x)
  f = func(x)
  g = grad(x)
  g2 = copy(g)
  for i=1:ndim
    x[i] += eps
    f2 = func(x)
    g2[i] = (f2-f)/eps
    x[i] -= eps
  end
  maximum(abs(g2-g))
end

function checkgrad_highdim(x::Array{Float64},func::Function,grad::Function; eps=1e-6)
  x = copy(x)
  ndim = length(x)
  f = func(x)
  g = grad(x)

  #randomly check 100 dimensions
  dimensions = sort(sample(1:ndim,100,replace=false))
  g2 = zeros(100)
  for i=1:100
    x[dimensions[i]] += eps
    f2 = func(x)
    g2[i] = (f2-f)/eps
    x[dimensions[i]] -= eps
  end
  maximum(abs(g2-g[dimensions]))
end

end