diff --git a/tests/testthat/test-twdtw_classify.R b/tests/testthat/test-twdtw_classify.R
index 1e591a3..7ba8852 100644
--- a/tests/testthat/test-twdtw_classify.R
+++ b/tests/testthat/test-twdtw_classify.R
@@ -1,28 +1,28 @@
+cat("Loading mgcv...\n")
 library(mgcv)
 
-# Read training samples
-samples <- st_read(system.file("mato_grosso_brazil/samples.gpkg", package = "dtwSat"), quiet = TRUE)
+cat("Reading training samples...\n")
+samples <- st_read(system.file("mato_grosso_brazil/samples.gpkg", package = "dtwSat"))
 
-# Satellite image time sereis files
+cat("Geting satellite image time series files...\n")
 tif_files <- dir(system.file("mato_grosso_brazil", package = "dtwSat"), pattern = "\\.tif$", full.names = TRUE)
 
-# The acquisition date is in the file name are not the true acquisition date
-# of each pixel. MOD13Q1 is a 16-day composite product, so the acquisition date
-# is the first day of the 16-day period
+cat("Geting acquisition dates from file names...\n")
 acquisition_date <- as.Date(regmatches(tif_files, regexpr("[0-9]{8}", tif_files)), format = "%Y%m%d")
 
-# Read the data as a stars object setting the time/date for each observation
-# using along. This will prodcue a 4D array (data-cube) which will then be converted
-# to a 3D array by spliting the 'band' dimension
+cat("Reading data as stars object...\n")
 dc <- read_stars(tif_files,
                  proxy = FALSE,
                  along = list(time = acquisition_date),
                  RasterIO = list(bands = 1:6))
 
+cat("Setting band names...\n")
 dc <- st_set_dimensions(dc, 3, c("EVI", "NDVI", "RED", "BLUE", "NIR", "MIR"))
+
+cat("Converting 4D into a 3D datacube...\n")
 dc <- split(dc, c("band"))
 
-# Create a knn1-twdtw model
+cat("Creating twdtw_knn1 model using GAM...\n")
 system.time(
   m <- twdtw_knn1(x = dc,
                   y = samples,
@@ -32,22 +32,24 @@ system.time(
                   time_weight = c(steepness = 0.1, midpoint = 50))
 )
 
+cat("Printing model...\n")
 print(m)
 
-# Visualize model patterns
+cat("Visualizing model patterns...\n")
 plot(m)
 
-# Classify satellite image time series
+cat("Classifying satellite image time series...\n")
 system.time(lu <- predict(dc, model = m))
 
+cat("Printing land use classification...\n")
 print(lu)
 
-# # Visualise land use classification
+cat("Visualizing land use classification...\n")
 ggplot() +
   geom_stars(data = lu) +
   theme_minimal()
 
-# ### OTHER TESTS
+cat("Testing model with resampling frequency of 60 and GAM...\n")
 m <- twdtw_knn1(x = dc,
                 y = samples,
                 cycle_length = 'year',
@@ -56,8 +58,10 @@ m <- twdtw_knn1(x = dc,
                 smooth_fun = function(x, y) gam(y ~ s(x), data = data.frame(x = x, y = y)),
                 resampling_freq = 60)
 
+cat("Visualizing resampled patterns...\n")
 plot(m)
 
+cat("Classifying satellite image time series with resampled patterns...\n")
 system.time(
   lu <- predict(dc,
                 model = m,
@@ -66,14 +70,15 @@ system.time(
                 time_weight = c(steepness = 0.1, midpoint = 50))
 )
 
+cat("Printing land use classification of resampled patterns...\n")
 print(lu)
 
-# Visualise land use classification
+cat("Visualizing land use classification of resampled patterns...\n")
 ggplot() +
   geom_stars(data = lu) +
   theme_minimal()
 
-# Test model without samples reduction
+cat("Testing model without sample reduction...\n")
 m <- twdtw_knn1(x = dc,
                 y = samples,
                 cycle_length = 'year',
@@ -81,13 +86,16 @@ m <- twdtw_knn1(x = dc,
                 time_weight = c(steepness = 0.1, midpoint = 50),
                 smooth_fun = NULL)
 
+cat("Printing model without sample reduction...\n")
 print(m)
 
+cat("Visualizing model patterns without sample reduction...\n")
 plot(m)
 
+cat("Visualizing model patterns without sample reduction for EVI and NDVI...\n")
 plot(m, bands = c('EVI', 'NDVI'))
 
-# Test custom smooth function -- mean value for each time in dc
+cat("Testing custom smooth function using the average for each observation date...\n")
 m <- twdtw_knn1(x = dc,
                 y = samples,
                 cycle_length = 'year',
@@ -95,8 +103,24 @@ m <- twdtw_knn1(x = dc,
                 time_weight = c(steepness = 0.1, midpoint = 50),
                 smooth_fun = function(x, y) lm(y ~ factor(x), data = data.frame(x = x, y = y)))
 
+cat("Printing model with custom smooth...\n")
 print(m)
 
+cat("Visualizing model patterns with custom smooth...\n")
 plot(m)
 
+cat("Visualizing model patterns with custom smooth for EVI and NDVI...\n")
 plot(m, bands = c('EVI', 'NDVI'))
+
+cat("Classifying satellite image time series with custom smooth...\n")
+system.time(lu <- predict(dc, model = m))
+
+cat("Printing land use classification with custom smooth...\n")
+print(lu)
+
+cat("Visualizing land use classification with custom smooth...\n")
+ggplot() +
+  geom_stars(data = lu) +
+  theme_minimal()
+
+cat("Script completed.\n")