text:
"""
    likelihood(DM::DataModel, θ::AbstractVector) -> Real
Calculates the likelihood ``L(\\mathrm{data} \\, | \\, \\theta)`` a `DataModel` and a parameter configuration ``\\theta``.
"""
likelihood(args...; kwargs...) = exp(loglikelihood(args...; kwargs...))

## Prefix underscore for likelihood, Score and FisherMetric indicates that Prior has already been accounted for upstream
loglikelihood(DM::AbstractDataModel; kwargs...) = θ::AbstractVector{<:Number} -> loglikelihood(DM, θ; kwargs...)

# import Distributions.loglikelihood
"""
    loglikelihood(DM::DataModel, θ::AbstractVector) -> Real
Calculates the logarithm of the likelihood ``L``, i.e. ``\\ell(\\mathrm{data} \\, | \\, \\theta) \\coloneqq \\mathrm{ln} \\big( L(\\mathrm{data} \\, | \\, \\theta) \\big)`` given a `DataModel` and a parameter configuration ``\\theta``.
"""
loglikelihood(DM::AbstractDataModel, θ::AbstractVector{<:Number}, LogPriorFn::Union{Nothing, Function}=LogPrior(DM); kwargs...) = loglikelihood(Data(DM), Predictor(DM), θ, LogPriorFn; kwargs...)

loglikelihood(DS::AbstractDataSet, model::ModelOrFunction, θ::AbstractVector{<:Number}, LogPriorFn::Nothing; kwargs...) = _loglikelihood(DS, model, θ; kwargs...)
loglikelihood(DS::AbstractDataSet, model::ModelOrFunction, θ::AbstractVector{<:Number}, LogPriorFn::Function; kwargs...) = _loglikelihood(DS, model, θ; kwargs...) + EvalLogPrior(LogPriorFn, θ)

# Specialize this for different DataSet types
@inline function _loglikelihood(DS::AbstractDataSet, model::ModelOrFunction, θ::AbstractVector{<:Number}; kwargs...)
    -0.5*(DataspaceDim(DS)*log(2π) - logdetInvCov(DS) + InnerProduct(yInvCov(DS), ydata(DS)-EmbeddingMap(DS, model, θ; kwargs...)))
end

AutoScore(DM::AbstractDataModel, θ::AbstractVector{<:Number}; ADmode::Union{Symbol,Val}=Val(:ForwardDiff), kwargs...) = _AutoScore(Data(DM), Predictor(DM), θ; ADmode=ADmode, kwargs...) + EvalLogPriorGrad(LogPrior(DM), θ; ADmode=ADmode, kwargs...)
AutoMetric(DM::AbstractDataModel, θ::AbstractVector{<:Number}; ADmode::Union{Symbol,Val}=Val(:ForwardDiff), kwargs...) = _AutoMetric(Data(DM), Predictor(DM), θ; ADmode=ADmode, kwargs...) - EvalLogPriorHess(LogPrior(DM), θ; ADmode=ADmode, kwargs...)
_AutoScore(DS::AbstractDataSet, model::ModelOrFunction, θ::AbstractVector{<:Number}; ADmode::Union{Symbol,Val}=Val(:ForwardDiff), kwargs...) = GetGrad(ADmode, x->_loglikelihood(DS, model, x; kwargs...))(θ)
_AutoMetric(DS::AbstractDataSet, model::ModelOrFunction, θ::AbstractVector{<:Number}; ADmode::Union{Symbol,Val}=Val(:ForwardDiff), kwargs...) = GetHess(ADmode, x->-_loglikelihood(DS, model, x; kwargs...))(θ)

Score(DM::AbstractDataModel; kwargs...) = θ::AbstractVector{<:Number} -> Score(DM, θ; kwargs...)
NegScore(DM::AbstractDataModel; kwargs...) = θ::AbstractVector{<:Number} -> -Score(DM, θ; kwargs...)

"""
    Score(DM::DataModel, θ::AbstractVector{<:Number}; Auto::Val=Val(false))
Calculates the gradient of the log-likelihood ``\\ell`` with respect to a set of parameters ``\\theta``. `Auto=Val(true)` uses automatic differentiation.
"""
Score(DM::AbstractDataModel, θ::AbstractVector{<:Number}, LogPriorFn::Union{Nothing,Function}=LogPrior(DM); kwargs...) = Score(Data(DM), Predictor(DM), dPredictor(DM), θ, LogPriorFn; kwargs...)

Score(DS::AbstractDataSet, model::ModelOrFunction, dmodel::ModelOrFunction, θ::AbstractVector{<:Number}, LogPriorFn::Nothing; kwargs...) = _Score(DS, model, dmodel, θ; kwargs...)
Score(DS::AbstractDataSet, model::ModelOrFunction, dmodel::ModelOrFunction, θ::AbstractVector{<:Number}, LogPriorFn::Function; kwargs...) = _Score(DS, model, dmodel, θ; kwargs...) + EvalLogPriorGrad(LogPriorFn, θ)

# convert Auto from kwarg to arg
_Score(DS::AbstractDataSet, model::ModelOrFunction, dmodel::ModelOrFunction, θ::AbstractVector{<:Number}; Auto::Val=Val(false), kwargs...) = _Score(DS, model, dmodel, θ, Auto; kwargs...)
# Delegate to AutoScore
_Score(DS::AbstractDataSet, model::ModelOrFunction, dmodel::ModelOrFunction, θ::AbstractVector{<:Number}, Auto::Val{true}; kwargs...) = _AutoScore(DS, model, θ; ADmode=Auto, kwargs...)
_Score(DS::AbstractDataSet, model::ModelOrFunction, dmodel::ModelOrFunction, θ::AbstractVector{<:Number}, Auto::Union{Symbol,Val}; kwargs...) = _AutoScore(DS, model, θ; ADmode=Auto, kwargs...)

# Specialize this for different DataSet types
@inline function _Score(DS::AbstractDataSet, model::ModelOrFunction, dmodel::ModelOrFunction, θ::AbstractVector{<:Number}, Auto::Val{false}; kwargs...)
    transpose(EmbeddingMatrix(DS,dmodel,θ; kwargs...)) * (yInvCov(DS) * (ydata(DS) - EmbeddingMap(DS,model,θ; kwargs...)))
end

"""
Point θ lies outside confidence region of level `Confvol` if this function > 0.
"""
WilksCriterion(DM::AbstractDataModel, θ::AbstractVector{<:BigFloat}, Confvol::BigFloat=ConfVol(BigFloat(1.)); dof::Int=pdim(DM), kwargs...) = ChisqCDF(dof, 2(LogLikeMLE(DM) - loglikelihood(DM,θ; kwargs...))) - Confvol
WilksCriterion(DM::AbstractDataModel, θ::AbstractVector{<:Number}, Confvol::Real=ConfVol(1.); dof::Int=pdim(DM), kwargs...) = cdf(Chisq(dof), 2(LogLikeMLE(DM) - loglikelihood(DM, θ; kwargs...))) - Confvol

# Do not give default to third argument here such as to not overrule the defaults from above
WilksCriterion(DM::AbstractDataModel, θ::AbstractVector{<:Float64}, Confvol::BigFloat; kwargs...) = WilksCriterion(DM, BigFloat.(θ), Confvol; kwargs...)
WilksCriterion(DM::AbstractDataModel, θ::AbstractVector{<:BigFloat}, Confvol::Float64; kwargs...) = WilksCriterion(DM, θ, BigFloat(Confvol); kwargs...)

"""
    WilksTest(DM::DataModel, θ::AbstractVector{<:Number}, Confvol=ConfVol(1)) -> Bool
Checks whether a given parameter configuration `θ` is within a confidence interval of level `Confvol` using Wilks' theorem.
This makes the assumption, that the likelihood has the form of a normal distribution, which is asymptotically correct in the limit that the number of datapoints is infinite.
The keyword `dof` can be used to manually specify the degrees of freedom.
"""
WilksTest(DM::AbstractDataModel, θ::AbstractVector{<:Number}, Confvol::Real=ConfVol(one(suff(θ))); kwargs...)::Bool = WilksCriterion(DM, θ, Confvol; kwargs...) < 0.

function FindConfBoundary(DM::AbstractDataModel, Confnum::Real; tol::Real=4e-15, maxiter::Int=10000, dof::Int=pdim(DM))
    CF = tol < 2e-15 ? ConfVol(BigFloat(Confnum)) : ConfVol(Confnum)
    mle = if CF isa BigFloat
        suff(MLE(DM)) != BigFloat && println("FindConfBoundary: Promoting MLE to BigFloat and continuing. However, it is advisable to promote the entire DataModel object via DM = BigFloat(DM) instead.")
        BigFloat.(MLE(DM))
    else
        MLE(DM)
    end
    Test(x) = WilksTest(DM, x .* BasisVector(1, length(mle)) + mle, CF; dof=dof)
    FindConfBoundary(Test, mle; tol=tol, maxiter=maxiter)
end
function FindConfBoundary(Test::Function, mle::AbstractVector{<:Number}; tol::Real=4e-15, maxiter::Int=10000)
    res = mle + LineSearch(Test, zero(suff(mle)); tol=tol, maxiter=maxiter) .* BasisVector(1, length(mle))
    tol < 2e-15 ? res : convert(Vector{Float64}, res)
end

# function FtestPrepared(DM::DataModel, θ::AbstractVector, S_MLE::Real, ConfVol=ConfVol(1))::Bool
#     n = length(ydata(DM));  p = length(θ);    S(P) = sum(((ydata(DM) .- map(x->DM.model(x,P),xdata(DM)))./ysigma(DM)).^2)
#     S(θ) ≤ S_MLE * (1. + p/(n-p)) * quantile(FDist(p, n-p),ConfVol)
# end
# Ftest(DM::DataModel, θ::AbstractVector, MLE::AbstractVector, Conf=ConfVol(1))::Bool = FtestPrepared(DM,θ,sum((ydata(DM) .- map(x->DM.model(x,MLE),xdata(DM))).^2),Conf)

# equivalent to ResidualSquares(DM,MLE(DM))
RS_MLE(DM::AbstractDataModel) = logdetInvCov(DM) - Npoints(DM)*ydim(DM)*log(2π) - 2LogLikeMLE(DM)
ResidualSquares(DM::AbstractDataModel, θ::AbstractVector{<:Number}) = InnerProduct(yInvCov(DM), ydata(DM) - EmbeddingMap(DM,θ))
function FCriterion(DM::AbstractDataModel, θ::AbstractVector{<:Number}, Confvol::Real=ConfVol(one(suff(θ))))
    n = length(ydata(DM));  p = length(θ)
    ResidualSquares(DM,θ) - RS_MLE(DM) * (1. + length(θ)/(n - p)) * quantile(FDist(p, n-p),Confvol)
end
function FTest(DM::AbstractDataModel, θ::AbstractVector{<:Number}, Confvol::Real=ConfVol(one(suff(θ))))::Bool
    FCriterion(DM,θ,Confvol) < 0
end
function FindFBoundary(DM::AbstractDataModel, Confnum::Real; tol::Real=4e-15, maxiter::Int=10000)
    ((suff(MLE(DM)) != BigFloat) && tol < 2e-15) && throw("For tol < 2e-15, MLE(DM) must first be promoted to BigFloat via DM = DataModel(Data(DM),DM.model,DM.dmodel,BigFloat.(MLE(DM))).")
    CF = ConfVol(Confnum)
    Test(x) = FTest(DM,x .* BasisVector(1,pdim(DM)) + MLE(DM), CF)
    res = MLE(DM) .+ LineSearch(Test, zero(suff(CF)); tol=tol, maxiter=maxiter) .* BasisVector(1,pdim(DM))
    tol < 2e-15 ? res : convert(Vector{Float64}, res)
end

FDistCDF(x,d1,d2) = beta_inc(d1/2.,d2/2.,d1*x/(d1*x + d2)) #, 1 .-d1*BigFloat(x)/(d1*BigFloat(x) + d2))[1]
# function Ftest2(DM::DataModel, point::AbstractVector{T}, MLE::AbstractVector{T}, ConfVol::T=ConfVol(1))::Bool where {T<:BigFloat}
#     n = length(ydata(DM));  p = length(point);    S(P) = sum(((ydata(DM) .- map(x->DM.model(x,P),xdata(DM)))./ysigma(DM)).^2)
#     FDistCDF(S(point) / (S(MLE) * (1 + p/(n-p))),p,n-p) ≤ ConfVol
# end

inversefactor(m) = 1. / sqrt((m - 1.) + (m - 1.)^2)
GetAlpha(x::AbstractVector{<:Number}) = GetAlpha(length(x))
@inline function GetAlpha(n::Int)
    V = Vector{Float64}(undef,n)
    fill!(V,-inversefactor(n))
    V[end] = (n-1) * inversefactor(n)
    V
end

"""
    OrthVF(DM::DataModel, θ::AbstractVector{<:Number}; Auto::Val=Val(false)) -> Vector
Calculates a direction (in parameter space) in which the value of the log-likelihood does not change, given a parameter configuration ``\\theta``.
`Auto=Val(true)` uses automatic differentiation to calculate the score.
"""
function OrthVF(DM::AbstractDataModel, θ::AbstractVector{<:Number}; alpha::AbstractVector=GetAlpha(length(θ)), Auto::Val=Val(false), kwargs...)
    length(θ) < 2 && throw(ArgumentError("dim(Parameter Space) < 2  --> No orthogonal VF possible."))
    S = -Score(DM, θ; Auto=Auto, kwargs...);    P = prod(S);    VF = P ./ S
    normalize(alpha .* VF)
end

