text:
# RecipesBase.@recipe f(DM::AbstractDataModel) = DM, MLE(DM)
RecipesBase.@recipe function f(DM::AbstractDataModel, mle::AbstractVector{<:Number}=MLE(DM), xpositions::AbstractVector{<:Number}=xdata(DM))
    (xdim(DM) != 1 && Npoints(DM) > 1) && throw("Not programmed for plotting xdim != 1 yet.")
    xguide -->              (ydim(DM) > Npoints(DM) ? "Positions" : xnames(DM)[1])
    yguide -->              (ydim(DM) ==1 ? ynames(DM)[1] : "Observations")
    @series begin
        Data(DM), xpositions
    end
    markeralpha :=      0.
    linewidth -->       2
    seriescolor :=     (ydim(DM) == 1 ? get(plotattributes, :seriescolor, :red) : (:auto))
    linestyle -->       :solid
    RSEs = ResidualStandardError(DM, mle)
    RSEs = !isnothing(RSEs) ? convert.(Float64, RSEs) : RSEs
    label -->  if ydim(DM) == 1
        # "Fit with RSE≈$(RSEs[1])"
        "Fit" * (isnothing(RSEs) ? "" : " with RSE≈$(round(RSEs[1]; sigdigits=3))")
    elseif ydim(DM) ≤ Npoints(DM)
        # reshape([ynames(DM)[i] * " Fit with RSE≈$(RSEs[i])" for i in 1:ydim(DM)], 1, ydim(DM))
        reshape([ynames(DM)[i] * " Fit"*(isnothing(RSEs) ? "" : " with RSE≈$(round(RSEs[i]; sigdigits=3))")  for i in 1:ydim(DM)], 1, ydim(DM))
    else
        reshape("Fit for $(xnames(DM)[1])=" .* string.(round.(xdata(DM); sigdigits=3)), 1, length(xdata(DM)))
    end
    # ydim(DS) > Npoints(DS) && length(xpositions) != ydim(DS)
    X = ydim(DM) ≤ Npoints(DM) ? DomainSamples(extrema(xdata(DM)); N=500) : xdata(DM)
    Y = predictedY(DM, mle, X)
    # Y = EmbeddingMap(Data(DM), Predictor(DM), mle, X)
    # Y = ydim(DM) == 1 ? Y : (ydim(DM) ≤ Npoints(DM) ? Unpack(Windup(Y, ydim(DM))) : transpose(Unpack(Windup(Y, ydim(DM)))))
    if ydim(DM) ≤ Npoints(DM)
        return X, Y
    elseif xpositions != xdata(DM)
        return xpositions, Y
    else
        return Y
    end
end
function predictedY(DM::AbstractDataModel, mle::AbstractVector{<:Number}=MLE(DM), X::AbstractVector=xdata(DM))
    predictedY(Data(DM), Predictor(DM), mle, X)
end
function predictedY(DS::AbstractDataSet, model::ModelOrFunction, mle::AbstractVector{<:Number}, X::AbstractVector=xdata(DS))
    # Ignore structure of missing values for CompositeDataSet for dense prediction curve
    Y = DS isa CompositeDataSet ? EmbeddingMap(Val(true), model, mle, X) : EmbeddingMap(DS, model, mle, X)
    ydim(DS) == 1 ? Y : (ydim(DS) ≤ Npoints(DS) ? Unpack(Windup(Y, ydim(DS))) : transpose(Unpack(Windup(Y, ydim(DS)))))
end

# xpositions for PDE Datasets
RecipesBase.@recipe function f(DS::AbstractDataSet, xpositions::AbstractVector{<:Number}=xdata(DS))
    xdim(DS) != 1 && throw("Not programmed for plotting xdim != 1 yet.")
    Σ_y = typeof(ysigma(DS)) <: AbstractVector ? ysigma(DS) : sqrt.(Diagonal(ysigma(DS)).diag)
    Σ_x = DS isa DataSetExact ? (xsigma(DS) isa AbstractVector ? xsigma(DS) : sqrt.(Diagonal(xsigma(DS)).diag)) : nothing
    line -->                (:scatter, 0.8)
    xguide -->              (ydim(DS) > Npoints(DS) ? "Positions" : xnames(DS)[1])
    yguide -->              (ydim(DS) == 1 ? ynames(DS)[1] : "Observations")
    seriescolor := :auto
    if ydim(DS) == 1
        label --> "Data"
        yerror --> Σ_y
        xerror --> Σ_x
    elseif ydim(DS) ≤ Npoints(DS)       # Series per y-component
        label --> reshape("Data: " .* ynames(DS), 1, ydim(DS))
        yerror --> Unpack(Windup(Σ_y, ydim(DS)))
        xerror --> Σ_x
    else                                # Series per x-component
        # Use of xdata instead of xpositions deliberate here!
        label --> reshape("Data for $(xnames(DS)[1])=" .* string.(round.(xdata(DS); sigdigits=3)), 1, length(xdata(DS)))
        yerror --> transpose(Unpack(Windup(Σ_y, ydim(DS))))
        # No way to incorporate errors in xpositions here....
    end
    Y = if ydim(DS) == 1
        ydata(DS)
    elseif ydim(DS) ≤ Npoints(DS)
        Unpack(Windup(ydata(DS), ydim(DS)))
    else
        Unpack(Windup(ydata(DS), ydim(DS))) |> transpose
    end
    if ydim(DS) > Npoints(DS) && length(xpositions) != ydim(DS)
        # Discard xpositions if xpositions == xdata but ydim > Npoints
        return Y
    else
        return xpositions, Y
    end
end

RecipesBase.@recipe function f(DSs::AbstractVector{<:AbstractDataSet})
    layout := length(DSs)
    for x in DSs
        @series begin x end
    end
end

RecipesBase.@recipe function f(LU::HyperCube)
    length(LU.U) != 2 && throw("Cube not Planar, cannot Plot Box.")
    rectangle(LU)[:,1], rectangle(LU)[:,2]
end

RecipesBase.@recipe function f(X::AbstractVector{<:AbstractVector{<:Number}})
    marker --> :hex
    linealpha --> 0
    markersize --> 1.8
    markerstrokewidth --> 0.5
    ToCols(Unpack(X))
end

"""
    Rsquared(DM::DataModel) -> Real
Calculates the R² value associated with the maximum likelihood estimate of a `DataModel`. It should be noted that the R² value is only a valid measure for the goodness of a fit for linear relationships.
"""
function Rsquared(DM::DataModel, MLE::AbstractVector{<:Number})
    !(xdim(DM) == ydim(DM) == 1) && return -1
    mean = sum(ydata(DM)) / length(ydata(DM))
    Stot = sum(abs2, ydata(DM) .- mean)
    Sres = sum(abs2, ydata(DM) - EmbeddingMap(DM,MLE))
    1 - Sres / Stot
end
Rsquared(DM::DataModel) = Rsquared(DM, MLE(DM))

ResidualStandardError(DM::AbstractDataModel, mle::AbstractVector{<:Number}=MLE(DM)) = ResidualStandardError(Data(DM), Predictor(DM), mle)
function ResidualStandardError(DS::AbstractDataSet, model::ModelOrFunction, MLE::AbstractVector{<:Number})
    Npoints(DS) ≤ length(MLE) && (println("Too few data points to compute RSE"); return nothing)
    ydiff = ydata(DS) - EmbeddingMap(DS, model, MLE)
    Res = map(i->sqrt(sum(abs2, view(ydiff, i:ydim(DS):length(ydiff))) / (Npoints(DS) - length(MLE))), 1:ydim(DS))
    ydim(DS) == 1 ? Res[1] : Res
end

TotalRSE(DS::AbstractDataSet, model::ModelOrFunction, MLE::AbstractVector{<:Number}) = norm(ResidualStandardError(DS, model, MLE))

FittedPlot(DM::AbstractDataModel, args...; kwargs...) = Plots.plot(DM, args...; kwargs...)

ResidualPlot(DM::AbstractDataModel; kwargs...) = ResidualPlot(Data(DM), Predictor(DM), MLE(DM); kwargs...)
function ResidualPlot(DS::DataSet, model::ModelOrFunction, mle::AbstractVector{<:Number}; kwargs...)
    Plots.plot(DataModel(DataSet(xdata(DS), ydata(DS)-EmbeddingMap(DS,model,mle), ysigma(DS), dims(DS)),
                ((x,p)->(ydim(DS) == 1 ? 0.0 : zeros(ydim(DS)))),
                (x,p)->zeros(ydim(DS), length(mle)),
                mle, _loglikelihood(DS, model, mle), true); kwargs...)
end
function ResidualPlot(DS::DataSetExact, model::ModelOrFunction, mle::AbstractVector{<:Number}; kwargs...)
    Plots.plot(DataModel(DataSetExact(xdata(DS), xsigma(DS), ydata(DS)-EmbeddingMap(DS,model,mle), ysigma(DS), dims(DS)),
                ((x,p)->(ydim(DS) == 1 ? 0.0 : zeros(ydim(DS)))),
                (x,p)->zeros(ydim(DS), length(mle)),
                mle, _loglikelihood(DS, model, mle), true); kwargs...)
end
# function ResidualPlot(DM::AbstractDataModel; kwargs...)
#     !(xdim(DS) == ydim(DS) == 1) && throw("Not programmed for plotting xdim != 1 or ydim != 1 yet.")
#     Plots.plot(DataSetExact(xdata(DM),resid,ysigma(DM));kwargs...)
#     Plots.plot!(x->0,[xdata(DM)[1],xdata(DM)[end]],label="Fit")
#     Plots.plot!(legendtitle="R² ≈ $(round(Rsquared(DM),sigdigits=3))")
# end

function _PrepareEllipsePlot(pos::AbstractVector{<:Number}, cov::AbstractMatrix{<:Number}; N::Int=100)
    @assert length(pos) == size(cov,1) == size(cov,2)
    C = cholesky(collect(cov)).U
    if length(pos) == 2
        ran = range(0, 2π; length=N)
        M = Unpack([pos + C*[cos(α),sin(α)] for α in ran])
        return view(M,:,1), view(M,:,2)
    elseif length(pos) == 3
        θran = range(0, π; length=N);  ϕran = range(0, 2π; length=N)
        # M = mapreduce(x->transpose(pos + C * x), vcat, eachrow([cos.(ϕran).*sin.(θran) sin.(ϕran).*sin.(θran) cos.(θran)]))
        M = Unpack([pos + C*[cos(ϕ)*sin(θ),sin(ϕ)*sin(θ),cos(θ)] for θ in θran for ϕ in ϕran])
        #return view(M,:,1), view(M,:,2), view(M,:,3)
        M[:,1], M[:,2], M[:,3]
    else
        throw("Cannot plot Ellipses for dim > 3.")
    end
end

function PlotEllipse(pos::AbstractVector{<:Number}, cov::AbstractMatrix{<:Number}; OverWrite::Bool=false, N::Int=100, c=:blue, kwargs...)
    @assert length(pos) == size(cov,1) == size(cov,2)
    @assert 2 ≤ length(pos) ≤ 3
    F = OverWrite ? Plots.plot : Plots.plot!
    F([pos]; marker=:hex, markersize=1.5, markeralpha=1, c=c, label="")
    M = _PrepareEllipsePlot(pos, cov; N=N)
    Plots.plot!(M...; label="", fillalpha=0.2, lw=0, seriestype=[:shape,], c=c, kwargs...)
end

function ToEllipsoidTuples(DS::AbstractDataSet)
    X, Σ = if xdim(DS) > 0
        vcat(xdata(DS), ydata(DS)), BlockMatrix(HealthyCovariance(xsigma(DS)), HealthyCovariance(ysigma(DS)))
    elseif xdim(DS) == 0 && ydim(DS) > 1
        throw("Not programmed yet.")
    else
        throw("Error.")
    end
    ToEllipsoidTuples(X, Σ, dims(DS))
end

function ToEllipsoidTuples(X::AbstractVector, Σ::AbstractMatrix)
    @assert length(X) % 2 == 0
    @info "ToEllipsoidTuples: No dims given, assuming xdim = ydim = 1"
    ToEllipsoidTuples(X, Σ, (length(X)÷2, 1, 1))
end

"""
Expects [xdata;ydata] as format for X (and Σ).
"""
function ToEllipsoidTuples(X::AbstractVector, Σ::AbstractMatrix, dims::Tuple{Int,Int,Int})
    @assert length(X) == size(Σ,1) == size(Σ,2) == dims[1]*(dims[2] + dims[3])
    Xit = Iterators.partition(1:dims[1]*dims[2], dims[2])
    Yit = Iterators.partition(dims[1]*dims[2]+1:dims[1]*(dims[2] + dims[3]), dims[3])
    function _Sub(X,Σ,tup)
        inds = vcat(tup[1],tup[2])
        view(X,inds), view(Σ, inds, inds)
    end
    [_Sub(X,Σ,tup) for tup in Iterators.zip(Xit,Yit)]
end

"""
Returns [xdata;ydata] as format for X (and Σ).
"""
function FromEllipsoidTuples(M::AbstractVector{<:Tuple{AbstractVector{<:Number}, AbstractMatrix{<:Number}}})
    @info "FromEllipsoidTuples: No dims given, assuming xdim = 1, ydim = $(length(M[1][1])-1)"
    FromEllipsoidTuples(M, (length(M), 1, length(M[1][1])-1))
end

function FromEllipsoidTuples(M::AbstractVector{<:Tuple{AbstractVector{<:Number}, AbstractMatrix{<:Number}}}, dims::Tuple{Int,Int,Int})
    @assert length(M) == dims[1] && ConsistentElDims(getindex.(M,1)) == dims[2] + dims[3]
    X = vcat(mapreduce(x->x[1][1:dims[2]], vcat, M), mapreduce(x->x[1][1+dims[2]:end], vcat, M))
    Σ = Matrix{suff(M)}(undef,dims[1]*(dims[2] + dims[3]), dims[1]*(dims[2] + dims[3]))
    for (i,tup) in enumerate(Iterators.zip(Iterators.partition(1:dims[1]*dims[2], dims[2]),Iterators.partition(dims[1]*dims[2]+1:dims[1]*(dims[2] + dims[3]), dims[3])))
        inds = vcat(tup[1],tup[2])
        view(Σ,inds,inds) .= M[i][2]
    end;    X, HealthyCovariance(Σ)
end

