diff --git a/Project.toml b/Project.toml
index c76d0865..aefc030d 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,7 +1,7 @@
 name = "InformationGeometry"
 uuid = "ee9d1287-bbdf-432c-9920-c447cf97a828"
 authors = ["Rafael Arutjunjan"]
-version = "1.13.3"
+version = "1.14.0"
 
 [deps]
 BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf"
diff --git a/src/ConfidenceRegions.jl b/src/ConfidenceRegions.jl
index 56c4db49..5dfb5a14 100644
--- a/src/ConfidenceRegions.jl
+++ b/src/ConfidenceRegions.jl
@@ -1,78 +1,5 @@
 
 
-"""
-    likelihood(DM::DataModel, θ::AbstractVector) -> Real
-Calculates the likelihood ``L(\\mathrm{data} \\, | \\, \\theta)`` a `DataModel` and a parameter configuration ``\\theta``.
-"""
-likelihood(args...; kwargs...) = exp(loglikelihood(args...; kwargs...))
-
-
-
-## Prefix underscore for likelihood, Score and FisherMetric indicates that Prior has already been accounted for upstream
-loglikelihood(DM::AbstractDataModel; kwargs...) = LogLikelihood(θ::AbstractVector{<:Number}; Kwargs...) = loglikelihood(DM, θ; kwargs..., Kwargs...)
-Negloglikelihood(DM::AbstractDataModel; kwargs...) = NegativeLogLikelihood(θ::AbstractVector{<:Number}; Kwargs...) = -loglikelihood(DM, θ; kwargs..., Kwargs...)
-
-# import Distributions.loglikelihood
-"""
-    loglikelihood(DM::DataModel, θ::AbstractVector) -> Real
-Calculates the logarithm of the likelihood ``L``, i.e. ``\\ell(\\mathrm{data} \\, | \\, \\theta) \\coloneqq \\mathrm{ln} \\big( L(\\mathrm{data} \\, | \\, \\theta) \\big)`` given a `DataModel` and a parameter configuration ``\\theta``.
-"""
-loglikelihood(DM::AbstractDataModel, θ::AbstractVector{<:Number}, LogPriorFn::Union{Nothing,Function}=LogPrior(DM); kwargs...) = loglikelihood(Data(DM), Predictor(DM), θ, LogPriorFn; kwargs...)
-
-loglikelihood(DS::AbstractDataSet, model::ModelOrFunction, θ::AbstractVector{<:Number}, LogPriorFn::Nothing; kwargs...) = _loglikelihood(DS, model, θ; kwargs...)
-loglikelihood(DS::AbstractDataSet, model::ModelOrFunction, θ::AbstractVector{<:Number}, LogPriorFn::Function; kwargs...) = _loglikelihood(DS, model, θ; kwargs...) + EvalLogPrior(LogPriorFn, θ)
-
-
-# Specialize this for different DataSet types
-@inline function _loglikelihood(DS::AbstractDataSet, model::ModelOrFunction, θ::AbstractVector{<:Number}; kwargs...)
-    -0.5*(DataspaceDim(DS)*log(2π) - logdetInvCov(DS) + InnerProduct(yInvCov(DS), ydata(DS)-EmbeddingMap(DS, model, θ; kwargs...)))
-end
-
-
-
-AutoScore(DM::AbstractDataModel, θ::AbstractVector{<:Number}, LogPriorFn::Union{Nothing, Function}=LogPrior(DM); kwargs...) = AutoScore(Data(DM), Predictor(DM), θ, LogPriorFn; kwargs...)
-AutoScore(DS::AbstractDataSet, model::ModelOrFunction, θ::AbstractVector{<:Number}, LogPriorFn::Nothing; kwargs...) = _AutoScore(DS, model, θ; kwargs...)
-function AutoScore(DS::AbstractDataSet, model::ModelOrFunction, θ::AbstractVector{<:Number}, LogPriorFn::Function; ADmode::Union{Symbol,Val}=Val(:ForwardDiff), kwargs...)
-    _AutoScore(DS, model, θ; ADmode=ADmode, kwargs...) + EvalLogPriorGrad(LogPriorFn, θ; ADmode=ADmode)
-end
-_AutoScore(DS::AbstractDataSet, model::ModelOrFunction, θ::AbstractVector{<:Number}; ADmode::Union{Symbol,Val}=Val(:ForwardDiff), kwargs...) = GetGrad(ADmode, x->_loglikelihood(DS, model, x; kwargs...))(θ)
-
-
-AutoMetric(DM::AbstractDataModel, θ::AbstractVector{<:Number}, LogPriorFn::Union{Nothing, Function}=LogPrior(DM); kwargs...) = AutoMetric(Data(DM), Predictor(DM), θ, LogPriorFn; kwargs...)
-AutoMetric(DS::AbstractDataSet, model::ModelOrFunction, θ::AbstractVector{<:Number}, LogPriorFn::Nothing; kwargs...) = _AutoMetric(DS, model, θ; kwargs...)
-function AutoMetric(DS::AbstractDataSet, model::ModelOrFunction, θ::AbstractVector{<:Number}, LogPriorFn::Function; ADmode::Union{Symbol,Val}=Val(:ForwardDiff), kwargs...)
-    _AutoMetric(DS, model, θ; ADmode=ADmode, kwargs...) - EvalLogPriorHess(LogPriorFn, θ; ADmode=ADmode)
-end
-_AutoMetric(DS::AbstractDataSet, model::ModelOrFunction, θ::AbstractVector{<:Number}; ADmode::Union{Symbol,Val}=Val(:ForwardDiff), kwargs...) = GetHess(ADmode, x->-_loglikelihood(DS, model, x; kwargs...))(θ)
-
-
-
-Score(DM::AbstractDataModel; kwargs...) = LogLikelihoodGradient(θ::AbstractVector{<:Number}; Kwargs...) = Score(DM, θ; kwargs..., Kwargs...)
-NegScore(DM::AbstractDataModel; kwargs...) = NegativeLogLikelihoodGradient(θ::AbstractVector{<:Number}; Kwargs...) = -Score(DM, θ; kwargs..., Kwargs...)
-
-"""
-    Score(DM::DataModel, θ::AbstractVector{<:Number}; ADmode::Val=Val(:ForwardDiff))
-Calculates the gradient of the log-likelihood ``\\ell`` with respect to a set of parameters ``\\theta``.
-`ADmode=Val(false)` computes the Score by separately evaluating the `model` as well as the Jacobian `dmodel` provided in `DM`.
-Other choices of `ADmode` directly compute the Score by differentiating the formula the log-likelihood, i.e. only one evaluation on a dual variable is performed.
-"""
-Score(DM::AbstractDataModel, θ::AbstractVector{<:Number}, LogPriorFn::Union{Nothing,Function}=LogPrior(DM); kwargs...) = Score(Data(DM), Predictor(DM), dPredictor(DM), θ, LogPriorFn; kwargs...)
-
-Score(DS::AbstractDataSet, model::ModelOrFunction, dmodel::ModelOrFunction, θ::AbstractVector{<:Number}, LogPriorFn::Nothing; kwargs...) = _Score(DS, model, dmodel, θ; kwargs...)
-Score(DS::AbstractDataSet, model::ModelOrFunction, dmodel::ModelOrFunction, θ::AbstractVector{<:Number}, LogPriorFn::Function; kwargs...) = _Score(DS, model, dmodel, θ; kwargs...) + EvalLogPriorGrad(LogPriorFn, θ)
-
-
-
-_Score(DS::AbstractDataSet, model::ModelOrFunction, dmodel::ModelOrFunction, θ::AbstractVector{<:Number}; ADmode::Val=Val(:ForwardDiff), kwargs...) = _Score(DS, model, dmodel, θ, ADmode; kwargs...)
-# Delegate to AutoScore
-_Score(DS::AbstractDataSet, model::ModelOrFunction, dmodel::ModelOrFunction, θ::AbstractVector{<:Number}, ADmode::Val{true}; kwargs...) = _AutoScore(DS, model, θ; ADmode=ADmode, kwargs...)
-_Score(DS::AbstractDataSet, model::ModelOrFunction, dmodel::ModelOrFunction, θ::AbstractVector{<:Number}, ADmode::Union{Symbol,Val}; kwargs...) = _AutoScore(DS, model, θ; ADmode=ADmode, kwargs...)
-
-# Specialize this for different DataSet types
-@inline function _Score(DS::AbstractDataSet, model::ModelOrFunction, dmodel::ModelOrFunction, θ::AbstractVector{<:Number}, ADmode::Val{false}; kwargs...)
-    transpose(EmbeddingMatrix(DS,dmodel,θ; kwargs...)) * (yInvCov(DS) * (ydata(DS) - EmbeddingMap(DS,model,θ; kwargs...)))
-end
-
 
 """
 Point θ lies outside confidence region of level `Confvol` if this function > 0.
@@ -783,108 +710,6 @@ function Sensitivity(DM::AbstractDataModel, Confnum::Real=1; Approx::Bool=false,
 end
 
 
-### for-loop typically slower than reduce(vcat, ...)
-### Apparently curve_fit() throws an error in conjuction with ForwardDiff when reinterpret() is used
-# Reduction(X::AbstractVector{<:SVector{Len,T}}) where Len where T = reinterpret(T, X)
-Reduction(X::AbstractVector{<:AbstractVector}) = reduce(vcat, X)
-Reduction(X::AbstractVector{<:Number}) = X
-Reduction(X::AbstractVector{<:SubArray{<:Number, 0}}) = [@inbounds X[i][1] for i in 1:length(X)]
-
-
-# h(θ) ∈ Dataspace
-"""
-    EmbeddingMap(DM::AbstractDataModel, θ::AbstractVector{<:Number}) -> Vector
-Returns a vector of the collective predictions of the `model` as evaluated at the x-values and the parameter configuration ``\\theta``.
-```
-h(\\theta) \\coloneqq \\big(y_\\mathrm{model}(x_1;\\theta),...,y_\\mathrm{model}(x_N;\\theta)\\big) \\in \\mathcal{D}
-```
-"""
-EmbeddingMap(DM::AbstractDataModel, θ::AbstractVector{<:Number}, woundX::AbstractVector=WoundX(DM); kwargs...) = EmbeddingMap(Data(DM), Predictor(DM), θ, woundX; kwargs...)
-EmbeddingMap(DS::AbstractDataSet, model::ModelOrFunction, θ::AbstractVector{<:Number}, woundX::AbstractVector=WoundX(DS); kwargs...) = _CustomOrNot(DS, model, θ, woundX; kwargs...)
-EmbeddingMap(DS::Val, model::ModelOrFunction, θ::AbstractVector{<:Number}, woundX::AbstractVector; kwargs...) = _CustomOrNot(DS, model, θ, woundX; kwargs...)
-
-_CustomOrNot(DS::Union{Val,AbstractDataSet}, model::Function, θ::AbstractVector{<:Number}, woundX::AbstractVector; kwargs...) = _CustomOrNot(DS, model, θ, woundX, Val(false), Val(false); kwargs...)
-_CustomOrNot(DS::Union{Val,AbstractDataSet}, M::ModelMap, θ::AbstractVector{<:Number}, woundX::AbstractVector; kwargs...) = _CustomOrNot(DS, M.Map, θ, woundX, M.CustomEmbedding, M.inplace; kwargs...)
-
-# Specialize this for different Dataset types
-_CustomOrNot(::Union{Val,AbstractDataSet}, model::ModelOrFunction, θ::AbstractVector{<:Number}, woundX::AbstractVector, custom::Val{false}, inplace::Val{false}; kwargs...) = Reduction(map(x->model(x,θ; kwargs...), woundX))
-_CustomOrNot(::Union{Val,AbstractDataSet}, model::ModelOrFunction, θ::AbstractVector{<:Number}, woundX::AbstractVector, custom::Val{true}, inplace::Val{false}; kwargs...) = model(woundX, θ; kwargs...)
-
-
-function _CustomOrNot(DS::AbstractDataSet, model!::ModelOrFunction, θ::AbstractVector{<:Number}, woundX::AbstractVector, custom::Val{true}, inplace::Val{true}; kwargs...)
-    Y = Vector{suff(θ)}(undef, length(woundX)*ydim(DS))
-    model!(Y, woundX, θ; kwargs...);    Y
-end
-function _CustomOrNot(DS::AbstractDataSet, model!::ModelOrFunction, θ::AbstractVector{<:Number}, woundX::AbstractVector, custom::Val{false}, inplace::Val{true}; kwargs...)
-    Y = Vector{suff(θ)}(undef, length(woundX)*ydim(DS))
-    EmbeddingMap!(Y, DS, model!, θ, woundX; kwargs...);     Y
-end
-
-
-function EmbeddingMap!(Y::AbstractVector{<:Number}, DS::AbstractDataSet, model!::ModelOrFunction, θ::AbstractVector{<:Number}, woundX::AbstractVector=WoundX(DS); kwargs...)
-    EmbeddingMap!(Y, model!, θ, woundX, Val(ydim(DS)); kwargs...)
-end
-function EmbeddingMap!(Y::AbstractVector{<:Number}, model!::ModelOrFunction, θ::AbstractVector{<:Number}, woundX::AbstractVector, Ydim::Val{1}; kwargs...)
-    @inbounds for i in Base.OneTo(length(Y))
-        model!(Y[i], woundX[i], θ; kwargs...)
-    end
-end
-function EmbeddingMap!(Y::AbstractVector{<:Number}, model!::ModelOrFunction, θ::AbstractVector{<:Number}, woundX::AbstractVector, Ydim::Val{T}; kwargs...) where T
-    @inbounds for (i, row) in enumerate(Iterators.partition(1:length(Y), T))
-        model!(view(Y,row), woundX[i], θ; kwargs...)
-    end
-end
-# Fallback for out-of-place models
-function EmbeddingMap!(Y::AbstractVector{<:Number}, DS::AbstractDataSet, model!::ModelMap{false}, θ::AbstractVector{<:Number}, woundX::AbstractVector=WoundX(DS); kwargs...)
-    copyto!(Y, EmbeddingMap(DS, model!, θ, woundX; kwargs...))
-end
-
-
-
-"""
-    EmbeddingMatrix(DM::AbstractDataModel, θ::AbstractVector{<:Number}) -> Matrix
-Returns the jacobian of the embedding map as evaluated at the x-values and the parameter configuration ``\\theta``.
-"""
-EmbeddingMatrix(DM::AbstractDataModel, θ::AbstractVector{<:Number}, woundX::AbstractVector=WoundX(DM); kwargs...) = EmbeddingMatrix(Data(DM), dPredictor(DM), θ, woundX; kwargs...)
-EmbeddingMatrix(DS::AbstractDataSet, dmodel::ModelOrFunction, θ::AbstractVector{<:Number}, woundX::AbstractVector=WoundX(DS); kwargs...) = _CustomOrNotdM(DS, dmodel, θ, woundX; kwargs...)
-EmbeddingMatrix(DS::Val, dmodel::ModelOrFunction, θ::AbstractVector{<:Number}, woundX::AbstractVector; kwargs...) = _CustomOrNotdM(DS, dmodel, θ, woundX; kwargs...)
-
-_CustomOrNotdM(DS::Union{Val,AbstractDataSet}, dmodel::Function, θ::AbstractVector{<:Number}, woundX::AbstractVector; kwargs...) = _CustomOrNotdM(DS, dmodel, floatify(θ), woundX, Val(false), Val(false); kwargs...)
-_CustomOrNotdM(DS::Union{Val,AbstractDataSet}, dM::ModelMap, θ::AbstractVector{<:Number}, woundX::AbstractVector; kwargs...) = _CustomOrNotdM(DS, dM.Map, floatify(θ), woundX, dM.CustomEmbedding, dM.inplace; kwargs...)
-
-# Specialize this for different Dataset types
-_CustomOrNotdM(::Union{Val,AbstractDataSet}, dmodel::ModelOrFunction, θ::AbstractVector{<:Number}, woundX::AbstractVector, custom::Val{false}, inplace::Val{false}; kwargs...) = reduce(vcat, map(x->dmodel(x,θ; kwargs...), woundX))
-_CustomOrNotdM(::Union{Val,AbstractDataSet}, dmodel::ModelOrFunction, θ::AbstractVector{<:Number}, woundX::AbstractVector, custom::Val{true}, inplace::Val{false}; kwargs...) = dmodel(woundX, θ; kwargs...)
-
-
-function _CustomOrNotdM(DS::AbstractDataSet, dmodel!::ModelOrFunction, θ::AbstractVector{<:Number}, woundX::AbstractVector, custom::Val{true}, inplace::Val{true}; kwargs...)
-    J = Matrix{suff(θ)}(undef, length(woundX)*ydim(DS), length(θ))
-    dmodel!(J, woundX, θ; kwargs...);   J
-end
-function _CustomOrNotdM(DS::AbstractDataSet, dmodel!::ModelOrFunction, θ::AbstractVector{<:Number}, woundX::AbstractVector, custom::Val{false}, inplace::Val{true}; kwargs...)
-    J = Matrix{suff(θ)}(undef, length(woundX)*ydim(DS), length(θ))
-    EmbeddingMatrix!(J, DS, dmodel!, θ, woundX; kwargs...);     J
-end
-
-
-function EmbeddingMatrix!(J::AbstractMatrix{<:Number}, DS::AbstractDataSet, dmodel!::ModelOrFunction, θ::AbstractVector{<:Number}, woundX::AbstractVector=WoundX(DS); kwargs...)
-    EmbeddingMatrix!(J, dmodel!, θ, woundX, Val(ydim(DS)); kwargs...)
-end
-function EmbeddingMatrix!(J::AbstractMatrix{<:Number}, dmodel!::ModelOrFunction, θ::AbstractVector{<:Number}, woundX::AbstractVector, Ydim::Val{1}; kwargs...)
-    @inbounds for row in Base.OneTo(size(J,1))
-        dmodel!(view(J,row,:), woundX[row], θ; kwargs...)
-    end
-end
-function EmbeddingMatrix!(J::AbstractMatrix{<:Number}, dmodel!::ModelOrFunction, θ::AbstractVector{<:Number}, woundX::AbstractVector, Ydim::Val{T}; kwargs...) where T
-    @inbounds for (i, row) in enumerate(Iterators.partition(1:size(J,1), T))
-        dmodel!(view(J,row,:), woundX[i], θ; kwargs...)
-    end
-end
-# Fallback for out-of-place models
-function EmbeddingMatrix!(J::AbstractMatrix{<:Number}, DS::AbstractDataSet, dmodel!::ModelMap{false}, θ::AbstractVector{<:Number}, woundX::AbstractVector=WoundX(DS); kwargs...)
-    copyto!(J, EmbeddingMatrix(DS, dmodel!, θ, woundX; kwargs...))
-end
-
 
 
 # M ⟵ D
@@ -1283,43 +1108,6 @@ end
 
 
 
-"""
-    LiftedLogLikelihood(DM::AbstractDataModel) -> Function
-Computes log-likelihood on the extended data space ``\\hat{\\ell} : \\mathcal{X}^N \\times\\mathcal{Y}^N \\longrightarrow \\mathbb{R}``.
-Should be maximized. Does not account for priors.
-"""
-LiftedLogLikelihood(DM::Union{AbstractDataModel,AbstractDataSet}) = (G = dist(DM);  ℓ(Z::AbstractVector{<:Number}) = logpdf(G, Z))
-"""
-    LiftedCost(DM::AbstractDataModel) -> Function
-Computes negative log-likelihood as cost function on the extended data space ``C : \\mathcal{X}^N \\times\\mathcal{Y}^N \\longrightarrow \\mathbb{R}``.
-Should be minimized. Does not account for priors.
-"""
-LiftedCost(DM::Union{AbstractDataModel,AbstractDataSet}) = (G = dist(DM);  Negativeℓ(Z::AbstractVector{<:Number}) = -logpdf(G, Z))
-
-"""
-    LiftedEmbedding(DM::AbstractDataModel) -> Function
-Constructs lifted embedding map from initial space into extended dataspace ``\\hat{h} : \\mathcal{X}^N \\times \\mathcal{M} \\longrightarrow \\mathcal{X}^N \\times\\mathcal{Y}^N`` effecting
-``\\xi = (x_\\text{opt}, \\theta) \\longmapsto \\hat{h}(\\xi) = (x_\\text{opt}, h(x_\\text{opt}, \\theta))``.
-"""
-LiftedEmbedding(DM::AbstractDataModel) = LiftedEmbedding(Data(DM), Predictor(DM), pdim(DM))
-function LiftedEmbedding(DS::AbstractDataSet, Model::ModelOrFunction, pd::Int)
-    ĥ(ξ::AbstractVector; kwargs...) = ĥ(view(ξ,1:length(ξ)-pd), view(ξ,length(ξ)-pd+1:length(ξ)); kwargs...)
-    ĥ(xdat::AbstractVector, θ::AbstractVector{<:Number}; kwargs...) = [xdat; EmbeddingMap(DS, Model, θ, Windup(xdat, xdim(DS)); kwargs...)]
-end
-
-"""
-    FullLiftedLogLikelihood(DM::AbstractDataModel) -> Function
-Computes the full likelihood given Xθ INCLUDING PRIOR.
-"""
-FullLiftedLogLikelihood(DM::AbstractDataModel) = FullLiftedLogLikelihood(Data(DM), Predictor(DM), LogPrior(DM), pdim(DM))
-FullLiftedLogLikelihood(DS::AbstractDataSet, model::ModelOrFunction, LogPriorFn::Nothing, pd::Int) = LiftedLogLikelihood(DS)∘LiftedEmbedding(DS, model, pd)
-function FullLiftedLogLikelihood(DS::AbstractDataSet, model::ModelOrFunction, LogPriorFn::Function, pd::Int)
-    L = LiftedLogLikelihood(DS)∘LiftedEmbedding(DS, model, pd)
-    ℓ(Xθ::AbstractVector{<:Number}; kwargs...) = L(Xθ; kwargs...) + EvalLogPrior(LogPriorFn, view(Xθ, length(Xθ)-pd+1:length(Xθ)))
-end
-FullLiftedNegLogLikelihood(args...; kwargs...) = (L=FullLiftedLogLikelihood(args...; kwargs...); Xθ::AbstractVector{<:Number}->-L(Xθ))
-
-
 abstract type AbstractBoundarySlice end
 
 struct ConfidenceBoundarySlice <: AbstractBoundarySlice
diff --git a/src/InformationGeometry.jl b/src/InformationGeometry.jl
index a6e12b1e..0eaadf05 100644
--- a/src/InformationGeometry.jl
+++ b/src/InformationGeometry.jl
@@ -181,8 +181,16 @@ export LineSearch, MonteCarloArea
 export curve_fit, RobustFit, TotalLeastSquares, BlockMatrix
 
 
+include("Likelihoods.jl")
+export likelihood, loglikelihood, Score
+
+
+include("ModelPredictions.jl")
+export EmbeddingMap, EmbeddingMatrix, EmbeddingMap!, EmbeddingMatrix!
+
+
 include("ConfidenceRegions.jl")
-export likelihood, loglikelihood, Score, WilksCriterion, WilksTest, OrthVF, OrthVF!, FindMLE
+export WilksCriterion, WilksTest, OrthVF, OrthVF!, FindMLE
 export AutoScore, AutoMetric
 export FindConfBoundary, FCriterion, FTest, FindFBoundary
 export GenerateBoundary, ConfidenceRegion, ConfidenceRegions
@@ -192,7 +200,7 @@ export FisherMetric, GeometricDensity
 export ConfidenceRegionVolume, CoordinateVolume
 export ExpectedInvariantVolume, GeodesicRadius, CoordinateDistortion, Sensitivity
 
-export EmbeddingMap, EmbeddingMatrix, EmbeddingMap!, EmbeddingMatrix!, Pullback, Pushforward
+export Pullback, Pushforward
 export AIC, AICc, BIC, ModelComparison, IsLinearParameter, IsLinear, LeastInformativeDirection
 
 export FindConfBoundaryOnPlane, LinearCuboid, IntersectCube, IntersectRegion, MincedBoundaries, ConfidenceBoundary
diff --git a/src/Likelihoods.jl b/src/Likelihoods.jl
new file mode 100644
index 00000000..109083aa
--- /dev/null
+++ b/src/Likelihoods.jl
@@ -0,0 +1,152 @@
+
+
+### Likelihoods
+
+"""
+    likelihood(DM::DataModel, θ::AbstractVector) -> Real
+Calculates the likelihood ``L(\\mathrm{data} \\, | \\, \\theta)`` a `DataModel` and a parameter configuration ``\\theta``.
+"""
+likelihood(args...; kwargs...) = exp(loglikelihood(args...; kwargs...))
+
+
+
+## Prefix underscore for likelihood, Score and FisherMetric indicates that Prior has already been accounted for upstream
+loglikelihood(DM::AbstractDataModel; kwargs...) = LogLikelihood(θ::AbstractVector{<:Number}; Kwargs...) = loglikelihood(DM, θ; kwargs..., Kwargs...)
+Negloglikelihood(DM::AbstractDataModel; kwargs...) = NegativeLogLikelihood(θ::AbstractVector{<:Number}; Kwargs...) = -loglikelihood(DM, θ; kwargs..., Kwargs...)
+
+# import Distributions.loglikelihood
+"""
+    loglikelihood(DM::DataModel, θ::AbstractVector) -> Real
+Calculates the logarithm of the likelihood ``L``, i.e. ``\\ell(\\mathrm{data} \\, | \\, \\theta) \\coloneqq \\mathrm{ln} \\big( L(\\mathrm{data} \\, | \\, \\theta) \\big)`` given a `DataModel` and a parameter configuration ``\\theta``.
+"""
+loglikelihood(DM::AbstractDataModel, θ::AbstractVector{<:Number}, LogPriorFn::Union{Nothing,Function}=LogPrior(DM); kwargs...) = loglikelihood(Data(DM), Predictor(DM), θ, LogPriorFn; kwargs...)
+
+loglikelihood(DS::AbstractDataSet, model::ModelOrFunction, θ::AbstractVector{<:Number}, LogPriorFn::Nothing; kwargs...) = _loglikelihood(DS, model, θ; kwargs...)
+loglikelihood(DS::AbstractDataSet, model::ModelOrFunction, θ::AbstractVector{<:Number}, LogPriorFn::Function; kwargs...) = _loglikelihood(DS, model, θ; kwargs...) + EvalLogPrior(LogPriorFn, θ)
+
+
+# Specialize this for different DataSet types
+@inline function _loglikelihood(DS::AbstractDataSet, model::ModelOrFunction, θ::AbstractVector{<:Number}; kwargs...)
+    -0.5*(DataspaceDim(DS)*log(2π) - logdetInvCov(DS) + InnerProduct(yInvCov(DS), ydata(DS)-EmbeddingMap(DS, model, θ; kwargs...)))
+end
+
+
+
+InnerProduct(Mat::AbstractMatrix, Y::AbstractVector) = transpose(Y) * Mat * Y
+# InnerProduct(Mat::PDMats.PDMat, Y::AbstractVector) = (R = Mat.chol.U * Y;  dot(R,R))
+
+InnerProductV(Mat::AbstractMatrix, Y::AbstractVector) = @tullio Res := Y[i] * Mat[i,j] * Y[j]
+InnerProductV(Mat::Diagonal, Y::AbstractVector) = @tullio Res := Mat.diag[j] * Y[j]^2
+
+# Does not hit BLAS, sadly
+"""
+    InnerProductChol(Mat::AbstractMatrix, Y::AbstractVector{T}) -> T
+Computes ``|| Mat * Y ||^2``, i.e. ``Y^t \\, * (Mat^t * Mat) * Y``.
+"""
+function InnerProductChol(Mat::UpperTriangular, Y::AbstractVector{T})::T where T <: Number
+    @assert size(Mat,2) == length(Y)
+    Res = zero(T);    temp = zero(T);    n = size(Mat,2)
+    @inbounds for i in 1:n
+        temp = dot(view(Mat,i,i:n), view(Y,i:n))
+        Res += temp^2
+    end;    Res
+end
+function InnerProductChol(Mat::Diagonal, Y::AbstractVector{T})::T where T <: Number
+    @assert length(Mat.diag) == length(Y)
+    # sum(abs2, Mat.diag .* Y)
+    sum((Mat.diag .* Y).^2) # faster for differentiation
+end
+function InnerProductChol(Mat::AbstractMatrix, Y::AbstractVector{T})::T where T <: Number
+    @assert size(Mat,1) == size(Mat,2) == length(Y)
+    sum(abs2, Mat*Y)
+end
+
+
+
+
+AutoScore(DM::AbstractDataModel, θ::AbstractVector{<:Number}, LogPriorFn::Union{Nothing, Function}=LogPrior(DM); kwargs...) = AutoScore(Data(DM), Predictor(DM), θ, LogPriorFn; kwargs...)
+AutoScore(DS::AbstractDataSet, model::ModelOrFunction, θ::AbstractVector{<:Number}, LogPriorFn::Nothing; kwargs...) = _AutoScore(DS, model, θ; kwargs...)
+function AutoScore(DS::AbstractDataSet, model::ModelOrFunction, θ::AbstractVector{<:Number}, LogPriorFn::Function; ADmode::Union{Symbol,Val}=Val(:ForwardDiff), kwargs...)
+    _AutoScore(DS, model, θ; ADmode=ADmode, kwargs...) + EvalLogPriorGrad(LogPriorFn, θ; ADmode=ADmode)
+end
+_AutoScore(DS::AbstractDataSet, model::ModelOrFunction, θ::AbstractVector{<:Number}; ADmode::Union{Symbol,Val}=Val(:ForwardDiff), kwargs...) = GetGrad(ADmode, x->_loglikelihood(DS, model, x; kwargs...))(θ)
+
+
+AutoMetric(DM::AbstractDataModel, θ::AbstractVector{<:Number}, LogPriorFn::Union{Nothing, Function}=LogPrior(DM); kwargs...) = AutoMetric(Data(DM), Predictor(DM), θ, LogPriorFn; kwargs...)
+AutoMetric(DS::AbstractDataSet, model::ModelOrFunction, θ::AbstractVector{<:Number}, LogPriorFn::Nothing; kwargs...) = _AutoMetric(DS, model, θ; kwargs...)
+function AutoMetric(DS::AbstractDataSet, model::ModelOrFunction, θ::AbstractVector{<:Number}, LogPriorFn::Function; ADmode::Union{Symbol,Val}=Val(:ForwardDiff), kwargs...)
+    _AutoMetric(DS, model, θ; ADmode=ADmode, kwargs...) - EvalLogPriorHess(LogPriorFn, θ; ADmode=ADmode)
+end
+_AutoMetric(DS::AbstractDataSet, model::ModelOrFunction, θ::AbstractVector{<:Number}; ADmode::Union{Symbol,Val}=Val(:ForwardDiff), kwargs...) = GetHess(ADmode, x->-_loglikelihood(DS, model, x; kwargs...))(θ)
+
+
+
+Score(DM::AbstractDataModel; kwargs...) = LogLikelihoodGradient(θ::AbstractVector{<:Number}; Kwargs...) = Score(DM, θ; kwargs..., Kwargs...)
+NegScore(DM::AbstractDataModel; kwargs...) = NegativeLogLikelihoodGradient(θ::AbstractVector{<:Number}; Kwargs...) = -Score(DM, θ; kwargs..., Kwargs...)
+
+"""
+    Score(DM::DataModel, θ::AbstractVector{<:Number}; ADmode::Val=Val(:ForwardDiff))
+Calculates the gradient of the log-likelihood ``\\ell`` with respect to a set of parameters ``\\theta``.
+`ADmode=Val(false)` computes the Score by separately evaluating the `model` as well as the Jacobian `dmodel` provided in `DM`.
+Other choices of `ADmode` directly compute the Score by differentiating the formula the log-likelihood, i.e. only one evaluation on a dual variable is performed.
+"""
+Score(DM::AbstractDataModel, θ::AbstractVector{<:Number}, LogPriorFn::Union{Nothing,Function}=LogPrior(DM); kwargs...) = Score(Data(DM), Predictor(DM), dPredictor(DM), θ, LogPriorFn; kwargs...)
+
+Score(DS::AbstractDataSet, model::ModelOrFunction, dmodel::ModelOrFunction, θ::AbstractVector{<:Number}, LogPriorFn::Nothing; kwargs...) = _Score(DS, model, dmodel, θ; kwargs...)
+Score(DS::AbstractDataSet, model::ModelOrFunction, dmodel::ModelOrFunction, θ::AbstractVector{<:Number}, LogPriorFn::Function; kwargs...) = _Score(DS, model, dmodel, θ; kwargs...) + EvalLogPriorGrad(LogPriorFn, θ)
+
+
+
+_Score(DS::AbstractDataSet, model::ModelOrFunction, dmodel::ModelOrFunction, θ::AbstractVector{<:Number}; ADmode::Val=Val(:ForwardDiff), kwargs...) = _Score(DS, model, dmodel, θ, ADmode; kwargs...)
+# Delegate to AutoScore
+_Score(DS::AbstractDataSet, model::ModelOrFunction, dmodel::ModelOrFunction, θ::AbstractVector{<:Number}, ADmode::Val{true}; kwargs...) = _AutoScore(DS, model, θ; ADmode=ADmode, kwargs...)
+_Score(DS::AbstractDataSet, model::ModelOrFunction, dmodel::ModelOrFunction, θ::AbstractVector{<:Number}, ADmode::Union{Symbol,Val}; kwargs...) = _AutoScore(DS, model, θ; ADmode=ADmode, kwargs...)
+
+# Specialize this for different DataSet types
+@inline function _Score(DS::AbstractDataSet, model::ModelOrFunction, dmodel::ModelOrFunction, θ::AbstractVector{<:Number}, ADmode::Val{false}; kwargs...)
+    transpose(EmbeddingMatrix(DS,dmodel,θ; kwargs...)) * (yInvCov(DS) * (ydata(DS) - EmbeddingMap(DS,model,θ; kwargs...)))
+end
+
+
+
+
+
+"""
+    LiftedLogLikelihood(DM::AbstractDataModel) -> Function
+Computes log-likelihood on the extended data space ``\\hat{\\ell} : \\mathcal{X}^N \\times\\mathcal{Y}^N \\longrightarrow \\mathbb{R}``.
+Should be maximized.
+!!! note
+    CANNOT ACCOUNT FOR PRIORS.
+"""
+LiftedLogLikelihood(DM::Union{AbstractDataModel,AbstractDataSet}) = (G = dist(DM);  ℓ(Z::AbstractVector{<:Number}) = logpdf(G, Z))
+"""
+    LiftedCost(DM::AbstractDataModel) -> Function
+Computes negative log-likelihood as cost function on the extended data space ``C : \\mathcal{X}^N \\times\\mathcal{Y}^N \\longrightarrow \\mathbb{R}``.
+Should be minimized.
+!!! note
+    CANNOT ACCOUNT FOR PRIORS.
+"""
+LiftedCost(DM::Union{AbstractDataModel,AbstractDataSet}) = (G = dist(DM);  Negativeℓ(Z::AbstractVector{<:Number}) = -logpdf(G, Z))
+
+"""
+    LiftedEmbedding(DM::AbstractDataModel) -> Function
+Constructs lifted embedding map from initial space into extended dataspace ``\\hat{h} : \\mathcal{X}^N \\times \\mathcal{M} \\longrightarrow \\mathcal{X}^N \\times\\mathcal{Y}^N`` effecting
+``\\xi = (x_\\text{opt}, \\theta) \\longmapsto \\hat{h}(\\xi) = (x_\\text{opt}, h(x_\\text{opt}, \\theta))``.
+"""
+LiftedEmbedding(DM::AbstractDataModel) = LiftedEmbedding(Data(DM), Predictor(DM), pdim(DM))
+function LiftedEmbedding(DS::AbstractDataSet, Model::ModelOrFunction, pd::Int)
+    ĥ(ξ::AbstractVector; kwargs...) = ĥ(view(ξ,1:length(ξ)-pd), view(ξ,length(ξ)-pd+1:length(ξ)); kwargs...)
+    ĥ(xdat::AbstractVector, θ::AbstractVector{<:Number}; kwargs...) = [xdat; EmbeddingMap(DS, Model, θ, Windup(xdat, xdim(DS)); kwargs...)]
+end
+
+"""
+    FullLiftedLogLikelihood(DM::AbstractDataModel) -> Function
+Computes the full likelihood ``\\hat{\\ell} : \\mathcal{X}^N \\times\\mathcal{M}^N \\longrightarrow \\mathbb{R}`` given `Xθ` from initial space INCLUDING PRIOR.
+"""
+FullLiftedLogLikelihood(DM::AbstractDataModel) = FullLiftedLogLikelihood(Data(DM), Predictor(DM), LogPrior(DM), pdim(DM))
+FullLiftedLogLikelihood(DS::AbstractDataSet, model::ModelOrFunction, LogPriorFn::Nothing, pd::Int) = LiftedLogLikelihood(DS)∘LiftedEmbedding(DS, model, pd)
+function FullLiftedLogLikelihood(DS::AbstractDataSet, model::ModelOrFunction, LogPriorFn::Function, pd::Int)
+    L = LiftedLogLikelihood(DS)∘LiftedEmbedding(DS, model, pd)
+    ℓ(Xθ::AbstractVector{<:Number}; kwargs...) = L(Xθ; kwargs...) + EvalLogPrior(LogPriorFn, view(Xθ, length(Xθ)-pd+1:length(Xθ)))
+end
+FullLiftedNegLogLikelihood(args...; kwargs...) = (L=FullLiftedLogLikelihood(args...; kwargs...); Xθ::AbstractVector{<:Number}->-L(Xθ))
diff --git a/src/ModelPredictions.jl b/src/ModelPredictions.jl
new file mode 100644
index 00000000..9c3f42c3
--- /dev/null
+++ b/src/ModelPredictions.jl
@@ -0,0 +1,104 @@
+
+
+
+### for-loop typically slower than reduce(vcat, ...)
+### Apparently curve_fit() throws an error in conjuction with ForwardDiff when reinterpret() is used
+# Reduction(X::AbstractVector{<:SVector{Len,T}}) where Len where T = reinterpret(T, X)
+Reduction(X::AbstractVector{<:AbstractVector}) = reduce(vcat, X)
+Reduction(X::AbstractVector{<:Number}) = X
+Reduction(X::AbstractVector{<:SubArray{<:Number, 0}}) = [@inbounds X[i][1] for i in 1:length(X)]
+
+
+# h(θ) ∈ Dataspace
+"""
+    EmbeddingMap(DM::AbstractDataModel, θ::AbstractVector{<:Number}) -> Vector
+Returns a vector of the collective predictions of the `model` as evaluated at the x-values and the parameter configuration ``\\theta``.
+```
+h(\\theta) \\coloneqq \\big(y_\\mathrm{model}(x_1;\\theta),...,y_\\mathrm{model}(x_N;\\theta)\\big) \\in \\mathcal{D}
+```
+"""
+EmbeddingMap(DM::AbstractDataModel, θ::AbstractVector{<:Number}, woundX::AbstractVector=WoundX(DM); kwargs...) = EmbeddingMap(Data(DM), Predictor(DM), θ, woundX; kwargs...)
+EmbeddingMap(DS::AbstractDataSet, model::ModelOrFunction, θ::AbstractVector{<:Number}, woundX::AbstractVector=WoundX(DS); kwargs...) = _CustomOrNot(DS, model, θ, woundX; kwargs...)
+EmbeddingMap(DS::Val, model::ModelOrFunction, θ::AbstractVector{<:Number}, woundX::AbstractVector; kwargs...) = _CustomOrNot(DS, model, θ, woundX; kwargs...)
+
+_CustomOrNot(DS::Union{Val,AbstractDataSet}, model::Function, θ::AbstractVector{<:Number}, woundX::AbstractVector; kwargs...) = _CustomOrNot(DS, model, θ, woundX, Val(false), Val(false); kwargs...)
+_CustomOrNot(DS::Union{Val,AbstractDataSet}, M::ModelMap, θ::AbstractVector{<:Number}, woundX::AbstractVector; kwargs...) = _CustomOrNot(DS, M.Map, θ, woundX, M.CustomEmbedding, M.inplace; kwargs...)
+
+# Specialize this for different Dataset types
+_CustomOrNot(::Union{Val,AbstractDataSet}, model::ModelOrFunction, θ::AbstractVector{<:Number}, woundX::AbstractVector, custom::Val{false}, inplace::Val{false}; kwargs...) = Reduction(map(x->model(x,θ; kwargs...), woundX))
+_CustomOrNot(::Union{Val,AbstractDataSet}, model::ModelOrFunction, θ::AbstractVector{<:Number}, woundX::AbstractVector, custom::Val{true}, inplace::Val{false}; kwargs...) = model(woundX, θ; kwargs...)
+
+
+function _CustomOrNot(DS::AbstractDataSet, model!::ModelOrFunction, θ::AbstractVector{<:Number}, woundX::AbstractVector, custom::Val{true}, inplace::Val{true}; kwargs...)
+    Y = Vector{suff(θ)}(undef, length(woundX)*ydim(DS))
+    model!(Y, woundX, θ; kwargs...);    Y
+end
+function _CustomOrNot(DS::AbstractDataSet, model!::ModelOrFunction, θ::AbstractVector{<:Number}, woundX::AbstractVector, custom::Val{false}, inplace::Val{true}; kwargs...)
+    Y = Vector{suff(θ)}(undef, length(woundX)*ydim(DS))
+    EmbeddingMap!(Y, DS, model!, θ, woundX; kwargs...);     Y
+end
+
+
+function EmbeddingMap!(Y::AbstractVector{<:Number}, DS::AbstractDataSet, model!::ModelOrFunction, θ::AbstractVector{<:Number}, woundX::AbstractVector=WoundX(DS); kwargs...)
+    EmbeddingMap!(Y, model!, θ, woundX, Val(ydim(DS)); kwargs...)
+end
+function EmbeddingMap!(Y::AbstractVector{<:Number}, model!::ModelOrFunction, θ::AbstractVector{<:Number}, woundX::AbstractVector, Ydim::Val{1}; kwargs...)
+    @inbounds for i in Base.OneTo(length(Y))
+        model!(Y[i], woundX[i], θ; kwargs...)
+    end
+end
+function EmbeddingMap!(Y::AbstractVector{<:Number}, model!::ModelOrFunction, θ::AbstractVector{<:Number}, woundX::AbstractVector, Ydim::Val{T}; kwargs...) where T
+    @inbounds for (i, row) in enumerate(Iterators.partition(1:length(Y), T))
+        model!(view(Y,row), woundX[i], θ; kwargs...)
+    end
+end
+# Fallback for out-of-place models
+function EmbeddingMap!(Y::AbstractVector{<:Number}, DS::AbstractDataSet, model!::ModelMap{false}, θ::AbstractVector{<:Number}, woundX::AbstractVector=WoundX(DS); kwargs...)
+    copyto!(Y, EmbeddingMap(DS, model!, θ, woundX; kwargs...))
+end
+
+
+
+"""
+    EmbeddingMatrix(DM::AbstractDataModel, θ::AbstractVector{<:Number}) -> Matrix
+Returns the jacobian of the embedding map as evaluated at the x-values and the parameter configuration ``\\theta``.
+"""
+EmbeddingMatrix(DM::AbstractDataModel, θ::AbstractVector{<:Number}, woundX::AbstractVector=WoundX(DM); kwargs...) = EmbeddingMatrix(Data(DM), dPredictor(DM), θ, woundX; kwargs...)
+EmbeddingMatrix(DS::AbstractDataSet, dmodel::ModelOrFunction, θ::AbstractVector{<:Number}, woundX::AbstractVector=WoundX(DS); kwargs...) = _CustomOrNotdM(DS, dmodel, θ, woundX; kwargs...)
+EmbeddingMatrix(DS::Val, dmodel::ModelOrFunction, θ::AbstractVector{<:Number}, woundX::AbstractVector; kwargs...) = _CustomOrNotdM(DS, dmodel, θ, woundX; kwargs...)
+
+_CustomOrNotdM(DS::Union{Val,AbstractDataSet}, dmodel::Function, θ::AbstractVector{<:Number}, woundX::AbstractVector; kwargs...) = _CustomOrNotdM(DS, dmodel, floatify(θ), woundX, Val(false), Val(false); kwargs...)
+_CustomOrNotdM(DS::Union{Val,AbstractDataSet}, dM::ModelMap, θ::AbstractVector{<:Number}, woundX::AbstractVector; kwargs...) = _CustomOrNotdM(DS, dM.Map, floatify(θ), woundX, dM.CustomEmbedding, dM.inplace; kwargs...)
+
+# Specialize this for different Dataset types
+_CustomOrNotdM(::Union{Val,AbstractDataSet}, dmodel::ModelOrFunction, θ::AbstractVector{<:Number}, woundX::AbstractVector, custom::Val{false}, inplace::Val{false}; kwargs...) = reduce(vcat, map(x->dmodel(x,θ; kwargs...), woundX))
+_CustomOrNotdM(::Union{Val,AbstractDataSet}, dmodel::ModelOrFunction, θ::AbstractVector{<:Number}, woundX::AbstractVector, custom::Val{true}, inplace::Val{false}; kwargs...) = dmodel(woundX, θ; kwargs...)
+
+
+function _CustomOrNotdM(DS::AbstractDataSet, dmodel!::ModelOrFunction, θ::AbstractVector{<:Number}, woundX::AbstractVector, custom::Val{true}, inplace::Val{true}; kwargs...)
+    J = Matrix{suff(θ)}(undef, length(woundX)*ydim(DS), length(θ))
+    dmodel!(J, woundX, θ; kwargs...);   J
+end
+function _CustomOrNotdM(DS::AbstractDataSet, dmodel!::ModelOrFunction, θ::AbstractVector{<:Number}, woundX::AbstractVector, custom::Val{false}, inplace::Val{true}; kwargs...)
+    J = Matrix{suff(θ)}(undef, length(woundX)*ydim(DS), length(θ))
+    EmbeddingMatrix!(J, DS, dmodel!, θ, woundX; kwargs...);     J
+end
+
+
+function EmbeddingMatrix!(J::AbstractMatrix{<:Number}, DS::AbstractDataSet, dmodel!::ModelOrFunction, θ::AbstractVector{<:Number}, woundX::AbstractVector=WoundX(DS); kwargs...)
+    EmbeddingMatrix!(J, dmodel!, θ, woundX, Val(ydim(DS)); kwargs...)
+end
+function EmbeddingMatrix!(J::AbstractMatrix{<:Number}, dmodel!::ModelOrFunction, θ::AbstractVector{<:Number}, woundX::AbstractVector, Ydim::Val{1}; kwargs...)
+    @inbounds for row in Base.OneTo(size(J,1))
+        dmodel!(view(J,row,:), woundX[row], θ; kwargs...)
+    end
+end
+function EmbeddingMatrix!(J::AbstractMatrix{<:Number}, dmodel!::ModelOrFunction, θ::AbstractVector{<:Number}, woundX::AbstractVector, Ydim::Val{T}; kwargs...) where T
+    @inbounds for (i, row) in enumerate(Iterators.partition(1:size(J,1), T))
+        dmodel!(view(J,row,:), woundX[i], θ; kwargs...)
+    end
+end
+# Fallback for out-of-place models
+function EmbeddingMatrix!(J::AbstractMatrix{<:Number}, DS::AbstractDataSet, dmodel!::ModelMap{false}, θ::AbstractVector{<:Number}, woundX::AbstractVector=WoundX(DS); kwargs...)
+    copyto!(J, EmbeddingMatrix(DS, dmodel!, θ, woundX; kwargs...))
+end
diff --git a/src/NumericalTools/Miscellaneous.jl b/src/NumericalTools/Miscellaneous.jl
index 00bf5e42..17cd8925 100644
--- a/src/NumericalTools/Miscellaneous.jl
+++ b/src/NumericalTools/Miscellaneous.jl
@@ -142,36 +142,6 @@ InvChisqCDF(k::Int, p::Real; kwargs...) = 2gamma_inc_inv(k/2., p, 1-p)
 InvChisqCDF(k::Int, p::BigFloat; tol::Real=GetH(p)) = invert(x->ChisqCDF(k, x), p; tol=tol)
 
 
-InnerProduct(Mat::AbstractMatrix, Y::AbstractVector) = transpose(Y) * Mat * Y
-# InnerProduct(Mat::PDMats.PDMat, Y::AbstractVector) = (R = Mat.chol.U * Y;  dot(R,R))
-
-InnerProductV(Mat::AbstractMatrix, Y::AbstractVector) = @tullio Res := Y[i] * Mat[i,j] * Y[j]
-InnerProductV(Mat::Diagonal, Y::AbstractVector) = @tullio Res := Mat.diag[j] * Y[j]^2
-
-
-
-# Does not hit BLAS, sadly
-"""
-    InnerProductChol(Mat::AbstractMatrix, Y::AbstractVector{T}) -> T
-Computes ``|| Mat * Y ||^2``, i.e. ``Y^t \\, * (Mat^t * Mat) * Y``.
-"""
-function InnerProductChol(Mat::UpperTriangular, Y::AbstractVector{T})::T where T <: Number
-    @assert size(Mat,2) == length(Y)
-    Res = zero(T);    temp = zero(T);    n = size(Mat,2)
-    @inbounds for i in 1:n
-        temp = dot(view(Mat,i,i:n), view(Y,i:n))
-        Res += temp^2
-    end;    Res
-end
-function InnerProductChol(Mat::Diagonal, Y::AbstractVector{T})::T where T <: Number
-    @assert length(Mat.diag) == length(Y)
-    # sum(abs2, Mat.diag .* Y)
-    sum((Mat.diag .* Y).^2) # faster for differentiation
-end
-function InnerProductChol(Mat::AbstractMatrix, Y::AbstractVector{T})::T where T <: Number
-    @assert size(Mat,1) == size(Mat,2) == length(Y)
-    sum(abs2, Mat*Y)
-end
 
 
 import Base.==