# Faster Method for Planar OrthVF
"""
    OrthVF(DM::DataModel, PL::Plane, θ::AbstractVector{<:Number}; Auto::Val=Val(false)) -> Vector
Calculates a direction (in parameter space) in which the value of the log-likelihood does not change, given a parameter configuration ``\\theta``.
Since a 2D `Plane` is specified, both the input `θ` as well as well as the output have 2 components.
`Auto=Val(true)` uses automatic differentiation to calculate the score.
"""
function OrthVF(DM::AbstractDataModel, PL::Plane, θ::AbstractVector{<:Number}, PlanarLogPrior::Union{Nothing,Function}=(isnothing(LogPrior(DM)) ? nothing : LogPrior(DM)∘PlaneCoordinates(PL)); Auto::Val=Val(false), kwargs...)
    S = transpose([PL.Vx PL.Vy]) * (-Score(Data(DM), Predictor(DM), dPredictor(DM), PlaneCoordinates(PL,θ), PlanarLogPrior; Auto=Auto, kwargs...))
    P = prod(S)
    normalize(SA[-P/S[1],P/S[2]])
end

# Method for general functions F
function OrthVF(F::Function, θ::AbstractVector{<:Number}; ADmode::Union{Symbol,Val}=Val(true), alpha::AbstractVector=GetAlpha(length(θ)))
    S = GetGrad(ADmode,F)(θ);    P = prod(S);    VF = P ./ S
    normalize(alpha .* VF)
end

function OrthVF(F::Function; ADmode::Union{Symbol,Val}=Val(true), alpha::AbstractVector=GetAlpha(length(θ)))
    OrthogonalVectorField(θ::AbstractVector; alpha::AbstractVector=alpha) = _OrthVF(GetGrad(ADmode, F), θ; alpha=alpha)
end

"""
    _OrthVF(dF::Function, θ::AbstractVector{<:Number}; alpha::AbstractVector=GetAlpha(length(θ)))
Computes OrthVF by evaluating the GRADIENT dF.
"""
_OrthVF(dF::Function, θ::AbstractVector{<:Number}; alpha::AbstractVector=GetAlpha(length(θ))) = _turn(dF(θ), alpha)
_OrthVF!(Res::AbstractVector{<:Number}, dF::Function, θ::AbstractVector{<:Number}; alpha::AbstractVector=GetAlpha(length(θ))) = _turn!(Res, dF(θ), alpha)

_turn!(Res::AbstractVector, S::AbstractVector, α::AbstractVector) = (map!(/, Res, α, S);    normalize!(Res))

# Out-of-place versions
function _turn(S::AbstractVector, α::AbstractVector)
    # normalize!(map(/, α, S))
    normalize!(map!(/, S, α, S))
    # normalize!(α ./ S)
end
_turn(S::AbstractVector, α::AbstractVector, normalize::Val{true}) = _turn(S::AbstractVector, α::AbstractVector)

function _turn(S::AbstractVector, α::AbstractVector, normalize::Val{false})
    P = prod(S);    VF = P ./ S
    alpha .* P ./ S
end

GetStartP(DM::AbstractDataModel) = GetStartP(Data(DM), Predictor(DM), pdim(DM))
GetStartP(DS::AbstractDataSet, model::Function, hint::Int=pdim(DS,model)) = GetStartP(hint)
GetStartP(DS::AbstractDataSet, model::ModelMap, hint::Int=-42) = ElaborateGetStartP(model)
GetStartP(hint::Int) = ones(hint) .+ 0.05*(rand(hint) .- 0.5)

ElaborateGetStartP(M::ModelMap; maxiters::Int=5000) = ElaborateGetStartP(M.Domain, M.InDomain; maxiters=maxiters)
function ElaborateGetStartP(C::HyperCube, InDom::Union{Nothing,Function}; maxiters::Int=5000)
    naivetry = GetStartP(length(C))
    IsInDomain(InDom, C, naivetry) ? naivetry : SobolStartP(C, InDom; maxiters=maxiters)
end
function SobolStartP(C::HyperCube, InDom::Union{Nothing,Function}; maxiters::Int=5000)
    X = rand(length(C));    i = 0
    S = Sobol.skip(SobolSeq(clamp(C.L, -1e5ones(length(C)), 1e5ones(length(C))), clamp(C.U, -1e5ones(length(C)), 1e5ones(length(C)))), rand(1:10*maxiters); exact=true)
    while i < maxiters
        Sobol.next!(S, X);    (_TestInDomain(InDom, X) && break);    i += 1
    end
    i == maxiters && throw("Unable to find point p satisfying InDomain(p) > 0 inside HyperCube within $maxiters iterations.")
    X
end

FindMLEBig(DM::AbstractDataModel,start::AbstractVector{<:Number}=MLE(DM),LogPriorFn::Union{Function,Nothing}=LogPrior(DM); kwargs...) = FindMLEBig(Data(DM), Predictor(DM), convert(Vector,start), LogPriorFn; kwargs...)
function FindMLEBig(DS::AbstractDataSet,model::ModelOrFunction,start::AbstractVector{<:Number}=GetStartP(DS,model),LogPriorFn::Union{Function,Nothing}=nothing;
                                    tol::Real=convert(BigFloat,10^(-precision(BigFloat)/30)), meth::Optim.AbstractOptimizer=BFGS(), kwargs...)
    NegEll(p::AbstractVector{<:Number}) = -loglikelihood(DS,model,p,LogPriorFn; kwargs...)
    sum(abs, xsigma(DS)) != 0.0 && @warn "Ignoring x-uncertainties in maximum likelihood estimation. Can be incorporated using the TotalLeastSquares() method."
    InformationGeometry.minimize(NegEll, BigFloat.(convert(Vector,start)), (model isa ModelMap ? model.Domain : nothing); meth=meth, tol=tol, kwargs...)
end

function FindMLE(DM::AbstractDataModel, start::AbstractVector{<:Number}=MLE(DM), LogPriorFn::Union{Function,Nothing}=LogPrior(DM); Big::Bool=false, tol::Real=1e-14, kwargs...)
    FindMLE(Data(DM), Predictor(DM), start, LogPriorFn; Big=Big, tol=tol, kwargs...)
end
function FindMLE(DS::AbstractDataSet, model::ModelOrFunction, start::AbstractVector{<:Number}=GetStartP(DS,model), LogPriorFn::Union{Function,Nothing}=nothing; Big::Bool=false, tol::Real=1e-14, kwargs...)
    (Big || tol < 2.3e-15 || suff(start) == BigFloat) && return FindMLEBig(DS, model, start, LogPriorFn)
    sum(abs, xsigma(DS)) != 0.0 && @warn "Ignoring x-uncertainties in maximum likelihood estimation. Can be incorporated using the TotalLeastSquares() method."
    # NegEll(p::AbstractVector{<:Number}) = -loglikelihood(DS,model,p)
    if isnothing(LogPriorFn) && !(DS isa GeneralizedDataSet)
        return curve_fit(DS, model, start; tol=tol).param
        # return InformationGeometry.minimize(NegEll, GetStartP(DS,model); meth=NelderMead(), tol=tol)
        # return optimize(NegEll, ones(pdim(DS,model)), BFGS(), Optim.Options(g_tol=tol), autodiff = :forward) |> Optim.minimizer
    else
        NegEll(p::AbstractVector{<:Number}) = -loglikelihood(DS,model,p,LogPriorFn; kwargs...)
        return InformationGeometry.minimize(NegEll, start, (model isa ModelMap ? model.Domain : nothing); tol=tol, kwargs...)
        # return curve_fit(DS, model, start; tol=tol).param
        # return optimize(NegEll, convert(Vector,start), BFGS(), Optim.Options(g_tol=tol), autodiff = :forward) |> Optim.minimizer
    end
end

# Slower than using curve_fit(; autodiff = :forward) but less prone to errors.
function FindMLE(DS::AbstractDataSet, model::ModelOrFunction, dmodel::ModelOrFunction, start::AbstractVector{<:Number}=GetStartP(DS,model), LogPriorFn::Union{Function,Nothing}=nothing; Big::Bool=false, tol::Real=1e-14, kwargs...)
    (Big || tol < 2.3e-15 || suff(start) == BigFloat) && return FindMLEBig(DS, model, start, LogPriorFn)
    sum(abs, xsigma(DS)) != 0.0 && @warn "Ignoring x-uncertainties in maximum likelihood estimation. Can be incorporated using the TotalLeastSquares() method."
    if isnothing(LogPriorFn) && !(DS isa GeneralizedDataSet)
        return curve_fit(DS, model, dmodel, start; tol=tol).param
        # return InformationGeometry.minimize(NegEll, GetStartP(DS,model); meth=NelderMead(), tol=tol)
        # return optimize(NegEll, ones(pdim(DS,model)), BFGS(), Optim.Options(g_tol=tol), autodiff = :forward) |> Optim.minimizer
    else
        NegEll(p::AbstractVector{<:Number}) = -loglikelihood(DS,model,p,LogPriorFn; kwargs...)
        NegScore(p::AbstractVector{<:Number}) = -Score(DS,model,dmodel,p,LogPriorFn; kwargs...)
        return InformationGeometry.minimize(NegEll, NegScore, start, (model isa ModelMap ? model.Domain : nothing); tol=tol, kwargs...)
        # return curve_fit(DS, model, start; tol=tol).param
        # return optimize(NegEll, convert(Vector,start), BFGS(), Optim.Options(g_tol=tol), autodiff = :forward) |> Optim.minimizer
    end
end

"""
    ConfidenceInterval1D(DM::AbstractDataModel, Confnum::Real=1.; tol::Real=1e-14) -> Tuple{Number,Number}
Returns the confidence interval associated with confidence level `Confnum` in the case of one-dimensional parameter spaces.
"""
function ConfidenceInterval1D(DM::AbstractDataModel, Confnum::Real=1.; tol::Real=1e-13)
    (tol < 2e-15 || Confnum > 8) && throw("ConfidenceInterval1D not programmed for BigFloat yet.")
    pdim(DM) != 1 && throw("ConfidenceInterval1D not defined for p != 1.")
    A = LogLikeMLE(DM) - (1/2)*InvChisqCDF(pdim(DM),ConfVol(Confnum))
    Func(p::Number) = loglikelihood(DM,MLE(DM) + p*BasisVector(1,pdim(DM))) - A
    D(f) = x->ForwardDiff.derivative(f,x);  NegFunc(x) = Func(-x)
    B = find_zero((Func,D(Func)),0.1,Roots.Order1(); xatol=tol)
    A = find_zero((Func,D(Func)),-B,Roots.Order1(); xatol=tol)
    rts = (MLE(DM)[1]+A, MLE(DM)[1]+B)
    rts[1] ≥ rts[2] ? throw("ConfidenceInterval1D errored...") : return rts
end

function SpatialBoundaryFunction(M::ModelMap)
    function ModelMapBoundaries(u,p,t)
        S = !IsInDomain(M, u)
        S && @warn "Curve ran into boundaries specified by ModelMap at $u."
        return S
    end
end
function SpatialBoundaryFunction(M::ModelMap, PL::Plane)
    function ModelMapBoundaries(u,p,t)
        S = !IsInDomain(M, PlaneCoordinates(PL,u))
        S && @warn "Curve ran into boundaries specified by ModelMap at $u."
        return S
    end
end

"""
    GenerateBoundary(DM::DataModel, u0::AbstractVector{<:Number}; tol::Real=1e-9, meth=Tsit5(), mfd::Bool=true) -> ODESolution
Basic method for constructing a curve lying on the confidence region associated with the initial configuration `u0`.
"""
function GenerateBoundary(DM::AbstractDataModel, u0::AbstractVector{<:Number}; tol::Real=1e-9, Boundaries::Union{Function,Nothing}=nothing,
                            meth::OrdinaryDiffEqAlgorithm=GetMethod(tol), mfd::Bool=false, Auto::Val=Val(false), kwargs...)
    u0 = !mfd ? PromoteStatic(u0, true) : u0
    LogLikeOnBoundary = loglikelihood(DM, u0)
    IntCurveODE!(du,u,p,t)  =  (du .= 0.1 * OrthVF(DM,u; Auto=Auto))
    g!(resid,u,p,t)  =  (resid[1] = LogLikeOnBoundary - loglikelihood(DM,u))
    terminatecondition(u,t,integrator) = u[2] - u0[2]
    # TerminateCondition only on upwards crossing --> supply two different affect functions, leave second free I
    CB = ContinuousCallback(terminatecondition,terminate!,nothing)
    CB = !isnothing(Boundaries) ? CallbackSet(CB, DiscreteCallback(Boundaries,terminate!)) : CB
    CB = Predictor(DM) isa ModelMap ? CallbackSet(CB, DiscreteCallback(SpatialBoundaryFunction(Predictor(DM)),terminate!)) : CB
    prob = ODEProblem(IntCurveODE!,u0,(0.,1e5))
    if mfd
        solve(prob, meth; reltol=tol, abstol=tol, callback=CallbackSet(CB,ManifoldProjection(g!)), kwargs...)
    else
        solve(prob, meth; reltol=tol, abstol=tol, callback=CB, kwargs...)
    end
end

