ENH: more complete options for simlr.search

R/multiscaleSVDxpts.R

@@ -3178,13 +3178,20 @@ simlr <- function(
     jointInitialization = TRUE,
     sparsenessAlg = NA,
     verbose = FALSE) {
+  parse_constraint <- function(x) {
+    num1=num2=NA
+    temp=unlist(strsplit( x, "x"))
+    if ( length(temp) > 1 ) num1=as.numeric(temp[2])
+    if ( length(temp) > 2 ) num2=as.numeric(temp[3])
+    return(c(temp[1], as.numeric(num1), as.numeric(num2)))
+    }
   if (missing(scale)) scale <- c("centerAndScale")
   if (missing(energyType)) energyType <- "cca"
   if (missing(mixAlg)) mixAlg <- "svd"
   if (missing(optimizationStyle)) optimizationStyle <- "lineSearch"
   if (!missing("randomSeed")) set.seed(randomSeed) #  else set.seed( 0 )
   energyType <- match.arg(energyType)
-  constraint <- match.arg(constraint)
+  constraint <- parse_constraint( constraint[1] )
   optimizationStyle <- match.arg(optimizationStyle)
   scalechoices = c(
       "sqrtnp", "np", "centerAndScale",
@@ -3402,7 +3409,8 @@ simlr <- function(
   nc <- ncol(initialUMatrix[[1]])
   myw <- matrix(rnorm(nc^2), nc, nc) # initialization for fastICA
   getSyME2 <- function(lineSearch, gradient, myw, mixAlg,
-                       avgU, whichModality, last_energy=0, verbose = FALSE ) {
+                       avgU, whichModality, last_energy=0, 
+                       constraint=c('ortho',0.5,2.0), verbose = FALSE ) {
     prediction <- 0
     myenergysearchv <- (vmats[[whichModality]] + gradient * lineSearch) # update the i^th v matrix
     if (verbose) {
@@ -3420,9 +3428,11 @@ simlr <- function(
         sparsenessAlg = sparsenessAlg
       )
     }
-    myorthEnergy = invariant_orthogonality_defect( myenergysearchv )
-    if ( last_energy > .Machine$double.eps & myorthEnergy > .Machine$double.eps)
-      myorthEnergy = 0.5 * myorthEnergy*(last_energy/myorthEnergy)
+    if ( constraint[1] == 'ortho' ) {
+      myorthEnergy = invariant_orthogonality_defect( myenergysearchv )
+      if ( last_energy > .Machine$double.eps & myorthEnergy > .Machine$double.eps)
+        myorthEnergy = as.numeric(constraint[2]) * myorthEnergy*(last_energy/myorthEnergy)
+      } else myorthEnergy = 0.0
     if (ccaEnergy) {
       # ( v'*X'*Y )/( norm2(X*v ) * norm2( u ) )
       t0 <- norm(voxmats[[whichModality]] %*% myenergysearchv, "F")
@@ -3471,7 +3481,8 @@ simlr <- function(
     loki <- getSyME2(0, 0,
       myw = myw, mixAlg = mixAlg,
       avgU = initialUMatrix[[i]],
-      whichModality = i
+      whichModality = i,
+      constraint=constraint
     )
     initialEnergy <- initialEnergy + loki / nModalities
   }
@@ -3640,7 +3651,7 @@ simlr <- function(
       }
       # initialize gradient line search
       temperv <- getSyMG(vmats[[i]], i, myw = myw, mixAlg = mixAlg)
-      temperv <- constrainG(temperv, i, constraint = constraint)
+      temperv <- constrainG(temperv, i, constraint = constraint[1] )
 
       useAdam <- FALSE
       if (useAdam) { 
@@ -3666,16 +3677,19 @@ simlr <- function(
         temperv <- temperv * (1.0 - expBeta) + lastG[[i]] * (expBeta)
         lastG[[i]] <- temperv
       }
-      orthgrad = gradient_invariant_orthogonality_defect( vmats[[i]] )
-      temperv = temperv - orthgrad * norm(orthgrad,"F")/norm(temperv,"F")*2.0
+      if ( constraint[1] == 'ortho' ) {
+        orthgrad = gradient_invariant_orthogonality_defect( vmats[[i]] )
+        temperv = temperv - orthgrad * norm(orthgrad,"F")/norm(temperv,"F")*as.numeric(constraint[3])
+      }
       if ( myit > 1 ) laste = energyPath[ myit - 1 ] else laste = 1e9
       if (optimizationLogic(energyPath, myit, i)) {
         temp <- optimize(getSyME2, # computes the energy
           interval = lineSearchRange,
           tol = lineSearchTolerance,
           gradient = temperv,
           myw = myw, mixAlg = mixAlg,
-          avgU = initialUMatrix[[i]], whichModality = i, last_energy=laste
+          avgU = initialUMatrix[[i]], whichModality = i, 
+          last_energy=laste, constraint=constraint
         )
         errterm[i] <- temp$objective
         gamma[i] <- temp$minimum
@@ -3685,7 +3699,8 @@ simlr <- function(
           gamma[i], temperv,
           myw = myw, mixAlg = mixAlg,
           avgU = initialUMatrix[[i]],
-          whichModality = i, last_energy=laste
+          whichModality = i, last_energy=laste,
+          constraint=constraint
         )
       }
       if (errterm[i] <= min(energyPath[, i], na.rm = T) |
@@ -3737,7 +3752,9 @@ simlr <- function(
       loki <- getSyME2(0, 0,
         myw = myw, mixAlg = mixAlg,
         avgU = initialUMatrix[[jj]],
-        whichModality = jj, verbose = FALSE
+        whichModality = jj, 
+        constraint=constraint,
+        verbose = FALSE
       )
       energyPath[myit, jj] <- loki
     } # matrix loop
@@ -4411,6 +4428,8 @@ simlr.parameters <- function(
 #' @param regs The regularization options for SIMLR.
 #' @param options_df A data frame of parameter combinations generated by `simlr.parameters`.
 #' @param maxits The maximum number of iterations for SIMLR. Default is 100.
+#' @param connectors a list ( length of projections or number of modalities )
+#' that indicates which modalities should be paired with current modality
 #' @param nperms The number of permutations for the significance test. Default is 1.
 #' @param verbose The verbosity level. Default is 0.
 #' @param FUN The function to use for the SIMLR evaluation. Default is `rvcoef`.
@@ -4421,6 +4440,7 @@ simlr.search <- function(
   regs,
   options_df,
   maxits = 100,
+  connectors = NULL,
   nperms = 1,
   verbose = 0,
   FUN = rvcoef
@@ -4500,7 +4520,7 @@ simlr.search <- function(
       optimizationStyle = optimus,
       initialUMatrix = if ("lowrank" %in% prescaling) lowrankRowMatrix(initu, nsimlr * 2) else initu,
       constraint = constraint,
-      connectors = simlr_path_models(length(mats), 0),
+      connectors = connectors,
       verbose = verbose > 2,
       nperms = nperms,
       FUN = FUN