From 3017401bfdeb2789d77d0a95b7b0072c17cd89f0 Mon Sep 17 00:00:00 2001
From: Sara Moreno <sara.moreno@bsc.es>
Date: Fri, 14 Apr 2023 11:38:43 +0200
Subject: [PATCH] Fixed ncores

---
 R/Analogs.R       |  52 ++++++-------
 R/Intbc.R         |  42 ++++++-----
 R/Interpolation.R |  82 +++++++++++++-------
 R/Intlr.R         | 186 +++++++++++++++++++++++++---------------------
 R/LogisticReg.R   |  73 +++++++++---------
 5 files changed, 242 insertions(+), 193 deletions(-)

diff --git a/R/Analogs.R b/R/Analogs.R
index a69a66d..f0ebd61 100644
--- a/R/Analogs.R
+++ b/R/Analogs.R
@@ -66,8 +66,8 @@
 #'@param loocv_window a logical vector only to be used if 'obs' does not have the dimension
 #''window'. It indicates whether to apply leave-one-out cross-validation in the creation 
 #'of the window. It is recommended to be set to TRUE. Default to TRUE.
-#'@param ncores an integer indicating the number of cores to use in parallel computation.
-#'
+#'@param ncores an integer indicating the number of cores to use in parallel computation. 
+#'The default value is NULL.
 #'@return An 's2dv' object. The element 'data' contains the dowscaled field, 'lat' the
 #'downscaled latitudes, and 'lon' the downscaled longitudes. If fun_analog is set to NULL
 #'(default), the output array in 'data' also contains the dimension 'analog' with the best 
@@ -87,7 +87,7 @@
 #'@export
 CST_Analogs <- function(exp, obs, grid_exp, obs2 = NULL, nanalogs = 3, fun_analog = NULL, lat_dim = "lat", 
                         lon_dim = "lon", sdate_dim = "sdate", time_dim = "time", member_dim = "member",
-                        region = NULL, return_indices = FALSE, loocv_window = TRUE, ncores = 1) {
+                        region = NULL, return_indices = FALSE, loocv_window = TRUE, ncores = NULL) {
 
   # input exp and obs must be s2dv_cube objects
   if (!inherits(exp,'s2dv_cube')) {
@@ -204,7 +204,8 @@ CST_Analogs <- function(exp, obs, grid_exp, obs2 = NULL, nanalogs = 3, fun_analo
 #'@param loocv_window a logical vector only to be used if 'obs' does not have the dimension
 #''window'. It indicates whether to apply leave-one-out cross-validation in the creation 
 #'of the window. It is recommended to be set to TRUE. Default to TRUE.
-#'@param ncores an integer indicating the number of cores to use in parallel computation.
+#'@param ncores an integer indicating the number of cores to use in parallel computation. 
+#'The default value is NULL.
 #'@import multiApply
 #'@import CSTools
 #'@importFrom s2dv InsertDim CDORemap
@@ -232,7 +233,7 @@ Analogs <- function(exp, obs, exp_lats, exp_lons, obs_lats, obs_lons, grid_exp,
                     obs2_lats = NULL, obs2_lons = NULL, nanalogs = 3, fun_analog = NULL, 
                     lat_dim = "lat", lon_dim = "lon", sdate_dim = "sdate", time_dim = "time", 
                     member_dim = "member", region = NULL, return_indices = FALSE, 
-                    loocv_window = TRUE, ncores = 1) {
+                    loocv_window = TRUE, ncores = NULL) {
   #-----------------------------------
   # Checkings
   #-----------------------------------
@@ -336,7 +337,14 @@ Analogs <- function(exp, obs, exp_lats, exp_lons, obs_lats, obs_lons, grid_exp,
       stop("Missing time dimension in 'obs2', or does not match the parameter 'time_dim'")
     }
   }
-
+  ## ncores
+  if (!is.null(ncores)) {
+    if (!is.numeric(ncores) | any(ncores %% 1 != 0) | any(ncores < 0) |
+        length(ncores) > 1) {
+      stop("Parameter 'ncores' must be a positive integer.")
+    }
+  }
+  
   # Select a function to apply to the analogs selected for a given observation
   if (!is.null(fun_analog)) {
     stopifnot(fun_analog %in% c("mean", "wmean", "max", "min", "median"))
@@ -370,22 +378,22 @@ Analogs <- function(exp, obs, exp_lats, exp_lons, obs_lats, obs_lons, grid_exp,
 
   obs_interpolated <- Interpolation(exp = obs_train, lats = obs_train_lats, lons = obs_train_lons, 
                                     target_grid = grid_exp, lat_dim = lat_dim, lon_dim = lon_dim, 
-                                    method_remap = "conservative", region = region)
+                                    method_remap = "conservative", region = region, ncores = ncores)
   # If after interpolating 'obs' data the coordinates do not match, the exp data is interpolated to
   # the same grid to force the matching
   if (!.check_coords(lat1 = obs_interpolated$lat, lat2 = exp_lats, lon1 = obs_interpolated$lon, lon2 = exp_lons)) {
     exp_interpolated <- Interpolation(exp = exp, lats = exp_lats, lons = exp_lons, target_grid = grid_exp, 
                                       lat_dim = lat_dim, lon_dim = lon_dim, method_remap = "conservative", 
-                                      region = region)$data
+                                      region = region, ncores = ncores)$data
   } else {
     exp_interpolated <- exp
   }
 
   # Create window if user does not have it in the training observations
   if ( !("window" %in% names(dim(obs_interpolated$data))) ) {
-    obs_train_interpolated <- generate_window(obj = obs_interpolated$data, sdate_dim = sdate_dim, 
-                                              time_dim = time_dim, loocv = loocv_window)
-    obs_hres <- generate_window(obj = obs, sdate_dim = sdate_dim, time_dim = time_dim, loocv = loocv_window)
+    obs_train_interpolated <- .generate_window(obj = obs_interpolated$data, sdate_dim = sdate_dim, 
+                                              time_dim = time_dim, loocv = loocv_window, ncores = ncores)
+    obs_hres <- .generate_window(obj = obs, sdate_dim = sdate_dim, time_dim = time_dim, loocv = loocv_window, ncores = ncores)
   }
 
   #-----------------------------------
@@ -411,7 +419,7 @@ Analogs <- function(exp, obs, exp_lats, exp_lons, obs_lats, obs_lons, grid_exp,
       ntimes <- dim(res.ind)[names(dim(res.ind)) == time_dim]
       res.ind <- Apply(res.ind, target_dims = c("ind", sdate_dim), function(x)   
         sapply(1:nsdates, function(s) seq(ntimes * nsdates)[ - (ntimes * (s - 1) + 1:ntimes)][x[, s]]),
-        output_dims = c("ind", sdate_dim))$output1
+        output_dims = c("ind", sdate_dim), ncores = ncores)$output1
     }
 
     # restore ensemble dimension in observations if it existed originally
@@ -484,7 +492,7 @@ Analogs <- function(exp, obs, exp_lats, exp_lons, obs_lats, obs_lons, grid_exp,
 # Add the dimension window to an array that contains, at least, the start date and time
 # dimensions
 # object has at least dimensions sdate and time
-generate_window <- function(obj, sdate_dim, time_dim, loocv, size = NULL) {
+.generate_window <- function(obj, sdate_dim, time_dim, loocv, size = NULL, ncores = NULL) {
 
   rsdates <- 1:dim(obj)[names(dim(obj)) == sdate_dim]
   ntimes <- dim(obj)[names(dim(obj)) == time_dim] 
@@ -497,11 +505,11 @@ generate_window <- function(obj, sdate_dim, time_dim, loocv, size = NULL) {
     if (loocv) {
         obj_window <- Apply(obj, target_dims = c(time_dim, sdate_dim),
           fun = function(x) sapply(rsdates, function(s) as.vector(x[ rtimes, -s])), 
-          output_dims = c('window', sdate_dim))$output1
+          output_dims = c('window', sdate_dim), ncores = ncores)$output1
     # Generate window without cross-validation        
     } else {
         obj_window <- Apply(obj, target_dims = c(time_dim, sdate_dim),
-          fun = as.vector, output_dims = 'window')$output1
+          fun = as.vector, output_dims = 'window', ncores = ncores)$output1
     }
   }
   # Generate window of the size specified by the user. Only applied with CV
@@ -513,28 +521,20 @@ generate_window <- function(obj, sdate_dim, time_dim, loocv, size = NULL) {
    
     # Concatenate data from previous, target and posterior months
     obj_new <- Apply(obj, target_dims = list(c("time", "smonth")),
-      fun = as.vector, output_dims = "time")$output1
+      fun = as.vector, output_dims = "time", ncores = ncores )$output1
 
     if (loocv) {
       obj_window <- Apply(list(obj_new, rtimes, rsdates), target_dims = list(c(time_dim, sdate_dim), NULL, NULL),
         fun = function(x, t, s) as.vector(x[(ntimes + t - size):(ntimes + t + size), -s]), 
-        output_dims = 'window')$output1  
+        output_dims = 'window', ncores = ncores)$output1  
       names(dim(obj_window))[(length(names(dim(obj_window))) - 1):length(names(dim(obj_window)))] <- c(time_dim, sdate_dim) 
     } else {
       obj_window <- Apply(obj_new, target_dims = c(time_dim, sdate_dim),
         fun = function(x) sapply(rtimes, function(t) as.vector(x[(ntimes + t - size):(ntimes + t + size), ])),
-        output_dims = c('window', time_dim))$output1
+        output_dims = c('window', time_dim), ncores = ncores)$output1
 
     }
   } 
 
   return(obj_window)
 }