function GenerateBoundary(DM::AbstractDataModel, PL::Plane, u0::AbstractVector{<:Number}; tol::Real=1e-9, mfd::Bool=false,
                            Boundaries::Union{Function,Nothing}=nothing, meth::OrdinaryDiffEqAlgorithm=GetMethod(tol), Auto::Val=Val(false), kwargs...)
    @assert length(u0) == 2
    u0 = !mfd ? PromoteStatic(u0, true) : u0
    PlanarLogPrior = isnothing(LogPrior(DM)) ? nothing : LogPrior(DM)∘PlaneCoordinates(PL)
    LogLikeOnBoundary = loglikelihood(DM, PlaneCoordinates(PL,u0), PlanarLogPrior)
    IntCurveODE!(du,u,p,t)  =  du .= 0.1 * OrthVF(DM, PL, u, PlanarLogPrior; Auto=Auto)
    g!(resid,u,p,t)  =  resid[1] = LogLikeOnBoundary - loglikelihood(DM, PlaneCoordinates(PL,u), PlanarLogPrior)
    terminatecondition(u,t,integrator) = u[2] - u0[2]
    CB = ContinuousCallback(terminatecondition,terminate!,nothing)
    CB = !isnothing(Boundaries) ? CallbackSet(CB, DiscreteCallback(Boundaries,terminate!)) : CB
    CB = Predictor(DM) isa ModelMap ? CallbackSet(CB, DiscreteCallback(SpatialBoundaryFunction(Predictor(DM),PL),terminate!)) : CB
    prob = ODEProblem(IntCurveODE!,u0,(0.,1e5))
    if mfd
        solve(prob, meth; reltol=tol, abstol=tol, callback=CallbackSet(CB, ManifoldProjection(g!)), kwargs...)
    else
        solve(prob, meth; reltol=tol, abstol=tol, callback=CB, kwargs...)
    end
end

"""
General function `F` with 2D domain whose Hessian should be negative-definite everywhere, i.e. negative cost function.
"""
GenerateBoundary(F::Function, u0::AbstractVector{<:Number}; ADmode::Union{Val,Symbol}=Val(:ForwardDiff), kwargs...) = GenerateBoundary(F, GetGrad(ADmode,F), u0; kwargs...)
function GenerateBoundary(F::Function, dF::Function, u0::AbstractVector{<:Number}; tol::Real=1e-9, mfd::Bool=false,
                            Boundaries::Union{Function,Nothing}=nothing, meth::OrdinaryDiffEqAlgorithm=GetMethod(tol), kwargs...)
    @assert length(u0) == 2
    !mfd && (u0 = PromoteStatic(u0, true))
    CheatingOrth!(du::AbstractVector, dF::AbstractVector) = (mul!(du, SA[0 1; -1 0.], dF);  normalize!(du); nothing)
    IntCurveODE!(du,u,p,t) = CheatingOrth!(du, dF(u))
    solve(ODEProblem(IntCurveODE!,u0,(0.,1e5)), meth; reltol=tol, abstol=tol, callback=_CallbackConstructor(F, u0, F(u0); Boundaries=Boundaries, mfd=mfd), kwargs...)
end

"""
Constructs appropriate callback kwarg for 2D cost function consisting of
- the termination condition when the integral curve has closed,
- a domain check based on `Boundaries` kwarg,
- optional manifold projection.
"""
function _CallbackConstructor(F::Function, u0::AbstractVector{<:Number}, FuncOnBoundary::Real; Boundaries::Union{Function,Nothing}=nothing, mfd::Bool=false)
    terminatecondition(u,t,integrator) = u[2] - u0[2]
    CB = ContinuousCallback(terminatecondition, terminate!, nothing)
    !isnothing(Boundaries) && (CB = CallbackSet(CB, DiscreteCallback(Boundaries,terminate!)))
    g!(resid,u,p,t) = (resid[1] = FuncOnBoundary - F(u))
    mfd && (CB = CallbackSet(ManifoldProjection(g!), CB)) # Eval projection first
    CB
end

"""
    ConfidenceRegion(DM::AbstractDataModel, Confnum::Real=1.; tol::Real=1e-9, meth=Tsit5(), mfd::Bool=false, Auto::Val=Val(false), parallel::Bool=false, Dirs::Tuple{Int,Int,Int}=(1,2,3), N::Int=30)
Computes confidence region of level `Confnum`. For `pdim(DM) > 2`, the confidence region is intersected by a family of `Plane`s in the directions specified by the keyword `Dirs`.
The `Plane`s and their embedded 2D confidence boundaries are returned as the respective first and second arguments in this case.
"""
function ConfidenceRegion(DM::AbstractDataModel, Confnum::Real=1.; tol::Real=1e-9, meth::OrdinaryDiffEqAlgorithm=GetMethod(tol), mfd::Bool=false,
                            Boundaries::Union{Function,Nothing}=nothing, Auto::Val=Val(false), parallel::Bool=false, Dirs::Tuple{Int,Int,Int}=(1,2,3), N::Int=30, dof::Int=pdim(DM), kwargs...)
    if pdim(DM) == 1
        ConfidenceInterval1D(DM, Confnum; tol=tol)
    elseif pdim(DM) == 2
        GenerateBoundary(DM, FindConfBoundary(DM, Confnum; tol=tol, dof=dof); tol=tol, Boundaries=Boundaries, meth=meth, mfd=mfd, Auto=Auto, kwargs...)
    else
        # println("ConfidenceRegion() computes solutions in the θ[1]-θ[2] plane which are separated in the θ[3] direction. For more explicit control, call MincedBoundaries() and set options manually.")
        Cube = LinearCuboid(DM, Confnum)
        Planes = IntersectCube(DM, Cube, Confnum; Dirs=Dirs, N=N, tol=tol)
        Planes, MincedBoundaries(DM, Planes, Confnum; tol=tol, Boundaries=Boundaries, Auto=Auto, meth=meth, mfd=mfd, parallel=parallel, kwargs...)
    end
end

IsStructurallyIdentifiable(DM::AbstractDataModel, sol::AbstractODESolution; kwargs...)::Bool = length(StructurallyIdentifiable(DM,sol; kwargs...)) == 0

function StructurallyIdentifiable(DM::AbstractDataModel, sol::AbstractODESolution; kwargs...)
    find_zeros(t->GeometricDensity(DM,sol(t); kwargs...), sol.t[1], sol.t[end])
end
function StructurallyIdentifiable(DM::AbstractDataModel, sols::AbstractVector{<:AbstractODESolution}; parallel::Bool=false, kwargs...)
    (parallel ? pmap : map)(x->StructurallyIdentifiable(DM,x; kwargs...), sols)
end

"""
    ConfidenceRegions(DM::DataModel, Range::AbstractVector)
Computes the boundaries of confidence regions for two-dimensional parameter spaces given a vector or range of confidence levels.
A convenient interface which extends this to higher dimensions is currently still under development.

For example,
```julia
ConfidenceRegions(DM, 1:3; tol=1e-9)
```
computes the ``1\\sigma``, ``2\\sigma`` and ``3\\sigma`` confidence regions associated with a given `DataModel` using a solver tolerance of ``10^{-9}``.

Keyword arguments:
* `IsConfVol = true` can be used to specify the desired confidence level directly in terms of a probability ``p \\in [0,1]`` instead of in units of standard deviations ``\\sigma``,
* `tol` can be used to quantify the tolerance with which the ODE which defines the confidence boundary is solved (default `tol = 1e-9`),
* `meth` can be used to specify the solver algorithm (default `meth = Tsit5()`),
* `Auto = Val(true)` can be chosen to compute the derivatives of the likelihood using automatic differentiation,
* `parallel = true` parallelizes the computations of the separate confidence regions provided each process has access to the necessary objects,
* `dof` can be used to manually specify the degrees of freedom.
"""
function ConfidenceRegions(DM::AbstractDataModel, Confnums::AbstractVector{<:Real}=1:1; IsConfVol::Bool=false,
                        tol::Real=1e-9, meth::OrdinaryDiffEqAlgorithm=GetMethod(tol), mfd::Bool=false, Auto::Val=Val(false),
                        Boundaries::Union{Function,Nothing}=nothing, tests::Bool=!(Predictor(DM) isa ModelMap), parallel::Bool=false, dof::Int=pdim(DM), kwargs...)
    det(FisherMetric(DM,MLE(DM))) < 1e-14 && throw("It appears as though the given model is not structurally identifiable.")
    Range = IsConfVol ? InvConfVol.(Confnums) : Confnums
    if pdim(DM) == 1
        return (parallel ? pmap : map)(x->ConfidenceRegion(DM,x; tol=tol, dof=dof), Range)
    elseif pdim(DM) == 2
        sols = if parallel
            @showprogress 1 "Computing boundaries... " pmap(x->ConfidenceRegion(DM, x; tol=tol, dof=dof, Boundaries=Boundaries, meth=meth, mfd=mfd, Auto=Auto, kwargs...), Range)
        else
            @showprogress 1 "Computing boundaries... " map(x->ConfidenceRegion(DM, x; tol=tol, dof=dof, Boundaries=Boundaries, meth=meth, mfd=mfd, Auto=Auto, kwargs...), Range)
        end
        if tests
            NotTerminated = map(x->(x.retcode != :Terminated), sols)
            sum(NotTerminated) != 0 && println("Solutions $((1:length(sols))[NotTerminated]) did not exit properly.")
            roots = StructurallyIdentifiable(DM, sols; parallel=parallel)
            Unidentifiables = map(x->(length(x) != 0), roots)
            for i in 1:length(roots)
                length(roots[i]) != 0 && println("Solution $i hits chart boundary at t = $(roots[i]) and should therefore be considered invalid.")
            end
        end
        return sols
    else
        throw("This functionality is still under construction. Use ConfidenceRegion() instead.")
    end
end

"""
    InterruptedConfidenceRegion(DM::AbstractDataModel, Confnum::Real; Boundaries::Union{Function,Nothing}=nothing, tol::Real=1e-9,
                                redo::Bool=true, meth::OrdinaryDiffEqAlgorithm=Tsit5(), mfd::Bool=true, Auto::Val=Val(false), kwargs...) -> ODESolution
Integrates along the level lines of the log-likelihood in the counter-clockwise direction until the model becomes either
1. structurally identifiable via `det(g) < tol`
2. the given `Boundaries(u,t,int)` method evaluates to `true`.
It then integrates from where this obstruction was met in the clockwise direction until said obstruction is hit again, resulting in a half-open confidence region.
"""
function InterruptedConfidenceRegion(DM::AbstractDataModel, Confnum::Real; Boundaries::Union{Function,Nothing}=nothing, tol::Real=1e-9,
                                redo::Bool=true, meth::OrdinaryDiffEqAlgorithm=GetMethod(tol), mfd::Bool=false, Auto::Val=Val(false), kwargs...)
    GenerateInterruptedBoundary(DM, FindConfBoundary(DM, Confnum; tol=tol); Boundaries=Boundaries, tol=tol, meth=meth, mfd=mfd, Auto=Auto, kwargs...)
end

"""
    GenerateInterruptedBoundary(DM::AbstractDataModel, u0::AbstractVector{<:Number}; Boundaries::Union{Function,Nothing}=nothing, tol::Real=1e-9,
                                redo::Bool=true, meth::OrdinaryDiffEqAlgorithm=Tsit5(), mfd::Bool=true, Auto::Val=Val(false), kwargs...) -> ODESolution
Integrates along the level lines of the log-likelihood in the counter-clockwise direction until the model becomes either
1. structurally identifiable via `det(g) < tol`
2. the given `Boundaries(u,t,int)` method evaluates to `true`.
It then integrates from where this obstruction was met in the clockwise direction until said obstruction is hit again, resulting in a half-open confidence region.
"""
function GenerateInterruptedBoundary(DM::AbstractDataModel, u0::AbstractVector{<:Number}; Boundaries::Union{Function,Nothing}=nothing, tol::Real=1e-9,
                                redo::Bool=true, meth::OrdinaryDiffEqAlgorithm=GetMethod(tol), mfd::Bool=false, Auto::Val=Val(false), kwargs...)
    u0 = !mfd ? PromoteStatic(u0, true) : u0
    LogLikeOnBoundary = loglikelihood(DM,u0)
    IntCurveODE!(du,u,p,t)  =  (du .= 0.1 .* OrthVF(DM, u; Auto=Auto))
    BackwardsIntCurveODE!(du,u,p,t)  =  (du .= -0.1 .* OrthVF(DM, u; Auto=Auto))
    g!(resid,u,p,t)  =  resid[1] = LogLikeOnBoundary - loglikelihood(DM, u)

    terminatecondition(u,t,integrator) = u[2] - u0[2]
    Singularity(u,t,integrator) = det(FisherMetric(DM, u)) - tol

    ForwardsTerminate = ContinuousCallback(terminatecondition,terminate!,nothing)

    CB = ContinuousCallback(Singularity,terminate!)
    CB = !isnothing(Boundaries) ? CallbackSet(CB, DiscreteCallback(Boundaries,terminate!)) : CB
    CB = Predictor(DM) isa ModelMap ? CallbackSet(CB, DiscreteCallback(SpatialBoundaryFunction(Predictor(DM)),terminate!)) : CB
    CB = mfd ? CallbackSet(CB, ManifoldProjection(g!)) : CB

    Forwardprob = ODEProblem(IntCurveODE!, u0, (0., 1e5))
    sol1 = solve(Forwardprob, meth; reltol=tol, abstol=tol, callback=CallbackSet(ForwardsTerminate, CB), kwargs...)

    if norm(sol1.u[end] - sol1.u[1]) < 10tol
        # closed loop, no apparent interruption or direct termination in worst case
        return sol1
    else
        Backprob = redo ? ODEProblem(BackwardsIntCurveODE!, sol1.u[end], (0., 4e5)) : ODEProblem(BackwardsIntCurveODE!, u0, (0., 4e5))
        sol2 = solve(Backprob, meth; reltol=tol, abstol=tol, callback=CB, kwargs...)
        return redo ? sol2 : [sol1, sol2]
    end
end

# Assume that sums from Fisher metric defined with first derivatives of loglikelihood pull out

