diff --git a/NAMESPACE b/NAMESPACE
index 29dc0811bb0d77faefcd9658edd5343a3bb2b944..91e3447127b00a81a074d786c35a9a165b08082c 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -25,6 +25,7 @@ export(CST_SaveExp)
 export(CST_SplitDim)
 export(CST_WeatherRegimes)
 export(Calibration)
+export(CategoricalEnsCombination)
 export(EnsClustering)
 export(MergeDims)
 export(MultiEOF)
@@ -101,7 +102,6 @@ importFrom(s2dverification,Ano_CrossValid)
 importFrom(s2dverification,Clim)
 importFrom(s2dverification,EOF)
 importFrom(s2dverification,Eno)
-importFrom(s2dverification,InsertDim)
 importFrom(s2dverification,Load)
 importFrom(s2dverification,Mean1Dim)
 importFrom(s2dverification,Subset)
diff --git a/NEWS.md b/NEWS.md
index 3a8dedb39807f2c7ed79009bffc864d8f54b1b8c..7e06bfcd2e321bbbeed066abcb7fba4a1b0ca1a4 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -8,6 +8,8 @@
     + CST_RainFARM accepts multi-dimensional slopes and weights and handless missing values in sample dimensions.
     + CST_MultiMetric includes 'rpss' metric and it is addapted to s2dv.
     + PlotMostLikelyQuantileMap vignette
+    + PlotTriangles4Categories includes two parameters to adjust axis and margins.
+    + CategoricalEnsCombination is exposed to users
 
 - Fixes:
     + PlotForecastPDF correctly displays terciles labels 
diff --git a/R/CST_Anomaly.R b/R/CST_Anomaly.R
index ef616225fe3d8bf827325982937a3c13428f2283..52a786fb577d5a4f397a9a205abe65e1ece46f41 100644
--- a/R/CST_Anomaly.R
+++ b/R/CST_Anomaly.R
@@ -18,7 +18,8 @@
 #'
 #' @return A list with two S3 objects, 'exp' and 'obs', of the class 's2dv_cube', containing experimental and date-corresponding observational anomalies, respectively. These 's2dv_cube's can be ingested by other functions in CSTools.
 #'
-#'@importFrom s2dverification Clim Ano_CrossValid InsertDim
+#'@importFrom s2dverification Clim Ano_CrossValid
+#'@importFrom s2dv InsertDim
 #'
 #'@examples
 #'# Example 1:
@@ -165,12 +166,12 @@ CST_Anomaly <- function(exp = NULL, obs = NULL, cross = FALSE, memb = TRUE,
       clim_exp <- tmp$clim_exp
       clim_obs <- tmp$clim_obs
     } else {
-      clim_exp <- s2dverification::InsertDim(tmp$clim_exp, 2, dim_exp[2]) 
-      clim_obs <- s2dverification::InsertDim(tmp$clim_obs, 2, dim_obs[2]) 
+      clim_exp <- InsertDim(tmp$clim_exp, 2, dim_exp[2]) 
+      clim_obs <- InsertDim(tmp$clim_obs, 2, dim_obs[2]) 
     }
    
-    clim_exp <- s2dverification::InsertDim(clim_exp, 3, dim_exp[3]) 
-    clim_obs <- s2dverification::InsertDim(clim_obs, 3, dim_obs[3]) 
+    clim_exp <- InsertDim(clim_exp, 3, dim_exp[3]) 
+    clim_obs <- InsertDim(clim_obs, 3, dim_obs[3]) 
     ano <- NULL    
     ano$ano_exp <- exp$data - clim_exp
     ano$ano_obs <- obs$data - clim_obs 
diff --git a/R/CST_Calibration.R b/R/CST_Calibration.R
index fd381790f936697fff109766e69b4aa847a77b73..874c351eab65459da83bd6e705e1902b8cbb91bf 100644
--- a/R/CST_Calibration.R
+++ b/R/CST_Calibration.R
@@ -13,7 +13,7 @@
 #'@param ncores is an integer that indicates the number of cores for parallel computations using multiApply function. The default value is one.
 #'@return an object of class \code{s2dv_cube} containing the calibrated forecasts in the element \code{$data} with the same dimensions as the one in the exp object.
 #'
