R version 4.0.5 (2021-03-31) -- "Shake and Throw"
Copyright (C) 2021 The R Foundation for Statistical Computing
Platform: i386-w64-mingw32/i386 (32-bit)

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.

> library("aroma.light")
aroma.light v3.20.0 successfully loaded. See ?aroma.light for help.
> 
> X <- matrix(1:8, nrow=4, ncol=2)
> X[2,2] <- NA_integer_
> 
> print(X)
     [,1] [,2]
[1,]    1    5
[2,]    2   NA
[3,]    3    7
[4,]    4    8
> 
> # Returns a 4x2 matrix
> print(backtransformAffine(X, a=c(1,5)))
     [,1] [,2]
[1,]    0    0
[2,]    1   NA
[3,]    2    2
[4,]    3    3
> 
> # Returns a 4x2 matrix
> print(backtransformAffine(X, b=c(1,1/2)))
     [,1] [,2]
[1,]    1   10
[2,]    2   NA
[3,]    3   14
[4,]    4   16
> 
> # Returns a 4x2 matrix
> print(backtransformAffine(X, a=matrix(1:4,ncol=1)))
     [,1] [,2]
[1,]    0    4
[2,]    0   NA
[3,]    0    4
[4,]    0    4
> 
> # Returns a 4x2 matrix
> print(backtransformAffine(X, a=matrix(1:3,ncol=1)))
     [,1] [,2]
[1,]    0    4
[2,]    0   NA
[3,]    0    4
[4,]    3    7
> 
> # Returns a 4x2 matrix
> print(backtransformAffine(X, a=matrix(1:2,ncol=1), b=c(1,2)))
     [,1] [,2]
[1,]    0    2
[2,]    0   NA
[3,]    2    3
[4,]    2    3
> 
> # Returns a 4x1 matrix
> print(backtransformAffine(X, b=c(1,1/2), project=TRUE))
     [,1]
[1,]  2.8
[2,]  1.6
[3,]  5.2
[4,]  6.4
> 
> # If the columns of X are identical, and a identity
> # backtransformation is applied and projected, the
> # same matrix is returned.
> X <- matrix(1:4, nrow=4, ncol=3)
> Y <- backtransformAffine(X, b=c(1,1,1), project=TRUE)
> print(X)
     [,1] [,2] [,3]
[1,]    1    1    1
[2,]    2    2    2
[3,]    3    3    3
[4,]    4    4    4
> print(Y)
     [,1]
[1,]    1
[2,]    2
[3,]    3
[4,]    4
> stopifnot(sum(X[,1]-Y) <= .Machine$double.eps)
> 
> 
> # If the columns of X are identical, and a identity
> # backtransformation is applied and projected, the
> # same matrix is returned.
> X <- matrix(1:4, nrow=4, ncol=3)
> X[,2] <- X[,2]*2; X[,3] <- X[,3]*3
> print(X)
     [,1] [,2] [,3]
[1,]    1    2    3
[2,]    2    4    6
[3,]    3    6    9
[4,]    4    8   12
> Y <- backtransformAffine(X, b=c(1,2,3))
> print(Y)
     [,1] [,2] [,3]
[1,]    1    1    1
[2,]    2    2    2
[3,]    3    3    3
[4,]    4    4    4
> Y <- backtransformAffine(X, b=c(1,2,3), project=TRUE)
> print(Y)
     [,1]
[1,]    1
[2,]    2
[3,]    3
[4,]    4
> stopifnot(sum(X[,1]-Y) <= .Machine$double.eps)
> 
> proc.time()
   user  system elapsed 
   0.29    0.01    0.29 

R version 4.0.5 (2021-03-31) -- "Shake and Throw"
Copyright (C) 2021 The R Foundation for Statistical Computing
Platform: x86_64-w64-mingw32/x64 (64-bit)

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.

> library("aroma.light")
aroma.light v3.20.0 successfully loaded. See ?aroma.light for help.
> 
> X <- matrix(1:8, nrow=4, ncol=2)
> X[2,2] <- NA_integer_
> 
> print(X)
     [,1] [,2]
[1,]    1    5
[2,]    2   NA
[3,]    3    7
[4,]    4    8
> 
> # Returns a 4x2 matrix
> print(backtransformAffine(X, a=c(1,5)))
     [,1] [,2]
[1,]    0    0
[2,]    1   NA
[3,]    2    2
[4,]    3    3
> 
> # Returns a 4x2 matrix
> print(backtransformAffine(X, b=c(1,1/2)))
     [,1] [,2]
[1,]    1   10
[2,]    2   NA
[3,]    3   14
[4,]    4   16
> 
> # Returns a 4x2 matrix
> print(backtransformAffine(X, a=matrix(1:4,ncol=1)))
     [,1] [,2]
[1,]    0    4
[2,]    0   NA
[3,]    0    4
[4,]    0    4
> 
> # Returns a 4x2 matrix
> print(backtransformAffine(X, a=matrix(1:3,ncol=1)))
     [,1] [,2]
[1,]    0    4
[2,]    0   NA
[3,]    0    4
[4,]    3    7
> 
> # Returns a 4x2 matrix
> print(backtransformAffine(X, a=matrix(1:2,ncol=1), b=c(1,2)))
     [,1] [,2]
[1,]    0    2
[2,]    0   NA
[3,]    2    3
[4,]    2    3
> 
> # Returns a 4x1 matrix
> print(backtransformAffine(X, b=c(1,1/2), project=TRUE))
     [,1]
[1,]  2.8
[2,]  1.6
[3,]  5.2
[4,]  6.4
> 
> # If the columns of X are identical, and a identity
> # backtransformation is applied and projected, the
> # same matrix is returned.
> X <- matrix(1:4, nrow=4, ncol=3)
> Y <- backtransformAffine(X, b=c(1,1,1), project=TRUE)
> print(X)
     [,1] [,2] [,3]
[1,]    1    1    1
[2,]    2    2    2
[3,]    3    3    3
[4,]    4    4    4
> print(Y)
     [,1]
[1,]    1
[2,]    2
[3,]    3
[4,]    4
> stopifnot(sum(X[,1]-Y) <= .Machine$double.eps)
> 
> 
> # If the columns of X are identical, and a identity
> # backtransformation is applied and projected, the
> # same matrix is returned.
> X <- matrix(1:4, nrow=4, ncol=3)
> X[,2] <- X[,2]*2; X[,3] <- X[,3]*3
> print(X)
     [,1] [,2] [,3]
[1,]    1    2    3
[2,]    2    4    6
[3,]    3    6    9
[4,]    4    8   12
> Y <- backtransformAffine(X, b=c(1,2,3))
> print(Y)
     [,1] [,2] [,3]
[1,]    1    1    1
[2,]    2    2    2
[3,]    3    3    3
[4,]    4    4    4
> Y <- backtransformAffine(X, b=c(1,2,3), project=TRUE)
> print(Y)
     [,1]
[1,]    1
[2,]    2
[3,]    3
[4,]    4
> stopifnot(sum(X[,1]-Y) <= .Machine$double.eps)
> 
> proc.time()
   user  system elapsed 
   0.31    0.07    0.37 

R version 4.0.5 (2021-03-31) -- "Shake and Throw"
Copyright (C) 2021 The R Foundation for Statistical Computing
Platform: i386-w64-mingw32/i386 (32-bit)

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.

> library("aroma.light")
aroma.light v3.20.0 successfully loaded. See ?aroma.light for help.
> 
> # Consider the case where K=4 measurements have been done
> # for the same underlying signals 'x'.  The different measurements
> # have different systematic variation
> #
> #   y_k = f(x_k) + eps_k; k = 1,...,K.
> #
> # In this example, we assume non-linear measurement functions
> #
> #   f(x) = a + b*x + x^c + eps(b*x)
> #
> # where 'a' is an offset, 'b' a scale factor, and 'c' an exponential.
> # We also assume heteroscedastic zero-mean noise with standard
> # deviation proportional to the rescaled underlying signal 'x'.
> #
> # Furthermore, we assume that measurements k=2 and k=3 undergo the
> # same transformation, which may illustrate that the come from
> # the same batch. However, when *fitting* the model below we
> # will assume they are independent.
> 
> # Transforms
> a <- c(2, 15, 15,   3)
> b <- c(2,  3,  3,   4)
> c <- c(1,  2,  2, 1/2)
> K <- length(a)
> 
> # The true signal
> N <- 1000
> x <- rexp(N)
> 
> # The noise
> bX <- outer(b,x)
> E <- apply(bX, MARGIN=2, FUN=function(x) rnorm(K, mean=0, sd=0.1*x))
> 
> # The transformed signals with noise
> Xc <- t(sapply(c, FUN=function(c) x^c))
> Y <- a + bX + Xc + E
> Y <- t(Y)
> 
> 
> 
> # - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
> # Fit principal curve
> # - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
> # Fit principal curve through Y = (y_1, y_2, ..., y_K)
> fit <- fitPrincipalCurve(Y)
> 
> # Flip direction of 'lambda'?
> rho <- cor(fit$lambda, Y[,1], use="complete.obs")
> flip <- (rho < 0)
> if (flip) {
+   fit$lambda <- max(fit$lambda, na.rm=TRUE)-fit$lambda
+ }
> 
> L <- ncol(fit$s)
> 
> # - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
> # Backtransform data according to model fit
> # - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
> # Backtransform toward the principal curve (the "common scale")
> YN1 <- backtransformPrincipalCurve(Y, fit=fit)
> stopifnot(ncol(YN1) == K)
> 
> 
> # Backtransform toward the first dimension
> YN2 <- backtransformPrincipalCurve(Y, fit=fit, targetDimension=1)
> stopifnot(ncol(YN2) == K)
> 
> 
> # Backtransform toward the last (fitted) dimension
> YN3 <- backtransformPrincipalCurve(Y, fit=fit, targetDimension=L)
> stopifnot(ncol(YN3) == K)
> 
> 
> # Backtransform toward the third dimension (dimension by dimension)
> # Note, this assumes that K == L.
> YN4 <- Y
> for (cc in 1:L) {
+   YN4[,cc] <- backtransformPrincipalCurve(Y, fit=fit,
+                                   targetDimension=1, dimensions=cc)
+ }
> stopifnot(identical(YN4, YN2))
> 
> 
> # Backtransform a subset toward the first dimension
> # Note, this assumes that K == L.
> YN5 <- backtransformPrincipalCurve(Y, fit=fit,
+                                targetDimension=1, dimensions=2:3)
> stopifnot(identical(YN5, YN2[,2:3]))
> stopifnot(ncol(YN5) == 2)
> 
> 
> # Extract signals from measurement #2 and backtransform according
> # its model fit.  Signals are standardized to target dimension 1.
> y6 <- Y[,2,drop=FALSE]
> yN6 <- backtransformPrincipalCurve(y6, fit=fit, dimensions=2,
+                                                targetDimension=1)
> stopifnot(identical(yN6, YN2[,2,drop=FALSE]))
> stopifnot(ncol(yN6) == 1)
> 
> 
> # Extract signals from measurement #2 and backtransform according
> # the the model fit of measurement #3 (because we believe these
> # two have undergone very similar transformations.
> # Signals are standardized to target dimension 1.
> y7 <- Y[,2,drop=FALSE]
> yN7 <- backtransformPrincipalCurve(y7, fit=fit, dimensions=3,
+                                                targetDimension=1)
> stopifnot(ncol(yN7) == 1)
> 
> rho <- cor(yN7, yN6)
> print(rho)
          [,1]
[1,] 0.9999635
> stopifnot(rho > 0.999)
> 
> proc.time()
   user  system elapsed 
   0.75    0.10    0.85 

R version 4.0.5 (2021-03-31) -- "Shake and Throw"
Copyright (C) 2021 The R Foundation for Statistical Computing
Platform: x86_64-w64-mingw32/x64 (64-bit)

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.

> library("aroma.light")
aroma.light v3.20.0 successfully loaded. See ?aroma.light for help.
> 
> # Consider the case where K=4 measurements have been done
> # for the same underlying signals 'x'.  The different measurements
> # have different systematic variation
> #
> #   y_k = f(x_k) + eps_k; k = 1,...,K.
> #
> # In this example, we assume non-linear measurement functions
> #
> #   f(x) = a + b*x + x^c + eps(b*x)
> #
> # where 'a' is an offset, 'b' a scale factor, and 'c' an exponential.
> # We also assume heteroscedastic zero-mean noise with standard
> # deviation proportional to the rescaled underlying signal 'x'.
> #
> # Furthermore, we assume that measurements k=2 and k=3 undergo the
> # same transformation, which may illustrate that the come from
> # the same batch. However, when *fitting* the model below we
> # will assume they are independent.
> 
> # Transforms
> a <- c(2, 15, 15,   3)
> b <- c(2,  3,  3,   4)
> c <- c(1,  2,  2, 1/2)
> K <- length(a)
> 
> # The true signal
> N <- 1000
> x <- rexp(N)
> 
> # The noise
> bX <- outer(b,x)
> E <- apply(bX, MARGIN=2, FUN=function(x) rnorm(K, mean=0, sd=0.1*x))
> 
> # The transformed signals with noise
> Xc <- t(sapply(c, FUN=function(c) x^c))
> Y <- a + bX + Xc + E
> Y <- t(Y)
> 
> 
> 
> # - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
> # Fit principal curve
> # - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
> # Fit principal curve through Y = (y_1, y_2, ..., y_K)
> fit <- fitPrincipalCurve(Y)
> 
> # Flip direction of 'lambda'?
> rho <- cor(fit$lambda, Y[,1], use="complete.obs")
> flip <- (rho < 0)
> if (flip) {
+   fit$lambda <- max(fit$lambda, na.rm=TRUE)-fit$lambda
+ }
> 
> L <- ncol(fit$s)
> 
> # - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
> # Backtransform data according to model fit
> # - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
> # Backtransform toward the principal curve (the "common scale")
> YN1 <- backtransformPrincipalCurve(Y, fit=fit)
> stopifnot(ncol(YN1) == K)
> 
> 
> # Backtransform toward the first dimension
> YN2 <- backtransformPrincipalCurve(Y, fit=fit, targetDimension=1)
> stopifnot(ncol(YN2) == K)
> 
> 
> # Backtransform toward the last (fitted) dimension
> YN3 <- backtransformPrincipalCurve(Y, fit=fit, targetDimension=L)
> stopifnot(ncol(YN3) == K)
> 
> 
> # Backtransform toward the third dimension (dimension by dimension)
> # Note, this assumes that K == L.
> YN4 <- Y
> for (cc in 1:L) {
+   YN4[,cc] <- backtransformPrincipalCurve(Y, fit=fit,
+                                   targetDimension=1, dimensions=cc)
+ }
> stopifnot(identical(YN4, YN2))
> 
> 
> # Backtransform a subset toward the first dimension
> # Note, this assumes that K == L.
> YN5 <- backtransformPrincipalCurve(Y, fit=fit,
+                                targetDimension=1, dimensions=2:3)
> stopifnot(identical(YN5, YN2[,2:3]))
> stopifnot(ncol(YN5) == 2)
> 
> 
> # Extract signals from measurement #2 and backtransform according
> # its model fit.  Signals are standardized to target dimension 1.
> y6 <- Y[,2,drop=FALSE]
> yN6 <- backtransformPrincipalCurve(y6, fit=fit, dimensions=2,
+                                                targetDimension=1)
> stopifnot(identical(yN6, YN2[,2,drop=FALSE]))
> stopifnot(ncol(yN6) == 1)
> 
> 
> # Extract signals from measurement #2 and backtransform according
> # the the model fit of measurement #3 (because we believe these
> # two have undergone very similar transformations.
> # Signals are standardized to target dimension 1.
> y7 <- Y[,2,drop=FALSE]
> yN7 <- backtransformPrincipalCurve(y7, fit=fit, dimensions=3,
+                                                targetDimension=1)
> stopifnot(ncol(yN7) == 1)
> 
> rho <- cor(yN7, yN6)
> print(rho)
          [,1]
[1,] 0.9999822
> stopifnot(rho > 0.999)
> 
> proc.time()
   user  system elapsed 
   0.89    0.04    0.92 

R version 4.0.5 (2021-03-31) -- "Shake and Throw"
Copyright (C) 2021 The R Foundation for Statistical Computing
Platform: i386-w64-mingw32/i386 (32-bit)

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.

> library("aroma.light")
aroma.light v3.20.0 successfully loaded. See ?aroma.light for help.
> 
> layout(matrix(1:3, ncol=1))
> par(mar=c(2,4,4,1)+0.1)
> 
> # - - - - - - - - - - - - - - - - - - - - - - - - - - - -
> # A bimodal distribution
> # - - - - - - - - - - - - - - - - - - - - - - - - - - - -
> xAA <- rnorm(n=10000, mean=0, sd=0.1)
> xBB <- rnorm(n=10000, mean=1, sd=0.1)
> x <- c(xAA,xBB)
> fit <- findPeaksAndValleys(x)
> print(fit)
    type           x      density
1   peak -0.00873527 1.7140702200
2 valley  0.48718983 0.0002852856
3   peak  0.99572319 1.6736759944
> calls <- callNaiveGenotypes(x, cn=rep(1,length(x)), verbose=-20)
Calling genotypes from allele B fractions (BAFs)...
 Fitting naive genotype model...
  Fitting naive genotype model from normal allele B fractions (BAFs)...
   Flavor: density
   Censoring BAFs...
    Before:
          Min.    1st Qu.     Median       Mean    3rd Qu.       Max. 
    -0.3649774 -0.0009633  0.4949453  0.5003704  1.0004612  1.4023984 
    [1] 20000
    After:
          Min.    1st Qu.     Median       Mean    3rd Qu.       Max. 
          -Inf -0.0009633  0.4949453             1.0004612        Inf 
    [1] 16848
   Censoring BAFs...done
   Copy number level #1 (C=1) of 1...
    Identified extreme points in density of BAF:
        type          x     density
    1   peak 0.01101669 1.652971537
    2 valley 0.49488869 0.003803567
    3   peak 0.97876069 1.624315408
    Local minimas ("valleys") in BAF:
        type         x     density
    2 valley 0.4948887 0.003803567
   Copy number level #1 (C=1) of 1...done
  Fitting naive genotype model from normal allele B fractions (BAFs)...done
  [[1]]
  [[1]]$flavor
  [1] "density"
  
  [[1]]$cn
  [1] 1
  
  [[1]]$nbrOfGenotypeGroups
  [1] 2
  
  [[1]]$tau
  [1] 0.4948887
  
  [[1]]$n
  [1] 16848
  
  [[1]]$fit
      type          x     density
  1   peak 0.01101669 1.652971537
  2 valley 0.49488869 0.003803567
  3   peak 0.97876069 1.624315408
  
  [[1]]$fitValleys
      type         x     density
  2 valley 0.4948887 0.003803567
  
  
  attr(,"class")
  [1] "NaiveGenotypeModelFit" "list"                 
 Fitting naive genotype model...done
 Copy number level #1 (C=1) of 1...
  Model fit:
  $flavor
  [1] "density"
  
  $cn
  [1] 1
  
  $nbrOfGenotypeGroups
  [1] 2
  
  $tau
  [1] 0.4948887
  
  $n
  [1] 16848
  
  $fit
      type          x     density
  1   peak 0.01101669 1.652971537
  2 valley 0.49488869 0.003803567
  3   peak 0.97876069 1.624315408
  
  $fitValleys
      type         x     density
  2 valley 0.4948887 0.003803567
  
  Genotype threshholds [1]: 0.494888688269728
  TCN=1 => BAF in {0,1}.
  Call regions: A = (-Inf,0.495], B = (0.495,+Inf)
 Copy number level #1 (C=1) of 1...done
Calling genotypes from allele B fractions (BAFs)...done
> xc <- split(x, calls)
> print(table(calls))
calls
    0     1 
10000 10000 
> xx <- c(list(x),xc)
> plotDensity(xx, adjust=1.5, lwd=2, col=seq_along(xx), main="(AA,BB)")
> abline(v=fit$x)
> 
> # - - - - - - - - - - - - - - - - - - - - - - - - - - - -
> # A trimodal distribution with missing values
> # - - - - - - - - - - - - - - - - - - - - - - - - - - - -
> xAB <- rnorm(n=10000, mean=1/2, sd=0.1)
> x <- c(xAA,xAB,xBB)
> x[sample(length(x), size=0.05*length(x))] <- NA_real_
> x[sample(length(x), size=0.01*length(x))] <- -Inf
> x[sample(length(x), size=0.01*length(x))] <- +Inf
> fit <- findPeaksAndValleys(x)
> print(fit)
    type            x   density