FisherMetric(DM::AbstractDataModel; kwargs...) = θ::AbstractVector{<:Number} -> FisherMetric(DM, θ; kwargs...)

"""
    FisherMetric(DM::DataModel, θ::AbstractVector{<:Number})
Computes the Fisher metric ``g`` given a `DataModel` and a parameter configuration ``\\theta`` under the assumption that the likelihood ``L(\\mathrm{data} \\, | \\, \\theta)`` is a multivariate normal distribution.
```math
g_{ab}(\\theta) \\coloneqq -\\int_{\\mathcal{D}} \\mathrm{d}^m y_{\\mathrm{data}} \\, L(y_{\\mathrm{data}} \\,|\\, \\theta) \\, \\frac{\\partial^2 \\, \\mathrm{ln}(L)}{\\partial \\theta^a \\, \\partial \\theta^b} = -\\mathbb{E} \\bigg( \\frac{\\partial^2 \\, \\mathrm{ln}(L)}{\\partial \\theta^a \\, \\partial \\theta^b} \\bigg)
```
"""
FisherMetric(DM::AbstractDataModel, θ::AbstractVector{<:Number}, LogPriorFn::Union{Nothing,Function}=LogPrior(DM); kwargs...) = FisherMetric(Data(DM), dPredictor(DM), θ, LogPriorFn; kwargs...)

FisherMetric(DS::AbstractDataSet, dmodel::ModelOrFunction, θ::AbstractVector{<:Number}, LogPriorFn::Nothing; kwargs...) = _FisherMetric(DS, dmodel, θ; kwargs...)
# ADD MINUS SIGN FOR LogPrior TERM TO ACCOUNT FOR NEGATIVE SIGN IN DEFINTION OF FISHER METRIC
FisherMetric(DS::AbstractDataSet, dmodel::ModelOrFunction, θ::AbstractVector{<:Number}, LogPriorFn::Function; kwargs...) = _FisherMetric(DS, dmodel, θ; kwargs...) - EvalLogPriorHess(LogPriorFn, θ)
# FisherMetric(DS::AbstractDataSet, dmodel::ModelOrFunction, θ::AbstractVector{<:Number}, LogPriorFn::Function; kwargs...) = exp(EvalLogPrior(LogPriorFn, θ)) * (_FisherMetric(DS, dmodel, θ; kwargs...) - EvalLogPriorHess(LogPriorFn, θ))

# Specialize this for other DataSet types
_FisherMetric(DS::AbstractDataSet, dmodel::ModelOrFunction, θ::AbstractVector{<:Number}; kwargs...) = Pullback(DS, dmodel, yInvCov(DS), θ; kwargs...)

"""
    GeometricDensity(DM::AbstractDataModel, θ::AbstractVector) -> Real
Computes the square root of the determinant of the Fisher metric ``\\sqrt{\\mathrm{det}\\big(g(\\theta)\\big)}`` at the point ``\\theta``.
"""
GeometricDensity(DM::AbstractDataModel, θ::AbstractVector{<:Number}; kwargs...) = FisherMetric(DM, θ; kwargs...) |> det |> sqrt
# GeometricDensity(DS::AbstractDataSet, dmodel::ModelOrFunction, θ::AbstractVector{<:Number}; kwargs...) = FisherMetric(DS, dmodel, θ; kwargs...) |> det |> sqrt
GeometricDensity(Metric::Function, θ::AbstractVector{<:Number}; kwargs...) = sqrt(det(Metric(θ; kwargs...)))
GeometricDensity(DM::AbstractDataModel; kwargs...) = θ::AbstractVector{<:Number} -> GeometricDensity(DM, θ; kwargs...)

"""
    ConfidenceRegionVolume(DM::AbstractDataModel, Confnum::Real; N::Int=Int(1e5), WE::Bool=true, Approx::Bool=false, kwargs...) -> Real
Computes coordinate-invariant volume of confidence region associated with level `Confnum` via Monte Carlo by integrating the geometric density factor.
For likelihoods which are particularly expensive to evaluate, `Approx=true` can improve the performance by approximating the confidence region via polygons.
"""
function ConfidenceRegionVolume(DM::AbstractDataModel, Confnum::Real; N::Int=Int(1e5), WE::Bool=true, Approx::Bool=false, kwargs...)
    if Approx
        return ConfidenceRegionVolume(DM, ConfidenceRegion(DM,Confnum;tol=1e-6); N=N, WE=WE, Approx=Approx, kwargs...)
    else
        # Might not need to compute ConfidenceRegion if pdim > 2
        Domain = if pdim(DM) == 2
            # For pdim == 2, Bounding box from confidence region more performant than ProfileLikelihood
            ConstructCube(ConfidenceRegion(DM,Confnum;tol=1e-6); Padding=1e-2)
        else
            ProfileBox(DM, InterpolatedProfiles(ProfileLikelihood(DM, Confnum+2; plot=false)), Confnum; Padding=1e-2)
        end
        return IntegrateOverConfidenceRegion(DM, Domain, Confnum, z->GeometricDensity(DM,z; kwargs...); N=N, WE=WE, kwargs...)
    end
end
function ConfidenceRegionVolume(DM::AbstractDataModel, sol::AbstractODESolution; N::Int=Int(1e5), WE::Bool=true, Approx::Bool=false, kwargs...)
    @assert pdim(DM) == length(sol.u[1]) == 2
    Domain = ConstructCube(sol; Padding=1e-2)
    if Approx
        IntegrateOverApproxConfidenceRegion(DM, Domain, sol, z->GeometricDensity(DM,z;kwargs...); N=N, WE=WE)
    else
        IntegrateOverConfidenceRegion(DM, Domain, GetConfnum(DM, sol), z->GeometricDensity(DM,z;kwargs...); N=N, WE=WE, kwargs...)
    end
end
function ConfidenceRegionVolume(DM::AbstractDataModel, Tup::Tuple{<:AbstractVector{<:Plane},<:AbstractVector{<:AbstractODESolution}}; N::Int=Int(1e5), WE::Bool=true, Approx::Bool=false, kwargs...)
    ConfidenceRegionVolume(DM, Tup[1], Tup[2]; N=N, WE=WE, Approx=Approx, kwargs...)
end
function ConfidenceRegionVolume(DM::AbstractDataModel, Planes::AbstractVector{<:Plane}, sols::AbstractVector{<:AbstractODESolution}, Confnum::Real=GetConfnum(DM,Planes,sols); N::Int=Int(1e5), WE::Bool=true, Approx::Bool=false, kwargs...)
    Domain = ProfileBox(DM, InterpolatedProfiles(ProfileLikelihood(DM, Confnum+2; plot=false)), Confnum; Padding=1e-2)
    if Approx
        IntegrateOverApproxConfidenceRegion(DM, Domain, Planes, sols, z->GeometricDensity(DM,z;kwargs...); N=N, WE=WE)
    else
        IntegrateOverConfidenceRegion(DM, Domain, Confnum, z->GeometricDensity(DM,z;kwargs...); N=N, WE=WE, kwargs...)
    end
end

"""
    CoordinateVolume(DM::AbstractDataModel, Confnum::Real; N::Int=Int(1e5), WE::Bool=true, Approx::Bool=false, kwargs...) -> Real
Computes coordinate-dependent apparent volume of confidence region associated with level `Confnum` via Monte Carlo integration.
For likelihoods which are particularly expensive to evaluate, `Approx=true` can improve the performance by approximating the confidence region via polygons.
"""
function CoordinateVolume(DM::AbstractDataModel, Confnum::Real; N::Int=Int(1e5), WE::Bool=true, Approx::Bool=false, kwargs...)
    if Approx
        return CoordinateVolume(DM, ConfidenceRegion(DM,Confnum;tol=1e-6); N=N, WE=WE, Approx=Approx, kwargs...)
    else
        # Might not need to compute ConfidenceRegion if pdim > 2
        Domain = if pdim(DM) == 2
            # For pdim == 2, Bounding box from confidence region more performant than ProfileLikelihood
            ConstructCube(ConfidenceRegion(DM,Confnum;tol=1e-6); Padding=1e-2)
        else
            ProfileBox(DM, InterpolatedProfiles(ProfileLikelihood(DM, Confnum+2; plot=false)), Confnum; Padding=1e-2)
        end
        return IntegrateOverConfidenceRegion(DM, Domain, Confnum, z->one(suff(z)); N=N, WE=WE, kwargs...)
    end
end
function CoordinateVolume(DM::AbstractDataModel, sol::AbstractODESolution; N::Int=Int(1e5), WE::Bool=true, Approx::Bool=false, kwargs...)
    @assert pdim(DM) == length(sol.u[1]) == 2
    Domain = ConstructCube(sol; Padding=1e-2)
    if Approx
        IntegrateOverApproxConfidenceRegion(DM, Domain, sol, z->one(suff(z)); N=N, WE=WE)
    else
        IntegrateOverConfidenceRegion(DM, Domain, GetConfnum(DM, sol), z->one(suff(z)); N=N, WE=WE, kwargs...)
    end
end
function CoordinateVolume(DM::AbstractDataModel, Tup::Tuple{<:AbstractVector{<:Plane},<:AbstractVector{<:AbstractODESolution}}; N::Int=Int(1e5), WE::Bool=true, Approx::Bool=false, kwargs...)
    CoordinateVolume(DM, Tup[1], Tup[2]; N=N, WE=WE, Approx=Approx, kwargs...)
end
function CoordinateVolume(DM::AbstractDataModel, Planes::AbstractVector{<:Plane}, sols::AbstractVector{<:AbstractODESolution}, Confnum::Real=GetConfnum(DM,Planes,sols); N::Int=Int(1e5), WE::Bool=true, Approx::Bool=false, kwargs...)
    Domain = ProfileBox(DM, InterpolatedProfiles(ProfileLikelihood(DM, Confnum+2; plot=false)), Confnum; Padding=1e-2)
    if Approx
        IntegrateOverApproxConfidenceRegion(DM, Domain, Planes, sols, z->one(suff(z)); N=N, WE=WE)
    else
        IntegrateOverConfidenceRegion(DM, Domain, Confnum, z->one(suff(z)); N=N, WE=WE, kwargs...)
    end
end

SphereVolumeFactor(n::Int) = π^(n/2) / gamma(n/2 + 1)
ExpectedInvariantVolume(DM::AbstractDataModel, Confnum::Real) = SphereVolumeFactor(pdim(DM)) * GeodesicRadius(DM, Confnum)^pdim(DM)

GeodesicRadius(DM::AbstractDataModel, Confnum::Real) = GeodesicRadius(Confnum, pdim(DM))
GeodesicRadius(Confnum::Real, dim::Int) = sqrt(InvChisqCDF(dim, ConfVol(Confnum)))

"""
    CoordinateDistortion(DM::AbstractDataModel, Confnum::Real=1) -> Real
For CoordinateDistortions ≪ 1, the model predictions are extremely sensitive with respect to the parameters.
For CoordinateDistortion ⪎ 1, the model is comparatively insensitive towards the parameters.

This quantity is computed from the ratio of the coordinate-dependent apparent volume of a confidence region compared with the coordinate-invariant volume, which is obtained from integrating over the appropriate volume form / geometric density factor.
The unit of this quantity is ``[L^n]`` where ``L`` is the unit of length of each of the components.
"""
function CoordinateDistortion(DM::AbstractDataModel, Confnum::Real=1; Approx::Bool=false, WE::Bool=true, N::Int=Int(1e5), kwargs...)
    CoordinateVolume(DM, Confnum; N=N, Approx=Approx, WE=WE, kwargs...) / ExpectedInvariantVolume(DM, Confnum)
end

# Sensitivity independent of quality of measured datapoints (number and uncertainties), roughly independent of Confnum
function Sensitivity(DM::AbstractDataModel, Confnum::Real=1; Approx::Bool=false, WE::Bool=true, N::Int=Int(1e5), kwargs...)
    1 / CoordinateDistortion(DM, Confnum; Approx=Approx, WE=WE, N=N, kwargs...)
end

### for-loop typically slower than reduce(vcat, ...)
### Apparently curve_fit() throws an error in conjuction with ForwardDiff when reinterpret() is used
# Reduction(X::AbstractVector{<:SVector}) = reinterpret(suff(X), X)
Reduction(X::AbstractVector{<:AbstractVector}) = reduce(vcat, X)
Reduction(X::AbstractVector{<:Number}) = X

# h(θ) ∈ Dataspace
"""
    EmbeddingMap(DM::AbstractDataModel, θ::AbstractVector{<:Number}) -> Vector
Returns a vector of the collective predictions of the `model` as evaluated at the x-values and the parameter configuration ``\\theta``.
```
h(\\theta) \\coloneqq \\big(y_\\mathrm{model}(x_1;\\theta),...,y_\\mathrm{model}(x_N;\\theta)\\big) \\in \\mathcal{D}
```
"""
EmbeddingMap(DM::AbstractDataModel, θ::AbstractVector{<:Number}, woundX::AbstractVector=WoundX(DM); kwargs...) = EmbeddingMap(Data(DM), Predictor(DM), θ, woundX; kwargs...)
EmbeddingMap(DS::AbstractDataSet, model::ModelOrFunction, θ::AbstractVector{<:Number}, woundX::AbstractVector=WoundX(DS); kwargs...) = _CustomOrNot(DS, model, θ, woundX; kwargs...)
EmbeddingMap(DS::Val, model::ModelOrFunction, θ::AbstractVector{<:Number}, woundX::AbstractVector; kwargs...) = _CustomOrNot(DS, model, θ, woundX; kwargs...)

