diff --git a/.JuliaFormatter.toml b/.JuliaFormatter.toml
new file mode 100644
index 000000000..ada0e45ce
--- /dev/null
+++ b/.JuliaFormatter.toml
@@ -0,0 +1 @@
+style="sciml"
diff --git a/.github/actions/install-julia/action.yaml b/.github/actions/install-julia/action.yaml
deleted file mode 100644
index 5456288f2..000000000
--- a/.github/actions/install-julia/action.yaml
+++ /dev/null
@@ -1,40 +0,0 @@
-name: Install Julia
-description: Installs a user-specified version of Julia
-inputs:
-  version:
-    description: 'The version of Julia to install (e.g., 1.10.0)'
-    required: true
-    default: '1.10.0'
-runs:
-  using: composite
-  steps:
-    - name: Cache Julia binary
-      id: cache-julia
-      uses: actions/cache@v3
-      with:
-        path: ~/julia
-        key: julia-binary-${{ inputs.version }}
-        restore-keys: |
-          julia-binary-${{ inputs.version }}
-
-    - name: Install Julia
-      if: steps.cache-julia.outputs.cache-hit != 'true'
-      run: |
-        mkdir -p ~/julia
-        version=$(echo ${{ inputs.version }} | cut -d. -f1,2)
-        wget https://julialang-s3.julialang.org/bin/linux/x64/$version/julia-${{ inputs.version }}-linux-x86_64.tar.gz -O ~/julia/julia.tar.gz
-      shell: bash
-
-    - name: Unzip Julia
-      run: tar -xzf ~/julia/julia.tar.gz -C ~/julia --strip-components=1
-      shell: bash
-
-    - run: echo "$HOME/julia/bin" >> $GITHUB_PATH
-      shell: bash
-
-    - uses: actions/cache@v3
-      with:
-        path: |
-          ~/.julia
-          ~/.cache/julia
-        key: julia-${{ inputs.version }}-${{ hashFiles('**/*.toml') }}
diff --git a/.github/workflows/CompatHelper-EpiAware.yaml b/.github/workflows/CompatHelper-EpiAware.yaml
new file mode 100644
index 000000000..83f8af07f
--- /dev/null
+++ b/.github/workflows/CompatHelper-EpiAware.yaml
@@ -0,0 +1,30 @@
+name: Compatibility Helper
+
+on:
+  schedule:
+    - cron: '00 * * * *'
+  issues:
+    types: [opened, reopened]
+
+concurrency:
+  group: ${{ github.workflow }}-${{ github.ref }}
+  cancel-in-progress: true
+
+jobs:
+  compatibility:
+    runs-on: ubuntu-latest
+    steps:
+      - uses: actions/checkout@v4
+        with:
+          sparse-checkout: EpiAware
+          sparse-checkout-cone-mode: false
+      - name: Move EpiAware to root
+        run: mv EpiAware/* .
+      - uses: julia-actions/setup-julia@v1
+      - uses: julia-actions/cache@v1
+      - name: Pkg.add("CompatHelper")
+        run: julia -e 'using Pkg; Pkg.add("CompatHelper")'
+      - name: CompatHelper.main()
+        env:
+          GITHUB_TOKEN: ${{ secrets.GITHUB_TOKEN }}
+        run: julia -e 'using CompatHelper; CompatHelper.main()'
diff --git a/.github/workflows/codecoverage-EpiAware.yaml b/.github/workflows/codecoverage-EpiAware.yaml
new file mode 100644
index 000000000..5a9e263af
--- /dev/null
+++ b/.github/workflows/codecoverage-EpiAware.yaml
@@ -0,0 +1,34 @@
+name: Code coverage
+
+on:
+  push:
+    branches:
+      - main
+  pull_request:
+
+concurrency:
+  group: ${{ github.workflow }}-${{ github.ref }}
+  cancel-in-progress: true
+
+jobs:
+  coverage:
+    runs-on: ubuntu-latest
+    steps:
+      - uses: actions/checkout@v4
+        with:
+          sparse-checkout: EpiAware
+          sparse-checkout-cone-mode: false
+      - name: Move EpiAware to root
+        run: mv EpiAware/* .
+      - uses: julia-actions/setup-julia@v1
+      - uses: julia-actions/cache@v1
+      - uses: julia-actions/julia-buildpkg@v1
+      - uses: julia-actions/julia-runtest@v1
+        with:
+          annotate: true
+      - uses: julia-actions/julia-processcoverage@v1
+      - uses: codecov/codecov-action@v4
+        with:
+          file: lcov.info
+          token: ${{ secrets.CODECOV_TOKEN }}
+          fail_ci_if_error: true
diff --git a/.github/workflows/pre-commit.yaml b/.github/workflows/pre-commit.yaml
index 7aa43557e..3deb88555 100644
--- a/.github/workflows/pre-commit.yaml
+++ b/.github/workflows/pre-commit.yaml
@@ -11,8 +11,7 @@ jobs:
     steps:
     - uses: actions/checkout@v3
     - uses: actions/setup-python@v3
-    - uses: ./.github/actions/install-julia
-      with:
-        version: '1.10.0'
+    - uses: julia-actions/setup-julia@v1
+    - uses: julia-actions/cache@v1
     - run: julia -e 'using Pkg; Pkg.add("JuliaFormatter")'
     - uses: ./.github/actions/pre-commit
diff --git a/.github/workflows/test-EpiAware.yaml b/.github/workflows/test-EpiAware.yaml
index 4e615c2c9..fa84c6389 100644
--- a/.github/workflows/test-EpiAware.yaml
+++ b/.github/workflows/test-EpiAware.yaml
@@ -1,17 +1,33 @@
 name: Test EpiAware
 
 on:
-  pull_request:
   push:
-    branches: [main]
+    branches:
+      - main
+  pull_request:
+
+concurrency:
+  group: ${{ github.workflow }}-${{ github.ref }}
+  cancel-in-progress: true
 
 jobs:
-  test-EpiAware:
-    runs-on: ubuntu-latest
+  test:
+    runs-on: ${{ matrix.os }}
+    strategy:
+      fail-fast: true
+      matrix:
+        julia-version: ['1']
+        os: [ubuntu-latest, windows-latest, macOS-latest]
     steps:
-    - uses: actions/checkout@v3
-    - uses: ./.github/actions/install-julia
-      with:
-        version: '1.10.0'
-    - name: Run unit tests for EpiAware
-      run: julia --project=EpiAware -e 'using Pkg; Pkg.test(coverage = true)'
+      - uses: actions/checkout@v4
+        with:
+          sparse-checkout: EpiAware
+          sparse-checkout-cone-mode: false
+      - name: Move EpiAware to root
+        run: mv EpiAware/* .
+      - uses: julia-actions/setup-julia@v1
+        with:
+          version: ${{ matrix.julia-version }}
+      - uses: julia-actions/cache@v1
+      - uses: julia-actions/julia-buildpkg@v1
+      - uses: julia-actions/julia-runtest@v1
diff --git a/EpiAware/.JuliaFormatter.toml b/EpiAware/.JuliaFormatter.toml
new file mode 100644
index 000000000..ada0e45ce
--- /dev/null
+++ b/EpiAware/.JuliaFormatter.toml
@@ -0,0 +1 @@
+style="sciml"
diff --git a/EpiAware/README.md b/EpiAware/README.md
index 5e699c1b5..32d6b8d6d 100644
--- a/EpiAware/README.md
+++ b/EpiAware/README.md
@@ -1,5 +1,8 @@
 # EpiAware
 
+[![SciML Code Style](https://img.shields.io/static/v1?label=code%20style&message=SciML&color=9558b2&labelColor=389826)](https://github.com/SciML/SciMLStyle)
+[![Test EpiAware](https://github.com/CDCgov/Rt-without-renewal/actions/workflows/test-EpiAware.yaml/badge.svg)](https://github.com/CDCgov/Rt-without-renewal/actions/workflows/test-EpiAware.yaml)
+[![codecov](https://codecov.io/gh/CDCgov/Rt-without-renewal/graph/badge.svg?token=IX4GIA8F0H)](https://codecov.io/gh/CDCgov/Rt-without-renewal)
 
 ## Model Diagram
 
diff --git a/EpiAware/src/EpiAware.jl b/EpiAware/src/EpiAware.jl
index 01dc5bb67..f9235431e 100644
--- a/EpiAware/src/EpiAware.jl
+++ b/EpiAware/src/EpiAware.jl
@@ -20,16 +20,16 @@ This module provides functionality for calculating Rt (effective reproduction nu
 module EpiAware
 
 using Distributions,
-    Turing,
-    LogExpFunctions,
-    LinearAlgebra,
-    SparseArrays,
-    Random,
-    ReverseDiff,
-    Optim,
-    Parameters,
-    QuadGK,
-    DataFramesMeta
+      Turing,
+      LogExpFunctions,
+      LinearAlgebra,
+      SparseArrays,
+      Random,
+      ReverseDiff,
+      Optim,
+      Parameters,
+      QuadGK,
+      DataFramesMeta
 
 # Exported utilities
 export create_discrete_pmf, default_rw_priors, default_delay_obs_priors, spread_draws
diff --git a/EpiAware/src/epimodel.jl b/EpiAware/src/epimodel.jl
index 7a39c88f6..67929d821 100644
--- a/EpiAware/src/epimodel.jl
+++ b/EpiAware/src/epimodel.jl
@@ -1,9 +1,9 @@
 abstract type AbstractEpiModel end
 
-struct EpiData{T<:Real,F<:Function}
+struct EpiData{T <: Real, F <: Function}
     gen_int::Vector{T}
     delay_int::Vector{T}
-    delay_kernel::SparseMatrixCSC{T,Integer}
+    delay_kernel::SparseMatrixCSC{T, Integer}
     cluster_coeff::T
     len_gen_int::Integer
     len_delay_int::Integer
@@ -12,20 +12,20 @@ struct EpiData{T<:Real,F<:Function}
 
     #Inner constructors for EpiData object
     function EpiData(
-        gen_int,
-        delay_int,
-        cluster_coeff,
-        time_horizon::Integer,
-        transformation::Function,
+            gen_int,
+            delay_int,
+            cluster_coeff,
+            time_horizon::Integer,
+            transformation::Function
     )
         @assert all(gen_int .>= 0) "Generation interval must be non-negative"
         @assert all(delay_int .>= 0) "Delay interval must be non-negative"
-        @assert sum(gen_int) ≈ 1 "Generation interval must sum to 1"
-        @assert sum(delay_int) ≈ 1 "Delay interval must sum to 1"
+        @assert sum(gen_int)≈1 "Generation interval must sum to 1"
+        @assert sum(delay_int)≈1 "Delay interval must sum to 1"
 
         K = generate_observation_kernel(delay_int, time_horizon)
 
-        new{eltype(gen_int),typeof(transformation)}(
+        new{eltype(gen_int), typeof(transformation)}(
             gen_int,
             delay_int,
             K,
@@ -33,23 +33,22 @@ struct EpiData{T<:Real,F<:Function}
             length(gen_int),
             length(delay_int),
             time_horizon,
-            transformation,
+            transformation
         )
     end
 
     function EpiData(
-        gen_distribution::ContinuousDistribution,
-        delay_distribution::ContinuousDistribution,
-        cluster_coeff,
-        time_horizon::Integer;
-        D_gen,
-        D_delay,
-        Δd = 1.0,
-        transformation::Function = exp,
+            gen_distribution::ContinuousDistribution,
+            delay_distribution::ContinuousDistribution,
+            cluster_coeff,
+            time_horizon::Integer;
+            D_gen,
+            D_delay,
+            Δd = 1.0,
+            transformation::Function = exp
     )
-        gen_int =
-            create_discrete_pmf(gen_distribution, Δd = Δd, D = D_gen) |>
-            p -> p[2:end] ./ sum(p[2:end])
+        gen_int = create_discrete_pmf(gen_distribution, Δd = Δd, D = D_gen) |>
+                  p -> p[2:end] ./ sum(p[2:end])
         delay_int = create_discrete_pmf(delay_distribution, Δd = Δd, D = D_delay)
 
         return EpiData(gen_int, delay_int, cluster_coeff, time_horizon, transformation)
@@ -81,11 +80,11 @@ function (epimodel::Renewal)(_Rt, init)
     Rt = epimodel.data.transformation.(_Rt)
 
     r_approx = R_to_r(Rt[1], epimodel)
-    init = I₀ * [exp(-r_approx * t) for t = 0:(epimodel.data.len_gen_int-1)]
+    init = I₀ * [exp(-r_approx * t) for t in 0:(epimodel.data.len_gen_int - 1)]
 
     function generate_infs(recent_incidence, Rt)
         new_incidence = Rt * dot(recent_incidence, epimodel.data.gen_int)
-        [new_incidence; recent_incidence[1:(epimodel.data.len_gen_int-1)]], new_incidence
+        [new_incidence; recent_incidence[1:(epimodel.data.len_gen_int - 1)]], new_incidence
     end
 
     I_t, _ = scan(generate_infs, init, Rt)
diff --git a/EpiAware/src/latent-processes.jl b/EpiAware/src/latent-processes.jl
index 2ac2de53b..eb43bd332 100644
--- a/EpiAware/src/latent-processes.jl
+++ b/EpiAware/src/latent-processes.jl
@@ -1,7 +1,7 @@
 function default_rw_priors()
     return (
         var_RW_dist = truncated(Normal(0.0, 0.05), 0.0, Inf),
-        init_rw_value_dist = Normal(),
+        init_rw_value_dist = Normal()
     )
 end
 
@@ -13,8 +13,8 @@ end
     rw = Vector{eltype(ϵ_t)}(undef, n)
 
     rw[1] = σ_RW * ϵ_t[1]
-    for t = 2:n
-        rw[t] = rw[t-1] + σ_RW * ϵ_t[t]
+    for t in 2:n
+        rw[t] = rw[t - 1] + σ_RW * ϵ_t[t]
     end
     return rw, init, (; σ_RW,)
 end
diff --git a/EpiAware/src/models.jl b/EpiAware/src/models.jl
index 24bf036fc..9396e6fd3 100644
--- a/EpiAware/src/models.jl
+++ b/EpiAware/src/models.jl
@@ -1,15 +1,15 @@
 @model function make_epi_inference_model(
-    y_t,
-    epimodel::AbstractEpiModel,
-    latent_process,
-    observation_process;
-    process_priors,
-    pos_shift = 1e-6,
+        y_t,
+        epimodel::AbstractEpiModel,
+        latent_process,
+        observation_process;
+        process_priors,
+        pos_shift = 1e-6
 )
     #Latent process
     time_steps = epimodel.data.time_horizon
-    @submodel latent_process, init, latent_process_aux =
-        latent_process(time_steps; latent_process_priors = process_priors)
+    @submodel latent_process, init, latent_process_aux = latent_process(
+        time_steps; latent_process_priors = process_priors)
 
     #Transform into infections
     I_t = epimodel(latent_process, init)
@@ -20,7 +20,7 @@
         I_t,
         epimodel::AbstractEpiModel;
         observation_process_priors = process_priors,
-        pos_shift = pos_shift,
+        pos_shift = pos_shift
     )
 
     #Generate quantities
@@ -28,6 +28,6 @@
         generated_y_t,
         I_t,
         latent_process,
-        process_aux = merge(latent_process_aux, generated_y_t_aux),
+        process_aux = merge(latent_process_aux, generated_y_t_aux)
     )
 end
diff --git a/EpiAware/src/observation-processes.jl b/EpiAware/src/observation-processes.jl
index 90d904449..379e0701b 100644
--- a/EpiAware/src/observation-processes.jl
+++ b/EpiAware/src/observation-processes.jl
@@ -3,19 +3,18 @@ function default_delay_obs_priors()
 end
 
 @model function delay_observations(
-    y_t,
-    I_t,
-    epimodel::AbstractEpiModel;
-    observation_process_priors = default_delay_obs_priors(),
-    pos_shift = 1e-6,
+        y_t,
+        I_t,
+        epimodel::AbstractEpiModel;
+        observation_process_priors = default_delay_obs_priors(),
+        pos_shift = 1e-6
 )
     #Parameters
     neg_bin_cluster_factor ~ observation_process_priors.neg_bin_cluster_factor_prior
 
     #Predictive distribution
-    case_pred_dists =
-        (epimodel.data.delay_kernel * I_t) .+ pos_shift .|>
-        μ -> mean_cc_neg_bin(μ, neg_bin_cluster_factor)
+    case_pred_dists = (epimodel.data.delay_kernel * I_t) .+ pos_shift .|>
+                      μ -> mean_cc_neg_bin(μ, neg_bin_cluster_factor)
 
     #Likelihood
     y_t ~ arraydist(case_pred_dists)
diff --git a/EpiAware/src/utilities.jl b/EpiAware/src/utilities.jl
index b27eaca75..5a2f9d035 100644
--- a/EpiAware/src/utilities.jl
+++ b/EpiAware/src/utilities.jl
@@ -17,7 +17,7 @@ value. This is similar to the JAX function `jax.lax.scan`.
 - `ys`: An array containing the result values of applying `f` to each element of `xs`.
 - `carry`: The final accumulator value.
 """
-function scan(f::Function, init, xs::Vector{T}) where {T<:Union{Integer,AbstractFloat}}
+function scan(f::Function, init, xs::Vector{T}) where {T <: Union{Integer, AbstractFloat}}
     carry = init
     ys = similar(xs)
     for (i, x) in enumerate(xs)
@@ -53,19 +53,19 @@ Raises:
 - `AssertionError` if `D` is not greater than `Δd`.
 """
 function create_discrete_pmf(
-    dist::Distribution,
-    ::Val{:single_censored};
-    primary_approximation_point = 0.5,
-    Δd = 1.0,
-    D,
+        dist::Distribution,
+        ::Val{:single_censored};
+        primary_approximation_point = 0.5,
+        Δd = 1.0,
+        D
 )
-    @assert minimum(dist) >= 0.0 "Distribution must be non-negative"
-    @assert Δd > 0.0 "Δd must be positive"
-    @assert D > Δd "D must be greater than Δd"
-    @assert primary_approximation_point >= 0.0 && primary_approximation_point <= 1.0 "`primary_approximation_point` must be in [0,1]."
+    @assert minimum(dist)>=0.0 "Distribution must be non-negative"
+    @assert Δd>0.0 "Δd must be positive"
+    @assert D>Δd "D must be greater than Δd"
+    @assert primary_approximation_point >= 0.0&&primary_approximation_point <= 1.0 "`primary_approximation_point` must be in [0,1]."
 
     ts = Δd:Δd:D |> collect
-    @assert ts[end] == D "D must be a multiple of Δd."
+    @assert ts[end]==D "D must be a multiple of Δd."
     ts = [primary_approximation_point * Δd; ts] #This covers situation where primary_approximation_point == 1
 
     ts .|> (t -> cdf(dist, t)) |> diff |> p -> p ./ sum(p)
@@ -96,13 +96,13 @@ Raises:
 - `AssertionError` if `D` is not greater than `Δd`.
 """
 function create_discrete_pmf(dist::Distribution; Δd = 1.0, D)
-    @assert minimum(dist) >= 0.0 "Distribution must be non-negative."
-    @assert Δd > 0.0 "Δd must be positive."
-    @assert D > Δd "D must be greater than Δd."
+    @assert minimum(dist)>=0.0 "Distribution must be non-negative."
+    @assert Δd>0.0 "Δd must be positive."
+    @assert D>Δd "D must be greater than Δd."
 
     ts = 0.0:Δd:D |> collect
 
-    @assert ts[end] == D "D must be a multiple of Δd."
+    @assert ts[end]==D "D must be a multiple of Δd."
 
     ∫F(dist, t, Δd) = quadgk(u -> cdf(dist, t - u) / Δd, 0.0, Δd)[1]
 
@@ -123,11 +123,11 @@ The value of the negative MGF.
 
 """
 function neg_MGF(r, w::AbstractVector)
-    return sum([w[i] * exp(-r * i) for i = 1:length(w)])
+    return sum([w[i] * exp(-r * i) for i in 1:length(w)])
 end
 
 function dneg_MGF_dr(r, w::AbstractVector)
-    return -sum([w[i] * i * exp(-r * i) for i = 1:length(w)])
+    return -sum([w[i] * i * exp(-r * i) for i in 1:length(w)])
 end
 
 """
@@ -156,12 +156,12 @@ The two step approximation is based on:
 Returns:
 - The approximate value of `r`.
 """
-function R_to_r(R₀, w::Vector{T}; newton_steps = 2, Δd = 1.0) where {T<:AbstractFloat}
+function R_to_r(R₀, w::Vector{T}; newton_steps = 2, Δd = 1.0) where {T <: AbstractFloat}
     mean_gen_time = dot(w, 1:length(w)) * Δd
     # Small r approximation as initial guess
     r_approx = (R₀ - 1) / (R₀ * mean_gen_time)
     # Newton's method
-    for _ = 1:newton_steps
+    for _ in 1:newton_steps
         r_approx -= (R₀ * neg_MGF(r_approx, w) - 1) / (R₀ * dneg_MGF_dr(r_approx, w))
     end
     return r_approx
@@ -171,7 +171,6 @@ function R_to_r(R₀, epimodel::AbstractEpiModel; newton_steps = 2, Δd = 1.0)
     R_to_r(R₀, epimodel.data.gen_int; newton_steps = newton_steps, Δd = Δd)
 end
 
-
 """
     r_to_R(r, w)
 
@@ -189,7 +188,6 @@ function r_to_R(r, w::AbstractVector)
     return 1 / neg_MGF(r, w::AbstractVector)
 end
 
-
 """
     mean_cc_neg_bin(μ, α)
 
@@ -209,7 +207,6 @@ function mean_cc_neg_bin(μ, α)
     return NegativeBinomial(r, p)
 end
 
-
 """
     generate_observation_kernel(delay_int, time_horizon)
 
@@ -224,10 +221,10 @@ Generate an observation kernel matrix based on the given delay interval and time
 """
 function generate_observation_kernel(delay_int, time_horizon)
     K = zeros(eltype(delay_int), time_horizon, time_horizon) |> SparseMatrixCSC
-    for i = 1:time_horizon, j = 1:time_horizon
+    for i in 1:time_horizon, j in 1:time_horizon
         m = i - j
         if m >= 0 && m <= (length(delay_int) - 1)
-            K[i, j] = delay_int[m+1]
+            K[i, j] = delay_int[m + 1]
         end
     end
     return K
@@ -247,9 +244,9 @@ Converts a `Chains` object into a DataFrame in `tidybayes` format.
 function spread_draws(chn::Chains)
     df = DataFrame(chn)
     df = hcat(DataFrame(draw = 1:size(df, 1)), df)
-    @rename!(df, $(".draw") = :draw)
-    @rename!(df, $(".chain") = :chain)
-    @rename!(df, $(".iteration") = :iteration)
+    @rename!(df, $(".draw")=:draw)
+    @rename!(df, $(".chain")=:chain)
+    @rename!(df, $(".iteration")=:iteration)
 
     return df
 end
diff --git a/EpiAware/test/predictive_checking/discretized_pmfs.jl b/EpiAware/test/predictive_checking/discretized_pmfs.jl
index bc4fb47ab..da6e71ef6 100644
--- a/EpiAware/test/predictive_checking/discretized_pmfs.jl
+++ b/EpiAware/test/predictive_checking/discretized_pmfs.jl
@@ -68,7 +68,7 @@ D = 21.0
 
 plt1 = let
     Δd = 1
-    ts = (0.0:Δd:(D-Δd)) |> collect
+    ts = (0.0:Δd:(D - Δd)) |> collect
     pmf1 = create_discrete_pmf(cont_dist, Val(:single_censored); Δd = Δd, D = D)
     pmf2 = create_discrete_pmf(cont_dist; Δd = Δd, D = D)
 
@@ -80,7 +80,7 @@ plt1 = let
         title = "Discrete PMF with Δd = 1 day",
         label = ["Single censoring (midpoint primary)" "Double Censoring"],
         xlabel = "Days",
-        ylabel = "Probability",
+        ylabel = "Probability"
     )
 end
 savefig(plt1, joinpath(@__DIR__(), "assets/", "discrete_pmf_daily.png"))
@@ -89,7 +89,7 @@ savefig(plt1, joinpath(@__DIR__(), "assets/", "discrete_pmf_daily.png"))
 
 plt2 = let
     Δd = 1 / 24
-    ts = (0.0:Δd:(D-Δd)) |> collect
+    ts = (0.0:Δd:(D - Δd)) |> collect
     pmf1 = create_discrete_pmf(cont_dist, Val(:single_censored); Δd = Δd, D = D)
     pmf2 = create_discrete_pmf(cont_dist; Δd = Δd, D = D)
 
@@ -101,7 +101,7 @@ plt2 = let
         title = "Discrete PMF with Δd = 1 hour",
         label = ["Single censoring (midpoint primary)" "Double Censoring"],
         xlabel = "Days",
-        ylabel = "Probability",
+        ylabel = "Probability"
     )
 end
 savefig(plt2, joinpath(@__DIR__(), "assets/", "discrete_pmf_hourly.png"))
diff --git a/EpiAware/test/predictive_checking/fast_approx_for_r.jl b/EpiAware/test/predictive_checking/fast_approx_for_r.jl
index 2646065ef..a36f89dd0 100644
--- a/EpiAware/test/predictive_checking/fast_approx_for_r.jl
+++ b/EpiAware/test/predictive_checking/fast_approx_for_r.jl
@@ -39,7 +39,6 @@ TestEnv.activate()
 ```
 =#
 
-
 using EpiAware
 using Distributions
 using StatsPlots
@@ -47,8 +46,7 @@ using StatsPlots
 # Create a discrete probability mass function (PMF) for a negative binomial distribution
 # with left truncation at 1.
 
-w =
-    create_discrete_pmf(NegativeBinomial(2, 0.5), D = 20.0) |>
+w = create_discrete_pmf(NegativeBinomial(2, 0.5), D = 20.0) |>
     p -> p[2:end] ./ sum(p[2:end])
 
 ##
@@ -76,5 +74,5 @@ plot(
     title = "Fast approximation for r",
     lab = ["T_2 = 1.0" "T_2 = 3.5" "T_2 = 7.0" "T_2 = 14.0"],
     yticks = [0.0, 1e-15, 1e-10, 1e-5, 1e0] |> x -> (x .+ jitter, string.(x)),
-    xticks = 0:2:10,
+    xticks = 0:2:10
 )
diff --git a/EpiAware/test/predictive_checking/toy_model_log_infs_RW.jl b/EpiAware/test/predictive_checking/toy_model_log_infs_RW.jl
index 41f0993c4..d2a46b632 100644
--- a/EpiAware/test/predictive_checking/toy_model_log_infs_RW.jl
+++ b/EpiAware/test/predictive_checking/toy_model_log_infs_RW.jl
@@ -91,7 +91,7 @@ model_data = EpiData(
     neg_bin_cluster_factor,
     time_horizon,
     D_gen = 10.0,
-    D_delay = 10.0,
+    D_delay = 10.0
 )
 
 #=
@@ -113,7 +113,7 @@ log_infs_model = make_epi_inference_model(
     random_walk,
     delay_observations;
     process_priors = merge(default_rw_priors(), default_delay_obs_priors()),
-    pos_shift = 1e-6,
+    pos_shift = 1e-6
 )
 
 #=
@@ -135,7 +135,7 @@ plot(
     label = "I_t",
     xlabel = "Time",
     ylabel = "Infections",
-    title = "Generated Infections",
+    title = "Generated Infections"
 )
 scatter!(random_epidemic.y_t, lab = "generated cases")
 
@@ -153,7 +153,7 @@ model = make_epi_inference_model(
     random_walk,
     delay_observations;
     process_priors = merge(default_rw_priors(), default_delay_obs_priors()),
-    pos_shift = 1e-6,
+    pos_shift = 1e-6
 )
 
 @time chn = sample(
@@ -162,7 +162,7 @@ model = make_epi_inference_model(
     MCMCThreads(),
     250,
     4;
-    drop_warmup = true,
+    drop_warmup = true
 )
 
 #=
@@ -182,7 +182,7 @@ scatter!(
     xlabel = "Time",
     ylabel = "Cases",
     title = "Posterior Predictive Checking",
-    ylims = (-0.5, maximum(truth_data) * 2.5),
+    ylims = (-0.5, maximum(truth_data) * 2.5)
 )
 
 #=
@@ -200,7 +200,7 @@ scatter!(
     xlabel = "Time",
     ylabel = "Unobserved Infections",
     title = "Posterior Predictive Checking",
-    ylims = (-0.5, maximum(gen.I_t) * 1.5),
+    ylims = (-0.5, maximum(gen.I_t) * 1.5)
 )
 
 #=
diff --git a/EpiAware/test/prior_predictive_checking/ppc-latent-processes.jl b/EpiAware/test/prior_predictive_checking/ppc-latent-processes.jl
index 7f212d3ac..4ec65f2e8 100644
--- a/EpiAware/test/prior_predictive_checking/ppc-latent-processes.jl
+++ b/EpiAware/test/prior_predictive_checking/ppc-latent-processes.jl
@@ -20,13 +20,11 @@ sampled_walks = prior_chn |> chn -> mapreduce(hcat, generated_quantities(model,
     gen[1]
 end
 ## From law of total variance and known mean of HalfNormal distribution
-theoretical_std =
-    [
-        t * latent_process_priors.var_RW_dist.untruncated.σ * sqrt(2) / sqrt(π) for t = 1:n
-    ] .|> sqrt
+theoretical_std = [t * latent_process_priors.var_RW_dist.untruncated.σ * sqrt(2) / sqrt(π)
+                   for t in 1:n] .|> sqrt
 
-plt_ppc_rw =
-    plot(sampled_walks, lab = "", ylabel = "RW", xlabel = "t", c = :grey, alpha = 0.1)
+plt_ppc_rw = plot(
+    sampled_walks, lab = "", ylabel = "RW", xlabel = "t", c = :grey, alpha = 0.1)
 plot!(
     plt_ppc_rw,
     zeros(n),
@@ -34,7 +32,7 @@ plot!(
     c = :red,
     lab = "Theoretical 3 sigma spread",
     ribbon = 3 * theoretical_std,
-    fillalpha = 0.2,
+    fillalpha = 0.2
 )
 
 σ_hist = histogram(
@@ -44,7 +42,7 @@ plot!(
     ylabel = "Density",
     xlabel = "σ²_RW",
     c = :grey,
-    alpha = 0.5,
+    alpha = 0.5
 )
 plot!(
     σ_hist,
@@ -53,7 +51,7 @@ plot!(
     c = :red,
     alpha = 0.5,
     lab = "Prior",
-    bins = 100,
+    bins = 100
 )
 
 plt_rw = plot(
@@ -62,6 +60,6 @@ plt_rw = plot(
     layout = (1, 2),
     size = (800, 400),
     left_margin = 3mm,
-    bottom_margin = 3mm,
+    bottom_margin = 3mm
 )
 savefig(plt_rw, joinpath(@__DIR__(), "assets", "ppc_rw.png"))
diff --git a/EpiAware/test/test_epimodel.jl b/EpiAware/test/test_epimodel.jl
index b2d9493ce..c975ef350 100644
--- a/EpiAware/test/test_epimodel.jl
+++ b/EpiAware/test/test_epimodel.jl
@@ -38,7 +38,7 @@ end
         cluster_coeff,
         time_horizon;
         D_gen = 10.0,
-        D_delay = 10.0,
+        D_delay = 10.0
     )
 
     @test data.cluster_coeff == 0.8
@@ -51,7 +51,6 @@ end
     @test size(data.delay_kernel) == (time_horizon, time_horizon)
 end
 
-
 @testitem "Renewal function: internal generate infs" begin
     using LinearAlgebra
     gen_int = [0.2, 0.3, 0.5]
@@ -65,16 +64,15 @@ end
 
     function generate_infs(recent_incidence, Rt)
         new_incidence = Rt * dot(recent_incidence, epimodel.data.gen_int)
-        [new_incidence; recent_incidence[1:(epimodel.data.len_gen_int-1)]], new_incidence
+        [new_incidence; recent_incidence[1:(epimodel.data.len_gen_int - 1)]], new_incidence
     end
 
     recent_incidence = [10, 20, 30]
     Rt = 1.5
 
     expected_new_incidence = Rt * dot(recent_incidence, [0.2, 0.3, 0.5])
-    expected_output =
-        [expected_new_incidence; recent_incidence[1:2]], expected_new_incidence
-
+    expected_output = [expected_new_incidence; recent_incidence[1:2]],
+    expected_new_incidence
 
     @test generate_infs(recent_incidence, Rt) == expected_output
 end
diff --git a/EpiAware/test/test_latent-processes.jl b/EpiAware/test/test_latent-processes.jl
index 9fb85a236..b8c7f7f00 100644
--- a/EpiAware/test/test_latent-processes.jl
+++ b/EpiAware/test/test_latent-processes.jl
@@ -5,11 +5,10 @@
     model = random_walk(n)
     fixed_model = fix(model, (σ²_RW = 1.0, init_rw_value = 0.0)) #Fixing the standard deviation of the random walk process
     n_samples = 1000
-    samples_day_5 =
-        sample(fixed_model, Prior(), n_samples) |>
-        chn -> mapreduce(vcat, generated_quantities(fixed_model, chn)) do gen
-            gen[1][5] #Extracting day 5 samples
-        end
+    samples_day_5 = sample(fixed_model, Prior(), n_samples) |>
+                    chn -> mapreduce(vcat, generated_quantities(fixed_model, chn)) do gen
+        gen[1][5] #Extracting day 5 samples
+    end
     #Check statistics are within 5 sigma
     #Theoretically, after 5 steps distribution is N(0, var = 5)
     theoretical_std_of_empiral_mean = sqrt(5 / n_samples)
@@ -23,10 +22,8 @@
     @info "var = $(var(samples_day_5)); theoretical_std_of_empiral_var = $(theoretical_std_of_empiral_var)"
     @test (var(samples_day_5) - 5) < 5 * theoretical_std_of_empiral_var &&
           (var(samples_day_5) - 5) > -5 * theoretical_std_of_empiral_var
-
 end
 @testitem "Testing default_rw_priors" begin
-
     @testset "var_RW_dist" begin
         priors = default_rw_priors()
         var_RW = rand(priors.var_RW_dist)
diff --git a/EpiAware/test/test_models.jl b/EpiAware/test/test_models.jl
index 61d4493b0..7e94bccb8 100644
--- a/EpiAware/test/test_models.jl
+++ b/EpiAware/test/test_models.jl
@@ -7,7 +7,6 @@
     process_priors = merge(default_rw_priors(), default_delay_obs_priors())
     pos_shift = 1e-6
 
-
     epimodel = DirectInfections(data)
 
     # Call the function
@@ -17,17 +16,17 @@
         random_walk,
         delay_observations;
         process_priors,
-        pos_shift,
+        pos_shift
     )
 
     # Define expected outputs for a conditional model
     # Underlying log-infections are const value 1 for all time steps and
     # any other unfixed parameters
 
-    fixed_test_mdl =
-        fix(test_mdl, (init = log(1.0), σ²_RW = 0.0, neg_bin_cluster_factor = 0.05))
+    fixed_test_mdl = fix(
+        test_mdl, (init = log(1.0), σ²_RW = 0.0, neg_bin_cluster_factor = 0.05))
     X = rand(fixed_test_mdl)
-    expected_I_t = [1.0 for _ = 1:epimodel.data.time_horizon]
+    expected_I_t = [1.0 for _ in 1:(epimodel.data.time_horizon)]
     gen = generated_quantities(fixed_test_mdl, rand(fixed_test_mdl))
 
     # Perform tests
@@ -51,17 +50,17 @@ end
         random_walk,
         delay_observations;
         process_priors,
-        pos_shift,
+        pos_shift
     )
 
     # Define expected outputs for a conditional model
     # Underlying log-infections are const value 1 for all time steps and
     # any other unfixed parameters
 
-    fixed_test_mdl =
-        fix(test_mdl, (init = log(1.0), σ²_RW = 0.0, neg_bin_cluster_factor = 0.05))
+    fixed_test_mdl = fix(
+        test_mdl, (init = log(1.0), σ²_RW = 0.0, neg_bin_cluster_factor = 0.05))
     X = rand(fixed_test_mdl)
-    expected_I_t = [1.0 for _ = 1:epimodel.data.time_horizon]
+    expected_I_t = [1.0 for _ in 1:(epimodel.data.time_horizon)]
     gen = generated_quantities(fixed_test_mdl, rand(fixed_test_mdl))
 
     # # Perform tests
@@ -85,17 +84,17 @@ end
         random_walk,
         delay_observations;
         process_priors,
-        pos_shift,
+        pos_shift
     )
 
     # Define expected outputs for a conditional model
     # Underlying log-infections are const value 1 for all time steps and
     # any other unfixed parameters
 
-    fixed_test_mdl =
-        fix(test_mdl, (init = log(1.0), σ²_RW = 0.0, neg_bin_cluster_factor = 0.05))
+    fixed_test_mdl = fix(
+        test_mdl, (init = log(1.0), σ²_RW = 0.0, neg_bin_cluster_factor = 0.05))
     X = rand(fixed_test_mdl)
-    expected_I_t = [1.0 for _ = 1:epimodel.data.time_horizon]
+    expected_I_t = [1.0 for _ in 1:(epimodel.data.time_horizon)]
     gen = generated_quantities(fixed_test_mdl, rand(fixed_test_mdl))
 
     # # Perform tests
diff --git a/EpiAware/test/test_observation-processes.jl b/EpiAware/test/test_observation-processes.jl
index 238bfd5fa..b30266f65 100644
--- a/EpiAware/test/test_observation-processes.jl
+++ b/EpiAware/test/test_observation-processes.jl
@@ -15,17 +15,16 @@
         missing,
         I_t,
         epimodel;
-        observation_process_priors = observation_process_priors,
+        observation_process_priors = observation_process_priors
     )
     fix_mdl = fix(mdl, neg_bin_cluster_factor = 0.00001) # Effectively Poisson sampling
 
     n_samples = 1000
-    mean_first_obs =
-        sample(fix_mdl, Prior(), n_samples) |>
-        chn -> generated_quantities(fix_mdl, chn) .|> (gen -> gen[1][1]) |> mean
+    mean_first_obs = sample(fix_mdl, Prior(), n_samples) |>
+                     chn -> generated_quantities(fix_mdl, chn) .|> (gen -> gen[1][1]) |>
+                            mean
 
     theoretical_std_of_empiral_mean = sqrt(I_t[1]) / sqrt(n_samples)
     @test mean(mean_first_obs) - I_t[1] < 5 * theoretical_std_of_empiral_mean &&
           mean(mean_first_obs) - I_t[1] > -5 * theoretical_std_of_empiral_mean
-
 end
diff --git a/EpiAware/test/test_utilities.jl b/EpiAware/test/test_utilities.jl
index f1fc6173d..47e869ca6 100644
--- a/EpiAware/test/test_utilities.jl
+++ b/EpiAware/test/test_utilities.jl
@@ -50,35 +50,32 @@ end
     # interval censoring with left hand approx.
     @testset "Test case 4" begin
         dist = Exponential(1.0)
-        expected_pmf = [(exp(-(t - 1)) - exp(-t)) / (1 - exp(-5)) for t = 1:5]
+        expected_pmf = [(exp(-(t - 1)) - exp(-t)) / (1 - exp(-5)) for t in 1:5]
         pmf = create_discrete_pmf(
             dist,
             Val(:single_censored);
             primary_approximation_point = 0.0,
             Δd = 1.0,
-            D = 5.0,
+            D = 5.0
         )
-        @test pmf ≈ expected_pmf atol = 1e-15
+        @test pmf≈expected_pmf atol=1e-15
     end
 
     # Test case 5: Testing output against expected PMF basic version - double
     # interval censoring
     @testset "Test case 5" begin
         dist = Exponential(1.0)
-        expected_pmf_uncond = [
-            exp(-1)
-            [(1 - exp(-1)) * (exp(1) - 1) * exp(-s) for s = 1:9]
-        ]
+        expected_pmf_uncond = [exp(-1)
+                               [(1 - exp(-1)) * (exp(1) - 1) * exp(-s) for s in 1:9]]
         expected_pmf = expected_pmf_uncond ./ sum(expected_pmf_uncond)
         pmf = create_discrete_pmf(dist; Δd = 1.0, D = 10.0)
-        @test expected_pmf ≈ pmf atol = 1e-15
+        @test expected_pmf≈pmf atol=1e-15
     end
 
     @testset "Test case 6" begin
         dist = Exponential(1.0)
         @test_throws AssertionError create_discrete_pmf(dist, Δd = 1.0, D = 3.5)
     end
-
 end
 
 @testitem "Testing r_to_R function" begin
@@ -88,7 +85,7 @@ end
         w = ones(5) |> x -> x ./ sum(x)
         expected_ratio = 1
         ratio = EpiAware.r_to_R(r, w)
-        @test ratio ≈ expected_ratio atol = 1e-15
+        @test ratio≈expected_ratio atol=1e-15
     end
 
     #Test MethodError when w is not a vector
@@ -97,7 +94,6 @@ end
         w = 1
         @test_throws MethodError EpiAware.r_to_R(r, w)
     end
-
 end
 
 @testitem "Testing generate_observation_kernel function" begin
@@ -106,18 +102,15 @@ end
         delay_int = [0.2, 0.5, 0.3]
         time_horizon = 5
         expected_K = SparseMatrixCSC(
-            [
-                0.2 0 0 0 0
-                0.5 0.2 0 0 0
-                0.3 0.5 0.2 0 0
-                0 0.3 0.5 0.2 0
-                0 0 0.3 0.5 0.2
-            ],
+            [0.2 0 0 0 0
+             0.5 0.2 0 0 0
+             0.3 0.5 0.2 0 0
+             0 0.3 0.5 0.2 0
+             0 0 0.3 0.5 0.2],
         )
         K = EpiAware.generate_observation_kernel(delay_int, time_horizon)
         @test K == expected_K
     end
-
 end
 
 @testitem "Testing neg_MGF function" begin
@@ -127,19 +120,18 @@ end
         w = [0.2, 0.3, 0.5]
         expected_result = 0.2 * exp(-0.5 * 1) + 0.3 * exp(-0.5 * 2) + 0.5 * exp(-0.5 * 3)
         result = EpiAware.neg_MGF(r, w)
-        @test result ≈ expected_result atol = 1e-15
+        @test result≈expected_result atol=1e-15
     end
 
     # Test case 2: Testing with zero r and non-empty weight vector
     @testset "Test case 2" begin
         r = 0
         w = [0.1, 0.2, 0.3, 0.4]
-        expected_result =
-            0.1 * exp(-0 * 1) + 0.2 * exp(-0 * 2) + 0.3 * exp(-0 * 3) + 0.4 * exp(-0 * 4)
+        expected_result = 0.1 * exp(-0 * 1) + 0.2 * exp(-0 * 2) + 0.3 * exp(-0 * 3) +
+                          0.4 * exp(-0 * 4)
         result = EpiAware.neg_MGF(r, w)
-        @test result ≈ expected_result atol = 1e-15
+        @test result≈expected_result atol=1e-15
     end
-
 end
 
 @testitem "Testing dneg_MGF_dr function" begin
@@ -148,10 +140,10 @@ end
     @testset "Test case 1" begin
         r = 0.5
         w = [0.2, 0.3, 0.5]
-        expected_result =
-            -(0.2 * 1 * exp(-0.5 * 1) + 0.3 * 2 * exp(-0.5 * 2) + 0.5 * 3 * exp(-0.5 * 3))
+        expected_result = -(0.2 * 1 * exp(-0.5 * 1) + 0.3 * 2 * exp(-0.5 * 2) +
+                            0.5 * 3 * exp(-0.5 * 3))
         result = EpiAware.dneg_MGF_dr(r, w)
-        @test result ≈ expected_result atol = 1e-15
+        @test result≈expected_result atol=1e-15
     end
 
     # Test case 2: Testing with zero r and non-empty weight vector
@@ -165,9 +157,8 @@ end
             0.4 * 4 * exp(-0 * 4)
         )
         result = EpiAware.dneg_MGF_dr(r, w)
-        @test result ≈ expected_result atol = 1e-15
+        @test result≈expected_result atol=1e-15
     end
-
 end
 @testitem "Testing spread_draws function" begin
     using DataFramesMeta, Turing