diff --git a/.github/workflows/test-package.yml b/.github/workflows/test-package.yml
new file mode 100644
index 0000000..a65e9af
--- /dev/null
+++ b/.github/workflows/test-package.yml
@@ -0,0 +1,82 @@
+on:
+  push:
+    branches:
+      - main
+      - master
+      - develop
+  pull_request:
+    branches:
+      - main
+      - master
+      - develop
+
+name: R-CMD-check
+
+jobs:
+  R-CMD-check:
+    runs-on: ${{ matrix.config.os }}
+    name: ${{ matrix.config.os }} (${{ matrix.config.r }})
+
+    strategy:
+      fail-fast: false
+      matrix:
+        config:
+          - {os: windows-latest, r: 'release'}
+          - {os: macOS-latest, r: 'release'}
+          - {os: ubuntu-20.04, r: 'release', rspm: "https://packagemanager.rstudio.com/cran/__linux__/focal/latest"}
+          - {os: ubuntu-20.04, r: 'devel', rspm: "https://packagemanager.rstudio.com/cran/__linux__/focal/latest"}
+
+    env:
+      R_REMOTES_NO_ERRORS_FROM_WARNINGS: true
+      RSPM: ${{ matrix.config.rspm }}
+      GITHUB_PAT: ${{ secrets.GITHUB_TOKEN }}
+
+    steps:
+      - uses: actions/checkout@v2
+      - uses: r-lib/actions/setup-r@v1
+        with:
+          r-version: ${{ matrix.config.r }}
+
+      - uses: r-lib/actions/setup-pandoc@v1
+
+      - name: Query dependencies
+        run: |
+          install.packages('remotes')
+          saveRDS(remotes::dev_package_deps(dependencies = TRUE), ".github/depends.Rds", version = 2)
+          writeLines(sprintf("R-%i.%i", getRversion()$major, getRversion()$minor), ".github/R-version")
+        shell: Rscript {0}
+
+      - name: Cache R packages
+        if: runner.os != 'Windows'
+        uses: actions/cache@v2
+        with:
+          path: ${{ env.R_LIBS_USER }}
+          key: ${{ runner.os }}-${{ hashFiles('.github/R-version') }}-1-${{ hashFiles('.github/depends.Rds') }}
+          restore-keys: ${{ runner.os }}-${{ hashFiles('.github/R-version') }}-1-
+
+      - name: Install system dependencies
+        if: runner.os == 'Linux'
+        run: |
+          while read -r cmd
+          do
+            eval sudo $cmd
+          done < <(Rscript -e 'writeLines(remotes::system_requirements("ubuntu", "20.04"))')
+
+      - name: Install dependencies
+        run: |
+          remotes::install_deps(dependencies = TRUE)
+          remotes::install_cran("rcmdcheck")
+        shell: Rscript {0}
+
+      - name: Check
+        env:
+          _R_CHECK_CRAN_INCOMING_REMOTE_: false
+        run: rcmdcheck::rcmdcheck(args = c("--no-manual", "--as-cran"), error_on = "warning", check_dir = "check")
+        shell: Rscript {0}
+
+      - name: Upload check results
+        if: failure()
+        uses: actions/upload-artifact@main
+        with:
+          name: ${{ runner.os }}-r${{ matrix.config.r }}-results
+          path: check
\ No newline at end of file
diff --git a/.gitignore b/.gitignore
index 1b8dad1..7cb37b6 100755
--- a/.gitignore
+++ b/.gitignore
@@ -6,6 +6,7 @@
 .DS_Store
 .lintr
 .vscode/*
+LICENSE
 
 tests/plots/*.pdf
 tests/plots/*.txt
diff --git a/.travis.yml b/.travis.yml
deleted file mode 100644
index 1ef7cc5..0000000
--- a/.travis.yml
+++ /dev/null
@@ -1,15 +0,0 @@
-language: r
-
-r:
-  - release
-  - devel
-
-r_packages:
-  - covr
-
-after_success:
-  - travis_wait 20 Rscript -e 'library(covr); codecov()'
-
-script:
-  - R CMD build .
-  - travis_wait 20 R CMD check --as-cran *.tar.gz
\ No newline at end of file
diff --git a/DESCRIPTION b/DESCRIPTION
index 69a9f0d..a058af2 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -1,8 +1,8 @@
 Package: mdatools
 Title: Multivariate Data Analysis for Chemometrics
-Version: 0.11.2
-Date: 2020-10-22
-Author: Sergey Kucheryavskiy
+Version: 0.11.3
+Date: 2021-01-21
+Author: Sergey Kucheryavskiy (<https://orcid.org/0000-0002-3145-7244>)
 Maintainer: Sergey Kucheryavskiy <svkucheryavski@gmail.com>
 Description: Projection based methods for preprocessing,
     exploring and analysis of multivariate data used in chemometrics.
diff --git a/LICENSE.md b/LICENSE.md
new file mode 100644
index 0000000..7cf9f2d
--- /dev/null
+++ b/LICENSE.md
@@ -0,0 +1,20 @@
+Copyright (c) 2016-2021 Sergey Kucheryavskiy
+
+Permission is hereby granted, free of charge, to any person obtaining
+a copy of this software and associated documentation files (the
+"Software"), to deal in the Software without restriction, including
+without limitation the rights to use, copy, modify, merge, publish,
+distribute, sublicense, and/or sell copies of the Software, and to
+permit persons to whom the Software is furnished to do so, subject to
+the following conditions:
+
+The above copyright notice and this permission notice shall be
+included in all copies or substantial portions of the Software.
+
+THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
+EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
+MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
+NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE
+LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION
+OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION
+WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
\ No newline at end of file
diff --git a/NAMESPACE b/NAMESPACE
index 8e4bec5..808c243 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -172,6 +172,7 @@ export(ddmoments.param)
 export(ddrobust.param)
 export(ellipse)
 export(employ)
+export(eye)
 export(fprintf)
 export(getCalibrationData)
 export(getConfusionMatrix)
@@ -233,6 +234,7 @@ export(pca)
 export(pca.mvreplace)
 export(pca.run)
 export(pcares)
+export(pcv)
 export(pinv)
 export(plotBars)
 export(plotBiplot)
@@ -297,6 +299,7 @@ export(prep.alsbasecorr)
 export(prep.autoscale)
 export(prep.msc)
 export(prep.norm)
+export(prep.ref2km)
 export(prep.savgol)
 export(prep.snv)
 export(prepCalData)
diff --git a/NEWS.md b/NEWS.md
index 175db33..ec834ba 100755
--- a/NEWS.md
+++ b/NEWS.md
@@ -1,3 +1,19 @@
+v.0.11.3
+========
+
+* added [Procrustes Cross-Validation method](https://doi.org/10.1021/acs.analchem.0c02175), `pcv()` (it is also available as a [separate project](https://github.com/svkucheryavski/pcv)).
+
+* added Kubelka-Munk transformation for diffuse reflectance spectra (`prep.ref2km()`).
+
+* fixed bug [#94](https://github.com/svkucheryavski/mdatools/issues/94) which caused wrong limits in PCA distance plot when outliers are present but excluded.
+
+* fixed bug [#95](https://github.com/svkucheryavski/mdatools/issues/95) which lead to issues when PLS regression methods (e.g. `plotRMSE()`) are used for PLS-DA model object.
+
+* added additional check that parameter `cgroup` for plotting functions is provided as a vector or as a factor to avoid confusion.
+
+* added link to [YouTube channel](https://www.youtube.com/channel/UCox0H4utfMq4FIu2kymuyTA) with Chemometric course based on *mdatools* package.
+
+
 v.0.11.2
 ========
 
diff --git a/R/ldecomp.R b/R/ldecomp.R
index f84cae5..d87c51b 100755
--- a/R/ldecomp.R
+++ b/R/ldecomp.R
@@ -964,6 +964,8 @@ ldecomp.getLimitsCoordinates <- function(Qlim, T2lim, ncomp, norm, log,
 #' logical, show or not legend on the plot (if more than one result object)
 #' @param legend.position
 #' if legend must be shown, where it should be
+#' @param show.excluded
+#' logical, show or hide rows marked as excluded (attribute `exclrows`).
 #' @param ...
 #' other plot parameters (see \code{mdaplotg} for details)
 #'
@@ -988,12 +990,18 @@ ldecomp.getLimitsCoordinates <- function(Qlim, T2lim, ncomp, norm, log,
 ldecomp.plotResiduals <- function(res, Qlim, T2lim, ncomp, log = FALSE, norm = FALSE,
    cgroup = NULL, xlim = NULL, ylim = NULL, show.limits = c(TRUE, TRUE),
    lim.col = c("darkgray", "darkgray"), lim.lwd = c(1, 1), lim.lty = c(2, 3),
-   show.legend = TRUE, legend.position = "topright", ...) {
+   show.legend = TRUE, legend.position = "topright", show.excluded = FALSE, ...) {
 
-   getPlotLim <- function(lim, pd, ld, dim, show.limits) {
+   # return column with values either with or without excluded outliers
+   getValues <- function(x, dim) {
+      return(if (show.excluded) x[, dim] else mda.purgeRows(x)[, dim])
+   }
+
+   # compute limits fo axis depending on values and position of critical limits
+   getPlotLim <- function(lim, pd, ld, dim) {
       if (!is.null(lim) || all(!show.limits)) return(lim)
-      limits <- if (show.limits[[2]]) ld$outliers else ld$extremes
-      return(c(0, max(sapply(pd, function(x) max(x[, dim])), limits[, dim])) * 1.05)
+      limits <- if (show.limits[[2]]) max(ld$outliers[, dim]) else max(ld$extremes[, dim])
+      return( c(0, max(sapply(pd, function(x) { max(c(getValues(x, dim), limits)) * 1.05}))) )
    }
 
    # check that show.limits is logical
@@ -1015,15 +1023,16 @@ ldecomp.plotResiduals <- function(res, Qlim, T2lim, ncomp, log = FALSE, norm = F
 
    # get coordinates for critical limits
    lim_data <- ldecomp.getLimitsCoordinates(Qlim, T2lim, ncomp = ncomp, norm = norm, log = log)
-   xlim <- getPlotLim(xlim, plot_data, lim_data, 1, show.limits)
-   ylim <- getPlotLim(ylim, plot_data, lim_data, 2, show.limits)
+   xlim <- getPlotLim(xlim, plot_data, lim_data, 1)
+   ylim <- getPlotLim(ylim, plot_data, lim_data, 2)
 
    # make plot
    if (length(plot_data) == 1) {
-      mdaplot(plot_data[[1]], type = "p", xlim = xlim, ylim = ylim, cgroup = cgroup, ...)
+      mdaplot(plot_data[[1]], type = "p", xlim = xlim, ylim = ylim, cgroup = cgroup,
+         show.excluded = show.excluded, ...)
    } else {
       mdaplotg(plot_data, type = "p", xlim = xlim, ylim = ylim, show.legend = show.legend,
-         legend.position = legend.position, ...)
+         show.excluded = show.excluded, legend.position = legend.position, ...)
    }
 
    # show critical limits
diff --git a/R/mdaplot.R b/R/mdaplot.R
index ed7e6be..b11cb17 100755
--- a/R/mdaplot.R
+++ b/R/mdaplot.R
@@ -273,6 +273,10 @@ mdaplot.getColors <- function(ngroups = NULL, cgroup = NULL, colmap = "default",
       return(mdaplot.prepareColors(palette, ngroups, opacity))
    }
 
+   if (!is.null(dim(cgroup))) {
+      stop("Parameter 'cgroup' should be a vector of values or a factor.")
+   }
+
    # if cgroup is factor return vector with corresponding values
    if (is.factor(cgroup)) {
       ngroups <- length(attr(cgroup, "levels"))
diff --git a/R/pcv.R b/R/pcv.R
new file mode 100644
index 0000000..96c33ab
--- /dev/null
+++ b/R/pcv.R
@@ -0,0 +1,177 @@
+#' Compute matrix with pseudo-validation set
+#'
+#' @param x
+#' matrix with calibration set (IxJ)
+#' @param ncomp
+#' number of components for PCA decomposition
+#' @param nseg
+#' number of segments in cross-validation
+#' @param scale
+#' logical, standardize columns of X prior to decompositon or not
+#'
+#' @description
+#' The method computes pseudo-validation matrix Xpv, based on PCA decomposition of calibration set X
+#' and systematic (venetian blinds) cross-validation. It is assumed that data rows are ordered
+#' correctly, so systematic cross-validation can be applied.
+#'
+#' All details can be found in [1]
+#'
+#' @references
+#' 1. Kucheryavskiy, S., Zhilin, S., Rodionova, O., & Pomerantsev, A.
+#' Procrustes Cross-Validation—A Bridge between Cross-Validation and Independent Validation Sets.
+#' Analytical Chemistry, 92 (17), 2020. pp.11842–11850. DOI: 10.1021/acs.analchem.0c02175
+#'
+#' @return
+#' Pseudo-validation matrix (IxJ)
+#'
+#' @export
+pcv <- function(x, ncomp = min(round(nrow(x)/nseg) - 1, col(x), 20), nseg = 4, scale = FALSE) {
+
+   # keep names if any
+   attrs <- attributes(x)
+
+   # convert to matrix if necessary
+   x <- mda.df2mat(x)
+
+   mx <- apply(x, 2, mean)
+   sx <- if (scale) apply(x, 2, sd) else rep(1, ncol(x))
+
+   # autoscale the calibration set
+   x <- scale(x, center = mx, scale = sx)
+
+   # create a global model
+   P <- svd(x)$v[, seq_len(ncomp), drop = FALSE]
+
+   # create matrix with indices for cross-validation, so
+   # each column is number of rows to be taken as local validation set
+   # in corresponding segment
+   ind <- matrix(seq_len(nrow(x)), ncol = nseg, byrow = TRUE)
+
+   # prepare empty matrix for pseudo-validation set
+   x.pv <- matrix(0, nrow(x), ncol(x))
+   a <- NULL
+
+   # cv-loop
+   for (k in seq_len(nseg)) {
+
+      # split data to calibration and validation
+      x.c <- x[-ind[, k], , drop = FALSE]
+      x.k <- x[ ind[, k], , drop = FALSE]
+
+      # get loadings for local model and rotation matrix between global and local models
+      P.k <- svd(x.c, nv = ncomp)$v[, seq_len(ncomp), drop = FALSE]
+
+      # correct direction of loadings for local model
+      a <- acos(colSums(P * P.k)) < pi / 2
+      P.k <- P.k %*% diag(a * 2 - 1, ncol(P), ncol(P))
+
+      # get rotation matrix between the PC spaces
+      R <- getR(P.k, P)
+
+      # rotate the local validation set and save as a part of Xpv
+      x.pv[ind[, k], ] <- tcrossprod(x.k, R)
+   }
+
+   # uscenter and unscale the data
+   x.pv <- sweep(x.pv, 2, sx, "*")
+   x.pv <- sweep(x.pv, 2, mx, "+")
+
+   attributes(x.pv) <- attrs
+   return(x.pv)
+}
+
+#' Creates rotation matrix to map a set vectors \code{base1} to a set of vectors \code{base2}.
+#'
+#' @param base1
+#' Matrix (JxA) with A orthonormal vectors as columns to be rotated (A <= J)
+#' @param base2
+#' Matrix (JxA) with A orthonormal vectors as columns, \code{base1} should be aligned with
+#'
+#' @description
+#' In both sets vectors should be orthonormal.
+#'
+#' @return
+#' Rotation matrix (JxJ)
+getR <- function(base1, base2) {
+   base1 <- as.matrix(base1);
+   base2 <- as.matrix(base2);
+
+   R1 <- rotationMatrixToX1(base1[, 1])
+   R2 <- rotationMatrixToX1(base2[, 1])
+
+   if (ncol(base1) == 1) {
+      R <- t(R2) %*% R1
+   } else {
+      # Compute bases rotated to match their first vectors to [1 0 0 ... 0]'
+      base1_r <- as.matrix(R1 %*% base1)
+      base2_r <- as.matrix(R2 %*% base2)
+
+      # Get bases of subspaces of dimension n-1 (forget x1)
+      nr <- nrow(base1_r) # equal to nrow(base2_r)
+      nc <- ncol(base1_r) # equal to ncol(base2_r)
+      base1_rs <- base1_r[2:nr, 2:nc]
+      base2_rs <- base2_r[2:nr, 2:nc]
+
+      # Recursevely compute rotation matrix to map subspaces
+      Rs <- getR(base1_rs, base2_rs)
+
+      # Construct rotation matrix of the whole space (recall x1)
+      M <- eye(nr)
+      M[2:nr, 2:nr] <- Rs
+
+      R <- crossprod(R2, (M %*% R1))
+   }
+
+   return(R);
+}
+
+#' Creates a rotation matrix to map a vector x to [1 0 0 ... 0]
+#'
+#' @param x
+#' Vector (sequence with J coordinates)
+#'
+#' @return
+#' Rotation matrix (JxJ)
+rotationMatrixToX1 <- function(x) {
+   N <- length(x)
+   R <- eye(N)
+   step <- 1
+   while (step < N) {
+      A <- eye(N)
+      n <- 1
+      while (n <= N - step) {
+         r2 <- x[n]^2 + x[n + step]^2
+         if (r2 > 0) {
+            r <- sqrt(r2)
+            pcos <- x[n] / r
+            psin <- -x[n + step] / r
+            A[n, n] <- pcos
+            A[n, n + step] <- -psin
+            A[n + step, n] <- psin
+            A[n + step, n + step] <- pcos
+         }
+         n <- n + 2 * step
+      }
+      step <- 2 * step
+      x <- A %*% x
+      R <- A %*% R
+   }
+   return(R)
+}
+
+
+#' Create the identity matrix
+#'
+#' @param n
+#' Size of the matrix
+#'
+#' @return
+#' The identity matrix (n x n)
+#'
+#' @export
+eye <- function(n) {
+   X <- matrix(0, n, n)
+   diag(X) <- 1
+   return(X)
+}
+
diff --git a/R/plotseries.R b/R/plotseries.R
index 63d691b..988361f 100644
--- a/R/plotseries.R
+++ b/R/plotseries.R
@@ -640,10 +640,13 @@ plotBars <- function(ps, col = ps$col, bwd = 0.8, border = NA, force.x.values =
    y <- ps$y_values[1, ]
 
    if (length(x) > 1) {
-      bwd_left <- c(x[seq(2, length(x))] - x[seq(1, length(x) - 1)])
-      bwd_right <- -c(x[seq(1, length(x) - 1)] - x[seq(2, length(x))])
-      bwd_left <- c(bwd_left[1], bwd_left) * bwd / 2
-      bwd_right <- c(bwd_right, bwd_right[length(bwd_right)]) * bwd / 2
+      # this gives variable width for bars and does not work well
+      #bwd_left <- c(x[seq(2, length(x))] - x[seq(1, length(x) - 1)])
+      #bwd_right <- -c(x[seq(1, length(x) - 1)] - x[seq(2, length(x))])
+      #bwd_left <- c(bwd_left[1], bwd_left) * bwd / 2
+      #bwd_right <- c(bwd_right, bwd_right[length(bwd_right)]) * bwd / 2
+      dx <- min(diff(x))
+      bwd_right <- bwd_left <- dx * bwd / 2
    } else {
       bwd_left <- bwd_right <- bwd * x / 2
    }
diff --git a/R/plsdares.R b/R/plsdares.R
index 4540fb9..1548f06 100755
--- a/R/plsdares.R
+++ b/R/plsdares.R
@@ -129,7 +129,7 @@
 #' @export
 plsdares <- function(plsres, cres) {
    obj <- c(plsres, cres)
-   class(obj) <- c("plsdares", "classres", "plsres")
+   class(obj) <- c("plsdares", "classres", "plsres", "regres")
    obj$call <- match.call()
 
    return(obj)
diff --git a/R/prep.R b/R/prep.R
index c20474b..1024792 100755
--- a/R/prep.R
+++ b/R/prep.R
@@ -26,42 +26,40 @@
 #' @export
 prep.autoscale <- function(data, center = TRUE, scale = FALSE, max.cov = 0) {
 
-   attrs <- mda.getattr(data)
-   dimnames <- dimnames(data)
+   f <- function(data, center, scale, max.cov) {
+      # define values for centering
+      if (is.logical(center) && center) center <- apply(data, 2, mean)
 
-   # define values for centering
-   if (is.logical(center) && center) center <- apply(data, 2, mean)
-
-   if (is.numeric(center) && length(center) != ncol(data)) {
-      stop("Number of values in 'center' should be the same as number of columns in 'daata'")
-   }
+      if (is.numeric(center) && length(center) != ncol(data)) {
+         stop("Number of values in 'center' should be the same as number of columns in 'daata'")
+      }
 
-   # define values for weigting
-   if (is.logical(scale) && scale) scale <- apply(data, 2, sd)
+      # define values for weigting
+      if (is.logical(scale) && scale) scale <- apply(data, 2, sd)
 
-   if (is.numeric(scale) && length(scale) != ncol(data)) {
-      stop("Number of values in 'scale' should be the same as number of columns in 'daata'")
-   }
+      if (is.numeric(scale) && length(scale) != ncol(data)) {
+         stop("Number of values in 'scale' should be the same as number of columns in 'daata'")
+      }
 
-   # compute coefficient of variation and set scale to 1 if it is below
-   # a user defined threshold
-   if (is.numeric(scale)) {
-      m <- if (is.numeric(center)) center else apply(data, 2, mean)
-      cv <- scale / abs(m) * 100
-      scale[is.nan(cv) | cv <= max.cov] <- 1
-   }
+      # compute coefficient of variation and set scale to 1 if it is below
+      # a user defined threshold
+      if (is.numeric(scale)) {
+         m <- if (is.numeric(center)) center else apply(data, 2, mean)
+         cv <- scale / abs(m) * 100
+         scale[is.nan(cv) | cv <= max.cov] <- 1
+      }
 
-   # make autoscaling and attach preprocessing attributes
-   data <- scale(data, center = center, scale = scale)
+      # make autoscaling and attach preprocessing attributes
+      data <- scale(data, center = center, scale = scale)
+      attr(data, "scaled:center") <- NULL
+      attr(data, "scaled:scale") <- NULL
+      attr(data, "prep:center") <- center
+      attr(data, "prep:scale") <- scale
 
-   data <- mda.setattr(data, attrs)
-   attr(data, "scaled:center") <- NULL
-   attr(data, "scaled:scale") <- NULL
-   attr(data, "prep:center") <- center
-   attr(data, "prep:scale") <- scale
-   dimnames(data) <- dimnames
+      return(data)
+   }
 
-   return(data)
+   return(prep.generic(data, f, center = center, scale = scale, max.cov = max.cov))
 }
 
 #' Standard Normal Variate transformation
@@ -97,14 +95,9 @@ prep.autoscale <- function(data, center = TRUE, scale = FALSE, max.cov = 0) {
 #'
 #' @export
 prep.snv <- function(data) {
-   attrs <- mda.getattr(data)
-   dimnames <- dimnames(data)
-
-   data <- t(scale(t(data), center = TRUE, scale = TRUE))
-   data <- mda.setattr(data, attrs)
-   dimnames(data) <- dimnames
 
-   return(data)
+   f <- function(data) t(scale(t(data), center = TRUE, scale = TRUE))
+   return(prep.generic(data, f))
 }
 
 #' Normalization
@@ -122,25 +115,21 @@ prep.snv <- function(data) {
 #'
 #' @export
 prep.norm <- function(data, type = "area") {
-   attrs <- mda.getattr(data)
-   dimnames <- dimnames(data)
-
-   w <- switch(
-      type,
-      "area" = apply(abs(data), 1, sum),
-      "length" = sqrt(apply(data^2, 1, sum)),
-      "sum" = apply(data, 1, sum)
-   )
-
-   if (is.null(w)) {
-      stop("Wrong value for argument 'type'.")
-   }
 
-   data <- sweep(data, 1, w, "/")
-   data <- mda.setattr(data, attrs)
-   dimnames(data) <- dimnames
+   f <- function(data, type) {
+
+      w <- switch(
+         type,
+         "area" = apply(abs(data), 1, sum),
+         "length" = sqrt(apply(data^2, 1, sum)),
+         "sum" = apply(data, 1, sum)
+      )
 
-   return(data)
+      if (is.null(w)) stop("Wrong value for argument 'type'.")
+      return(sweep(data, 1, w, "/"))
+   }
+
+   return(prep.generic(data, f, type = type))
 }
 
 #' Savytzky-Golay filter
@@ -159,45 +148,30 @@ prep.norm <- function(data, type = "area") {
 #'
 #' @export
 prep.savgol <- function(data, width = 3, porder = 1, dorder = 0) {
-   attrs <- mda.getattr(data)
-   dimnames <- dimnames(data)
-
-   if (width < 3) {
-      stop("Filter width ('width') should be equal at least to 3.")
-   }
-
-   if (width %% 2 != 1) {
-      stop("Filter width ('width') should be an odd integer number.")
-   }
-
-   if (dorder < 0 || dorder > 2) {
-      stop("Wrong value for the derivative order (should be 0, 1, or 2).")
-   }
-
-   if (porder < 0 || porder > 4) {
-      stop("Wrong value for the polynomial degree (should be integer number between 0 and 4).")
-   }
-
-   if (porder < dorder) {
-      stop("Polynomal degree ('porder') should not be smaller the derivative order ('dorder').")
-   }
 
-   nobj <- nrow(data)
-   nvar <- ncol(data)
-   pdata <- matrix(0, ncol = nvar, nrow = nobj)
+   f <- function(data, width, porder, dorder) {
+      stopifnot("Filter width ('width') should be equal at least to 3." = width > 2)
+      stopifnot("Filter width ('width') should be an odd integer number." = width %% 2 == 1)
+      stopifnot("Wrong value for the derivative order (should be 0, 1, or 2)." = dorder %in% (0:2))
+      stopifnot("Wrong value for the polynomial degree (should be integer number between 0 and 4)." = porder %in% (0:4))
+      stopifnot("Polynomal degree ('porder') should not be smaller the derivative order ('dorder')." = porder >= dorder)
+
+      nobj <- nrow(data)
+      nvar <- ncol(data)
+      pdata <- matrix(0, ncol = nvar, nrow = nobj)
+
+      for (i in seq_len(nobj)) {
+         d <- data[i, ]
+         w <- (width - 1) / 2
+         f <- pinv(outer(-w:w, 0:porder, FUN = "^"))
+         d <- convolve(d, rev(f[dorder + 1, ]), type = "o")
+         pdata[i, ] <- d[(w + 1) : (length(d) - w)]
+      }
 
-   for (i in seq_len(nobj)) {
-      d <- data[i, ]
-      w <- (width - 1) / 2
-      f <- pinv(outer(-w:w, 0:porder, FUN = "^"))
-      d <- convolve(d, rev(f[dorder + 1, ]), type = "o")
-      pdata[i, ] <- d[(w + 1) : (length(d) - w)]
+      return(pdata)
    }
 
-   pdata <- mda.setattr(pdata, attrs)
-   dimnames(pdata) <- dimnames
-
-   return(pdata)
+   return(prep.generic(data, f, width = width, porder = porder, dorder = dorder))
 }
 
 #' Multiplicative Scatter Correction transformation
@@ -233,49 +207,57 @@ prep.savgol <- function(data, width = 3, porder = 1, dorder = 0) {
 #'
 #' @export
 prep.msc <- function(spectra, mspectrum = NULL) {
-   attrs <- mda.getattr(spectra)
-   dimnames <- dimnames(spectra)
 
-   if (is.null(mspectrum)) {
-      mspectrum <- apply(spectra, 2, mean)
-   }
+   f <- function(spectra, mspectrum) {
+      if (is.null(mspectrum)) {
+         mspectrum <- apply(spectra, 2, mean)
+      }
 
-   if (!is.null(mspectrum)) {
-      dim(mspectrum) <- NULL
-   }
+      if (!is.null(mspectrum)) {
+         dim(mspectrum) <- NULL
+      }
 
-   if (length(mspectrum) != ncol(spectra)) {
-      stop("Length of 'mspectrum' should be the same as number of columns in 'spectra'.")
-   }
+      if (length(mspectrum) != ncol(spectra)) {
+         stop("Length of 'mspectrum' should be the same as number of columns in 'spectra'.")
+      }
 
-   cspectra <- matrix(0, nrow = nrow(spectra), ncol = ncol(spectra))
-   for (i in seq_len(nrow(spectra))) {
-      coef <- coef(lm(spectra[i, ] ~ mspectrum))
-      cspectra[i, ] <- (spectra[i, ] - coef[1]) / coef[2]
-   }
+      pspectra <- matrix(0, nrow = nrow(spectra), ncol = ncol(spectra))
+      for (i in seq_len(nrow(spectra))) {
+         coef <- coef(lm(spectra[i, ] ~ mspectrum))
+         pspectra[i, ] <- (spectra[i, ] - coef[1]) / coef[2]
+      }
 
-   cspectra <- mda.setattr(cspectra, attrs)
-   attr(cspectra, "mspectrum") <- mspectrum
-   dimnames(cspectra) <- dimnames
+      attr(pspectra, "mspectrum") <- mspectrum
+      return(pspectra)
+   }
 
-   return(cspectra)
+   return(prep.generic(spectra, f, mspectrum = mspectrum))
 }
 
-#' Pseudo-inverse matrix
+
+#' Kubelka-Munk transformation
 #'
 #' @description
-#' Computes pseudo-inverse matrix using SVD
+#' Applies Kubelka-Munk (km) transformation to data matrix (spectra)
 #'
-#' @param data
-#' a matrix with data values to compute inverse for
+#' @param spectra
+#' a matrix with spectra values (absolute reflectance values)
+#'
+#' @return
+#' preprocessed spectra.
+#'
+#' @details
+#' Kubelka-Munk is useful preprocessing method for diffuse reflection spectra (e.g. taken for
+#' powders or rough surface). It transforms the reflectance spectra R to K/M units as follows:
+#' (1 - R)^2 / 2R
 #'
 #' @export
-pinv <- function(data) {
-   # Calculates pseudo-inverse of data matrix
-   s <- svd(data)
-   s$v %*% diag(1 / s$d) %*% t(s$u)
-}
+prep.ref2km <- function(spectra) {
+   stopifnot("Can't use Kubelka-Munk transformation as some of the values are zeros or negative." = all(spectra > 0))
+   f <- function(x) (1 - x)^2 / (2 * x)
 
+   return(prep.generic(spectra, f))
+}
 
 #' Baseline correction using assymetric least squares
 #'
@@ -347,3 +329,37 @@ prep.alsbasecorr <- function(spectra, plambda = 5, p = 0.1, max.niter = 10) {
    return(pspectra)
 }
 
+
+#' Pseudo-inverse matrix
+#'
+#' @description
+#' Computes pseudo-inverse matrix using SVD
+#'
+#' @param data
+#' a matrix with data values to compute inverse for
+#'
+#' @export
+pinv <- function(data) {
+   # Calculates pseudo-inverse of data matrix
+   s <- svd(data)
+   s$v %*% diag(1 / s$d) %*% t(s$u)
+}
+
+#' Generic function for preprocessing
+#'
+#' @param x
+#' data matrix to be preprocessed
+#' @param f
+#' function for preprocessing
+#' @param ...
+#' arguments for the function f
+#'
+#'
+prep.generic <- function(x, f, ...) {
+   attrs <- mda.getattr(x)
+   dimnames <- dimnames(x)
+   x.p <- f(x, ...)
+   x.p <- mda.setattr(x.p, attrs)
+   dimnames(x.p) <- dimnames
+   return(x.p)
+}
\ No newline at end of file
diff --git a/R/regmodel.R b/R/regmodel.R
index a261f9f..bc183f4 100755
--- a/R/regmodel.R
+++ b/R/regmodel.R
@@ -325,7 +325,7 @@ plotPredictions.regmodel <- function(obj, ncomp = obj$ncomp.selected, ny = 1,
       stop("Wrong value for 'ncomp' parameter.")
    }
 
-   plot_data <- lapply(res, plotPredictions, ny = ny, ncomp = ncomp, show.plot = FALSE)
+   plot_data <- lapply(res, plotPredictions.regres, ny = ny, ncomp = ncomp, show.plot = FALSE)
    attr(plot_data[[1]], "name") <- sprintf("Predictions (ncomp = %d)", ncomp)
    plots <- mdaplotg(plot_data, type = "p", legend.position = legend.position, ...)
 
diff --git a/README.md b/README.md
index 0598e52..8f9ad7f 100755
--- a/README.md
+++ b/README.md
@@ -1,8 +1,8 @@
 Multivariate Data Analysis Tools
 ===========================================
-![Travis (build status)](https://img.shields.io/travis/svkucheryavski/mdatools?color=blue&style=flat-square "Travis CI build status")
+![GitHub build status](https://github.com/svkucheryavski/mdatools/workflows/R-CMD-check/badge.svg)
+![GitHub All Releases](https://img.shields.io/github/downloads/svkucheryavski/mdatools/total?color=blue&logo=Github "Downloads from GitHub")
 ![Downloads (CRAN)](https://cranlogs.r-pkg.org/badges/grand-total/mdatools?color=blue&logo=R&style=flat-square "Downloads from CRAN")
-![GitHub All Releases](https://img.shields.io/github/downloads/svkucheryavski/mdatools/total?color=blue&logo=Github&style=flat-square "Downloads from GitHub")
 [![codecov](https://codecov.io/gh/svkucheryavski/mdatools/branch/0.10.0/graph/badge.svg?style=flat-square)](https://codecov.io/gh/svkucheryavski/mdatools)
 
 <img src="https://mdatools.com//images/logo.svg" align="left" style="top: 5px 10px 5px 0;" />
@@ -11,22 +11,24 @@ Multivariate Data Analysis Tools
 
 For more details and examples read a [Bookdown tutorial](https://mdatools.com/docs/). The project website, [mdatools.com](https://mdatools.com), contains additional information about supplementary materials and tools.
 
+You can also take video-lectures from [YouTube channel](https://www.youtube.com/channel/UCox0H4utfMq4FIu2kymuyTA) devoted to introductory Chemometric course I give to master students. The lectures explain theory behind basic Chemometric methods but also show how to use them in *mdatools*.
+
 If you want to cite the package, please use the following: Sergey Kucheryavskiy, *mdatools – R package for chemometrics*, Chemometrics and Intelligent Laboratory Systems, Volume 198,
 2020 (DOI: [10.1016/j.chemolab.2020.103937](https://doi.org/10.1016/j.chemolab.2020.103937)).
 
 What is new
 -----------
 
-Latest release (0.11.2) is available both from GitHub and CRAN. You can see the full list of changes [here](NEWS.md). The Bookdown tutorial has been also updated and contains the description of new methods added in the last release.
+Latest release (0.11.3) is available both from GitHub and CRAN. You can see the full list of changes [here](NEWS.md). The Bookdown tutorial has been also updated and contains the description of new methods added in the last release.
 
 
 How to install
 --------------
 
-The package is available from CRAN by usual installing procedure. However, due to restrictions in CRAN politics regarding number of submissions (one in 3-4 month), mostly major releases will be published there (with 2-3 weeks delay after GitHub release as more thorough testing is needed). You can [download](https://github.com/svkucheryavski/mdatools/releases) a zip-file with source package and install it using the `install.packages` command, e.g. if the downloaded file is `mdatools_0.11.2.tar.gz` and it is located in a current working directory, just run the following:
+The package is available from CRAN by usual installing procedure. However, due to restrictions in CRAN politics regarding number of submissions (one in 3-4 month), mostly major releases will be published there (with 2-3 weeks delay after GitHub release as more thorough testing is needed). You can [download](https://github.com/svkucheryavski/mdatools/releases) a zip-file with source package and install it using the `install.packages` command, e.g. if the downloaded file is `mdatools_0.11.3.tar.gz` and it is located in a current working directory, just run the following:
 
 ```
-install.packages("mdatools_0.11.2.tar.gz")
+install.packages("mdatools_0.11.3.tar.gz")
 ```
 
 If you have `devtools` package installed, the following command will install the current developer version from the master branch of GitHub repository (do not forget to load the `devtools` package first):
diff --git a/man/eye.Rd b/man/eye.Rd
new file mode 100644
index 0000000..1d257c9
--- /dev/null
+++ b/man/eye.Rd
@@ -0,0 +1,17 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/pcv.R
+\name{eye}
+\alias{eye}
+\title{Create the identity matrix}
+\usage{
+eye(n)
+}
+\arguments{
+\item{n}{Size of the matrix}
+}
+\value{
+The identity matrix (n x n)
+}
+\description{
+Create the identity matrix
+}
diff --git a/man/getR.Rd b/man/getR.Rd
new file mode 100644
index 0000000..045b7e7
--- /dev/null
+++ b/man/getR.Rd
@@ -0,0 +1,19 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/pcv.R
+\name{getR}
+\alias{getR}
+\title{Creates rotation matrix to map a set vectors \code{base1} to a set of vectors \code{base2}.}
+\usage{
+getR(base1, base2)
+}
+\arguments{
+\item{base1}{Matrix (JxA) with A orthonormal vectors as columns to be rotated (A <= J)}
+
+\item{base2}{Matrix (JxA) with A orthonormal vectors as columns, \code{base1} should be aligned with}
+}
+\value{
+Rotation matrix (JxJ)
+}
+\description{
+In both sets vectors should be orthonormal.
+}
diff --git a/man/ldecomp.plotResiduals.Rd b/man/ldecomp.plotResiduals.Rd
index df384fb..1351769 100644
--- a/man/ldecomp.plotResiduals.Rd
+++ b/man/ldecomp.plotResiduals.Rd
@@ -20,6 +20,7 @@ ldecomp.plotResiduals(
   lim.lty = c(2, 3),
   show.legend = TRUE,
   legend.position = "topright",
+  show.excluded = FALSE,
   ...
 )
 }
@@ -54,6 +55,8 @@ ldecomp.plotResiduals(
 
 \item{legend.position}{if legend must be shown, where it should be}
 
+\item{show.excluded}{logical, show or hide rows marked as excluded (attribute `exclrows`).}
+
 \item{...}{other plot parameters (see \code{mdaplotg} for details)}
 }
 \description{
diff --git a/man/mdatools.Rd b/man/mdatools.Rd
index 7271620..e715882 100644
--- a/man/mdatools.Rd
+++ b/man/mdatools.Rd
@@ -6,17 +6,17 @@ Package for Multivariate Data Analysis (Chemometrics)
 }
 
 \description{
-This package contains classes and functions for most common methods used in Chemometrics. For a complete list of functions, use \code{library(help = 'mdatools')}. 
+This package contains classes and functions for most common methods used in Chemometrics. For a complete list of functions, use \code{library(help = 'mdatools')}.
 }
 
 \details{
-The porject is hosted on GitHub, there you can also find a Bookdown user tutorial (https://svkucheryavski.github.io/mdatools/) explaining most important features of the package.
+The project is hosted on GitHub (https://svkucheryavski.github.io/mdatools/), there you can also find a Bookdown user tutorial explaining most important features of the package. There is also a dedicated YouTube channel (https://www.youtube.com/channel/UCox0H4utfMq4FIu2kymuyTA) with introductory Chemometric course with examples based on mdatools functionality.
 
 Every method is represented by two classes: a model class for keeping all parameters and information about the model, and a class for keeping and visualising results of applying the model to particular data values.
 
-Every model class, e.g. \code{\link{pls}}, has all needed functionality implemented as class methods, including model calibration, validation (test set and cross-validation), visualisation of the calibration and validation results with various plots and summary statistics. 
+Every model class, e.g. \code{\link{pls}}, has all needed functionality implemented as class methods, including model calibration, validation (test set and cross-validation), visualisation of the calibration and validation results with various plots and summary statistics.
 
-So far the following modelling methods are implemented:
+So far the following modelling and validation methods are implemented:
 \tabular{ll}{
       \code{\link{pca}}, \code{\link{pcares}} \tab Principal Component Analysis (PCA).\cr
       \code{\link{pls}}, \code{\link{plsres}} \tab Partial Least Squares regression (PLS).\cr
@@ -25,6 +25,9 @@ So far the following modelling methods are implemented:
       \code{\link{plsda}}, \code{\link{plsdares}} \tab Partial Least Squares Dscriminant Analysis (PLS-DA).\cr
       \code{\link{randtest}} \tab Randomization test for PLS-regression.\cr
       \code{\link{ipls}} \tab Interval PLS variable.\cr
+      \code{\link{mcrals}} \tab Multivariate Curve Resolution with Alternating Least Squares.\cr
+      \code{\link{mcrpure}} \tab Multivariate Curve Resolution with Purity approach.\cr
+      \code{\link{pcv}} \tab Procrustes Cross Validation.\cr
 }
 
 Methods for data preprocessing:
@@ -34,9 +37,10 @@ Methods for data preprocessing:
       \code{\link{prep.snv}} \tab Standard normal variate.\cr
       \code{\link{prep.msc}} \tab Multiplicative scatter correction.\cr
       \code{\link{prep.norm}} \tab Spectra normalization.\cr
+      \code{\link{prep.alsbasecorr}} \tab Baseline correction with Asymmetric Least Squares.\cr
 }
 
-All plotting methods are based on two functions, \code{\link{mdaplot}} and \code{\link{mdaplotg}}. The functions extend the basic functionality of R plots and allow to make automatic legend and color grouping of data points or lines with colorbar legend, automatically adjust axes limits when several data groups are plotted and so on. 
+All plotting methods are based on two functions, \code{\link{mdaplot}} and \code{\link{mdaplotg}}. The functions extend the basic functionality of R plots and allow to make automatic legend and color grouping of data points or lines with colorbar legend, automatically adjust axes limits when several data groups are plotted and so on.
 }
 
 \author{
diff --git a/man/pcv.Rd b/man/pcv.Rd
new file mode 100644
index 0000000..f5ae4b6
--- /dev/null
+++ b/man/pcv.Rd
@@ -0,0 +1,37 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/pcv.R
+\name{pcv}
+\alias{pcv}
+\title{Compute matrix with pseudo-validation set}
+\usage{
+pcv(
+  x,
+  ncomp = min(round(nrow(x)/nseg) - 1, col(x), 20),
+  nseg = 4,
+  scale = FALSE
+)
+}
+\arguments{
+\item{x}{matrix with calibration set (IxJ)}
+
+\item{ncomp}{number of components for PCA decomposition}
+
+\item{nseg}{number of segments in cross-validation}
+
+\item{scale}{logical, standardize columns of X prior to decompositon or not}
+}
+\value{
+Pseudo-validation matrix (IxJ)
+}
+\description{
+The method computes pseudo-validation matrix Xpv, based on PCA decomposition of calibration set X
+and systematic (venetian blinds) cross-validation. It is assumed that data rows are ordered
+correctly, so systematic cross-validation can be applied.
+
+All details can be found in [1]
+}
+\references{
+1. Kucheryavskiy, S., Zhilin, S., Rodionova, O., & Pomerantsev, A.
+Procrustes Cross-Validation—A Bridge between Cross-Validation and Independent Validation Sets.
+Analytical Chemistry, 92 (17), 2020. pp.11842–11850. DOI: 10.1021/acs.analchem.0c02175
+}
diff --git a/man/prep.generic.Rd b/man/prep.generic.Rd
new file mode 100644
index 0000000..400ca75
--- /dev/null
+++ b/man/prep.generic.Rd
@@ -0,0 +1,18 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/prep.R
+\name{prep.generic}
+\alias{prep.generic}
+\title{Generic function for preprocessing}
+\usage{
+prep.generic(x, f, ...)
+}
+\arguments{
+\item{x}{data matrix to be preprocessed}
+
+\item{f}{function for preprocessing}
+
+\item{...}{arguments for the function f}
+}
+\description{
+Generic function for preprocessing
+}
diff --git a/man/prep.ref2km.Rd b/man/prep.ref2km.Rd
new file mode 100644
index 0000000..3b16ea7
--- /dev/null
+++ b/man/prep.ref2km.Rd
@@ -0,0 +1,22 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/prep.R
+\name{prep.ref2km}
+\alias{prep.ref2km}
+\title{Kubelka-Munk transformation}
+\usage{
+prep.ref2km(spectra)
+}
+\arguments{
+\item{spectra}{a matrix with spectra values (absolute reflectance values)}
+}
+\value{
+preprocessed spectra.
+}
+\description{
+Applies Kubelka-Munk (km) transformation to data matrix (spectra)
+}
+\details{
+Kubelka-Munk is useful preprocessing method for diffuse reflection spectra (e.g. taken for
+powders or rough surface). It transforms the reflectance spectra R to K/M units as follows:
+(1 - R)^2 / 2R
+}
diff --git a/man/rotationMatrixToX1.Rd b/man/rotationMatrixToX1.Rd
new file mode 100644
index 0000000..8b3b2bf
--- /dev/null
+++ b/man/rotationMatrixToX1.Rd
@@ -0,0 +1,17 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/pcv.R
+\name{rotationMatrixToX1}
+\alias{rotationMatrixToX1}
+\title{Creates a rotation matrix to map a vector x to [1 0 0 ... 0]}
+\usage{
+rotationMatrixToX1(x)
+}
+\arguments{
+\item{x}{Vector (sequence with J coordinates)}
+}
+\value{
+Rotation matrix (JxJ)
+}
+\description{
+Creates a rotation matrix to map a vector x to [1 0 0 ... 0]
+}
diff --git a/tests/testthat/test-ldecompplots.R b/tests/testthat/test-ldecompplots.R
index ca74947..9039ef8 100644
--- a/tests/testthat/test-ldecompplots.R
+++ b/tests/testthat/test-ldecompplots.R
@@ -195,3 +195,41 @@ test_that("variance plot works fine.", {
       plotCumVariance(obj, type = "h", col = "red", show.labels = TRUE, labels = "names")
    })
 })
+
+test_that("Bug 94: limits for residuals plot are computed correctly", {
+   set.seed(42)
+   X <- matrix(rnorm(100 * 10), 100, 10)
+   X[20, ] <- X[20, ] * 10
+
+   m1 <- pca(X, 5)
+   m2 <- pca(X, 5, exclrows = 20)
+
+   ldecomp.plotResiduals(res = list(m1$res$cal), m1$Qlim, m1$T2lim, 5, norm = TRUE)
+   p1 <- par("usr")
+
+   ldecomp.plotResiduals(res = list(m2$res$cal), m2$Qlim, m2$T2lim, 5, norm = TRUE)
+   p2 <- par("usr")
+
+   ldecomp.plotResiduals(res = list(m2$res$cal), m2$Qlim, m2$T2lim, 5, norm = TRUE, show.limits = FALSE)
+   p3 <- par("usr")
+
+   ldecomp.plotResiduals(res = list(m2$res$cal), m2$Qlim, m2$T2lim, 5, norm = TRUE, show.excluded = TRUE)
+   p4 <- par("usr")
+
+   # check that limits without outlier are small enough
+   expect_lt(p1[2], 35)
+   expect_lt(p1[4], 7)
+
+   # check that limits with hidden outlier are small enough
+   expect_lt(p2[2], 8)
+   expect_lt(p2[4], 8)
+
+   # check that limits with shown outlier are larger
+   expect_lt(p2[2], p4[2])
+   expect_lt(p2[4], p4[4])
+
+   # now limits and hidden outliers should have smallest limits
+   expect_lt(p3[2], p2[2])
+   expect_lt(p3[4], p2[4])
+
+})
\ No newline at end of file
diff --git a/tests/testthat/test-mdaplot2.R b/tests/testthat/test-mdaplot2.R
index ae5ce86..51ba08e 100644
--- a/tests/testthat/test-mdaplot2.R
+++ b/tests/testthat/test-mdaplot2.R
@@ -60,6 +60,10 @@ test_that("excluded values and color grouping work fine together", {
    expect_silent(mdaplot(pd, cgroup = d, main = "Colormap: jet", colmap = "jet"))
    expect_silent(mdaplot(pd, cgroup = d, main = "Colormap: user", colmap = usr_cmap1))
    expect_silent(mdaplot(pd, cgroup = d, main = "Colormap: user", colmap = usr_cmap2))
+
+   # if cgroup is a matrix or is a vector - raise an error
+   expect_error(mdaplot(pd, cgroup = as.matrix(d)))
+   expect_error(mdaplot(pd, cgroup = as.data.frame(d)))
 })
 
 par(mfrow = c(2, 2))
@@ -75,7 +79,7 @@ test_that("opacity parameter works well together with color grouping", {
    expect_silent(tf(type = "p", cgroup = people[, "Wine"], opacity = 0.5))
    expect_silent(tf(type = "l", cgroup = people[, "Wine"], opacity = 0.5))
    expect_silent(tf(type = "b", cgroup = people[, "Wine"], opacity = 0.5))
-   expect_silent(tf(type = "h", cgroup = people[1, ], opacity = 0.5))
+   expect_silent(tf(type = "h", cgroup = as.numeric(people[1, ]), opacity = 0.5))
 })
 
 par(mfrow = c(4, 2))
diff --git a/tests/testthat/test-pcv.R b/tests/testthat/test-pcv.R
new file mode 100644
index 0000000..e22e140
--- /dev/null
+++ b/tests/testthat/test-pcv.R
@@ -0,0 +1,42 @@
+#####################################################
+# Tests for Procrustes Cross-Validation methods     #
+#####################################################
+
+setup({
+   pdf(file = tempfile("mdatools-test-pcv-", fileext = ".pdf"))
+   sink(tempfile("mdatools-test-pcv-", fileext = ".txt"), append = FALSE, split = FALSE)
+})
+
+teardown({
+   dev.off()
+   sink()
+})
+
+# Simple tests for the PCV R implementation
+I <- 100
+J <- 50
+A <- 10
+K <- 4
+set.seed(42)
+x <- matrix(rnorm(I * J), I, J)
+
+params <- list()
+params[[1]] <- list(x = x)
+params[[2]] <- list(x = x, ncomp = 1)
+params[[3]] <- list(x = x, ncomp = A)
+params[[4]] <- list(x = x, ncomp = A, nseg = 4)
+params[[5]] <- list(x = x, ncomp = A, nseg = 10)
+params[[6]] <- list(x = x, ncomp = A, nseg = nrow(x))
+params[[7]] <- list(x = x, ncomp = A, nseg = 4, scale = TRUE)
+params[[8]] <- list(x = x, ncomp = A, nseg = 10, scale = TRUE)
+params[[9]] <- list(x = x, ncomp = A, nseg = nrow(x), scale = TRUE)
+
+context("pcv: for PCA")
+for (i in seq_along(params)) {
+   test_that(paste0("pcv() for PCA works well with parameters set #", i), {
+      x.pv <- do.call(pcv, params[[i]])
+      expect_equal(nrow(x.pv), nrow(x))
+      expect_equal(ncol(x.pv), ncol(x))
+      expect_gt(ks.test(x, x.pv)$p.value, 0.01)
+   })
+}
diff --git a/tests/testthat/test-plsda.R b/tests/testthat/test-plsda.R
index 5448b09..c8c4977 100644
--- a/tests/testthat/test-plsda.R
+++ b/tests/testthat/test-plsda.R
@@ -110,3 +110,16 @@ test_that("predictions work fine", {
    summary(res)
    print(res)
 })
+
+test_that("Bug #95: plots for PLS regression can be also used with PLS-DA objects", {
+   data(people)
+   X <- people[, -c(3, 9)]
+   c <- as.factor(people[, 9])
+   m <- plsda(X, c, 3, cv = 1)
+
+   expect_equal(class(m$res$cal), c("plsdares", "classres", "plsres", "regres"))
+
+   par(mfrow = c(1, 2))
+   expect_silent(plotRMSE(m))
+   expect_silent(plotPredictions.regmodel(m))
+})
\ No newline at end of file
diff --git a/tests/testthat/test-prep.R b/tests/testthat/test-prep.R
index a866d1f..3a33763 100644
--- a/tests/testthat/test-prep.R
+++ b/tests/testthat/test-prep.R
@@ -57,7 +57,7 @@ test_that("Full autoscaling with limits for max.cov is correct for negative data
 })
 
 
-context("prep: snv, msc, norm")
+context("prep: snv, msc, norm, km")
 
 test_that("SNV works correctly", {
    spectra <- simdata$spectra.c
@@ -110,6 +110,23 @@ test_that("Normalization to unit length works correctly", {
    expect_equivalent(apply(pspectra, 1, function(x) sqrt(sum(x^2))), rep(1, nrow(pspectra)))
 })
 
+
+
+test_that("Kubelka-Munk works correctly", {
+   spectra <- simdata$spectra.c
+   spectra <- spectra - min(spectra) + 0.01 * max(spectra)
+   expect_silent(pspectra <- prep.ref2km(spectra))
+   expect_equal(mda.getattr(spectra), mda.getattr(pspectra))
+   expect_equal(dim(spectra), dim(pspectra))
+   expect_true(all(pspectra > 0))
+
+   spectra[10, 10] <- 0
+   expect_error(prep.km(spectra))
+
+   spectra[10, 10] <- -1
+   expect_error(prep.km(spectra))
+})
+
 context("prep: savgol")
 
 test_that("SavGol smoothing works correctly", {