_CustomOrNot(DS::Union{Val,AbstractDataSet}, model::Function, θ::AbstractVector{<:Number}, woundX::AbstractVector; kwargs...) = _CustomOrNot(DS, model, θ, woundX, Val(false); kwargs...)
_CustomOrNot(DS::Union{Val,AbstractDataSet}, M::ModelMap, θ::AbstractVector{<:Number}, woundX::AbstractVector; kwargs...) = _CustomOrNot(DS, M.Map, θ, woundX, M.CustomEmbedding; kwargs...)

# Specialize this for different Dataset types e.g. via ::Val{CompositeDataSet}
_CustomOrNot(::Union{Val,AbstractDataSet}, model::Function, θ::AbstractVector{<:Number}, woundX::AbstractVector, custom::Val{false}; kwargs...) = Reduction(map(x->model(x,θ; kwargs...), woundX))
_CustomOrNot(::Union{Val,AbstractDataSet}, model::Function, θ::AbstractVector{<:Number}, woundX::AbstractVector, custom::Val{true}; kwargs...) = model(woundX, θ; kwargs...)

"""
    EmbeddingMatrix(DM::AbstractDataModel, θ::AbstractVector{<:Number}) -> Matrix
Returns the jacobian of the embedding map as evaluated at the x-values and the parameter configuration ``\\theta``.
"""
EmbeddingMatrix(DM::AbstractDataModel, θ::AbstractVector{<:Number}, woundX::AbstractVector=WoundX(DM); kwargs...) = EmbeddingMatrix(Data(DM), dPredictor(DM), θ, woundX; kwargs...)
EmbeddingMatrix(DS::AbstractDataSet, dmodel::ModelOrFunction, θ::AbstractVector{<:Number}, woundX::AbstractVector=WoundX(DS); kwargs...) = _CustomOrNotdM(DS, dmodel, θ, woundX; kwargs...)
EmbeddingMatrix(DS::Val, dmodel::ModelOrFunction, θ::AbstractVector{<:Number}, woundX::AbstractVector; kwargs...) = _CustomOrNotdM(DS, dmodel, θ, woundX; kwargs...)

_CustomOrNotdM(DS::Union{Val,AbstractDataSet}, dmodel::Function, θ::AbstractVector{<:Number}, woundX::AbstractVector; kwargs...) = _CustomOrNotdM(DS, dmodel, floatify(θ), woundX, Val(false); kwargs...)
_CustomOrNotdM(DS::Union{Val,AbstractDataSet}, dM::ModelMap, θ::AbstractVector{<:Number}, woundX::AbstractVector; kwargs...) = _CustomOrNotdM(DS, dM.Map, floatify(θ), woundX, dM.CustomEmbedding; kwargs...)

# Specialize this for different Dataset types e.g. via ::Val{CompositeDataSet}
_CustomOrNotdM(::Union{Val,AbstractDataSet}, dmodel::Function, θ::AbstractVector{<:Number}, woundX::AbstractVector, custom::Val{false}; kwargs...) = reduce(vcat, map(x->dmodel(x,θ; kwargs...), woundX))
_CustomOrNotdM(::Union{Val,AbstractDataSet}, dmodel::Function, θ::AbstractVector{<:Number}, woundX::AbstractVector, custom::Val{true}; kwargs...) = dmodel(woundX, θ; kwargs...)

### very slightly faster apparently
# function performDMap2(DS::AbstractDataSet,dmodel::ModelOrFunction,θ::AbstractVector{<:Number},woundX::AbstractVector)
#     Res = Array{suff(θ)}(undef,Npoints(DS)*ydim(DS),length(θ))
#     for i in 1:Npoints(DS)
#         Res[1+(i-1)*ydim(DS):(i*ydim(DS)),:] = dmodel(woundX[i],θ)
#     end;    Res
# end

# M ⟵ D
Pullback(DM::AbstractDataModel, F::Function, θ::AbstractVector{<:Number}; kwargs...) = F(EmbeddingMap(DM, θ; kwargs...))

"""
    Pullback(DM::AbstractDataModel, ω::AbstractVector{<:Number}, θ::AbstractVector) -> Vector
Pull-back of a covector to the parameter manifold ``T^*\\mathcal{M} \\longleftarrow T^*\\mathcal{D}``.
"""
Pullback(DM::AbstractDataModel, ω::AbstractVector{<:Number}, θ::AbstractVector{<:Number}; kwargs...) = transpose(EmbeddingMatrix(DM, θ; kwargs...)) * ω

"""
    Pullback(DM::DataModel, G::AbstractArray{<:Number,2}, θ::AbstractVector) -> Matrix
Pull-back of a (0,2)-tensor `G` to the parameter manifold.
"""
Pullback(DM::AbstractDataModel, G::AbstractMatrix, θ::AbstractVector{<:Number}; kwargs...) = Pullback(Data(DM), dPredictor(DM), G, θ; kwargs...)
@inline function Pullback(DS::AbstractDataSet, dmodel::ModelOrFunction, G::AbstractMatrix, θ::AbstractVector{<:Number}; kwargs...)
    J = EmbeddingMatrix(DS, dmodel, θ; kwargs...)
    transpose(J) * G * J
end

# M ⟶ D
"""
    Pushforward(DM::DataModel, X::AbstractVector, θ::AbstractVector) -> Vector
Calculates the push-forward of a vector `X` from the parameter manifold to the data space ``T\\mathcal{M} \\longrightarrow T\\mathcal{D}``.
"""
Pushforward(DM::AbstractDataModel, X::AbstractVector, θ::AbstractVector{<:Number}; kwargs...) = EmbeddingMatrix(DM, θ; kwargs...) * X

"""
    AIC(DM::DataModel, θ::AbstractVector) -> Real
Calculates the Akaike Information Criterion given a parameter configuration ``\\theta`` defined by ``\\mathrm{AIC} = 2 \\, \\mathrm{length}(\\theta) -2 \\, \\ell(\\mathrm{data} \\, | \\, \\theta)``.
Lower values for the AIC indicate that the associated model function is more likely to be correct. For linearly parametrized models and small sample sizes, it is advisable to instead use the AICc which is more accurate.
"""
AIC(DM::AbstractDataModel, θ::AbstractVector{<:Number}=MLE(DM); kwargs...) = 2length(θ) - 2loglikelihood(DM, θ; kwargs...)

"""
    AICc(DM::DataModel, θ::AbstractVector) -> Real
Computes Akaike Information Criterion with an added correction term that prevents the AIC from selecting models with too many parameters (i.e. overfitting) in the case of small sample sizes.
``\\mathrm{AICc} = \\mathrm{AIC} + \\frac{2\\mathrm{length}(\\theta)^2 + 2 \\mathrm{length}(\\theta)}{N - \\mathrm{length}(\\theta) - 1}`` where ``N`` is the number of data points.
Whereas AIC constitutes a first order estimate of the information loss, the AICc constitutes a second order estimate. However, this particular correction term assumes that the model is **linearly parametrized**.
"""
function AICc(DM::AbstractDataModel, θ::AbstractVector{<:Number}=MLE(DM); kwargs...)
    (Npoints(DM) - length(θ) - 1) == 0 && throw("DataSet too small to appy AIC correction. Use AIC instead.")
    AIC(DM, θ; kwargs...) + (2length(θ)^2 + 2length(θ)) / (Npoints(DM) - length(θ) - 1)
end

"""
    BIC(DM::DataModel, θ::AbstractVector) -> Real
Calculates the Bayesian Information Criterion given a parameter configuration ``\\theta`` defined by ``\\mathrm{BIC} = \\mathrm{ln}(N) \\cdot \\mathrm{length}(\\theta) -2 \\, \\ell(\\mathrm{data} \\, | \\, \\theta)`` where ``N`` is the number of data points.
"""
BIC(DM::AbstractDataModel, θ::AbstractVector{<:Number}=MLE(DM); kwargs...) = length(θ)*log(Npoints(DM)) - 2loglikelihood(DM, θ; kwargs...)

"""
    ModelComparison(DM1::AbstractDataModel, DM2::AbstractDataModel) -> Tuple{Int,Real}
Compares the AICc values of both models at best fit and estimates probability that one model is more likely than the other.
First entry of tuple returns which model is more likely to be correct (1 or 2) whereas the second entry returns the ratio of probabilities.
"""
function ModelComparison(DM1::AbstractDataModel, DM2::AbstractDataModel; kwargs...)
    !(ydata(DM1) == ydata(DM2) && xdata(DM1) == xdata(DM2) && yInvCov(DM1) == yInvCov(DM2)) && throw("Not comparing against same data!")
    Mod1 = AICc(DM1,MLE(DM1); kwargs...);      Mod2 = AICc(DM2,MLE(DM2); kwargs...)
    res = (Int((Mod1 > Mod2) + 1), round(exp(0.5*abs(Mod2-Mod1)),sigdigits=5))
    println("Model $(res[1]) is estimated to be $(res[2]) times as likely to be correct from difference in AICc values.")
    res
end

"""
    IsLinearParameter(DM::DataModel) -> BitVector
Checks with respect to which parameters the model function `model(x,θ)` is linear and returns vector of booleans where `true` indicates linearity.
This test is performed by comparing the Jacobians of the model for two random configurations ``\\theta_1, \\theta_2 \\in \\mathcal{M}`` column by column.
"""
function IsLinearParameter(DM::AbstractDataModel; kwargs...)
    P = pdim(DM);    J1 = EmbeddingMatrix(DM,rand(P); kwargs...);    J2 = EmbeddingMatrix(DM,rand(P); kwargs...)
    BitArray(J1[:,i] == J2[:,i]  for i in 1:size(J1,2))
end

"""
    IsLinear(DM::DataModel) -> Bool
Checks whether the `model(x,θ)` function is linear with respect to all of its parameters ``\\theta \\in \\mathcal{M}``.
A componentwise check can be attained via the method `IsLinearParameter(DM)`.
"""
IsLinear(DM::AbstractDataModel; kwargs...) = all(IsLinearParameter(DM; kwargs...))

"""
    LeastInformativeDirection(DM::DataModel,θ::AbstractVector{<:Number}=MLE(DM)) -> Vector{Float64}
Returns a vector which points in the direction in which the likelihood decreases most slowly.
"""
function LeastInformativeDirection(DM::AbstractDataModel, θ::AbstractVector{<:Number}=MLE(DM); kwargs...)
    M = eigen(FisherMetric(DM,θ; kwargs...));  i = findmin(M.values)[2]
    M.vectors[:,i] / sqrt(M.values[i])
end

"""
    FindConfBoundaryOnPlane(DM::AbstractDataModel,PL::Plane,Confnum::Real=1.; tol::Real=1e-8) -> Union{AbstractVector{<:Number},Bool}
Computes point inside the plane `PL` which lies on the boundary of a confidence region of level `Confnum`.
If such a point cannot be found (i.e. does not seem to exist), the method returns `false`.
"""
function FindConfBoundaryOnPlane(DM::AbstractDataModel, PL::Plane, Confnum::Real=1.; tol::Real=1e-8, kwargs...)
    FindConfBoundaryOnPlane(DM, PL, MLEinPlane(DM, PL; tol=tol), Confnum; kwargs...)
end
function FindConfBoundaryOnPlane(DM::AbstractDataModel, PL::Plane, mle::AbstractVector{<:Number}, Confnum::Real=1.; dof::Int=length(PL), tol::Real=1e-8, maxiter::Int=10000)
    CF = ConfVol(Confnum);      model = Predictor(DM)
    PlanarLogPrior = isnothing(LogPrior(DM)) ? nothing : LogPrior(DM)∘PlaneCoordinates(PL)
    planarmod(x,p::AbstractVector{<:Number}) = model(x, PlaneCoordinates(PL,p))
    Test(x::Number) = ChisqCDF(dof, abs(2(LogLikeMLE(DM) - loglikelihood(Data(DM), planarmod, mle + [x,0.], PlanarLogPrior)))) - CF < 0.
    !Test(0.) && return false
    [LineSearch(Test, 0.; tol=tol, maxiter=maxiter), 0.] + mle
end

DoPruning(DM, Planes::AbstractVector{<:Plane}, Confnum; kwargs...) = Prune(DM, AntiPrune(DM, Planes, Confnum), Confnum; kwargs...)

function Prune(DM::AbstractDataModel, Pls::AbstractVector{<:Plane}, Confnum::Real=1.; tol::Real=1e-8)
    CF = ConfVol(Confnum)
    Planes = copy(Pls)
    while length(Planes) > 2
        !WilksTest(DM, PlaneCoordinates(Planes[1],MLEinPlane(DM,Planes[1];tol=tol)), CF) ? popfirst!(Planes) : break
    end
    while length(Planes) > 2
        !WilksTest(DM, PlaneCoordinates(Planes[end],MLEinPlane(DM,Planes[end];tol=tol)), CF) ? pop!(Planes) : break
    end
    length(Planes) == 2 && throw("For some reason, all Planes were pruned away?!")
    return Planes
end

function AntiPrune(DM::AbstractDataModel, Pls::AbstractVector{<:Plane}, Confnum::Real=1.; tol::Real=1e-8)
    Planes = copy(Pls)
    length(Planes) < 2 && throw("Not enough Planes to infer translation direction.")
    CF = ConfVol(Confnum)
    while true
        TestPlane = Shift(Planes[2], Planes[1])
        WilksTest(DM, PlaneCoordinates(TestPlane,MLEinPlane(DM,TestPlane;tol=tol)), CF) ? pushfirst!(Planes,TestPlane) : break
    end
    while true
        TestPlane = Shift(Planes[end-1], Planes[end])
        WilksTest(DM, PlaneCoordinates(TestPlane,MLEinPlane(DM,TestPlane;tol=tol)), CF) ? push!(Planes,TestPlane) : break
    end;    Planes