PlotEllipses(GDS::GeneralizedDataSet; kwargs...) = PlotEllipses(dist(GDS), dims(GDS); kwargs...)
PlotEllipses(dist::ContinuousMultivariateDistribution, dims::Tuple{Int,Int,Int}; kwargs...) = PlotEllipses(GetMean(dist), Sigma(dist), dims; kwargs...)
function PlotEllipses(X::AbstractVector, Σ::AbstractMatrix, dims::Tuple{Int,Int,Int}; OverWrite::Bool=true, c=:blue, kwargs...)
    @assert length(X) == size(Σ,1) == size(Σ,2) == dims[1]*(dims[2] + dims[3])
    @assert 2 ≤ dims[2] + dims[3] ≤ 3
    p = []
    OverWrite && Plots.plot()
    Ellipses = ToEllipsoidTuples(X,Σ,dims)
    for (x, σ) in Ellipses
        p = PlotEllipse(x, σ; OverWrite=false, c=c, kwargs...)
    end;    display(p)
end

RecipesBase.@recipe function f(GDS::GeneralizedDataSet, xpositions::AbstractVector{<:Number}=xdata(GDS))
    @assert xdim(GDS) + ydim(GDS) == 2
    names = vcat(xnames(GDS),ynames(GDS))
    xguide --> names[1]
    yguide --> names[2]
    N = xdim(GDS) + ydim(GDS) == 2 ? 200 : 80
    for (x,σ) in ToEllipsoidTuples(GetMean(dist(GDS)),Sigma(dist(GDS)),dims(GDS))
        A = cholesky(collect(σ)).U * [cos.(range(0, 2π; length=N))'; sin.(range(0, 2π; length=N))']
        @series begin
            seriesalpha --> 0.3
            label := ""
            fillcolor --> :blue
            Plots.Shape(x[1] .+ A[1,:], x[2] .+ A[2,:])
        end
    end
end

meshgrid(x, y) = (repeat(x, outer=length(y)), repeat(y, inner=length(x)))

function PlotScalar(F::Function, PlanarCube::HyperCube; N::Int=100, Save::Bool=false, parallel::Bool=false, OverWrite::Bool=true, kwargs...)
    length(PlanarCube) != 2 && throw(ArgumentError("Cube not Planar."))
    Lims = PlanarCube;    A = range(Lims.L[1], Lims.U[1], length=N);    B = range(Lims.L[2], Lims.U[2], length=N)
    func(args...) = F([args...])
    Map = parallel ? pmap : map
    if Save
        X,Y = meshgrid(A,B)
        Z = Map(func,X,Y)
        p = OverWrite ? contour(X,Y,Z; fill=true, nlevels=40, kwargs...) : contour!(X,Y,Z; fill=true, nlevels=40, kwargs...)
        display(p)
        return [X Y Z]
    else
        p = OverWrite ? contour(A,B,func; fill=true, nlevels=40, kwargs...) : contour(A,B,func; fill=true, nlevels=40, kwargs...)
        return p
    end
end

function PlotScalar(F::Function, PlotPlane::Plane, PlanarCube::HyperCube, N::Int=100; Save::Bool=true, parallel::Bool=false, kwargs...)
    length(PlanarCube) != 2 && throw(ArgumentError("Cube not Planar."))
    Lims = PlanarCube;    A = range(Lims.L[1], Lims.U[1], length=N);    B = range(Lims.L[2], Lims.U[2], length=N)
    Lcomp(x,y) = F(PlaneCoordinates(PlotPlane,[x,y]))
    Map = parallel ? pmap : map
    if Save
        X,Y = meshgrid(A,B)
        Z = Map(Lcomp,X,Y)
        p = contour(X,Y,Z; fill=true, leg=false, nlevels=40, title="Plot centered around: $(round.(PlotPlane.stütz,sigdigits=4))",
            xlabel="$(PlotPlane.Vx) direction", ylabel="$(PlotPlane.Vy) direction", kwargs...)
        p = scatter!([0],[0]; lab="Center", marker=:hex)
        display(p)
        return [X Y Z]
    else
        p = contour(A,B,Lcomp; fill=true, leg=false, nlevels=40, title="Plot centered around: $(round.(PlotPlane.stütz,sigdigits=4))",
            xlabel="$(PlotPlane.Vx) direction", ylabel="$(PlotPlane.Vy) direction", kwargs...)
        p = scatter!([0],[0]; lab="Center", marker=:hex)
        return p
    end
end

############################# Plotting
function PlotLoglikelihood(DM::DataModel, MLE::AbstractVector, PlanarCube::HyperCube, N::Int=100; Save::Bool=true, parallel::Bool=false)
    length(MLE) !=2 && throw(ArgumentError("Only 2D supported."))
    length(PlanarCube) != 2 && throw(ArgumentError("Cube not Planar."))
    Lcomp(args...) = loglikelihood(DM,[args...])
    Lims = TranslateCube(PlanarCube,MLE)
    A = range(Lims.L[1], Lims.U[1], length=N);  B = range(Lims.L[2], Lims.U[2], length=N)
    Map = parallel ? pmap : map
    if Save
        X,Y = meshgrid(A,B)
        Z = Map(Lcomp,X,Y)
        p = contour(X,Y,Z, fill=true, leg=false, nlevels=40)
        p = scatter!([MLE[1]],[MLE[2]], lab="MLE: [$(round(MLE[1],sigdigits=4)),$(round(MLE[2],sigdigits=4))]", marker=:hex)
        display(p)
        return [X Y Z]
    else
        p = contour(A,B,Lcomp, fill=true, leg=false, nlevels=40)
        p = scatter!([MLE[1]],[MLE[2]], lab="MLE: [$(round(MLE[1],sigdigits=4)),$(round(MLE[2],sigdigits=4))]", marker=:hex)
        return p
    end
end

PlotLoglikelihood(DM::DataModel,MLE::AbstractVector,size::Float64=0.5,N::Int=100) = PlotLoglikelihood(DM, MLE, HyperCube([[-size,size],[-size,size]]), N)

function PlotLoglikelihood(DM::DataModel, PlotPlane::Plane, PlanarCube::HyperCube, N::Int=100; Save::Bool=true, parallel::Bool=false)
    Lcomp(x,y) = loglikelihood(DM,PlaneCoordinates(PlotPlane,[x,y]))
    length(PlanarCube) != 2 && throw(ArgumentError("Cube not Planar."))
    Lims = PlanarCube
    A = range(Lims.L[1], Lims.U[1], length=N)
    B = range(Lims.L[2], Lims.U[2], length=N)
    Map = parallel ? pmap : map
    if Save
        X,Y = meshgrid(A,B)
        Z = Map(Lcomp,X,Y)
        p = contour(X,Y,Z, fill=true, leg=false, nlevels=40, title="Plot centered around: $(round.(PlotPlane.stütz,sigdigits=4))",
            xlabel="$(PlotPlane.Vx) direction", ylabel="$(PlotPlane.Vy) direction")
        p = scatter!([0],[0],lab="Center", marker=:hex)
        display(p)
        return [X Y Z]
    else
        p = contour(A,B,Lcomp, fill=true, leg=false, nlevels=40, title="Plot centered around: $(round.(PlotPlane.stütz,sigdigits=4))",
            xlabel="$(PlotPlane.Vx) direction", ylabel="$(PlotPlane.Vy) direction")
        p = scatter!([0],[0],lab="Center", marker=:hex)
        return p
    end
end

PlotLoglikelihood(DM::DataModel, PlotPlane::Plane, size::Float64=0.5,N::Int=100;Save::Bool=true) = PlotLoglikelihood(DM, PlotPlane, HyperCube([[-size,size],[-size,size]]), N, Save=Save)

function ConstructBox(fit::LsqFit.LsqFitResult, Confnum::Real; AxisCS::Bool=true)
    E = Confnum * stderror(fit)
    HyperCube(fit.param - E, fit.param + E)
end

# Choose a Plane?
VisualizeMC(Test::Function, Boundaries::AbstractVector,N::Int=2000) = VisualizeMC(Test, HyperCube(Boundaries), N)

function VisualizeMC(Test::Function, PlanarCube::HyperCube, N::Int=2000)
    length(PlanarCube) != 2 && throw(ArgumentError("Cube not Planar."))
    Lims = PlanarCube
    YesPoints = Vector{Vector{Float64}}(undef,0)
    NoPoints = Vector{Vector{Float64}}(undef,0)
    for i in 1:N
        num = rand.(Uniform.(Lims.L,Lims.U))
        if Test(num)
            push!(YesPoints,num)
        else
            push!(NoPoints,num)
        end
    end
    Plots.plot!(rectangle(Lims),lw=2,lab="Sample Space")
    p = scatter!(Unpack(YesPoints), marker=(0.7,:hex,:green), markersize=2,lab="Inside")
    p = scatter!(Unpack(NoPoints), marker=(0.7,:circle,:red), markersize=2,lab="Outside")
    p
end

function VisualizeMC(Test::Function, sol::AbstractODESolution, N::Int=2000)
    Cube = ConstructCube(sol)
    lowers, uppers = Cube.L, Cube.U
    YesPoints = Vector{Vector{Float64}}(undef,0)
    NoPoints = Vector{Vector{Float64}}(undef,0)
    for i in 1:N
        num = rand.(Uniform.(lowers,uppers))
        if Test(num)
            push!(YesPoints,num)
        else
            push!(NoPoints,num)
        end
    end
    p = Plots.plot(sol,vars=(1,2), lw=2,xlims=(lowers[1],uppers[1]), ylims=(lowers[2],uppers[2]))
        try
        Box = rectangle(lowers...,uppers...)
        Plots.plot!(Box[:,1],Box[:,2],lw=2,lab="Sample Space")
    catch x
        println("Could not plot sampling rectangle. Fix Error:")
        println(x)
    end
    p = scatter!(Unpack(YesPoints), marker=(0.5,:hex,:green), markersize=2,lab="Inside")
    p = scatter!(Unpack(NoPoints), marker=(0.5,:circle,:red), markersize=2,lab="Outside")
    p
end

rectangle(ax,ay,bx,by) = [ax ay; bx ay; bx by; ax by; ax ay]
function rectangle(LU::HyperCube)
    length(LU.L) != 2 && throw(ArgumentError("Cube not Planar."))
    rectangle((LU.L)...,(LU.U)...)
end

"""
    VFRescale(ZeilenVecs::Array{<:Number,2},C::HyperCube;scaling=0.85)
Rescale vector to look good in 2D plot.
"""
function VFRescale(ZeilenVecs::Array{<:Number,2},C::HyperCube;scaling=0.85)
    VecsPerLine = sqrt(size(ZeilenVecs)[1])
    SpacePerVec = scaling/VecsPerLine .* CubeWidths(C)
    for i in 1:size(ZeilenVecs)[1]
        if SpacePerVec[1] < abs(ZeilenVecs[i,1])
            rescale = SpacePerVec[1] / abs(ZeilenVecs[i,1])
            ZeilenVecs[i,:] .*= rescale
        end
        if SpacePerVec[2] < abs(ZeilenVecs[i,2])
            rescale = SpacePerVec[2] / abs(ZeilenVecs[i,2])
            ZeilenVecs[i,:] .*= rescale
        end
    end
    ZeilenVecs[:,1],ZeilenVecs[:,2]
end

# function Plot2DVF(DM::DataModel,V::Function,MLE::AbstractVector,PlanarCube::HyperCube,N::Int=25;scaling::Float64=0.85, OverWrite::Bool=false)
#     length(MLE) !=2 && throw(ArgumentError("Only 2D supported."))
#     length(PlanarCube) != 2 && throw(ArgumentError("Cube not Planar."))
#     Lims = TranslateCube(PlanarCube,MLE)
#     AV, BV  = meshgrid(range(Lims.L[1], Lims.U[1], length=N), range(Lims.L[2], Lims.U[2], length=N))
#     Vcomp(a,b) = V([a,b])
#     u,v = VFRescale(Unpack(Vcomp.(AV,BV)), PlanarCube; scaling=scaling)
#     # u,v = Vcomp.(AV,BV) |> Unpack
#     # u,v = VFRescale([u v],PlanarCube,scaling=scaling)
#     if OverWrite
#         quiver(AV,BV,quiver=(u,v)) |> display
#     else
#         quiver!(AV,BV,quiver=(u,v)) |> display
#     end
#     [AV BV u v]
# end
#
# function Plot2DVF(DM::DataModel,V::Function,MLE::AbstractVector,size::Float64=0.5,N::Int=25;scaling::Float64=0.85, OverWrite::Bool=false)
#     Plot2DVF(DM,V, MLE, HyperCube([[-size,size],[-size,size]]), N; scaling=scaling, OverWrite=OverWrite)
# end

# function Plot2DVF(DM::DataModel, V::Function, PlotPlane::Plane, PlanarCube::HyperCube, N::Int=25; scaling::Float64=0.85, OverWrite::Bool=false)
#     length(PlanarCube) != 2 && throw(ArgumentError("Cube not Planar."))
#     Lims = PlanarCube
#     AV, BV  = meshgrid(range(Lims.L[1], Lims.U[1], length=N), range(Lims.L[2], Lims.U[2], length=N))
#     Vcomp(a,b) = V(PlaneCoordinates(PlotPlane,[a,b]))
#     u,v = Vcomp.(AV,BV) |> Unpack
#     u,v = VFRescale([u v],PlanarCube,scaling=scaling)
#     if OverWrite
#         quiver(AV,BV,quiver=(u,v)) |> display
#     else
#         quiver!(AV,BV,quiver=(u,v)) |> display
#     end
#     [AV BV u v]
# end

# function Plot2DVF(DM::DataModel,V::Function, PlotPlane::Plane,size::Float64=0.5, N::Int=25; scaling::Float64=0.85, OverWrite::Bool=false)
#     Plot2DVF(DM,V, PlotPlane, HyperCube([[-size,size],[-size,size]]), N; scaling=scaling, OverWrite=OverWrite)
# end

function Plot2DVF(V::Function, Lims::HyperCube; N::Int=25, scaling::Float64=0.85, OverWrite::Bool=false)
    @assert length(Lims) == length(V(Center(Lims))) == 2
    AV, BV = meshgrid(range(Lims.L[1], Lims.U[1]; length=N), range(Lims.L[2], Lims.U[2]; length=N))
    Vcomp(a,b) = V([a,b])
    u, v = VFRescale(Unpack(Vcomp.(AV,BV)), Lims; scaling=scaling)
    if OverWrite
        quiver(AV, BV, quiver=(u,v)) |> display
    else
        quiver!(AV, BV, quiver=(u,v)) |> display
    end
    [AV BV u v]
end

"""
    Deplanarize(PL::Plane,sol::AbstractODESolution; N::Int=500) -> Matrix
    Deplanarize(PL::Plane,sol::AbstractODESolution, Ts::AbstractVector{<:Number}) -> Matrix
Converts the 2D outputs of `sol` from planar coordinates associated with `PL` to the coordinates of the ambient space of `PL`.
"""
Deplanarize(PL::Plane,sol::AbstractODESolution; N::Int=500) = Deplanarize(PL,sol,range(sol.t[1],sol.t[end]; length=N))
Deplanarize(PL::Plane,sol::AbstractODESolution,Ts::AbstractVector{<:Number}) = map(t->PlaneCoordinates(PL,sol(t)),Ts) |> Unpack