-#'@importFrom s2dverification InsertDim
+#'@importFrom s2dv InsertDim
 #'@import abind
 #'
 #'@seealso \code{\link{CST_Load}}
@@ -84,7 +84,7 @@ CST_Calibration <- function(exp, obs, cal.method = "mse_min",
 #'@param ncores is an integer that indicates the number of cores for parallel computations using multiApply function. The default value is one.
 #'@return an array containing the calibrated forecasts with the same dimensions as the \code{exp} array.
 #'
-#'@importFrom s2dverification InsertDim
+#'@importFrom s2dv InsertDim
 #'@import abind
 #'
 #'@seealso \code{\link{CST_Load}}
@@ -251,7 +251,7 @@ Calibration <- function(exp, obs, cal.method = "mse_min", eval.method = "leave-o
 
 .calc.obs.fc.quant <- function(obs, fc){ #function to calculate different quantities of a series of ensemble forecasts and corresponding observations
   amt.mbr <- dim(fc)[1]
-  obs.per.ens <- s2dverification::InsertDim(obs, posdim = 1, lendim = amt.mbr)
+  obs.per.ens <- InsertDim(obs, posdim = 1, lendim = amt.mbr)
   fc.ens.av <- apply(fc, c(2), mean, na.rm = TRUE)
   cor.obs.fc <- cor(fc.ens.av, obs, use = "complete.obs")
   obs.av <- mean(obs, na.rm = TRUE)
@@ -271,7 +271,7 @@ Calibration <- function(exp, obs, cal.method = "mse_min", eval.method = "leave-o
 
 .calc.obs.fc.quant.ext <- function(obs, fc){ #extended function to calculate different quantities of a series of ensemble forecasts and corresponding observations
   amt.mbr <- dim(fc)[1]
-  obs.per.ens <- s2dverification::InsertDim(obs, posdim = 1, lendim = amt.mbr)
+  obs.per.ens <- InsertDim(obs, posdim = 1, lendim = amt.mbr)
   fc.ens.av <- apply(fc, c(2), mean, na.rm = TRUE)
   cor.obs.fc <- cor(fc.ens.av, obs, use = "complete.obs")
   obs.av <- mean(obs, na.rm = TRUE)
@@ -296,7 +296,7 @@ Calibration <- function(exp, obs, cal.method = "mse_min", eval.method = "leave-o
   fc.ens.av <- apply(fc, c(2), mean, na.rm = TRUE)
   fc.ens.av.av <- mean(fc.ens.av, na.rm = TRUE)
   fc.ens.av.sd <- sd(fc.ens.av, na.rm = TRUE)
-  fc.ens.av.per.ens <- s2dverification::InsertDim(fc.ens.av, posdim = 1, lendim = amt.mbr)
+  fc.ens.av.per.ens <- InsertDim(fc.ens.av, posdim = 1, lendim = amt.mbr)
   fc.ens.sd <- apply(fc, c(2), sd, na.rm = TRUE)
   fc.ens.var.av.sqrt <- sqrt(mean(fc.ens.sd^2, na.rm = TRUE))
   fc.dev <- fc - fc.ens.av.per.ens
@@ -322,10 +322,10 @@ Calibration <- function(exp, obs, cal.method = "mse_min", eval.method = "leave-o
 .calc.fc.quant.ext <- function(fc){ #extended function to calculate different quantities of a series of ensemble forecasts
 
   amt.mbr <- dim(fc)[1]
-  repmat1.tmp <- s2dverification::InsertDim(fc, posdim = 1, lendim = amt.mbr)
+  repmat1.tmp <- InsertDim(fc, posdim = 1, lendim = amt.mbr)
   repmat2.tmp <- aperm(repmat1.tmp, c(2, 1, 3))
   spr.abs <- apply(abs(repmat1.tmp - repmat2.tmp), c(3), mean, na.rm = TRUE)
-  spr.abs.per.ens <- s2dverification::InsertDim(spr.abs, posdim = 1, lendim = amt.mbr)
+  spr.abs.per.ens <- InsertDim(spr.abs, posdim = 1, lendim = amt.mbr)
 
   return(
     append(.calc.fc.quant(fc),
diff --git a/R/CST_CategoricalEnsCombination.R b/R/CST_CategoricalEnsCombination.R
index 26ac54492112884c660ec3f4495aaed0205a8090..909792ee593606a4c2e3b9e5c5afb551cadff206 100644
--- a/R/CST_CategoricalEnsCombination.R
+++ b/R/CST_CategoricalEnsCombination.R
@@ -55,7 +55,7 @@
 #'@references Robertson, A. W., Lall, U., Zebiak, S. E., & Goddard, L. (2004). Improved combination of multiple atmospheric GCM ensembles for seasonal prediction. Monthly Weather Review, 132(12), 2732-2744.
 #'@references Van Schaeybroeck, B., & Vannitsem, S. (2019). Postprocessing of Long-Range Forecasts. In Statistical Postprocessing of Ensemble Forecasts (pp. 267-290).
 #'
-#'@importFrom s2dverification InsertDim
+#'@importFrom s2dv InsertDim
 #'@import abind
 #'@examples 
 #'
@@ -86,19 +86,45 @@ CST_CategoricalEnsCombination <- function(exp, obs, cat.method = "pool", eval.me
          "of the parameter 'obs' must be equal to 1.")
   }
   names.dim.tmp <- names(dim(exp$data))
-  exp$data <- .CategoricalEnsCombination.wrap(fc = exp$data, obs = obs$data, cat.method = cat.method, eval.method = eval.method, amt.cat = amt.cat, ...)
+  exp$data <- CategoricalEnsCombination(fc = exp$data, obs = obs$data, cat.method = cat.method,
+                                        eval.method = eval.method, amt.cat = amt.cat, ...)
   names.dim.tmp[which(names.dim.tmp == "member")] <- "category"
   names(dim(exp$data)) <- names.dim.tmp
-  exp$data <- s2dverification::InsertDim(exp$data, lendim = 1, posdim = 2)
+  exp$data <- suppressWarnings(InsertDim(exp$data, lendim = 1, posdim = 2))
   names(dim(exp$data))[2] <- "member"
   exp$Datasets <- c(exp$Datasets, obs$Datasets)
   exp$source_files <- c(exp$source_files, obs$source_files)
   return(exp)
 }
 
+#' Make categorical forecast based on a multi-model forecast with potential for calibrate
+#'
+#'@author Bert Van Schaeybroeck, \email{bertvs@meteo.be}
+#'@description This function converts a multi-model ensemble forecast 
+#' into a categorical forecast by giving the probability
+#' for each category. Different methods are available to combine 
+#' the different ensemble forecasting models into 
+#' probabilistic categorical forecasts. 
+#'
+#' See details in ?CST_CategoricalEnsCombination
+#'@param fc a multi-dimensional array with named dimensions containing the seasonal forecast experiment data in the element named \code{$data}. The amount of forecasting models is equal to the size of the \code{dataset} dimension of the data array. The amount of members per model may be different. The  size of the \code{member} dimension of the data array is equal to the maximum of the ensemble members among the models. Models with smaller ensemble sizes have residual indices of \code{member} dimension in the data array filled with NA values.
+#'@param obs a multidimensional array with named dimensions containing the observed data in the element named \code{$data}.
+#'@param amt.cat is the amount of categories. Equally-sized quantiles will be calculated based on the amount of categories.
+#'@param cat.method method used to produce the categorical forecast, can be either \code{pool}, \code{comb}, \code{mmw} or \code{obs}. The method pool assumes equal weight for all ensemble members while the method comb assumes equal weight for each model. The weighting method is descirbed in Rajagopalan et al. (2002), Robertson et al. (2004) and Van Schaeybroeck and Vannitsem (2019). Finally, the \code{obs} method classifies the observations into the different categories and therefore contains only 0 and 1 values. 
+#'@param eval.method is the sampling method used, can be either \code{"in-sample"} or \code{"leave-one-out"}. Default value is the \code{"leave-one-out"} cross validation. 
+#'@param ... other parameters to be passed on to the calibration procedure.
+#'
+#'@return an array containing the categorical forecasts in the element called \code{$data}. The first two dimensions of the returned object are named dataset and member and are both of size one. An additional dimension named category is introduced and is of size amt.cat.
+#'
+#'@references Rajagopalan, B., Lall, U., & Zebiak, S. E. (2002). Categorical climate forecasts through regularization and optimal combination of multiple GCM ensembles. Monthly Weather Review, 130(7), 1792-1811.
+#'@references Robertson, A. W., Lall, U., Zebiak, S. E., & Goddard, L. (2004). Improved combination of multiple atmospheric GCM ensembles for seasonal prediction. Monthly Weather Review, 132(12), 2732-2744.
+#'@references Van Schaeybroeck, B., & Vannitsem, S. (2019). Postprocessing of Long-Range Forecasts. In Statistical Postprocessing of Ensemble Forecasts (pp. 267-290).
+#'
+#'@importFrom s2dv InsertDim
+#'@import abind
+#'@export
 
-
-.CategoricalEnsCombination.wrap <- function (fc, obs, cat.method, eval.method, amt.cat, ...) {
+CategoricalEnsCombination <- function (fc, obs, cat.method, eval.method, amt.cat, ...) {
   
   if (!all(c("member", "sdate") %in% names(dim(fc)))) {
     stop("Parameter 'exp' must have the dimensions 'member' and 'sdate'.")
@@ -276,10 +302,10 @@ comb.dims <- function(arr.in, dims.to.combine){
       optim.tmp <- constrOptim(theta = init.par, f = .funct.optim, grad = .funct.optim.grad, 
                           ui = constr.mtrx, ci = constr.vec,
                           freq.per.mdl.at.obs = freq.per.mdl.at.obs) 
-      init.par <- optim.tmp$par * (1 - abs(rnorm(amt.coeff, 0, 0.01)))          
-      var.cat.fc[ , eval.dexes] <- apply(s2dverification::InsertDim( 
-        s2dverification::InsertDim(optim.tmp$par, lendim = amt.cat, posdim = 2),
-        lendim =  amt.sdate.ev, posdim = 3) *
+      init.par <- optim.tmp$par * (1 - abs(rnorm(amt.coeff, 0, 0.01)))
+      var.cat.fc[ , eval.dexes] <- apply(suppressWarnings(InsertDim( 
+        InsertDim(optim.tmp$par, lendim = amt.cat, posdim = 2),
+        lendim =  amt.sdate.ev, posdim = 3)) *
         freq.per.mdl.ev[ , , , drop = FALSE], c(2,3), sum, na.rm = TRUE)
     } else if (cat.method == "comb") {
       freq.per.mdl.ev <- .calc.freq.per.mdl(cat.fc = cat.fc.ev, mdl.feat = mdl.feat, amt.cat = amt.cat)
@@ -363,9 +389,12 @@ comb.dims <- function(arr.in, dims.to.combine){
 
 .funct.optim.grad <- function(par, freq.per.mdl.at.obs){
   amt.model <- dim(freq.per.mdl.at.obs)[1]
-  return(-apply(freq.per.mdl.at.obs/s2dverification::InsertDim(drop(par %*% 
-freq.per.mdl.at.obs),
-   lendim = amt.model, posdim = 1), c(1), mean, na.rm = TRUE))
+  preprocess <- drop(par %*% freq.per.mdl.at.obs)
+  if (is.null(dim(preprocess))) {
+    dim(preprocess) <- c(dim = length(preprocess))
+  }
+  return(-apply(freq.per.mdl.at.obs/suppressWarnings(InsertDim(preprocess,
+          lendim = as.numeric(amt.model), posdim = 1)), c(1), mean, na.rm = TRUE))
 }
 
 .calc.freq.per.mdl.at.obs <- function(cat.obs, cat.fc, amt.cat, mdl.feat){
@@ -374,7 +403,7 @@ freq.per.mdl.at.obs),
   amt.mdl <- mdl.feat$amt.mdl
   mdl.msk.tmp <- mdl.feat$mdl.msk
   amt.coeff <- amt.mdl + 1
-  msk.fc.obs <- (cat.fc == s2dverification::InsertDim(cat.obs, posdim = 1, lendim = amt.mbr))
+  msk.fc.obs <- (cat.fc == InsertDim(cat.obs, posdim = 1, lendim = amt.mbr))
   freq.per.mdl.at.obs <- array(NA, c(amt.coeff, amt.sdate))
   for (i.mdl in seq(1, amt.mdl)){
     freq.per.mdl.at.obs[i.mdl, ] <- apply(msk.fc.obs[mdl.msk.tmp[i.mdl, ], , drop = FALSE],
diff --git a/R/CST_RegimesAssign.R b/R/CST_RegimesAssign.R
index 72ae69f1e16ce301f194794d55ad5f9e55251201..8985d3355bb5c048aab525bfb8a9653ec1f182cb 100644
--- a/R/CST_RegimesAssign.R
+++ b/R/CST_RegimesAssign.R
@@ -27,7 +27,8 @@
 #'         that accounts for the serial dependence of the data with the same structure as Composite.)(only when composite = 'TRUE'),
 #'         \code{$cluster} (array with the same dimensions as data (except latitude and longitude which are removed) indicating the ref_maps to which each point is allocated.) ,
 #'         \code{$frequency} (A vector of integers (from k=1,...k n reference maps) indicating the percentage of assignations corresponding to each map.),
-#'@importFrom s2dverification ACC Mean1Dim InsertDim
+#'@importFrom s2dverification ACC Mean1Dim
+#'@importFrom s2dv InsertDim
 #'@import multiApply
 #'@examples
 #'\dontrun{
@@ -102,7 +103,8 @@ CST_RegimesAssign <- function(data,  ref_maps,
 #'         \code{$cluster} (array with the same dimensions as data (except latitude and longitude which are removed) indicating the ref_maps to which each point is allocated.) ,
 #'         \code{$frequency} (A vector of integers (from k = 1, ... k n reference maps) indicating the percentage of assignations corresponding to each map.),
 #'  
-#'@importFrom s2dverification ACC Mean1Dim Eno InsertDim
+#'@importFrom s2dverification ACC Mean1Dim Eno
+#'@importFrom s2dv InsertDim
 #'@import multiApply   
 #'@examples 
 #'\dontrun{
@@ -190,9 +192,9 @@ RegimesAssign <- function(data, ref_maps, lat, method = "distance", composite =
     
     if (any(is.na(index))) {
       recon <-list(
-          composite = s2dverification::InsertDim(array(NA, dim = c(dim(dataComp)[-postime])), 
+          composite = InsertDim(array(NA, dim = c(dim(dataComp)[-postime])), 
                                 postime, dim(ref_maps)['composite.cluster']),
-          pvalue = s2dverification::InsertDim(array(NA, dim = c(dim(dataComp)[-postime])), 
+          pvalue = InsertDim(array(NA, dim = c(dim(dataComp)[-postime])), 
                              postime, dim(ref_maps)['composite.cluster']))
     } else {
       if (memb) {
@@ -255,7 +257,7 @@ RegimesAssign <- function(data, ref_maps, lat, method = "distance", composite =
                     which(names(dim(target)) == 'lon')))
   
   # weights are defined
-  latWeights <- s2dverification::InsertDim(sqrt(cos(lat * pi / 180)), 2, dim(ref)[3])
+  latWeights <- InsertDim(sqrt(cos(lat * pi / 180)), 2, dim(ref)[3])
   
   
   rmsdiff <- function(x, y) {
@@ -278,12 +280,8 @@ RegimesAssign <- function(data, ref_maps, lat, method = "distance", composite =
     corr <- rep(NA, nclust)
     for (i in 1:nclust) {
       corr[i] <-
-        ACC(s2dverification::InsertDim(InsertDim(
-          s2dverification::InsertDim(ref[i, , ] * latWeights, 1, 1), 2, 1
-        ), 3, 1),
-        s2dverification::InsertDim(InsertDim(
-          s2dverification::InsertDim(target * latWeights, 1, 1), 2, 1
-        ), 3, 1))$ACC[2]
+        ACC(InsertDim(InsertDim(InsertDim(ref[i, , ] * latWeights, 1, 1), 2, 1), 3, 1),
+          InsertDim(InsertDim(InsertDim(target * latWeights, 1, 1), 2, 1), 3, 1))$ACC[2]
     }
     assign <- which(corr == max(corr))
   }
diff --git a/R/CST_SaveExp.R b/R/CST_SaveExp.R
index 4ae560bff5757c03a91a3e251a9442dca13ed2fe..ac377101538d846e2e154dbb7f0f9a2e710e231e 100644
--- a/R/CST_SaveExp.R
+++ b/R/CST_SaveExp.R
@@ -17,8 +17,7 @@
 #'@seealso \code{\link{CST_Load}}, \code{\link{as.s2dv_cube}} and \code{\link{s2dv_cube}}
 #'
 #'@import ncdf4
-#'@importFrom s2dv Reorder
-#'@importFrom s2dverification InsertDim
+#'@importFrom s2dv Reorder InsertDim
 #'@import multiApply
 #'
 #'@examples
@@ -79,8 +78,7 @@ CST_SaveExp <- function(data, destination = "./CST_Data") {
 #' The path will be created with the name of the variable and each Datasets.
 #' 
 #'@import ncdf4
-#'@importFrom s2dv Reorder
-#'@importFrom s2dverification InsertDim
+#'@importFrom s2dv Reorder InsertDim
 #'@import multiApply
 #'
 #'@examples
@@ -144,7 +142,7 @@ SaveExp <- function(data, lon, lat, Dataset, var_name, units, startdates, Dates,
   if (length(dataset_pos) == 0) {
     warning("Element 'data' in parameter 'data' hasn't 'dataset' dimension. ",
             "All data is stored in the same 'dataset' folder.")
-    data$data <- s2dverification::InsertDim(data, posdim = 1, lendim = 1)
+    data$data <- InsertDim(data, posdim = 1, lendim = 1)
     names(dim(data))[1] <- "dataset"
     dimname <- c("dataset", dimname)
     dataset_pos = 1
diff --git a/R/PlotForecastPDF.R b/R/PlotForecastPDF.R
index 8c4b1068b9e96893018f29c378c04ec30fdae71b..4cf6de6c82c370de6a370326f78873b3416edf2c 100644
--- a/R/PlotForecastPDF.R
+++ b/R/PlotForecastPDF.R
@@ -23,7 +23,7 @@
 #'@importFrom reshape2 melt
 #'@importFrom plyr .
 #'@importFrom plyr dlply
-#'@importFrom s2dverification InsertDim
+#'@importFrom s2dv InsertDim
 #'@examples
 #'fcsts <- data.frame(fcst1 = rnorm(10), fcst2 = rnorm(10, 0.5, 1.2), 
 #'                    fcst3 = rnorm(10, -0.5, 0.9))
@@ -108,7 +108,7 @@ PlotForecastPDF <- function(fcst, tercile.limits, extreme.limits = NULL, obs = N
         if (length(tercile.limits) != 2) {
             stop("Provide two tercile limits")
         } 
-        tercile.limits <- s2dverification::InsertDim(tercile.limits, 1, npanels)
+        tercile.limits <- InsertDim(tercile.limits, 1, npanels)
     } else if (is.array(tercile.limits)) {
         if (length(dim(tercile.limits)) == 2) {
             if (dim(tercile.limits)[2] != 2) {
@@ -135,7 +135,7 @@ PlotForecastPDF <- function(fcst, tercile.limits, extreme.limits = NULL, obs = N
             if (length(extreme.limits) != 2) {
                 stop("Provide two extreme limits")
             } 
-            extreme.limits <- s2dverification::InsertDim(extreme.limits, 1, npanels)
+            extreme.limits <- InsertDim(extreme.limits, 1, npanels)
         } else if (is.array(extreme.limits)) {
             if (length(dim(extreme.limits)) == 2) {
                 if (dim(extreme.limits)[2] != 2) {
diff --git a/R/PlotTriangles4Categories.R b/R/PlotTriangles4Categories.R
index 14d578a54b912c6a358cb1ac913ce63901f8c2dc..9e3d995df521f9b0a9f9b8a8a31a8df01272c1ed 100644
--- a/R/PlotTriangles4Categories.R
+++ b/R/PlotTriangles4Categories.R
@@ -40,8 +40,10 @@
 #'@param cex_leg a number to indicate the increase/reductuion of the lab_legend used 
 #'  to represent sig_data.
 #'@param col_leg color of the legend (triangles).
+#'@param cex_axis a number to indicate the increase/reduction of the axis labels.
 #'@param fileout A string of full directory path and file name indicating where 
 #'  to save the plot. If not specified (default), a graphics device will pop up. 
+#'@param mar A numerical vector of the form c(bottom, left, top, right) which gives the number of lines of margin to be specified on the four sides of the plot. 
 #'@param size_units A string indicating the units of the size of the device 
 #'  (file or window) to plot in. Set 'px' as default. See ?Devices and the 
 #'  creator function of the corresponding device.
@@ -82,6 +84,7 @@ PlotTriangles4Categories <- function(data, brks = NULL, cols = NULL,
                                      ylabels = NULL, ytitle = NULL,
                                      legend = TRUE, lab_legend = NULL,
                                      cex_leg = 1, col_leg = 'black',
+                                     cex_axis = 1.5, mar = c(5, 4, 0, 0),
                                      fileout = NULL, size_units = 'px',
                                      res = 100, figure.width = 1, ...) {
   # Checking the dimensions
@@ -101,7 +104,9 @@ PlotTriangles4Categories <- function(data, brks = NULL, cols = NULL,
       stop("Parameter 'data' should contain 'dimx', 'dimy' and 'dimcat' dimension names. ")
     }
   }
-  
+  if (!is.vector(mar) & length(mar) != 4) {
+    stop("Parameter 'mar' must be a vector of length 4.")
+  }
   if (!is.null(sig_data)) {
     if (!is.logical(sig_data)) {
       stop("Parameter 'sig_data' array must be logical.")}
@@ -180,7 +185,7 @@ PlotTriangles4Categories <- function(data, brks = NULL, cols = NULL,
   if(legend){
     layout(matrix(c(1, 2, 1, 3), 2, 2, byrow = T),
            widths = c(10, 3.4), heights = c(10, 3.5))
-    par(oma = c(1, 1, 1, 1), mar = c(5, 4, 0, 0))
+    par(oma = c(1, 1, 1, 1), mar = mar)
     if(is.null(lab_legend)) {
       lab_legend = 1:ncat
     }
@@ -203,10 +208,10 @@ PlotTriangles4Categories <- function(data, brks = NULL, cols = NULL,
   }
   
   if (xlab){
-    axis(1, at =(1:ncol) - 0.5, las = 2, labels = xlabels, cex.axis = 1.5)
+    axis(1, at =(1:ncol) - 0.5, las = 2, labels = xlabels, cex.axis = cex_axis)
   }
   if (ylab){
-    axis(2, at = (1:nrow) - 0.5, las = 2, labels = ylabels, cex.axis = 1.5)
+    axis(2, at = (1:nrow) - 0.5, las = 2, labels = ylabels, cex.axis = cex_axis)
   }
   
   
diff --git a/man/CST_MultiMetric.Rd b/man/CST_MultiMetric.Rd
index e8e047dc6f6fba4e4dcaf6f30af4be5e3e0f4636..dc2f756623acb51ff90e108d1bade669a6a98201 100644
--- a/man/CST_MultiMetric.Rd
+++ b/man/CST_MultiMetric.Rd
@@ -30,7 +30,7 @@ CST_MultiMetric(
 \item{sdate_dim}{name of the start date dimension or a dimension name identifiying the different forecast. It can be NULL, the default value is 'sdate'.}
 }
 \value{
-an object of class \code{s2dv_cube} containing the statistics of the selected metric in the element \code{$data} which is a list of arrays: for the metric requested and others for statistics about its signeificance. The arrays have two dataset dimensions equal to the 'dataset' dimension in the \code{exp$data} and \code{obs$data} inputs. If \code{multimodel} is TRUE, the greatest position in the first dimension correspons to the Multi-Model Mean.
+an object of class \code{s2dv_cube} containing the statistics of the selected metric in the element \code{$data} which is a list of arrays: for the metric requested and others for statistics about its signeificance. The arrays have two dataset dimensions equal to the 'dataset' dimension in the \code{exp$data} and \code{obs$data} inputs. If \code{multimodel} is TRUE, the first position in the first 'nexp' dimension correspons to the Multi-Model Mean.
 }
 \description{
 This function calculates correlation (Anomaly Correlation Coefficient; ACC), root mean square error (RMS) and the root mean square error skill score (RMSSS) of individual anomaly models and multi-models mean (if desired) with the observations.
diff --git a/man/CategoricalEnsCombination.Rd b/man/CategoricalEnsCombination.Rd
new file mode 100644
index 0000000000000000000000000000000000000000..2f5ad14d69ba9db0ea970da8bc340ec73f66e2a1
--- /dev/null
+++ b/man/CategoricalEnsCombination.Rd
@@ -0,0 +1,43 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/CST_CategoricalEnsCombination.R
+\name{CategoricalEnsCombination}
+\alias{CategoricalEnsCombination}
+\title{Make categorical forecast based on a multi-model forecast with potential for calibrate}
+\usage{
+CategoricalEnsCombination(fc, obs, cat.method, eval.method, amt.cat, ...)
+}
+\arguments{
+\item{fc}{a multi-dimensional array with named dimensions containing the seasonal forecast experiment data in the element named \code{$data}. The amount of forecasting models is equal to the size of the \code{dataset} dimension of the data array. The amount of members per model may be different. The  size of the \code{member} dimension of the data array is equal to the maximum of the ensemble members among the models. Models with smaller ensemble sizes have residual indices of \code{member} dimension in the data array filled with NA values.}
+
+\item{obs}{a multidimensional array with named dimensions containing the observed data in the element named \code{$data}.}
+
+\item{cat.method}{method used to produce the categorical forecast, can be either \code{pool}, \code{comb}, \code{mmw} or \code{obs}. The method pool assumes equal weight for all ensemble members while the method comb assumes equal weight for each model. The weighting method is descirbed in Rajagopalan et al. (2002), Robertson et al. (2004) and Van Schaeybroeck and Vannitsem (2019). Finally, the \code{obs} method classifies the observations into the different categories and therefore contains only 0 and 1 values.}
+
+\item{eval.method}{is the sampling method used, can be either \code{"in-sample"} or \code{"leave-one-out"}. Default value is the \code{"leave-one-out"} cross validation.}
+
+\item{amt.cat}{is the amount of categories. Equally-sized quantiles will be calculated based on the amount of categories.}
+
+\item{...}{other parameters to be passed on to the calibration procedure.}
+}
+\value{
+an array containing the categorical forecasts in the element called \code{$data}. The first two dimensions of the returned object are named dataset and member and are both of size one. An additional dimension named category is introduced and is of size amt.cat.
+}
+\description{
+This function converts a multi-model ensemble forecast 
+into a categorical forecast by giving the probability
+for each category. Different methods are available to combine 
+the different ensemble forecasting models into 
+probabilistic categorical forecasts. 
+
+See details in ?CST_CategoricalEnsCombination
+}
+\references{
+Rajagopalan, B., Lall, U., & Zebiak, S. E. (2002). Categorical climate forecasts through regularization and optimal combination of multiple GCM ensembles. Monthly Weather Review, 130(7), 1792-1811.
+
+Robertson, A. W., Lall, U., Zebiak, S. E., & Goddard, L. (2004). Improved combination of multiple atmospheric GCM ensembles for seasonal prediction. Monthly Weather Review, 132(12), 2732-2744.
+
+Van Schaeybroeck, B., & Vannitsem, S. (2019). Postprocessing of Long-Range Forecasts. In Statistical Postprocessing of Ensemble Forecasts (pp. 267-290).
+}
+\author{
+Bert Van Schaeybroeck, \email{bertvs@meteo.be}
+}
diff --git a/man/PlotTriangles4Categories.Rd b/man/PlotTriangles4Categories.Rd
index 6e95df38351cb9f4c7bfcfa749e36af265f0d160..abd3085eff6da0742189afb1ce5f3c7d16d00ff5 100644
--- a/man/PlotTriangles4Categories.Rd
+++ b/man/PlotTriangles4Categories.Rd
@@ -23,6 +23,8 @@ PlotTriangles4Categories(
   lab_legend = NULL,
   cex_leg = 1,
   col_leg = "black",
+  cex_axis = 1.5,
+  mar = c(5, 4, 0, 0),
   fileout = NULL,
   size_units = "px",
   res = 100,
@@ -84,6 +86,10 @@ to represent sig_data.}
 
 \item{col_leg}{color of the legend (triangles).}
 
+\item{cex_axis}{a number to indicate the increase/reduction of the axis labels.}
+
+\item{mar}{A numerical vector of the form c(bottom, left, top, right) which gives the number of lines of margin to be specified on the four sides of the plot.}
+
 \item{fileout}{A string of full directory path and file name indicating where 
 to save the plot. If not specified (default), a graphics device will pop up.}