end

"""
    LinearCuboid(DM::AbstractDataModel, Confnum::Real=1.; Padding::Number=1/30, N::Int=200) -> HyperCube
Returns `HyperCube` which bounds the linearized confidence region of level `Confnum` for a `DataModel`.
"""
function LinearCuboid(DM::AbstractDataModel, Confnum::Real=1.; Padding::Number=1/30, N::Int=200)
    # LinearCuboid(Symmetric(InvChisqCDF(pdim(DM),ConfVol(Confnum))*inv(FisherMetric(DM, MLE(DM)))), MLE(DM); Padding=Padding, N=N)
    BoundingBox(Symmetric(InvChisqCDF(pdim(DM),ConfVol(Confnum))*inv(FisherMetric(DM, MLE(DM)))), MLE(DM); Padding=Padding)
end

# Formula from https://tavianator.com/2014/ellipsoid_bounding_boxes.html
function BoundingBox(Σ::AbstractMatrix, μ::AbstractVector=zeros(size(Σ,1)); Padding::Number=1/30)
    @assert size(Σ,1) == size(Σ,2) == length(μ)
    E = eigen(Σ)
    offsets = [dot(E.values, E.vectors[i,:].^2) for i in 1:length(E.values)] .|> sqrt
    HyperCube(μ-offsets, μ+offsets; Padding=Padding)
end

"""
    IntersectCube(DM::AbstractDataModel,Cube::HyperCube,Confnum::Real=1.; Dirs::Tuple{Int,Int,Int}=(1,2,3), N::Int=31) -> Vector{Plane}
Returns a set of parallel 2D planes which intersect `Cube`. The planes span the directions corresponding to the basis vectors corresponding to the first two components of `Dirs`.
They are separated in the direction of the basis vector associated with the third component of `Dirs`.
The keyword `N` can be used to approximately control the number of planes which are returned.
This depends on whether more (or fewer) planes than `N` are necessary to cover the whole confidence region of level `Confnum`.
"""
function IntersectCube(DM::AbstractDataModel, Cube::HyperCube, Confnum::Real=1.; N::Int=31, Dirs::Tuple{Int,Int,Int}=(1,2,3), tol::Real=1e-8)
    (!allunique(Dirs) || !all(x->(1 ≤ x ≤ length(Cube)), Dirs)) && throw("Invalid choice of Dirs: $Dirs.")
    widths = CubeWidths(Cube)
    # PL = Plane(Center(Cube), BasisVector(Dirs[1], length(Cube)), BasisVector(Dirs[2],length(Cube)))
    PL = Plane(Center(Cube), 0.5widths[Dirs[1]]*BasisVector(Dirs[1], length(Cube)), 0.5widths[Dirs[2]]*BasisVector(Dirs[2],length(Cube)))
    IntersectRegion(DM, PL, widths[Dirs[3]] * BasisVector(Dirs[3],length(Cube)), Confnum; N=N, tol=tol)
end

"""
    IntersectRegion(DM::AbstractDataModel,PL::Plane,v::AbstractVector{<:Number},Confnum::Real=1.; N::Int=31) -> Vector{Plane}
Translates family of `N` planes which are translated approximately from `-v` to `+v` and intersect the confidence region of level `Confnum`.
If necessary, planes are removed or more planes added such that the maximal family of planes is found.
"""
function IntersectRegion(DM::AbstractDataModel, PL::Plane, v::AbstractVector{<:Number}, Confnum::Real=1.; N::Int=31, tol::Real=1e-8)
    IsOnPlane(Plane(zeros(length(v)), PL.Vx, PL.Vy),v) && throw("Translation vector v = $v lies in given Plane $PL.")
    # Planes = ParallelPlanes(PL, v, range(-0.5,0.5; length=N))
    # AntiPrune(DM, Prune(DM,Planes,Confnum;tol=tol), Confnum; tol=tol)
    DoPruning(DM, ParallelPlanes(PL, v, range(-0.5,0.5;length=N)), Confnum; tol=tol)
end

function GenerateEmbeddedBoundary(DM::AbstractDataModel, PL::Plane, Confnum::Real=1.; tol::Real=1e-8, Auto::Val=Val(false),
                                Boundaries::Union{Function,Nothing}=nothing, meth::OrdinaryDiffEqAlgorithm=GetMethod(tol), mfd::Bool=false, kwargs...)
    GenerateBoundary(DM, PL, FindConfBoundaryOnPlane(DM, PL, Confnum; tol=tol); tol=tol, Boundaries=Boundaries, meth=meth, mfd=mfd, Auto=Auto, kwargs...)
end

"""
    MincedBoundaries(DM::AbstractDataModel, Planes::AbstractVector{<:Plane}, Confnum::Real=1.; tol::Real=1e-9, Auto::Val=Val(false), meth=Tsit5(), mfd::Bool=false)
Intersects the confidence boundary of level `Confnum` with `Planes` and computes `ODESolution`s which parametrize this intersection.
"""
function MincedBoundaries(DM::AbstractDataModel, Planes::AbstractVector{<:Plane}, Confnum::Real=1.; tol::Real=1e-8, Auto::Val=Val(false),
                        Boundaries::Union{Function,Nothing}=nothing, meth::OrdinaryDiffEqAlgorithm=GetMethod(tol), mfd::Bool=false, parallel::Bool=false, kwargs...)
    Map = parallel ? pmap : map
    if parallel
        @showprogress 1 "Computing planar solutions... " pmap(X->GenerateEmbeddedBoundary(DM, X, Confnum; tol=tol, Boundaries=Boundaries, meth=meth, mfd=mfd, Auto=Auto, kwargs...), Planes)
    else
        @showprogress 1 "Computing planar solutions... " map(X->GenerateEmbeddedBoundary(DM, X, Confnum; tol=tol, Boundaries=Boundaries, meth=meth, mfd=mfd, Auto=Auto, kwargs...), Planes)
    end
end

function Thinner(S::AbstractVector{<:AbstractVector{<:Number}}; threshold::Real=0.2)
    M = S |> diff .|> norm;    b = threshold * median(M);    Res = [S[1]]
    for i in 2:length(M)
        M[i-1] + M[i] < b ? continue : push!(Res, S[i])
    end;    Res
end

CastShadow(DM::AbstractDataModel, Tup::Tuple{<:AbstractVector{<:Plane},<:AbstractVector{<:AbstractODESolution}}, args...; kwargs...) = CastShadow(DM, Tup[1], Tup[2], args...; kwargs...)
CastShadow(DM::DataModel, Planes::AbstractVector{<:Plane}, sols::AbstractVector{<:AbstractODESolution}, dirs::Tuple{<:Int,<:Int}; kwargs...) = CastShadow(DM, Planes, sols, dirs[1], dirs[2]; kwargs...)
function CastShadow(DM::DataModel, Planes::AbstractVector{<:Plane}, sols::AbstractVector{<:AbstractODESolution}, dir1::Int, dir2::Int; threshold::Real=0.2)
    @assert length(Planes) == length(sols)
    @assert pdim(DM) == length(Planes[1])
    @assert 0 < dir1 ≤ pdim(DM) && 0 < dir2 ≤ pdim(DM) && dir1 != dir2

    Project(p::AbstractVector{<:Number}, dir1::Int, dir2::Int) = SA[p[dir1], p[dir2]]

    poly = map(x->Project(PlaneCoordinates(Planes[1],x), dir1, dir2), sols[1])
    for i in 2:length(Planes)
        poly = UnionPolygons(poly, map(x->Project(PlaneCoordinates(Planes[i],x), dir1, dir2), sols[i]))
    end;    Thinner(poly; threshold=threshold)
end

function ToGeos(pointlist::AbstractVector{<:AbstractVector{<:Number}})
    @assert 2 == InformationGeometry.ConsistentElDims(pointlist)
    text = "POLYGON(("
    for point in pointlist
        text *= "$(point[1]) $(point[2]),"
    end
    text *= "$(pointlist[1][1]) $(pointlist[1][2])" * "))"
    LibGEOS.readgeom(text)
end
UnionPolygons(p1::AbstractVector{<:AbstractVector{<:Number}}, p2::AbstractVector{<:AbstractVector{<:Number}}) = LibGEOS.coordinates(UnionPolygons(ToGeos(p1), ToGeos(p2)))[1]
UnionPolygons(p1::LibGEOS.AbstractPolygon, p2::LibGEOS.AbstractPolygon) = LibGEOS.union(p1,p2)

ToAmbient(DM::AbstractDataModel, pointlist::AbstractVector{<:AbstractVector{<:Number}}, dirs::Tuple{<:Int, <:Int}) = ToAmbient(DM, pointlist, dirs[1], dirs[2])
function ToAmbient(DM::AbstractDataModel, pointlist::AbstractVector{<:AbstractVector{<:Number}}, dir1::Int, dir2::Int)
    @assert 2 == InformationGeometry.ConsistentElDims(pointlist)
    @assert 0 < dir1 ≤ pdim(DM) && 0 < dir2 ≤ pdim(DM) && dir1 != dir2
    mle = copy(MLE(DM));      mle[[dir1,dir2]] .= 0.0
    PL = Plane(mle, BasisVector(dir1, pdim(DM)), BasisVector(dir2, pdim(DM)))
    map(x->PlaneCoordinates(PL,x), pointlist)
end

function ShadowTheatre(DM::AbstractDataModel, Confnum::Real=1, dirs::Tuple{<:Int,<:Int}=(1,2); tol::Real=1e-7, N::Int=50)
    @assert (1 ≤ dirs[1] ≤ pdim(DM)) && (1 ≤ dirs[2] ≤ pdim(DM)) && dirs[1] != dirs[2] && pdim(DM) > 2
    keep = trues(pdim(DM));     keep[dirs[1]] = false;      keep[dirs[2]] = false
    translationdirs = collect(1:pdim(DM))[keep]

    Planes, sols = ConfidenceRegion(DM, Confnum; tol=tol, N=N, Dirs=(dirs[1],dirs[2],translationdirs[1]))
    list = CastShadow(DM, Planes, sols, dirs)
    if length(translationdirs) > 1
        for i in translationdirs[2:end]
            Planes, sols = ConfidenceRegion(DM, Confnum; tol=tol, N=N, Dirs=(dirs[1],dirs[2],i))
            list = UnionPolygons(list, CastShadow(DM, Planes, sols, dirs))
        end
    end
    ToAmbient(DM, list, dirs)
end

"""
    LeftOfLine(q₁::AbstractVector, q₂::AbstractVector, p::AbstractVector) -> Bool
Checks if point `p` is left of the line from `q₁` to `q₂` via `det([q₁-p  q₂-p]) > 0` for 2D points.
"""
function LeftOfLine(q₁::AbstractVector, q₂::AbstractVector, p::AbstractVector)::Bool
    # @assert length(q₁) == length(q₂) == length(p) == 2
    (q₁[1] - p[1]) * (q₂[2] - p[2]) - (q₂[1] - p[1]) * (q₁[2] - p[2]) > 0
end

# Copied from Luxor.jl
"""
    isinside(p, pointlist) -> Bool
Is a point `p` inside a polygon defined by a counterclockwise list of points.
"""
function isinside(p::AbstractVector{<:Number}, pointlist::AbstractVector{<:AbstractVector})
    @assert ConsistentElDims(pointlist) == length(p) == 2
    c = false
    @inbounds for counter in eachindex(pointlist)
        q1 = pointlist[counter]
        # if reached last point, set "next point" to first point
        if counter == length(pointlist)
            q2 = pointlist[1]
        else
            q2 = pointlist[counter + 1]
        end
        if (q1[2] < p[2]) != (q2[2] < p[2]) # crossing
            if q1[1] >= p[1]
                if q2[1] > p[1]
                    c = !c
                elseif (LeftOfLine(q1, q2, p) == (q2[2] > q1[2]))
                    c = !c
                end
            elseif q2[1] > p[1]
                if (LeftOfLine(q1, q2, p) == (q2[2] > q1[2]))
                    c = !c
                end
            end
        end
    end
    c
end

"""
    ApproxInRegion(sol::ODESolution, p::AbstractVector{<:Number}) -> Bool
Blazingly fast approximative test whether a point lies within the polygon defined by the base points of a 2D ODESolution.
Especially well-suited for hypothesis testing once a confidence boundary has been explicitly computed.
"""
ApproxInRegion(sol::AbstractODESolution, p::AbstractVector{<:Number}) = isinside(p, sol.u)

function ApproxInRegion(Planes::AbstractVector{<:Plane}, sols::AbstractVector{<:AbstractODESolution}, p::AbstractVector{<:Number})
    !(ConsistentElDims(Planes) == length(p) == 3) && throw("ApproxInRegion: Cannot determine for length(p) > 3.")      # Unclear how to do this for higher dimensions.
    @assert length(Planes) == length(sols) && all(x->length(x.u[1])==2, sols)
    # Assuming all planes parallel
    ProjectionOp = InformationGeometry.ProjectionOperator(Planes[1])
    minind = findmin([DistanceToPlane(Planes[i], p, ProjectionOp) for i in eachindex(Planes)])[2]
    ApproxInRegion(sols[minind], DecomposeWRTPlane(Planes[minind], ProjectOntoPlane(Planes[minind], p)))
end

GetConfnum(DM::AbstractDataModel, θ::AbstractVector{<:Number}; dof::Int=length(θ), kwargs...) = InvConfVol(ChisqCDF(dof, 2(LogLikeMLE(DM) - loglikelihood(DM, θ; kwargs...))))
GetConfnum(DM::AbstractDataModel, sol::AbstractODESolution; kwargs...) = GetConfnum(DM, sol.u[end]; kwargs...)
GetConfnum(DM::AbstractDataModel, PL::Plane, sol::AbstractODESolution; kwargs...) = GetConfnum(DM, PlaneCoordinates(PL, sol.u[end]); kwargs...)