"""
    VisualizeSols(sols::AbstractVector{<:AbstractODESolution}; OverWrite::Bool=true)
Visualizes vectors of type `ODESolution` using the `Plots.jl` package. If `OverWrite=false`, the solution is displayed on top of the previous plot object.
"""
function VisualizeSols(sols::AbstractVector{<:AbstractODESolution}; vars::Tuple=Tuple(1:length(sols[1].u[1])), OverWrite::Bool=true, leg::Bool=false, kwargs...)
    p = [];     OverWrite && Plots.plot()
    for sol in sols
        p = VisualizeSols(sol; vars=vars, leg=leg, kwargs...)
    end;    p
end
function VisualizeSols(sol::AbstractODESolution; vars::Tuple=Tuple(1:length(sol.u[1])), leg::Bool=false, OverWrite::Bool=false,
                                        ModelMapMeta::Union{ModelMap,Bool}=false, kwargs...)
    OverWrite && Plots.plot()
    if ModelMapMeta isa ModelMap
        names = pnames(ModelMapMeta)
        if length(vars) == 2
            return Plots.plot!(sol; xlabel=names[vars[1]], ylabel=names[vars[2]], vars=vars, leg=leg, kwargs...)
        elseif length(vars) == 3
            return Plots.plot!(sol; xlabel=names[vars[1]], ylabel=names[vars[2]], zlabel=names[vars[3]], vars=vars, leg=leg, kwargs...)
        end
        # What if vars > 3? Does Plots.jl throw an error?
    end
    Plots.plot!(sol; vars=vars, leg=leg, kwargs...)
end

function VisualizeSols(PL::Plane, sol::AbstractODESolution; vars::Tuple=Tuple(1:length(PL)), leg::Bool=false, N::Int=500,
                            ModelMapMeta::Union{ModelMap,Bool}=false, kwargs...)
    H = Deplanarize(PL, sol; N=N)
    if ModelMapMeta isa ModelMap
        names = pnames(ModelMapMeta)
        return Plots.plot!(H[:,vars[1]], H[:,vars[2]], H[:,vars[3]]; xlabel=names[vars[1]], ylabel=names[vars[2]], zlabel=names[vars[3]], leg=leg, kwargs...)
    else
        return Plots.plot!(H[:,vars[1]], H[:,vars[2]], H[:,vars[3]]; leg=leg, kwargs...)
    end
end
function VisualizeSols(PL::Plane, sols::AbstractVector{<:AbstractODESolution}; vars::Tuple=Tuple(1:length(PL)), N::Int=500, OverWrite::Bool=true, leg::Bool=false, kwargs...)
    p = [];     OverWrite && Plots.plot()
    for sol in sols
        p = VisualizeSols(PL, sol; N=N, vars=vars, leg=leg, kwargs...)
    end;    p
end

VisualizeSols(X::Tuple, args...; kwargs...) = VisualizeSols(X..., args...; kwargs...)
function VisualizeSols(PL::AbstractVector{<:Plane},sols::AbstractVector{<:AbstractODESolution}; vars::Tuple=Tuple(1:length(PL[1])), N::Int=500,
                OverWrite::Bool=true,leg::Bool=false, color=rand([:red,:blue,:green,:orange,:grey]), kwargs...)
    length(PL) != length(sols) && throw("VisualizeSols: Must receive same number of Planes and Solutions.")
    p = [];     OverWrite && Plots.plot()
    for i in 1:length(sols)
        p = VisualizeSols(PL[i], sols[i]; N=N, vars=vars, leg=leg, color=color, kwargs...)
    end;    p
end
function VisualizeSols(DM::AbstractDataModel, args...; OverWrite::Bool=true, kwargs...)
    OverWrite ? scatter([MLE(DM)]; label="MLE") : scatter!([MLE(DM)]; label="MLE")
    if Predictor(DM) isa ModelMap
        return VisualizeSols(args...; OverWrite=false, ModelMapMeta=Predictor(DM), kwargs...)
    else
        return VisualizeSols(args...; OverWrite=false, kwargs...)
    end
end

function VisualizeSols(CB::ConfidenceBoundary; vars::Tuple=Tuple(1:length(CB.MLE)), OverWrite::Bool=true, color=rand([:red,:blue,:green,:orange,:grey]), kwargs...)
    p = OverWrite ? scatter([CB.MLE]; label="MLE") : []
    if length(vars) == 2
        VisualizeSols(CB.sols[1]; vars=vars, color=color, OverWrite=false, label="$(round(CB.Confnum, sigdigits=3))σ Conf. Boundary",
                        xlabel=CB.pnames[vars[1]], ylabel=CB.pnames[vars[2]], leg=true, kwargs...)
    else
        VisualizeSols(CB.sols[1]; vars=vars, color=color, OverWrite=false, label="$(round(CB.Confnum, sigdigits=3))σ Conf. Boundary",
                        xlabel=CB.pnames[vars[1]], ylabel=CB.pnames[vars[2]], zlabel=CB.pnames[vars[3]],leg=true, kwargs...)
    end
    for sol in CB.sols[2:end]
        p = VisualizeSols(sol; vars=vars, color=color, OverWrite=false, label="", leg=true, kwargs...)
    end; p
end

function VisualizeSols(CBs::AbstractVector{<:ConfidenceBoundary}; vars::Tuple=Tuple(1:length(CBs[1].MLE)), OverWrite::Bool=true, kwargs...)
    @assert all(x->x.MLE==CBs[1].MLE, CBs)
    @assert allunique(map(x->x.Confnum,CBs))
    p = OverWrite ? Plots.scatter([CBs[1].MLE]; label="MLE") : []
    for CB in CBs
        VisualizeSols(CB; vars=vars, OverWrite=false, kwargs...)
    end; p
end

VisualizeGeos(sol::AbstractODESolution; kwargs...) = VisualizeGeos([sol]; kwargs...)
function VisualizeGeos(sols::AbstractVector{<:AbstractODESolution}; OverWrite::Bool=false, leg::Bool=false, kwargs...)
    VisualizeSols(sols; vars=Tuple(1:Int(length(sols[1].u[1])/2)), OverWrite=OverWrite, leg=leg, kwargs...)
end

function VisualizeSolPoints(sol::AbstractODESolution; kwargs...)
    Plots.plot!(sol.u; marker=:hex, line=:dash, linealpha=1, markersize=2, kwargs...)
    VisualizeSols(sol; OverWrite=false, kwargs...)
end
function VisualizeSolPoints(sols::AbstractVector{<:AbstractODESolution}; OverWrite::Bool=false, kwargs...)
    p = [];     OverWrite && Plots.plot()
    for sol in sols
        p = VisualizeSolPoints(sol; kwargs...)
    end;    p
end

XCube(DS::AbstractDataSet; Padding::Number=0.) = ConstructCube(Unpack(WoundX(DS)); Padding=Padding)
XCube(DM::AbstractDataModel; Padding::Number=0.) = XCube(Data(DM); Padding=Padding)
Grid(Cube::HyperCube, N::Int=5) = [range(Cube.L[i], Cube.U[i]; length=N) for i in 1:length(Cube)]

# function GetExtrema(DM::AbstractDataModel,sols::Union{AbstractODESolution,AbstractVector{<:AbstractODESolution}},X::Union{<:Number,AbstractVector{<:Number}}; N::Int=200)
#     low = Inf;   up = -Inf
#     for sol in sols
#         Y = map(Z->DM.model(X,sol(Z)), range(sol.t[1],sol.t[end]; length=N))
#         templow = minimum(Y);   tempup = maximum(Y)
#         if templow < low    low = templow       end
#         if up < tempup      up = tempup         end
#     end;    (low, up)
# end
# function GetExtrema(DM::AbstractDataModel,PL::Plane,sols::Union{AbstractODESolution,AbstractVector{<:AbstractODESolution}},X::Union{<:Number,AbstractVector{<:Number}}; N::Int=200)
#     low = Inf;   up = -Inf
#     for sol in sols
#         templow, tempup = map(t->DM.model(X,PlaneCoordinates(PL,sol(t))), range(sol.t[1],sol.t[end]; length=N)) |> extrema
#         if templow < low    low = templow       end
#         if up < tempup      up = tempup         end
#     end;    (low, up)
# end
# function GetExtrema(DM::AbstractDataModel,PL::AbstractVector{<:Plane},sols::AbstractVector{<:AbstractODESolution},X::Union{<:Number,AbstractVector{<:Number}}; N::Int=200)
#     length(PL) != length(sols) && throw("Dimensional Mismatch.")
#     low = Inf;   up = -Inf
#     for i in 1:length(sols)
#         templow, tempup = GetExtrema(DM,PL[i],sols[i],X; N=N)
#         if templow < low    low = templow       end
#         if up < tempup      up = tempup         end
#     end;    (low, up)
# end
#
# """
#     ConfidenceBands(DM::DataModel,sol::AbstractODESolution,domain::HyperCube; N::Int=200)
# Given a confidence interval `sol`, the pointwise confidence band around the model prediction is computed for x values in `domain` by evaluating the model on the boundary of the confidence region.
# """
# function ConfidenceBands(DM::AbstractDataModel,sols::Union{AbstractODESolution,AbstractVector{<:AbstractODESolution}},domain::HyperCube=XCube(DM); N::Int=200, Np::Int=200)
#     !(length(domain) == xdim(DM)) && throw("Dimensionality of domain inconsistent with xdim.")
#     low = Vector{Float64}(undef,N^xdim(DM));     up = Vector{Float64}(undef,N^xdim(DM))
#     X = xdim(DM) == 1 ? range(domain.L[1],domain.U[1]; length=N) : Iterators.product(Grid(domain,N)...)
#     for i in 1:length(X)
#         low[i], up[i] = GetExtrema(DM,sols,X[i]; N=Np)
#     end
#     col = rand([:red,:blue,:green,:orange,:grey])
#     Plots.plot!(X,low,color=col,label="Lower Conf. Band")
#     Plots.plot!(X,up,color=col,label="Upper Conf. Band") |> display
#     return [Unpack(collect(X)) low up]
# end
#
# function ConfidenceBands(DM::AbstractDataModel,Planes::Union{Plane,AbstractVector{<:Plane}},sols::Union{AbstractODESolution,AbstractVector{<:AbstractODESolution}},
#                                             domain::HyperCube=XCube(DM); N::Int=200, Np::Int=300)
#     length(domain) != xdim(DM) && throw("Dimensionality of domain inconsistent with xdim.")
#     !(length(Planes) == 1 || length(Planes) == length(sols)) && throw("Number of Planes inconsisten with number of ODESolutions.")
#     low = Vector{Float64}(undef,N^xdim(DM));     up = Vector{Float64}(undef,N^xdim(DM))
#     X = xdim(DM) == 1 ? range(domain.L[1],domain.U[1]; length=N) : Iterators.product(Grid(domain,N)...)
#     for i in 1:length(X)
#         low[i], up[i] = GetExtrema(DM,Planes,sols,X[i]; N=Np)
#     end
#     col = rand([:red,:blue,:green,:orange,:grey])
#     Plots.plot!(X,low,color=col,label="Lower Conf. Band")
#     Plots.plot!(X,up,color=col,label="Upper Conf. Band") |> display
#     return [Unpack(collect(X)) low up]
# end

function PlotConfidenceBands(DM::AbstractDataModel, M::AbstractMatrix{<:Number}, xpositions::Union{AbstractVector{<:Number},Nothing}=nothing;
                                        Confnum::Real=-1)
    lab = 0 < Confnum ? "$(Confnum)σ " : ""
    if size(M,2) == 3
        return Plots.plot!(view(M,:,1), view(M,:,2:3); label=[lab*"Conf. Band" ""], color=rand([:red,:blue,:green,:orange,:grey])) |> display
    else # Assume the FittedPlot splits every y-component into a separate series of points and have same number of rows as x-values
        @assert xdim(DM) == 1 && size(M,2) == 1 + 2ydim(DM)
        if xpositions isa Nothing
            for i in 1:(size(M,1)-1)
                Plots.plot!([M[i,2:2:end] M[i,3:2:end]]; color=rand([:red,:blue,:green,:orange,:grey]), label=["" ""])
            end
            Plots.plot!([M[end,2:2:end] M[end,3:2:end]]; color=rand([:red,:blue,:green,:orange,:grey]), label=["" lab*"Conf. Band"]) |> display
        elseif length(xpositions) == (size(M,2)-1) / 2
            for i in 1:(size(M,1)-1)
                Plots.plot!(xpositions, [M[i,2:2:end] M[i,3:2:end]]; color=rand([:red,:blue,:green,:orange,:grey]), label=["" ""])
            end
            Plots.plot!(xpositions, [M[end,2:2:end] M[end,3:2:end]]; color=rand([:red,:blue,:green,:orange,:grey]), label=["" lab*"Conf. Band"]) |> display
        else
            throw("Vector of xpositions wrong length.")
        end
    end
end

function ConfidenceBands(DM::AbstractDataModel, Confnum::Real, Xdomain::HyperCube=XCube(DM); N::Int=300, plot::Bool=true, samples::Int=200)
    ConfidenceBands(DM, ConfidenceRegion(DM,Confnum), Xdomain; N=N, plot=plot, samples=samples)
end

"""
    ConfidenceBands(DM::DataModel, sol::AbstractODESolution, Xdomain::HyperCube; N::Int=300, plot::Bool=true) -> Matrix
Given a confidence interval `sol`, the pointwise confidence band around the model prediction is computed for x values in `Xdomain`
by evaluating the model on the boundary of the confidence region.
"""
function ConfidenceBands(DM::AbstractDataModel, sols::Union{AbstractODESolution,AbstractVector{<:AbstractODESolution}}, Xdomain::HyperCube=XCube(DM);
                            N::Int=300, plot::Bool=true, samples::Int=max(2*length(sols),100))
    ConfidenceBands(DM, sols, DomainSamples(Xdomain; N=N); plot=plot, samples=samples)
end

function ConfidenceBands(DM::AbstractDataModel, Tup::Tuple{<:AbstractVector{<:Plane},AbstractVector{<:AbstractODESolution}}, woundX=XCube(DM); N::Int=300, plot::Bool=true, samples::Int=max(2*length(Tup[1]),100))
    ConfidenceBands(DM, Tup[1], Tup[2], woundX; plot=plot, samples=samples)
end
function ConfidenceBands(DM::AbstractDataModel, Planes::AbstractVector{<:Plane}, sols::AbstractVector{<:AbstractODESolution}, Xdomain::HyperCube=XCube(DM); N::Int=300, plot::Bool=true, samples::Int=max(2*length(sols),100))
    ConfidenceBands(DM, Planes, sols, DomainSamples(Xdomain; N=N); plot=plot, samples=samples)