1   peak -0.007562787 1.1944864
2 valley  0.246500174 0.1863059
3   peak  0.492497645 1.1727941
4 valley  0.742527860 0.1849325
5   peak  0.996590821 1.1578546
> calls <- callNaiveGenotypes(x)
> xc <- split(x, calls)
> print(table(calls))
calls
   0  0.5    1 
9605 9261 9669 
> xx <- c(list(x),xc)
> plotDensity(xx, adjust=1.5, lwd=2, col=seq_along(xx), main="(AA,AB,BB)")
> abline(v=fit$x)
> 
> # - - - - - - - - - - - - - - - - - - - - - - - - - - - -
> # A trimodal distribution with clear separation
> # - - - - - - - - - - - - - - - - - - - - - - - - - - - -
> xAA <- rnorm(n=10000, mean=0, sd=0.02)
> xAB <- rnorm(n=10000, mean=1/2, sd=0.02)
> xBB <- rnorm(n=10000, mean=1, sd=0.02)
> x <- c(xAA,xAB,xBB)
> fit <- findPeaksAndValleys(x)
> print(fit)
    type            x      density
1   peak -0.001930485 2.606041e+00
2 valley  0.248589778 3.066903e-05
3   peak  0.496295206 2.608020e+00
4 valley  0.746815468 3.319422e-05
5   peak  0.997335731 2.608589e+00
> calls <- callNaiveGenotypes(x)
> xc <- split(x, calls)
> print(table(calls))
calls
    0   0.5     1 
10000 10000 10000 
> xx <- c(list(x),xc)
> plotDensity(xx, adjust=1.5, lwd=2, col=seq_along(xx), main="(AA',AB',BB')")
> abline(v=fit$x)
> 
> proc.time()
   user  system elapsed 
   0.65    0.06    0.70 

R version 4.0.5 (2021-03-31) -- "Shake and Throw"
Copyright (C) 2021 The R Foundation for Statistical Computing
Platform: x86_64-w64-mingw32/x64 (64-bit)

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.

> library("aroma.light")
aroma.light v3.20.0 successfully loaded. See ?aroma.light for help.
> 
> layout(matrix(1:3, ncol=1))
> par(mar=c(2,4,4,1)+0.1)
> 
> # - - - - - - - - - - - - - - - - - - - - - - - - - - - -
> # A bimodal distribution
> # - - - - - - - - - - - - - - - - - - - - - - - - - - - -
> xAA <- rnorm(n=10000, mean=0, sd=0.1)
> xBB <- rnorm(n=10000, mean=1, sd=0.1)
> x <- c(xAA,xBB)
> fit <- findPeaksAndValleys(x)
> print(fit)
    type           x      density
1   peak -0.00550228 1.6746286836
2 valley  0.49770694 0.0005825472
3   peak  0.99679149 1.6834880834
> calls <- callNaiveGenotypes(x, cn=rep(1,length(x)), verbose=-20)
Calling genotypes from allele B fractions (BAFs)...
 Fitting naive genotype model...
  Fitting naive genotype model from normal allele B fractions (BAFs)...
   Flavor: density
   Censoring BAFs...
    Before:
          Min.    1st Qu.     Median       Mean    3rd Qu.       Max. 
    -0.3679160 -0.0003755  0.4807680  0.4998221  0.9996986  1.3592052 
    [1] 20000
    After:
          Min.    1st Qu.     Median       Mean    3rd Qu.       Max. 
          -Inf -0.0003755  0.4807680             0.9996986        Inf 
    [1] 16758
   Censoring BAFs...done
   Copy number level #1 (C=1) of 1...
    Identified extreme points in density of BAF:
        type          x     density
    1   peak 0.01101243 1.634459850
    2 valley 0.49824430 0.004428433
    3   peak 0.97861374 1.637361239
    Local minimas ("valleys") in BAF:
        type         x     density
    2 valley 0.4982443 0.004428433
   Copy number level #1 (C=1) of 1...done
  Fitting naive genotype model from normal allele B fractions (BAFs)...done
  [[1]]
  [[1]]$flavor
  [1] "density"
  
  [[1]]$cn
  [1] 1
  
  [[1]]$nbrOfGenotypeGroups
  [1] 2
  
  [[1]]$tau
  [1] 0.4982443
  
  [[1]]$n
  [1] 16758
  
  [[1]]$fit
      type          x     density
  1   peak 0.01101243 1.634459850
  2 valley 0.49824430 0.004428433
  3   peak 0.97861374 1.637361239
  
  [[1]]$fitValleys
      type         x     density
  2 valley 0.4982443 0.004428433
  
  
  attr(,"class")
  [1] "NaiveGenotypeModelFit" "list"                 
 Fitting naive genotype model...done
 Copy number level #1 (C=1) of 1...
  Model fit:
  $flavor
  [1] "density"
  
  $cn
  [1] 1
  
  $nbrOfGenotypeGroups
  [1] 2
  
  $tau
  [1] 0.4982443
  
  $n
  [1] 16758
  
  $fit
      type          x     density
  1   peak 0.01101243 1.634459850
  2 valley 0.49824430 0.004428433
  3   peak 0.97861374 1.637361239
  
  $fitValleys
      type         x     density
  2 valley 0.4982443 0.004428433
  
  Genotype threshholds [1]: 0.498244295187841
  TCN=1 => BAF in {0,1}.
  Call regions: A = (-Inf,0.498], B = (0.498,+Inf)
 Copy number level #1 (C=1) of 1...done
Calling genotypes from allele B fractions (BAFs)...done
> xc <- split(x, calls)
> print(table(calls))
calls
    0     1 
10000 10000 
> xx <- c(list(x),xc)
> plotDensity(xx, adjust=1.5, lwd=2, col=seq_along(xx), main="(AA,BB)")
> abline(v=fit$x)
> 
> # - - - - - - - - - - - - - - - - - - - - - - - - - - - -
> # A trimodal distribution with missing values
> # - - - - - - - - - - - - - - - - - - - - - - - - - - - -
> xAB <- rnorm(n=10000, mean=1/2, sd=0.1)
> x <- c(xAA,xAB,xBB)
> x[sample(length(x), size=0.05*length(x))] <- NA_real_
> x[sample(length(x), size=0.01*length(x))] <- -Inf
> x[sample(length(x), size=0.01*length(x))] <- +Inf
> fit <- findPeaksAndValleys(x)
> print(fit)
    type           x   density
1   peak -0.00459096 1.1581739
2 valley  0.24849264 0.1832683
3   peak  0.48971296 1.1903049
4 valley  0.74279656 0.1828331
5   peak  0.99588017 1.1806170
> calls <- callNaiveGenotypes(x)
> xc <- split(x, calls)
> print(table(calls))
calls
   0  0.5    1 
9580 9278 9665 
> xx <- c(list(x),xc)
> plotDensity(xx, adjust=1.5, lwd=2, col=seq_along(xx), main="(AA,AB,BB)")
> abline(v=fit$x)
> 
> # - - - - - - - - - - - - - - - - - - - - - - - - - - - -
> # A trimodal distribution with clear separation
> # - - - - - - - - - - - - - - - - - - - - - - - - - - - -
> xAA <- rnorm(n=10000, mean=0, sd=0.02)
> xAB <- rnorm(n=10000, mean=1/2, sd=0.02)
> xBB <- rnorm(n=10000, mean=1, sd=0.02)
> x <- c(xAA,xAB,xBB)
> fit <- findPeaksAndValleys(x)
> print(fit)
    type            x      density
1   peak -0.003861662 2.606226e+00
2 valley  0.245916144 3.420654e-05
3   peak  0.495693951 2.602026e+00
4 valley  0.745471758 3.503801e-05
5   peak  0.998056057 2.606695e+00
> calls <- callNaiveGenotypes(x)
> xc <- split(x, calls)
> print(table(calls))
calls
    0   0.5     1 
10000 10000 10000 
> xx <- c(list(x),xc)
> plotDensity(xx, adjust=1.5, lwd=2, col=seq_along(xx), main="(AA',AB',BB')")
> abline(v=fit$x)
> 
> proc.time()
   user  system elapsed 
   0.71    0.09    0.79 

R version 4.0.5 (2021-03-31) -- "Shake and Throw"
Copyright (C) 2021 The R Foundation for Statistical Computing
Platform: i386-w64-mingw32/i386 (32-bit)

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.