function GetConfnum(DM::AbstractDataModel, Planes::AbstractVector{<:Plane}, sols::AbstractVector{<:AbstractODESolution}; kwargs...)
    Nums = [GetConfnum(DM, Planes[i], sols[i]; kwargs...) for i in eachindex(Planes)]
    mean = sum(Nums) / length(Nums)
    !all(x->abs(x-mean) < 1e-5, Nums) && @warn "High Variance in given Confnums, continuing anyway with arithmetic mean."
    return mean
end

struct ConfidenceBoundary
    sols::AbstractVector{<:AbstractODESolution}
    Confnum::Real
    MLE::AbstractVector{<:Number}
    pnames::AbstractVector{String}
end

function ConfidenceBoundary(DM::AbstractDataModel, sol::AbstractODESolution)
    @assert pdim(DM) == length(sol.u[1]) == 2
    ConfidenceBoundary([sol], GetConfnum(DM, sol), MLE(DM), pnames(DM))
end
function ConfidenceBoundary(DM::AbstractDataModel, Planes::AbstractVector{<:Plane}, sols::AbstractVector{<:AbstractODESolution})
    @assert length(Planes) == length(sols)
    ConfidenceBoundary(EmbeddedODESolution(sols, Planes), GetConfnum(DM, Planes, sols), MLE(DM), pnames(DM))
end

# function isplanar(sol::AbstractODESolution)::Bool
#     p1 = sol.u[1];      p2 = sol.u[Int(ceil(length(sol.t)/3))];     p3 = sol.u[Int(ceil(2length(sol.t)/3))]
#     PL = Plane(p1, p2-p1, Make2ndOrthogonal(p2-p1,p3-p1));    all(x->DistanceToPlane(PL,x) < 1e-12, sol.u)
# end
#
# GetPlane(CB::ConfidenceBoundary) = GetPlane(CB.sols[1], CB.MLE)
# GetPlane(DM::AbstractDataModel, sol::AbstractODESolution) = GetPlane(sol, MLE(DM))
# function GetPlane(sol::AbstractODESolution, MLE::AbstractVector{<:Number})
#     @assert isplanar(sol)
#     # Assuming that the initial point was located at [a,0,0...,0] relative to MLE
#     # return sol, MLE
#     Plane(MLE, sol.u[1] .- MLE, sol.u[end÷4] - MLE)
# end

# LinearPredictionUncertainties(DM::AbstractDataModel, F::Function, Cube::HyperCube) = LinearPredictionUncertainties(DM, F, FaceCenters(Cube))
# function LinearPredictionUncertainties(DM::AbstractDataModel, F::Function, points::AbstractVector{<:AbstractVector{<:Number}})
#     best = F(MLE(DM))
#     bounds = if typeof(best) <: AbstractVector
#         extrema(map(F, points))
#     else
#         map(extrema, eachcol(Unpack(map(F, points))))
#     end
#     best, bounds
# end


text:
# This file was generated, do not modify it. # hide
figure(figsize=(8,6))
plot(curve.parameter_values, curve.measurements)
xlabel("Number of rounds", fontsize=14)
ylabel("Cross entropy", fontsize=14)
xticks([10, 100, 250, 500], fontsize=12)
yticks(0.8:0.05:1, fontsize=12)

savefig("assets/literate/EX-crabs-xgb-curve1.svg") # hide

ext: jl

max_stars_repo_name: vdayanand/MLJTutorials

max_issues_repo_name: vdayanand/MLJTutorials

max_forks_repo_name: vdayanand/MLJTutorials


text:
using ForwardDiff
function cournot_newton(eta, c1, c2, guess, tolerance) # note that guess is a 1*2 array
    profit_1(guess) = (sum(guess))^(-1/eta)*guess[1] - 0.5*c1*guess[1]^2
    profit_2(guess) = (sum(guess))^(-1/eta)*guess[2] - 0.5*c2*guess[2]^2

    diff = 1e3
    while diff > tolerance #here diff is the distance matrix
        println("Now the guess for q1 is $(guess[1]), the guess for q2 is $(guess[2]).")
        f1 = guess -> ForwardDiff.gradient(profit_1,guess) # calculate FOC w.r.t q1
        f2 = guess -> ForwardDiff.gradient(profit_2,guess) # calculate FOC w.r.t q2
        f1_prime = guess -> ForwardDiff.hessian(profit_1,guess) # calculate SOC w.r.t q1
        f2_prime = guess -> ForwardDiff.hessian(profit_2,guess) # calculate SOC w.r.t q2
        new_guess_q1 = guess[1] - f1(guess)[1] / f1_prime(guess)[1,1]
        new_guess_q2 = guess[2] - f2(guess)[2] / f2_prime(guess)[2,2]
        diff = sqrt((guess[1]-new_guess_q1)^2 + (guess[2]-new_guess_q2)^2) # eclidean distance from old guess
        guess = [new_guess_q1 new_guess_q2] # form new guess
    end
    println("The optimal quantity for firm 1 is $(guess[1]), for firm 2 is $(guess[2])")
end


text:
export download_data, get_scedc_files, query_name, get_dict_name, query_FDSN, evaluate_done, get_seisdata

function download_data(station, source::String, dday::Date, yr::Int64, path::String, aws::AWSConfig, rootdir::String="")
    """ Stream data from IRIS and NCEDC servers to local and then to seisbasin
        Works only for IRIS and NCEDC. For data from SCEDE, refer to func get_seisdata().
    """
    try
        data = get_data("FDSN", station, src=source, s=string(dday), t = string(dday+Day(1)))
        println("Downloaded data for $station on $dday")
        # save data
        for ind in 1:3
            try
                datum = data[ind]
                println(datum.id)
                if (source == "NCEDC") && (datum !=nothing)
                    if length(datum.x)>500000 # check second since NCEDC and IRIS send different formats

                        # filepathing and directories
                        fpath, upath = joinpath(rootdir, "ncedc_waveforms/$yr/$path/$(datum.id)"),"ncedc_waveforms/$yr/$path/$(datum.id)"
                        if !isdir(dirname(fpath))
                            mkpath(dirname(fpath))
                        end

                        # write, upload, and cleanup!
                        wseis(fpath, datum)
                        s3_put(aws, "seisbasin", upath, read(fpath))
                        rm(fpath)
                    end
                else # data is from IRIS
                    println(typeof(datum), "  ", length(SeisData(datum).x[1]))
                    if (datum !=nothing) && (length(SeisData(datum).x[1])>500000)

                        # filepathing and directories
                        fpath, upath = joinpath(rootdir, "iris_waveforms/$yr/$path/$(datum.id)"), "iris_waveforms/$yr/$path/$(datum.id)"
                        if !isdir(dirname(fpath))
                            mkpath(dirname(fpath))
                        end

                        # write, upload, and cleanup!
                        wseis(fpath, SeisData(datum))
                        s3_put(aws, "seisbasin", upath, read(fpath))
                        rm(fpath)
                    end
                end
            catch
                println("Difficutly processing $station.")
            end
        end
    catch
        println("Cannot download file!")
    end
end

function get_scedc_files(dd::Date, aws::AWSConfig)
    """ Get filelist of all available SCEDC ms data from SCEDC at a given date dd.
        For now, the function only query CI network.

        dd::Date
        aws:AWSConfig
        return::Array{String, 1}
    """
    ar_filelist = map(x -> s3query(aws, dd, enddate = dd, network="CI", channel=x),["BH?", "HH?"])
    filelist_scedc_BH = ar_filelist[1]
    filelist_scedc_HH = ar_filelist[2]
    # create dictionary and overwrite HH keys with available BH data

    BH_keys = [get_dict_name(file) for file in filelist_scedc_BH]
    HH_keys = [get_dict_name(file) for file in filelist_scedc_HH]

    # Convert to dictionary
    HH_dict = Dict([(name, file) for (name, file) in zip(HH_keys, filelist_scedc_HH)])
    BH_dict = Dict([(name, file) for (name, file) in zip(BH_keys, filelist_scedc_BH)])
    filelist_dict = merge(HH_dict, BH_dict) # BH dict overwrite HH_dict. This is essentually the union

    #As Julia grammer indicates, the rightmost channel will replace the leftmost one.
    #But why do we need to overwrite BH with HH ?

    filelist_scedc = collect(values(filelist_dict)) # return values as array for download
    return filelist_scedc
end

function query_name(s::String, full::Bool)
    """ Split AWS SCEDC-PDS style key into net.cha.sta.comp style.
        s::String
            e.g. "continuous_waveforms/2019/2019_034/AZBZN__BHE___2019034.ms"
        full::Bool
            Decide whether to return station's code
        return::String
            e.g. "AZ.BZN..E"
    """
    name = convert(String, split(s,"/")[end])
    net, cha, sta, comp = name[1:2], strip(convert(String, name[3:6]),'_'),
                          strip(convert(String, name[11:12]),'_'), name[10]
    if full
        return join([net, cha, sta, comp],".")
    else
        return join([net, cha, comp],".") #comp - add back in for type
    end
end

function get_dict_name(file::String)
    """ Helper function for get_scedc_files"""
    station = convert(String, split(split(file,"/")[end],"_")[1])
    component = split(file,"/")[end][10:12]
    return string(station, "_", component)
end

function query_FDSN(file_str::String, src::String, dday::Date)
    """ Query SCEDC/NCEDC and IRIS for metadata. Try twice

        file_str::String
            e.g. "UW.*.*.BH*"
        string::String
            e.g. "IRIS", "SCEDC", "NCEDC"
        dday::Date
            e.g. Date("2021-01-01"), Date("2021-01-01T12:29:00")

        return::Array{String, 1}:
            "CI.PER..BHE"
            "CI.PER..BHN"
            "CI.PER..BHZ"
    """
    println("Beginning metadata download for $src ($file_str).")
    try
        if src in ["SCEDC","NCEDC"]
            """
                As most SCEDC/NCEDC network could be found in a relatively small region,
                we restrict the region for faster query.
            """
            if file_str[1:2] ∉ ["CJ", "MX", "NP", "NR", "PB", "RB", "SB", "TA", "UL"]
                # then the region works for network from SCEDC/NCEDC.
                return FDSNsta(file_str, src=src, s = string(dday), t = string(dday+Day(1)), reg = [30., 45., -126., -112.]).id
            else
                # not much stations are beyonds the region.
                return FDSNsta(file_str, src=src, s = string(dday), t = string(dday+Day(1))).id
            end
        else # then data is IRIS so we filter to region
             # not much can be done for general stations from IRIS...
            return FDSNsta(file_str,src=src, s = string(dday), t = string(dday+Day(1))).id
        end
    catch e
        println(e)
        println("Could not download metadata for $file_str from $src on $dday.")
        println("Second trail.")
        try
            if src in ["SCEDC","NCEDC"]
                if file_str[1:2] ∉ ["CJ", "MX", "NP", "NR", "PB", "RB", "SB", "TA", "UL"]
                    return FDSNsta(file_str, src=src, s = string(dday), t = string(dday+Day(1)), reg = [30., 45., -126., -112.]).id
                else
                    return FDSNsta(file_str, src=src, s = string(dday), t = string(dday+Day(1))).id
                end
            else
                return FDSNsta(file_str,src=src, s = string(dday), t = string(dday+Day(1))).id
            end
        catch
            println("Could not download metadata for $file_str from $src on $dday.")
            println("Error on both attempts.")
        end
    end
end

function evaluate_done(yr::Int64, path::String, aws::AWSConfig)
    try
        #iris_query = [elt["Key"] for elt in collect(s3_list_objects(aws, "seisbasin", "iris_waveforms/$yr/$path/", max_items=1000))]
        #ncedc_query = [elt["Key"] for elt in collect(s3_list_objects(aws, "seisbasin", "ncedc_waveforms/$yr/$path/", max_items=1000))]
        iris_path, ncedc_path = S3Path("s3://seisbasin/iris_waveforms/$yr/$path/", config=aws), S3Path("s3://seisbasin/ncedc_waveforms/$yr/$path/", config=aws)
        iris_query, ncedc_query = convert.(String, readdir(iris_path)), convert.(String, readdir(ncedc_path))
        println("$(length(iris_query)-1) iris waveforms and $(length(ncedc_query)-1) ncedc waveforms already downloaded for $path.")
        if ((length(iris_query) < 20) && (yr <=2004)) || ((length(iris_query)<100) && (yr >=2005)) || ((length(ncedc_query) < 5) && yr <=  2003) || ((length(ncedc_query) < 80) && (yr >= 2004))
            println("Continuing to download phase for $path...")
            return false # we need to attempt to reprocess this day!
        else
            println("Data is already in seisbasin for $path. Continuing to next day...")
            return true # day is already processed
        end
    catch e
        println(e)
        println("Error in evaluating day. Continuing to download phase")
        return false
    end
end