end

function ConfidenceBands(DM::AbstractDataModel, sols::Union{AbstractODESolution,AbstractVector{<:AbstractODESolution}}, woundX::AbstractVector{<:Number}; plot::Bool=true, samples::Int=max(2*length(sols),100))
    Res = Array{Float64,2}(undef, length(woundX), 2*ydim(DM))
    for col in 1:2:2ydim(DM)    fill!(view(Res,:,col), Inf);     fill!(view(Res,:,col+1), -Inf)    end
    # gradually refine Res for each solution to avoid having to allocate a huge list of points
    for sol in sols
        _ConfidenceBands!(Res, DM, map(sol, range(sol.t[1], sol.t[end]; length=samples)), woundX)
    end
    M = hcat(woundX, Res)
    if plot
        Confnum = round(GetConfnum(DM, sols[1]); sigdigits=2)
        PlotConfidenceBands(DM, M; Confnum=Confnum)
    end
    return M
end

function ConfidenceBands(DM::AbstractDataModel, Planes::AbstractVector{<:Plane}, sols::AbstractVector{<:AbstractODESolution}, woundX::AbstractVector{<:Number}; plot::Bool=true, samples::Int=max(2*length(sols),100))
    @assert length(Planes) == length(sols)
    Res = Array{Float64,2}(undef, length(woundX), 2*ydim(DM))
    for col in 1:2:2ydim(DM)    fill!(view(Res,:,col), Inf);     fill!(view(Res,:,col+1), -Inf)    end
    # gradually refine Res for each solution to avoid having to allocate a huge list of points
    for i in 1:length(sols)
        _ConfidenceBands!(Res, DM, map(x->PlaneCoordinates(Planes[i], sols[i](x)), range(sols[i].t[1], sols[i].t[end]; length=samples)), woundX)
    end
    M = hcat(woundX, Res)
    if plot
        Confnum = round(GetConfnum(DM, PlaneCoordinates(Planes[1], sols[1].u[1])); sigdigits=2)
        PlotConfidenceBands(DM, M; Confnum=Confnum)
    end
    return M
end

"""
ApproxConfidenceBands(DM::AbstractDataModel, Confnum::Real, Xdomain=XCube(DM); N::Int=300, plot::Bool=true, add::Real=1.5)
"""
function ApproxConfidenceBands(DM::AbstractDataModel, Confnum::Real, Xdomain=XCube(DM); N::Int=300, plot::Bool=true, add::Real=1.5)
    higherConfnum = Confnum + add
    @warn "Trying to establish HyperCube which circumscribes confidence region by passing $higherConfnum to ProfileLikelihood(). If this errors, establish a suitable box manually."
    Box = ProfileBox(DM, InterpolatedProfiles(ProfileLikelihood(DM, higherConfnum; plot=false)), Confnum)
    ApproxConfidenceBands(DM, Box, Xdomain; N=N, plot=plot)
end

"""
    ApproxConfidenceBands(DM::AbstractDataModel, ParameterCube::HyperCube, Xdomain=XCube(DM); N::Int=300, plot::Bool=true)
Computes confidence bands associated with the face centers of the `ParameterCube`.
If the `ParameterCube` circumscribes a given confidence region, this will typically result in a gross and asymmetric overestimation of the true pointwise confidence bands associated with this confidence level.
"""
function ApproxConfidenceBands(DM::AbstractDataModel, ParameterCube::HyperCube, Xdomain::HyperCube=XCube(DM); N::Int=300, plot::Bool=true)
    length(Xdomain) != xdim(DM) && throw("Dimensionality of domain inconsistent with xdim.")
    ApproxConfidenceBands(DM, ParameterCube, DomainSamples(Xdomain; N=N); plot=plot)
end

function ApproxConfidenceBands(DM::AbstractDataModel, ParameterCube::HyperCube, woundX::AbstractVector{<:Number}; plot::Bool=true)
    ConfidenceBands(DM, FaceCenters(ParameterCube), woundX; plot=plot)
end

# Devise version with woundX::AbstractVector{<:AbstractVector{<:Number}} for xdim > 1
function ConfidenceBands(DM::AbstractDataModel, points::AbstractVector{<:AbstractVector{<:Number}}, woundX::AbstractVector{<:Number}; plot::Bool=true)
    xdim(DM) != 1 && throw("Not programmed for xdim != 1 yet.")
    Res = Array{Float64,2}(undef, length(woundX), 2*ydim(DM))
    for col in 1:2:2ydim(DM)    fill!(view(Res,:,col), Inf);     fill!(view(Res,:,col+1), -Inf)    end
    _ConfidenceBands!(Res, DM, points, woundX)
    M = hcat(woundX, Res)
    if plot
        Confnum = round(GetConfnum(DM, points[1]); sigdigits=2)
        PlotConfidenceBands(DM, M; Confnum=Confnum)
    end
    return M
end

# Does the computations
function _ConfidenceBands!(Res::AbstractMatrix, DM::AbstractDataModel, points::AbstractVector{<:AbstractVector{<:Number}}, woundX::AbstractVector{<:Number})
    @assert pdim(DM) == ConsistentElDims(points)
    xdim(DM) != 1 && throw("Not programmed for xdim != 1 yet.")
    @assert size(Res) == (length(woundX), 2ydim(DM))
    # BELOW LOOP WILL NOT WORK FOR CompositeDataSet!!!!!!
    @assert Data(DM) isa DataSet || Data(DM) isa DataSetExact
    for point in points
        # Do it like this to exploit CustomEmbeddings
        Y = Windup(EmbeddingMap(Data(DM), Predictor(DM), point, woundX), ydim(DM)) |> Unpack
        for col in 1:2:2ydim(DM)
            Ycol = Int(ceil(col/2))
            for row in 1:size(Res,1)
                Res[row, col] = min(Res[row, col], Y[row,Ycol])
                Res[row, col+1] = max(Res[row, col+1], Y[row,Ycol])
            end
        end
    end
end

"""
Computes width of confidence bands.
"""
function ConfidenceBandWidth(args...; plot::Bool=true, OverWrite::Bool=true, kwargs...)
    band = ConfidenceBands(args...; plot=false, kwargs...)
    Res = hcat(view(band,:,1), view(band,:,3)-view(band,:,2))
    F = OverWrite ? Plots.plot : Plots.plot!
    plot && display(F(view(Res,:,1), view(Res,:,2), label="Conf. band width"))
    Res
end

"""
    PredictionEnsemble(DM::AbstractDataModel, pDomain::HyperCube, Xs::AbstractVector{<:Number}=DomainSamples(XCube(DM),300); N::Int=50, uniform::Bool=true, MaxConfnum::Real=3, plot::Bool=true, kwargs...)
Plots `N` model predictions which are randomly chosen from a confidence region of level `MaxConfnum`.
"""
function PredictionEnsemble(DM::AbstractDataModel, pDomain::HyperCube, Xs::AbstractVector{<:Number}=DomainSamples(XCube(DM),300); N::Int=50, uniform::Bool=true, MaxConfnum::Real=3, plot::Bool=true, kwargs...)
    @assert xdim(DM) == 1
    function GenerateUniformPoint(DM::AbstractDataModel, pDomain::HyperCube, MaxConfnum::Real)
        p = SVector{length(pDomain)}(rand(pDomain));        Conf = GetConfnum(DM, p)
        Conf ≤ MaxConfnum ? (p, Conf) : GenerateUniformPoint(DM, pDomain, MaxConfnum)
    end
    iF = inv(FisherMetric(DM, MLE(DM))) |> Symmetric
    function GenerateGaussianPoint(DM::AbstractDataModel, MaxConfnum::Real)
        p = rand(MvNormal(MLE(DM), iF));    Conf = GetConfnum(DM, p)
        Conf ≤ MaxConfnum ? (p, Conf) : GenerateGaussianPoint(DM, MaxConfnum)
    end
    function MakePrediction(DM::AbstractDataModel, pDomain::HyperCube, Xs::AbstractVector{<:Number}; MaxConfnum::Real=3)
        point, Conf = uniform ? GenerateUniformPoint(DM, pDomain, MaxConfnum) : GenerateGaussianPoint(DM, MaxConfnum)
        MakePrediction(DM, point, Xs), Conf
    end
    function MakePrediction(DM::AbstractDataModel, point::AbstractVector{<:Number}, Xs::AbstractVector{<:Number})
        Unpack(Windup(EmbeddingMap(Data(DM), Predictor(DM), point, Xs), ydim(DM)))
    end
    Preds = [MakePrediction(DM, pDomain, Xs; MaxConfnum=MaxConfnum) for i in 1:N]
    plot && display(PlotEnsemble(Xs, getindex.(Preds,1), getindex.(Preds,2); kwargs...))
    Xs, getindex.(Preds,1), getindex.(Preds,2)
end
PredictionEnsemble(DM::AbstractDataModel; MaxConfnum::Real=3, kwargs...) = PredictionEnsemble(DM, ProfileBox(DM,MaxConfnum); MaxConfnum=MaxConfnum, kwargs...)

function PlotEnsemble(Xs::AbstractVector{<:Number}, Preds::Union{AbstractVector{<:AbstractArray{<:Number}}}, Confs::AbstractVector{<:Real}; palette::Symbol=:Oranges_5, OverWrite::Bool=false, kwargs...)
    extr = extrema(Confs);    p = OverWrite ? Plots.plot() : Plots.plot!()
    for i in eachindex(Confs)
        p = Plots.plot!(Xs, Preds[i]; color=get(cgrad(palette), Confs[i], extr), label="", alpha=0.3, kwargs...)
    end;    p
end

PointwiseConfidenceBandFULL(DM::DataModel,sol::AbstractODESolution,Cube::HyperCube,Confnum::Real=1; N::Int=500) = PointwiseConfidenceBandFULL(DM,sol,FindMLE(DM),Cube,Confnum; N=N)
function PointwiseConfidenceBandFULL(DM::DataModel,sol::AbstractODESolution,MLE::AbstractVector,Cube::HyperCube,Confnum::Real=1; N::Int=500)
    !(length(Cube) == xdim(DM)) && throw("PWConfBand: Wrong Cube dim.")
    if ydim(DM) == 1
        Lims = ConstructCube(sol)
        low = Vector{Float64}(undef,N); up = Vector{Float64}(undef,N)
        X = DomainSamples(Cube; N=N)
        for i in 1:length(X)
            Y = Predictor(DM)(X[i],MLE)
            up[i] = maximum(Y); low[i] = minimum(Y)
        end
        LogLikeMLE = loglikelihood(DM,MLE)
        Confvol = ConfVol(Confnum)
        for i in 1:N
            num = rand.(Uniform.(Lims.L,Lims.U))
            if WilksTestPrepared(DM,num,LogLikeMLE,Confvol)
                for i in 1:length(X)
                    Y = Predictor(DM)(X[i],num)
                    Y > up[i] && (up[i] = Y)
                    Y < low[i] && (low[i] = Y)
                end
            else i = i-1 end
        end
        Plots.plot!(X,low)
        Plots.plot!(X,up) |> display
        return [X low up]
    else
        throw("Not programmed yet.")
        # Evaluate on boundary of cube
    end
end

"""
    PlotMatrix(Mat::AbstractMatrix, mle::AbstractVector; Confnum::Real=0., dims::Tuple{Int,Int}=(1,2), N::Int=400, plot::Bool=true, OverWrite::Bool=true, kwargs...)
Plots ellipse corresponding to a given covariance matrix which may additionally be offset by a vector `mle`.
By providing information on the confidence level of the given matrix via the `Confnum` kwarg, a correction factor is computed to rescale the given matrix appropriately.

Example:
```
PlotMatrix(inv(FisherMetric(DM,mle)),mle)
```
"""
function PlotMatrix(Mat::AbstractMatrix, MLE::AbstractVector{<:Number}=zeros(size(Mat,1)); Confnum::Real=0., dims::Tuple{Int,Int}=(1,2), N::Int=400, plot::Bool=true, OverWrite::Bool=true, kwargs...)
    !(length(MLE) == size(Mat,1) == size(Mat,2)) && throw("PlotMatrix: Dimensional mismatch.")
    corr = Confnum != 0. ? sqrt(quantile(Chisq(length(MLE)),ConfVol(Confnum))) : 1.0
    C = corr .* cholesky(Symmetric(Mat)).L
    angles = range(0, 2π; length=N)
    F(α::Number) = MLE + C * RotatedVector(α, dims[1], dims[2], length(MLE))
    Data = Unpack(F.(angles))
    Pl = OverWrite ? Plots.plot : Plots.plot!
    if plot   display(Pl(ToCols(Data)...; label="Matrix", kwargs...))  end
    Data
end

function PlotCurves(Curves::AbstractVector{<:AbstractODESolution}; N::Int=100)
    p = [];    A = Array{Float64,2}(undef,N,2)
    for sol in Curves
        ran = range(sol.t[1],sol.t[end],length=N)
        for i in 1:length(ran)    A[i,:] = sol(ran[i])[1:2]  end
        p = Plots.plot!(A[:,1],A[:,2])
        # p = Plots.plot!(sol,vars=(1,2))
    end
    p
end

EvaluateAlongGeodesic(F::Function, sol::AbstractODESolution, Interval::Tuple{<:Number,<:Number}=(sol.t[1],sol.t[end]); N::Int=300) = [F(sol(t)[1:Int(length(sol.u[1])/2)]) for t in range(Interval[1],Interval[2],length=N)]
function PlotAlongGeodesic(F::Function,sol::AbstractODESolution, Interval::Tuple{<:Number,<:Number}=(sol.t[1],sol.t[end]); N::Int=300, OverWrite::Bool=false)
    Z = EvaluateAlongGeodesic(F, sol, Interval; N=N)
    @assert ConsistentElDims(Z) == 1
    X = DomainSamples(Interval; N=N)
    OverWrite ? display(Plots.plot(X, Z)) : display(Plots.plot!(X, Z))
    [X Z]
end
EvaluateAlongGeodesicLength(DM::AbstractDataModel, F::Function, sol::AbstractODESolution, Interval::Tuple{<:Number,<:Number}=(sol.t[1],sol.t[end]); N::Int=300) = EvaluateAlongGeodesic(F,sol,Interval, N=N)
function PlotAlongGeodesicLength(DM::AbstractDataModel, F::Function, sol::AbstractODESolution, Interval::Tuple{<:Number,<:Number}=(sol.t[1],sol.t[end]); N::Int=300, OverWrite::Bool=false)
    Z = EvaluateAlongGeodesic(F, sol, Interval; N=N)
    @assert ConsistentElDims(Z) == 1
    X = DomainSamples(Interval; N=N)
    Geo = GeodesicLength(x->FisherMetric(DM,x), sol, sol.t[end]; FullSol=true, Auto=true, tol=1e-14)
    Ls = map(Geo, X)
    OverWrite ? display(Plots.plot(Ls, Z)) : display(Plots.plot!(Ls, Z))
    [Ls Z]