-
-
-
-
-
-
-
-
diff --git a/R/Intbc.R b/R/Intbc.R
index 27de373..1cb558d 100644
--- a/R/Intbc.R
+++ b/R/Intbc.R
@@ -49,7 +49,8 @@
 #'to the left border, while lonmax refers to the right border. latmin indicates the lower
 #'border, whereas latmax indicates the upper border. If set to NULL (default), the function
 #'takes the first and last elements of the latitudes and longitudes.
-#'@param ncores an integer indicating the number of cores to use in parallel computation.
+#'@param ncores an integer indicating the number of cores to use in parallel computation. 
+#'The default value is NULL.
 #'@return An 's2dv' object. The element 'data' contains the dowscaled field, 'lat' the
 #'downscaled latitudes, and 'lon' the downscaled longitudes.
 #'@examples
@@ -68,7 +69,7 @@
 
 CST_Intbc <- function(exp, obs, target_grid, bc_method, int_method = NULL, points = NULL, 
                       method_point_interp = NULL, lat_dim = "lat", lon_dim = "lon", 
-                      sdate_dim = "sdate", member_dim = "member", region = NULL, ncores = 1,...)
+                      sdate_dim = "sdate", member_dim = "member", region = NULL, ncores = NULL, ...)
 {
   if (!inherits(exp,'s2dv_cube')) {
     stop("Parameter 'exp' must be of the class 's2dv_cube'")
@@ -83,7 +84,7 @@ CST_Intbc <- function(exp, obs, target_grid, bc_method, int_method = NULL, point
                int_method = int_method, bc_method = bc_method, points = points, 
                source_file = exp$attrs$source_files[1], method_point_interp = method_point_interp, 
                lat_dim = lat_dim, lon_dim = lon_dim, sdate_dim = sdate_dim, member_dim = member_dim, 
-               region = region, ncores = ncores,...)
+               region = region, ncores = ncores, ...)
 
   # Modify data, lat and lon in the origina s2dv_cube, adding the downscaled data
   exp$data <- res$data
@@ -162,8 +163,8 @@ CST_Intbc <- function(exp, obs, target_grid, bc_method, int_method = NULL, point
 #'to the left border, while lonmax refers to the right border. latmin indicates the lower
 #'border, whereas latmax indicates the upper border. If set to NULL (default), the function
 #'takes the first and last elements of the latitudes and longitudes.
-#'@param ncores an integer indicating the number of cores to use in parallel computation.
-#'
+#'@param ncores an integer indicating the number of cores to use in parallel computation. 
+#'The default value is NULL.
 #'@import CSTools 
 #'
 #'@seealso \code{\link[CSTools]{BiasCorrection}}
@@ -186,7 +187,7 @@ CST_Intbc <- function(exp, obs, target_grid, bc_method, int_method = NULL, point
 #'@export
 Intbc <- function(exp, obs, exp_lats, exp_lons, obs_lats, obs_lons, target_grid, bc_method, int_method = NULL,
                   points = NULL, method_point_interp = NULL, lat_dim = "lat", lon_dim = "lon", sdate_dim = "sdate", 
-                  time_dim = "time", member_dim = "member", source_file = NULL, region = NULL, ncores = 1, ...) {
+                  time_dim = "time", member_dim = "member", source_file = NULL, region = NULL, ncores = NULL, ...) {
 
   if (!inherits(bc_method, 'character')) {
     stop("Parameter 'bc_method' must be of the class 'character'")
@@ -250,11 +251,18 @@ Intbc <- function(exp, obs, exp_lats, exp_lons, obs_lats, obs_lons, target_grid,
             "'obs_lats' and 'obs_lons'.")
     region <- c(obs_lons[1], obs_lons[length(obs_lons)], obs_lats[1], obs_lats[length(obs_lats)])
   }
-
+  ## ncores
+  if (!is.null(ncores)) {
+    if (!is.numeric(ncores) | any(ncores %% 1 != 0) | any(ncores < 0) |
+        length(ncores) > 1) {
+      stop("Parameter 'ncores' must be a positive integer.")
+    }
+  }
+  
   exp_interpolated <- Interpolation(exp = exp, lats = exp_lats, lons = exp_lons, target_grid = target_grid, 
                                     method_remap = int_method, points = points, source_file = source_file,
                                     lat_dim = lat_dim, lon_dim = lon_dim, method_point_interp = method_point_interp, 
-                                    region = region)
+                                    region = region, ncores = ncores)
  
   # If after interpolating 'exp' data the coordinates do not match, the obs data is interpolated to
   # the same grid to force the matching
@@ -263,7 +271,7 @@ Intbc <- function(exp, obs, exp_lats, exp_lons, obs_lats, obs_lons, target_grid,
     obs_interpolated <- Interpolation(exp = obs, lats = obs_lats, lons = obs_lons, target_grid = target_grid, 
                                       method_remap = int_method, points = points, source_file = source_file,
                                       lat_dim = lat_dim, lon_dim = lon_dim, 
-                                      method_point_interp = method_point_interp, region = region)
+                                      method_point_interp = method_point_interp, region = region, ncores = ncores)
     obs_ref <- obs_interpolated$data
   } else {
     obs_ref <- obs
@@ -289,32 +297,32 @@ Intbc <- function(exp, obs, exp_lats, exp_lons, obs_lats, obs_lons, target_grid,
     }
 
     res <- BiasCorrection(exp = exp_interpolated$data, obs = obs_ref, memb_dim = member_dim,
-                          sdate_dim = sdate_dim, ...) 
+                          sdate_dim = sdate_dim, ncores = ncores, ...) 
   } 
   else if (bc_method == 'cal' | bc_method == 'calibration') {
     if (dim(exp_interpolated$data)[member_dim] == 1) {
       stop('Calibration must not be used with only one ensemble member.')
     }
     res <- Calibration(exp = exp_interpolated$data, obs = obs_ref, memb_dim = member_dim, 
-                       sdate_dim = sdate_dim, ...)
+                       sdate_dim = sdate_dim, ncores = ncores, ...)
   }
   else if (bc_method == 'qm' | bc_method == 'quantile_mapping') {
    
     res <- QuantileMapping(exp = exp_interpolated$data, obs = obs_ref, na.rm = TRUE, 
-                           memb_dim = member_dim, sdate_dim = sdate_dim, ...)
+                           memb_dim = member_dim, sdate_dim = sdate_dim, ncores = ncores, ...)
   }
   else if (bc_method == 'dbc' | bc_method == 'dynamical_bias') {
     # the temporal dimension must be only one dimension called "time"
     if (all(c(time_dim, sdate_dim) %in% names(dim(exp_interpolated$data)))) {
       exp_interpolated$data <- Apply(exp_interpolated$data, target_dims = c(time_dim, sdate_dim), 
-        fun = as.vector, output_dims = "time")$output1
+        fun = as.vector, output_dims = "time", ncores = ncores)$output1
     } 
     if (all(c(time_dim, sdate_dim) %in% names(dim(obs_ref)))) {
       obs_ref <- Apply(obs_ref, target_dims = c(time_dim, sdate_dim), fun = as.vector, 
-        output_dims = "time")$output1
+        output_dims = "time", ncores = ncores)$output1
     }
     # REMEMBER to add na.rm = T in colMeans in .proxiesattractor 