function get_seisdata(date::Date, aws::AWSConfig, data_sources::Array{Array{String,1},1}=[["*.*.*.HH*","NCEDC"],["*.*.*.BH*","NCEDC"],
            ["*.*.*.HH*","IRIS"],["*.*.*.BH*","IRIS"]], rootdir::String="")
    """ Downloads available seismic data from NCEDC and IRIS for given date. Uploads data to seisbasin
        Also update csv (only) info from SCEDC to seisbasin.

    """
    try
        yr = Dates.year(date)
        path = join([Dates.year(date),lpad(Dates.dayofyear(date),3,"0")],"_") # Yeilds "YEAR_JDY"
        @eval @everywhere yr, path, rootdir = $yr, $path, $rootdir
        if evaluate_done(yr, path, aws)
            return 0
        else
            ##################### Query SCEDC S3 for easy stations ########################
            filelist_scedc = get_scedc_files(date, aws)

            scedc_list = map(x -> query_name(x, true), filelist_scedc)
            scedc_brief = map(x -> query_name(x, false), filelist_scedc)
            scedc_sta = [split(elt, ".")[2] for elt in scedc_brief]

            save_sources = DataFrame(Bucket = vec(fill("SCEDC", (length(scedc_brief),1))), Source = scedc_list, Station = scedc_sta, ID = scedc_brief, Path = filelist_scedc)
            # write scedc name files immediately to bucket - donwload during correlate routine
            if !isdir(joinpath(rootdir, "scedc_path/$yr/"))
                mkpath(joinpath(rootdir, "scedc_path/$yr/"))
            end
            CSV.write(joinpath(rootdir, "scedc_path/$yr/$path.csv"), save_sources)
            s3_put(aws, "seisbasin", "scedc_path/$yr/$path.csv", read(joinpath(rootdir, "scedc_path/$yr/$path.csv")))
            println("There are $(length(save_sources.ID)) SCEDC stations for $path.")

            ####################### Query NCEDC and IRIS for extra missing stations ######################
            T_list = @elapsed list_queries = map(x -> query_FDSN(x[1],x[2], date), data_sources) # pmap doesn't work bc server ....
            for list in list_queries
                try
                    filter!(x -> x[end] in ['E','N','Z'], list)
                catch e
                    println(e)
                    println("Likely no metadata from that source")
                end
            end
            data_sources_ch = [[join([split(elt,".")[1], split(elt,".")[2], elt[end]],".") for elt in array] for array in list_queries]

            for (ind, source) in enumerate(data_sources_ch) # add new stations to dataframe
                for (j, station) in enumerate(source)
                    if station ∉ save_sources.ID # check if already contained in list (lower priority stations won't be added)
                        push!(save_sources, [data_sources[ind][2], list_queries[ind][j], split(station,".")[2], station, list_queries[ind][j]])
                    end
                end
            end

            ####################### Download NCEDC and IRIS Data ######################
            # Write data and then upload to seisbasin
            if !isdir(joinpath(rootdir, "ncedc_waveforms/$yr/$path/"))
                mkpath(joinpath(rootdir, "ncedc_waveforms/$yr/$path/"))
            end
            if !isdir(joinpath(rootdir, "iris_waveforms/$yr/$path/"))
                mkpath(joinpath(rootdir, "iris_waveforms/$yr/$path/"))
            end

            NCEDC_list = filter(x -> x.Bucket =="NCEDC", save_sources).Path
            IRIS_list = filter(x -> x.Bucket =="IRIS", save_sources).Path
            println("Now trying download for $(length(NCEDC_list)) NCEDC stations and $(length(IRIS_list)) IRIS stations.")

            IRIS_stations = unique([string(x[1:end-1],"*") for x in IRIS_list])
            NCEDC_stations = unique([string(x[1:end-1],"*") for x in NCEDC_list])
            map(x -> download_data(x, "IRIS", date, yr, path, aws), IRIS_stations)
            map(x -> download_data(x, "NCEDC", date, yr, path, aws), NCEDC_stations)
            ####################### Filter and Upload data to seisbasin ######################
            # Can use glob on ec2 - doesn't like "/" on local
            # iris_filepaths = joinpath.("iris_waveforms/$yr/$path/", readdir(joinpath(rootdir, "iris_waveforms/$yr/$path/"))) # duplicates contain "[NAME]" so we filter these
            # ncedc_filepaths = joinpath.("ncedc_waveforms/$yr/$path/", readdir(joinpath(rootdir, "ncedc_waveforms/$yr/$path/")))
            # filter!(x -> !occursin(".DS_Store", x), iris_filepaths)
            # filter!(x -> !occursin(".DS_Store", x), ncedc_filepaths)

            # map(x -> s3_put(aws, "seisbasin", x, read(joinpath(rootdir, x))), iris_filepaths);
            # map(x -> s3_put(aws, "seisbasin", x, read(joinpath(rootdir, x))), ncedc_filepaths);

            # println("Transfered $(length(iris_filepaths)+length(ncedc_filepaths)) raw seisdata objects from NCEDC and IRIS to seisbasin.")

            # # Cleanup - doesn't delete directories because AWS batch gets mad
            # rm.(joinpath.(rootdir, "iris_waveforms/$yr/$path", readdir(joinpath(rootdir, "iris_waveforms/$yr/$path"))), recursive=true)
            # rm.(joinpath.(rootdir, "ncedc_waveforms/$yr/$path", readdir(joinpath(rootdir, "ncedc_waveforms/$yr/$path"))), recursive=true)
            return 0
        end
    catch e
        println(e)
        path = join([Dates.year(date),lpad(Dates.dayofyear(date),3,"0")],"_") # Yeilds "YEAR_JDY"
        println("Unable to process data for $path. Continuing ...")
        return 1
    end
end


text:
module TestExamples3DEuler

using Test
using Trixi

include("test_trixi.jl")

# pathof(Trixi) returns /path/to/Trixi/src/Trixi.jl, dirname gives the parent directory
EXAMPLES_DIR = joinpath(pathof(Trixi) |> dirname |> dirname, "examples", "3d")

@testset "Compressible Euler" begin
  @testset "elixir_euler_source_terms.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_source_terms.jl"),
      l2   = [0.01032310150257373, 0.009728768969448439, 0.009728768969448494, 0.009728768969448388, 0.015080412597559597],
      linf = [0.034894790428615874, 0.033835365548322116, 0.033835365548322116, 0.03383536554832034, 0.06785765131417065])
  end

  @testset "elixir_euler_convergence_test_pure_fv.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_convergence_test_pure_fv.jl"),
      l2   = [0.037182410351406,  0.032062252638283974, 0.032062252638283974, 0.03206225263828395,  0.12228177813586687],
      linf = [0.0693648413632646, 0.0622101894740843,   0.06221018947408474,  0.062210189474084965, 0.24196451799555962])
  end

  @testset "elixir_euler_source_terms.jl with split_form" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_source_terms.jl"),
      l2   = [0.010323101502574773, 0.009728768969449282, 0.009728768969449344, 0.009728768969449401, 0.015080412597560888],
      linf = [0.0348947904286212, 0.03383536554832034, 0.03383536554831723, 0.03383536554831679, 0.0678576513141671],
      volume_integral=VolumeIntegralFluxDifferencing(flux_central))
  end

  @testset "elixir_euler_eoc_test.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_eoc_test.jl"),
      l2   = [0.0003637241020254405, 0.0003955570866382718, 0.0003955570866383613, 0.00039555708663834417, 0.0007811613481640202],
      linf = [0.0024000660244674066, 0.0029635410025339315, 0.0029635410025292686, 0.002963541002525938, 0.007191437359396424])
  end

  @testset "elixir_euler_mortar.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_mortar.jl"),
      l2   = [0.0019011097544691046, 0.0018289464161846331, 0.0018289464161847266, 0.0018289464161847851, 0.0033547668596639966],
      linf = [0.011918626829790169, 0.011808582902362641, 0.01180858290237552, 0.011808582902357312, 0.024648094686513744])
  end

  @testset "elixir_euler_amr.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_amr.jl"),
      l2   = [0.0038281919675174224, 0.003828191967517416, 0.003828191967517415, 0.003828191967517416, 0.0057422879512759525],
      linf = [0.07390148817126874, 0.07390148817126896, 0.07390148817126896, 0.07390148817126874, 0.110852232256903],
      tspan=(0.0, 0.1))
  end

  @testset "elixir_euler_taylor_green_vortex.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_taylor_green_vortex.jl"),
      l2   = [0.000349497725802709, 0.031333864159293845, 0.03133386415929371, 0.04378595044851934, 0.015796652303527357],
      linf = [0.0013934750216293423, 0.07242017454880123, 0.07242017454880156, 0.12796560115483002, 0.07680757651078807],
      tspan = (0.0, 0.5))
  end

  @testset "elixir_euler_shockcapturing.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_shockcapturing.jl"),
      l2   = [0.02570137197844877, 0.016179934130642552, 0.01617993413064253, 0.016172648598753545, 0.09261669328795467],
      linf = [0.3954458125573179, 0.26876916180359345, 0.26876916180359345, 0.26933123042178553, 1.3724137121660251])
  end

  @testset "elixir_euler_shockcapturing.jl with initial_condition_sedov_blast_wave" begin
    # OBS! This setup does not run longer but crashes (also the parameters do not make sense) -> only for testing the IC!
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_shockcapturing.jl"),
      l2   = [0.03627061709701234, 0.05178777768691398, 0.05178777768691404, 0.051787777686913984, 0.2308467334945545],
      linf = [0.9542995612960261, 1.2654683412445416, 1.2654683412445418, 1.2654683412445418, 12.805752164787227],
      initial_condition=initial_condition_sedov_blast_wave, cfl=0.7, alpha_max=1.0, tspan=(0.0, 0.1))
  end

  @testset "elixir_euler_shockcapturing.jl with initial_condition_sedov_self_gravity" begin
    # OBS! This setup does not run longer but crashes (also the parameters do not make sense) -> only for testing the IC!
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_shockcapturing.jl"),
      l2   = [0.048222481449477446, 0.0517877776783649, 0.051787777678364795, 0.05178777767836493, 0.230846733492145],
      linf = [0.954296651568088, 1.2654683413909855, 1.2654683413909855, 1.2654683413909857, 12.805752164822744],
      initial_condition=initial_condition_sedov_self_gravity, cfl=0.7, alpha_max=1.0, tspan=(0.0, 0.1))
  end

  @testset "elixir_euler_shockcapturing_amr.jl" begin
    # OBS! This setup does not make much practical sense. It is only added to exercise the
    # `sedov_self_gravity` AMR indicator, which in its original configuration is too expensive for
    # CI testing
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_shockcapturing_amr.jl"),
      l2   = [0.02217299067704248, 0.012771561294571411, 0.01277156129457143, 0.012770635779336643, 0.08091898488262424],
      linf = [0.4047819603427084, 0.27493532130155474, 0.2749353213015551, 0.2749304638368023, 1.4053942765487641],
      maxiters=10)
  end

  @testset "elixir_euler_density_pulse.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_density_pulse.jl"),
      l2   = [0.057196526814004715, 0.057196526814004715, 0.05719652681400473, 0.057196526814004736, 0.08579479022100575],
      linf = [0.27415246703018203, 0.2741524670301829, 0.2741524670301827, 0.27415246703018226, 0.41122870054527816])
  end

  @testset "elixir_euler_ec.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_ec.jl"),
      l2   = [0.02526341317987378, 0.016632068583699623, 0.016632068583699623, 0.01662548715216875, 0.0913477018048886],
      linf = [0.4372549540810414, 0.28613118232798984, 0.28613118232799006, 0.28796686065271876, 1.5072828647309124])
  end

  @testset "elixir_euler_ec.jl with initial_condition=initial_condition_constant" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_ec.jl"),
      l2   = [4.183721551616214e-16, 6.059779958716338e-16, 4.916596221090319e-16, 9.739943366304456e-16, 3.7485908743251566e-15],
      linf = [2.4424906541753444e-15, 3.733124920302089e-15, 4.440892098500626e-15, 5.329070518200751e-15, 2.4868995751603507e-14],
      initial_condition=initial_condition_constant)
  end

  @testset "elixir_euler_ec.jl with flux_chandrashekar" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_ec.jl"),
      l2   = [0.025265721172813106, 0.016649800693500427, 0.01664980069350042, 0.01664379306708522, 0.09137248646784184],
      linf = [0.4373399329742198, 0.28434487167605427, 0.28434487167605427, 0.28522678968890774, 1.532471676033761],
      surface_flux=flux_chandrashekar, volume_flux=flux_chandrashekar)
  end

  @testset "elixir_euler_ec.jl with flux_kennedy_gruber" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_ec.jl"),
      l2   = [0.025280033869871984, 0.016675487948639846, 0.016675487948639853, 0.016668992714991282, 0.091455613470441],
      linf = [0.43348628145015766, 0.28853549062014217, 0.28853549062014217, 0.2903943042772536, 1.5236557526482426],
      surface_flux=flux_kennedy_gruber, volume_flux=flux_kennedy_gruber)
  end

  @testset "elixir_euler_ec.jl with flux_shima_etal" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_ec.jl"),
      l2   = [0.025261716925811403, 0.016637655557848952, 0.01663765555784895, 0.01663105921013437, 0.09136239054024566],
      linf = [0.43692416928732536, 0.28622033209064734, 0.28622033209064746, 0.2881197143457632, 1.506534270303663],
      surface_flux=flux_shima_etal, volume_flux=flux_shima_etal)
  end

  @testset "elixir_euler_blob_amr.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_blob_amr.jl"),
      l2   = [0.04867856452253151, 0.2640486962336911, 0.0354927658652858, 0.03549276586528571, 1.0777274757408568],
      linf = [9.558543313792217, 49.4518309553356, 10.319859082570309, 10.319859082570487, 195.1066220797401],
      tspan = (0.0, 0.2))
  end

  @testset "elixir_euler_sedov_blast_wave.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_sedov_blast_wave.jl"),
      l2   = [0.0007127163978031706, 0.0023166296394624025, 0.002316629639462401, 0.0023166296394624038, 0.010200581509653256],
      linf = [0.06344190883105805, 0.6292607955969378, 0.6292607955969377, 0.6292607955969377, 2.397746252817731],
      maxiters=5, max_level=6)
  end
end

end # module