end
EvaluateAlongCurve(F::Function, sol::AbstractODESolution, Interval::Tuple{<:Number,<:Number}=(sol.t[1],sol.t[end]); N::Int=300) = [F(sol(t)) for t in range(Interval[1],Interval[2],length=N)]
function PlotAlongCurve(F::Function, sol::AbstractODESolution, Interval::Tuple{<:Number,<:Number}=(sol.t[1],sol.t[end]); N::Int=300, OverWrite::Bool=false)
    Z = EvaluateAlongCurve(F, sol, Interval, N=N)
    @assert ConsistentElDims(Z) == 1
    X = DomainSamples(Interval; N=N)
    OverWrite ? display(Plots.plot(X, Z)) : display(Plots.plot!(X, Z))
    [X Z]
end

PhaseSpacePlot(DM::AbstractDataModel; kwargs...) = PhaseSpacePlot(DM, (C=InformationGeometry.XCube(DM); range(C.L[1], C.U[1]; length=300)); kwargs...)
function PhaseSpacePlot(DM::AbstractDataModel, ts::AbstractVector{<:Number}, mle::AbstractVector{<:Number}=MLE(DM); OverWrite::Bool=true, kwargs...)
    OverWrite && plot()
    p = scatter!(collect(Iterators.partition(ydata(DM),ydim(DM))); label="Observed Data")
    p = if ydim(DM) == 2
        plot!(collect(Iterators.partition(EmbeddingMap(DM, mle, ts),ydim(DM))); label="Phase Space Trajectory", xlabel="$(ynames(DM)[1])", ylabel="$(ynames(DM)[2])", markeralpha=0, linealpha=1, kwargs...)
    elseif ydim(DM) == 3
        plot!(collect(Iterators.partition(EmbeddingMap(DM, mle, ts),ydim(DM))); label="Phase Space Trajectory", xlabel="$(ynames(DM)[1])", ylabel="$(ynames(DM)[2])", zlabel="$(ynames(DM)[3])", markeralpha=0, linealpha=1, kwargs...)
    else
        throw("Cannot display phase space for ydim=$(ydim(DM)).")
    end;    display(p);    p
end

# Helper methods for creating a rectangular or triangular mesh from the outputs of MincedBoundaries()

# function RectangularFaceIndices(N::Int=20, zerobased::Bool=false)
#     G = Matrix{Int64}(undef,N,4)
#     for i in 1:(N-1)
#         G[i,:] = SA[i, i+1, N+i+1, N+i]
#     end;    G[N,:] = SA[N, 1, N+1, 2N]
#     zerobased && (G .-= 1)
#     G
# end
"""
Returns Array of vertex numbers where every row constitutes a trapezoid for two adjacent curves from which N samples have been drawn.
"""
RectangularFaceIndices(N::Int, zerobased::Bool=false) = (G=Matrix{Int64}(undef,N,4); RectangularFaceIndices!(G, zerobased); G)
function RectangularFaceIndices!(G::AbstractMatrix{<:Int}, zerobased::Bool=false)
    N, cols = size(G);    @assert cols == 4
    @inbounds for i in Base.OneTo(N-1)
        G[i,:] = SA[i, i+1, N+i+1, N+i]
    end;    G[N,:] = SA[N, 1, N+1, 2N]
    zerobased && (G .-= 1)
    nothing
end

"""
Turns Array which specifies rectangular faces into triangular faces.
"""
function RectToTriangFaces(M::AbstractMatrix{<:Int})
    G = Matrix{Int64}(undef, 2size(M,1), 3)
    @inbounds for i in Base.OneTo(size(M,1))
        G[2i-1,:], G[2i,:] = SA[M[i,1], M[i,2], M[i,4]], SA[M[i,2], M[i,3], M[i,4]]
    end;    G
end
"""
    CreateMesh(Planes::AbstractVector{<:Plane}, Sols::AbstractVector{<:AbstractODESolution}; N::Int=2*length(Sols), rectangular::Bool=true) -> (Matrix{Float64}, Matrix{Int64})
Returns a N×3 matrix whose rows correspond to the coordinates of various points in 3D space as the first argument.
The second Matrix is either N×4 or N×3 depending on the value of `rectangular` and enumerates the points which are to be connected up to a rectangular or triangular face in counter-clockwise fashion. The indices of the points correspond to the lines in the first Matrix.
"""
function CreateMesh(Planes::AbstractVector{<:Plane}, Sols::AbstractVector{<:AbstractODESolution}, fact::Real=1.4; N::Int=2*length(Sols), rectangular::Bool=true, pointy::Bool=true)
    Vertices = reduce(vcat, [Deplanarize(Planes[i], Sols[i]; N=N) for i in 1:length(Sols)])
    M = RectangularFaceIndices(N)
    Faces = reduce(vcat, [M .+ (i-1)*N for i in 1:(length(Sols)-1)])
    !rectangular && (Faces = RectToTriangFaces(Faces))
    pointy && return AddCaps(Planes, Sols, Vertices, Faces, fact; N=N)
    Vertices, Faces
end
function AddCaps(Planes::AbstractVector{<:Plane}, Sols::AbstractVector{<:AbstractODESolution}, Vertices::AbstractMatrix, Faces::AbstractMatrix, fact::Real=1.4; N::Int=2*length(Sols))
    @assert length(Planes[1]) == 3 # Can relax this if cross() is extended to higher dims
    @assert fact ≥ 1
    linep1 = size(Vertices,1) + 1
    function rowmean(PL::Plane, sol::AbstractODESolution, N::Int)
        PlaneCoordinates(PL, sum(sol(t) for t in range(sol.t[1], sol.t[end]; length=N))/N)
    end
    # add two points on the top and bottom of the confidence region
    n = cross(Planes[1].Vx, Planes[1].Vy) |> normalize
    height = dot(n, Planes[1].stütz - Planes[2].stütz)
    p1 = rowmean(Planes[1], Sols[1], N) + (fact-1)*height*n
    p2 = rowmean(Planes[end], Sols[end], N) - (fact-1)*height*n
    # p1 = (fact-1) * (Planes[1].stütz - Planes[2].stütz) + Planes[1].stütz
    # p2 = (fact-1) * (Planes[end].stütz - Planes[end-1].stütz) + Planes[end].stütz

    connectp1, connectp2 = if size(Faces,2) == 4
        vcat(reduce(vcat, [[i i+1 linep1 linep1] for i in 1:(N-1)]), [N 1 linep1 linep1]),
        vcat(reduce(vcat, [[i i+1 linep1+1 linep1+1] for i in ((length(Sols)-1)*N + 1):(linep1-2)]), [linep1-1 ((length(Sols)-1)*N+1) linep1+1 linep1+1])
    else
        vcat(reduce(vcat, [[i i+1 linep1] for i in 1:(N-1)]), [N 1 linep1]),
        vcat(reduce(vcat, [[i i+1 linep1+1] for i in ((length(Sols)-1)*N + 1):(linep1-2)]), [linep1-1 ((length(Sols)-1)*N+1) linep1+1])
    end
    vcat(Vertices, vcat(transpose(p1), transpose(p2))), vcat(Faces, connectp1, connectp2)
end

function ToObj(Vertices::AbstractMatrix, Faces::AbstractMatrix)
    text = ""
    for i in 1:size(Vertices,1)
        text *= "v"
        for j in 1:size(Vertices,2)
            text *= " $(Vertices[i,j])"
        end
        text *= '\n'
    end
    # ONE-BASED facet indices for .obj files
    for i in 1:size(Faces,1)
        text *= "f"
        for j in 1:size(Faces,2)
            text *= " $(Faces[i,j])"
        end
        text *= '\n'
    end;    text
end
function WriteObj(Vertices::AbstractMatrix, Faces::AbstractMatrix, path::String="D:/Boundary.obj")
    open(path,"w") do f
        write(f,ToObj(Vertices,Faces))
    end;    return
end

"""
    PolarSolution(Planes::AbstractVector{<:Plane}, sols::AbstractVector{<:AbstractODESolution}) -> Function
Returns spherical function `f: [0,π]×[0,2π] ⟶ R^3` `(θ, ϕ) ⟼ [x(θ, ϕ), y(θ, ϕ), z(θ, ϕ)]` which returns closest point on confidence boundary corresponding to the given angular coordinates.
"""
function PolarSolution(Planes::AbstractVector{<:Plane}, sols::AbstractVector{<:AbstractODESolution})
    @assert 3 == length(Planes[1])
    len = length(Planes)
    PolarToInd(θ::Real, n::Int=len) = ceil(Int, θ/π * (n-1)) + 1
    function SurfacePoint(θ::Real, ϕ::Real)
        @assert 0 ≤ θ ≤ π "θ not in range."
        ind = PolarToInd(θ)
        SurfacePoint(Planes[ind], sols[ind], ϕ)
    end
    function SurfacePoint(PL::Plane, sol::AbstractODESolution, ϕ::Real)
        # @assert 0 ≤ ϕ ≤ 2π "ϕ not in range."
        PlaneCoordinates(PL, sol(sol.t[end] * ϕ / (2π)))
    end
end

function AddCaps(X::M, Y::M, Z::M, fact::Real=1.4) where M <: AbstractMatrix{<:Number}
    # thetas varies along rows, phis varies along columns
    @assert size(X) == size(Y) == size(Z)
    @assert fact ≥ 1
    function rowmean(X::M, Y::M, Z::M, i::Int) where M <: AbstractMatrix{<:Number}
        point = zeros(3)
        for j in 1:size(X,2)
            point += SA[X[i,j], Y[i,j], Z[i,j]]
        end;    point / size(X,2)
    end
    # cap1 = (z = rowmean(X,Y,Z,2);   z + fact*(rowmean(X,Y,Z,1) - z))
    # cap2 = (z = rowmean(X,Y,Z, size(X,1)-1);   z + fact*(rowmean(X,Y,Z, size(X,1)) - z))
    n = (i=1; normalize(cross([X[i,3], Y[i,3], Z[i,3]]-[X[i,2], Y[i,2], Z[i,2]], [X[i,1], Y[i,1], Z[i,1]]-[X[i,2], Y[i,2], Z[i,2]])))
    height = dot(n, [X[1,1], Y[1,1], Z[1,1]]-[X[2,1], Y[2,1], Z[2,1]])
    cap1 = rowmean(X, Y, Z, 1) + (fact-1)*height*n
    cap2 = rowmean(X, Y, Z, size(X,1)) - (fact-1)*height*n
    vcat(transpose(cap1[1]*ones(size(X,2))), X, transpose(cap2[1]*ones(size(X,2)))),
    vcat(transpose(cap1[2]*ones(size(X,2))), Y, transpose(cap2[2]*ones(size(X,2)))),
    vcat(transpose(cap1[3]*ones(size(X,2))), Z, transpose(cap2[3]*ones(size(X,2))))
end

"""
    ToMatrices(Planes::AbstractVector{<:Plane}, sols::AbstractVector{<:AbstractODESolution}; N::Int=2*length(sols)) -> Tuple{Matrix,Matrix,Matrix}
Returns tuple of three `length(Planes) × N` matrices `X, Y, Z` for plotting of confidence boundary surfaces in 3D with tools such as Makie.
"""
function ToMatrices(Planes::AbstractVector{<:Plane}, sols::AbstractVector{<:AbstractODESolution}, fact::Real=1.4; N::Int=2*length(sols), pointy::Bool=true)
    θs = range(0, π; length=size(Planes,1));    ϕs = range(0, 2π; length=N);    P = PolarSolution(Planes, sols)
    X = Array{Float64}(undef, length(θs), length(ϕs));    Y = Array{Float64}(undef, length(θs), length(ϕs));    Z = Array{Float64}(undef, length(θs), length(ϕs))
    for i in 1:length(θs)
        for j in 1:length(ϕs)
            # X[i,j], Y[i,j], Z[i,j] = P(θs[i], ϕs[j])
            X[i,j], Y[i,j], Z[i,j] = P(Planes[i], sols[i], ϕs[j])
        end
    end
    pointy && return AddCaps(X, Y, Z, fact)
    (X, Y, Z)
end


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule argp_standalone_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("argp_standalone")
JLLWrappers.@generate_main_file("argp_standalone", UUID("c53206cc-00f7-50bf-ad1e-3ae1f6e49bc3"))
end  # module argp_standalone_jll

ext: jl


text:
# Autogenerated wrapper script for AMGX_jll for x86_64-linux-gnu-cxx11
export libamgxsh

