diff --git a/Project.toml b/Project.toml
index 7b32ca1..259399e 100644
--- a/Project.toml
+++ b/Project.toml
@@ -15,6 +15,7 @@ Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 CategoricalArrays = "0.10"
 MLJBase = "0.20, 0.21"
 MLJModelInterface = "1"
+NaturalSort = "1"
 Plots = "1"
 julia = "1.7"
 
diff --git a/dev/quick_tour/notebook.jl b/dev/quick_tour/notebook.jl
index a392be5..8b3dd33 100644
--- a/dev/quick_tour/notebook.jl
+++ b/dev/quick_tour/notebook.jl
@@ -7,7 +7,14 @@ using InteractiveUtils
 # This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error).
 macro bind(def, element)
     quote
-        local iv = try Base.loaded_modules[Base.PkgId(Base.UUID("6e696c72-6542-2067-7265-42206c756150"), "AbstractPlutoDingetjes")].Bonds.initial_value catch; b -> missing; end
+        local iv = try
+            Base.loaded_modules[Base.PkgId(
+                Base.UUID("6e696c72-6542-2067-7265-42206c756150"),
+                "AbstractPlutoDingetjes",
+            )].Bonds.initial_value
+        catch
+            b -> missing
+        end
         local el = $(esc(element))
         global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : iv(el)
         el
@@ -16,17 +23,17 @@ end
 
 # ╔═╡ aad62ef1-4136-4732-a9e6-3746524978ee
 begin
-	using ConformalPrediction
-	using Distributions
-	using EvoTrees
-	using LightGBM
-	using MLJ
-	using MLJDecisionTreeInterface
-	using MLJLinearModels
-	using NearestNeighborModels
-	using Plots
-	using PlutoUI
-	include("utils.jl")
+    using ConformalPrediction
+    using Distributions
+    using EvoTrees
+    using LightGBM
+    using MLJ
+    using MLJDecisionTreeInterface
+    using MLJLinearModels
+    using NearestNeighborModels
+    using Plots
+    using PlutoUI
+    include("utils.jl")
 end;
 
 # ╔═╡ bc0d7575-dabd-472d-a0ce-db69d242ced8
@@ -49,18 +56,18 @@ First, we create a simple helper function that generates our data:
 
 # ╔═╡ 2f1c8da3-77dc-4bd7-8fa4-7669c2861aaa
 begin
-	function get_data(N=600, xmax=3.0, noise=0.5; fun::Function=fun(X) = X * sin(X))
-		# Inputs:
-		d = Distributions.Uniform(-xmax, xmax)
-		X = rand(d, N)
-		X = MLJ.table(reshape(X, :, 1))
-		
-		# Outputs:
-		ε = randn(N) .* noise
-		y = @.(fun(X.x1)) + ε
-		y = vec(y)
-		return X, y
-	end
+    function get_data(N = 600, xmax = 3.0, noise = 0.5; fun::Function = fun(X) = X * sin(X))
+        # Inputs:
+        d = Distributions.Uniform(-xmax, xmax)
+        X = rand(d, N)
+        X = MLJ.table(reshape(X, :, 1))
+
+        # Outputs:
+        ε = randn(N) .* noise
+        y = @.(fun(X.x1)) + ε
+        y = vec(y)
+        return X, y
+    end
 end;
 
 # ╔═╡ eb251479-ce0f-4158-8627-099da3516c73