-    res <- DynBiasCorrection(exp = exp_interpolated$data, obs = obs_ref, ...)
+    res <- DynBiasCorrection(exp = exp_interpolated$data, obs = obs_ref, ncores = ncores, ...)
   } 
  
   # Return a list of three elements
@@ -322,7 +330,3 @@ Intbc <- function(exp, obs, exp_lats, exp_lons, obs_lats, obs_lons, target_grid,
 
   return(res)
 }
-
-
-
-
diff --git a/R/Interpolation.R b/R/Interpolation.R
index 1599bf3..0a53afa 100644
--- a/R/Interpolation.R
+++ b/R/Interpolation.R
@@ -37,7 +37,8 @@
 #'@param method_point_interp a character vector indicating the interpolation method to
 #'interpolate model gridded data into the point locations. Accepted methods are "nearest",
 #'"bilinear", "9point", "invdist4nn", "NE", "NW", "SE", "SW".
-#'
+#'@param ncores an integer indicating the number of cores to use in parallel computation. 
+#'The default value is NULL.
 #'@seealso \code{\link[s2dverification]{CDORemap}}
 #'
 #'@return An s2dv object containing the dowscaled field.
@@ -52,7 +53,7 @@
 #'@export
 CST_Interpolation <- function(exp, points = NULL, method_remap = NULL, target_grid = NULL, 
                               lat_dim = "lat", lon_dim = "lon", region = NULL, 
-                              method_point_interp = NULL)
+                              method_point_interp = NULL, ncores = NULL)
 {
   if (!inherits(exp,'s2dv_cube')) {
     stop("Parameter 'exp' must be of the class 's2dv_cube'")
@@ -70,7 +71,7 @@ CST_Interpolation <- function(exp, points = NULL, method_remap = NULL, target_gr
   res <- Interpolation(exp = exp$data, lats = exp$coords[[lat_dim]], lons = exp$coords[[lon_dim]], 
                        source_file = exp$attrs$source_files[1], points = points,
                        method_remap = method_remap, target_grid = target_grid, lat_dim = lat_dim,
-                       lon_dim = lon_dim, region = region, method_point_interp = method_point_interp)
+                       lon_dim = lon_dim, region = region, method_point_interp = method_point_interp, ncores = ncores)
 
   # Modify data, lat and lon in the origina s2dv_cube, adding the downscaled data
   exp$data <- res$data
@@ -129,6 +130,8 @@ CST_Interpolation <- function(exp, points = NULL, method_remap = NULL, target_gr
 #'interpolate model gridded data into the point locations. Accepted methods are "nearest",
 #'"bilinear", "9point", "invdist4nn", "NE", "NW", "SE", "SW". Only needed if the downscaling 
 #'is to a point location.
+#'@param ncores an integer indicating the number of cores to use in parallel computation. 
+#'The default value is NULL.
 #'@import multiApply
 #'@import plyr 
 #'@importFrom s2dv CDORemap
@@ -147,7 +150,7 @@ CST_Interpolation <- function(exp, points = NULL, method_remap = NULL, target_gr
 #'@export
 Interpolation <- function(exp, lats, lons, points = NULL, source_file = NULL, method_remap = NULL,  
                           target_grid = NULL, lat_dim = "lat", lon_dim = "lon", region = NULL, 
-                          method_point_interp = NULL) 
+                          method_point_interp = NULL, ncores = NULL) 
 {
   if (!is.null(method_remap)) {
     if (!inherits(method_remap, 'character')) {
@@ -232,7 +235,14 @@ Interpolation <- function(exp, lats, lons, points = NULL, source_file = NULL, me
            "be passed to CDO, and it must be a grid recognised by CDO or a NetCDF file.")
     }
   }
-
+  ## ncores
+  if (!is.null(ncores)) {
+    if (!is.numeric(ncores) | any(ncores %% 1 != 0) | any(ncores < 0) |
+        length(ncores) > 1) {
+      stop("Parameter 'ncores' must be a positive integer.")
+    }
+  }
+  
   #----------------------------------
   # Limits of the region defined by the model data
   #----------------------------------
@@ -256,12 +266,27 @@ Interpolation <- function(exp, lats, lons, points = NULL, source_file = NULL, me
   # Map regrid with CDO
   #----------------------------------
   if (is.null(points)) {
-    res <- CDORemap(data_array = exp, 
-                    lats = lats, 
-                    lons = lons, 
-                    grid = target_grid, 
-                    method = method_remap, 
-                    crop = region)
+    
+    .KnownLonNames <- s2dv:::.KnownLonNames
+    .KnownLatNames <- s2dv:::.KnownLatNames
+    .warning <- s2dv:::.warning
+    
+    res <- multiApply::Apply(data = list(data_array = exp), 
+                             target_dims = c(lon_dim,lat_dim),
+                             fun = CDORemap, 
+                             lats = lats, 
+                             lons = lons,
+                             grid = target_grid, 
+                             method = method_remap,
+                             crop = region, 
+                             ncores = ncores)
+
+    # res <- CDORemap(data_array = exp,
+    #                 lats = lats,
+    #                 lons = lons,
+    #                 grid = target_grid,
+    #                 method = method_remap,
+    #                 crop = region)
 
     # Return a list
     res <- list(data = res$data_array, obs = NULL, lon = res$lons, lat = res$lats)
@@ -271,16 +296,16 @@ Interpolation <- function(exp, lats, lons, points = NULL, source_file = NULL, me
   #----------------------------------
   } else {  
     # First create interpolation weights, depending on the chosen method
-    weights <- create_interp_weights(ncfile = source_file, locids = 1:unique(lengths(points)),
+    weights <- .create_interp_weights(ncfile = source_file, locids = 1:unique(lengths(points)),
                                      lats = points[[lat_dim]], lons = points[[lon_dim]], 
                                      method = method_point_interp, region = list(lat_min = region[3], 
                                      lat_max = region[4], lon_min = region[1], lon_max = region[2])) 
 
     # Select coarse-scale data to be interpolated
-    model_data_gridpoints <- get_model_data(weights.df = weights, mdata = exp)
+    model_data_gridpoints <- .get_model_data(weights.df = weights, mdata = exp, ncores = ncores)
 
     # Interpolate model data to point locations
-    res <- interpolate_data(model_data_gridpoints, weights)
+    res <- .interpolate_data(model_data_gridpoints, weights, ncores = ncores)
   
     # Return a list
     res <- list(data = res, obs = NULL, lon = points[[lon_dim]], lat = points[[lat_dim]])
@@ -293,7 +318,7 @@ Interpolation <- function(exp, lats, lons, points = NULL, source_file = NULL, me
 # Compute weights for interpolation at several (lat,lon) positions
 # We assume that grid boxes are centered in the grid point.
 #======================
-create_interp_weights <- function(ncfile, locids, lats, lons, region = NULL,
+.create_interp_weights <- function(ncfile, locids, lats, lons, region = NULL,
                                   method = c("nearest", "bilinear", "9point", "invdist4nn", "NE", 
                                              "NW", "SE", "SW")) 
 { 
@@ -307,7 +332,7 @@ create_interp_weights <- function(ncfile, locids, lats, lons, region = NULL,
   #----------------
   # Read grid description and compute (i,j) of requested locations (including decimals)
   #----------------
-  griddes <- get_griddes(nc_cropped1)
+  griddes <- .get_griddes(nc_cropped1)
 
   if (is.null(griddes$yinc)) {
     system(paste0('rm ', nc_cropped1))
@@ -318,7 +343,7 @@ create_interp_weights <- function(ncfile, locids, lats, lons, region = NULL,
   # If latitudes are decreasingly ordered, revert them
   if (griddes$yinc < 0) {
     system(paste0('cdo invertlat ', nc_cropped1, ' ', nc_cropped2))
-    griddes <- get_griddes(nc_cropped2)
+    griddes <- .get_griddes(nc_cropped2)
   }
   # remove temporary files
   system(paste0('rm ', nc_cropped1))
@@ -328,7 +353,7 @@ create_interp_weights <- function(ncfile, locids, lats, lons, region = NULL,
     stop("'griddes' not found in the NetCDF source files")
   } 
   
-  gridpoints <- latlon2ij(griddes, lats, lons)
+  gridpoints <- .latlon2ij(griddes, lats, lons)
 
   #----------------
   # Compute the weights according to the selected method
@@ -573,7 +598,7 @@ create_interp_weights <- function(ncfile, locids, lats, lons, region = NULL,
 # Works only for 'lonlat' and 'gaussian' grids.
 # Grids are supposed to cover whole globe.
 #======================
-latlon2ij <- function(griddes, lats, lons) {
+.latlon2ij <- function(griddes, lats, lons) {
   #------------
   # Check input params
   #------------
@@ -623,7 +648,7 @@ latlon2ij <- function(griddes, lats, lons) {
 #======================
 # Use cdo griddes to obtain grid information
 #======================
-get_griddes <- function(ncfile) {
+.get_griddes <- function(ncfile) {
   tmp <- system(paste0("cdo griddes ", ncfile,
                        " 2>/dev/null | egrep 'gridtype|xsize|ysize|xfirst|xinc|yfirst|yinc'"), intern = T)
   arr <- do.call(rbind, strsplit(tmp,"\\s+= ", perl = T))
@@ -631,7 +656,7 @@ get_griddes <- function(ncfile) {
   names(griddes) <- arr[,1]
 
   if(griddes$gridtype == "gaussian") {
-    griddes$yvals <- get_lats(ncfile)
+    griddes$yvals <- .get_lats(ncfile)
   }
 
   # Convert some fields to numeric. Ensures all fields are present.
@@ -645,7 +670,7 @@ get_griddes <- function(ncfile) {
 #======================
 # Use nco to obtain latitudes. Latitudes shall be named "lat" or "latitude".
 #======================
-get_lats <- function(ncfile) {
+.get_lats <- function(ncfile) {
 
   tmp <- system(paste0('ncks -H -s "%f " -v latitude ',ncfile),intern=T)
 
@@ -663,7 +688,7 @@ get_lats <- function(ncfile) {
 # Uses StartR. All ... parameters go to Start (i.e. specify dat, var, 
 # sdate, time, ensemble, num_procs, etc)
 #======================
-get_model_data <- function(weights.df, mdata) {
+.get_model_data <- function(weights.df, mdata, ncores = NULL) {
  
   #-----------------
   # Get data for all combinations of i and j.
@@ -705,7 +730,9 @@ get_model_data <- function(weights.df, mdata) {
   #-----------------
   # Retrieve with multiApply
   #-----------------
-  sub_mdata <- Apply(mdata, target_dims = list(c(latdim, londim)), fun = function(x) {laply(1:length(is),function(k) { x[js[k],is[k]] }) })$output1
+  sub_mdata <- Apply(mdata, target_dims = list(c(latdim, londim)), 
+                     fun = function(x) {laply(1:length(is),function(k) { x[js[k],is[k]] }) }, 
+                     ncores = ncores)$output1
   names(dim(sub_mdata))[1] <- "gridpoint"
 
   #-----------------
@@ -718,7 +745,7 @@ get_model_data <- function(weights.df, mdata) {
 # Multiply the grid-point series by the weights, 
 # to obtain the desired interpolations
 #======================
-interpolate_data <- function(model_data, weights.df) {
+.interpolate_data <- function(model_data, weights.df, ncores) {
   #-----------------
   # Multiply each gridpoint matrix by its corresponding weight
   #-----------------
@@ -732,9 +759,8 @@ interpolate_data <- function(model_data, weights.df) {
   #interp_data <- apply(weighted_data, -gpdim, function(x) { rowsum(x, weights.df$locid) }) 
   #names(dim(interp_data))[1] <- "location"
   interp_data <- Apply(weighted_data, target_dims = gpdim, fun = function(x) { 
-    rowsum(x, weights.df$locid)}, output_dims = c("location", "aux"))$output1
+    rowsum(x, weights.df$locid)}, output_dims = c("location", "aux"),
+    ncores = ncores)$output1
   
   return(interp_data)
 }
-
-
diff --git a/R/Intlr.R b/R/Intlr.R
index 24c909f..e170138 100644
--- a/R/Intlr.R
+++ b/R/Intlr.R
@@ -68,7 +68,7 @@
 #'border, whereas latmax indicates the upper border. If set to NULL (default), the function
 #'takes the first and last elements of the latitudes and longitudes.
 #'@param ncores an integer indicating the number of cores to use in parallel computation.
-#'
+#'The default value is NULL.
 #'@import multiApply 
 #'
 #'@return A list of three elements. 'data' contains the dowscaled field, 'lat' the
@@ -89,16 +89,16 @@
 CST_Intlr <- function(exp, obs, lr_method, target_grid = NULL, points = NULL, int_method = NULL, 
                       method_point_interp = NULL, predictors = NULL, lat_dim = "lat", lon_dim = "lon", 
                       sdate_dim = "sdate", time_dim = "time", member_dim = "member",
-                      large_scale_predictor_dimname = 'vars', loocv = FALSE, region = NULL, ncores = 1) {
-
+                      large_scale_predictor_dimname = 'vars', loocv = FALSE, region = NULL, ncores = NULL) {
+  
   if (!inherits(exp,'s2dv_cube')) {
     stop("Parameter 'exp' must be of the class 's2dv_cube'")
   }
-
+  
   if (!inherits(obs,'s2dv_cube')) {
     stop("Parameter 'obs' must be of the class 's2dv_cube'")
   }
-
+  
   res <- Intlr(exp = exp$data, obs = obs$data, exp_lats = exp$coords[[lat_dim]], exp_lons = exp$coords[[lon_dim]],
                obs_lats = obs$coords[[lat_dim]], obs_lons = obs$coords[[lon_dim]], points = points, 
                source_file_exp = exp$attrs$source_files[1], source_file_obs = obs$attrs$source_files[1], 
@@ -107,18 +107,18 @@ CST_Intlr <- function(exp, obs, lr_method, target_grid = NULL, points = NULL, in
                lat_dim = lat_dim, lon_dim = lon_dim, sdate_dim = sdate_dim, time_dim = time_dim, 
                member_dim = member_dim, large_scale_predictor_dimname = large_scale_predictor_dimname, 
                loocv = loocv, region = region, ncores = ncores) 
-
+  
   # Modify data, lat and lon in the origina s2dv_cube, adding the downscaled data
   exp$data <- res$data
   exp$dims <- dim(exp$data)
   exp$coords[[lon_dim]] <- res$lon
   exp$coords[[lat_dim]] <- res$lat
-
+  
   obs$data <- res$obs
   obs$dims <- dim(obs$data)
   obs$coords[[lon_dim]] <- res$lon
   obs$coords[[lat_dim]] <- res$lat
-
+  
   res_s2dv <- list(exp = exp, obs = obs)
   return(res_s2dv)
 }
@@ -205,7 +205,7 @@ CST_Intlr <- function(exp, obs, lr_method, target_grid = NULL, points = NULL, in
 #'@param loocv a logical indicating whether to apply leave-one-out cross-validation when
 #'generating the linear regressions. Default to FALSE.
 #'@param ncores an integer indicating the number of cores to use in parallel computation.
-#'
+#'The default value is NULL.
 #'@import multiApply 
 #'
 #'@return A list of three elements. 'data' contains the dowscaled field, 'lat' the
@@ -226,68 +226,68 @@ Intlr <- function(exp, obs, exp_lats, exp_lons, obs_lats, obs_lons, lr_method, t
                   int_method = NULL, method_point_interp = NULL, source_file_exp = NULL, source_file_obs = NULL, 
                   predictors = NULL, lat_dim = "lat", lon_dim = "lon", sdate_dim = "sdate", time_dim = "time", 
                   member_dim = "member", region = NULL, large_scale_predictor_dimname = 'vars', loocv = FALSE, 
-                  ncores = 1) {
-
+                  ncores = NULL) {
+  
   #-----------------------------------
   # Checkings
   #-----------------------------------
   if (!inherits(lr_method, 'character')) {
     stop("Parameter 'lr_method' must be of the class 'character'")
   }
-
+  
   if (!inherits(large_scale_predictor_dimname, 'character')) {
     stop("Parameter 'large_scale_predictor_dimname' must be of the class 'character'")
   }
-
+  
   if (!inherits(loocv, 'logical')) {
     stop("Parameter 'loocv' must be set to TRUE or FALSE")
   }
-
+  
   if (!inherits(lat_dim, 'character')) {
     stop("Parameter 'lat_dim' must be of the class 'character'")
   }
-
+  
   if (!inherits(lon_dim, 'character')) {
     stop("Parameter 'lon_dim' must be of the class 'character'")
   }
-
+  
   if (!inherits(sdate_dim, 'character')) {
     stop("Parameter 'sdate_dim' must be of the class 'character'")
   }
-
+  
   if (!inherits(large_scale_predictor_dimname, 'character')) {
     stop("Parameter 'large_scale_predictor_dimname' must be of the class 'character'")
   }
-
+  
   if (is.na(match(lon_dim, names(dim(exp)))) | is.na(match(lon_dim, names(dim(obs))))) {
     stop("Missing longitude dimension in 'exp' and/or 'obs', or does not match the parameter ",
          "'lon_dim'")
   }
-
+  
   if (is.na(match(lat_dim, names(dim(exp)))) | is.na(match(lat_dim, names(dim(obs))))) {
     stop("Missing latitude dimension in 'exp' and/or 'obs', or does not match the parameter ",
          "'lat_dim'")
   }
-
+  
   if (is.na(match(sdate_dim, names(dim(exp)))) | is.na(match(sdate_dim, names(dim(obs))))) {
     stop("Missing start date dimension in 'exp' and/or 'obs', or does not match the parameter ",
          "'sdate_dim'")
   }
-
+  
   if (!is.null(points) & (is.null(source_file_exp) | is.null(source_file_obs))) {
     stop("No source files found. Source files for exp and obs must be provided in the parameters ",
          "'source_file_exp' and 'source_file_obs', respectively.")
   }
-
+  
   if (!is.null(points) & is.null(method_point_interp)) {
     stop("Please provide the interpolation method to interpolate gridded data to point locations ",
          "through the parameter 'method_point_interp'.")
   }
-
+  
   # sdate must be the time dimension in the input data
   stopifnot(sdate_dim %in% names(dim(exp)))
   stopifnot(sdate_dim %in% names(dim(obs))) 
-
+  
   # the code is not yet prepared to handle members in the observations
   restore_ens <- FALSE
   if (member_dim %in% names(dim(obs))) {       
@@ -299,7 +299,7 @@ Intlr <- function(exp, obs, exp_lats, exp_lons, obs_lats, obs_lons, lr_method, t
            "but it should be of length = 1).")       
     }     
   }
- 
+  
   # checkings for the parametre 'predictors' 
   if (!is.null(predictors)) {
     if (!is.array(predictors)) {
@@ -311,29 +311,36 @@ Intlr <- function(exp, obs, exp_lats, exp_lons, obs_lats, obs_lons, lr_method, t
       stopifnot(dim(predictors)[sdate_dim] == dim(exp)[sdate_dim]) 
     }
   }
-
+  ## ncores
+  if (!is.null(ncores)) {
+    if (!is.numeric(ncores) | any(ncores %% 1 != 0) | any(ncores < 0) |
+        length(ncores) > 1) {
+      stop("Parameter 'ncores' must be a positive integer.")
+    }
+  }
+  
   #-----------------------------------
   # Interpolation
   #-----------------------------------
   if (lr_method != '4nn') {
- 
+    
     if (is.null(int_method)) {
       stop("Parameter 'int_method' must be a character vector indicating the interpolation method. ",
            "Accepted methods are con, bil, bic, nn, con2")
     }
- 
+    
     if (is.null(region)) {
       warning("The borders of the downscaling region have not been provided. Assuming the ",
               "four borders of the downscaling region are defined by the first and last ",
               "elements of the parametres 'obs_lats' and 'obs_lons'.")      
       region <- c(obs_lons[1], obs_lons[length(obs_lons)], obs_lats[1], obs_lats[length(obs_lats)])
     }
-
+    
     exp_interpolated <- Interpolation(exp = exp, lats = exp_lats, lons = exp_lons, target_grid = target_grid, 
                                       points = points, method_point_interp = method_point_interp, 
                                       source_file = source_file_exp, lat_dim = lat_dim, lon_dim = lon_dim, 
-                                      method_remap = int_method, region = region)
-
+                                      method_remap = int_method, region = region, ncores = ncores)
+    
     # If after interpolating 'exp' data the coordinates do not match, the obs data is interpolated to
     # the same grid to force the matching
     if ((!.check_coords(lat1 = exp_interpolated$lat, lat2 = obs_lats, 
@@ -341,8 +348,8 @@ Intlr <- function(exp, obs, exp_lats, exp_lons, obs_lats, obs_lons, lr_method, t
       obs_interpolated <- Interpolation(exp = obs, lats = obs_lats, lons = obs_lons, target_grid = target_grid, 
                                         points = points, method_point_interp = method_point_interp, 
                                         source_file = source_file_obs, lat_dim = lat_dim, lon_dim = lon_dim, 
-                                        method_remap = int_method, region = region)
-
+                                        method_remap = int_method, region = region, ncores = ncores)
+      
       lats <- obs_interpolated$lat
       lons <- obs_interpolated$lon 
       obs_interpolated <- obs_interpolated$data
@@ -352,7 +359,7 @@ Intlr <- function(exp, obs, exp_lats, exp_lons, obs_lats, obs_lons, lr_method, t
       lons <- obs_lons
     }
   }
- 
+  
   #-----------------------------------
   # Linear regressions
   #-----------------------------------
@@ -366,7 +373,7 @@ Intlr <- function(exp, obs, exp_lats, exp_lons, obs_lats, obs_lons, lr_method, t
     target_dims_predictor <- sdate_dim
     target_dims_predictand <- sdate_dim
   } 
-
+  
   # (Multi) linear regression with large-scale predictors
   # Predictor: passed through the parameter 'predictors' by the user. Can be model or observations
   # Predictand: model data
@@ -374,27 +381,27 @@ Intlr <- function(exp, obs, exp_lats, exp_lons, obs_lats, obs_lons, lr_method, t
     if (is.null(predictors)) {
       stop("The large-scale predictors must be passed through the parametre 'predictors'")
     }  
-
+    
     predictand <- obs_interpolated
     predictor <- predictors
     
     target_dims_predictor <- c(sdate_dim, large_scale_predictor_dimname)
     target_dims_predictand <- sdate_dim
   } 
-
+  
   # Multi-linear regression with the four nearest neighbours
   # Predictors: model data
   # Predictand: observations
   else if (lr_method == '4nn') {
-    predictor <- find_nn(coar = exp, lats_hres = obs_lats, lons_hres = obs_lons, lats_coar = exp_lats, 
-                         lons_coar = exp_lons, lat_dim = lat_dim, lon_dim = lon_dim, nn = 4)
-     
+    predictor <- .find_nn(coar = exp, lats_hres = obs_lats, lons_hres = obs_lons, lats_coar = exp_lats, 
+                         lons_coar = exp_lons, lat_dim = lat_dim, lon_dim = lon_dim, nn = 4, ncores = ncores)
+    
     if (is.null(points)) {
       if (!is.null(target_grid)) {
         warning("Interpolating to the 'obs' grid")
       }
       predictand <- obs
-    
+      
       lats <- obs_lats
       lons <- obs_lons
     }
@@ -403,40 +410,40 @@ Intlr <- function(exp, obs, exp_lats, exp_lons, obs_lats, obs_lons, lr_method, t
     else {
       predictor <- Interpolation(exp = predictor, lats = obs_lats, lons = obs_lons, target_grid = NULL,
                                  points = points, method_point_interp = method_point_interp, 
-                                 source_file = source_file_obs, method_remap = NULL, region = region)
-
+                                 source_file = source_file_obs, method_remap = NULL, region = region, ncores = ncores)
+      
       predictand <- Interpolation(exp = obs, lats = obs_lats, lons = obs_lons, target_grid = NULL,
                                   points = points, method_point_interp = method_point_interp, 
-                                  source_file = source_file_obs, method_remap = NULL, region = region)
-           
+                                  source_file = source_file_obs, method_remap = NULL, region = region, ncores = ncores)
+      
       lats <- predictor$lat
       lons <- predictor$lon 
       predictor <- predictor$data
       predictand <- predictand$data
     }
-  
+    
     target_dims_predictor <- c(sdate_dim,'nn')
     target_dims_predictand <- sdate_dim
   }
-
+  
   else {
     stop(paste0(lr_method, " method is not implemented yet"))
   }
-
+  
   # Apply the linear regressions
   res <- Apply(list(predictor, predictand), target_dims = list(target_dims_predictor, target_dims_predictand), 
                fun = .intlr, loocv = loocv, ncores = ncores)$output1
- 
+  
   names(dim(res))[1] <- sdate_dim
-
+  
   # restore ensemble dimension in observations if it existed originally
   if (restore_ens) {
     predictand <- s2dv::InsertDim(predictand, posdim = 1, lendim = 1, name = member_dim)
   }
-
+  
   # Return a list of three elements
   res <- list(data = res, obs = predictand, lon = lons, lat = lats)
-
+  
   return(res)
 }
 
@@ -444,63 +451,70 @@ Intlr <- function(exp, obs, exp_lats, exp_lons, obs_lats, obs_lons, lr_method, t
 # Atomic function to generate and apply the linear regressions
 #-----------------------------------
 .intlr <- function(x, y, loocv) {
-
+  
   tmp_df <- data.frame(x = x, y = y)  
-
   # if the data is all NA, force return return NA
   if (all(is.na(tmp_df)) | (sum(apply(tmp_df, 2, function(x) !all(is.na(x)))) == 1)) {
-
+    
     n <- nrow(tmp_df)
     res <- rep(NA, n)
-  
+    
   } else {
     # training
-    lm1 <- train_lm(df = tmp_df, loocv = loocv)  
-
+    lm1 <- .train_lm(df = tmp_df, loocv = loocv)  
+    
     # prediction
-    res <- pred_lm(lm1 = lm1, df = tmp_df, loocv = loocv)
+    res <- .pred_lm(lm1 = lm1, df = tmp_df, loocv = loocv)
   }
   
   return(res)
- 
 }
 
 #-----------------------------------
 # Function to generate the linear regressions.
 # Returns a list
 #-----------------------------------
-train_lm <- function(df, loocv) {
-
-  # Remove columns containing only NA's
-  df <- df[ , apply(df, 2, function(x) !all(is.na(x)))]
-
+.train_lm <- function(df, loocv) {
+  
+  # Remove predictor columns containing only NA's
+  df <- df[ , apply(df[,colnames(df) != 'y'], 2, function(x) !all(is.na(x)))]
+  
   if (loocv) {
-   
-    lm1 <- lapply(1:nrow(df), function(j) lm(df[-j,], formula = y ~ .))
- 
+    
+    lm1 <- lapply(1:nrow(df), function(j) {
+      if (all(is.na(df[-j,]$y))) {
+        return(NA)
+      } else {
+        return(lm(df[-j,], formula = y ~ .)) 
+      }})
   } else {
-
-    lm1 <- list(lm(data = df, formula = y ~ .))
+    
+    lm1 <- ifelse(all(is.na(df$y)), NA, list(lm(data = df, formula = y ~ .)))
   }
-
+  
   return(lm1)
 }
 
 #-----------------------------------
 # Function to apply the linear regressions.
 #-----------------------------------
-pred_lm <- function(df, lm1, loocv) {
-
-  if (loocv) {
-
-    pred_vals <- sapply(1:nrow(df), function(j) predict(lm1[[j]], df[j,]))
+.pred_lm <- function(df, lm1, loocv) {
   
+  if (loocv) {
+    pred_vals <- sapply(1:nrow(df), function(j) {
+      if (all(is.na(lm1[[j]]))) {
+        return(NA)
+      } else {
+        return(predict(lm1[[j]], df[j,]))
+      }})
   } else {
-  
-    pred_vals_ls <- lapply(lm1, predict, data = df)    
-    pred_vals <- unlist(pred_vals_ls)
+    if (!is.na(lm1)) { 
+      pred_vals_ls <- lapply(lm1, predict, data = df)    
+      pred_vals <- unlist(pred_vals_ls)
+    } else {
+      pred_vals <- rep(NA, nrow(df))
+    }
   } 
- 
   return(pred_vals)
 }
 
@@ -508,20 +522,20 @@ pred_lm <- function(df, lm1, loocv) {
 # Function to find N nearest neighbours.
 # 'coar' is an array with named dimensions
 #-----------------------------------
-find_nn <- function(coar, lats_hres, lons_hres, lats_coar, lons_coar, lat_dim, lon_dim, nn = 4) {
-
+.find_nn <- function(coar, lats_hres, lons_hres, lats_coar, lons_coar, lat_dim, lon_dim, nn = 4, ncores = NULL) {
+  
   # Sort the distances from closest to furthest
   idx_lat <- as.array(sapply(lats_hres, function(x) order(abs(lats_coar - x))[1:nn]))
   idx_lon <- as.array(sapply(lons_hres, function(x) order(abs(lons_coar - x))[1:nn]))
-
+  
   names(dim(idx_lat)) <- c('nn', lat_dim)
   names(dim(idx_lon)) <- c('nn', lon_dim)
-
+  
   # obtain the values of the nearest neighbours 
   nearest <- Apply(list(coar, idx_lat, idx_lon), 
                    target_dims = list(c(lat_dim, lon_dim), lat_dim, lon_dim), 
-                   fun = function(x, y, z) x[y, z])$output1
-
+                   fun = function(x, y, z) x[y, z], ncores = ncores)$output1
+  
   return(nearest)
 }
 
diff --git a/R/LogisticReg.R b/R/LogisticReg.R
index c514d25..4581fc8 100644
--- a/R/LogisticReg.R
+++ b/R/LogisticReg.R
@@ -62,8 +62,8 @@
 #'takes the first and last elements of the latitudes and longitudes.
 #'@param loocv a logical vector indicating whether to perform leave-one-out cross-validation
 #'in the fitting of the logistic regression. Default to FALSE.
-#'@param ncores an integer indicating the number of cores to use in parallel computation.
-#'
+#'@param ncores an integer indicating the number of cores to use in parallel computation. 
+#'The default value is NULL.
 #'@import multiApply
 #'@import nnet
 #'@importFrom laply plyr  
@@ -91,7 +91,7 @@
 CST_LogisticReg <- function(exp, obs, target_grid, int_method = NULL, log_reg_method = "ens_mean",
                   probs_cat = c(1/3,2/3), return_most_likely_cat = FALSE, points = NULL, 
                   method_point_interp = NULL, lat_dim = "lat", lon_dim = "lon", sdate_dim = "sdate", 
-                  member_dim = "member", region = NULL, loocv = FALSE, ncores = 1) {
+                  member_dim = "member", region = NULL, loocv = FALSE, ncores = NULL) {
 
   if (!inherits(exp,'s2dv_cube')) {
     stop("Parameter 'exp' must be of the class 's2dv_cube'")
@@ -198,8 +198,8 @@ CST_LogisticReg <- function(exp, obs, target_grid, int_method = NULL, log_reg_me
 #'takes the first and last elements of the latitudes and longitudes.
 #'@param loocv a logical vector indicating whether to perform leave-one-out cross-validation
 #'in the fitting of the logistic regression. Default to FALSE.
-#'@param ncores an integer indicating the number of cores to use in parallel computation.
-#'
+#'@param ncores an integer indicating the number of cores to use in parallel computation. 
+#'The default value is NULL.
 #'@import multiApply
 #'@import nnet
 #'@importFrom laply plyr  
@@ -227,7 +227,7 @@ LogisticReg <- function(exp, obs, exp_lats, exp_lons, obs_lats, obs_lons, target
                   int_method = NULL, log_reg_method = "ens_mean", probs_cat = c(1/3,2/3), 
                   return_most_likely_cat = FALSE, points = NULL, method_point_interp = NULL, 
                   lat_dim = "lat", lon_dim = "lon", sdate_dim = "sdate", member_dim = "member", 
-                  source_file = NULL, region = NULL, loocv = FALSE, ncores = 1) {
+                  source_file = NULL, region = NULL, loocv = FALSE, ncores = NULL) {
 
   #-----------------------------------
   # Checkings
@@ -303,7 +303,14 @@ LogisticReg <- function(exp, obs, exp_lats, exp_lons, obs_lats, obs_lons, target
             "'obs_lats' and 'obs_lons'.")
     region <- c(obs_lons[1], obs_lons[length(obs_lons)], obs_lats[1], obs_lats[length(obs_lats)])
   }
-
+  ## ncores
+  if (!is.null(ncores)) {
+    if (!is.numeric(ncores) | any(ncores %% 1 != 0) | any(ncores < 0) |
+        length(ncores) > 1) {
+      stop("Parameter 'ncores' must be a positive integer.")
+    }
+  }
+  
   # the code is not yet prepared to handle members in the observations
   restore_ens <- FALSE
   if (member_dim %in% names(dim(obs))) {
@@ -319,18 +326,18 @@ LogisticReg <- function(exp, obs, exp_lats, exp_lons, obs_lats, obs_lons, target
   exp_interpolated <- Interpolation(exp = exp, lats = exp_lats, lons = exp_lons, target_grid = target_grid,
                                     method_remap = int_method, points = points, source_file = source_file,
                                     lat_dim = lat_dim, lon_dim = lon_dim, method_point_interp = method_point_interp,
-                                    region = region)
+                                    region = region, ncores = ncores)
 
   # compute ensemble mean anomalies
   if (log_reg_method == "ens_mean") {
-    predictor <- get_ens_mean_anom(obj_ens = exp_interpolated$data, member_dim = member_dim, sdate_dim = sdate_dim)
+    predictor <- .get_ens_mean_anom(obj_ens = exp_interpolated$data, member_dim = member_dim, sdate_dim = sdate_dim, ncores = ncores)
 
     target_dims_predictor <- sdate_dim
   } 
   else if (log_reg_method == "ens_mean_sd")  {
-    ens_mean_anom <- get_ens_mean_anom(obj_ens = exp_interpolated$data, member_dim = member_dim, 
-                                       sdate_dim = sdate_dim)
-    ens_sd <- get_ens_sd(obj_ens = exp_interpolated$data, member_dim = member_dim)
+    ens_mean_anom <- .get_ens_mean_anom(obj_ens = exp_interpolated$data, member_dim = member_dim, 
+                                       sdate_dim = sdate_dim, ncores = ncores)
+    ens_sd <- .get_ens_sd(obj_ens = exp_interpolated$data, member_dim = member_dim, ncores = ncores)
 
     #merge two arrays into one array of predictors
     predictor <- abind(ens_mean_anom, ens_sd, along = 1/2)
@@ -338,7 +345,7 @@ LogisticReg <- function(exp, obs, exp_lats, exp_lons, obs_lats, obs_lons, target
      
     target_dims_predictor <- c(sdate_dim, "pred")
   } else if (log_reg_method == "sorted_members") {   
-    predictor <- sort_members(obj_ens = exp_interpolated$data, member_dim = member_dim)
+    predictor <- .sort_members(obj_ens = exp_interpolated$data, member_dim = member_dim, ncores = ncores)
 
     target_dims_predictor <- c(sdate_dim, member_dim)
   } else {
@@ -352,7 +359,7 @@ LogisticReg <- function(exp, obs, exp_lats, exp_lons, obs_lats, obs_lons, target
     obs_interpolated <- Interpolation(exp = obs, lats = obs_lats, lons = obs_lons, target_grid = target_grid,
                                       method_remap = int_method, points = points, source_file = source_file,
                                       lat_dim = lat_dim, lon_dim = lon_dim,
-                                      method_point_interp = method_point_interp, region = region)
+                                      method_point_interp = method_point_interp, region = region, ncores = ncores)
     obs_ref <- obs_interpolated$data
   } else {
     obs_ref <- obs
@@ -362,19 +369,19 @@ LogisticReg <- function(exp, obs, exp_lats, exp_lons, obs_lats, obs_lons, target
   obs_cat <- Apply(obs_ref, target_dims = sdate_dim, function(x) {
                    terc <- convert2prob(as.vector(x), prob = probs_cat) 
                    apply(terc, 1, function(r) which (r == 1))},          
-               output_dims = sdate_dim)$output1
+               output_dims = sdate_dim, ncores = ncores)$output1
 
   res <- Apply(list(predictor, obs_cat), target_dims = list(target_dims_predictor, sdate_dim), fun = function(x, y)
-               .log_reg(x = x, y = y, loocv = loocv), output_dims = c(sdate_dim, "category"))$output1
+               .log_reg(x = x, y = y, loocv = loocv), output_dims = c(sdate_dim, "category"), ncores = ncores)$output1
 
   if (return_most_likely_cat) {
-    res <- Apply(res, target_dims = c(sdate_dim, "category"), most_likely_category,
-                 output_dims = sdate_dim)$output1
+    res <- Apply(res, target_dims = c(sdate_dim, "category"), .most_likely_category,
+                 output_dims = sdate_dim, ncores = ncores)$output1
   }
 
   # restore ensemble dimension in observations if it existed originally
   if (restore_ens) {
-    obs_ref <- s2dv::InsertDim(obs_ref, posdim = 1, lendim = 1, name = member_dim)
+    obs_ref <- s2dv::InsertDim(obs_ref, posdim = 1, lendim = 1, name = member_dim, ncores = ncores)
   }
 
   res <- list(data = res, obs = obs_ref, lon = exp_interpolated$lon, lat = exp_interpolated$lat)
@@ -382,7 +389,7 @@ LogisticReg <- function(exp, obs, exp_lats, exp_lons, obs_lats, obs_lons, target
   return(res)
 }
 
-most_likely_category <- function(data) {
+.most_likely_category <- function(data) {
 # data, expected dims: start date, category (in this order)
 
   if (all(is.na(data))) {
@@ -393,33 +400,33 @@ most_likely_category <- function(data) {
   return(mlc)
 }
 
-sort_members <- function(obj_ens, member_dim) {
+.sort_members <- function(obj_ens, member_dim, ncores = NULL) {
 
-  sorted <- Apply(obj_ens, target_dims = member_dim, sort, decreasing = TRUE, na.last = TRUE)$output1  
+  sorted <- Apply(obj_ens, target_dims = member_dim, sort, decreasing = TRUE, na.last = TRUE, ncores = ncores)$output1  
 
   return(sorted)
 }
 
-get_ens_sd <- function(obj_ens, member_dim) {
+.get_ens_sd <- function(obj_ens, member_dim, ncores = NULL) {
 
   # compute ensemble spread
-  ens_sd <- Apply(obj_ens, target_dims = member_dim, sd, na.rm = TRUE)$output1
+  ens_sd <- Apply(obj_ens, target_dims = member_dim, sd, na.rm = TRUE, ncores = ncores)$output1
 
   return(ens_sd)
 }
 
-get_ens_mean_anom <- function(obj_ens, member_dim, sdate_dim) {
+.get_ens_mean_anom <- function(obj_ens, member_dim, sdate_dim, ncores = NULL) {
 
   require(s2dv)
 
   # compute climatology
-  clim <- Apply(obj_ens, target_dims = c(member_dim, sdate_dim), mean)$output1
+  clim <- Apply(obj_ens, target_dims = c(member_dim, sdate_dim), mean, ncores = ncores)$output1
 
   # compute ensemble mean
-  ens_mean <- Apply(obj_ens, target_dims = member_dim, mean, na.rm = TRUE)$output1
+  ens_mean <- Apply(obj_ens, target_dims = member_dim, mean, na.rm = TRUE, ncores = ncores)$output1
 
   # compute ensemble mean anomalies
-  anom <- Ano(ens_mean, clim)
+  anom <- Ano(ens_mean, clim, ncores = ncores)
 
   return(anom)
 }
@@ -437,10 +444,10 @@ get_ens_mean_anom <- function(obj_ens, member_dim, sdate_dim) {
 
   } else {
     # training
-    lm1 <- train_lr(df = tmp_df, loocv = loocv)
+    lm1 <- .train_lr(df = tmp_df, loocv = loocv)
 
     # prediction
-    res <- pred_lr(lm1 = lm1, df = tmp_df, loocv = loocv)
+    res <- .pred_lr(lm1 = lm1, df = tmp_df, loocv = loocv)
   }
   return(res)
 }
@@ -448,7 +455,7 @@ get_ens_mean_anom <- function(obj_ens, member_dim, sdate_dim) {
 #-----------------------------------
 # Function to train the logistic regressions.
 #-----------------------------------
-train_lr <- function(df, loocv) {
+.train_lr <- function(df, loocv) {
 
   require(nnet)
   
@@ -471,7 +478,7 @@ train_lr <- function(df, loocv) {
 #-----------------------------------
 # Function to apply the logistic regressions.
 #-----------------------------------
-pred_lr <- function(df, lm1, loocv) {
+.pred_lr <- function(df, lm1, loocv) {
 
   require(plyr)
 
@@ -496,5 +503,3 @@ pred_lr <- function(df, lm1, loocv) {
 
   return(pred_vals)
 }
-
-
-- 
GitLab