using CUDA_jll
JLLWrappers.@generate_wrapper_header("AMGX")
JLLWrappers.@declare_library_product(libamgxsh, "libamgxsh.so")
function __init__()
    JLLWrappers.@generate_init_header(CUDA_jll)
    JLLWrappers.@init_library_product(
        libamgxsh,
        "lib/libamgxsh.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
mutable struct SimulationState
    g::Global
    population::Population
    setup::Evaluation
    evaluate!::Function
    function SimulationState(g::Global)
        p = Population(g, g.cf.checkpoint_filename)
        setup = Evaluation(g.cf.dataset_id, g.cf.fitness_function_id, normalise=g.cf.normalise_input,
                           balance=g.cf.balance_dataset)
        fitness_function = (f::Function) -> f.(p.population)
        new(g, p, setup, fitness_function)
    end
    function SimulationState(g::Global, p::Population)
        setup = Evaluation(g.cf.dataset_id, g.cf.fitness_function_id, normalise=g.cf.normalise_input,
                           balance=g.cf.balance_dataset)
        fitness_function = (f::Function) -> f.(p.population)
        new(g, p, setup, fitness_function)
    end
end

function Base.show(io::IO, state::SimulationState)
    print(io, "State: \n")
    println(io, state.population)
    print(io, state.setup.dataset_setup)
end

max_stars_repo_name: Lucas-Muniz/MNEAT

max_issues_repo_name: Lucas-Muniz/MNEAT

max_forks_repo_name: Lucas-Muniz/MNEAT


text:
if VERSION >= v"0.7.0"
    using Printf
    using LinearAlgebra
end

function descent( z::Array{T,1}, alpha::Array{T,1}, c::SLiSeFilters.Config, lb::R=R(0.0) )

    project = function( z::Array{T,1} )
        mask = abs.(imag(z)) .< lb
        z[mask] = real(z[mask]) + lb*im*sign.(imag(z[mask]));
        z
    end

    maxiter = 20000;
    i = 1;

    while i < maxiter
        dz,da = Fgradient( z, alpha, c );

        a = 0.5;
        b = 0.5;
        t::Float64 = 1*b;

        resd = residual(z, alpha, c);
        #emm = dot(-([da; dz]), [da; dz]);
        emm = -norm([da; dz])^2;

        while( residual( project(z-t .* dz), alpha - t.*da, c) > real(resd + a * t * emm) )
            t = t*b;
            if 1e-16 > t
                println("done...")
                i=maxiter;
                break;
            end
        end

        alpha -= t*da;
        z = project(z - t*dz);

        if mod(i,500) == 0
            @printf "%.16e\t%d\n" resd i
        end
        i+=1
    end

    return z,alpha
end


text:
Tables.istable(::Type{<:ResultSet}) = true

Tables.rowaccess(::Type{<:ResultSet}) = true
Tables.columnaccess(::Type{<:ResultSet}) = true

Tables.rows(rs::ResultSet) = rs
Tables.columns(rs::ResultSet) = rs

Tables.schema(rs::ResultSet) = Tables.Schema(names(rs), eltype.(columns(rs)))


text:
module TestAqua

using Aqua
using BangBang
using Setfield

Aqua.test_all(
    BangBang;
    ambiguities=(exclude=[Base.get, Setfield.set, Setfield.modify],),
)

end  # module


text:
export first, FirstNotFoundException

import Base: first

struct FirstNotFoundException <: Exception end

"""
    first(; default = nothing)

Creates a first operator, which returns an Observable
that emits only the first value emitted by the source Observable.
Sends `FirstNotFoundException` error message if a given source completes without emitting a single value.

# Arguments
- `default`: an optional default value to provide if no values were emitted

# Producing

Stream of type `<: Subscribable{L}` where `L` refers to type of source stream

# Examples
```jldoctest
using Rocket

source       = from(1:100)
subscription = subscribe!(source |> first(), logger())
;

# output

[LogActor] Data: 1
[LogActor] Completed
```

See also: [`take`](@ref), [`logger`](@ref)
"""
first(; default = nothing) = default === nothing ? take(1) + error_if_empty(FirstNotFoundException()) : take(1) + default_if_empty(default)


text:
mutable struct ZeroMean{T<:Real} <: PriorMean{T}
end

"""
ZeroMean
```julia
ZeroMean()
```
Construct a mean prior set to 0 and cannot be changed.
"""
function ZeroMean()
    ZeroMean{Float64}()
end

update!(μ::ZeroMean{T},grad::AbstractVector{T}) where {T<:Real} = nothing

get_opt(μ::ZeroMean) = nothing

array(μ::ZeroMean{T},length::Int) where {T<:Real} = zeros(T,length)

(μ::ZeroMean{T})(x::Real) where {T<:Real} = zero(T)
(μ::ZeroMean{T})(x::AbstractMatrix) where {T<:Real} = zeros(T,size(x,1))


text:
include("../aimajulia.jl");

using Test;

using Main.aimajulia;

#The following search tests are from the aima-python doctest

@test depth_first_tree_search(NQueensProblem(8)) == Node{Array{Int64, 1}}([8, 4, 1, 3, 6, 2, 7, 5]);

#Specify wordlist path for travis-ci testing.
filename = "./aima-data/EN-text/wordlist.txt";

@test length(BoggleFinder(board=collect("SARTELNID"), fn=filename)) == 206;

ab = GraphProblem("A", "B", aimajulia.romania);

@test solution(breadth_first_tree_search(ab)) == ["S", "F", "B"];

@test solution(breadth_first_search(ab)) == ["S", "F", "B"];

@test solution(uniform_cost_search(ab)) == ["S", "R", "P", "B"];

@test solution(depth_first_graph_search(ab)) == ["T", "L", "M", "D", "C", "P", "B"];

@test solution(iterative_deepening_search(ab)) == ["S", "F", "B"];

@test length(solution(depth_limited_search(ab))) == 50;

@test solution(astar_search(ab)) == ["S", "R", "P", "B"];

@test solution(recursive_best_first_search(ab)) == ["S", "R", "P", "B"];

@test compare_searchers([GraphProblem("A", "B", aimajulia.romania),
                        GraphProblem("O", "N", aimajulia.romania),
                        GraphProblem("Q", "WA", aimajulia.australia)],
                        ["Searcher", "Romania(A, B)", "Romania(O, N)", "Australia"]) ==
    ["Searcher"                     "Romania(A, B)"         "Romania(O, N)"         "Australia";
    "Main.aimajulia.breadth_first_tree_search"     "<  23/  24/  63/B>"    "<1191/1192/3378/N>"    "<   9/  10/  32/WA>";
    "Main.aimajulia.breadth_first_search"          "<   7/  11/  18/B>"    "<  18/  20/  44/N>"    "<   3/   6/   9/WA>";
    "Main.aimajulia.depth_first_graph_search"      "<   8/   9/  20/B>"    "<  16/  17/  37/N>"    "<   2/   3/   8/WA>";
    "Main.aimajulia.iterative_deepening_search"    "<  13/  36/  36/B>"    "< 683/1874/1875/N>"    "<   4/  13/  12/WA>";
    "Main.aimajulia.depth_limited_search"          "<  64/  94/ 167/B>"    "< 948/2629/2701/N>"    "<  51/  57/ 153/WA>";
    "Main.aimajulia.recursive_best_first_search"   "<  11/  12/  35/B>"    "<8481/8482/23788/N>"   "<  10/  11/  38/WA>"];

# Initialize LRTAStarAgentProgram with an OnlineSearchProblem.

lrtastar_program = OnlineSearchProblem("State_3", "State_5", aimajulia.one_dim_state_space, aimajulia.one_dim_state_space_least_costs);
lrtastar_agentprogram = LRTAStarAgentProgram(lrtastar_program);

@test execute(lrtastar_agentprogram, "State_3") == "Right";

@test execute(lrtastar_agentprogram, "State_4") == "Left";

@test execute(lrtastar_agentprogram, "State_3") == "Right";

@test execute(lrtastar_agentprogram, "State_4") == "Right";

@test execute(lrtastar_agentprogram, "State_5") == nothing;

lrtastar_agentprogram = LRTAStarAgentProgram(lrtastar_program);

@test execute(lrtastar_agentprogram, "State_4") == "Left";

lrtastar_agentprogram = LRTAStarAgentProgram(lrtastar_program);

@test execute(lrtastar_agentprogram, "State_5") == nothing;

max_stars_repo_name: mikhail-j/aima-julia

max_issues_repo_name: mikhail-j/aima-julia

max_forks_repo_name: mikhail-j/aima-julia


text:
module EMLheat_flow
	using DanaTypes
	using DotPlusInheritance
	using Reexport
	@reexport using ...streams.EMLstreams
	import EMLtypes.length
	include("heat_flow/heat_flow.jl")
end


text:
using ImageSegmentation
using ImageSegmentation.Colors
using ImageSegmentation.FixedPointNumbers
using FileIO
using Statistics
using SparseArrays
using Test

@testset "flood_fill" begin
    # 0d
    a = reshape([true])
    @test flood(identity, a, CartesianIndex()) == a
    @test_throws ArgumentError flood(!, a, CartesianIndex())
    # 1d
    a = 1:7
    @test flood(==(2), a, CartesianIndex(2)) == (a .== 2)
    @test_throws ArgumentError flood(==(2), a, CartesianIndex(3))
    @test flood(x -> 1 < x < 4, a, CartesianIndex(2)) == [false, true, true, false, false, false, false]
    @test flood(isinteger, a, CartesianIndex(2)) == trues(7)
    # 2d
    ab = [true false false false;
         true true false false;
         true false false true;
         true true true true]
    an0f8 = N0f8.(ab)
    agray = Gray.(an0f8)
        for (f, a) in ((identity, ab), (==(1), an0f8), (==(1), agray))
        for idx in CartesianIndices(a)
            if f(a[idx])
                @test flood(f, a, idx) == a
            else
                @test_throws ArgumentError flood(f, a, idx)
            end
        end
    end
    @test flood(identity, ab, Int16(1)) == ab
    # 3d
    k = 10
    a = falses(k, k, k)
    idx = CartesianIndex(1,1,1)
    incs = [CartesianIndex(1,0,0), CartesianIndex(0,1,0), CartesianIndex(0,0,1)]
    a[idx] = true
    while any(<(k), Tuple(idx))
        d = rand(1:3)
        idx += incs[d]
        idx = min(idx, CartesianIndex(k,k,k))
        a[idx] = true
    end
    for idx in eachindex(a)
        if a[idx]
            @test flood(identity, a, idx) == a
        else
            @test_throws ArgumentError flood(identity, a, idx)
        end
    end
    # Colors
    path = download("https://github.com/JuliaImages/juliaimages.github.io/raw/source/docs/src/pkgs/segmentation/assets/flower.jpg")
    img = load(path)
    seg = flood(img, CartesianIndex(87,280); thresh=0.3*sqrt(3))   # TODO: eliminate the sqrt(3) when we transition to `abs2(c) = c ⋅ c`
    @test 0.2*length(seg) <= sum(seg) <= 0.25*length(seg)
    c = mean(img[seg])
    # N0f8 makes for easier approximate testing
    @test N0f8(red(c)) ≈ N0f8(0.855)
    @test N0f8(green(c)) ≈ N0f8(0.161)
    @test N0f8(blue(c)) ≈ N0f8(0.439)

    # flood_fill!
    near3(x) = round(Int, x) == 3
    a0 = [range(2, 4, length=9);]
    a = copy(a0)
    idx = (length(a)+1)÷2
    dest = fill!(similar(a, Bool), false)
    @test flood_fill!(near3, dest, a, idx) == (round.(a) .== 3)
    a = copy(a0)
    flood_fill!(near3, a, idx; fillvalue=-1)
    @test a == [near3(a0[i]) ? -1 : a[i] for i in eachindex(a)]
    a = copy(a0)
    @test_throws ArgumentError flood_fill!(near3, a, idx; fillvalue=-1, isfilled=near3)
    # warning
    a = [1:7;]
    @test_logs (:warn, r"distinct.*incomplete") flood_fill!(<(5), a, 1; fillvalue=3)
    @test a == [3,3,3,4,5,6,7]
    a = [1:7;]
    dest = fill(-1, size(a))
    @test_logs flood_fill!(<(5), dest, a, 1; fillvalue=3)   # no warnings
    @test dest == [3,3,3,3,-1,-1,-1]
    a = [1:7;]
    @test_logs flood_fill!(<(5), a, 1; fillvalue=11)
    @test a == [11,11,11,11,5,6,7]

    # This mimics a "big data" application in which we have several structures we want
    # to label with different segment numbers, and the `src` array is too big to fit
    # in memory.
    # It would be better to use a package like SparseArrayKit, which allows efficient
    # insertions and supports arbitrary dimensions.
    a = Bool[0 0 0 0 0 0 1 1;
             1 1 0 0 0 0 0 0]
    dest = spzeros(Int, size(a)...)   # stores the nonzero indexes in a Dict
    flood_fill!(identity, dest, a, CartesianIndex(2, 1); fillvalue=1)
    flood_fill!(identity, dest, a, CartesianIndex(1, 7); fillvalue=2)
    @test dest == [0 0 0 0 0 0 2 2;
                   1 1 0 0 0 0 0 0]
end


text:
const strip_patterns            = UInt32(0)
# Flags that activate function-based processors
const strip_corrupt_utf8        = UInt32(0x1) << 0
const strip_case                = UInt32(0x1) << 1
const strip_accents             = UInt32(0x1) << 2
const strip_html_tags           = UInt32(0x1) << 3
# Flags that activate function-based processors (external to this file)
const stem_words                = UInt32(0x1) << 7
# Flags that activate Regex based processors
const strip_punctuation         = UInt32(0x1) << 9
const strip_whitespace          = UInt32(0x1) << 10
const strip_numbers             = UInt32(0x1) << 11
const strip_non_ascii           = UInt32(0x1) << 12
const strip_single_chars        = UInt32(0x1) << 13
# Word list based
const strip_indefinite_articles = UInt32(0x1) << 20
const strip_definite_articles   = UInt32(0x1) << 21
const strip_prepositions        = UInt32(0x1) << 22
const strip_pronouns            = UInt32(0x1) << 23
const strip_stopwords           = UInt32(0x1) << 24
const strip_sparse_terms        = UInt32(0x1) << 25
const strip_frequent_terms      = UInt32(0x1) << 26

# Generate custom flag combinations (bit-shift)
function flag_generate(bs::Vector{<:Integer})
    b = UInt32(1)
    n = UInt32(0)
    for s in bs
        n|= (b << s)
    end
    return n
end

# Generate custom flag combinations (bit-toogle)
flag_generate(flags::UInt32...) = reduce(|, flags)

# Compound stripping flags
const strip_articles = flag_generate(
                        strip_indefinite_articles,
                        strip_definite_articles)

const strip_everything = flag_generate(
                            [0,1,2,3,
                            9,10,11,12,13,
                            20,21,22,23,24])
const strip_everything_stem = flag_generate(
                                [0,1,2,3,
                                 7,
                                 9,10,11,12,13,
                                 20,21,22,23,24])

# RegEx Expressions for various stripping flags
# Format: flag => (match=>replacement)
const strip2regex = Dict{UInt32,Regex}(
    strip_whitespace => r"[\s]+",
    strip_numbers => r"\d+",
    strip_non_ascii => r"[^a-zA-Z\s]",
    strip_single_chars => r"(\s|\b)[\w]{1}(\b|\s)",
    #strip_html_tags => r"(<script\b[^>]*>([\s\S]*?)</script>|<[^>]*>)",
    #strip_punctuation =>r"[^\d\w\s\b]+"
    strip_punctuation => r"[!\"#$%&\'()*+,-./:;<=>?@\[\\\]^_`\{\|\}~]+"
)

# Basic string processing functions
# Remove corrupt UTF8 characters
remove_corrupt_utf8(s::AbstractString) = begin
    return map(x->isvalid(x) ? x : ' ', s)
end

# Conversion to lowercase
remove_case(s::T) where T<:AbstractString = lowercase(s)

# Removing accents
remove_accents(s::T) where T<:AbstractString =
    Unicode.normalize(s, stripmark=true)

# Remove HTML tags
remove_html_tags(s::T) where T<:AbstractString =
    replace(s, r"(<script\b[^>]*>([\s\S]*?)</script>|<[^>]*>)" => ' ')

# Generate automatically functions for various Document types and Corpus
# Note: One has to add a simple method for `AbstractString` and the name
#       of the function in the `for` container to generate all needed
#       methods
for fname in [:remove_corrupt_utf8, :remove_case, :remove_accents, :remove_html_tags]
    # File document
    definition = """
        $(fname)!(d::FileDocument) = error("FileDocument cannot be modified.")
        """
    eval(Meta.parse(definition))
    # String Document
    definition = """
        function $(fname)!(d::StringDocument)
            d.text = $(fname)(d.text)
            return nothing
        end
        """
    eval(Meta.parse(definition))
    # Token Document
    definition = """
        function $(fname)!(d::TokenDocument)
            @inbounds for i in 1:length(d.tokens)
                d.tokens[i] = $(fname)(d.tokens[i])
            end
            filter!(t->(!isempty(t) && isvalid(t[1])), d.tokens)
        end
        """
    eval(Meta.parse(definition))
    # NGramDocument
    definition = """
        function $(fname)!(d::NGramDocument{S}) where S
            _ngrams = Dict{S, Int}()
            for token in keys(d.ngrams)
                _token = $(fname)(token)
                if !isempty(_token)
                    _ngrams[_token] = get(_ngrams, _token, 0) + 1
                end
            end
            filter!(p->isvalid(p.first[1]), _ngrams)
            d.ngrams = _ngrams
            return nothing
        end
        """
    eval(Meta.parse(definition))
    # Corpus
    definition = """
        function $(fname)!(crps::Corpus)
            for doc in crps
                $(fname)!(doc)
            end
        end
        """
    eval(Meta.parse(definition))
end

function write_sub(to::IOBuffer, a::AbstractArray{UInt8}, offs, nel)
     if offs+nel-1 > length(a) || offs < 1 || nel < 0
         throw(BoundsError())
     end
     GC.@preserve a unsafe_write(to, pointer(a, offs), UInt(nel))
 end

"""
    remove_patterns(s, rex)

Removes from the string `s` the text matching the pattern described
by the regular expression `rex`.
"""
function remove_patterns(s::AbstractString, rex::Regex)
    iob = IOBuffer()
    ibegin = 1
    v=codeunits(s)
    for m in eachmatch(rex, s, overlap=true)
        len = m.match.offset-ibegin+1
        if len > 0
            write_sub(iob, v, ibegin, len)
            write(iob, ' ')
        end
        ibegin = nextind(s, lastindex(m.match)+m.match.offset)
    end
    len = length(v) - ibegin + 1
    (len > 0) && write_sub(iob, v, ibegin, len)
    String(take!(iob))
end

function remove_patterns(s::SubString{T}, rex::Regex) where T <: String
    iob = IOBuffer()
    ioffset = s.offset
    data = codeunits(s.string)
    ibegin = 1
    for m in eachmatch(rex, s, overlap=true)
        len = m.match.offset-ibegin+1
        if len > 0
            write_sub(iob, data, ibegin+ioffset, len)
            write(iob, ' ')
        end
        ibegin = nextind(s, lastindex(m.match)+m.match.offset)
    end
    len = lastindex(s) - ibegin + 1
    (len > 0) && write_sub(iob, data, ibegin+ioffset, len)
    String(take!(iob))
end

"""
    remove_patterns!(d, rex)

Removes from the document or corpus `d` the text matching the pattern described
by the regular expression `rex`.
"""
remove_patterns!(d::FileDocument, rex::Regex) = error("FileDocument cannot be modified.")

remove_patterns!(d::StringDocument, rex::Regex) = begin
    d.text = remove_patterns(d.text, rex)
    nothing
end

remove_patterns!(d::TokenDocument, rex::Regex) = begin
    @inbounds for i in 1:length(d.tokens)
        d.tokens[i] = remove_patterns(d.tokens[i], rex)
    end
    filter!(t->(!isempty(t) && isvalid(t[1])), d.tokens)
end

remove_patterns!(d::NGramDocument{S}, rex::Regex) where S = begin
    _ngrams = Dict{S, Int}()
    for token in keys(d.ngrams)
        _token = remove_patterns(token, rex)
        if !isempty(_token)
            _ngrams[_token] = get(_ngrams, _token, 0) + 1
        end
    end
    filter!(p->isvalid(p.first[1]), _ngrams)
    d.ngrams = _ngrams
    return nothing
end

function remove_patterns!(crps::Corpus, rex::Regex)
    for doc in crps
        remove_patterns!(doc, rex)
    end
end

# Remove specified words
function remove_words!(entity, words::Vector{T}) where T<: AbstractString
    skipwords = Vector{T}()
    union!(skipwords, words)
    prepare!(entity, strip_patterns, skip_words = skipwords)
end

"""
    sparse_terms(crps::Corpus, alpha)

Returns a vector with rare terms among all documents. The parameter
`alpha` indicates the sparsity threshold (a frequency <= alpha means sparse).
"""
function sparse_terms(crps::Corpus, alpha=DEFAULT_CORPUS_SPARSITY)
    isempty(crps.lexicon) && update_lexicon!(crps)
    isempty(crps.inverse_index) && update_inverse_index!(crps)
    res = Vector{String}(undef, 0)
    ndocs = length(crps.documents)
    for term in keys(crps.lexicon)
        f = length(crps.inverse_index[term]) / ndocs
        if f <= alpha
            push!(res, String(term))
        end
    end
    return res
end

"""
    frequent_terms(crps::Corpus, alpha)

Returns a vector with frequent terms among all documents. The parameter
`alpha` indicates the sparsity threshold (a frequency <= alpha means sparse).
"""
function frequent_terms(crps::Corpus, alpha=1.0-DEFAULT_CORPUS_SPARSITY)
    isempty(crps.lexicon) && update_lexicon!(crps)
    isempty(crps.inverse_index) && update_inverse_index!(crps)
    res = Vector{String}(undef, 0)
    ndocs = length(crps.documents)
    for term in keys(crps.lexicon)
        f = length(crps.inverse_index[term]) / ndocs
        if f > alpha
            push!(res, String(term))
        end
    end
    return res
end

"""
    sparse_terms(doc, alpha)

Returns a vector with rare terms in the document `doc`. The parameter
`alpha` indicates the sparsity threshold (a frequency <= alpha means sparse).
"""
function sparse_terms(doc, alpha=DEFAULT_DOC_SPARSITY)
    ng = ngrams(doc)
    n = sum(values(ng))
    res = Vector{String}(undef, 0)
    for (term, count) in ng
        if count/n <= alpha
            push!(res, String(term))
        end
    end
    return res
end

"""
    frequent_terms(doc, alpha)

Returns a vector with frequent terms in the document `doc`. The parameter
`alpha` indicates the sparsity threshold (a frequency <= alpha means sparse).
"""
function frequent_terms(doc, alpha=1.0-DEFAULT_DOC_SPARSITY)
    ng = ngrams(doc)
    n = sum(values(ng))
    res = Vector{String}(undef, 0)
    for (term, count) in ng
        if count/n > alpha
            push!(res, String(term))
        end
    end
    return res
end

# Function that builds a regex out of a set of strings
_build_words_pattern(words::Vector{T}) where T<:AbstractString = begin
    Regex(ifelse(isempty(words), "", "\\b("* join(words,"|","|") *")\\b"))
end

# Function that builds a big regex out of a set of regexes
_build_regex_pattern(regexes::Vector{T}) where T<:Regex = begin
    l = length(regexes)
    if l == 0
        return r""
    elseif l == 1
        return pop!(regexes)
    else
        iob = IOBuffer()
        write(iob, "($(pop!(regexes).pattern))")
        for re in regexes
            write(iob, "|($(re.pattern))")
        end
        return Regex(String(take!(iob)))
    end
end

# Get the language
_language(doc::AbstractDocument) = doc.metadata.language
_language(crps::Corpus) = begin
    length(crps) < 1 && return DEFAULT_LANGUAGE
    return _language(crps[1])
end

function prepare!(entity,  # can be an AbstractDocument or Corpus
                  flags::UInt32;
                  skip_patterns=Vector{Regex}(),
                  skip_words=Vector{String}(),
                  alpha_sparse=0.05,
                  alpha_frequent=0.95)
    # Do function-based stripping
    ((flags & strip_corrupt_utf8) > 0) && remove_corrupt_utf8!(entity)
    ((flags & strip_html_tags) > 0) && remove_html_tags!(entity)
    ((flags & strip_case) > 0) && remove_case!(entity)
    ((flags & strip_accents) > 0) && remove_accents!(entity)
    # regex
    rpatterns = Vector{Regex}(undef, 0)  # patterns to remove
    ((flags & strip_whitespace) > 0) && push!(rpatterns, strip2regex[strip_whitespace])
    if (flags & strip_non_ascii) > 0
        push!(rpatterns, strip2regex[strip_non_ascii])
    else
        ((flags & strip_numbers) > 0) && push!(rpatterns, strip2regex[strip_numbers])
        ((flags & strip_punctuation) > 0) && push!(rpatterns, strip2regex[strip_punctuation])
        ((flags & strip_single_chars) > 0) && push!(rpatterns, strip2regex[strip_single_chars])
    end
    # known words
    language = _language(entity)
    if (flags & strip_articles) > 0
        union!(skip_words, articles(language))
    else
        ((flags & strip_indefinite_articles) > 0) && union!(skip_words, indefinite_articles(language))
        ((flags & strip_definite_articles) > 0) && union!(skip_words, definite_articles(language))
    end
    ((flags & strip_prepositions) > 0) && union!(skip_words, prepositions(language))
    ((flags & strip_pronouns) > 0) && union!(skip_words, pronouns(language))
    ((flags & strip_stopwords) > 0) && union!(skip_words, stopwords(language))
    # sparse, frequent terms
    ((flags & strip_sparse_terms) > 0) && union!(skip_words, sparse_terms(entity, alpha_sparse))
    ((flags & strip_frequent_terms) > 0) && union!(skip_words, frequent_terms(entity, alpha_frequent))
    if !isempty(skip_words)
        push!(rpatterns, _build_words_pattern(skip_words))
    end
    # custom regex
    if !isempty(skip_patterns)
        push!(rpatterns, _build_regex_pattern(skip_patterns))
    end
    # Do regex-based stripping
    if !isempty(rpatterns)
        r = _build_regex_pattern(rpatterns)
        remove_patterns!(entity, r)
    end
    # Stemming
    ((flags & stem_words) > 0) && stem!(entity)
    nothing
end

function prepare(s::AbstractString,
                 flags::UInt32;
                 language::Language = DEFAULT_LANGUAGE,
                 skip_patterns = Vector{Regex}(),
                 skip_words = Vector{String}(),
                 alpha_sparse=0.05,
                 alpha_frequent=0.95)
    os = s  # Initialize output string
    # Do function-based stripping
    ((flags & strip_corrupt_utf8) > 0) && (os = remove_corrupt_utf8(os))
    ((flags & strip_html_tags) > 0)    && (os = remove_html_tags(os))
    ((flags & strip_case) > 0)         && (os = remove_case(os))
    ((flags & strip_accents) > 0)      && (os = remove_accents(os))
    # regex
    rpatterns = Vector{Regex}(undef, 0)  # patterns to remove
    ((flags & strip_whitespace) > 0) && push!(rpatterns, strip2regex[strip_whitespace])
    if (flags & strip_non_ascii) > 0
        push!(rpatterns, strip2regex[strip_non_ascii])
    else
        ((flags & strip_numbers) > 0) && push!(rpatterns, strip2regex[strip_numbers])
        ((flags & strip_punctuation) > 0) && push!(rpatterns, strip2regex[strip_punctuation])
        ((flags & strip_single_chars) > 0) && push!(rpatterns, strip2regex[strip_single_chars])
    end
    # known words
    if (flags & strip_articles) > 0
        union!(skip_words, articles(language))
    else
        ((flags & strip_indefinite_articles) > 0) && union!(skip_words, indefinite_articles(language))
        ((flags & strip_definite_articles) > 0) && union!(skip_words, definite_articles(language))
    end
    ((flags & strip_prepositions) > 0) && union!(skip_words, prepositions(language))
    ((flags & strip_pronouns) > 0) && union!(skip_words, pronouns(language))
    ((flags & strip_stopwords) > 0) && union!(skip_words, stopwords(language))
    # sparse, frequent terms
    ((flags & strip_sparse_terms) > 0) && union!(skip_words, sparse_terms(os, alpha_sparse))
    ((flags & strip_frequent_terms) > 0) && union!(skip_words, frequent_terms(os, alpha_frequent))
    if !isempty(skip_words)
        push!(rpatterns, _build_words_pattern(skip_words))
    end
    # custom regex
    if !isempty(skip_patterns)
        push!(rpatterns, _build_regex_pattern(skip_patterns))
    end
    # Do regex-based stripping
    if !isempty(rpatterns)
        r = _build_regex_pattern(rpatterns)
        os = remove_patterns(os, r)
    end
    # Stemming
    ((flags & stem_words) > 0) && (os = stem(os, language=language))
    return os
end


text:
module Repl

const kliodir = joinpath(homedir(), "Klio")
const repldir = joinpath(kliodir, "repl.it")

using Pkg
Pkg.activate(kliodir)

include(joinpath(repldir, "install.jl"))

ENV["PATH"] = joinpath(homedir(), "opt/maxima/usr/bin") * ":" * ENV["PATH"]

include(joinpath(repldir, "run.jl"))

end


text:
using TuringModels, StatsFuns

Turing.setadbackend(:reversediff)

# outcome and predictor almost perfectly associated

x = repeat([-1], 9); append!(x, repeat([1],11))
y = repeat([0], 10); append!(y, repeat([1],10))

@model m_good_stan(x,y) = begin
    α ~ Normal(0,10)
    β ~ Normal(0,10)

    logits = α .+ β * x

    y .~ BinomialLogit.(1, logits)
end

chns = sample(m_good_stan(x,y), Turing.NUTS(0.65), 1000)

# Stan results

m_10_x,_results = "
    mean   sd   5.5% 94.5% n_eff Rhat
 a -5.09 4.08 -12.62 -0.25   100 1.00
 b  7.86 4.09   2.96 15.75   104 1.01
";

# Look at the proper draws (in corrected chn2)

describe(chns)

# End of 10/m10.xxt.jl


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

@doc raw"""
    ObjsReminder(;
        complete_ts=nothing,
        creator=nothing,
        id=nothing,
        recurring=nothing,
        text=nothing,
        time=nothing,
        user=nothing,
    )

    - complete_ts::Int32
    - creator::DefsUserId
    - id::DefsReminderId
    - recurring::Bool
    - text::String
    - time::Int32
    - user::DefsUserId
"""
mutable struct ObjsReminder <: SwaggerModel
    complete_ts::Any # spec type: Union{ Nothing, Int32 } # spec name: complete_ts
    creator::Any # spec type: Union{ Nothing, DefsUserId } # spec name: creator
    id::Any # spec type: Union{ Nothing, DefsReminderId } # spec name: id
    recurring::Any # spec type: Union{ Nothing, Bool } # spec name: recurring
    text::Any # spec type: Union{ Nothing, String } # spec name: text
    time::Any # spec type: Union{ Nothing, Int32 } # spec name: time
    user::Any # spec type: Union{ Nothing, DefsUserId } # spec name: user

    function ObjsReminder(;complete_ts=nothing, creator=nothing, id=nothing, recurring=nothing, text=nothing, time=nothing, user=nothing)
        o = new()
        validate_property(ObjsReminder, Symbol("complete_ts"), complete_ts)
        setfield!(o, Symbol("complete_ts"), complete_ts)
        validate_property(ObjsReminder, Symbol("creator"), creator)
        setfield!(o, Symbol("creator"), creator)
        validate_property(ObjsReminder, Symbol("id"), id)
        setfield!(o, Symbol("id"), id)
        validate_property(ObjsReminder, Symbol("recurring"), recurring)
        setfield!(o, Symbol("recurring"), recurring)
        validate_property(ObjsReminder, Symbol("text"), text)
        setfield!(o, Symbol("text"), text)
        validate_property(ObjsReminder, Symbol("time"), time)
        setfield!(o, Symbol("time"), time)
        validate_property(ObjsReminder, Symbol("user"), user)
        setfield!(o, Symbol("user"), user)
        o
    end
end # type ObjsReminder

const _property_map_ObjsReminder = Dict{Symbol,Symbol}(Symbol("complete_ts")=>Symbol("complete_ts"), Symbol("creator")=>Symbol("creator"), Symbol("id")=>Symbol("id"), Symbol("recurring")=>Symbol("recurring"), Symbol("text")=>Symbol("text"), Symbol("time")=>Symbol("time"), Symbol("user")=>Symbol("user"))
const _property_types_ObjsReminder = Dict{Symbol,String}(Symbol("complete_ts")=>"Int32", Symbol("creator")=>"DefsUserId", Symbol("id")=>"DefsReminderId", Symbol("recurring")=>"Bool", Symbol("text")=>"String", Symbol("time")=>"Int32", Symbol("user")=>"DefsUserId")
Base.propertynames(::Type{ ObjsReminder }) = collect(keys(_property_map_ObjsReminder))
Swagger.property_type(::Type{ ObjsReminder }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_ObjsReminder[name]))}
Swagger.field_name(::Type{ ObjsReminder }, property_name::Symbol) =  _property_map_ObjsReminder[property_name]