> library("aroma.light")
aroma.light v3.20.0 successfully loaded. See ?aroma.light for help.
> 
> for (zzz in 0) {
+ 
+ # This example requires plot3d() in R.basic [http://www.braju.com/R/]
+ if (!require(pkgName <- "R.basic", character.only=TRUE)) break
+ 
+ layout(matrix(1:4, nrow=2, ncol=2, byrow=TRUE))
+ 
+ ############################################################
+ # Lines in two-dimensions
+ ############################################################
+ x <- list(a=c(1,0), b=c(1,2))
+ y <- list(a=c(0,2), b=c(1,1))
+ fit <- distanceBetweenLines(ax=x$a, bx=x$b, ay=y$a, by=y$b)
+ 
+ xlim <- ylim <- c(-1,8)
+ plot(NA, xlab="", ylab="", xlim=ylim, ylim=ylim)
+ 
+ # Highlight the offset coordinates for both lines
+ points(t(x$a), pch="+", col="red")
+ text(t(x$a), label=expression(a[x]), adj=c(-1,0.5))
+ points(t(y$a), pch="+", col="blue")
+ text(t(y$a), label=expression(a[y]), adj=c(-1,0.5))
+ 
+ v <- c(-1,1)*10
+ xv <- list(x=x$a[1]+x$b[1]*v, y=x$a[2]+x$b[2]*v)
+ yv <- list(x=y$a[1]+y$b[1]*v, y=y$a[2]+y$b[2]*v)
+ 
+ lines(xv, col="red")
+ lines(yv, col="blue")
+ 
+ points(t(fit$xs), cex=2.0, col="red")
+ text(t(fit$xs), label=expression(x(s)), adj=c(+2,0.5))
+ points(t(fit$yt), cex=1.5, col="blue")
+ text(t(fit$yt), label=expression(y(t)), adj=c(-1,0.5))
+ print(fit)
+ 
+ 
+ ############################################################
+ # Lines in three-dimensions
+ ############################################################
+ x <- list(a=c(0,0,0), b=c(1,1,1))  # The 'diagonal'
+ y <- list(a=c(2,1,2), b=c(2,1,3))  # A 'fitted' line
+ fit <- distanceBetweenLines(ax=x$a, bx=x$b, ay=y$a, by=y$b)
+ 
+ xlim <- ylim <- zlim <- c(-1,3)
+ dummy <- t(c(1,1,1))*100
+ 
+ # Coordinates for the lines in 3d
+ v <- seq(-10,10, by=1)
+ xv <- list(x=x$a[1]+x$b[1]*v, y=x$a[2]+x$b[2]*v, z=x$a[3]+x$b[3]*v)
+ yv <- list(x=y$a[1]+y$b[1]*v, y=y$a[2]+y$b[2]*v, z=y$a[3]+y$b[3]*v)
+ 
+ for (theta in seq(30,140,length.out=3)) {
+   plot3d(dummy, theta=theta, phi=30, xlab="", ylab="", zlab="",
+                              xlim=ylim, ylim=ylim, zlim=zlim)
+ 
+   # Highlight the offset coordinates for both lines
+   points3d(t(x$a), pch="+", col="red")
+   text3d(t(x$a), label=expression(a[x]), adj=c(-1,0.5))
+   points3d(t(y$a), pch="+", col="blue")
+   text3d(t(y$a), label=expression(a[y]), adj=c(-1,0.5))
+ 
+   # Draw the lines
+   lines3d(xv, col="red")
+   lines3d(yv, col="blue")
+ 
+   # Draw the two points that are closest to each other
+   points3d(t(fit$xs), cex=2.0, col="red")
+   text3d(t(fit$xs), label=expression(x(s)), adj=c(+2,0.5))
+   points3d(t(fit$yt), cex=1.5, col="blue")
+   text3d(t(fit$yt), label=expression(y(t)), adj=c(-1,0.5))
+ 
+   # Draw the distance between the two points
+   lines3d(rbind(fit$xs,fit$yt), col="purple", lwd=2)
+ }
+ 
+ print(fit)
+ 
+ } # for (zzz in 0)
Loading required package: R.basic
Warning message:
In library(package, lib.loc = lib.loc, character.only = TRUE, logical.return = TRUE,  :
  there is no package called 'R.basic'
> rm(zzz)
> 
> proc.time()
   user  system elapsed 
   0.39    0.03    0.40 

R version 4.0.5 (2021-03-31) -- "Shake and Throw"
Copyright (C) 2021 The R Foundation for Statistical Computing
Platform: x86_64-w64-mingw32/x64 (64-bit)

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.

> library("aroma.light")
aroma.light v3.20.0 successfully loaded. See ?aroma.light for help.
> 
> for (zzz in 0) {
+ 
+ # This example requires plot3d() in R.basic [http://www.braju.com/R/]
+ if (!require(pkgName <- "R.basic", character.only=TRUE)) break
+ 
+ layout(matrix(1:4, nrow=2, ncol=2, byrow=TRUE))
+ 
+ ############################################################
+ # Lines in two-dimensions
+ ############################################################
+ x <- list(a=c(1,0), b=c(1,2))
+ y <- list(a=c(0,2), b=c(1,1))
+ fit <- distanceBetweenLines(ax=x$a, bx=x$b, ay=y$a, by=y$b)
+ 
+ xlim <- ylim <- c(-1,8)
+ plot(NA, xlab="", ylab="", xlim=ylim, ylim=ylim)
+ 
+ # Highlight the offset coordinates for both lines
+ points(t(x$a), pch="+", col="red")
+ text(t(x$a), label=expression(a[x]), adj=c(-1,0.5))
+ points(t(y$a), pch="+", col="blue")
+ text(t(y$a), label=expression(a[y]), adj=c(-1,0.5))
+ 
+ v <- c(-1,1)*10
+ xv <- list(x=x$a[1]+x$b[1]*v, y=x$a[2]+x$b[2]*v)
+ yv <- list(x=y$a[1]+y$b[1]*v, y=y$a[2]+y$b[2]*v)
+ 
+ lines(xv, col="red")
+ lines(yv, col="blue")
+ 
+ points(t(fit$xs), cex=2.0, col="red")
+ text(t(fit$xs), label=expression(x(s)), adj=c(+2,0.5))
+ points(t(fit$yt), cex=1.5, col="blue")
+ text(t(fit$yt), label=expression(y(t)), adj=c(-1,0.5))
+ print(fit)
+ 
+ 
+ ############################################################
+ # Lines in three-dimensions
+ ############################################################
+ x <- list(a=c(0,0,0), b=c(1,1,1))  # The 'diagonal'
+ y <- list(a=c(2,1,2), b=c(2,1,3))  # A 'fitted' line
+ fit <- distanceBetweenLines(ax=x$a, bx=x$b, ay=y$a, by=y$b)
+ 
+ xlim <- ylim <- zlim <- c(-1,3)
+ dummy <- t(c(1,1,1))*100
+ 
+ # Coordinates for the lines in 3d
+ v <- seq(-10,10, by=1)
+ xv <- list(x=x$a[1]+x$b[1]*v, y=x$a[2]+x$b[2]*v, z=x$a[3]+x$b[3]*v)
+ yv <- list(x=y$a[1]+y$b[1]*v, y=y$a[2]+y$b[2]*v, z=y$a[3]+y$b[3]*v)
+ 
+ for (theta in seq(30,140,length.out=3)) {
+   plot3d(dummy, theta=theta, phi=30, xlab="", ylab="", zlab="",
+                              xlim=ylim, ylim=ylim, zlim=zlim)
+ 
+   # Highlight the offset coordinates for both lines
+   points3d(t(x$a), pch="+", col="red")
+   text3d(t(x$a), label=expression(a[x]), adj=c(-1,0.5))
+   points3d(t(y$a), pch="+", col="blue")
+   text3d(t(y$a), label=expression(a[y]), adj=c(-1,0.5))
+ 
+   # Draw the lines
+   lines3d(xv, col="red")
+   lines3d(yv, col="blue")
+ 
+   # Draw the two points that are closest to each other
+   points3d(t(fit$xs), cex=2.0, col="red")
+   text3d(t(fit$xs), label=expression(x(s)), adj=c(+2,0.5))
+   points3d(t(fit$yt), cex=1.5, col="blue")
+   text3d(t(fit$yt), label=expression(y(t)), adj=c(-1,0.5))
+ 
+   # Draw the distance between the two points
+   lines3d(rbind(fit$xs,fit$yt), col="purple", lwd=2)
+ }
+ 
+ print(fit)
+ 
+ } # for (zzz in 0)
Loading required package: R.basic
Warning message:
In library(package, lib.loc = lib.loc, character.only = TRUE, logical.return = TRUE,  :
  there is no package called 'R.basic'
> rm(zzz)
> 
> proc.time()
   user  system elapsed 
   0.40    0.01    0.40 

R version 4.0.5 (2021-03-31) -- "Shake and Throw"
Copyright (C) 2021 The R Foundation for Statistical Computing
Platform: i386-w64-mingw32/i386 (32-bit)

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.

> library("aroma.light")
aroma.light v3.20.0 successfully loaded. See ?aroma.light for help.
> 
> layout(matrix(1:3, ncol=1))
> par(mar=c(2,4,4,1)+0.1)
> 
> # - - - - - - - - - - - - - - - - - - - - - - - - - - - -
> # A unimodal distribution
> # - - - - - - - - - - - - - - - - - - - - - - - - - - - -
> x1 <- rnorm(n=10000, mean=0, sd=1)
> x <- x1
> fit <- findPeaksAndValleys(x)
> print(fit)
    type          x      density
1   peak -4.6866555 2.758458e-04
2 valley -4.0878218 9.177323e-08
3   peak  0.1584533 3.888324e-01
4 valley  3.4066723 1.036337e-03
5   peak  3.4429652 1.036740e-03
> plot(density(x), lwd=2, main="x1")
> abline(v=fit$x)
> 
> # - - - - - - - - - - - - - - - - - - - - - - - - - - - -
> # A trimodal distribution
> # - - - - - - - - - - - - - - - - - - - - - - - - - - - -
> x2 <- rnorm(n=10000, mean=4, sd=1)
> x3 <- rnorm(n=10000, mean=8, sd=1)
> x <- c(x1,x2,x3)
> fit <- findPeaksAndValleys(x)
> print(fit)
    type           x      density
1   peak -4.67621934 3.466620e-05
2 valley -4.41982328 3.282096e-05
3   peak  0.01216584 1.222382e-01
4 valley  1.91682232 4.467742e-02
5   peak  4.00461885 1.218104e-01
6 valley  5.98253135 4.426994e-02
7   peak  7.99707185 1.250449e-01
> plot(density(x), lwd=2, main="c(x1,x2,x3)")
> abline(v=fit$x)
> 
> 
> # - - - - - - - - - - - - - - - - - - - - - - - - - - - -
> # A trimodal distribution with clear separation
> # - - - - - - - - - - - - - - - - - - - - - - - - - - - -
> x1b <- rnorm(n=10000, mean=0, sd=0.1)
> x2b <- rnorm(n=10000, mean=4, sd=0.1)
> x3b <- rnorm(n=10000, mean=8, sd=0.1)
> x <- c(x1b,x2b,x3b)
> 
> # Illustrating explicit usage of density()
> d <- density(x)
> fit <- findPeaksAndValleys(d, tol=0)
> print(fit)
    type           x      density
1   peak -0.03127426 3.418703e-01
2 valley  1.98912525 1.174349e-06
3   peak  3.98780004 3.431323e-01
4 valley  5.98647482 1.187113e-06
5   peak  7.98514961 3.431534e-01
> plot(d, lwd=2, main="c(x1b,x2b,x3b)")
> abline(v=fit$x)
> 
> proc.time()
   user  system elapsed 
   0.59    0.04    0.62 

R version 4.0.5 (2021-03-31) -- "Shake and Throw"
Copyright (C) 2021 The R Foundation for Statistical Computing
Platform: x86_64-w64-mingw32/x64 (64-bit)

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.

> library("aroma.light")
aroma.light v3.20.0 successfully loaded. See ?aroma.light for help.
> 
> layout(matrix(1:3, ncol=1))
> par(mar=c(2,4,4,1)+0.1)
> 
> # - - - - - - - - - - - - - - - - - - - - - - - - - - - -
> # A unimodal distribution
> # - - - - - - - - - - - - - - - - - - - - - - - - - - - -
> x1 <- rnorm(n=10000, mean=0, sd=1)
> x <- x1
> fit <- findPeaksAndValleys(x)
> print(fit)
    type           x   density
1   peak -0.19560152 0.3793843
2 valley -0.03605744 0.3757474
3   peak  0.23516749 0.3879747
> plot(density(x), lwd=2, main="x1")
> abline(v=fit$x)
> 
> # - - - - - - - - - - - - - - - - - - - - - - - - - - - -
> # A trimodal distribution
> # - - - - - - - - - - - - - - - - - - - - - - - - - - - -
> x2 <- rnorm(n=10000, mean=4, sd=1)
> x3 <- rnorm(n=10000, mean=8, sd=1)
> x <- c(x1,x2,x3)
> fit <- findPeaksAndValleys(x)
> print(fit)
    type          x    density
1   peak 0.01948255 0.12235959
2 valley 1.95187826 0.04467602
3   peak 3.95328811 0.12236843
4 valley 5.95469795 0.04420204
5   peak 7.99061486 0.12417660
> plot(density(x), lwd=2, main="c(x1,x2,x3)")
> abline(v=fit$x)
> 
> 
> # - - - - - - - - - - - - - - - - - - - - - - - - - - - -
> # A trimodal distribution with clear separation
> # - - - - - - - - - - - - - - - - - - - - - - - - - - - -
> x1b <- rnorm(n=10000, mean=0, sd=0.1)
> x2b <- rnorm(n=10000, mean=4, sd=0.1)
> x3b <- rnorm(n=10000, mean=8, sd=0.1)
> x <- c(x1b,x2b,x3b)
> 
> # Illustrating explicit usage of density()
> d <- density(x)
> fit <- findPeaksAndValleys(d, tol=0)
> print(fit)
    type          x      density
1   peak -0.0293411 3.419062e-01
2 valley  1.9692961 1.266080e-06
3   peak  3.9679332 3.419840e-01
4 valley  5.9665704 1.297371e-06
5   peak  7.9866983 3.426214e-01
> plot(d, lwd=2, main="c(x1b,x2b,x3b)")
> abline(v=fit$x)
> 
> proc.time()
   user  system elapsed 
   0.35    0.04    0.39 

R version 4.0.5 (2021-03-31) -- "Shake and Throw"
Copyright (C) 2021 The R Foundation for Statistical Computing
Platform: i386-w64-mingw32/i386 (32-bit)

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.

> library("aroma.light")
aroma.light v3.20.0 successfully loaded. See ?aroma.light for help.
> 
> # Simulate data from the model y <- a + bx + x^c + eps(bx)
> J <- 1000
> x <- rexp(J)
> a <- c(2,15,3)
> b <- c(2,3,4)
> c <- c(1,2,1/2)
> bx <- outer(b,x)
> xc <- t(sapply(c, FUN=function(c) x^c))
> eps <- apply(bx, MARGIN=2, FUN=function(x) rnorm(length(b), mean=0, sd=0.1*x))
> y <- a + bx + xc + eps
> y <- t(y)
> 
> # Fit principal curve through (y_1, y_2, y_3)
> fit <- fitPrincipalCurve(y, verbose=TRUE)
Fitting principal curve...
 Data size: 1000x3
 Identifying missing values...
 Identifying missing values...done
 Data size after removing non-finite data points: 1000x3
 Calling principal_curve()...
Starting curve---distance^2: 1514522
Iteration 1---distance^2: 400.6155
Iteration 2---distance^2: 400.1743
Iteration 3---distance^2: 400.1808
  Converged: TRUE
  Number of iterations: 3
  Processing time/iteration: 0.1s (0.0s/iteration)
 Calling principal_curve()...done
Fitting principal curve...done
> 
> # Flip direction of 'lambda'?
> rho <- cor(fit$lambda, y[,1], use="complete.obs")
> flip <- (rho < 0)
> if (flip) {
+   fit$lambda <- max(fit$lambda, na.rm=TRUE)-fit$lambda
+ }
> 
> 
> # Backtransform (y_1, y_2, y_3) to be proportional to each other
> yN <- backtransformPrincipalCurve(y, fit=fit)
> 
> # Same backtransformation dimension by dimension
> yN2 <- y
> for (cc in 1:ncol(y)) {
+   yN2[,cc] <- backtransformPrincipalCurve(y, fit=fit, dimensions=cc)
+ }
> stopifnot(identical(yN2, yN))
> 
> 
> xlim <- c(0, 1.04*max(x))
> ylim <- range(c(y,yN), na.rm=TRUE)
> 
> 
> # Pairwise signals vs x before and after transform
> layout(matrix(1:4, nrow=2, byrow=TRUE))
> par(mar=c(4,4,3,2)+0.1)
> for (cc in 1:3) {
+   ylab <- substitute(y[c], env=list(c=cc))
+   plot(NA, xlim=xlim, ylim=ylim, xlab="x", ylab=ylab)
+   abline(h=a[cc], lty=3)
+   mtext(side=4, at=a[cc], sprintf("a=%g", a[cc]),
+         cex=0.8, las=2, line=0, adj=1.1, padj=-0.2)
+   points(x, y[,cc])
+   points(x, yN[,cc], col="tomato")
+   legend("topleft", col=c("black", "tomato"), pch=19,
+                     c("orignal", "transformed"), bty="n")
+ }
> title(main="Pairwise signals vs x before and after transform", outer=TRUE, line=-2)
> 
> 
> # Pairwise signals before and after transform
> layout(matrix(1:4, nrow=2, byrow=TRUE))
> par(mar=c(4,4,3,2)+0.1)
> for (rr in 3:2) {
+   ylab <- substitute(y[c], env=list(c=rr))
+   for (cc in 1:2) {
+     if (cc == rr) {
+       plot.new()
+       next
+     }
+     xlab <- substitute(y[c], env=list(c=cc))
+     plot(NA, xlim=ylim, ylim=ylim, xlab=xlab, ylab=ylab)
+     abline(a=0, b=1, lty=2)
+     points(y[,c(cc,rr)])
+     points(yN[,c(cc,rr)], col="tomato")
+     legend("topleft", col=c("black", "tomato"), pch=19,
+                       c("orignal", "transformed"), bty="n")
+   }
+ }
> title(main="Pairwise signals before and after transform", outer=TRUE, line=-2)
> 
> proc.time()
   user  system elapsed 
   1.14    0.06    1.18 

R version 4.0.5 (2021-03-31) -- "Shake and Throw"
Copyright (C) 2021 The R Foundation for Statistical Computing
Platform: x86_64-w64-mingw32/x64 (64-bit)

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.

> library("aroma.light")
aroma.light v3.20.0 successfully loaded. See ?aroma.light for help.
> 
> # Simulate data from the model y <- a + bx + x^c + eps(bx)
> J <- 1000
> x <- rexp(J)
> a <- c(2,15,3)
> b <- c(2,3,4)
> c <- c(1,2,1/2)
> bx <- outer(b,x)
> xc <- t(sapply(c, FUN=function(c) x^c))
> eps <- apply(bx, MARGIN=2, FUN=function(x) rnorm(length(b), mean=0, sd=0.1*x))
> y <- a + bx + xc + eps
> y <- t(y)
> 
> # Fit principal curve through (y_1, y_2, y_3)
> fit <- fitPrincipalCurve(y, verbose=TRUE)
Fitting principal curve...
 Data size: 1000x3
 Identifying missing values...
 Identifying missing values...done
 Data size after removing non-finite data points: 1000x3
 Calling principal_curve()...
Starting curve---distance^2: 1800863
Iteration 1---distance^2: 398.5954
Iteration 2---distance^2: 397.9932
Iteration 3---distance^2: 398.011
  Converged: TRUE
  Number of iterations: 3
  Processing time/iteration: 0.1s (0.0s/iteration)
 Calling principal_curve()...done
Fitting principal curve...done
> 
> # Flip direction of 'lambda'?
> rho <- cor(fit$lambda, y[,1], use="complete.obs")
> flip <- (rho < 0)
> if (flip) {
+   fit$lambda <- max(fit$lambda, na.rm=TRUE)-fit$lambda
+ }
> 
> 
> # Backtransform (y_1, y_2, y_3) to be proportional to each other
> yN <- backtransformPrincipalCurve(y, fit=fit)
> 
> # Same backtransformation dimension by dimension
> yN2 <- y
> for (cc in 1:ncol(y)) {
+   yN2[,cc] <- backtransformPrincipalCurve(y, fit=fit, dimensions=cc)
+ }
> stopifnot(identical(yN2, yN))
> 
> 
> xlim <- c(0, 1.04*max(x))
> ylim <- range(c(y,yN), na.rm=TRUE)
> 
> 
> # Pairwise signals vs x before and after transform
> layout(matrix(1:4, nrow=2, byrow=TRUE))
> par(mar=c(4,4,3,2)+0.1)
> for (cc in 1:3) {
+   ylab <- substitute(y[c], env=list(c=cc))
+   plot(NA, xlim=xlim, ylim=ylim, xlab="x", ylab=ylab)
+   abline(h=a[cc], lty=3)
+   mtext(side=4, at=a[cc], sprintf("a=%g", a[cc]),
+         cex=0.8, las=2, line=0, adj=1.1, padj=-0.2)
+   points(x, y[,cc])
+   points(x, yN[,cc], col="tomato")
+   legend("topleft", col=c("black", "tomato"), pch=19,
+                     c("orignal", "transformed"), bty="n")
+ }
> title(main="Pairwise signals vs x before and after transform", outer=TRUE, line=-2)
> 
> 
> # Pairwise signals before and after transform
> layout(matrix(1:4, nrow=2, byrow=TRUE))
> par(mar=c(4,4,3,2)+0.1)
> for (rr in 3:2) {
+   ylab <- substitute(y[c], env=list(c=rr))
+   for (cc in 1:2) {
+     if (cc == rr) {
+       plot.new()
+       next
+     }
+     xlab <- substitute(y[c], env=list(c=cc))
+     plot(NA, xlim=ylim, ylim=ylim, xlab=xlab, ylab=ylab)
+     abline(a=0, b=1, lty=2)
+     points(y[,c(cc,rr)])
+     points(yN[,c(cc,rr)], col="tomato")
+     legend("topleft", col=c("black", "tomato"), pch=19,
+                       c("orignal", "transformed"), bty="n")
+   }
+ }
> title(main="Pairwise signals before and after transform", outer=TRUE, line=-2)
> 
> proc.time()
   user  system elapsed 
   1.01    0.14    1.14 

R version 4.0.5 (2021-03-31) -- "Shake and Throw"
Copyright (C) 2021 The R Foundation for Statistical Computing
Platform: i386-w64-mingw32/i386 (32-bit)

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.

> library("aroma.light")
aroma.light v3.20.0 successfully loaded. See ?aroma.light for help.
> 
> # Simulate data from the model y <- a + bx + x^c + eps(bx)
> x <- rexp(1000)
> a <- c(2,15)
> b <- c(2,1)
> c <- c(1,2)
> bx <- outer(b,x)
> xc <- t(sapply(c, FUN=function(c) x^c))
> eps <- apply(bx, MARGIN=2, FUN=function(x) rnorm(length(x), mean=0, sd=0.1*x))
> Y <- a + bx + xc + eps
> Y <- t(Y)
> 
> lim <- c(0,70)
> plot(Y, xlim=lim, ylim=lim)
> 
> # Fit principal curve through a subset of (y_1, y_2)
> subset <- sample(nrow(Y), size=0.3*nrow(Y))
> fit <- fitXYCurve(Y[subset,], bandwidth=0.2)
> 
> lines(fit, col="red", lwd=2)
> 
> # Backtransform (y_1, y_2) keeping y_1 unchanged
> YN <- backtransformXYCurve(Y, fit=fit)
> points(YN, col="blue")
> abline(a=0, b=1, col="red", lwd=2)
> 
> proc.time()
   user  system elapsed 
   0.34    0.12    0.45 

R version 4.0.5 (2021-03-31) -- "Shake and Throw"
Copyright (C) 2021 The R Foundation for Statistical Computing
Platform: x86_64-w64-mingw32/x64 (64-bit)

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.

> library("aroma.light")
aroma.light v3.20.0 successfully loaded. See ?aroma.light for help.
> 
> # Simulate data from the model y <- a + bx + x^c + eps(bx)
> x <- rexp(1000)
> a <- c(2,15)
> b <- c(2,1)
> c <- c(1,2)
> bx <- outer(b,x)
> xc <- t(sapply(c, FUN=function(c) x^c))
> eps <- apply(bx, MARGIN=2, FUN=function(x) rnorm(length(x), mean=0, sd=0.1*x))
> Y <- a + bx + xc + eps
> Y <- t(Y)
> 
> lim <- c(0,70)
> plot(Y, xlim=lim, ylim=lim)
> 
> # Fit principal curve through a subset of (y_1, y_2)
> subset <- sample(nrow(Y), size=0.3*nrow(Y))
> fit <- fitXYCurve(Y[subset,], bandwidth=0.2)
> 
> lines(fit, col="red", lwd=2)
> 
> # Backtransform (y_1, y_2) keeping y_1 unchanged
> YN <- backtransformXYCurve(Y, fit=fit)
> points(YN, col="blue")
> abline(a=0, b=1, col="red", lwd=2)
> 
> proc.time()
   user  system elapsed 
   0.42    0.09    0.50 

R version 4.0.5 (2021-03-31) -- "Shake and Throw"
Copyright (C) 2021 The R Foundation for Statistical Computing
Platform: i386-w64-mingw32/i386 (32-bit)

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.

> library("aroma.light")
aroma.light v3.20.0 successfully loaded. See ?aroma.light for help.
> 
> for (zzz in 0) {
+ 
+ # This example requires plot3d() in R.basic [http://www.braju.com/R/]
+ if (!require(pkgName <- "R.basic", character.only=TRUE)) break
+ 
+ # Simulate data from the model y <- a + bx + eps(bx)
+ x <- rexp(1000)
+ a <- c(2,15,3)
+ b <- c(2,3,4)
+ bx <- outer(b,x)
+ eps <- apply(bx, MARGIN=2, FUN=function(x) rnorm(length(x), mean=0, sd=0.1*x))
+ y <- a + bx + eps
+ y <- t(y)
+ 
+ # Add some outliers by permuting the dimensions for 1/10 of the observations
+ idx <- sample(1:nrow(y), size=1/10*nrow(y))
+ y[idx,] <- y[idx,c(2,3,1)]
+ 
+ # Plot the data with fitted lines at four different view points
+ opar <- par(mar=c(1,1,1,1)+0.1)
+ N <- 4
+ layout(matrix(1:N, nrow=2, byrow=TRUE))
+ theta <- seq(0,270,length.out=N)
+ phi <- rep(20, length.out=N)
+ xlim <- ylim <- zlim <- c(0,45)
+ persp <- list()
+ for (kk in seq_along(theta)) {
+   # Plot the data
+   persp[[kk]] <- plot3d(y, theta=theta[kk], phi=phi[kk], xlim=xlim, ylim=ylim, zlim=zlim)
+ }
+ 
+ # Weights on the observations
+ # Example a: Equal weights
+ w <- NULL
+ # Example b: More weight on the outliers (uncomment to test)
+ w <- rep(1, length(x)); w[idx] <- 0.8
+ 
+ # ...and show all iterations too with different colors.
+ maxIter <- c(seq(1,20,length.out=10),Inf)
+ col <- topo.colors(length(maxIter))
+ # Show the fitted value for every iteration
+ for (ii in seq_along(maxIter)) {
+   # Fit a line using IWPCA through data
+   fit <- iwpca(y, w=w, maxIter=maxIter[ii], swapDirections=TRUE)
+ 
+   ymid <- fit$xMean
+   d0 <- apply(y, MARGIN=2, FUN=min) - ymid
+   d1 <- apply(y, MARGIN=2, FUN=max) - ymid
+   b <- fit$vt[1,]
+   y0 <- -b * max(abs(d0))
+   y1 <-  b * max(abs(d1))
+   yline <- matrix(c(y0,y1), nrow=length(b), ncol=2)
+   yline <- yline + ymid
+ 
+   for (kk in seq_along(theta)) {
+     # Set pane to draw in
+     par(mfg=c((kk-1) %/% 2, (kk-1) %% 2) + 1)
+     # Set the viewpoint of the pane
+     options(persp.matrix=persp[[kk]])
+ 
+     # Get the first principal component
+     points3d(t(ymid), col=col[ii])
+     lines3d(t(yline), col=col[ii])
+ 
+     # Highlight the last one
+     if (ii == length(maxIter))
+       lines3d(t(yline), col="red", lwd=3)
+   }
+ }
+ 
+ par(opar)
+ 
+ } # for (zzz in 0)
Loading required package: R.basic
Warning message:
In library(package, lib.loc = lib.loc, character.only = TRUE, logical.return = TRUE,  :
  there is no package called 'R.basic'
> rm(zzz)
> 
> proc.time()
   user  system elapsed 
   0.29    0.07    0.35 

R version 4.0.5 (2021-03-31) -- "Shake and Throw"
Copyright (C) 2021 The R Foundation for Statistical Computing
Platform: x86_64-w64-mingw32/x64 (64-bit)

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.

> library("aroma.light")
aroma.light v3.20.0 successfully loaded. See ?aroma.light for help.
> 
> for (zzz in 0) {
+ 
+ # This example requires plot3d() in R.basic [http://www.braju.com/R/]
+ if (!require(pkgName <- "R.basic", character.only=TRUE)) break
+ 
+ # Simulate data from the model y <- a + bx + eps(bx)
+ x <- rexp(1000)
+ a <- c(2,15,3)
+ b <- c(2,3,4)
+ bx <- outer(b,x)
+ eps <- apply(bx, MARGIN=2, FUN=function(x) rnorm(length(x), mean=0, sd=0.1*x))
+ y <- a + bx + eps
+ y <- t(y)
+ 
+ # Add some outliers by permuting the dimensions for 1/10 of the observations
+ idx <- sample(1:nrow(y), size=1/10*nrow(y))
+ y[idx,] <- y[idx,c(2,3,1)]
+ 
+ # Plot the data with fitted lines at four different view points
+ opar <- par(mar=c(1,1,1,1)+0.1)
+ N <- 4
+ layout(matrix(1:N, nrow=2, byrow=TRUE))
+ theta <- seq(0,270,length.out=N)
+ phi <- rep(20, length.out=N)
+ xlim <- ylim <- zlim <- c(0,45)
+ persp <- list()
+ for (kk in seq_along(theta)) {
+   # Plot the data
+   persp[[kk]] <- plot3d(y, theta=theta[kk], phi=phi[kk], xlim=xlim, ylim=ylim, zlim=zlim)
+ }
+ 
+ # Weights on the observations
+ # Example a: Equal weights
+ w <- NULL
+ # Example b: More weight on the outliers (uncomment to test)
+ w <- rep(1, length(x)); w[idx] <- 0.8
+ 
+ # ...and show all iterations too with different colors.
+ maxIter <- c(seq(1,20,length.out=10),Inf)
+ col <- topo.colors(length(maxIter))
+ # Show the fitted value for every iteration
+ for (ii in seq_along(maxIter)) {
+   # Fit a line using IWPCA through data
+   fit <- iwpca(y, w=w, maxIter=maxIter[ii], swapDirections=TRUE)
+ 
+   ymid <- fit$xMean
+   d0 <- apply(y, MARGIN=2, FUN=min) - ymid
+   d1 <- apply(y, MARGIN=2, FUN=max) - ymid
+   b <- fit$vt[1,]
+   y0 <- -b * max(abs(d0))
+   y1 <-  b * max(abs(d1))
+   yline <- matrix(c(y0,y1), nrow=length(b), ncol=2)
+   yline <- yline + ymid
+ 
+   for (kk in seq_along(theta)) {
+     # Set pane to draw in
+     par(mfg=c((kk-1) %/% 2, (kk-1) %% 2) + 1)
+     # Set the viewpoint of the pane
+     options(persp.matrix=persp[[kk]])
+ 
+     # Get the first principal component
+     points3d(t(ymid), col=col[ii])
+     lines3d(t(yline), col=col[ii])
+ 
+     # Highlight the last one
+     if (ii == length(maxIter))
+       lines3d(t(yline), col="red", lwd=3)
+   }
+ }
+ 
+ par(opar)
+ 
+ } # for (zzz in 0)
Loading required package: R.basic
Warning message:
In library(package, lib.loc = lib.loc, character.only = TRUE, logical.return = TRUE,  :
  there is no package called 'R.basic'
> rm(zzz)
> 
> proc.time()
   user  system elapsed 
   0.34    0.03    0.37 

R version 4.0.5 (2021-03-31) -- "Shake and Throw"
Copyright (C) 2021 The R Foundation for Statistical Computing
Platform: i386-w64-mingw32/i386 (32-bit)

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.

> library("aroma.light")
aroma.light v3.20.0 successfully loaded. See ?aroma.light for help.
> 
> # Define f(x)
> f <- expression(0.1*x^4 + 1*x^3 + 2*x^2 + x + 10*sin(2*x))
> 
> # Simulate data from this function in the range [a,b]
> a <- -2; b <- 5
> x <- seq(a, b, length.out=3000)
> y <- eval(f)
> 
> # Add some noise to the data
> y <- y + rnorm(length(y), 0, 10)
> 
> # Plot the function and its second derivative
> plot(x,y, type="l", lwd=4)
> 
> # Fit a cubic smoothing spline and plot it
> g <- smooth.spline(x,y, df=16)
> lines(g, col="yellow", lwd=2, lty=2)
> 
> # Calculating the (log) likelihood of the fitted spline
> l <- likelihood(g)
> 
> cat("Log likelihood with unique x values:\n")
Log likelihood with unique x values:
> print(l)
Likelihood of smoothing spline: -298948.9 
 Log base: 2.718282 
 Weighted residuals sum of square: 298949.1 
 Penalty: -0.1164618 
 Smoothing parameter lambda: 0.0009257147 
 Roughness score: 125.8075 
> 
> # Note that this is not the same as the log likelihood of the
> # data on the fitted spline iff the x values are non-unique
> x[1:5] <- x[1]  # Non-unique x values
> g <- smooth.spline(x,y, df=16)
> l <- likelihood(g)
> 
> cat("\nLog likelihood of the *spline* data set:\n")

Log likelihood of the *spline* data set:
> print(l)
Likelihood of smoothing spline: -298228.6 
 Log base: 2.718282 
 Weighted residuals sum of square: 298228.7 
 Penalty: -0.1163556 
 Smoothing parameter lambda: 0.0009261969 
 Roughness score: 125.6273 
> 
> # In cases with non unique x values one has to proceed as
> # below if one want to get the log likelihood for the original
> # data.
> l <- likelihood(g, x=x, y=y)
> cat("\nLog likelihood of the *original* data set:\n")

Log likelihood of the *original* data set:
> print(l)
Likelihood of smoothing spline: -298954.7 
 Log base: 2.718282 
 Weighted residuals sum of square: 298954.8 
 Penalty: -0.1163556 
 Smoothing parameter lambda: 0.0009261969 
 Roughness score: 125.6272 
> 
> 
> 
> 
> 
> 
> proc.time()
   user  system elapsed 
   0.67    0.04    0.70 

R version 4.0.5 (2021-03-31) -- "Shake and Throw"
Copyright (C) 2021 The R Foundation for Statistical Computing
Platform: x86_64-w64-mingw32/x64 (64-bit)

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.

> library("aroma.light")
aroma.light v3.20.0 successfully loaded. See ?aroma.light for help.
> 
> # Define f(x)
> f <- expression(0.1*x^4 + 1*x^3 + 2*x^2 + x + 10*sin(2*x))
> 
> # Simulate data from this function in the range [a,b]
> a <- -2; b <- 5
> x <- seq(a, b, length.out=3000)
> y <- eval(f)
> 
> # Add some noise to the data
> y <- y + rnorm(length(y), 0, 10)
> 
> # Plot the function and its second derivative
> plot(x,y, type="l", lwd=4)
> 
> # Fit a cubic smoothing spline and plot it
> g <- smooth.spline(x,y, df=16)
> lines(g, col="yellow", lwd=2, lty=2)
> 
> # Calculating the (log) likelihood of the fitted spline
> l <- likelihood(g)
> 
> cat("Log likelihood with unique x values:\n")
Log likelihood with unique x values:
> print(l)
Likelihood of smoothing spline: -312905.5 
 Log base: 2.718282 
 Weighted residuals sum of square: 312905.6 
 Penalty: -0.1138784 
 Smoothing parameter lambda: 0.0009257147 
 Roughness score: 123.0167 
> 
> # Note that this is not the same as the log likelihood of the
> # data on the fitted spline iff the x values are non-unique
> x[1:5] <- x[1]  # Non-unique x values
> g <- smooth.spline(x,y, df=16)
> l <- likelihood(g)
> 
> cat("\nLog likelihood of the *spline* data set:\n")

Log likelihood of the *spline* data set:
> print(l)
Likelihood of smoothing spline: -311635.5 
 Log base: 2.718282 
 Weighted residuals sum of square: 311635.6 
 Penalty: -0.1138943 
 Smoothing parameter lambda: 0.0009261969 
 Roughness score: 122.9698 
> 
> # In cases with non unique x values one has to proceed as
> # below if one want to get the log likelihood for the original
> # data.
> l <- likelihood(g, x=x, y=y)
> cat("\nLog likelihood of the *original* data set:\n")

Log likelihood of the *original* data set:
> print(l)
Likelihood of smoothing spline: -312906.3 
 Log base: 2.718282 
 Weighted residuals sum of square: 312906.4 
 Penalty: -0.1138942 
 Smoothing parameter lambda: 0.0009261969 
 Roughness score: 122.9697 
> 
> 
> 
> 
> 
> 
> proc.time()
   user  system elapsed 
   0.42    0.09    0.50 

R version 4.0.5 (2021-03-31) -- "Shake and Throw"
Copyright (C) 2021 The R Foundation for Statistical Computing
Platform: i386-w64-mingw32/i386 (32-bit)

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.

> library("aroma.light")
aroma.light v3.20.0 successfully loaded. See ?aroma.light for help.
> 
> # Deaths from sport parachuting;  from ABC of EDA, p.224:
> deaths <- matrix(c(14,15,14, 7,4,7, 8,2,10, 15,9,10, 0,2,0), ncol=3, byrow=TRUE)
> rownames(deaths) <- c("1-24", "25-74", "75-199", "200++", "NA")
> colnames(deaths) <- 1973:1975
> 
> print(deaths)
       1973 1974 1975
1-24     14   15   14
25-74     7    4    7
75-199    8    2   10
200++    15    9   10
NA        0    2    0
> 
> mp <- medianPolish(deaths)
> mp1 <- medpolish(deaths, trace=FALSE)
> print(mp)

Median Polish Results (Dataset: "deaths")

Overall: 8

Row Effects:
  1-24  25-74 75-199  200++     NA 
     6     -1      0      2     -8 

Column Effects:
1973 1974 1975 
   0   -1    0 

Residuals:
       1973 1974 1975
1-24      0    2    0
25-74     0   -2    0
75-199    0   -5    2
200++     5    0    0
NA        0    3    0

> 
> ff <- c("overall", "row", "col", "residuals")
> stopifnot(all.equal(mp[ff], mp1[ff]))
> 
> # Validate decomposition:
> stopifnot(all.equal(deaths, mp$overall+outer(mp$row,mp$col,"+")+mp$resid))
> 
> proc.time()
   user  system elapsed 
   0.43    0.03    0.45 

R version 4.0.5 (2021-03-31) -- "Shake and Throw"
Copyright (C) 2021 The R Foundation for Statistical Computing
Platform: x86_64-w64-mingw32/x64 (64-bit)

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.

> library("aroma.light")
aroma.light v3.20.0 successfully loaded. See ?aroma.light for help.
> 
> # Deaths from sport parachuting;  from ABC of EDA, p.224:
> deaths <- matrix(c(14,15,14, 7,4,7, 8,2,10, 15,9,10, 0,2,0), ncol=3, byrow=TRUE)
> rownames(deaths) <- c("1-24", "25-74", "75-199", "200++", "NA")
> colnames(deaths) <- 1973:1975
> 
> print(deaths)
       1973 1974 1975
1-24     14   15   14
25-74     7    4    7
75-199    8    2   10
200++    15    9   10
NA        0    2    0
> 
> mp <- medianPolish(deaths)
> mp1 <- medpolish(deaths, trace=FALSE)
> print(mp)

Median Polish Results (Dataset: "deaths")

Overall: 8

Row Effects:
  1-24  25-74 75-199  200++     NA 
     6     -1      0      2     -8 

Column Effects:
1973 1974 1975 
   0   -1    0 

Residuals:
       1973 1974 1975
1-24      0    2    0
25-74     0   -2    0
75-199    0   -5    2
200++     5    0    0
NA        0    3    0

> 
> ff <- c("overall", "row", "col", "residuals")
> stopifnot(all.equal(mp[ff], mp1[ff]))
> 
> # Validate decomposition:
> stopifnot(all.equal(deaths, mp$overall+outer(mp$row,mp$col,"+")+mp$resid))
> 
> proc.time()
   user  system elapsed 
   0.23    0.06    0.29 

R version 4.0.5 (2021-03-31) -- "Shake and Throw"
Copyright (C) 2021 The R Foundation for Statistical Computing
Platform: i386-w64-mingw32/i386 (32-bit)

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.

> library("aroma.light")
aroma.light v3.20.0 successfully loaded. See ?aroma.light for help.
> 
> pathname <- system.file("data-ex", "PMT-RGData.dat", package="aroma.light")
> rg <- read.table(pathname, header=TRUE, sep="\t")
> nbrOfScans <- max(rg$slide)
> 
> rg <- as.list(rg)
> for (field in c("R", "G"))
+   rg[[field]] <- matrix(as.double(rg[[field]]), ncol=nbrOfScans)
> rg$slide <- rg$spot <- NULL
> rg <- as.matrix(as.data.frame(rg))
> colnames(rg) <- rep(c("R", "G"), each=nbrOfScans)
> 
> rgC <- rg
> 
> layout(matrix(c(1,2,0,3,4,0,5,6,7), ncol=3, byrow=TRUE))
> 
> for (channel in c("R", "G")) {
+   sidx <- which(colnames(rg) == channel)
+   channelColor <- switch(channel, R="red", G="green")
+ 
+   # - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
+   # The raw data
+   # - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
+   plotMvsAPairs(rg, channel=channel)
+   title(main=paste("Observed", channel))
+   box(col=channelColor)
+ 
+   # - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
+   # The calibrated data
+   # - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
+   rgC[,sidx] <- calibrateMultiscan(rg[,sidx], average=NULL)
+ 
+   plotMvsAPairs(rgC, channel=channel)
+   title(main=paste("Calibrated", channel))
+   box(col=channelColor)
+ } # for (channel ...)
There were 50 or more warnings (use warnings() to see the first 50)
> 
> 
> # - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
> # The average calibrated data
> #
> # Note how the red signals are weaker than the green. The reason
> # for this can be that the scale factor in the green channel is
> # greater than in the red channel, but it can also be that there
> # is a remaining relative difference in bias between the green
> # and the red channel, a bias that precedes the scanning.
> # - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
> rgCA <- matrix(NA_real_, nrow=nrow(rg), ncol=2)
> colnames(rgCA) <- c("R", "G")
> for (channel in c("R", "G")) {
+   sidx <- which(colnames(rg) == channel)
+   rgCA[,channel] <- calibrateMultiscan(rg[,sidx])
+ }
> 
> plotMvsA(rgCA)
> title(main="Average calibrated")
> 
> 
> # - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
> # The affine normalized average calibrated data
> # - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
> # Create a matrix where the columns represent the channels
> # to be normalized.
> rgCAN <- rgCA
> # Affine normalization of channels
> rgCAN <- normalizeAffine(rgCAN)
> 
> plotMvsA(rgCAN)
> title(main="Affine normalized A.C.")
> 
> 
> # - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
> # It is always ok to rescale the affine normalized data if its
> # done on (R,G); not on (A,M)! However, this is only needed for
> # esthetic purposes.
> # - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
> rgCAN <- rgCAN * 2^5
> plotMvsA(rgCAN)
> title(main="Rescaled normalized")
> 
> 
> 
> proc.time()
   user  system elapsed 
   4.00    0.12    4.11 

R version 4.0.5 (2021-03-31) -- "Shake and Throw"
Copyright (C) 2021 The R Foundation for Statistical Computing
Platform: x86_64-w64-mingw32/x64 (64-bit)

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.

> library("aroma.light")
aroma.light v3.20.0 successfully loaded. See ?aroma.light for help.
> 
> pathname <- system.file("data-ex", "PMT-RGData.dat", package="aroma.light")
> rg <- read.table(pathname, header=TRUE, sep="\t")
> nbrOfScans <- max(rg$slide)
> 
> rg <- as.list(rg)
> for (field in c("R", "G"))
+   rg[[field]] <- matrix(as.double(rg[[field]]), ncol=nbrOfScans)
> rg$slide <- rg$spot <- NULL
> rg <- as.matrix(as.data.frame(rg))
> colnames(rg) <- rep(c("R", "G"), each=nbrOfScans)
> 
> rgC <- rg
> 
> layout(matrix(c(1,2,0,3,4,0,5,6,7), ncol=3, byrow=TRUE))
> 
> for (channel in c("R", "G")) {
+   sidx <- which(colnames(rg) == channel)
+   channelColor <- switch(channel, R="red", G="green")
+ 
+   # - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
+   # The raw data
+   # - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
+   plotMvsAPairs(rg, channel=channel)
+   title(main=paste("Observed", channel))
+   box(col=channelColor)
+ 
+   # - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
+   # The calibrated data
+   # - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
+   rgC[,sidx] <- calibrateMultiscan(rg[,sidx], average=NULL)
+ 
+   plotMvsAPairs(rgC, channel=channel)
+   title(main=paste("Calibrated", channel))
+   box(col=channelColor)
+ } # for (channel ...)
There were 50 or more warnings (use warnings() to see the first 50)
> 
> 
> # - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
> # The average calibrated data
> #
> # Note how the red signals are weaker than the green. The reason
> # for this can be that the scale factor in the green channel is
> # greater than in the red channel, but it can also be that there
> # is a remaining relative difference in bias between the green
> # and the red channel, a bias that precedes the scanning.
> # - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
> rgCA <- matrix(NA_real_, nrow=nrow(rg), ncol=2)
> colnames(rgCA) <- c("R", "G")
> for (channel in c("R", "G")) {
+   sidx <- which(colnames(rg) == channel)
+   rgCA[,channel] <- calibrateMultiscan(rg[,sidx])
+ }
> 
> plotMvsA(rgCA)
> title(main="Average calibrated")
> 
> 
> # - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
> # The affine normalized average calibrated data
> # - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
> # Create a matrix where the columns represent the channels
> # to be normalized.
> rgCAN <- rgCA
> # Affine normalization of channels
> rgCAN <- normalizeAffine(rgCAN)
> 
> plotMvsA(rgCAN)
> title(main="Affine normalized A.C.")
> 
> 
> # - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
> # It is always ok to rescale the affine normalized data if its
> # done on (R,G); not on (A,M)! However, this is only needed for
> # esthetic purposes.
> # - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
> rgCAN <- rgCAN * 2^5
> plotMvsA(rgCAN)
> title(main="Rescaled normalized")
> 
> 
> 
> proc.time()
   user  system elapsed 
   3.25    0.10    3.34 

R version 4.0.5 (2021-03-31) -- "Shake and Throw"
Copyright (C) 2021 The R Foundation for Statistical Computing
Platform: i386-w64-mingw32/i386 (32-bit)

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.

> library("aroma.light")
aroma.light v3.20.0 successfully loaded. See ?aroma.light for help.
> 
> # Simulate ten samples of different lengths
> N <- 10000
> X <- list()
> for (kk in 1:8) {
+   rfcn <- list(rnorm, rgamma)[[sample(2, size=1)]]
+   size <- runif(1, min=0.3, max=1)
+   a <- rgamma(1, shape=20, rate=10)
+   b <- rgamma(1, shape=10, rate=10)
+   values <- rfcn(size*N, a, b)
+ 
+   # "Censor" values
+   values[values < 0 | values > 8] <- NA_real_
+ 
+   X[[kk]] <- values
+ }
> 
> # Add 20% missing values
> X <- lapply(X, FUN=function(x) {
+   x[sample(length(x), size=0.20*length(x))] <- NA_real_
+   x
+ })
> 
> # Normalize quantiles
> Xn <- normalizeAverage(X, na.rm=TRUE, targetAvg=median(unlist(X), na.rm=TRUE))
> 
> # Plot the data
> layout(matrix(1:2, ncol=1))
> xlim <- range(X, Xn, na.rm=TRUE)
> plotDensity(X, lwd=2, xlim=xlim, main="The original distributions")
> plotDensity(Xn, lwd=2, xlim=xlim, main="The normalized distributions")
> 
> proc.time()
   user  system elapsed 
   0.37    0.09    0.45 

R version 4.0.5 (2021-03-31) -- "Shake and Throw"
Copyright (C) 2021 The R Foundation for Statistical Computing
Platform: x86_64-w64-mingw32/x64 (64-bit)

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.

> library("aroma.light")
aroma.light v3.20.0 successfully loaded. See ?aroma.light for help.
> 
> # Simulate ten samples of different lengths
> N <- 10000
> X <- list()
> for (kk in 1:8) {
+   rfcn <- list(rnorm, rgamma)[[sample(2, size=1)]]
+   size <- runif(1, min=0.3, max=1)
+   a <- rgamma(1, shape=20, rate=10)
+   b <- rgamma(1, shape=10, rate=10)
+   values <- rfcn(size*N, a, b)
+ 
+   # "Censor" values
+   values[values < 0 | values > 8] <- NA_real_
+ 
+   X[[kk]] <- values
+ }
> 
> # Add 20% missing values
> X <- lapply(X, FUN=function(x) {
+   x[sample(length(x), size=0.20*length(x))] <- NA_real_
+   x
+ })
> 
> # Normalize quantiles
> Xn <- normalizeAverage(X, na.rm=TRUE, targetAvg=median(unlist(X), na.rm=TRUE))
> 
> # Plot the data
> layout(matrix(1:2, ncol=1))
> xlim <- range(X, Xn, na.rm=TRUE)
> plotDensity(X, lwd=2, xlim=xlim, main="The original distributions")
> plotDensity(Xn, lwd=2, xlim=xlim, main="The normalized distributions")
> 
> proc.time()
   user  system elapsed 
   0.42    0.03    0.45 

R version 4.0.5 (2021-03-31) -- "Shake and Throw"
Copyright (C) 2021 The R Foundation for Statistical Computing
Platform: i386-w64-mingw32/i386 (32-bit)

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.

> library("aroma.light")
aroma.light v3.20.0 successfully loaded. See ?aroma.light for help.
> 
> # Simulate three samples with on average 20% missing values
> N <- 10000
> X <- cbind(rnorm(N, mean=3, sd=1),
+            rnorm(N, mean=4, sd=2),
+            rgamma(N, shape=2, rate=1))
> X[sample(3*N, size=0.20*3*N)] <- NA_real_
> 
> # Normalize quantiles
> Xn <- normalizeAverage(X, na.rm=TRUE, targetAvg=median(X, na.rm=TRUE))
> 
> # Plot the data
> layout(matrix(1:2, ncol=1))
> xlim <- range(X, Xn, na.rm=TRUE)
> plotDensity(X, lwd=2, xlim=xlim, main="The three original distributions")
> plotDensity(Xn, lwd=2, xlim=xlim, main="The three normalized distributions")
> 
> proc.time()
   user  system elapsed 
   0.37    0.03    0.39 

R version 4.0.5 (2021-03-31) -- "Shake and Throw"
Copyright (C) 2021 The R Foundation for Statistical Computing
Platform: x86_64-w64-mingw32/x64 (64-bit)

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.

> library("aroma.light")
aroma.light v3.20.0 successfully loaded. See ?aroma.light for help.
> 
> # Simulate three samples with on average 20% missing values
> N <- 10000
> X <- cbind(rnorm(N, mean=3, sd=1),
+            rnorm(N, mean=4, sd=2),
+            rgamma(N, shape=2, rate=1))
> X[sample(3*N, size=0.20*3*N)] <- NA_real_
> 
> # Normalize quantiles
> Xn <- normalizeAverage(X, na.rm=TRUE, targetAvg=median(X, na.rm=TRUE))
> 
> # Plot the data
> layout(matrix(1:2, ncol=1))
> xlim <- range(X, Xn, na.rm=TRUE)
> plotDensity(X, lwd=2, xlim=xlim, main="The three original distributions")
> plotDensity(Xn, lwd=2, xlim=xlim, main="The three normalized distributions")
> 
> proc.time()
   user  system elapsed 
   0.32    0.06    0.39 

R version 4.0.5 (2021-03-31) -- "Shake and Throw"
Copyright (C) 2021 The R Foundation for Statistical Computing
Platform: i386-w64-mingw32/i386 (32-bit)

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.

> library("aroma.light")
aroma.light v3.20.0 successfully loaded. See ?aroma.light for help.
> 
> pathname <- system.file("data-ex", "PMT-RGData.dat", package="aroma.light")
> rg <- read.table(pathname, header=TRUE, sep="\t")
> nbrOfScans <- max(rg$slide)
> 
> rg <- as.list(rg)
> for (field in c("R", "G"))
+   rg[[field]] <- matrix(as.double(rg[[field]]), ncol=nbrOfScans)
> rg$slide <- rg$spot <- NULL
> rg <- as.matrix(as.data.frame(rg))
> colnames(rg) <- rep(c("R", "G"), each=nbrOfScans)
> 
> layout(matrix(c(1,2,0,3,4,0,5,6,7), ncol=3, byrow=TRUE))
> 
> rgC <- rg
> for (channel in c("R", "G")) {
+   sidx <- which(colnames(rg) == channel)
+   channelColor <- switch(channel, R="red", G="green")
+ 
+   # - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
+   # The raw data
+   # - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
+   plotMvsAPairs(rg[,sidx])
+   title(main=paste("Observed", channel))
+   box(col=channelColor)
+ 
+   # - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
+   # The calibrated data
+   # - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
+   rgC[,sidx] <- calibrateMultiscan(rg[,sidx], average=NULL)
+ 
+   plotMvsAPairs(rgC[,sidx])
+   title(main=paste("Calibrated", channel))
+   box(col=channelColor)
+ } # for (channel ...)
> 
> 
> # - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
> # The average calibrated data
> #
> # Note how the red signals are weaker than the green. The reason
> # for this can be that the scale factor in the green channel is
> # greater than in the red channel, but it can also be that there
> # is a remaining relative difference in bias between the green
> # and the red channel, a bias that precedes the scanning.
> # - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
> rgCA <- rg
> for (channel in c("R", "G")) {
+   sidx <- which(colnames(rg) == channel)
+   rgCA[,sidx] <- calibrateMultiscan(rg[,sidx])
+ }
> 
> rgCAavg <- matrix(NA_real_, nrow=nrow(rgCA), ncol=2)
> colnames(rgCAavg) <- c("R", "G")
> for (channel in c("R", "G")) {
+   sidx <- which(colnames(rg) == channel)
+   rgCAavg[,channel] <- apply(rgCA[,sidx], MARGIN=1, FUN=median, na.rm=TRUE)
+ }
> 
> # Add some "fake" outliers
> outliers <- 1:600
> rgCAavg[outliers,"G"] <- 50000
> 
> plotMvsA(rgCAavg)
> title(main="Average calibrated (AC)")
> 
> # - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
> # Normalize data
> # - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
> # Weight-down outliers when normalizing
> weights <- rep(1, nrow(rgCAavg))
> weights[outliers] <- 0.001
> 
> # Affine normalization of channels
> rgCANa <- normalizeAffine(rgCAavg, weights=weights)
> # It is always ok to rescale the affine normalized data if its
> # done on (R,G); not on (A,M)! However, this is only needed for
> # esthetic purposes.
> rgCANa <- rgCANa *2^1.4
> plotMvsA(rgCANa)
> title(main="Normalized AC")
> 
> # Curve-fit (lowess) normalization
> rgCANlw <- normalizeLowess(rgCAavg, weights=weights)
Warning message:
In normalizeCurveFit.matrix(X, method = "lowess", ...) :
  Weights were rounded to {0,1} since 'lowess' normalization supports only zero-one weights.
> plotMvsA(rgCANlw, col="orange", add=TRUE)
> 
> # Curve-fit (loess) normalization
> rgCANl <- normalizeLoess(rgCAavg, weights=weights)
> plotMvsA(rgCANl, col="red", add=TRUE)
> 
> # Curve-fit (robust spline) normalization
> rgCANrs <- normalizeRobustSpline(rgCAavg, weights=weights)
> plotMvsA(rgCANrs, col="blue", add=TRUE)
> 
> legend(x=0,y=16, legend=c("affine", "lowess", "loess", "r. spline"), pch=19,
+        col=c("black", "orange", "red", "blue"), ncol=2, x.intersp=0.3, bty="n")
> 
> 
> plotMvsMPairs(cbind(rgCANa, rgCANlw), col="orange", xlab=expression(M[affine]))
> title(main="Normalized AC")
> plotMvsMPairs(cbind(rgCANa, rgCANl), col="red", add=TRUE)
> plotMvsMPairs(cbind(rgCANa, rgCANrs), col="blue", add=TRUE)
> abline(a=0, b=1, lty=2)
> legend(x=-6,y=6, legend=c("lowess", "loess", "r. spline"), pch=19,
+        col=c("orange", "red", "blue"), ncol=2, x.intersp=0.3, bty="n")
> 
> 
> proc.time()
   user  system elapsed 
   7.14    0.10    7.25 

R version 4.0.5 (2021-03-31) -- "Shake and Throw"
Copyright (C) 2021 The R Foundation for Statistical Computing
Platform: x86_64-w64-mingw32/x64 (64-bit)

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.

> library("aroma.light")
aroma.light v3.20.0 successfully loaded. See ?aroma.light for help.
> 
> pathname <- system.file("data-ex", "PMT-RGData.dat", package="aroma.light")
> rg <- read.table(pathname, header=TRUE, sep="\t")
> nbrOfScans <- max(rg$slide)
> 
> rg <- as.list(rg)
> for (field in c("R", "G"))
+   rg[[field]] <- matrix(as.double(rg[[field]]), ncol=nbrOfScans)
> rg$slide <- rg$spot <- NULL
> rg <- as.matrix(as.data.frame(rg))
> colnames(rg) <- rep(c("R", "G"), each=nbrOfScans)
> 
> layout(matrix(c(1,2,0,3,4,0,5,6,7), ncol=3, byrow=TRUE))
> 
> rgC <- rg
> for (channel in c("R", "G")) {
+   sidx <- which(colnames(rg) == channel)
+   channelColor <- switch(channel, R="red", G="green")
+ 
+   # - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
+   # The raw data
+   # - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
+   plotMvsAPairs(rg[,sidx])
+   title(main=paste("Observed", channel))
+   box(col=channelColor)
+ 
+   # - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
+   # The calibrated data
+   # - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
+   rgC[,sidx] <- calibrateMultiscan(rg[,sidx], average=NULL)
+ 
+   plotMvsAPairs(rgC[,sidx])
+   title(main=paste("Calibrated", channel))
+   box(col=channelColor)
+ } # for (channel ...)
> 
> 
> # - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
> # The average calibrated data
> #
> # Note how the red signals are weaker than the green. The reason
> # for this can be that the scale factor in the green channel is
> # greater than in the red channel, but it can also be that there
> # is a remaining relative difference in bias between the green
> # and the red channel, a bias that precedes the scanning.
> # - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
> rgCA <- rg
> for (channel in c("R", "G")) {
+   sidx <- which(colnames(rg) == channel)
+   rgCA[,sidx] <- calibrateMultiscan(rg[,sidx])
+ }
> 
> rgCAavg <- matrix(NA_real_, nrow=nrow(rgCA), ncol=2)
> colnames(rgCAavg) <- c("R", "G")
> for (channel in c("R", "G")) {
+   sidx <- which(colnames(rg) == channel)
+   rgCAavg[,channel] <- apply(rgCA[,sidx], MARGIN=1, FUN=median, na.rm=TRUE)
+ }
> 
> # Add some "fake" outliers
> outliers <- 1:600
> rgCAavg[outliers,"G"] <- 50000
> 
> plotMvsA(rgCAavg)
> title(main="Average calibrated (AC)")
> 
> # - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
> # Normalize data
> # - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
> # Weight-down outliers when normalizing
> weights <- rep(1, nrow(rgCAavg))
> weights[outliers] <- 0.001
> 
> # Affine normalization of channels
> rgCANa <- normalizeAffine(rgCAavg, weights=weights)
> # It is always ok to rescale the affine normalized data if its
> # done on (R,G); not on (A,M)! However, this is only needed for
> # esthetic purposes.
> rgCANa <- rgCANa *2^1.4
> plotMvsA(rgCANa)
> title(main="Normalized AC")
> 
> # Curve-fit (lowess) normalization
> rgCANlw <- normalizeLowess(rgCAavg, weights=weights)
Warning message:
In normalizeCurveFit.matrix(X, method = "lowess", ...) :
  Weights were rounded to {0,1} since 'lowess' normalization supports only zero-one weights.
> plotMvsA(rgCANlw, col="orange", add=TRUE)
> 
> # Curve-fit (loess) normalization
> rgCANl <- normalizeLoess(rgCAavg, weights=weights)
> plotMvsA(rgCANl, col="red", add=TRUE)
> 
> # Curve-fit (robust spline) normalization
> rgCANrs <- normalizeRobustSpline(rgCAavg, weights=weights)
> plotMvsA(rgCANrs, col="blue", add=TRUE)
> 
> legend(x=0,y=16, legend=c("affine", "lowess", "loess", "r. spline"), pch=19,
+        col=c("black", "orange", "red", "blue"), ncol=2, x.intersp=0.3, bty="n")
> 
> 
> plotMvsMPairs(cbind(rgCANa, rgCANlw), col="orange", xlab=expression(M[affine]))
> title(main="Normalized AC")
> plotMvsMPairs(cbind(rgCANa, rgCANl), col="red", add=TRUE)
> plotMvsMPairs(cbind(rgCANa, rgCANrs), col="blue", add=TRUE)
> abline(a=0, b=1, lty=2)
> legend(x=-6,y=6, legend=c("lowess", "loess", "r. spline"), pch=19,
+        col=c("orange", "red", "blue"), ncol=2, x.intersp=0.3, bty="n")
> 
> 
> proc.time()
   user  system elapsed 
   7.73    0.06    7.78 

R version 4.0.5 (2021-03-31) -- "Shake and Throw"
Copyright (C) 2021 The R Foundation for Statistical Computing
Platform: i386-w64-mingw32/i386 (32-bit)

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.

> library("aroma.light")
aroma.light v3.20.0 successfully loaded. See ?aroma.light for help.
> 
> # Simulate three shifted tracks of different lengths with same profiles
> ns <- c(A=2, B=1, C=0.25)*1000
> xx <- lapply(ns, FUN=function(n) { seq(from=1, to=max(ns), length.out=n) })
> zz <- mapply(seq_along(ns), ns, FUN=function(z,n) rep(z,n))
> 
> yy <- list(
+   A = rnorm(ns["A"], mean=0, sd=0.5),
+   B = rnorm(ns["B"], mean=5, sd=0.4),
+   C = rnorm(ns["C"], mean=-5, sd=1.1)
+ )
> yy <- lapply(yy, FUN=function(y) {
+   n <- length(y)
+   y[1:(n/2)] <- y[1:(n/2)] + 2
+   y[1:(n/4)] <- y[1:(n/4)] - 4
+   y
+ })
> 
> # Shift all tracks toward the first track
> yyN <- normalizeDifferencesToAverage(yy, baseline=1)
> 
> # The baseline channel is not changed
> stopifnot(identical(yy[[1]], yyN[[1]]))
> 
> # Get the estimated parameters
> fit <- attr(yyN, "fit")
> 
> # Plot the tracks
> layout(matrix(1:2, ncol=1))
> x <- unlist(xx)
> col <- unlist(zz)
> y <- unlist(yy)
> yN <- unlist(yyN)
> plot(x, y, col=col, ylim=c(-10,10))
> plot(x, yN, col=col, ylim=c(-10,10))
> 
> proc.time()
   user  system elapsed 
   0.64    0.06    0.68 

R version 4.0.5 (2021-03-31) -- "Shake and Throw"
Copyright (C) 2021 The R Foundation for Statistical Computing
Platform: x86_64-w64-mingw32/x64 (64-bit)

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.

> library("aroma.light")
aroma.light v3.20.0 successfully loaded. See ?aroma.light for help.
> 
> # Simulate three shifted tracks of different lengths with same profiles
> ns <- c(A=2, B=1, C=0.25)*1000
> xx <- lapply(ns, FUN=function(n) { seq(from=1, to=max(ns), length.out=n) })
> zz <- mapply(seq_along(ns), ns, FUN=function(z,n) rep(z,n))
> 
> yy <- list(
+   A = rnorm(ns["A"], mean=0, sd=0.5),
+   B = rnorm(ns["B"], mean=5, sd=0.4),
+   C = rnorm(ns["C"], mean=-5, sd=1.1)
+ )
> yy <- lapply(yy, FUN=function(y) {
+   n <- length(y)
+   y[1:(n/2)] <- y[1:(n/2)] + 2
+   y[1:(n/4)] <- y[1:(n/4)] - 4
+   y
+ })
> 
> # Shift all tracks toward the first track
> yyN <- normalizeDifferencesToAverage(yy, baseline=1)
> 
> # The baseline channel is not changed
> stopifnot(identical(yy[[1]], yyN[[1]]))
> 
> # Get the estimated parameters
> fit <- attr(yyN, "fit")
> 
> # Plot the tracks
> layout(matrix(1:2, ncol=1))
> x <- unlist(xx)
> col <- unlist(zz)
> y <- unlist(yy)
> yN <- unlist(yyN)
> plot(x, y, col=col, ylim=c(-10,10))
> plot(x, yN, col=col, ylim=c(-10,10))
> 
> proc.time()
   user  system elapsed 
   0.45    0.03    0.48 

R version 4.0.5 (2021-03-31) -- "Shake and Throw"
Copyright (C) 2021 The R Foundation for Statistical Computing
Platform: i386-w64-mingw32/i386 (32-bit)

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.

> library("aroma.light")
aroma.light v3.20.0 successfully loaded. See ?aroma.light for help.
> 
> # - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
> # Example 1: Single-enzyme fragment-length normalization of 6 arrays
> # - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
> # Number samples
> I <- 9
> 
> # Number of loci
> J <- 1000
> 
> # Fragment lengths
> fl <- seq(from=100, to=1000, length.out=J)
> 
> # Simulate data points with unknown fragment lengths
> hasUnknownFL <- seq(from=1, to=J, by=50)
> fl[hasUnknownFL] <- NA_real_
> 
> # Simulate data
> y <- matrix(0, nrow=J, ncol=I)
> maxY <- 12
> for (kk in 1:I) {
+   k <- runif(n=1, min=3, max=5)
+   mu <- function(fl) {
+     mu <- rep(maxY, length(fl))
+     ok <- !is.na(fl)
+     mu[ok] <- mu[ok] - fl[ok]^{1/k}
+     mu
+   }
+   eps <- rnorm(J, mean=0, sd=1)
+   y[,kk] <- mu(fl) + eps
+ }
> 
> # Normalize data (to a zero baseline)
> yN <- apply(y, MARGIN=2, FUN=function(y) {
+   normalizeFragmentLength(y, fragmentLengths=fl, onMissing="median")
+ })
> 
> # The correction factors
> rho <- y-yN
> print(summary(rho))
       V1              V2              V3              V4       
 Min.   :6.278   Min.   :7.509   Min.   :6.276   Min.   :7.874  
 1st Qu.:6.645   1st Qu.:7.824   1st Qu.:6.748   1st Qu.:8.100  
 Median :7.084   Median :8.149   Median :7.185   Median :8.343  
 Mean   :7.185   Mean   :8.206   Mean   :7.296   Mean   :8.400  
 3rd Qu.:7.689   3rd Qu.:8.565   3rd Qu.:7.826   3rd Qu.:8.679  
 Max.   :8.436   Max.   :9.107   Max.   :8.635   Max.   :9.116  
       V5              V6              V7              V8       
 Min.   :7.902   Min.   :6.037   Min.   :7.188   Min.   :3.977  
 1st Qu.:8.111   1st Qu.:6.536   1st Qu.:7.555   1st Qu.:4.607  
 Median :8.416   Median :7.080   Median :7.878   Median :5.315  
 Mean   :8.474   Mean   :7.146   Mean   :7.957   Mean   :5.531  
 3rd Qu.:8.798   3rd Qu.:7.718   3rd Qu.:8.343   3rd Qu.:6.388  
 Max.   :9.306   Max.   :8.550   Max.   :8.965   Max.   :7.814  
       V9       
 Min.   :7.885  
 1st Qu.:8.149  
 Median :8.412  
 Mean   :8.460  
 3rd Qu.:8.763  
 Max.   :9.182  
> # The correction for units with unknown fragment lengths
> # equals the median correction factor of all other units
> print(summary(rho[hasUnknownFL,]))
       V1              V2              V3              V4       
 Min.   :7.084   Min.   :8.149   Min.   :7.185   Min.   :8.343  
 1st Qu.:7.084   1st Qu.:8.149   1st Qu.:7.185   1st Qu.:8.343  
 Median :7.084   Median :8.149   Median :7.185   Median :8.343  
 Mean   :7.084   Mean   :8.149   Mean   :7.185   Mean   :8.343  
 3rd Qu.:7.084   3rd Qu.:8.149   3rd Qu.:7.185   3rd Qu.:8.343  
 Max.   :7.084   Max.   :8.149   Max.   :7.185   Max.   :8.343  
       V5              V6             V7              V8              V9       
 Min.   :8.416   Min.   :7.08   Min.   :7.878   Min.   :5.315   Min.   :8.412  
 1st Qu.:8.416   1st Qu.:7.08   1st Qu.:7.878   1st Qu.:5.315   1st Qu.:8.412  
 Median :8.416   Median :7.08   Median :7.878   Median :5.315   Median :8.412  
 Mean   :8.416   Mean   :7.08   Mean   :7.878   Mean   :5.315   Mean   :8.412  
 3rd Qu.:8.416   3rd Qu.:7.08   3rd Qu.:7.878   3rd Qu.:5.315   3rd Qu.:8.412  
 Max.   :8.416   Max.   :7.08   Max.   :7.878   Max.   :5.315   Max.   :8.412  
> 
> # Plot raw data
> layout(matrix(1:9, ncol=3))
> xlim <- c(0,max(fl, na.rm=TRUE))
> ylim <- c(0,max(y, na.rm=TRUE))
> xlab <- "Fragment length"
> ylab <- expression(log2(theta))
> for (kk in 1:I) {
+   plot(fl, y[,kk], xlim=xlim, ylim=ylim, xlab=xlab, ylab=ylab)
+   ok <- (is.finite(fl) & is.finite(y[,kk]))
+   lines(lowess(fl[ok], y[ok,kk]), col="red", lwd=2)
+ }
> 
> # Plot normalized data
> layout(matrix(1:9, ncol=3))
> ylim <- c(-1,1)*max(y, na.rm=TRUE)/2
> for (kk in 1:I) {
+   plot(fl, yN[,kk], xlim=xlim, ylim=ylim, xlab=xlab, ylab=ylab)
+   ok <- (is.finite(fl) & is.finite(y[,kk]))
+   lines(lowess(fl[ok], yN[ok,kk]), col="blue", lwd=2)
+ }
> 
> proc.time()
   user  system elapsed 
   1.64    0.09    1.71 

R version 4.0.5 (2021-03-31) -- "Shake and Throw"
Copyright (C) 2021 The R Foundation for Statistical Computing
Platform: x86_64-w64-mingw32/x64 (64-bit)

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.

> library("aroma.light")
aroma.light v3.20.0 successfully loaded. See ?aroma.light for help.
> 
> # - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
> # Example 1: Single-enzyme fragment-length normalization of 6 arrays
> # - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
> # Number samples
> I <- 9
> 
> # Number of loci
> J <- 1000
> 
> # Fragment lengths
> fl <- seq(from=100, to=1000, length.out=J)
> 
> # Simulate data points with unknown fragment lengths
> hasUnknownFL <- seq(from=1, to=J, by=50)
> fl[hasUnknownFL] <- NA_real_
> 
> # Simulate data
> y <- matrix(0, nrow=J, ncol=I)
> maxY <- 12
> for (kk in 1:I) {
+   k <- runif(n=1, min=3, max=5)
+   mu <- function(fl) {
+     mu <- rep(maxY, length(fl))
+     ok <- !is.na(fl)
+     mu[ok] <- mu[ok] - fl[ok]^{1/k}
+     mu
+   }
+   eps <- rnorm(J, mean=0, sd=1)
+   y[,kk] <- mu(fl) + eps
+ }
> 
> # Normalize data (to a zero baseline)
> yN <- apply(y, MARGIN=2, FUN=function(y) {
+   normalizeFragmentLength(y, fragmentLengths=fl, onMissing="median")
+ })
> 
> # The correction factors
> rho <- y-yN
> print(summary(rho))
       V1              V2              V3              V4       
 Min.   :7.900   Min.   :5.002   Min.   :3.815   Min.   :7.406  
 1st Qu.:8.171   1st Qu.:5.469   1st Qu.:4.515   1st Qu.:7.583  
 Median :8.434   Median :6.033   Median :5.223   Median :7.861  
 Mean   :8.514   Mean   :6.204   Mean   :5.426   Mean   :8.007  
 3rd Qu.:8.823   3rd Qu.:6.903   3rd Qu.:6.289   3rd Qu.:8.396  
 Max.   :9.418   Max.   :7.924   Max.   :7.665   Max.   :9.070  
       V5              V6              V7              V8       
 Min.   :5.661   Min.   :5.563   Min.   :7.770   Min.   :7.884  
 1st Qu.:6.164   1st Qu.:5.956   1st Qu.:8.009   1st Qu.:8.203  
 Median :6.616   Median :6.451   Median :8.305   Median :8.463  
 Mean   :6.747   Mean   :6.606   Mean   :8.370   Mean   :8.528  
 3rd Qu.:7.292   3rd Qu.:7.202   3rd Qu.:8.692   3rd Qu.:8.849  
 Max.   :8.288   Max.   :8.189   Max.   :9.258   Max.   :9.326  
       V9       
 Min.   :6.883  
 1st Qu.:7.187  
 Median :7.493  
 Mean   :7.618  
 3rd Qu.:8.026  
 Max.   :8.725  
> # The correction for units with unknown fragment lengths
> # equals the median correction factor of all other units
> print(summary(rho[hasUnknownFL,]))
       V1              V2              V3              V4       
 Min.   :8.434   Min.   :6.033   Min.   :5.223   Min.   :7.861  
 1st Qu.:8.434   1st Qu.:6.033   1st Qu.:5.223   1st Qu.:7.861  
 Median :8.434   Median :6.033   Median :5.223   Median :7.861  
 Mean   :8.434   Mean   :6.033   Mean   :5.223   Mean   :7.861  
 3rd Qu.:8.434   3rd Qu.:6.033   3rd Qu.:5.223   3rd Qu.:7.861  
 Max.   :8.434   Max.   :6.033   Max.   :5.223   Max.   :7.861  
       V5              V6              V7              V8       
 Min.   :6.616   Min.   :6.451   Min.   :8.305   Min.   :8.463  
 1st Qu.:6.616   1st Qu.:6.451   1st Qu.:8.305   1st Qu.:8.463  
 Median :6.616   Median :6.451   Median :8.305   Median :8.463  
 Mean   :6.616   Mean   :6.451   Mean   :8.305   Mean   :8.463  
 3rd Qu.:6.616   3rd Qu.:6.451   3rd Qu.:8.305   3rd Qu.:8.463  
 Max.   :6.616   Max.   :6.451   Max.   :8.305   Max.   :8.463  
       V9       
 Min.   :7.493  
 1st Qu.:7.493  
 Median :7.493  
 Mean   :7.493  
 3rd Qu.:7.493  
 Max.   :7.493  
> 
> # Plot raw data
> layout(matrix(1:9, ncol=3))
> xlim <- c(0,max(fl, na.rm=TRUE))
> ylim <- c(0,max(y, na.rm=TRUE))
> xlab <- "Fragment length"
> ylab <- expression(log2(theta))
> for (kk in 1:I) {
+   plot(fl, y[,kk], xlim=xlim, ylim=ylim, xlab=xlab, ylab=ylab)
+   ok <- (is.finite(fl) & is.finite(y[,kk]))
+   lines(lowess(fl[ok], y[ok,kk]), col="red", lwd=2)
+ }
> 
> # Plot normalized data
> layout(matrix(1:9, ncol=3))
> ylim <- c(-1,1)*max(y, na.rm=TRUE)/2
> for (kk in 1:I) {
+   plot(fl, yN[,kk], xlim=xlim, ylim=ylim, xlab=xlab, ylab=ylab)
+   ok <- (is.finite(fl) & is.finite(y[,kk]))
+   lines(lowess(fl[ok], yN[ok,kk]), col="blue", lwd=2)
+ }
> 
> proc.time()
   user  system elapsed 
   1.46    0.04    1.50 

R version 4.0.5 (2021-03-31) -- "Shake and Throw"
Copyright (C) 2021 The R Foundation for Statistical Computing
Platform: i386-w64-mingw32/i386 (32-bit)

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.

> library("aroma.light")
aroma.light v3.20.0 successfully loaded. See ?aroma.light for help.
> 
> # - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
> # Example 2: Two-enzyme fragment-length normalization of 6 arrays
> # - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
> set.seed(0xbeef)
> 
> # Number samples
> I <- 5
> 
> # Number of loci
> J <- 3000
> 
> # Fragment lengths (two enzymes)
> fl <- matrix(0, nrow=J, ncol=2)
> fl[,1] <- seq(from=100, to=1000, length.out=J)
> fl[,2] <- seq(from=1000, to=100, length.out=J)
> 
> # Let 1/2 of the units be on both enzymes
> fl[seq(from=1, to=J, by=4),1] <- NA_real_
> fl[seq(from=2, to=J, by=4),2] <- NA_real_
> 
> # Let some have unknown fragment lengths
> hasUnknownFL <- seq(from=1, to=J, by=15)
> fl[hasUnknownFL,] <- NA_real_
> 
> # Sty/Nsp mixing proportions:
> rho <- rep(1, I)
> rho[1] <- 1/3;  # Less Sty in 1st sample
> rho[3] <- 3/2;  # More Sty in 3rd sample
> 
> 
> # Simulate data
> z <- array(0, dim=c(J,2,I))
> maxLog2Theta <- 12
> for (ii in 1:I) {
+   # Common effect for both enzymes
+   mu <- function(fl) {
+     k <- runif(n=1, min=3, max=5)
+     mu <- rep(maxLog2Theta, length(fl))
+     ok <- is.finite(fl)
+     mu[ok] <- mu[ok] - fl[ok]^{1/k}
+     mu
+   }
+ 
+   # Calculate the effect for each data point
+   for (ee in 1:2) {
+     z[,ee,ii] <- mu(fl[,ee])
+   }
+ 
+   # Update the Sty/Nsp mixing proportions
+   ee <- 2
+   z[,ee,ii] <- rho[ii]*z[,ee,ii]
+ 
+   # Add random errors
+   for (ee in 1:2) {
+     eps <- rnorm(J, mean=0, sd=1/sqrt(2))
+     z[,ee,ii] <- z[,ee,ii] + eps
+   }
+ }
> 
> 
> hasFl <- is.finite(fl)
> 
> unitSets <- list(
+   nsp  = which( hasFl[,1] & !hasFl[,2]),
+   sty  = which(!hasFl[,1] &  hasFl[,2]),
+   both = which( hasFl[,1] &  hasFl[,2]),
+   none = which(!hasFl[,1] & !hasFl[,2])
+ )
> 
> # The observed data is a mix of two enzymes
> theta <- matrix(NA_real_, nrow=J, ncol=I)
> 
> # Single-enzyme units
> for (ee in 1:2) {
+   uu <- unitSets[[ee]]
+   theta[uu,] <- 2^z[uu,ee,]
+ }
> 
> # Both-enzyme units (sum on intensity scale)
> uu <- unitSets$both
> theta[uu,] <- (2^z[uu,1,]+2^z[uu,2,])/2
> 
> # Missing units (sample from the others)
> uu <- unitSets$none
> theta[uu,] <- apply(theta, MARGIN=2, sample, size=length(uu))
> 
> # Calculate target array
> thetaT <- rowMeans(theta, na.rm=TRUE)
> targetFcns <- list()
> for (ee in 1:2) {
+   uu <- unitSets[[ee]]
+   fit <- lowess(fl[uu,ee], log2(thetaT[uu]))
+   class(fit) <- "lowess"
+   targetFcns[[ee]] <- function(fl, ...) {
+     predict(fit, newdata=fl)
+   }
+ }
> 
> 
> # Fit model only to a subset of the data
> subsetToFit <- setdiff(1:J, seq(from=1, to=J, by=10))
> 
> # Normalize data (to a target baseline)
> thetaN <- matrix(NA_real_, nrow=J, ncol=I)
> fits <- vector("list", I)
> for (ii in 1:I) {
+   lthetaNi <- normalizeFragmentLength(log2(theta[,ii]), targetFcns=targetFcns,
+                      fragmentLengths=fl, onMissing="median",
+                      subsetToFit=subsetToFit, .returnFit=TRUE)
+   fits[[ii]] <- attr(lthetaNi, "modelFit")
+   thetaN[,ii] <- 2^lthetaNi
+ }
> 
> 
> # Plot raw data
> xlim <- c(0, max(fl, na.rm=TRUE))
> ylim <- c(0, max(log2(theta), na.rm=TRUE))
> Mlim <- c(-1,1)*4
> xlab <- "Fragment length"
> ylab <- expression(log2(theta))
> Mlab <- expression(M==log[2](theta/theta[R]))
> 
> layout(matrix(1:(3*I), ncol=I, byrow=TRUE))
> for (ii in 1:I) {
+   plot(NA, xlim=xlim, ylim=ylim, xlab=xlab, ylab=ylab, main="raw")
+ 
+   # Single-enzyme units
+   for (ee in 1:2) {
+     # The raw data
+     uu <- unitSets[[ee]]
+     points(fl[uu,ee], log2(theta[uu,ii]), col=ee+1)
+   }
+ 
+   # Both-enzyme units (use fragment-length for enzyme #1)
+   uu <- unitSets$both
+   points(fl[uu,1], log2(theta[uu,ii]), col=3+1)
+ 
+   for (ee in 1:2) {
+     # The true effects
+     uu <- unitSets[[ee]]
+     lines(lowess(fl[uu,ee], log2(theta[uu,ii])), col="black", lwd=4, lty=3)
+ 
+     # The estimated effects
+     fit <- fits[[ii]][[ee]]$fit
+     lines(fit, col="orange", lwd=3)
+ 
+     muT <- targetFcns[[ee]](fl[uu,ee])
+     lines(fl[uu,ee], muT, col="cyan", lwd=1)
+   }
+ }
> 
> # Calculate log-ratios
> thetaR <- rowMeans(thetaN, na.rm=TRUE)
> M <- log2(thetaN/thetaR)
> 
> # Plot normalized data
> for (ii in 1:I) {
+   plot(NA, xlim=xlim, ylim=Mlim, xlab=xlab, ylab=Mlab, main="normalized")
+   # Single-enzyme units
+   for (ee in 1:2) {
+     # The normalized data
+     uu <- unitSets[[ee]]
+     points(fl[uu,ee], M[uu,ii], col=ee+1)
+   }
+   # Both-enzyme units (use fragment-length for enzyme #1)
+   uu <- unitSets$both
+   points(fl[uu,1], M[uu,ii], col=3+1)
+ }
> 
> ylim <- c(0,1.5)
> for (ii in 1:I) {
+   data <- list()
+   for (ee in 1:2) {
+     # The normalized data
+     uu <- unitSets[[ee]]
+     data[[ee]] <- M[uu,ii]
+   }
+   uu <- unitSets$both
+   if (length(uu) > 0)
+     data[[3]] <- M[uu,ii]
+ 
+   uu <- unitSets$none
+   if (length(uu) > 0)
+     data[[4]] <- M[uu,ii]
+ 
+   cols <- seq_along(data)+1
+   plotDensity(data, col=cols, xlim=Mlim, xlab=Mlab, main="normalized")
+ 
+   abline(v=0, lty=2)
+ }
> 
> 
> proc.time()
   user  system elapsed 
   1.59    0.10    1.68 

R version 4.0.5 (2021-03-31) -- "Shake and Throw"
Copyright (C) 2021 The R Foundation for Statistical Computing
Platform: x86_64-w64-mingw32/x64 (64-bit)

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.

> library("aroma.light")
aroma.light v3.20.0 successfully loaded. See ?aroma.light for help.
> 
> # - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
> # Example 2: Two-enzyme fragment-length normalization of 6 arrays
> # - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
> set.seed(0xbeef)
> 
> # Number samples
> I <- 5
> 
> # Number of loci
> J <- 3000
> 
> # Fragment lengths (two enzymes)
> fl <- matrix(0, nrow=J, ncol=2)
> fl[,1] <- seq(from=100, to=1000, length.out=J)
> fl[,2] <- seq(from=1000, to=100, length.out=J)
> 
> # Let 1/2 of the units be on both enzymes
> fl[seq(from=1, to=J, by=4),1] <- NA_real_
> fl[seq(from=2, to=J, by=4),2] <- NA_real_
> 
> # Let some have unknown fragment lengths
> hasUnknownFL <- seq(from=1, to=J, by=15)
> fl[hasUnknownFL,] <- NA_real_
> 
> # Sty/Nsp mixing proportions:
> rho <- rep(1, I)
> rho[1] <- 1/3;  # Less Sty in 1st sample
> rho[3] <- 3/2;  # More Sty in 3rd sample
> 
> 
> # Simulate data
> z <- array(0, dim=c(J,2,I))
> maxLog2Theta <- 12
> for (ii in 1:I) {
+   # Common effect for both enzymes
+   mu <- function(fl) {
+     k <- runif(n=1, min=3, max=5)
+     mu <- rep(maxLog2Theta, length(fl))
+     ok <- is.finite(fl)
+     mu[ok] <- mu[ok] - fl[ok]^{1/k}
+     mu
+   }
+ 
+   # Calculate the effect for each data point
+   for (ee in 1:2) {
+     z[,ee,ii] <- mu(fl[,ee])
+   }
+ 
+   # Update the Sty/Nsp mixing proportions
+   ee <- 2
+   z[,ee,ii] <- rho[ii]*z[,ee,ii]
+ 
+   # Add random errors
+   for (ee in 1:2) {
+     eps <- rnorm(J, mean=0, sd=1/sqrt(2))
+     z[,ee,ii] <- z[,ee,ii] + eps
+   }
+ }
> 
> 
> hasFl <- is.finite(fl)
> 
> unitSets <- list(
+   nsp  = which( hasFl[,1] & !hasFl[,2]),
+   sty  = which(!hasFl[,1] &  hasFl[,2]),
+   both = which( hasFl[,1] &  hasFl[,2]),
+   none = which(!hasFl[,1] & !hasFl[,2])
+ )
> 
> # The observed data is a mix of two enzymes
> theta <- matrix(NA_real_, nrow=J, ncol=I)
> 
> # Single-enzyme units
> for (ee in 1:2) {
+   uu <- unitSets[[ee]]
+   theta[uu,] <- 2^z[uu,ee,]
+ }
> 
> # Both-enzyme units (sum on intensity scale)
> uu <- unitSets$both
> theta[uu,] <- (2^z[uu,1,]+2^z[uu,2,])/2
> 
> # Missing units (sample from the others)
> uu <- unitSets$none
> theta[uu,] <- apply(theta, MARGIN=2, sample, size=length(uu))
> 
> # Calculate target array
> thetaT <- rowMeans(theta, na.rm=TRUE)
> targetFcns <- list()
> for (ee in 1:2) {
+   uu <- unitSets[[ee]]
+   fit <- lowess(fl[uu,ee], log2(thetaT[uu]))
+   class(fit) <- "lowess"
+   targetFcns[[ee]] <- function(fl, ...) {
+     predict(fit, newdata=fl)
+   }
+ }
> 
> 
> # Fit model only to a subset of the data
> subsetToFit <- setdiff(1:J, seq(from=1, to=J, by=10))
> 
> # Normalize data (to a target baseline)
> thetaN <- matrix(NA_real_, nrow=J, ncol=I)
> fits <- vector("list", I)
> for (ii in 1:I) {
+   lthetaNi <- normalizeFragmentLength(log2(theta[,ii]), targetFcns=targetFcns,
+                      fragmentLengths=fl, onMissing="median",
+                      subsetToFit=subsetToFit, .returnFit=TRUE)
+   fits[[ii]] <- attr(lthetaNi, "modelFit")
+   thetaN[,ii] <- 2^lthetaNi
+ }
> 
> 
> # Plot raw data
> xlim <- c(0, max(fl, na.rm=TRUE))
> ylim <- c(0, max(log2(theta), na.rm=TRUE))
> Mlim <- c(-1,1)*4
> xlab <- "Fragment length"
> ylab <- expression(log2(theta))
> Mlab <- expression(M==log[2](theta/theta[R]))
> 
> layout(matrix(1:(3*I), ncol=I, byrow=TRUE))
> for (ii in 1:I) {
+   plot(NA, xlim=xlim, ylim=ylim, xlab=xlab, ylab=ylab, main="raw")
+ 
+   # Single-enzyme units
+   for (ee in 1:2) {
+     # The raw data
+     uu <- unitSets[[ee]]
+     points(fl[uu,ee], log2(theta[uu,ii]), col=ee+1)
+   }
+ 
+   # Both-enzyme units (use fragment-length for enzyme #1)
+   uu <- unitSets$both
+   points(fl[uu,1], log2(theta[uu,ii]), col=3+1)
+ 
+   for (ee in 1:2) {
+     # The true effects
+     uu <- unitSets[[ee]]
+     lines(lowess(fl[uu,ee], log2(theta[uu,ii])), col="black", lwd=4, lty=3)
+ 
+     # The estimated effects
+     fit <- fits[[ii]][[ee]]$fit
+     lines(fit, col="orange", lwd=3)
+ 
+     muT <- targetFcns[[ee]](fl[uu,ee])
+     lines(fl[uu,ee], muT, col="cyan", lwd=1)
+   }
+ }
> 
> # Calculate log-ratios
> thetaR <- rowMeans(thetaN, na.rm=TRUE)
> M <- log2(thetaN/thetaR)
> 
> # Plot normalized data
> for (ii in 1:I) {
+   plot(NA, xlim=xlim, ylim=Mlim, xlab=xlab, ylab=Mlab, main="normalized")
+   # Single-enzyme units
+   for (ee in 1:2) {
+     # The normalized data
+     uu <- unitSets[[ee]]
+     points(fl[uu,ee], M[uu,ii], col=ee+1)
+   }
+   # Both-enzyme units (use fragment-length for enzyme #1)
+   uu <- unitSets$both
+   points(fl[uu,1], M[uu,ii], col=3+1)
+ }
> 
> ylim <- c(0,1.5)
> for (ii in 1:I) {
+   data <- list()
+   for (ee in 1:2) {
+     # The normalized data
+     uu <- unitSets[[ee]]
+     data[[ee]] <- M[uu,ii]
+   }
+   uu <- unitSets$both
+   if (length(uu) > 0)
+     data[[3]] <- M[uu,ii]
+ 
+   uu <- unitSets$none
+   if (length(uu) > 0)
+     data[[4]] <- M[uu,ii]
+ 
+   cols <- seq_along(data)+1
+   plotDensity(data, col=cols, xlim=Mlim, xlab=Mlab, main="normalized")
+ 
+   abline(v=0, lty=2)
+ }
> 
> 
> proc.time()
   user  system elapsed 
   1.04    0.10    1.12 

R version 4.0.5 (2021-03-31) -- "Shake and Throw"
Copyright (C) 2021 The R Foundation for Statistical Computing
Platform: i386-w64-mingw32/i386 (32-bit)

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.

> library("aroma.light")
aroma.light v3.20.0 successfully loaded. See ?aroma.light for help.
> 
> # Simulate ten samples of different lengths
> N <- 10000
> X <- list()
> for (kk in 1:8) {
+   rfcn <- list(rnorm, rgamma)[[sample(2, size=1)]]
+   size <- runif(1, min=0.3, max=1)
+   a <- rgamma(1, shape=20, rate=10)
+   b <- rgamma(1, shape=10, rate=10)
+   values <- rfcn(size*N, a, b)
+ 
+   # "Censor" values
+   values[values < 0 | values > 8] <- NA_real_
+ 
+   X[[kk]] <- values
+ }
> 
> # Add 20% missing values
> X <- lapply(X, FUN=function(x) {
+   x[sample(length(x), size=0.20*length(x))] <- NA_real_
+   x
+ })
> 
> # Normalize quantiles
> Xn <- normalizeQuantile(X)
> 
> # Plot the data
> layout(matrix(1:2, ncol=1))
> xlim <- range(X, na.rm=TRUE)
> plotDensity(X, lwd=2, xlim=xlim, main="The original distributions")
> plotDensity(Xn, lwd=2, xlim=xlim, main="The normalized distributions")
> 
> proc.time()
   user  system elapsed 
   0.71    0.04    0.75 

R version 4.0.5 (2021-03-31) -- "Shake and Throw"
Copyright (C) 2021 The R Foundation for Statistical Computing
Platform: x86_64-w64-mingw32/x64 (64-bit)

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.

> library("aroma.light")
aroma.light v3.20.0 successfully loaded. See ?aroma.light for help.
> 
> # Simulate ten samples of different lengths
> N <- 10000
> X <- list()
> for (kk in 1:8) {
+   rfcn <- list(rnorm, rgamma)[[sample(2, size=1)]]
+   size <- runif(1, min=0.3, max=1)
+   a <- rgamma(1, shape=20, rate=10)
+   b <- rgamma(1, shape=10, rate=10)
+   values <- rfcn(size*N, a, b)
+ 
+   # "Censor" values
+   values[values < 0 | values > 8] <- NA_real_
+ 
+   X[[kk]] <- values
+ }
> 
> # Add 20% missing values
> X <- lapply(X, FUN=function(x) {
+   x[sample(length(x), size=0.20*length(x))] <- NA_real_
+   x
+ })
> 
> # Normalize quantiles
> Xn <- normalizeQuantile(X)
> 
> # Plot the data
> layout(matrix(1:2, ncol=1))
> xlim <- range(X, na.rm=TRUE)
> plotDensity(X, lwd=2, xlim=xlim, main="The original distributions")
> plotDensity(Xn, lwd=2, xlim=xlim, main="The normalized distributions")
> 
> proc.time()
   user  system elapsed 
   0.51    0.07    0.57 

R version 4.0.5 (2021-03-31) -- "Shake and Throw"
Copyright (C) 2021 The R Foundation for Statistical Computing
Platform: i386-w64-mingw32/i386 (32-bit)

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.

> library("aroma.light")
aroma.light v3.20.0 successfully loaded. See ?aroma.light for help.
> 
> # Simulate three samples with on average 20% missing values
> N <- 10000
> X <- cbind(rnorm(N, mean=3, sd=1),
+            rnorm(N, mean=4, sd=2),
+            rgamma(N, shape=2, rate=1))
> X[sample(3*N, size=0.20*3*N)] <- NA_real_
> 
> # Normalize quantiles
> Xn <- normalizeQuantile(X)
> 
> # Plot the data
> layout(matrix(1:2, ncol=1))
> xlim <- range(X, Xn, na.rm=TRUE)
> plotDensity(X, lwd=2, xlim=xlim, main="The three original distributions")
> plotDensity(Xn, lwd=2, xlim=xlim, main="The three normalized distributions")
> 
> proc.time()
   user  system elapsed 
   0.43    0.06    0.46 

R version 4.0.5 (2021-03-31) -- "Shake and Throw"
Copyright (C) 2021 The R Foundation for Statistical Computing
Platform: x86_64-w64-mingw32/x64 (64-bit)

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.

> library("aroma.light")
aroma.light v3.20.0 successfully loaded. See ?aroma.light for help.
> 
> # Simulate three samples with on average 20% missing values
> N <- 10000
> X <- cbind(rnorm(N, mean=3, sd=1),
+            rnorm(N, mean=4, sd=2),
+            rgamma(N, shape=2, rate=1))
> X[sample(3*N, size=0.20*3*N)] <- NA_real_
> 
> # Normalize quantiles
> Xn <- normalizeQuantile(X)
> 
> # Plot the data
> layout(matrix(1:2, ncol=1))
> xlim <- range(X, Xn, na.rm=TRUE)
> plotDensity(X, lwd=2, xlim=xlim, main="The three original distributions")
> plotDensity(Xn, lwd=2, xlim=xlim, main="The three normalized distributions")
> 
> proc.time()
   user  system elapsed 
   0.57    0.06    0.62 

R version 4.0.5 (2021-03-31) -- "Shake and Throw"
Copyright (C) 2021 The R Foundation for Statistical Computing
Platform: i386-w64-mingw32/i386 (32-bit)

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.

> library("aroma.light")
aroma.light v3.20.0 successfully loaded. See ?aroma.light for help.
> 
> # Simulate three samples with on average 20% missing values
> N <- 10000
> X <- cbind(rnorm(N, mean=3, sd=1),
+            rnorm(N, mean=4, sd=2),
+            rgamma(N, shape=2, rate=1))
> X[sample(3*N, size=0.20*3*N)] <- NA_real_
> 
> # Plot the data
> layout(matrix(c(1,0,2:5), ncol=2, byrow=TRUE))
> xlim <- range(X, na.rm=TRUE)
> plotDensity(X, lwd=2, xlim=xlim, main="The three original distributions")
> 
> Xn <- normalizeQuantile(X)
> plotDensity(Xn, lwd=2, xlim=xlim, main="The three normalized distributions")
> plotXYCurve(X, Xn, xlim=xlim, main="The three normalized distributions")
> 
> Xn2 <- normalizeQuantileSpline(X, xTarget=Xn[,1], spar=0.99)
> plotDensity(Xn2, lwd=2, xlim=xlim, main="The three normalized distributions")
> plotXYCurve(X, Xn2, xlim=xlim, main="The three normalized distributions")
> 
> proc.time()
   user  system elapsed 
   1.09    0.10    1.20 

R version 4.0.5 (2021-03-31) -- "Shake and Throw"
Copyright (C) 2021 The R Foundation for Statistical Computing
Platform: x86_64-w64-mingw32/x64 (64-bit)

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.

> library("aroma.light")
aroma.light v3.20.0 successfully loaded. See ?aroma.light for help.
> 
> # Simulate three samples with on average 20% missing values
> N <- 10000
> X <- cbind(rnorm(N, mean=3, sd=1),
+            rnorm(N, mean=4, sd=2),
+            rgamma(N, shape=2, rate=1))
> X[sample(3*N, size=0.20*3*N)] <- NA_real_
> 
> # Plot the data
> layout(matrix(c(1,0,2:5), ncol=2, byrow=TRUE))
> xlim <- range(X, na.rm=TRUE)
> plotDensity(X, lwd=2, xlim=xlim, main="The three original distributions")
> 
> Xn <- normalizeQuantile(X)
> plotDensity(Xn, lwd=2, xlim=xlim, main="The three normalized distributions")
> plotXYCurve(X, Xn, xlim=xlim, main="The three normalized distributions")
> 
> Xn2 <- normalizeQuantileSpline(X, xTarget=Xn[,1], spar=0.99)
> plotDensity(Xn2, lwd=2, xlim=xlim, main="The three normalized distributions")
> plotXYCurve(X, Xn2, xlim=xlim, main="The three normalized distributions")
> 
> proc.time()
   user  system elapsed 
   1.46    0.07    1.53 

R version 4.0.5 (2021-03-31) -- "Shake and Throw"
Copyright (C) 2021 The R Foundation for Statistical Computing
Platform: i386-w64-mingw32/i386 (32-bit)

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.

> library("aroma.light")
aroma.light v3.20.0 successfully loaded. See ?aroma.light for help.
> library("R.utils")
Loading required package: R.oo
Loading required package: R.methodsS3
R.methodsS3 v1.8.1 (2020-08-26 16:20:06 UTC) successfully loaded. See ?R.methodsS3 for help.
R.oo v1.24.0 (2020-08-26 16:11:58 UTC) successfully loaded. See ?R.oo for help.

Attaching package: 'R.oo'

The following object is masked from 'package:R.methodsS3':

    throw

The following objects are masked from 'package:methods':

    getClasses, getMethods

The following objects are masked from 'package:base':

    attach, detach, load, save

R.utils v2.10.1 (2020-08-26 22:50:31 UTC) successfully loaded. See ?R.utils for help.

Attaching package: 'R.utils'

The following object is masked from 'package:utils':

    timestamp

The following objects are masked from 'package:base':

    cat, commandArgs, getOption, inherits, isOpen, nullfile, parse,
    warnings

> 
> # Load data
> pathname <- system.file("data-ex/TumorBoost,fracB,exampleData.Rbin", package="aroma.light")
> data <- loadObject(pathname)
> attachLocally(data)
> pos <- position/1e6
> muN <- genotypeN
> 
> layout(matrix(1:4, ncol=1))
> par(mar=c(2.5,4,0.5,1)+0.1)
> ylim <- c(-0.05, 1.05)
> col <- rep("#999999", length(muN))
> col[muN == 1/2] <- "#000000"
> 
> # Allele B fractions for the normal sample
> plot(pos, betaN, col=col, ylim=ylim)
> 
> # Allele B fractions for the tumor sample
> plot(pos, betaT, col=col, ylim=ylim)
> 
> # TumorBoost w/ naive genotype calls
> betaTN <- normalizeTumorBoost(betaT=betaT, betaN=betaN, preserveScale=FALSE)
> plot(pos, betaTN, col=col, ylim=ylim)
> 
> # TumorBoost w/ external multi-sample genotype calls
> betaTNx <- normalizeTumorBoost(betaT=betaT, betaN=betaN, muN=muN, preserveScale=FALSE)
> plot(pos, betaTNx, col=col, ylim=ylim)
> 
> proc.time()
   user  system elapsed 
   1.12    0.06    1.18 

R version 4.0.5 (2021-03-31) -- "Shake and Throw"
Copyright (C) 2021 The R Foundation for Statistical Computing
Platform: x86_64-w64-mingw32/x64 (64-bit)

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.

> library("aroma.light")
aroma.light v3.20.0 successfully loaded. See ?aroma.light for help.
> library("R.utils")
Loading required package: R.oo
Loading required package: R.methodsS3
R.methodsS3 v1.8.1 (2020-08-26 16:20:06 UTC) successfully loaded. See ?R.methodsS3 for help.
R.oo v1.24.0 (2020-08-26 16:11:58 UTC) successfully loaded. See ?R.oo for help.

Attaching package: 'R.oo'

The following object is masked from 'package:R.methodsS3':

    throw

The following objects are masked from 'package:methods':

    getClasses, getMethods

The following objects are masked from 'package:base':

    attach, detach, load, save

R.utils v2.10.1 (2020-08-26 22:50:31 UTC) successfully loaded. See ?R.utils for help.

Attaching package: 'R.utils'

The following object is masked from 'package:utils':

    timestamp

The following objects are masked from 'package:base':

    cat, commandArgs, getOption, inherits, isOpen, nullfile, parse,
    warnings

> 
> # Load data
> pathname <- system.file("data-ex/TumorBoost,fracB,exampleData.Rbin", package="aroma.light")
> data <- loadObject(pathname)
> attachLocally(data)
> pos <- position/1e6
> muN <- genotypeN
> 
> layout(matrix(1:4, ncol=1))
> par(mar=c(2.5,4,0.5,1)+0.1)
> ylim <- c(-0.05, 1.05)
> col <- rep("#999999", length(muN))
> col[muN == 1/2] <- "#000000"
> 
> # Allele B fractions for the normal sample
> plot(pos, betaN, col=col, ylim=ylim)
> 
> # Allele B fractions for the tumor sample
> plot(pos, betaT, col=col, ylim=ylim)
> 
> # TumorBoost w/ naive genotype calls
> betaTN <- normalizeTumorBoost(betaT=betaT, betaN=betaN, preserveScale=FALSE)
> plot(pos, betaTN, col=col, ylim=ylim)
> 
> # TumorBoost w/ external multi-sample genotype calls
> betaTNx <- normalizeTumorBoost(betaT=betaT, betaN=betaN, muN=muN, preserveScale=FALSE)
> plot(pos, betaTNx, col=col, ylim=ylim)
> 
> proc.time()
   user  system elapsed 
   1.00    0.06    1.04 

R version 4.0.5 (2021-03-31) -- "Shake and Throw"
Copyright (C) 2021 The R Foundation for Statistical Computing
Platform: i386-w64-mingw32/i386 (32-bit)

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.

> library("aroma.light")
aroma.light v3.20.0 successfully loaded. See ?aroma.light for help.
> library("R.utils")
Loading required package: R.oo
Loading required package: R.methodsS3
R.methodsS3 v1.8.1 (2020-08-26 16:20:06 UTC) successfully loaded. See ?R.methodsS3 for help.
R.oo v1.24.0 (2020-08-26 16:11:58 UTC) successfully loaded. See ?R.oo for help.

Attaching package: 'R.oo'

The following object is masked from 'package:R.methodsS3':

    throw

The following objects are masked from 'package:methods':

    getClasses, getMethods

The following objects are masked from 'package:base':

    attach, detach, load, save

R.utils v2.10.1 (2020-08-26 22:50:31 UTC) successfully loaded. See ?R.utils for help.

Attaching package: 'R.utils'

The following object is masked from 'package:utils':

    timestamp

The following objects are masked from 'package:base':

    cat, commandArgs, getOption, inherits, isOpen, nullfile, parse,
    warnings

> 
> # Load data
> pathname <- system.file("data-ex/TumorBoost,fracB,exampleData.Rbin", package="aroma.light")
> data <- loadObject(pathname)
> 
> # Drop loci with missing values
> data <- na.omit(data)
> 
> attachLocally(data)
> pos <- position/1e6
> 
> # Call naive genotypes
> muN <- callNaiveGenotypes(betaN)
> 
> # Genotype classes
> isAA <- (muN == 0)
> isAB <- (muN == 1/2)
> isBB <- (muN == 1)
> 
> # Sanity checks
> stopifnot(all(muN[isAA] == 0))
> stopifnot(all(muN[isAB] == 1/2))
> stopifnot(all(muN[isBB] == 1))
> 
> # TumorBoost normalization with different flavors
> betaTNs <- list()
> for (flavor in c("v1", "v2", "v3", "v4")) {
+   betaTN <- normalizeTumorBoost(betaT=betaT, betaN=betaN, preserveScale=FALSE, flavor=flavor)
+ 
+   # Assert that no non-finite values are introduced
+   stopifnot(all(is.finite(betaTN)))
+ 
+   # Assert that nothing is flipped
+   stopifnot(all(betaTN[isAA] < 1/2))
+   stopifnot(all(betaTN[isBB] > 1/2))
+ 
+   betaTNs[[flavor]] <- betaTN
+ }
> 
> # Plot
> layout(matrix(1:4, ncol=1))
> par(mar=c(2.5,4,0.5,1)+0.1)
> ylim <- c(-0.05, 1.05)
> col <- rep("#999999", length(muN))
> col[muN == 1/2] <- "#000000"
> for (flavor in names(betaTNs)) {
+   betaTN <- betaTNs[[flavor]]
+   ylab <- sprintf("betaTN[%s]", flavor)
+   plot(pos, betaTN, col=col, ylim=ylim, ylab=ylab)
+ }
> 
> proc.time()
   user  system elapsed 
   0.75    0.14    0.87 

R version 4.0.5 (2021-03-31) -- "Shake and Throw"
Copyright (C) 2021 The R Foundation for Statistical Computing
Platform: x86_64-w64-mingw32/x64 (64-bit)

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.

> library("aroma.light")
aroma.light v3.20.0 successfully loaded. See ?aroma.light for help.
> library("R.utils")
Loading required package: R.oo
Loading required package: R.methodsS3
R.methodsS3 v1.8.1 (2020-08-26 16:20:06 UTC) successfully loaded. See ?R.methodsS3 for help.
R.oo v1.24.0 (2020-08-26 16:11:58 UTC) successfully loaded. See ?R.oo for help.

Attaching package: 'R.oo'

The following object is masked from 'package:R.methodsS3':

    throw

The following objects are masked from 'package:methods':

    getClasses, getMethods

The following objects are masked from 'package:base':

    attach, detach, load, save

R.utils v2.10.1 (2020-08-26 22:50:31 UTC) successfully loaded. See ?R.utils for help.

Attaching package: 'R.utils'

The following object is masked from 'package:utils':

    timestamp

The following objects are masked from 'package:base':

    cat, commandArgs, getOption, inherits, isOpen, nullfile, parse,
    warnings

> 
> # Load data
> pathname <- system.file("data-ex/TumorBoost,fracB,exampleData.Rbin", package="aroma.light")
> data <- loadObject(pathname)
> 
> # Drop loci with missing values
> data <- na.omit(data)
> 
> attachLocally(data)
> pos <- position/1e6
> 
> # Call naive genotypes
> muN <- callNaiveGenotypes(betaN)
> 
> # Genotype classes
> isAA <- (muN == 0)
> isAB <- (muN == 1/2)
> isBB <- (muN == 1)
> 
> # Sanity checks
> stopifnot(all(muN[isAA] == 0))
> stopifnot(all(muN[isAB] == 1/2))
> stopifnot(all(muN[isBB] == 1))
> 
> # TumorBoost normalization with different flavors
> betaTNs <- list()
> for (flavor in c("v1", "v2", "v3", "v4")) {
+   betaTN <- normalizeTumorBoost(betaT=betaT, betaN=betaN, preserveScale=FALSE, flavor=flavor)
+ 
+   # Assert that no non-finite values are introduced
+   stopifnot(all(is.finite(betaTN)))
+ 
+   # Assert that nothing is flipped
+   stopifnot(all(betaTN[isAA] < 1/2))
+   stopifnot(all(betaTN[isBB] > 1/2))
+ 
+   betaTNs[[flavor]] <- betaTN
+ }
> 
> # Plot
> layout(matrix(1:4, ncol=1))
> par(mar=c(2.5,4,0.5,1)+0.1)
> ylim <- c(-0.05, 1.05)
> col <- rep("#999999", length(muN))
> col[muN == 1/2] <- "#000000"
> for (flavor in names(betaTNs)) {
+   betaTN <- betaTNs[[flavor]]
+   ylab <- sprintf("betaTN[%s]", flavor)
+   plot(pos, betaTN, col=col, ylim=ylim, ylab=ylab)
+ }
> 
> proc.time()
   user  system elapsed 
   0.79    0.01    0.79 

R version 4.0.5 (2021-03-31) -- "Shake and Throw"
Copyright (C) 2021 The R Foundation for Statistical Computing
Platform: i386-w64-mingw32/i386 (32-bit)

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.

> library("aroma.light")
aroma.light v3.20.0 successfully loaded. See ?aroma.light for help.
> 
> data(cars)
> attach(cars)
> plot(speed, dist, main = "data(cars)  &  robust smoothing splines")
> 
> # Fit a smoothing spline using L_2 norm
> cars.spl <- smooth.spline(speed, dist)
> lines(cars.spl, col = "blue")
> 
> # Fit a smoothing spline using L_1 norm
> cars.rspl <- robustSmoothSpline(speed, dist)
> lines(cars.rspl, col = "red")
> 
> # Fit a smoothing spline using L_2 norm with 10 degrees of freedom
> lines(smooth.spline(speed, dist, df=10), lty=2, col = "blue")
> 
> # Fit a smoothing spline using L_1 norm with 10 degrees of freedom
> lines(robustSmoothSpline(speed, dist, df=10), lty=2, col = "red")
> 
> # Fit a smoothing spline using Tukey's biweight norm
> cars.rspl <- robustSmoothSpline(speed, dist, method = "symmetric")
> lines(cars.rspl, col = "purple")
> 
> legend(5,120, c(
+       paste("smooth.spline [C.V.] => df =",round(cars.spl$df,1)),
+       paste("robustSmoothSpline L1 [C.V.] => df =",round(cars.rspl$df,1)),
+       paste("robustSmoothSpline symmetric [C.V.] => df =",round(cars.rspl$df,1)),
+       "standard with s( * , df = 10)", "robust with s( * , df = 10)"
+     ),
+     col = c("blue","red","purple","blue","red"), lty = c(1,1,1,2,2),
+     bg='bisque')
> 
> proc.time()
   user  system elapsed 
   0.59    0.01    0.59 

R version 4.0.5 (2021-03-31) -- "Shake and Throw"
Copyright (C) 2021 The R Foundation for Statistical Computing
Platform: x86_64-w64-mingw32/x64 (64-bit)

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.

> library("aroma.light")
aroma.light v3.20.0 successfully loaded. See ?aroma.light for help.
> 
> data(cars)
> attach(cars)
> plot(speed, dist, main = "data(cars)  &  robust smoothing splines")
> 
> # Fit a smoothing spline using L_2 norm
> cars.spl <- smooth.spline(speed, dist)
> lines(cars.spl, col = "blue")
> 
> # Fit a smoothing spline using L_1 norm
> cars.rspl <- robustSmoothSpline(speed, dist)
> lines(cars.rspl, col = "red")
> 
> # Fit a smoothing spline using L_2 norm with 10 degrees of freedom
> lines(smooth.spline(speed, dist, df=10), lty=2, col = "blue")
> 
> # Fit a smoothing spline using L_1 norm with 10 degrees of freedom
> lines(robustSmoothSpline(speed, dist, df=10), lty=2, col = "red")
> 
> # Fit a smoothing spline using Tukey's biweight norm
> cars.rspl <- robustSmoothSpline(speed, dist, method = "symmetric")
> lines(cars.rspl, col = "purple")
> 
> legend(5,120, c(
+       paste("smooth.spline [C.V.] => df =",round(cars.spl$df,1)),
+       paste("robustSmoothSpline L1 [C.V.] => df =",round(cars.rspl$df,1)),
+       paste("robustSmoothSpline symmetric [C.V.] => df =",round(cars.rspl$df,1)),
+       "standard with s( * , df = 10)", "robust with s( * , df = 10)"
+     ),
+     col = c("blue","red","purple","blue","red"), lty = c(1,1,1,2,2),
+     bg='bisque')
> 
> proc.time()
   user  system elapsed 
   0.42    0.09    0.51 

R version 4.0.5 (2021-03-31) -- "Shake and Throw"
Copyright (C) 2021 The R Foundation for Statistical Computing
Platform: i386-w64-mingw32/i386 (32-bit)

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.

> library("aroma.light")
aroma.light v3.20.0 successfully loaded. See ?aroma.light for help.
> 
> X <- matrix(1:30, nrow=5L, ncol=6L)
> mu <- rowMeans(X)
> sd <- apply(X, MARGIN=1L, FUN=sd)
> 
> y <- rowAverages(X)
> stopifnot(all(y == mu))
> stopifnot(all(attr(y,"deviance") == sd))
> stopifnot(all(attr(y,"df") == ncol(X)))
> 
> proc.time()
   user  system elapsed 
   0.31    0.04    0.34 

R version 4.0.5 (2021-03-31) -- "Shake and Throw"
Copyright (C) 2021 The R Foundation for Statistical Computing
Platform: x86_64-w64-mingw32/x64 (64-bit)

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.

> library("aroma.light")
aroma.light v3.20.0 successfully loaded. See ?aroma.light for help.
> 
> X <- matrix(1:30, nrow=5L, ncol=6L)
> mu <- rowMeans(X)
> sd <- apply(X, MARGIN=1L, FUN=sd)
> 
> y <- rowAverages(X)
> stopifnot(all(y == mu))
> stopifnot(all(attr(y,"deviance") == sd))
> stopifnot(all(attr(y,"df") == ncol(X)))
> 
> proc.time()
   user  system elapsed 
   0.34    0.06    0.40 

R version 4.0.5 (2021-03-31) -- "Shake and Throw"
Copyright (C) 2021 The R Foundation for Statistical Computing
Platform: i386-w64-mingw32/i386 (32-bit)

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.

> library("aroma.light")
aroma.light v3.20.0 successfully loaded. See ?aroma.light for help.
> 
> # Simulate 20000 genes with 10 observations each
> X <- matrix(rnorm(n=20000), ncol=10)
> 
> # Calculate the correlation for 5000 random gene pairs
> cor <- sampleCorrelations(X, npairs=5000)
> print(summary(cor))
      Min.    1st Qu.     Median       Mean    3rd Qu.       Max. 
-0.9444409 -0.2545422  0.0089420  0.0009194  0.2464251  0.9084933 
> 
> 
> proc.time()
   user  system elapsed 
   0.57    0.03    0.59 

R version 4.0.5 (2021-03-31) -- "Shake and Throw"
Copyright (C) 2021 The R Foundation for Statistical Computing
Platform: x86_64-w64-mingw32/x64 (64-bit)

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.

> library("aroma.light")
aroma.light v3.20.0 successfully loaded. See ?aroma.light for help.
> 
> # Simulate 20000 genes with 10 observations each
> X <- matrix(rnorm(n=20000), ncol=10)
> 
> # Calculate the correlation for 5000 random gene pairs
> cor <- sampleCorrelations(X, npairs=5000)
> print(summary(cor))
      Min.    1st Qu.     Median       Mean    3rd Qu.       Max. 
-0.8653253 -0.2512924  0.0008791 -0.0021794  0.2377650  0.9026140 
> 
> 
> proc.time()
   user  system elapsed 
   0.53    0.07    0.59 

R version 4.0.5 (2021-03-31) -- "Shake and Throw"
Copyright (C) 2021 The R Foundation for Statistical Computing
Platform: i386-w64-mingw32/i386 (32-bit)

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.

> library("aroma.light")
aroma.light v3.20.0 successfully loaded. See ?aroma.light for help.
> 
> pairs <- sampleTuples(1:10, size=5, length=2)
> print(pairs)
     [,1] [,2]
[1,]    6   10
[2,]    7    8
[3,]    7    1
[4,]    8    5
[5,]    6    5
> 
> triples <- sampleTuples(1:10, size=5, length=3)
> print(triples)
     [,1] [,2] [,3]
[1,]    9    1    8
[2,]    5    3    6
[3,]    6    4    9
[4,]    9    6    5
[5,]    9    3    2
> 
> # Allow tuples with repeated elements
> quadruples <- sampleTuples(1:3, size=5, length=4, replace=TRUE)
> print(quadruples)
     [,1] [,2] [,3] [,4]
[1,]    3    2    2    1
[2,]    1    2    1    3
[3,]    1    3    2    3
[4,]    1    2    2    3
[5,]    3    3    3    1
> 
> proc.time()
   user  system elapsed 
   0.37    0.01    0.37 

R version 4.0.5 (2021-03-31) -- "Shake and Throw"
Copyright (C) 2021 The R Foundation for Statistical Computing
Platform: x86_64-w64-mingw32/x64 (64-bit)

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.

> library("aroma.light")
aroma.light v3.20.0 successfully loaded. See ?aroma.light for help.
> 
> pairs <- sampleTuples(1:10, size=5, length=2)
> print(pairs)
     [,1] [,2]
[1,]    5    9
[2,]    3    5
[3,]    6    2
[4,]    2    8
[5,]    3    2
> 
> triples <- sampleTuples(1:10, size=5, length=3)
> print(triples)
     [,1] [,2] [,3]
[1,]    4    2    8
[2,]    4    3    9
[3,]   10    2    8
[4,]    3    9    6
[5,]    2    9    4
> 
> # Allow tuples with repeated elements
> quadruples <- sampleTuples(1:3, size=5, length=4, replace=TRUE)
> print(quadruples)
     [,1] [,2] [,3] [,4]
[1,]    2    3    3    2
[2,]    2    1    3    2
[3,]    2    2    3    3
[4,]    2    1    1    3
[5,]    3    3    1    3
> 
> proc.time()
   user  system elapsed 
   0.26    0.07    0.32 

R version 4.0.5 (2021-03-31) -- "Shake and Throw"
Copyright (C) 2021 The R Foundation for Statistical Computing
Platform: i386-w64-mingw32/i386 (32-bit)

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.

> library("aroma.light")
aroma.light v3.20.0 successfully loaded. See ?aroma.light for help.
> 
> for (zzz in 0) {
+ 
+ # This example requires plot3d() in R.basic [http://www.braju.com/R/]
+ if (!require(pkgName <- "R.basic", character.only=TRUE)) break
+ 
+ # -------------------------------------------------------------
+ # A first example
+ # -------------------------------------------------------------
+ # Simulate data from the model y <- a + bx + eps(bx)
+ x <- rexp(1000)
+ a <- c(2,15,3)
+ b <- c(2,3,15)
+ bx <- outer(b,x)
+ eps <- apply(bx, MARGIN=2, FUN=function(x) rnorm(length(x), mean=0, sd=0.1*x))
+ y <- a + bx + eps
+ y <- t(y)
+ 
+ # Add some outliers by permuting the dimensions for 1/3 of the observations
+ idx <- sample(1:nrow(y), size=1/3*nrow(y))
+ y[idx,] <- y[idx,c(2,3,1)]
+ 
+ # Down-weight the outliers W times to demonstrate how weights are used
+ W <- 10
+ 
+ # Plot the data with fitted lines at four different view points
+ N <- 4
+ theta <- seq(0,180,length.out=N)
+ phi <- rep(30, length.out=N)
+ 
+ # Use a different color for each set of weights
+ col <- topo.colors(W)
+ 
+ opar <- par(mar=c(1,1,1,1)+0.1)
+ layout(matrix(1:N, nrow=2, byrow=TRUE))
+ for (kk in seq(theta)) {
+   # Plot the data
+   plot3d(y, theta=theta[kk], phi=phi[kk])
+ 
+   # First, same weights for all observations
+   w <- rep(1, length=nrow(y))
+ 
+   for (ww in 1:W) {
+     # Fit a line using IWPCA through data
+     fit <- wpca(y, w=w, swapDirections=TRUE)
+ 
+     # Get the first principal component
+     ymid <- fit$xMean
+     d0 <- apply(y, MARGIN=2, FUN=min) - ymid
+     d1 <- apply(y, MARGIN=2, FUN=max) - ymid
+     b <- fit$vt[1,]
+     y0 <- -b * max(abs(d0))
+     y1 <-  b * max(abs(d1))
+     yline <- matrix(c(y0,y1), nrow=length(b), ncol=2)
+     yline <- yline + ymid
+ 
+     points3d(t(ymid), col=col)
+     lines3d(t(yline), col=col)
+ 
+     # Down-weight outliers only, because here we know which they are.
+     w[idx] <- w[idx]/2
+   }
+ 
+   # Highlight the last one
+   lines3d(t(yline), col="red", lwd=3)
+ }
+ 
+ par(opar)
+ 
+ } # for (zzz in 0)
Loading required package: R.basic
Warning message:
In library(package, lib.loc = lib.loc, character.only = TRUE, logical.return = TRUE,  :
  there is no package called 'R.basic'
> rm(zzz)
> 
> proc.time()
   user  system elapsed 
   0.51    0.09    0.57 

R version 4.0.5 (2021-03-31) -- "Shake and Throw"
Copyright (C) 2021 The R Foundation for Statistical Computing
Platform: x86_64-w64-mingw32/x64 (64-bit)

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.

> library("aroma.light")
aroma.light v3.20.0 successfully loaded. See ?aroma.light for help.
> 
> for (zzz in 0) {
+ 
+ # This example requires plot3d() in R.basic [http://www.braju.com/R/]
+ if (!require(pkgName <- "R.basic", character.only=TRUE)) break
+ 
+ # -------------------------------------------------------------
+ # A first example
+ # -------------------------------------------------------------
+ # Simulate data from the model y <- a + bx + eps(bx)
+ x <- rexp(1000)
+ a <- c(2,15,3)
+ b <- c(2,3,15)
+ bx <- outer(b,x)
+ eps <- apply(bx, MARGIN=2, FUN=function(x) rnorm(length(x), mean=0, sd=0.1*x))
+ y <- a + bx + eps
+ y <- t(y)
+ 
+ # Add some outliers by permuting the dimensions for 1/3 of the observations
+ idx <- sample(1:nrow(y), size=1/3*nrow(y))
+ y[idx,] <- y[idx,c(2,3,1)]
+ 
+ # Down-weight the outliers W times to demonstrate how weights are used
+ W <- 10
+ 
+ # Plot the data with fitted lines at four different view points
+ N <- 4
+ theta <- seq(0,180,length.out=N)
+ phi <- rep(30, length.out=N)
+ 
+ # Use a different color for each set of weights
+ col <- topo.colors(W)
+ 
+ opar <- par(mar=c(1,1,1,1)+0.1)
+ layout(matrix(1:N, nrow=2, byrow=TRUE))
+ for (kk in seq(theta)) {
+   # Plot the data
+   plot3d(y, theta=theta[kk], phi=phi[kk])
+ 
+   # First, same weights for all observations
+   w <- rep(1, length=nrow(y))
+ 
+   for (ww in 1:W) {
+     # Fit a line using IWPCA through data
+     fit <- wpca(y, w=w, swapDirections=TRUE)
+ 
+     # Get the first principal component
+     ymid <- fit$xMean
+     d0 <- apply(y, MARGIN=2, FUN=min) - ymid
+     d1 <- apply(y, MARGIN=2, FUN=max) - ymid
+     b <- fit$vt[1,]
+     y0 <- -b * max(abs(d0))
+     y1 <-  b * max(abs(d1))
+     yline <- matrix(c(y0,y1), nrow=length(b), ncol=2)
+     yline <- yline + ymid
+ 
+     points3d(t(ymid), col=col)
+     lines3d(t(yline), col=col)
+ 
+     # Down-weight outliers only, because here we know which they are.
+     w[idx] <- w[idx]/2
+   }
+ 
+   # Highlight the last one
+   lines3d(t(yline), col="red", lwd=3)
+ }
+ 
+ par(opar)
+ 
+ } # for (zzz in 0)
Loading required package: R.basic
Warning message:
In library(package, lib.loc = lib.loc, character.only = TRUE, logical.return = TRUE,  :
  there is no package called 'R.basic'
> rm(zzz)
> 
> proc.time()
   user  system elapsed 
   0.50    0.03    0.53 

R version 4.0.5 (2021-03-31) -- "Shake and Throw"
Copyright (C) 2021 The R Foundation for Statistical Computing
Platform: i386-w64-mingw32/i386 (32-bit)

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.

> library("aroma.light")
aroma.light v3.20.0 successfully loaded. See ?aroma.light for help.
> 
> # -------------------------------------------------------------
> # A second example
> # -------------------------------------------------------------
> # Data
> x <- c(1,2,3,4,5)
> y <- c(2,4,3,3,6)
> 
> opar <- par(bty="L")
> opalette <- palette(c("blue", "red", "black"))
> xlim <- ylim <- c(0,6)
> 
> # Plot the data and the center mass
> plot(x,y, pch=16, cex=1.5, xlim=xlim, ylim=ylim)
> points(mean(x), mean(y), cex=2, lwd=2, col="blue")
> 
> 
> # Linear regression y ~ x
> fit <- lm(y ~ x)
> abline(fit, lty=1, col=1)
> 
> # Linear regression y ~ x through without intercept
> fit <- lm(y ~ x - 1)
> abline(fit, lty=2, col=1)
> 
> 
> # Linear regression x ~ y
> fit <- lm(x ~ y)
> c <- coefficients(fit)
> b <- 1/c[2]
> a <- -b*c[1]
> abline(a=a, b=b, lty=1, col=2)
> 
> # Linear regression x ~ y through without intercept
> fit <- lm(x ~ y - 1)
> b <- 1/coefficients(fit)
> abline(a=0, b=b, lty=2, col=2)
> 
> 
> # Orthogonal linear "regression"
> fit <- wpca(cbind(x,y))
> 
> b <- fit$vt[1,2]/fit$vt[1,1]
> a <- fit$xMean[2]-b*fit$xMean[1]
> abline(a=a, b=b, lwd=2, col=3)
> 
> # Orthogonal linear "regression" without intercept
> fit <- wpca(cbind(x,y), center=FALSE)
> b <- fit$vt[1,2]/fit$vt[1,1]
> a <- fit$xMean[2]-b*fit$xMean[1]
> abline(a=a, b=b, lty=2, lwd=2, col=3)
> 
> legend(xlim[1],ylim[2], legend=c("lm(y~x)", "lm(y~x-1)", "lm(x~y)",
+           "lm(x~y-1)", "pca", "pca w/o intercept"), lty=rep(1:2,3),
+                      lwd=rep(c(1,1,2),each=2), col=rep(1:3,each=2))
> 
> palette(opalette)
> par(opar)
> 
> proc.time()
   user  system elapsed 
   0.51    0.06    0.54 

R version 4.0.5 (2021-03-31) -- "Shake and Throw"
Copyright (C) 2021 The R Foundation for Statistical Computing
Platform: x86_64-w64-mingw32/x64 (64-bit)

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.

> library("aroma.light")
aroma.light v3.20.0 successfully loaded. See ?aroma.light for help.
> 
> # -------------------------------------------------------------
> # A second example
> # -------------------------------------------------------------
> # Data
> x <- c(1,2,3,4,5)
> y <- c(2,4,3,3,6)
> 
> opar <- par(bty="L")
> opalette <- palette(c("blue", "red", "black"))
> xlim <- ylim <- c(0,6)
> 
> # Plot the data and the center mass
> plot(x,y, pch=16, cex=1.5, xlim=xlim, ylim=ylim)
> points(mean(x), mean(y), cex=2, lwd=2, col="blue")
> 
> 
> # Linear regression y ~ x
> fit <- lm(y ~ x)
> abline(fit, lty=1, col=1)
> 
> # Linear regression y ~ x through without intercept
> fit <- lm(y ~ x - 1)
> abline(fit, lty=2, col=1)
> 
> 
> # Linear regression x ~ y
> fit <- lm(x ~ y)
> c <- coefficients(fit)
> b <- 1/c[2]
> a <- -b*c[1]
> abline(a=a, b=b, lty=1, col=2)
> 
> # Linear regression x ~ y through without intercept
> fit <- lm(x ~ y - 1)
> b <- 1/coefficients(fit)
> abline(a=0, b=b, lty=2, col=2)
> 
> 
> # Orthogonal linear "regression"
> fit <- wpca(cbind(x,y))
> 
> b <- fit$vt[1,2]/fit$vt[1,1]
> a <- fit$xMean[2]-b*fit$xMean[1]
> abline(a=a, b=b, lwd=2, col=3)
> 
> # Orthogonal linear "regression" without intercept
> fit <- wpca(cbind(x,y), center=FALSE)
> b <- fit$vt[1,2]/fit$vt[1,1]
> a <- fit$xMean[2]-b*fit$xMean[1]
> abline(a=a, b=b, lty=2, lwd=2, col=3)
> 
> legend(xlim[1],ylim[2], legend=c("lm(y~x)", "lm(y~x-1)", "lm(x~y)",
+           "lm(x~y-1)", "pca", "pca w/o intercept"), lty=rep(1:2,3),
+                      lwd=rep(c(1,1,2),each=2), col=rep(1:3,each=2))
> 
> palette(opalette)
> par(opar)
> 
> proc.time()
   user  system elapsed 
   0.54    0.03    0.56