struct DeepGPModel{T,L}
    layers::T
    nLayers::Int
    likelihood::L
end

function DeepGPModel(layers...; likelihood=Normal())
    # Here do some checks with likelihood and last layer
    # Check dimensions between layers
    return DeepGPModel{typeof(layers),typeof(likelihood)}(
        layers, length(layers), likelihood
    )
end

# Base.:∘(l1::AbstractGPLayer,l2::AbstractGPLayer) = DeepGPModel([l1,l2])
Base.getindex(m::DeepGPModel, i::Int) = m.layers[i]
Base.length(m::DeepGPModel) = m.length
Base.iterate(m::DeepGPModel) = iterate(m.layers)

Base.iterate(m::DeepGPModel, state) = iterate(m.layers, state)

params(m::DeepGPModel) = params.(m.layers) #Eventually add parameters from likelihood

function ELBO(m::DeepGPModel, X, y; nsamples=1)
    return expec_log_like(m, X, y; nsamples=nsamples) - KL(m)
end

function KL(m::DeepGPModel)
    return sum(KL, m)
end

propagate(m::DeepGPModel, f) = propagate(Base.tail(m.layers), propagate(first(m), f))
propagate(::Tuple{}, f) = first.(f) # Final output
propagate(layers::Tuple, f) = propagate(Base.tail(layers), propagate(first(layers), f))

"""
    predict(m, X; nsamples)

Create `nsamples` samples with model `m` given samples `X`
"""
function predict(m::DeepGPModel, X::AbstractVector; nsamples=1)
    return reshape(propagate(m, repeat(X, nsamples)), :, nsamples)
end

function expec_log_like(m::DeepGPModel, X::AbstractVector, y; nsamples=1)
    f = predict(m, X; nsamples=nsamples)
    loss = mean(_logpdf.(Ref(m), eachcol(f), Ref(y)))
    return loss
end

function _logpdf(m::DeepGPModel, f, y)
    return sum(_logpdf.(Ref(m.likelihood), f, y))
end

function _logpdf(l::Normal, f::Real, y::Real)
    return -0.5 * abs2(y - f) / l.σ^2
    # Distributions.logpdf(l,f-y)
end