function check_required(o::ObjsReminder)
    (getproperty(o, Symbol("creator")) === nothing) && (return false)
    (getproperty(o, Symbol("id")) === nothing) && (return false)
    (getproperty(o, Symbol("recurring")) === nothing) && (return false)
    (getproperty(o, Symbol("text")) === nothing) && (return false)
    (getproperty(o, Symbol("user")) === nothing) && (return false)
    true
end

function validate_property(::Type{ ObjsReminder }, name::Symbol, val)
end


text:
function meanColumn(A)
  return mean(A,1)
end

function stdColumn(A)
  return std(A,1)
end

function sumColumn(A)
  return sum(A,1)
end

function getMinIndices(A,K=1)
  return sortperm(A)[1:K]
end

function getMajority(A)
  B=Dict()
  for i in A
    if haskey(B,i)
      B[i]+=1
    else
      B[i]=1
    end

  end

  max= maximum(values(B))
  # println(B)
  # println(max)
  filteredresults = filter(((k,v)->v==max),B)
  indices = collect(keys(filteredresults))
  return indices[rand(1:length(indices))]
end

max_stars_repo_name: ajermaky/CustomMLLibrary

max_issues_repo_name: ajermaky/CustomMLLibrary

max_forks_repo_name: ajermaky/CustomMLLibrary