@@ -78,20 +85,27 @@ The slides can be used to change the number of observations `N`, the maximum (an
 
 # ╔═╡ 931ce259-d5fb-4a56-beb8-61a69a2fc09e
 begin
-	data_dict = Dict(
-		"N" => (500:100:5000,1000),
-		"noise" => (0.1:0.1:1.0,0.5),
-		"xmax" => (1:10,5),
-	)
-	@bind data_specs multi_slider(data_dict, title="Parameters")
+    data_dict = Dict(
+        "N" => (500:100:5000, 1000),
+        "noise" => (0.1:0.1:1.0, 0.5),
+        "xmax" => (1:10, 5),
+    )
+    @bind data_specs multi_slider(data_dict, title = "Parameters")
 end
 
 # ╔═╡ f0106aa5-b1c5-4857-af94-2711f80d25a8
 begin
-	X, y = get_data(data_specs.N, data_specs.xmax, data_specs.noise; fun=f)
-	scatter(X.x1, y, label="Observed data")
-	xrange = range(-data_specs.xmax,data_specs.xmax,length=50)
-	plot!(xrange, @.(f(xrange)), lw=4, label="Ground truth", ls=:dash, colour=:black)
+    X, y = get_data(data_specs.N, data_specs.xmax, data_specs.noise; fun = f)
+    scatter(X.x1, y, label = "Observed data")
+    xrange = range(-data_specs.xmax, data_specs.xmax, length = 50)
+    plot!(
+        xrange,
+        @.(f(xrange)),
+        lw = 4,
+        label = "Ground truth",
+        ls = :dash,
+        colour = :black,
+    )
 end
 
 # ╔═╡ 2fe1065e-d1b8-4e3c-930c-654f50349222
@@ -111,7 +125,7 @@ To start with, let's split our data into a training and test set:
 """
 
 # ╔═╡ 3a4fe2bc-387c-4d7e-b45f-292075a01bcd
-train, test = partition(eachindex(y), 0.4, 0.4, shuffle=true);
+train, test = partition(eachindex(y), 0.4, 0.4, shuffle = true);
 
 # ╔═╡ a34b8c07-08e0-4a0e-a0f9-8054b41b038b
 md"Now let's choose a model for our regression task:"
@@ -121,8 +135,8 @@ md"Now let's choose a model for our regression task:"
 
 # ╔═╡ 292978a2-1941-44d3-af5b-13456d16b656
 begin
-	Model = eval(tested_atomic_models[:regression][model_name])
-	model = Model()
+    Model = eval(tested_atomic_models[:regression][model_name])
+    model = Model()
 end;
 
 # ╔═╡ 10340f3f-7981-42da-846a-7599a9edb7f3
@@ -137,7 +151,7 @@ mach_raw = machine(model, X, y);
 md"Then we fit the machine to the training data:"
 
 # ╔═╡ aabfbbfb-7fb0-4f37-9a05-b96207636232
-MLJ.fit!(mach_raw, rows=train, verbosity=0);
+MLJ.fit!(mach_raw, rows = train, verbosity = 0);
 
 # ╔═╡ 5506e1b5-5f2f-4972-a845-9c0434d4b31c
 md"""
@@ -146,13 +160,20 @@ The chart below shows the resulting point predictions for the test data set:
 
 # ╔═╡ 9bb977fe-d7e0-4420-b472-a50e8bd6d94f
 begin
-	Xtest = MLJ.matrix(selectrows(X, test))
-	ytest = y[test]
-	ŷ = MLJ.predict(mach_raw, Xtest)
-	scatter(vec(Xtest), vec(ytest), label="Observed")
-	_order = sortperm(vec(Xtest))
-	plot!(vec(Xtest)[_order], vec(ŷ)[_order], lw=4, label="Predicted")
-	plot!(xrange, @.(f(xrange)), lw=2, ls=:dash, colour=:black, label="Ground truth")
+    Xtest = MLJ.matrix(selectrows(X, test))
+    ytest = y[test]
+    ŷ = MLJ.predict(mach_raw, Xtest)
+    scatter(vec(Xtest), vec(ytest), label = "Observed")
+    _order = sortperm(vec(Xtest))
+    plot!(vec(Xtest)[_order], vec(ŷ)[_order], lw = 4, label = "Predicted")
+    plot!(
+        xrange,
+        @.(f(xrange)),
+        lw = 2,
+        ls = :dash,
+        colour = :black,
+        label = "Ground truth",
+    )
 end
 
 # ╔═╡ 36eef47f-ad55-49be-ac60-7aa1cf50e61a
@@ -186,7 +207,7 @@ Then we fit the machine to the data:
 """
 
 # ╔═╡ 6b574688-ff3c-441a-a616-169685731883
-MLJ.fit!(mach, rows=train, verbosity=0);
+MLJ.fit!(mach, rows = train, verbosity = 0);
 
 # ╔═╡ da6e8f90-a3f9-4d06-86ab-b0f6705bbf54
 md"""
@@ -196,15 +217,22 @@ Now let us look at the predictions for our test data again. The chart below show
 """
 
 # ╔═╡ 797746e9-235f-4fb1-8cdb-9be295b54bbe
-@bind coverage Slider(0.1:0.1:1.0, default=0.8, show_value=true)
+@bind coverage Slider(0.1:0.1:1.0, default = 0.8, show_value = true)
 
 # ╔═╡ ad3e290b-c1f5-4008-81c7-a1a56ab10563
 begin
-	_conf_model = conformal_model(model, coverage=coverage)
-	_mach = machine(_conf_model, X, y)
-	MLJ.fit!(_mach, rows=train, verbosity=0)
-	plot(_mach.model, _mach.fitresult, Xtest, ytest, zoom=0, observed_lab="Test points")
-	plot!(xrange, @.(f(xrange)), lw=2, ls=:dash, colour=:black, label="Ground truth")
+    _conf_model = conformal_model(model, coverage = coverage)
+    _mach = machine(_conf_model, X, y)
+    MLJ.fit!(_mach, rows = train, verbosity = 0)
+    plot(_mach.model, _mach.fitresult, Xtest, ytest, zoom = 0, observed_lab = "Test points")
+    plot!(
+        xrange,
+        @.(f(xrange)),
+        lw = 2,
+        ls = :dash,
+        colour = :black,
+        label = "Ground truth",
+    )
 end
 
 # ╔═╡ b3a88859-0442-41ff-bfea-313437042830
@@ -225,29 +253,32 @@ To verify the marginal coverage property empirically we can look at the empirica
 
 # ╔═╡ d1140af9-608a-4669-9595-aee72ffbaa46
 begin
-	model_evaluation = evaluate!(_mach,operation=MLJ.predict,measure=emp_coverage, verbosity=0);
-	println("Empirical coverage: $(round(model_evaluation.measurement[1], digits=3))")
-	println("Coverage per fold: $(round.(model_evaluation.per_fold[1], digits=3))")
+    model_evaluation =
+        evaluate!(_mach, operation = MLJ.predict, measure = emp_coverage, verbosity = 0)
+    println("Empirical coverage: $(round(model_evaluation.measurement[1], digits=3))")
+    println("Coverage per fold: $(round.(model_evaluation.per_fold[1], digits=3))")
 end
 
 # ╔═╡ f742440b-258e-488a-9c8b-c9267cf1fb99
 begin
-	ncal = Int(conf_model.train_ratio * data_specs.N)
-	Markdown.parse("""
-	The empirical coverage rate should be close to the desired level of coverage. In most cases it will be slightly higher, since ``(1-\\alpha)`` is a lower bound.
-	
-	> Found an empirical coverage rate that is slightly lower than desired? The coverage property is "marginal" in the sense that the probability averaged over the randomness in the data. For most purposes a large enough calibration set size (``n>1000``) mitigates that randomness enough. Depending on your choices above, the calibration set may be quite small (currently $ncal), which can lead to **coverage slack** (see Section 3 in the [tutorial](https://arxiv.org/pdf/2107.07511.pdf)).
-	
-	### *So what's happening under the hood?*
-	
-	Inductive Conformal Prediction (also referred to as Split Conformal Prediction) broadly speaking works as follows:
-	
-	1. Partition the training into a proper training set and a separate calibration set
-	2. Train the machine learning model on the proper training set.
-	3. Using some heuristic notion of uncertainty (e.g. absolute error in the regression case) compute nonconformity scores using the calibration data and the fitted model. 
-	4. For the given coverage ratio compute the corresponding quantile of the empirical distribution of nonconformity scores.
-	5. For the given quantile and test sample ``X_{\\text{test}}``, form the corresponding conformal prediction set like so: ``C(X_{\\text{test}})=\\{y:s(X_{\\text{test}},y) \\le \\hat{q}\\}``
-	""")
+    ncal = Int(conf_model.train_ratio * data_specs.N)
+    Markdown.parse(
+        """
+The empirical coverage rate should be close to the desired level of coverage. In most cases it will be slightly higher, since ``(1-\\alpha)`` is a lower bound.
+
+> Found an empirical coverage rate that is slightly lower than desired? The coverage property is "marginal" in the sense that the probability averaged over the randomness in the data. For most purposes a large enough calibration set size (``n>1000``) mitigates that randomness enough. Depending on your choices above, the calibration set may be quite small (currently $ncal), which can lead to **coverage slack** (see Section 3 in the [tutorial](https://arxiv.org/pdf/2107.07511.pdf)).
+
+### *So what's happening under the hood?*
+
+Inductive Conformal Prediction (also referred to as Split Conformal Prediction) broadly speaking works as follows:
+
+1. Partition the training into a proper training set and a separate calibration set
+2. Train the machine learning model on the proper training set.
+3. Using some heuristic notion of uncertainty (e.g. absolute error in the regression case) compute nonconformity scores using the calibration data and the fitted model. 
+4. For the given coverage ratio compute the corresponding quantile of the empirical distribution of nonconformity scores.
+5. For the given quantile and test sample ``X_{\\text{test}}``, form the corresponding conformal prediction set like so: ``C(X_{\\text{test}})=\\{y:s(X_{\\text{test}},y) \\le \\hat{q}\\}``
+""",
+    )
 end
 
 # ╔═╡ 74444c01-1a0a-47a7-9b14-749946614f07
@@ -267,14 +298,25 @@ Quite cool, right? Using a single API call we are able to generate rigorous pred
 
 
 # ╔═╡ 824bd383-2fcb-4888-8ad1-260c85333edf
-@bind xmax_ood Slider(data_specs.xmax:(data_specs.xmax+5), default=(data_specs.xmax), show_value=true)
+@bind xmax_ood Slider(
+    data_specs.xmax:(data_specs.xmax+5),
+    default = (data_specs.xmax),
+    show_value = true,
+)
 
 # ╔═╡ 072cc72d-20a2-4ee9-954c-7ea70dfb8eea
 begin
-	Xood, yood = get_data(data_specs.N, xmax_ood, data_specs.noise; fun=f)
-	plot(_mach.model, _mach.fitresult, Xood, yood, zoom=0, observed_lab="Test points")
-	xood_range = range(-xmax_ood,xmax_ood,length=50)
-	plot!(xood_range, @.(f(xood_range)), lw=2, ls=:dash, colour=:black, label="Ground truth")
+    Xood, yood = get_data(data_specs.N, xmax_ood, data_specs.noise; fun = f)
+    plot(_mach.model, _mach.fitresult, Xood, yood, zoom = 0, observed_lab = "Test points")
+    xood_range = range(-xmax_ood, xmax_ood, length = 50)
+    plot!(
+        xood_range,
+        @.(f(xood_range)),
+        lw = 2,
+        ls = :dash,
+        colour = :black,
+        label = "Ground truth",
+    )
 end
 
 # ╔═╡ 4f41ec7c-aedd-475f-942d-33e2d1174902
diff --git a/dev/quick_tour/utils.jl b/dev/quick_tour/utils.jl
index ce9350b..979579f 100644
--- a/dev/quick_tour/utils.jl
+++ b/dev/quick_tour/utils.jl
@@ -1,19 +1,19 @@
-function multi_slider(vals::Dict; title="")
-    
+function multi_slider(vals::Dict; title = "")
+
     return PlutoUI.combine() do Child
-        
+
         inputs = [
             md""" $(_name): $(
                 Child(_name, Slider(_vals[1], default=_vals[2], show_value=true))
             )"""
-            
+
             for (_name, _vals) in vals
         ]
-        
+
         md"""
         #### $title
         $(inputs)
         """
     end
-    
-end
\ No newline at end of file
+
+end
diff --git a/src/conformal_models/plotting.jl b/src/conformal_models/plotting.jl
index 0dc09d0..248d407 100644
--- a/src/conformal_models/plotting.jl
+++ b/src/conformal_models/plotting.jl
@@ -59,13 +59,13 @@ function Plots.contourf(
     fitresult,
     X,
     y;
-    target::Union{Nothing,Real}=nothing,
-    ntest=50,
-    zoom=-1,
-    xlims=nothing,
-    ylims=nothing,
-    plot_set_size=false,
-    kwargs...
+    target::Union{Nothing,Real} = nothing,
+    ntest = 50,
+    zoom = -1,
+    xlims = nothing,
+    ylims = nothing,
+    plot_set_size = false,
+    kwargs...,
 )
 
     # Setup:
@@ -80,8 +80,8 @@ function Plots.contourf(
     xlims, ylims = generate_lims(x1, x2, xlims, ylims, zoom)
 
     # Surface range:
-    x1range = range(xlims[1], stop=xlims[2], length=ntest)
-    x2range = range(ylims[1], stop=ylims[2], length=ntest)
+    x1range = range(xlims[1], stop = xlims[2], length = ntest)
+    x2range = range(ylims[1], stop = ylims[2], length = ntest)
 
     # Target
     if !isnothing(target)
@@ -106,7 +106,7 @@ function Plots.contourf(
         if plot_set_size
             z = ismissing(p̂) ? 0 : sum(pdf.(p̂, p̂.decoder.classes) .> 0)
         else
-            z = ismissing(p̂) ? [missing for i in 1:length(levels(y))] : pdf.(p̂, levels(y))
+            z = ismissing(p̂) ? [missing for i = 1:length(levels(y))] : pdf.(p̂, levels(y))
             z = replace(z, 0 => missing)
         end
         push!(Z, z)
@@ -122,12 +122,12 @@ function Plots.contourf(
             x1range,
             x2range,
             Z;
-            title=title,
-            xlims=xlims,
-            ylims=ylims,
-            c=cgrad(:blues, _n + 1, categorical=true),
-            clim=clim,
-            kwargs...
+            title = title,
+            xlims = xlims,
+            ylims = ylims,
+            c = cgrad(:blues, _n + 1, categorical = true),
+            clim = clim,
+            kwargs...,
         )
     else
         clim = @isdefined(clim) ? clim : (0, 1)
@@ -135,17 +135,17 @@ function Plots.contourf(
             x1range,
             x2range,
             Z;
-            title=title,
-            xlims=xlims,
-            ylims=ylims,
-            clim=clim,
-            kwargs...
+            title = title,
+            xlims = xlims,
+            ylims = ylims,
+            clim = clim,
+            kwargs...,
         )
     end
 
     # Samples:
     y = typeof(y) <: CategoricalArrays.CategoricalArray ? y : Int.(y)
-    scatter!(plt, x1, x2, group=y; kwargs...)
+    scatter!(plt, x1, x2, group = y; kwargs...)
 
 end
 
@@ -159,9 +159,12 @@ end
 A `Plots.jl` recipe/method extension that can be used to visualize the conformal predictions of any fitted conformal classifier. Using a stacked area chart, this function plots the softmax output(s) contained the the conformal predictions set on the vertical axis against an input variable `X` on the horizontal axis. In the case of multiple input variables, the `input_var` argument can be used to specify the desired input variable.
 """
 function Plots.areaplot(
-    conf_model::ConformalProbabilisticSet, fitresult, X, y;
-    input_var::Union{Nothing,Int,Symbol}=nothing,
-    kwargs...
+    conf_model::ConformalProbabilisticSet,
+    fitresult,
+    X,
+    y;
+    input_var::Union{Nothing,Int,Symbol} = nothing,
+    kwargs...,
 )
 
     # Setup:
@@ -191,7 +194,8 @@ function Plots.areaplot(
     # Predictions:
     ŷ = predict(conf_model, fitresult, Xraw)
     nout = length(levels(y))
-    ŷ = map(_y -> ismissing(_y) ? [0 for i in 1:nout] : pdf.(_y, levels(y)), ŷ) |> _y -> reduce(hcat, _y)
+    ŷ =
+        map(_y -> ismissing(_y) ? [0 for i = 1:nout] : pdf.(_y, levels(y)), ŷ) |> _y -> reduce(hcat, _y)
     ŷ = permutedims(ŷ)
 
     areaplot(x, ŷ; kwargs...)
@@ -207,15 +211,18 @@ end
 A `Plots.jl` recipe/method extension that can be used to visualize the conformal predictions of a fitted conformal regressor. Data (`X`,`y`) are plotted as dots and overlaid with predictions intervals. `y` is plotted on the vertical axis against a single variable `X` on the horizontal axis. A shaded area indicates the prediction interval. The line in the center of the interval is the midpoint of the interval and can be interpreted as the point estimate of the conformal regressor. In case `X` is multi-dimensional, `input_var` can be used to specify the input variable of interest that will be used for the horizontal axis. If unspecified, the first variable will be plotting by default.
 """
 function Plots.plot(
-    conf_model::ConformalInterval, fitresult, X, y;
-    input_var::Union{Nothing,Int,Symbol}=nothing,
-    xlims::Union{Nothing,Tuple}=nothing,
-    ylims::Union{Nothing,Tuple}=nothing,
-    zoom::Real=-0.5,
-    train_lab::Union{Nothing,String}=nothing,
-    test_lab::Union{Nothing,String}=nothing,
-    ymid_lw::Int=1,
-    kwargs...
+    conf_model::ConformalInterval,
+    fitresult,
+    X,
+    y;
+    input_var::Union{Nothing,Int,Symbol} = nothing,
+    xlims::Union{Nothing,Tuple} = nothing,
+    ylims::Union{Nothing,Tuple} = nothing,
+    zoom::Real = -0.5,
+    train_lab::Union{Nothing,String} = nothing,
+    test_lab::Union{Nothing,String} = nothing,
+    ymid_lw::Int = 1,
+    kwargs...,
 )
 
     # Setup
@@ -253,11 +260,11 @@ function Plots.plot(
     plt = scatter(
         vec(x),
         vec(y),
-        label=train_lab,
-        xlim=xlims,
-        ylim=ylims,
-        title=title;
-        kwargs...
+        label = train_lab,
+        xlim = xlims,
+        ylim = ylims,
+        title = title;
+        kwargs...,
     )
 
     # Plot predictions:
@@ -267,7 +274,15 @@ function Plots.plot(
     yerror = (ub .- lb) ./ 2
     xplot = vec(x)
     _idx = sortperm(xplot)
-    plot!(plt, xplot[_idx], ymid[_idx], label=test_lab, ribbon=(yerror, yerror), lw=ymid_lw; kwargs...)
+    plot!(
+        plt,
+        xplot[_idx],
+        ymid[_idx],
+        label = test_lab,
+        ribbon = (yerror, yerror),
+        lw = ymid_lw;
+        kwargs...,
+    )
 
 end
 
@@ -276,13 +291,17 @@ end
 
 A `Plots.jl` recipe/method extension that can be used to visualize the set size distribution of a conformal predictor. In the regression case, prediction interval widths are stratified into discrete bins. It can be useful to plot the distribution of set sizes in order to visually asses how adaptive a conformal predictor is. For more adaptive predictors the distribution of set sizes is typically spread out more widely, which reflects that “the procedure is effectively distinguishing between easy and hard inputs”. This is desirable: when for a given sample it is difficult to make predictions, this should be reflected in the set size (or interval width in the regression case). Since ‘difficult’ lies on some spectrum that ranges from ‘very easy’ to ‘very difficult’ the set size should very across the spectrum of ‘empty set’ to ‘all labels included’.
 """
-function Plots.bar(conf_model::ConformalModel, fitresult, X; label="", xtickfontsize=6, kwrgs...)
+function Plots.bar(
+    conf_model::ConformalModel,
+    fitresult,
+    X;
+    label = "",
+    xtickfontsize = 6,
+    kwrgs...,
+)
     ŷ = predict(conf_model, fitresult, X)
     idx = size_indicator(ŷ)
-    x = sort(levels(idx), lt=natural)
+    x = sort(levels(idx), lt = natural)
     y = [sum(idx .== _x) for _x in x]
-    Plots.bar(x, y; label=label, xtickfontsize=xtickfontsize, kwrgs...)
+    Plots.bar(x, y; label = label, xtickfontsize = xtickfontsize, kwrgs...)
 end
-
-
-
diff --git a/src/conformal_models/utils.jl b/src/conformal_models/utils.jl
index 0882566..af55cc9 100644
--- a/src/conformal_models/utils.jl
+++ b/src/conformal_models/utils.jl
@@ -53,7 +53,7 @@ function set_size(ŷ)
     return _size
 end
 
-function size_indicator(ŷ::AbstractVector; bins = 5, tol=1e-10)
+function size_indicator(ŷ::AbstractVector; bins = 5, tol = 1e-10)
 
     _sizes = set_size.(ŷ)
     unique_sizes = unique(_sizes)
@@ -63,22 +63,24 @@ function size_indicator(ŷ::AbstractVector; bins = 5, tol=1e-10)
         if abs.(diff(collect(extrema(set_size.(ŷ)))))[1] < tol
             idx = categorical(ones(length(_sizes)))
         else
-            bin_caps = collect(range(minimum(unique_sizes), maximum(unique_sizes), length=bins+1))[2:end]
+            bin_caps = collect(
+                range(minimum(unique_sizes), maximum(unique_sizes), length = bins + 1),
+            )[2:end]
             idx = map(_sizes) do s
                 # Check which is the largest bin cap that _size exceeds:
-                ub = argmax(x -> s-x <= 0 ? s-x : -Inf, bin_caps)
+                ub = argmax(x -> s - x <= 0 ? s - x : -Inf, bin_caps)
                 if ub == minimum(bin_caps)
-                    ub = round(ub, digits=2)
-                    lb = round(minimum(_sizes), digits=2)
+                    ub = round(ub, digits = 2)
+                    lb = round(minimum(_sizes), digits = 2)
                     _idx = "|C| ∈ ($lb,$ub]"
                 else
-                    ub = round(ub, digits=2)
-                    lb = round(argmin(x -> s-x > 0 ? s-x : Inf, bin_caps),digits=2)
+                    ub = round(ub, digits = 2)
+                    lb = round(argmin(x -> s - x > 0 ? s - x : Inf, bin_caps), digits = 2)
                     _idx = "|C| ∈ ($lb,$ub]"
                 end
-                return _idx 
+                return _idx
             end
-            idx = categorical(idx)          
+            idx = categorical(idx)
         end
     end
 
@@ -87,7 +89,7 @@ function size_indicator(ŷ::AbstractVector; bins = 5, tol=1e-10)
         bin_caps = collect(1:2:(maximum(unique_sizes)+1))
         idx = map(_sizes) do s
             # Check which is the largest bin cap that _size exceeds:
-            ub = bin_caps[sum(s .> bin_caps) + 1]
+            ub = bin_caps[sum(s .> bin_caps)+1]
             if ub > maximum(_sizes)
                 ub = ub - 1
                 _idx = "|C| ∈ [$ub]"
@@ -97,7 +99,7 @@ function size_indicator(ŷ::AbstractVector; bins = 5, tol=1e-10)
             end
             return _idx
         end
-        idx = categorical(idx)   
+        idx = categorical(idx)
     end
 
     return idx
diff --git a/src/evaluation/measures.jl b/src/evaluation/measures.jl
index 4a8512b..a730268 100644
--- a/src/evaluation/measures.jl
+++ b/src/evaluation/measures.jl
@@ -19,7 +19,7 @@ function size_stratified_coverage(ŷ, y)
     # Setup:
     stratum_indicator = size_indicator(ŷ) |> x -> x.refs
     unique_stratums = sort(unique(stratum_indicator))
-    unique_stratums = unique_stratums[unique_stratums .!= 0]
+    unique_stratums = unique_stratums[unique_stratums.!=0]
     _covs = []
 
     if length(unique_stratums) == 1 && is_regression(ŷ)
@@ -38,4 +38,4 @@ function size_stratified_coverage(ŷ, y)
 
     return C̄
 
-end
\ No newline at end of file
+end
diff --git a/src/evaluation/utils.jl b/src/evaluation/utils.jl
index 41b0b2e..8bcb101 100644
--- a/src/evaluation/utils.jl
+++ b/src/evaluation/utils.jl
@@ -43,4 +43,4 @@ function is_covered(ŷ, y)
         return _is_covered
     end
     return is_covered
-end
\ No newline at end of file
+end
diff --git a/test/classification.jl b/test/classification.jl
index 4b957ee..b77650a 100644
--- a/test/classification.jl
+++ b/test/classification.jl
@@ -5,18 +5,14 @@ using Plots
 data_specs = (
     "Single Input - Binary" => (1, 2),
     "Two Inputs - Binary" => (2, 2),
-    "Multiple Inputs - Multiclass" => (5,5),
+    "Multiple Inputs - Multiclass" => (5, 5),
 )
-data_sets = Dict{String, Any}()
-for (k,v) in data_specs
+data_sets = Dict{String,Any}()
+for (k, v) in data_specs
     X, y = MLJ.make_blobs(500, v[1], centers = v[2])
     X = MLJ.table(MLJ.matrix(X))
     train, test = partition(eachindex(y), 0.8)
-    _set = Dict(
-        :data => (X,y),
-        :split => (train, test),
-        :specs => v
-    )
+    _set = Dict(:data => (X, y), :split => (train, test), :specs => v)
     data_sets[k] = _set
 end
 
@@ -63,26 +59,47 @@ conformal_models = merge(values(available_models[:classification])...)
                             # Plotting:
                             @test isplot(bar(mach.model, mach.fitresult, X))
                             @test isplot(areaplot(mach.model, mach.fitresult, X, y))
-                            @test isplot(areaplot(mach.model, mach.fitresult, X, y; input_var=1))
-                            @test isplot(areaplot(mach.model, mach.fitresult, X, y; input_var=:x1))
+                            @test isplot(
+                                areaplot(mach.model, mach.fitresult, X, y; input_var = 1),
+                            )
+                            @test isplot(
+                                areaplot(mach.model, mach.fitresult, X, y; input_var = :x1),
+                            )
                             if data_set[:specs][1] != 2
-                                @test_throws AssertionError contourf(mach.model, mach.fitresult, X, y)
+                                @test_throws AssertionError contourf(
+                                    mach.model,
+                                    mach.fitresult,
+                                    X,
+                                    y,
+                                )
                             else
                                 @test isplot(contourf(mach.model, mach.fitresult, X, y))
-                                @test isplot(contourf(mach.model, mach.fitresult, X, y; zoom = -1, plot_set_size = true))
-                                @test isplot(contourf(mach.model, mach.fitresult, X, y; target=1))
+                                @test isplot(
+                                    contourf(
+                                        mach.model,
+                                        mach.fitresult,
+                                        X,
+                                        y;
+                                        zoom = -1,
+                                        plot_set_size = true,
+                                    ),
+                                )
+                                @test isplot(
+                                    contourf(mach.model, mach.fitresult, X, y; target = 1),
+                                )
                             end
-                            
+
                             # Evaluation:
                             # Evaluation takes some time, so only testing for one method.
                             if _method == :simple_inductive && data_set[:specs][1] > 1
                                 # Empirical coverage:
-                                _eval = evaluate!(mach; measure = emp_coverage, verbosity = 0)
+                                _eval =
+                                    evaluate!(mach; measure = emp_coverage, verbosity = 0)
                                 Δ = _eval.measurement[1] - _cov     # over-/under-coverage
                                 @test Δ >= -0.05                    # we don't undercover too much
                                 # Size-stratified coverage:
                                 _eval = evaluate!(mach; measure = ssc, verbosity = 0)
-                            end 
+                            end
 
                         end
 
diff --git a/test/regression.jl b/test/regression.jl
index 90b73cf..80db55c 100644
--- a/test/regression.jl
+++ b/test/regression.jl
@@ -2,20 +2,14 @@ using MLJ
 using Plots
 
 # Data:
-data_specs = (
-    "Single Input - Single Output" => (1, 1),
-    "Multiple Inputs - Single Output" => (5, 1),
-)
+data_specs =
+    ("Single Input - Single Output" => (1, 1), "Multiple Inputs - Single Output" => (5, 1))
 data_sets = Dict{String,Any}()
 for (k, v) in data_specs
     X, y = MLJ.make_regression(500, v[1])
     X = MLJ.table(MLJ.matrix(X))
     train, test = partition(eachindex(y), 0.8)
-    _set = Dict(
-        :data => (X, y),
-        :split => (train, test),
-        :specs => v
-    )
+    _set = Dict(:data => (X, y), :split => (train, test), :specs => v)
     data_sets[k] = _set
 end
 
@@ -40,7 +34,7 @@ conformal_models = merge(values(available_models[:regression])...)
 
                     # Instantiate conformal models:
                     _cov = 0.85
-                    conf_model = conformal_model(model; method=_method, coverage=_cov)
+                    conf_model = conformal_model(model; method = _method, coverage = _cov)
                     conf_model = conformal_models[_method](model)
                     @test isnothing(conf_model.scores)
 
@@ -54,25 +48,38 @@ conformal_models = merge(values(available_models[:regression])...)
 
                             # Fit/Predict:
                             mach = machine(conf_model, X, y)
-                            fit!(mach, rows=train)
+                            fit!(mach, rows = train)
                             @test !isnothing(conf_model.scores)
                             predict(mach, selectrows(X, test))
 
                             # Plotting:
                             @test isplot(plot(mach.model, mach.fitresult, X, y))
-                            @test isplot(plot(mach.model, mach.fitresult, X, y; input_var=1, xlims=(-1,1), ylims=(-1,1)))
-                            @test isplot(plot(mach.model, mach.fitresult, X, y; input_var=:x1))
+                            @test isplot(
+                                plot(
+                                    mach.model,
+                                    mach.fitresult,
+                                    X,
+                                    y;
+                                    input_var = 1,
+                                    xlims = (-1, 1),
+                                    ylims = (-1, 1),
+                                ),
+                            )
+                            @test isplot(
+                                plot(mach.model, mach.fitresult, X, y; input_var = :x1),
+                            )
                             @test isplot(bar(mach.model, mach.fitresult, X))
 
                             # Evaluation:
                             # Evaluation takes some time, so only testing for one method.
                             if _method == :simple_inductive
                                 # Empirical coverage:
-                                _eval = evaluate!(mach; measure=emp_coverage, verbosity=0)
+                                _eval =
+                                    evaluate!(mach; measure = emp_coverage, verbosity = 0)
                                 Δ = _eval.measurement[1] - _cov     # over-/under-coverage
                                 @test Δ >= -0.05                    # we don't undercover too much
                                 # Size-stratified coverage:
-                                _eval = evaluate!(mach; measure=ssc, verbosity=0)
+                                _eval = evaluate!(mach; measure = ssc, verbosity = 0)
                             end
 
                         end
diff --git a/test/utils.jl b/test/utils.jl
index 4876688..1c3a400 100644
--- a/test/utils.jl
+++ b/test/utils.jl
@@ -1,3 +1,3 @@
 using Plots
 
-isplot(plt) = typeof(plt) <: Plots.Plot
\ No newline at end of file
+isplot(plt) = typeof(plt) <: Plots.Plot