text:
is_nothing(cur_object) = false
is_nothing(cur_object::Void) = true

is_present(cur_object) = !is_nothing(cur_object)


text:
struct ModelSetup{T<:Real,Tprog<:Real}
    parameters::Parameter
    grid::Grid{T,Tprog}
    constants::Constants{T,Tprog}
    forcing::Forcing{T}
end


text:
#===============================================================================

    optical_test_dual_channel.jl

===============================================================================#

import DataStructures.OrderedDict
import StaticArrays.SVector
import Memento: debug, warn

#===============================================================================
    constants >>
===============================================================================#

## channel descriptors
const CHANNELS = DYES = SVector(1,2)
const DYE_SYMBOLS = SVector(:FAM, :HEX)
const OPTICAL_TEST_SYMBOLS = SVector(:baseline, :water, DYE_SYMBOLS...)

## preset values >>

## bounds of signal-to-noise ratio (SNR)
const SNR_FAM_CH1_MIN = 0.75
const SNR_FAM_CH2_MAX = 1
const SNR_HEX_CH1_MAX = 0.50
const SNR_HEX_CH2_MIN = 0.88

## signal-to-noise ratio discriminant functions for each well
dscrmnt_snr_fam(snr_2chs) = [snr_2chs[1] > SNR_FAM_CH1_MIN, snr_2chs[2] < SNR_FAM_CH2_MAX]
dscrmnt_snr_hex(snr_2chs) = [snr_2chs[1] < SNR_HEX_CH1_MAX, snr_2chs[2] > SNR_HEX_CH2_MIN]
const dscrmnts_snr = OrderedDict(map(1:2) do i
    DYE_SYMBOLS[i] => [dscrmnt_snr_fam, dscrmnt_snr_hex][i]
end) ## do i

## maximum and minimum fluorescence values: channel 1, channel 2
const WATER_MAX = [32000, 5000]
const WATER_MIN = [1000, -1000]

#===============================================================================
    function definition >>
===============================================================================#

## called by dispatch()
function act(
    ::Type{Val{optical_test_dual_channel}},
    req             ::Associative; ## keys: "baseline", "water", "FAM", "HEX"
    out_format      ::OutputFormat = pre_json_output
    ## remove MySqldependency
    #
    # db_conn ::MySQL.MySQLHandle,
    # exp_id ::Integer,
    # calib_info ::Union{Integer,OrderedDict}; ## keys: "baseline", "water", "channel_1", "channel_2". Each value's "calibration_id" value is the same as `exp_id`
    #
    # start: arguments that might be passed by upstream code
    # well_nums ::AbstractVector =[],
)
    function SNR_test(w, label)
        (baseline_2chs, water_fluo_2chs, signal_fluo_2chs) =
            map(key -> fluo_dict[key][w, :], [:baseline, :water, label])
        snr_2chs =
            (signal_fluo_2chs .- water_fluo_2chs) ./ (signal_fluo_2chs .- baseline_2chs)
        transpose(dscrmnts_snr[label](transpose(snr_2chs)))
    end

    debug(logger, "at act(::Type{Val{optical_test_dual_channel}})")

    ## remove MySql dependency
    #
    # fluo_qry_2b = "SELECT step_id, well_num, fluorescence_value, channel
    #     FROM fluorescence_data
    #     WHERE experiment_id = $exp_id AND cycle_num = 1 AND step_id is not NULL
    #     well_constraint
    #     ORDER BY well_num, channel
    # "
    # fluo_data, fluo_well_nums = get_mysql_data_well(
    #     well_nums, fluo_qry_2b, db_conn, false
    # )
    #
    # num_wells = length(fluo_well_nums)

    # fluo_dict = OrderedDict(map(old_labels) do label ## old
    fluo_dict =
        OrderedDict(
            map(OPTICAL_TEST_SYMBOLS) do label
                label =>
                    mapreduce(
                        ## old
                        # fluo_data[:fluorescence_value][
                        #     (fluo_data[:step_id] .== calib_info[label]["step_id"]) .& (fluo_data[:channel] .== channel)
                        # ]
                        channel -> req[string(label)][FLUORESCENCE_VALUE_KEY][channel],
                        hcat,
                        CHANNELS)
            end)
    num_wells = size(fluo_dict[:baseline], 1)
    #
    bool_dict = OrderedDict(:baseline => fill(true, num_wells, length(CHANNELS)))
    ## water test
    bool_dict[:water] =
        mapreduce(
            channel ->
                map(fluo_dict[:water][:, channel]) do fluo_pw
                    WATER_MIN[channel] < fluo_pw < WATER_MAX[channel]
                end,
            hcat,
            CHANNELS)
    ## FAM and HEX SNR test
    for label in DYE_SYMBOLS
        bool_dict[label] =
            mapreduce(
                well -> SNR_test(well, label),
                vcat,
                1:num_wells)
    end
    #
    ## organize "optical_data"
    optical_data =
        map(1:num_wells) do well
            OrderedDict(
                map(OPTICAL_TEST_SYMBOLS) do label
                    label =>
                        map(CHANNELS) do channel
                            (fluo_dict[label][well, channel],
                                bool_dict[label][well, channel])
                        end ## do channel

                end) ## do label
        end ## do well
    #
    ## FAM and HEX self-calibrated
    ## ((signal_of_dye_x_in_channel_k - water_in_channel_k) /
    ##     (signal_of_target_dye_in_channel_k - water_in_channel_k); x=FAM,HEX; k=1,2)
    ## ratio (self_calib_of_dye_x_where_k_equals_1 / self_calib_of_dye_x_where_k_equals_2) test.
    #
    ## Note:
    ## Baseline is not subtracted because it is not part of the calibration procedure.
    #
    ## Issues:
    ## Is it more sensible to calculate target:off-target ratio than channel_1:channel_2?
    ## Calculation of ratio may fail if water >= signal values,
    ## reporting negative or infinite values
    #
    ## substract water values from signal values
    swd_vec =
        map(DYE_SYMBOLS) do label
            map(CHANNELS) do channel
                fluo_dict[label][:, channel] .- fluo_dict[:water][:, channel]
            end ## do channel
        end ## do label
    ## calculate normalization values from data in target channels
    swd_normd =
        map(CHANNELS) do target
            sweep(mean)(/)(swd_vec[target][target]) ## dye == channel
        end
    ## normalize signal data
    self_calib_vec =
        map(swd_vec) do swd_dye
            map(CHANNELS) do channel
                swd_dye[channel] ./ swd_normd[channel]
            end ## do channel
        end ## do swd_dye
    #
    ## call as invalid analysis if there are negative or zero values in the normalized data
    ## that will cause the channel1:channel2 ratio to be zero, infinite, or negative
    ## devectorized
    error_msgs = Vector{String}()
    for dye in DYES, channel in CHANNELS, value in self_calib_vec[dye][channel]
        if value <= 0.0
            ## call as invalid analysis instead of raising an error
            push!(error_msgs, "zero or negative values in the self-calibrated fluorescence data")
            warn(logger, error_msgs[end])
            break ## exit nested loops
        end ## if
    end ## for dye, channel, value
    #
    ## calculate channel1:channel2 ratios
    ch12_ratios =
        OrderedDict(
            map(DYES) do dye
                sc_dye = self_calib_vec[dye]
                DYE_SYMBOLS[dye] => round.(sc_dye[1] ./ sc_dye[2], JSON_DIGITS)
            end) # do channel_i
    if !(ch12_ratios |> values |> mold(x -> all(isfinite.(x))) |> all) ||
        (ch12_ratios |> values |> mold(x -> any(x .<= 0))      |> any)
            push!(error_msgs, "zero, negative, or infinite values of channel 1 : channel 2 ratio")
            warn(logger, error_msgs[end])
    end

    ## return values
    output = OrderedDict(
        :optical_data       => optical_data,
        Symbol("Ch1:Ch2")   => ch12_ratios,
        :valid              => length(error_msgs) == 0,
        :error              => join(error_msgs, "; "))
    return output |> out(out_format)
end ## act(::Type{Val{optical_test_dual_channel}})

max_stars_repo_name: MakerButt/chaipcr

max_issues_repo_name: MakerButt/chaipcr

max_forks_repo_name: MakerButt/chaipcr


text:
using BranchedFlow
using Test

@testset "BranchedFlow.jl" begin
    # Write your tests here.
end


text:
###############################################################################
#
# Initialize SPACMono structure
#
###############################################################################
"""
    initialize_spac_canopy!(node::SPACMono{FT}) where {FT<:AbstractFloat}

Initialize the RT parameters for a given
- `node` [`SPACMono`](@ref) type struct
"""
function initialize_spac_canopy!(
            node::SPACMono{FT}
) where {FT<:AbstractFloat}
    # 0.1 create variables required
    @unpack angles, can_opt, can_rad, canopy_rt, envirs, in_rad, leaves_rt,
            n_canopy, plant_hs, plant_ps, photo_set, rt_con, soil_opt,
            wl_set = node;
    fraction_sl::Array{FT,1} = repeat(canopy_rt.lidf, outer=[canopy_rt.nAzi]) /
                               length(canopy_rt.lazitab);
    n_sl = length(canopy_rt.lidf) * length(canopy_rt.lazitab);

    # fluspect the canopy layers
    for ican in 1:n_canopy
        fluspect!(leaves_rt[ican], wl_set);
    end

    # Four Different steps to compute Short-Wave RT
    canopy_geometry!(canopy_rt, angles, can_opt, rt_con)
    canopy_matrices!(leaves_rt, can_opt);
    short_wave!(canopy_rt, can_opt, can_rad, in_rad,
                soil_opt, rt_con);
    canopy_fluxes!(canopy_rt, can_opt, can_rad, in_rad,
                   soil_opt, leaves_rt, wl_set, rt_con);

    # Compute Long Wave (Last term is LW incoming in W m^-2)
    thermal_fluxes!(leaves_rt, can_opt, can_rad, canopy_rt,
                    soil_opt, [FT(400.0)], wl_set);

    # update the canopy leaf area partition information
    for i_can in 1:n_canopy
        rt_layer = n_canopy + 1 - i_can;
        iPS      = plant_ps[i_can];

        # calculate the fraction of sunlit and shaded leaves
        f_view = (can_opt.Ps[rt_layer] + can_opt.Ps[rt_layer+1]) / 2;

        for iLF in 1:n_sl
            iPS.APAR[iLF] = can_rad.absPAR_sunCab[(rt_layer-1)*n_sl + iLF] *
                            FT(1e6);
            iPS.LAIx[iLF] = f_view * fraction_sl[iLF];
        end
        iPS.APAR[end] = can_rad.absPAR_shadeCab[rt_layer] * FT(1e6);
        iPS.LAIx[end] = 1 - f_view;

        update_leaf_TP!(photo_set, iPS, plant_hs.leaves[i_can], envirs[i_can]);

        envirs[i_can].t_air = K_25(FT);
        envirs[i_can].p_sat = saturation_vapor_pressure(K_25(FT));
        envirs[i_can].p_H₂O = envirs[i_can].p_sat / 2;
    end

    return nothing
end

max_stars_repo_name: bbuman/Land

max_issues_repo_name: bbuman/Land

max_forks_repo_name: bbuman/Land