PRIA coverage - 39.82%

Files
Source

#' Title
#'
#' @param x TBD
#' @param ... TBD
#'
#' @return TBD
#' @export
"print.pria_table" <- function(x, ...) {
  print(data.table::as.data.table(x), ...);
}

#' Title
#'
#' @param x TBD
#' @param pattern TBD
#'
#' @return TBD
#' @export
".DollarNames.pria_table" <- function(x, pattern = "") {
  columns <- colnames(x);
  choose_from <- c(columns, c("codes.maintained", "codes.removed", "ena.set"));
  grep(pattern = pattern, x = choose_from, value = TRUE)
}

#' Title
#'
#' @param x  TBD
#' @param i TBD
#' @param j  TBD
#' @param ...  TBD
#' @param drop  TBD
#'
#' @return TBD
#' @export
"[.pria_table" <- function(x, i, j, ..., drop = TRUE) {
  x_ <- data.table::copy(x);
  class(x_) <- class(x_)[!class(x_) %chin% "pria_table"]

  pria_row <- x_[i, ]

  class(pria_row) <- c("pria_table", class(pria_row));

  pria_row;
}

#' Override element extraction
#'
#' @param x TBD
#' @param i TBD
#' @param j TBD
#'
#' @return TBD
#' @export
"$.pria_table" <- function(x, i, j) {
  if(i %in% colnames(x)) {
    x[[i]]
  }
  else {
    res <- attr(x = x, which = i, exact = TRUE)
    res
  }
}

#' Title
#'
#' @param e1 TBD
#' @param e2 TBD
#'
#' @return TBD
#' @export
">.pria_table" <- function(e1, e2) {
  reduced <- e1[
    e1$goodness.of.fit.ci_lower > e2 &
    e1$points.ci_lower > e2 &
    e1$nodes.ci_lower > e2,
  ]

  reduced <- reduced[order(-reduced$num.removed.codes, -reduced$goodness.of.fit),]
  class(reduced) <- c("pria_table", class(reduced));

  reduced
}

#' Title
#'
#' @param ... TBD
#' @param na.rm TBD
#'
#' @return TBD
#' @export
"max.pria_table" <- function(..., na.rm = FALSE) {
  all <- list(...)[[1]];
  best_index <- order(-all$num.removed.codes, -all$goodness.of.fit)[1]
  best <- all[best_index,][1,]

  best
}

#' Title
#'
#' @param x  TBD
#' @param i TBD
#' @param ...  TBD
#'
#' @return TBD
#' @export
"[.ena.sets" <- function(x, i, ...) {
  x_ <- x;
  class(x_) <- class(x_)[!class(x_) %chin% "ena.sets"]

  res <- `[`(x_, i, ...);
  class(res) <- c("ena.sets", class(res))
  res
}

#' Title
#'
#' @param x TBD
#' @param pattern TBD
#'
#' @return TBD
#' @export
".DollarNames.ena.sets" <- function(x, pattern = "") {
  grep(pattern = pattern, x = names(x[[1]]), value = TRUE)
}

#' Override element extraction
#'
#' @param x TBD
#' @param i TBD
#' @param j TBD
#'
#' @return TBD
#' @export
"$.ena.sets" <- function(x, i, j) {
  res <- lapply(x, `[[`, i)
  class(res) <- c("ena.sets", class(res))
  res
}

#' Title
#'
#' @param e1 TBD
#' @param e2 TBD
#'
#' @return TBD
#' @export
"==.ena.sets" <- function(e1,e2) {
  lapply(e1, `==`, e2)
}

#' Internal Performance Calculation
#' @noRd
calc_perf_internal <- function(result, min_gof) {
  gof <- result$goodness.of.fit
  if (gof < min_gof) gof <- 0
  #return(gof * max(0, result$points) * max(0, result$nodes))
  return(gof * result$points * result$nodes)
}

#' Internal Model Scorer
#' @noRd
get_score_internal <- function(current_codes, x, correlation_fun, min_gof,use_dist_corr=TRUE) {
  full_codes <- x$rotation$codes
  removed <- setdiff(full_codes, current_codes)

  if (length(removed) == 0) {
    reduced.set <- x
  } else {
    reduced.accum <- remove.codes.from.accum(x, removed)
    if (is.matrix(reduced.accum$rotation$adjacency.key)) {
      wh <- if (is.null(x$`_function.params`$rotation.by)) rENA::ena.svd else x$`_function.params`$rotation.by
      reduced.set <- reduced.accum |>
        rENA::sphere_norm() |>
        rENA::center() |>
        rENA::rotate(wh = wh) |>
        rENA::project() |>
        rENA::optimize()
      reduced.set$rotation$adjacency.key <- reduced.accum$rotation$adjacency.key
    } else {
      reduced.set <- rENA::ena.make.set(
        reduced.accum,
        rotation.by = x$`_function.params`$rotation.by,
        rotation.params = x$`_function.params`$rotation.params
      )
    }
  }

  if (use_dist_corr) {
    res <- compare.two.models.by.dist(x, reduced.set, keep_sets = FALSE)
    return(res$score)
  } else {
    res <- compare.two.models(x, reduced.set, FUN = correlation_fun, keep_sets = FALSE)
  }
  return(calc_perf_internal(res, min_gof))
}
#' Compute all pria scores
#'
#' @param current_codes current codes with some codes removed
#' @param x ENA set with full codes
#' @param correlation_fun correlation function
#' @param min_gof GOF threshold
#' @param grouping_col a grouping columns for effect size computation, default=NULL
#' @param include_reduced_set whether or not include the reduced set in the return,default=FALSE
#' @param n_dim number of first dimensions to use in similarity computations, default=2
#'
#' @return a list of scores, and optionally the reduced ENA set
#' @export
get_pria_scores <- function(current_codes, x, correlation_fun, min_gof,grouping_col=NULL,include_reduced_set=FALSE,n_dim=2) {
  full_codes <- x$rotation$codes
  removed <- setdiff(full_codes, current_codes)
  reduced.set =NULL
  if (length(removed) == 0) {
    reduced.set <- x
  } else {
    reduced.accum <- remove.codes.from.accum(x, removed)
    reduced.set <- rENA::ena.make.set(
      reduced.accum,
      rotation.by = x$`_function.params`$rotation.by,
      rotation.params = x$`_function.params`$rotation.params
    )
  }

  res1<-compare.two.models.by.dist(x,reduced.set,keep_sets=FALSE,n_dim)
  space_sim=res1$score

  res2 <- compare.two.models(x, reduced.set, FUN = correlation_fun, keep_sets = FALSE)
  gof <-res2$goodness.of.fit
  point_corr<-res2$points
  node_corr<-res2$nodes
  eta2 = NULL
  if(!is.null(grouping_col))
  {
    first_dim <- as.matrix(reduced.set$points)[,1]
    cond_vector <- reduced.set$meta.data[[grouping_col]]

    # Run ANOVA
    fit <- aov(first_dim ~ cond_vector)
    sum_stats <- summary(fit)[[1]]

    ss_between <- sum_stats[1, "Sum Sq"]
    ss_residual <- sum_stats[2, "Sum Sq"]
    eta2 <- ss_between / (ss_between + ss_residual)
  }
  if(!include_reduced_set)
  {
    reduced.set=NULL
  }
  return(list(space_sim = space_sim,
              gof = gof,
              point_corr = point_corr,
              node_corr=node_corr,
              eta2=eta2,
              reduced_set=reduced.set))
}
#' Compute pria scores on first two dimensions
#'
#' @param x ENA set with full codes
#' @param removed_codes codes removed
#' @param threshold threshold, default=0.95
#' @param rebuild optional function taking a reduced accumulation and returning
#'   the corresponding ENA set; see \code{\link{pria}}.
#' @return variance explained on dim 1 of the reduced set
#' @export
get_pria_scores_2Ds <- function(x,removed_codes,threshold=0.95,rebuild=NULL) {
  # construct reduced set
  reduced.set =NULL
  if (length(removed_codes) == 0) {
    reduced.set <- x
  } else {
    reduced.accum <- remove.codes.from.accum(x, removed_codes)
    reduced.set <- if (is.function(rebuild)) {
      # Callers that know how the original set was built supply `rebuild`. The
      # default below cannot serve every model: ena.make.set only builds
      # undirected sets, so an ordered/directed accumulation fails outright,
      # and it reads the rotation off `_function.params`, which callers such
      # as the ENA webtool never populate -- silently scoring each candidate
      # against a default-SVD rotation instead of the model's own.
      rebuild(reduced.accum)
    } else {
      rENA::ena.make.set(
        reduced.accum,
        rotation.by = x$`_function.params`$rotation.by,
        rotation.params = x$`_function.params`$rotation.params
      )
    }
  }
  # compute relative variance
  vr1 = reduced.set$model$variance[1]

  # compute gof, original set
  pts0 = x$points
  cts0 = x$model$centroids
  pts20 = as.matrix(pts0)[,1:2,drop=FALSE]
  cts20 = as.matrix(cts0)[,1:2,drop=FALSE]
  cor_result0 <-rENA::ena_correlation(pts20,cts20)
  gof_1_o = cor_result0[1,1]
  gof_2_o = cor_result0[2,1]
  # compute gof, reduced set
  pts = reduced.set$points
  cts = reduced.set$model$centroids
  pts2 = as.matrix(pts)[,1:2,drop=FALSE]
  cts2 = as.matrix(cts)[,1:2,drop=FALSE]
  cor_result1 <-rENA::ena_correlation(pts2,cts2)
  gof_1 = cor_result1[1,1]
  gof_2 = cor_result1[2,1]

  if(min(gof_1/gof_1_o,gof_2/gof_2_o)<threshold)
  {
    return(NULL)
  }

  # compute point, node correlations
  full_points = as.matrix(x$points)[,1:2,drop=FALSE]
  reduce_points = pts2
  cor_result2 = rENA::ena_correlation(full_points,reduce_points)

  full_nodes = x$rotation$nodes
  reduced_nodes =reduced.set$rotation$nodes
  # Keep the full model's rows for the codes the reduced model retained. Align
  # on rotation$codes rather than a `code` column: node tables from some
  # pipelines (the ENA webtool among them) carry only dimension columns, and
  # `full_nodes$code` is then NULL, selecting zero rows and leaving the
  # correlation below to fail with "subscript out of bounds".
  keep_nodes <- if (!is.null(full_nodes$code) && !is.null(reduced_nodes$code)) {
    full_nodes$code %in% reduced_nodes$code
  } else {
    x$rotation$codes %in% reduced.set$rotation$codes
  }
  full_nodes <- full_nodes[keep_nodes, ]
  full_nodes<-as.matrix(full_nodes)[,1:2,drop=FALSE]
  reduced_nodes<-as.matrix(reduced_nodes)[,1:2,drop=FALSE]
  cor_result3 = rENA::ena_correlation(full_nodes,reduced_nodes)

  # check flipping
  point_corr_1 = cor_result2[1,1]
  point_corr_2 = cor_result2[2,1]
  node_corr_1 = cor_result3[1,1]
  node_corr_2 = cor_result3[2,1]
  if(point_corr_1<0)
  {
    point_corr_1 = -point_corr_1
    node_corr_1 = - node_corr_1
  }
  if(point_corr_2<0)
  {
    point_corr_2 = -point_corr_2
    node_corr_2 = - node_corr_2
  }

  if(min(point_corr_1,point_corr_2,node_corr_1,node_corr_2)<threshold)
  {
    return(NULL)
  }

  return(list(
              removed_codes = removed_codes,
              vr1=vr1,
              gof1 = gof_1/gof_1_o,
              gof2 = gof_2/gof_2_o,
              point1 = point_corr_1,
              point2 = point_corr_2,
              node1 = node_corr_1,
              node2 = node_corr_2))
}
#' Compute pria scores on first two dimensions
#'
#' @param current_codes current codes with some codes removed
#' @param x ENA set with full codes
#' @param correlation_fun correlation function
#' @param min_gof GOF threshold
#' @param grouping_col a grouping columns for effect size computation, default=NULL
#' @param include_reduced_set whether or not include the reduced set in the return,default=FALSE
#' @param n_dim number of first dimensions to use in similarity computations, default=2
#'
#' @return a list of scores, and optionally the reduced ENA set
#' @export
get_pria_scores_2D <- function(x,removed_codes ) {
  # construct reduced set
  reduced.set =NULL
  if (length(removed_codes) == 0) {
    reduced.set <- x
  } else {
    reduced.accum <- remove.codes.from.accum(x, removed_codes)
    reduced.set <- rENA::ena.make.set(
      reduced.accum,
      rotation.by = x$`_function.params`$rotation.by,
      rotation.params = x$`_function.params`$rotation.params
    )
  }
  # compute pair-wise distance similarity in 2D
  res1<-compare.two.models.by.dist(x,reduced.set,keep_sets=FALSE,2)

  #***
  space_sim=res1$score
  #***
  #
  # compute gof, original set
  pts0 = x$points
  cts0 = x$model$centroids
  pts20 = as.matrix(pts0)[,1:2,drop=FALSE]
  cts20 = as.matrix(cts0)[,1:2,drop=FALSE]
  cor_result0 <-rENA::ena_correlation(pts20,cts20)
  # compute gof, reduced set
  pts = reduced.set$points
  cts = reduced.set$model$centroids
  pts2 = as.matrix(pts)[,1:2,drop=FALSE]
  cts2 = as.matrix(cts)[,1:2,drop=FALSE]
  cor_result1 <-rENA::ena_correlation(pts2,cts2)

  # compute point, node correlations
  full_points = as.matrix(x$points)[,1:2,drop=FALSE]
  reduce_points = pts2
  cor_result2 = rENA::ena_correlation(full_points,reduce_points)

  full_nodes = x$rotation$nodes
  reduced_nodes =reduced.set$rotation$nodes
  # Keep the full model's rows for the codes the reduced model retained. Align
  # on rotation$codes rather than a `code` column: node tables from some
  # pipelines (the ENA webtool among them) carry only dimension columns, and
  # `full_nodes$code` is then NULL, selecting zero rows and leaving the
  # correlation below to fail with "subscript out of bounds".
  keep_nodes <- if (!is.null(full_nodes$code) && !is.null(reduced_nodes$code)) {
    full_nodes$code %in% reduced_nodes$code
  } else {
    x$rotation$codes %in% reduced.set$rotation$codes
  }
  full_nodes <- full_nodes[keep_nodes, ]
  full_nodes<-as.matrix(full_nodes)[,1:2,drop=FALSE]
  reduced_nodes<-as.matrix(reduced_nodes)[,1:2,drop=FALSE]
  cor_result3 = rENA::ena_correlation(full_nodes,reduced_nodes)


  # check flipping
  point_corr_1 = cor_result2[1,1]
  point_corr_2 = cor_result2[2,1]
  node_corr_1 = cor_result3[1,1]
  node_corr_2 = cor_result3[2,1]
  if(point_corr_1<0)
  {
    point_corr_1 = -point_corr_1
    node_corr_1 = - node_corr_1
  }
  if(point_corr_2<0)
  {
    point_corr_2 = -point_corr_2
    node_corr_2 = - node_corr_2
  }

  # compute cosine between two rotation hyperplanes
  full_rm = x$rotation.matrix
  reduced_rm = reduced.set$rotation.matrix
  # base-R subset (dplyr isn't imported, so `filter` fell through to stats::filter)
  full_rm <- full_rm[full_rm$codes %in% reduced_rm$codes, ]
  full_rm = as.matrix(full_rm)[,2:3,drop=FALSE]
  reduced_rm = as.matrix(reduced_rm)[,2:3,drop=FALSE]
  aa = t(full_rm)%*%reduced_rm
  cos_rotation = svd(aa)$d[2]
  # compute consine between 2D points hyperplanes

    A <- full_points
    B <- reduce_points

    # 2. Get orthonormal bases using QR decomposition
    Q1 <- qr.Q(qr(A))
    Q2 <- qr.Q(qr(B))

    # 3. Compute the inner product matrix
    M <- t(Q1) %*% Q2

    # 4. Get singular values (cosines of principal angles)
    cos_point <- svd(M)$d[2]


  return(list(space_sim = space_sim,
              gof_1_o = cor_result0[1,1],
              gof_2_o = cor_result0[2,1],
              gof_1 = cor_result1[1,1],
              gof_2 = cor_result1[2,1],
              point_corr_1 = point_corr_1,
              point_corr_2 = point_corr_2,
              node_corr_1 = node_corr_1 ,
              node_corr_2 = node_corr_2,
              cos_rotation = cos_rotation,
              cos_point = cos_point
              ))
}

#' Correlations
#'
#' @param enaset TBD
#' @param pts TBD
#' @param cts TBD
#' @param dims TBD
#'
#' @return TBD
#' @export
correlations_slow <- function(enaset, pts = NULL, cts = NULL, dims = c(1:2)) {
  # I rewrote the goodness of fit test since the current package doesn't return CI
  if(is.null(pts) || is.null(cts)) {
    pts <- enaset$points
    cts <- enaset$model$centroids
  }

  pComb = combn(nrow(pts), 2)
  # pComb = rENA::combn_c2(nrow(pts)) + 1
  point1 = pComb[1,]
  point2 = pComb[2,]

  points = as.matrix(pts)
  centroids = as.matrix(cts)

  # svdDiff = matrix(points[point1, dims] - points[point2, dims], ncol=length(dims), nrow=length(point1))
  svdDiff = points[point1, dims, drop = FALSE] - points[point2, dims, drop = FALSE];
  # optDiff = matrix(centroids[point1, dims] - centroids[point2, dims], ncol=length(dims), nrow=length(point1))
  optDiff = centroids[point1, dims, drop = FALSE] - centroids[point2, dims, drop = FALSE]

  cor_result <- cor.test(as.numeric(svdDiff[,1]), as.numeric(optDiff[,1]))
  corr_lower_bound <- cor_result$conf.int[1]
  corr_upper_bound <- cor_result$conf.int[2]
  corr_estimates   <- cor_result$estimate[1]

  corr_results <- list(
    "corr_lower_bound"=corr_lower_bound,
    "corr_upper_bound"=corr_upper_bound,
    "corr_estimates"=as.numeric(corr_estimates)
  )
  return(corr_results);
}

#' Correlations
#'
#' @param enaset TBD
#' @param pts TBD
#' @param cts TBD
#' @param dims TBD
#' @param conf.level TBD
#'
#' @return TBD
#' @export
correlations <- function(enaset, pts = NULL, cts = NULL, dims = c(1), conf.level = 0.95) {
  # I rewrote the goodness of fit test since the current package doesn't return CI
  if(is.null(pts) || is.null(cts)) {
    pts <- enaset$points
    cts <- enaset$model$centroids
  }

  points = as.matrix(pts)[, dims, drop = FALSE];
  centroids = as.matrix(cts)[, dims, drop = FALSE];

  # pc <- matrix(c(points[,1], centroids[,1]), ncol = 2)

  cor_result <- rENA::ena_correlation(points, centroids)

  corr_results <- list(
    "corr_lower_bound" = abs(cor_result[1,2]),
    "corr_upper_bound" = abs(cor_result[1,3]),
    "corr_estimates" = abs(cor_result[1,1])
  )

  return(corr_results);
}

#' Compute GOF and cors
#'
#' @param full.set TBD
#' @param reduced.set TBD
#' @param FUN TBD
#' @param keep_sets TBD
#'
#' @return TBD
#' @export
compare.two.models = function(full.set, reduced.set, FUN = correlations_slow, keep_sets = TRUE){
  ### Arguments:
  # dimension.type: whether you would like to run comparison for MR dimension ("MR") or regular first SVD dimension.
  # full.set (or model 0 set): you need to pass an ena.set object. It could be mean rotation or regular ENA.
  # reduced.set (or model alternative set): you need to pass an ena.set object. It could be mean rotation or regular ENA.

  ### Return:
  # upper or lower boundary of C.I.s for three criteria (1. goodneess of fit, 2. location of points, 3. location of nodes)
  full.points <- as.matrix(full.set$points)[,1]
  full.nodes <- as.matrix(full.set$rotation$nodes)[,1]
  reduced.points <- as.matrix(reduced.set$points)[,1]
  reduced.nodes <- as.matrix(reduced.set$rotation$nodes)[,1]
  full.codes <- full.set$rotation$codes
  reduced.codes <- reduced.set$rotation$codes

  ### 1. Compare goodness of fit
  cor.set.res <- FUN(reduced.set);
  cor.set.gof = round(cor.set.res$corr_estimates, digits = 4)
  lower.CI.cor.set.gof = round(cor.set.res$corr_lower_bound, digits = 4)
  upper.CI.cor.set.gof = round(cor.set.res$corr_upper_bound, digits = 4)

  ### 2. Compare point locations
  cor.points <- cor.test(full.points, reduced.points)
  lower.CI.cor.points <- round(as.numeric(cor.points$conf.int[1]), digits = 4)
  upper.CI.cor.points <- round(as.numeric(cor.points$conf.int[2]), digits = 4)

  ### 3. Compare node positions
  full.nodes.reduced <- full.nodes[full.codes %in% reduced.codes]
  cor.nodes <- cor.test(full.nodes.reduced, reduced.nodes)
  lower.CI.cor.nodes <- round(as.numeric(cor.nodes$conf.int[1]), digits = 4)
  upper.CI.cor.nodes <- round(as.numeric(cor.nodes$conf.int[2]), digits = 4)

  result = list(
    "num.removed.codes" = length(full.codes) - length(reduced.codes),
    "goodness.of.fit" = cor.set.gof,
    "goodness.of.fit.ci_lower" = lower.CI.cor.set.gof,
    "goodness.of.fit.ci_upper" = upper.CI.cor.set.gof,
    "points" = as.numeric(cor.points$estimate),
    "points.ci_lower" = lower.CI.cor.points,
    "points.ci_upper" = upper.CI.cor.points,
    "nodes" = as.numeric(cor.nodes$estimate),
    "nodes.ci_lower" = lower.CI.cor.nodes,
    "nodes.ci_upper" = upper.CI.cor.nodes
    ,"codes.used" = list(reduced.codes)
  )

  attr(result, "codes.removed") <- as.character(reduced.codes);
  attr(result, "codes.maintained") <- as.character(full.codes %in% reduced.codes);
  if (keep_sets == TRUE) {
    attr(result, "ena.set") <- reduced.set;
  }
  # class(result) <- c("pria_result", class(result))
  return(result)
}
#' Compute GOF and cors
#'
#' @param full.set TBD
#' @param reduced.set TBD
#' @param FUN TBD
#' @param keep_sets TBD
#'
#' @return TBD
#' @export
compare.two.models_cors = function(full.set, reduced.set, FUN = correlations, keep_sets = TRUE){
  ### Arguments:
  # dimension.type: whether you would like to run comparison for MR dimension ("MR") or regular first SVD dimension.
  # full.set (or model 0 set): you need to pass an ena.set object. It could be mean rotation or regular ENA.
  # reduced.set (or model alternative set): you need to pass an ena.set object. It could be mean rotation or regular ENA.

  ### Return:
  # upper or lower boundary of C.I.s for three criteria (1. goodneess of fit, 2. location of points, 3. location of nodes)
  full.points <- as.matrix(full.set$points)[,1]
  full.nodes <- as.matrix(full.set$rotation$nodes)[,1]
  reduced.points <- as.matrix(reduced.set$points)[,1]
  reduced.nodes <- as.matrix(reduced.set$rotation$nodes)[,1]
  full.codes <- full.set$rotation$codes
  reduced.codes <- reduced.set$rotation$codes

  ### 1. Compare goodness of fit
  cor.set.res <- FUN(reduced.set);
  cor.set.gof = round(cor.set.res$corr_estimates, digits = 4)
  lower.CI.cor.set.gof = round(cor.set.res$corr_lower_bound, digits = 4)
  upper.CI.cor.set.gof = round(cor.set.res$corr_upper_bound, digits = 4)

  ### 2. Compare point locations
  cor.points <- cor.test(full.points, reduced.points)
  lower.CI.cor.points <- round(as.numeric(cor.points$conf.int[1]), digits = 4)
  upper.CI.cor.points <- round(as.numeric(cor.points$conf.int[2]), digits = 4)

  ### 3. Compare node positions
  full.nodes.reduced <- full.nodes[full.codes %in% reduced.codes]
  cor.nodes <- cor.test(full.nodes.reduced, reduced.nodes)
  lower.CI.cor.nodes <- round(as.numeric(cor.nodes$conf.int[1]), digits = 4)
  upper.CI.cor.nodes <- round(as.numeric(cor.nodes$conf.int[2]), digits = 4)

  result = list(
    "num.removed.codes" = length(full.codes) - length(reduced.codes),
    "goodness.of.fit" = cor.set.gof,
    "goodness.of.fit.ci_lower" = lower.CI.cor.set.gof,
    "goodness.of.fit.ci_upper" = upper.CI.cor.set.gof,
    "points" = as.numeric(cor.points$estimate),
    "points.ci_lower" = lower.CI.cor.points,
    "points.ci_upper" = upper.CI.cor.points,
    "nodes" = as.numeric(cor.nodes$estimate),
    "nodes.ci_lower" = lower.CI.cor.nodes,
    "nodes.ci_upper" = upper.CI.cor.nodes
    ,"codes.used" = list(reduced.codes)
  )

  attr(result, "codes.removed") <- as.character(reduced.codes);
  attr(result, "codes.maintained") <- as.character(full.codes %in% reduced.codes);
  if (keep_sets == TRUE) {
    attr(result, "ena.set") <- reduced.set;
  }
  # class(result) <- c("pria_result", class(result))
  return(result)
}
#' Compare to models by distance
#'
#' @param full.set TBD
#' @param reduced.set TBD
#' @param FUN TBD
#' @param keep_sets TBD
#'
#' @return TBD
#' @export
compare.two.models.by.dist = function(full.set, reduced.set, keep_sets = TRUE,n_dim=2){
  ### Arguments:
  # dimension.type: whether you would like to run comparison for MR dimension ("MR") or regular first SVD dimension.
  # full.set (or model 0 set): you need to pass an ena.set object. It could be mean rotation or regular ENA.
  # reduced.set (or model alternative set): you need to pass an ena.set object. It could be mean rotation or regular ENA.

  ### Return:
  # upper or lower boundary of C.I.s for three criteria (1. goodneess of fit, 2. location of points, 3. location of nodes)
  full.points <- as.matrix(full.set$points)
  reduced.points <- as.matrix(reduced.set$points)[,1:n_dim]
  if(n_dim>0)
  {
    full.points=full.points[,1:min(n_dim,ncol(full.points))]
    reduced.points=reduced.points[,1:min(n_dim,ncol(reduced.points))]
  }

  full.codes <- full.set$rotation$codes
  reduced.codes <- reduced.set$rotation$codes
  dist_corr_score <-dist_corr(full.points,reduced.points)

  result = list(
    "num.removed.codes" = length(full.codes) - length(reduced.codes),
    "score" = dist_corr_score,
    "codes.used" = list(reduced.codes)
  )

  attr(result, "codes.removed") <- as.character(reduced.codes);
  attr(result, "codes.maintained") <- as.character(full.codes %in% reduced.codes);
  if (keep_sets == TRUE) {
    attr(result, "ena.set") <- reduced.set;
  }
  # class(result) <- c("pria_result", class(result))
  return(result)
}

dist_corr <- function(A, B, sample_size = 1000000, method = "pearson") {
  m <- nrow(A)
  total_possible_pairs <- choose(m, 2)

  if (total_possible_pairs <= sample_size) {
    # CASE 1: Small m - Compute exact correlation
    #message("Computing exact correlation for all pairs...")
    d_A <- dist(A)
    d_B <- dist(B)
    return(cor(d_A, d_B, method = method))

  } else {
    # CASE 2: Large m - Compute sampled correlation
    #message(paste("Sampling", sample_size, "pairs for estimation..."))

    # Generate random indices
    idx1 <- sample.int(m, sample_size, replace = TRUE)
    idx2 <- sample.int(m, sample_size, replace = TRUE)

    # Filter out identity pairs (where idx1 == idx2)
    diff_mask <- idx1 != idx2
    idx1 <- idx1[diff_mask]
    idx2 <- idx2[diff_mask]

    # Memory-efficient distance calculation (Vectorized)
    # Distance = sqrt(sum of squared differences)
    dist_A <- sqrt(rowSums((A[idx1, , drop=FALSE] - A[idx2, , drop=FALSE])^2))
    dist_B <- sqrt(rowSums((B[idx1, , drop=FALSE] - B[idx2, , drop=FALSE])^2))

    return(cor(dist_A, dist_B, method = method))
  }
}
#' Conmpare
#'
#' @param A matrix, m by n1
#' @param B matrix, m by n2
#' @param n_dims number of first dimensions
#'
#' @return space similarity in terms of weighted correlation
#' @export
compute_weighted_average_correlation <- function(A, B, n_dims = NULL) {
  # A: Original Space (m x n1)
  # B: Reduced Space (m x n2)
  # k: Number of dimensions to compare (defaults to all available)

  # 1. Determine how many dimensions we can actually compare
  k = n_dims
  max_k <- min(ncol(A), ncol(B))
  if (is.null(k)) {
    k <- max_k
  } else {
    k <- min(k, max_k)
  }

  # 2. Extract the first k dimensions
  A_k <- A[, 1:k, drop = FALSE]
  B_k <- B[, 1:k, drop = FALSE]

  # 3. Calculate weights based on variance (Singular Values squared)
  # This ensures Dim 1 is more important than Dim 2
  weights <- apply(A_k, 2, var)

  # 4. Compute 1-to-1 correlation for each dimension
  r2_values <- numeric(k)
  for (j in 1:k) {
    # We use cor()^2 to get R-squared
    # use = "complete.obs" handles potential NAs
    r2_values[j] <- cor(A_k[, j], B_k[, j], use = "complete.obs")^2
  }

  # 5. Return the Weighted Average
  weighted_r2 <- sum(r2_values * weights) / sum(weights)

  return(list(
    similarity = weighted_r2,
    dimension_r2 = r2_values,
    weights = weights / sum(weights) # Normalized weights for reporting
  ))
}

#' @noRd
pria_backward_greedy <- function(x, remove.num, min_gof, correlation_fun,use_dist_corr=TRUE) {
  full_codes <- x$rotation$codes
  current_codes <- full_codes
  best_valid_set <- full_codes
  best_valid_score <- get_score_internal(full_codes, x, correlation_fun, min_gof,use_dist_corr=use_dist_corr)

  for (step in 1:remove.num) {
    if (length(current_codes) <= 3) break
    iteration_scores <- sapply(seq_along(current_codes), function(i) {
      get_score_internal(current_codes[-i], x, correlation_fun, min_gof,use_dist_corr=use_dist_corr)
    })

    if (max(iteration_scores) <= 0) break

    best_idx <- which.max(iteration_scores)
    current_codes <- current_codes[-best_idx]
    best_valid_set <- current_codes
    best_valid_score <- max(iteration_scores)
  }
  list(best_code_set = best_valid_set, final_performance = best_valid_score)
}

#' @noRd
pria_backward_beam<- function(x, remove.num, width, min_gof, correlation_fun,use_dist_corr=TRUE) {
  # 1. INITIALIZATION
  current_codes <- x$rotation$codes
  beams <- list(current_codes)

  # Start with the full set
  best_overall_set <- current_codes
  best_overall_score <- get_score_internal(current_codes, x, correlation_fun, min_gof,use_dist_corr=use_dist_corr)

  for (step in 1:remove.num) {
    # Generate candidates (using lapply for speed/efficiency)
    new_candidates <- list()
    for (b in seq_along(beams)) {
      curr <- beams[[b]]
      if(length(curr) <= 3) next
      # Add all possible 'minus-one' subsets
      new_candidates <- c(new_candidates, lapply(seq_along(curr), function(i) curr[-i]))
    }

    if (length(new_candidates) == 0) break

    # 2. UNIQUE & PRE-FILTERING (Improvement 1)
    # Sort and remove duplicates BEFORE scoring to save time
    candidate_strings <- vapply(new_candidates, function(s) paste(sort(s), collapse="|"), character(1))
    unique_idx <- which(!duplicated(candidate_strings))
    new_candidates <- new_candidates[unique_idx]

    # 3. EFFICIENT SCORING (Improvement 2)
    # Using vapply instead of growing a vector with c()
    candidate_scores <- vapply(new_candidates, function(s) {
      get_score_internal(s, x, correlation_fun, min_gof,use_dist_corr=use_dist_corr)
    }, numeric(1))

    # Identify valid candidates (score > 0)
    valid_idx <- which(candidate_scores > 0)

    # If no valid candidates exist, we cannot remove any more codes
    if (length(valid_idx) == 0) break

    # Filter to only valid ones
    new_candidates <- new_candidates[valid_idx]
    candidate_scores <- candidate_scores[valid_idx]

    # 4. PRUNING THE BEAM
    # Even if scores are decreasing, we pick the 'width' best remaining scores
    top_idx <- order(candidate_scores, decreasing = TRUE)[1:min(width, length(new_candidates))]
    beams <- new_candidates[top_idx]

    # 5. ALWAYS UPDATE (The "Greedy-Style" Fix)
    # Since you want to remove as many as possible, we update the best set
    # to be the best performer of THIS step, regardless of the previous step's score.
    best_overall_set <- beams[[1]]
    best_overall_score <- candidate_scores[top_idx[1]]
  }

  list(best_code_set = best_overall_set, final_performance = best_overall_score)
}

#' @noRd
pria_backward_pairwise <- function(x, remove.num, min_gof, correlation_fun,use_dist_corr=TRUE) {
  # 1. INITIALIZATION
  current_codes <- x$rotation$codes
  best_score <- get_score_internal(current_codes, x, correlation_fun, min_gof,use_dist_corr=use_dist_corr)

  # Calculate how many pairs we can remove
  # e.g., if remove.num is 5, floor(5/2) = 2 pairs (4 codes total)
  max_pair_steps <- floor(remove.num / 2)

  # 2. PAIRWISE REMOVAL PHASE
  if (max_pair_steps > 0) {
    for (s in 1:max_pair_steps) {
      # Safety: Need at least 5 codes to remove 2 and leave 3
      if (length(current_codes) <= 4) break

      combs <- combn(current_codes, 2)

      # Evaluate removing each pair
      iteration_scores <- apply(combs, 2, function(pair) {
        get_score_internal(setdiff(current_codes, pair), x, correlation_fun, min_gof,use_dist_corr=use_dist_corr)
      })

      # If no pairs can be removed safely, stop the pairwise phase
      if (all(iteration_scores <= 0)) break

      # Update the set with the best pair removal found
      best_idx <- which.max(iteration_scores)
      current_codes <- setdiff(current_codes, combs[, best_idx])
      best_score <- iteration_scores[best_idx]
    }
  }

  # 3. GREEDY CLEANUP PHASE
  # Calculate how many codes we have actually removed so far
  removed_so_far <- length(x$rotation$codes) - length(current_codes)

  # If we haven't reached remove.num (usually happens if remove.num is odd),
  # try to remove the remaining required codes one by one.
  while (removed_so_far < remove.num) {
    if (length(current_codes) <= 3) break

    # Single-item evaluation (Greedy)
    step_scores <- vapply(seq_along(current_codes), function(i) {
      get_score_internal(current_codes[-i], x, correlation_fun, min_gof)
    }, numeric(1))

    # If no further single removals are valid, stop completely
    if (all(step_scores <= 0)) break

    # Remove the best single item
    best_step_idx <- which.max(step_scores)
    current_codes <- current_codes[-best_step_idx]
    best_score <- step_scores[best_step_idx]

    # Update counter
    removed_so_far <- length(x$rotation$codes) - length(current_codes)
  }

  # 4. FINAL OUTPUT
  list(best_code_set = current_codes, final_performance = best_score)
}

#' @noRd
pria_brute_force <- function(x, remove.num, min_gof, correlation_fun, code_weights = NULL, lower_n = 12,use_dist_corr=TRUE) {
  # 1. INITIALIZATION & FILTERING
  full_codes <- x$rotation$codes
  total_n <- length(full_codes)

  # Logic for selecting candidate codes for removal
  if (!is.null(code_weights)) {
    # Ensure lower_n isn't larger than the total codes available
    actual_lower_n <- min(lower_n, total_n)

    # Sort codes by weight and pick the 'actual_lower_n' smallest
    # Assumes code_weights is a data.frame with columns 'code' and 'weight'
    sorted_weights <- code_weights[order(code_weights$weight), ]
    candidate_to_remove <- sorted_weights$code[1:actual_lower_n]

    # The 'fixed' codes that we definitely want to keep
    fixed_codes <- setdiff(full_codes, candidate_to_remove)

    # We need to remove 'remove.num' from the 'candidate_to_remove' pool
    # Therefore, we need to choose (actual_lower_n - remove.num) to KEEP from that pool
    pool_to_keep_count <- actual_lower_n - remove.num

    if (pool_to_keep_count < 0) {
      stop("remove.num cannot be greater than the number of candidate codes (lower_n).")
    }

    search_pool <- candidate_to_remove
    target_from_pool <- pool_to_keep_count

  } else {
    # Fallback to original behavior if no weights provided
    fixed_codes <- c()
    search_pool <- full_codes
    target_from_pool <- total_n - remove.num
  }

  # 2. COMBINATORIAL SAFETY CHECK
  num_combs <- choose(length(search_pool), target_from_pool)

  if (num_combs > 1e6) {
    stop(paste("Too many combinations:", round(num_combs),
               ". Even with filtering, brute force is too slow."))
  }

  # 3. GENERATE SUBSETS FROM THE CANDIDATE POOL
  pool_combinations <- combn(search_pool, target_from_pool, simplify = FALSE)

  # 4. EVALUATE ALL SETS (Combining Fixed + Combinations)
  scores <- vapply(pool_combinations, function(pool_set) {
    # Merge the guaranteed codes with this specific combination
    current_set <- c(fixed_codes, pool_set)
    get_score_internal(current_set, x, correlation_fun, min_gof,use_dist_corr=use_dist_corr)
  }, numeric(1))

  # 5. FIND THE OPTIMUM
  best_idx <- which.max(scores)

  if (scores[best_idx] <= 0) {
    warning("No valid code sets (score > 0) were found.")
  }

  list(
    best_code_set = c(fixed_codes, pool_combinations[[best_idx]]),
    final_performance = scores[best_idx],
    combinations_tested = num_combs
  )
}

#' Remove Codes
#'
#' @param accum TBD
#' @param code_name_to_remove TBD
#'
#' @return TBD
#' @export
remove.codes.from.accum <- function(accum, code_name_to_remove){
  # browser(expr = { length(code_name_to_remove) > 1})
  # adj_column_has_code <- code_name_to_remove == accum$rotation$adjacency.key
  # adj_column_has_code_as_vector <- apply(adj_column_has_code, 2, any)
  adj <- accum$rotation$adjacency.key

  if (is.matrix(adj)) {
    cols_to_keep <- !apply(adj, 2, function(col) any(col %in% code_name_to_remove))
    accum$rotation$adjacency.key <- adj[, cols_to_keep, drop = FALSE]
  } else {
    adj_column_has_code_as_vector <- colSums(sapply(adj, `%in%`, code_name_to_remove)) > 0
    accum$rotation$adjacency.key <- accum$rotation$adjacency.key[, !..adj_column_has_code_as_vector]
  }

  accum$rotation$codes <- accum$rotation$codes[!accum$rotation$codes %in% code_name_to_remove]
  colnames_to_remove <- colnames(as.matrix(accum$connection.counts))[if(is.matrix(adj)) !cols_to_keep else adj_column_has_code_as_vector]
  accum$connection.counts <- accum$connection.counts[, !colnames(accum$connection.counts) %in% colnames_to_remove, with = F]
  return(accum)
}

#' PRIA: Proactive Reduced Interaction Analysis
#'
#' @description
#' Optimizes ENA models by identifying the best subset of codes to maintain
#' model performance while reducing complexity.
#'
#' @param x An \code{ena.set} object.
#' @param remove.num number of codes to remove.
#' @param threshold threshold for code selection, default = 0.95.
#' @param rebuild optional function taking a reduced accumulation and returning
#'   the corresponding ENA set. Supply this when \code{x} was not built by
#'   \code{rENA::ena.make.set} with rotation recorded on \code{_function.params}
#'   -- notably ordered/directed models, which the default cannot rebuild.
#'
#' @export
pria <- function(x,
                 remove.num=3,
                 threshold = 0.95,
                 rebuild = NULL) {
  codes = x$rotation$codes
  m = length(codes)
  remove.num = min(remove.num,m-3)
  # Guard remove.num, not m: a model at the 3-code floor gives remove.num = 0,
  # and `for(k in 1:0)` counts *down* through k = 1, 0, so combn(codes, 0) is
  # reached and fails with "length of 'dimnames' [2] not equal to array extent".
  if(remove.num < 1)
  {
    return(NULL)
  }
  vr = 0
  n = 0
  result = NULL
  for(k in 1:remove.num)
  {
    combos <- combn(codes,k)
    for(i in 1:ncol(combos)){
      current_removed <- combos[,i]
      score <- get_pria_scores_2Ds(x,
                          removed_codes = current_removed,
                          threshold = threshold,
                          rebuild = rebuild )
      vr1 = score$vr1
      if(!is.null(vr1))
      {
        if(k>n)
        {
          n = k
          vr = vr1
          result = score
        }else
        {
          if(vr1>vr)
          {
            vr = vr1
            result = score
          }
        }
      }
    }
  }
  return(result$removed_codes)
}


#' backup code
pria_backup2 <- function(x, remove.num,
                 method = c("greedy", "beam", "pairwise", "brute"),
                 width = 3,
                 min_gof = 0.5,
                 correlation_fun = correlations,
                 use_dist_corr = TRUE,
                 code_weight=NULL,
                 lower_n=12) {

  method <- match.arg(method)

  res <- switch(method,
                "greedy"   = pria_backward_greedy(x, remove.num, min_gof, correlation_fun,use_dist_corr=use_dist_corr),
                "beam"     = pria_backward_beam(x, remove.num, width, min_gof, correlation_fun,use_dist_corr=use_dist_corr),
                "pairwise" = pria_backward_pairwise(x, remove.num, min_gof, correlation_fun,use_dist_corr=use_dist_corr),
                "brute"    = pria_brute_force(x, remove.num, min_gof, correlation_fun,code_weight=NULL,lower_n=12,use_dist_corr=use_dist_corr)
  )

  # Attach common metadata
  res$method <- method
  res$removed_codes <- setdiff(x$rotation$codes, res$best_code_set)

  return(res)
}
#' PRIA
#'
#' @param x TBD
#' @param remove.num TBD
#' @param correlation_fun TBD
#' @param keep_sets TBD
#'
#' @return TBD
#' @export
pria_backup <- function(x, remove.num, correlation_fun = correlations, keep_sets = FALSE) {
  code.left.all <- list()
  full.code.num <- length(x$rotation$codes)
  full.codes <- x$rotation$codes

  if(remove.num > full.code.num - 3){
    stop("Try removing less codes.") ### clear comment
  }
  else{
    for (remove.num.each.iter in 1:remove.num) {
      code.combs <- combn(x$rotation$codes, remove.num.each.iter)
      for (code.comb.num in 1:ncol(code.combs)) {
        code.comb <- code.combs[,code.comb.num]
        num.code.comb <- length(code.comb)
        reduced.accum <- remove.codes.from.accum(x, code.comb)
        reduced.set <- rENA::ena.make.set(
          reduced.accum,
          rotation.by = x$`_function.params`$rotation.by,
          rotation.params = x$`_function.params`$rotation.params
        )

        # browser(expr = { length(code.left.all) + 1 == 367 })
        result <- compare.two.models(x, reduced.set, FUN = correlation_fun, keep_sets = keep_sets)
        code.left.all[[length(code.left.all) + 1]] <- result
      }
    }
  }

  attrs <- lapply(code.left.all, attributes)
  code.left.all <- data.table::rbindlist(code.left.all, idcol =  TRUE);
  # data.table::setorderv(code.left.all, "goodness.of.fit", order = -1L)
  class(code.left.all) <- c("pria_table", "data.frame"); #class(code.left.all));
  data.table::setattr(code.left.all, "codes.maintained", lapply(attrs, `[[`, "codes.maintained"))
  data.table::setattr(code.left.all, "codes.removed", lapply(attrs, `[[`, "codes.removed"))
  ena_sets <- lapply(attrs, `[[`, "ena.set");
  class(ena_sets) <- c("ena.sets", class(ena_sets));
  data.table::setattr(code.left.all, "ena.set", ena_sets)

  code.left.all
}

1		#' Remove Codes
2		#'
3		#' @param accum TBD
4		#' @param code_name_to_remove TBD
5		#'
6		#' @return TBD
7		#' @export
8		remove.codes.from.accum <- function(accum, code_name_to_remove){
9		# browser(expr = { length(code_name_to_remove) > 1})
10		# adj_column_has_code <- code_name_to_remove == accum$rotation$adjacency.key
11		# adj_column_has_code_as_vector <- apply(adj_column_has_code, 2, any)
12	69x	adj <- accum$rotation$adjacency.key
13
14	69x	if (is.matrix(adj)) {
15	!	cols_to_keep <- !apply(adj, 2, function(col) any(col %in% code_name_to_remove))
16	!	accum$rotation$adjacency.key <- adj[, cols_to_keep, drop = FALSE]
17		} else {
18	69x	adj_column_has_code_as_vector <- colSums(sapply(adj, `%in%`, code_name_to_remove)) > 0
19	69x	accum$rotation$adjacency.key <- accum$rotation$adjacency.key[, !..adj_column_has_code_as_vector]
20		}
21
22	69x	accum$rotation$codes <- accum$rotation$codes[!accum$rotation$codes %in% code_name_to_remove]
23	69x	colnames_to_remove <- colnames(as.matrix(accum$connection.counts))[if(is.matrix(adj)) !cols_to_keep else adj_column_has_code_as_vector]
24	69x	accum$connection.counts <- accum$connection.counts[, !colnames(accum$connection.counts) %in% colnames_to_remove, with = F]
25	69x	return(accum)
26		}
27
28		#' PRIA: Proactive Reduced Interaction Analysis
29		#'
30		#' @description
31		#' Optimizes ENA models by identifying the best subset of codes to maintain
32		#' model performance while reducing complexity.
33		#'
34		#' @param x An \code{ena.set} object.
35		#' @param remove.num number of codes to remove.
36		#' @param threshold threshold for code selection, default = 0.95.
37		#' @param rebuild optional function taking a reduced accumulation and returning
38		#' the corresponding ENA set. Supply this when \code{x} was not built by
39		#' \code{rENA::ena.make.set} with rotation recorded on \code{_function.params}
40		#' -- notably ordered/directed models, which the default cannot rebuild.
41		#'
42		#' @export
43		pria <- function(x,
44		remove.num=3,
45		threshold = 0.95,
46		rebuild = NULL) {
47	2x	codes = x$rotation$codes
48	2x	m = length(codes)
49	2x	remove.num = min(remove.num,m-3)
50		# Guard remove.num, not m: a model at the 3-code floor gives remove.num = 0,
51		# and `for(k in 1:0)` counts down through k = 1, 0, so combn(codes, 0) is
52		# reached and fails with "length of 'dimnames' [2] not equal to array extent".
53	2x	if(remove.num < 1)
54		{
55	!	return(NULL)
56		}
57	2x	vr = 0
58	2x	n = 0
59	2x	result = NULL
60	2x	for(k in 1:remove.num)
61		{
62	3x	combos <- combn(codes,k)
63	3x	for(i in 1:ncol(combos)){
64	44x	current_removed <- combos[,i]
65	44x	score <- get_pria_scores_2Ds(x,
66	44x	removed_codes = current_removed,
67	44x	threshold = threshold,
68	44x	rebuild = rebuild )
69	44x	vr1 = score$vr1
70	44x	if(!is.null(vr1))
71		{
72	5x	if(k>n)
73		{
74	3x	n = k
75	3x	vr = vr1
76	3x	result = score
77		}else
78		{
79	2x	if(vr1>vr)
80		{
81	!	vr = vr1
82	!	result = score
83		}
84		}
85		}
86		}
87		}
88	2x	return(result$removed_codes)
89		}
90
91
92		#' backup code
93		pria_backup2 <- function(x, remove.num,
94		method = c("greedy", "beam", "pairwise", "brute"),
95		width = 3,
96		min_gof = 0.5,
97		correlation_fun = correlations,
98		use_dist_corr = TRUE,
99		code_weight=NULL,
100		lower_n=12) {
101
102	!	method <- match.arg(method)
103
104	!	res <- switch(method,
105	!	"greedy" = pria_backward_greedy(x, remove.num, min_gof, correlation_fun,use_dist_corr=use_dist_corr),
106	!	"beam" = pria_backward_beam(x, remove.num, width, min_gof, correlation_fun,use_dist_corr=use_dist_corr),
107	!	"pairwise" = pria_backward_pairwise(x, remove.num, min_gof, correlation_fun,use_dist_corr=use_dist_corr),
108	!	"brute" = pria_brute_force(x, remove.num, min_gof, correlation_fun,code_weight=NULL,lower_n=12,use_dist_corr=use_dist_corr)
109		)
110
111		# Attach common metadata
112	!	res$method <- method
113	!	res$removed_codes <- setdiff(x$rotation$codes, res$best_code_set)
114
115	!	return(res)
116		}
117		#' PRIA
118		#'
119		#' @param x TBD
120		#' @param remove.num TBD
121		#' @param correlation_fun TBD
122		#' @param keep_sets TBD
123		#'
124		#' @return TBD
125		#' @export
126		pria_backup <- function(x, remove.num, correlation_fun = correlations, keep_sets = FALSE) {
127	3x	code.left.all <- list()
128	3x	full.code.num <- length(x$rotation$codes)
129	3x	full.codes <- x$rotation$codes
130
131	3x	if(remove.num > full.code.num - 3){
132	!	stop("Try removing less codes.") ### clear comment
133		}
134		else{
135	3x	for (remove.num.each.iter in 1:remove.num) {
136	3x	code.combs <- combn(x$rotation$codes, remove.num.each.iter)
137	3x	for (code.comb.num in 1:ncol(code.combs)) {
138	24x	code.comb <- code.combs[,code.comb.num]
139	24x	num.code.comb <- length(code.comb)
140	24x	reduced.accum <- remove.codes.from.accum(x, code.comb)
141	24x	reduced.set <- rENA::ena.make.set(
142	24x	reduced.accum,
143	24x	rotation.by = x$`_function.params`$rotation.by,
144	24x	rotation.params = x$`_function.params`$rotation.params
145		)
146
147		# browser(expr = { length(code.left.all) + 1 == 367 })
148	24x	result <- compare.two.models(x, reduced.set, FUN = correlation_fun, keep_sets = keep_sets)
149	24x	code.left.all[[length(code.left.all) + 1]] <- result
150		}
151		}
152		}
153
154	3x	attrs <- lapply(code.left.all, attributes)
155	3x	code.left.all <- data.table::rbindlist(code.left.all, idcol = TRUE);
156		# data.table::setorderv(code.left.all, "goodness.of.fit", order = -1L)
157	3x	class(code.left.all) <- c("pria_table", "data.frame"); #class(code.left.all));
158	3x	data.table::setattr(code.left.all, "codes.maintained", lapply(attrs, `[[`, "codes.maintained"))
159	3x	data.table::setattr(code.left.all, "codes.removed", lapply(attrs, `[[`, "codes.removed"))
160	3x	ena_sets <- lapply(attrs, `[[`, "ena.set");
161	3x	class(ena_sets) <- c("ena.sets", class(ena_sets));
162	3x	data.table::setattr(code.left.all, "ena.set", ena_sets)
163
164	3x	code.left.all
165		}

1		#' Title
2		#'
3		#' @param x TBD
4		#' @param ... TBD
5		#'
6		#' @return TBD
7		#' @export
8		"print.pria_table" <- function(x, ...) {
9	3x	print(data.table::as.data.table(x), ...);
10		}
11
12		#' Title
13		#'
14		#' @param x TBD
15		#' @param pattern TBD
16		#'
17		#' @return TBD
18		#' @export
19		".DollarNames.pria_table" <- function(x, pattern = "") {
20	2x	columns <- colnames(x);
21	2x	choose_from <- c(columns, c("codes.maintained", "codes.removed", "ena.set"));
22	2x	grep(pattern = pattern, x = choose_from, value = TRUE)
23		}
24
25		#' Title
26		#'
27		#' @param x TBD
28		#' @param i TBD
29		#' @param j TBD
30		#' @param ... TBD
31		#' @param drop TBD
32		#'
33		#' @return TBD
34		#' @export
35		"[.pria_table" <- function(x, i, j, ..., drop = TRUE) {
36	26x	x_ <- data.table::copy(x);
37	26x	class(x_) <- class(x_)[!class(x_) %chin% "pria_table"]
38
39	26x	pria_row <- x_[i, ]
40
41	26x	class(pria_row) <- c("pria_table", class(pria_row));
42
43	26x	pria_row;
44		}
45
46		#' Override element extraction
47		#'
48		#' @param x TBD
49		#' @param i TBD
50		#' @param j TBD
51		#'
52		#' @return TBD
53		#' @export
54		"$.pria_table" <- function(x, i, j) {
55	41x	if(i %in% colnames(x)) {
56	27x	x[[i]]
57		}
58		else {
59	14x	res <- attr(x = x, which = i, exact = TRUE)
60	14x	res
61		}
62		}
63
64		#' Title
65		#'
66		#' @param e1 TBD
67		#' @param e2 TBD
68		#'
69		#' @return TBD
70		#' @export
71		">.pria_table" <- function(e1, e2) {
72	3x	reduced <- e1[
73	3x	e1$goodness.of.fit.ci_lower > e2 &
74	3x	e1$points.ci_lower > e2 &
75	3x	e1$nodes.ci_lower > e2,
76		]
77
78	3x	reduced <- reduced[order(-reduced$num.removed.codes, -reduced$goodness.of.fit),]
79	3x	class(reduced) <- c("pria_table", class(reduced));
80
81	3x	reduced
82		}
83
84		#' Title
85		#'
86		#' @param ... TBD
87		#' @param na.rm TBD
88		#'
89		#' @return TBD
90		#' @export
91		"max.pria_table" <- function(..., na.rm = FALSE) {
92	4x	all <- list(...)[[1]];
93	4x	best_index <- order(-all$num.removed.codes, -all$goodness.of.fit)[1]
94	4x	best <- all[best_index,][1,]
95
96	4x	best
97		}
98
99		#' Title
100		#'
101		#' @param x TBD
102		#' @param i TBD
103		#' @param ... TBD
104		#'
105		#' @return TBD
106		#' @export
107		"[.ena.sets" <- function(x, i, ...) {
108	3x	x_ <- x;
109	3x	class(x_) <- class(x_)[!class(x_) %chin% "ena.sets"]
110
111	3x	res <- `[`(x_, i, ...);
112	3x	class(res) <- c("ena.sets", class(res))
113	3x	res
114		}
115
116		#' Title
117		#'
118		#' @param x TBD
119		#' @param pattern TBD
120		#'
121		#' @return TBD
122		#' @export
123		".DollarNames.ena.sets" <- function(x, pattern = "") {
124	2x	grep(pattern = pattern, x = names(x[[1]]), value = TRUE)
125		}
126
127		#' Override element extraction
128		#'
129		#' @param x TBD
130		#' @param i TBD
131		#' @param j TBD
132		#'
133		#' @return TBD
134		#' @export
135		"$.ena.sets" <- function(x, i, j) {
136	5x	res <- lapply(x, `[[`, i)
137	5x	class(res) <- c("ena.sets", class(res))
138	5x	res
139		}
140
141		#' Title
142		#'
143		#' @param e1 TBD
144		#' @param e2 TBD
145		#'
146		#' @return TBD
147		#' @export
148		"==.ena.sets" <- function(e1,e2) {
149	1x	lapply(e1, `==`, e2)
150		}
151
152		#' Internal Performance Calculation
153		#' @noRd
154		calc_perf_internal <- function(result, min_gof) {
155	!	gof <- result$goodness.of.fit
156	!	if (gof < min_gof) gof <- 0
157		#return(gof * max(0, result$points) * max(0, result$nodes))
158	!	return(gof * result$points * result$nodes)
159		}
160
161		#' Internal Model Scorer
162		#' @noRd
163		get_score_internal <- function(current_codes, x, correlation_fun, min_gof,use_dist_corr=TRUE) {
164	!	full_codes <- x$rotation$codes
165	!	removed <- setdiff(full_codes, current_codes)
166
167	!	if (length(removed) == 0) {
168	!	reduced.set <- x
169		} else {
170	!	reduced.accum <- remove.codes.from.accum(x, removed)
171	!	if (is.matrix(reduced.accum$rotation$adjacency.key)) {
172	!	wh <- if (is.null(x$`_function.params`$rotation.by)) rENA::ena.svd else x$`_function.params`$rotation.by
173	!	reduced.set <- reduced.accum \|>
174	!	rENA::sphere_norm() \|>
175	!	rENA::center() \|>
176	!	rENA::rotate(wh = wh) \|>
177	!	rENA::project() \|>
178	!	rENA::optimize()
179	!	reduced.set$rotation$adjacency.key <- reduced.accum$rotation$adjacency.key
180		} else {
181	!	reduced.set <- rENA::ena.make.set(
182	!	reduced.accum,
183	!	rotation.by = x$`_function.params`$rotation.by,
184	!	rotation.params = x$`_function.params`$rotation.params
185		)
186		}
187		}
188
189	!	if (use_dist_corr) {
190	!	res <- compare.two.models.by.dist(x, reduced.set, keep_sets = FALSE)
191	!	return(res$score)
192		} else {
193	!	res <- compare.two.models(x, reduced.set, FUN = correlation_fun, keep_sets = FALSE)
194		}
195	!	return(calc_perf_internal(res, min_gof))
196		}
197		#' Compute all pria scores
198		#'
199		#' @param current_codes current codes with some codes removed
200		#' @param x ENA set with full codes
201		#' @param correlation_fun correlation function
202		#' @param min_gof GOF threshold
203		#' @param grouping_col a grouping columns for effect size computation, default=NULL
204		#' @param include_reduced_set whether or not include the reduced set in the return,default=FALSE
205		#' @param n_dim number of first dimensions to use in similarity computations, default=2
206		#'
207		#' @return a list of scores, and optionally the reduced ENA set
208		#' @export
209		get_pria_scores <- function(current_codes, x, correlation_fun, min_gof,grouping_col=NULL,include_reduced_set=FALSE,n_dim=2) {
210	!	full_codes <- x$rotation$codes
211	!	removed <- setdiff(full_codes, current_codes)
212	!	reduced.set =NULL
213	!	if (length(removed) == 0) {
214	!	reduced.set <- x
215		} else {
216	!	reduced.accum <- remove.codes.from.accum(x, removed)
217	!	reduced.set <- rENA::ena.make.set(
218	!	reduced.accum,
219	!	rotation.by = x$`_function.params`$rotation.by,
220	!	rotation.params = x$`_function.params`$rotation.params
221		)
222		}
223
224	!	res1<-compare.two.models.by.dist(x,reduced.set,keep_sets=FALSE,n_dim)
225	!	space_sim=res1$score
226
227	!	res2 <- compare.two.models(x, reduced.set, FUN = correlation_fun, keep_sets = FALSE)
228	!	gof <-res2$goodness.of.fit
229	!	point_corr<-res2$points
230	!	node_corr<-res2$nodes
231	!	eta2 = NULL
232	!	if(!is.null(grouping_col))
233		{
234	!	first_dim <- as.matrix(reduced.set$points)[,1]
235	!	cond_vector <- reduced.set$meta.data[[grouping_col]]
236
237		# Run ANOVA
238	!	fit <- aov(first_dim ~ cond_vector)
239	!	sum_stats <- summary(fit)[[1]]
240
241	!	ss_between <- sum_stats[1, "Sum Sq"]
242	!	ss_residual <- sum_stats[2, "Sum Sq"]
243	!	eta2 <- ss_between / (ss_between + ss_residual)
244		}
245	!	if(!include_reduced_set)
246		{
247	!	reduced.set=NULL
248		}
249	!	return(list(space_sim = space_sim,
250	!	gof = gof,
251	!	point_corr = point_corr,
252	!	node_corr=node_corr,
253	!	eta2=eta2,
254	!	reduced_set=reduced.set))
255		}
256		#' Compute pria scores on first two dimensions
257		#'
258		#' @param x ENA set with full codes
259		#' @param removed_codes codes removed
260		#' @param threshold threshold, default=0.95
261		#' @param rebuild optional function taking a reduced accumulation and returning
262		#' the corresponding ENA set; see \code{\link{pria}}.
263		#' @return variance explained on dim 1 of the reduced set
264		#' @export
265		get_pria_scores_2Ds <- function(x,removed_codes,threshold=0.95,rebuild=NULL) {
266		# construct reduced set
267	44x	reduced.set =NULL
268	44x	if (length(removed_codes) == 0) {
269	!	reduced.set <- x
270		} else {
271	44x	reduced.accum <- remove.codes.from.accum(x, removed_codes)
272	44x	reduced.set <- if (is.function(rebuild)) {
273		# Callers that know how the original set was built supply `rebuild`. The
274		# default below cannot serve every model: ena.make.set only builds
275		# undirected sets, so an ordered/directed accumulation fails outright,
276		# and it reads the rotation off `_function.params`, which callers such
277		# as the ENA webtool never populate -- silently scoring each candidate
278		# against a default-SVD rotation instead of the model's own.
279	!	rebuild(reduced.accum)
280		} else {
281	44x	rENA::ena.make.set(
282	44x	reduced.accum,
283	44x	rotation.by = x$`_function.params`$rotation.by,
284	44x	rotation.params = x$`_function.params`$rotation.params
285		)
286		}
287		}
288		# compute relative variance
289	44x	vr1 = reduced.set$model$variance[1]
290
291		# compute gof, original set
292	44x	pts0 = x$points
293	44x	cts0 = x$model$centroids
294	44x	pts20 = as.matrix(pts0)[,1:2,drop=FALSE]
295	44x	cts20 = as.matrix(cts0)[,1:2,drop=FALSE]
296	44x	cor_result0 <-rENA::ena_correlation(pts20,cts20)
297	44x	gof_1_o = cor_result0[1,1]
298	44x	gof_2_o = cor_result0[2,1]
299		# compute gof, reduced set
300	44x	pts = reduced.set$points
301	44x	cts = reduced.set$model$centroids
302	44x	pts2 = as.matrix(pts)[,1:2,drop=FALSE]
303	44x	cts2 = as.matrix(cts)[,1:2,drop=FALSE]
304	44x	cor_result1 <-rENA::ena_correlation(pts2,cts2)
305	44x	gof_1 = cor_result1[1,1]
306	44x	gof_2 = cor_result1[2,1]
307
308	44x	if(min(gof_1/gof_1_o,gof_2/gof_2_o)<threshold)
309		{
310	5x	return(NULL)
311		}
312
313		# compute point, node correlations
314	39x	full_points = as.matrix(x$points)[,1:2,drop=FALSE]
315	39x	reduce_points = pts2
316	39x	cor_result2 = rENA::ena_correlation(full_points,reduce_points)
317
318	39x	full_nodes = x$rotation$nodes
319	39x	reduced_nodes =reduced.set$rotation$nodes
320		# Keep the full model's rows for the codes the reduced model retained. Align
321		# on rotation$codes rather than a `code` column: node tables from some
322		# pipelines (the ENA webtool among them) carry only dimension columns, and
323		# `full_nodes$code` is then NULL, selecting zero rows and leaving the
324		# correlation below to fail with "subscript out of bounds".
325	39x	keep_nodes <- if (!is.null(full_nodes$code) && !is.null(reduced_nodes$code)) {
326	39x	full_nodes$code %in% reduced_nodes$code
327		} else {
328	!	x$rotation$codes %in% reduced.set$rotation$codes
329		}
330	39x	full_nodes <- full_nodes[keep_nodes, ]
331	39x	full_nodes<-as.matrix(full_nodes)[,1:2,drop=FALSE]
332	39x	reduced_nodes<-as.matrix(reduced_nodes)[,1:2,drop=FALSE]
333	39x	cor_result3 = rENA::ena_correlation(full_nodes,reduced_nodes)
334
335		# check flipping
336	39x	point_corr_1 = cor_result2[1,1]
337	39x	point_corr_2 = cor_result2[2,1]
338	39x	node_corr_1 = cor_result3[1,1]
339	39x	node_corr_2 = cor_result3[2,1]
340	39x	if(point_corr_1<0)
341		{
342	10x	point_corr_1 = -point_corr_1
343	10x	node_corr_1 = - node_corr_1
344		}
345	39x	if(point_corr_2<0)
346		{
347	28x	point_corr_2 = -point_corr_2
348	28x	node_corr_2 = - node_corr_2
349		}
350
351	39x	if(min(point_corr_1,point_corr_2,node_corr_1,node_corr_2)<threshold)
352		{
353	34x	return(NULL)
354		}
355
356	5x	return(list(
357	5x	removed_codes = removed_codes,
358	5x	vr1=vr1,
359	5x	gof1 = gof_1/gof_1_o,
360	5x	gof2 = gof_2/gof_2_o,
361	5x	point1 = point_corr_1,
362	5x	point2 = point_corr_2,
363	5x	node1 = node_corr_1,
364	5x	node2 = node_corr_2))
365		}
366		#' Compute pria scores on first two dimensions
367		#'
368		#' @param current_codes current codes with some codes removed
369		#' @param x ENA set with full codes
370		#' @param correlation_fun correlation function
371		#' @param min_gof GOF threshold
372		#' @param grouping_col a grouping columns for effect size computation, default=NULL
373		#' @param include_reduced_set whether or not include the reduced set in the return,default=FALSE
374		#' @param n_dim number of first dimensions to use in similarity computations, default=2
375		#'
376		#' @return a list of scores, and optionally the reduced ENA set
377		#' @export
378		get_pria_scores_2D <- function(x,removed_codes ) {
379		# construct reduced set
380	!	reduced.set =NULL
381	!	if (length(removed_codes) == 0) {
382	!	reduced.set <- x
383		} else {
384	!	reduced.accum <- remove.codes.from.accum(x, removed_codes)
385	!	reduced.set <- rENA::ena.make.set(
386	!	reduced.accum,
387	!	rotation.by = x$`_function.params`$rotation.by,
388	!	rotation.params = x$`_function.params`$rotation.params
389		)
390		}
391		# compute pair-wise distance similarity in 2D
392	!	res1<-compare.two.models.by.dist(x,reduced.set,keep_sets=FALSE,2)
393
394		#***
395	!	space_sim=res1$score
396		#***
397		#
398		# compute gof, original set
399	!	pts0 = x$points
400	!	cts0 = x$model$centroids
401	!	pts20 = as.matrix(pts0)[,1:2,drop=FALSE]
402	!	cts20 = as.matrix(cts0)[,1:2,drop=FALSE]
403	!	cor_result0 <-rENA::ena_correlation(pts20,cts20)
404		# compute gof, reduced set
405	!	pts = reduced.set$points
406	!	cts = reduced.set$model$centroids
407	!	pts2 = as.matrix(pts)[,1:2,drop=FALSE]
408	!	cts2 = as.matrix(cts)[,1:2,drop=FALSE]
409	!	cor_result1 <-rENA::ena_correlation(pts2,cts2)
410
411		# compute point, node correlations
412	!	full_points = as.matrix(x$points)[,1:2,drop=FALSE]
413	!	reduce_points = pts2
414	!	cor_result2 = rENA::ena_correlation(full_points,reduce_points)
415
416	!	full_nodes = x$rotation$nodes
417	!	reduced_nodes =reduced.set$rotation$nodes
418		# Keep the full model's rows for the codes the reduced model retained. Align
419		# on rotation$codes rather than a `code` column: node tables from some
420		# pipelines (the ENA webtool among them) carry only dimension columns, and
421		# `full_nodes$code` is then NULL, selecting zero rows and leaving the
422		# correlation below to fail with "subscript out of bounds".
423	!	keep_nodes <- if (!is.null(full_nodes$code) && !is.null(reduced_nodes$code)) {
424	!	full_nodes$code %in% reduced_nodes$code
425		} else {
426	!	x$rotation$codes %in% reduced.set$rotation$codes
427		}
428	!	full_nodes <- full_nodes[keep_nodes, ]
429	!	full_nodes<-as.matrix(full_nodes)[,1:2,drop=FALSE]
430	!	reduced_nodes<-as.matrix(reduced_nodes)[,1:2,drop=FALSE]
431	!	cor_result3 = rENA::ena_correlation(full_nodes,reduced_nodes)
432
433
434		# check flipping
435	!	point_corr_1 = cor_result2[1,1]
436	!	point_corr_2 = cor_result2[2,1]
437	!	node_corr_1 = cor_result3[1,1]
438	!	node_corr_2 = cor_result3[2,1]
439	!	if(point_corr_1<0)
440		{
441	!	point_corr_1 = -point_corr_1
442	!	node_corr_1 = - node_corr_1
443		}
444	!	if(point_corr_2<0)
445		{
446	!	point_corr_2 = -point_corr_2
447	!	node_corr_2 = - node_corr_2
448		}
449
450		# compute cosine between two rotation hyperplanes
451	!	full_rm = x$rotation.matrix
452	!	reduced_rm = reduced.set$rotation.matrix
453		# base-R subset (dplyr isn't imported, so `filter` fell through to stats::filter)
454	!	full_rm <- full_rm[full_rm$codes %in% reduced_rm$codes, ]
455	!	full_rm = as.matrix(full_rm)[,2:3,drop=FALSE]
456	!	reduced_rm = as.matrix(reduced_rm)[,2:3,drop=FALSE]
457	!	aa = t(full_rm)%*%reduced_rm
458	!	cos_rotation = svd(aa)$d[2]
459		# compute consine between 2D points hyperplanes
460
461	!	A <- full_points
462	!	B <- reduce_points
463
464		# 2. Get orthonormal bases using QR decomposition
465	!	Q1 <- qr.Q(qr(A))
466	!	Q2 <- qr.Q(qr(B))
467
468		# 3. Compute the inner product matrix
469	!	M <- t(Q1) %*% Q2
470
471		# 4. Get singular values (cosines of principal angles)
472	!	cos_point <- svd(M)$d[2]
473
474
475	!	return(list(space_sim = space_sim,
476	!	gof_1_o = cor_result0[1,1],
477	!	gof_2_o = cor_result0[2,1],
478	!	gof_1 = cor_result1[1,1],
479	!	gof_2 = cor_result1[2,1],
480	!	point_corr_1 = point_corr_1,
481	!	point_corr_2 = point_corr_2,
482	!	node_corr_1 = node_corr_1 ,
483	!	node_corr_2 = node_corr_2,
484	!	cos_rotation = cos_rotation,
485	!	cos_point = cos_point
486		))
487		}

1		#' Correlations
2		#'
3		#' @param enaset TBD
4		#' @param pts TBD
5		#' @param cts TBD
6		#' @param dims TBD
7		#'
8		#' @return TBD
9		#' @export
10		correlations_slow <- function(enaset, pts = NULL, cts = NULL, dims = c(1:2)) {
11		# I rewrote the goodness of fit test since the current package doesn't return CI
12	9x	if(is.null(pts) \|\| is.null(cts)) {
13	9x	pts <- enaset$points
14	9x	cts <- enaset$model$centroids
15		}
16
17	9x	pComb = combn(nrow(pts), 2)
18		# pComb = rENA::combn_c2(nrow(pts)) + 1
19	9x	point1 = pComb[1,]
20	9x	point2 = pComb[2,]
21
22	9x	points = as.matrix(pts)
23	9x	centroids = as.matrix(cts)
24
25		# svdDiff = matrix(points[point1, dims] - points[point2, dims], ncol=length(dims), nrow=length(point1))
26	9x	svdDiff = points[point1, dims, drop = FALSE] - points[point2, dims, drop = FALSE];
27		# optDiff = matrix(centroids[point1, dims] - centroids[point2, dims], ncol=length(dims), nrow=length(point1))
28	9x	optDiff = centroids[point1, dims, drop = FALSE] - centroids[point2, dims, drop = FALSE]
29
30	9x	cor_result <- cor.test(as.numeric(svdDiff[,1]), as.numeric(optDiff[,1]))
31	9x	corr_lower_bound <- cor_result$conf.int[1]
32	9x	corr_upper_bound <- cor_result$conf.int[2]
33	9x	corr_estimates <- cor_result$estimate[1]
34
35	9x	corr_results <- list(
36	9x	"corr_lower_bound"=corr_lower_bound,
37	9x	"corr_upper_bound"=corr_upper_bound,
38	9x	"corr_estimates"=as.numeric(corr_estimates)
39		)
40	9x	return(corr_results);
41		}
42
43		#' Correlations
44		#'
45		#' @param enaset TBD
46		#' @param pts TBD
47		#' @param cts TBD
48		#' @param dims TBD
49		#' @param conf.level TBD
50		#'
51		#' @return TBD
52		#' @export
53		correlations <- function(enaset, pts = NULL, cts = NULL, dims = c(1), conf.level = 0.95) {
54		# I rewrote the goodness of fit test since the current package doesn't return CI
55	16x	if(is.null(pts) \|\| is.null(cts)) {
56	16x	pts <- enaset$points
57	16x	cts <- enaset$model$centroids
58		}
59
60	16x	points = as.matrix(pts)[, dims, drop = FALSE];
61	16x	centroids = as.matrix(cts)[, dims, drop = FALSE];
62
63		# pc <- matrix(c(points[,1], centroids[,1]), ncol = 2)
64
65	16x	cor_result <- rENA::ena_correlation(points, centroids)
66
67	16x	corr_results <- list(
68	16x	"corr_lower_bound" = abs(cor_result[1,2]),
69	16x	"corr_upper_bound" = abs(cor_result[1,3]),
70	16x	"corr_estimates" = abs(cor_result[1,1])
71		)
72
73	16x	return(corr_results);
74		}
75
76		#' Compute GOF and cors
77		#'
78		#' @param full.set TBD
79		#' @param reduced.set TBD
80		#' @param FUN TBD
81		#' @param keep_sets TBD
82		#'
83		#' @return TBD
84		#' @export
85		compare.two.models = function(full.set, reduced.set, FUN = correlations_slow, keep_sets = TRUE){
86		### Arguments:
87		# dimension.type: whether you would like to run comparison for MR dimension ("MR") or regular first SVD dimension.
88		# full.set (or model 0 set): you need to pass an ena.set object. It could be mean rotation or regular ENA.
89		# reduced.set (or model alternative set): you need to pass an ena.set object. It could be mean rotation or regular ENA.
90
91		### Return:
92		# upper or lower boundary of C.I.s for three criteria (1. goodneess of fit, 2. location of points, 3. location of nodes)
93	25x	full.points <- as.matrix(full.set$points)[,1]
94	25x	full.nodes <- as.matrix(full.set$rotation$nodes)[,1]
95	25x	reduced.points <- as.matrix(reduced.set$points)[,1]
96	25x	reduced.nodes <- as.matrix(reduced.set$rotation$nodes)[,1]
97	25x	full.codes <- full.set$rotation$codes
98	25x	reduced.codes <- reduced.set$rotation$codes
99
100		### 1. Compare goodness of fit
101	25x	cor.set.res <- FUN(reduced.set);
102	25x	cor.set.gof = round(cor.set.res$corr_estimates, digits = 4)
103	25x	lower.CI.cor.set.gof = round(cor.set.res$corr_lower_bound, digits = 4)
104	25x	upper.CI.cor.set.gof = round(cor.set.res$corr_upper_bound, digits = 4)
105
106		### 2. Compare point locations
107	25x	cor.points <- cor.test(full.points, reduced.points)
108	25x	lower.CI.cor.points <- round(as.numeric(cor.points$conf.int[1]), digits = 4)
109	25x	upper.CI.cor.points <- round(as.numeric(cor.points$conf.int[2]), digits = 4)
110
111		### 3. Compare node positions
112	25x	full.nodes.reduced <- full.nodes[full.codes %in% reduced.codes]
113	25x	cor.nodes <- cor.test(full.nodes.reduced, reduced.nodes)
114	25x	lower.CI.cor.nodes <- round(as.numeric(cor.nodes$conf.int[1]), digits = 4)
115	25x	upper.CI.cor.nodes <- round(as.numeric(cor.nodes$conf.int[2]), digits = 4)
116
117	25x	result = list(
118	25x	"num.removed.codes" = length(full.codes) - length(reduced.codes),
119	25x	"goodness.of.fit" = cor.set.gof,
120	25x	"goodness.of.fit.ci_lower" = lower.CI.cor.set.gof,
121	25x	"goodness.of.fit.ci_upper" = upper.CI.cor.set.gof,
122	25x	"points" = as.numeric(cor.points$estimate),
123	25x	"points.ci_lower" = lower.CI.cor.points,
124	25x	"points.ci_upper" = upper.CI.cor.points,
125	25x	"nodes" = as.numeric(cor.nodes$estimate),
126	25x	"nodes.ci_lower" = lower.CI.cor.nodes,
127	25x	"nodes.ci_upper" = upper.CI.cor.nodes
128	25x	,"codes.used" = list(reduced.codes)
129		)
130
131	25x	attr(result, "codes.removed") <- as.character(reduced.codes);
132	25x	attr(result, "codes.maintained") <- as.character(full.codes %in% reduced.codes);
133	25x	if (keep_sets == TRUE) {
134	9x	attr(result, "ena.set") <- reduced.set;
135		}
136		# class(result) <- c("pria_result", class(result))
137	25x	return(result)
138		}
139		#' Compute GOF and cors
140		#'
141		#' @param full.set TBD
142		#' @param reduced.set TBD
143		#' @param FUN TBD
144		#' @param keep_sets TBD
145		#'
146		#' @return TBD
147		#' @export
148		compare.two.models_cors = function(full.set, reduced.set, FUN = correlations, keep_sets = TRUE){
149		### Arguments:
150		# dimension.type: whether you would like to run comparison for MR dimension ("MR") or regular first SVD dimension.
151		# full.set (or model 0 set): you need to pass an ena.set object. It could be mean rotation or regular ENA.
152		# reduced.set (or model alternative set): you need to pass an ena.set object. It could be mean rotation or regular ENA.
153
154		### Return:
155		# upper or lower boundary of C.I.s for three criteria (1. goodneess of fit, 2. location of points, 3. location of nodes)
156	!	full.points <- as.matrix(full.set$points)[,1]
157	!	full.nodes <- as.matrix(full.set$rotation$nodes)[,1]
158	!	reduced.points <- as.matrix(reduced.set$points)[,1]
159	!	reduced.nodes <- as.matrix(reduced.set$rotation$nodes)[,1]
160	!	full.codes <- full.set$rotation$codes
161	!	reduced.codes <- reduced.set$rotation$codes
162
163		### 1. Compare goodness of fit
164	!	cor.set.res <- FUN(reduced.set);
165	!	cor.set.gof = round(cor.set.res$corr_estimates, digits = 4)
166	!	lower.CI.cor.set.gof = round(cor.set.res$corr_lower_bound, digits = 4)
167	!	upper.CI.cor.set.gof = round(cor.set.res$corr_upper_bound, digits = 4)
168
169		### 2. Compare point locations
170	!	cor.points <- cor.test(full.points, reduced.points)
171	!	lower.CI.cor.points <- round(as.numeric(cor.points$conf.int[1]), digits = 4)
172	!	upper.CI.cor.points <- round(as.numeric(cor.points$conf.int[2]), digits = 4)
173
174		### 3. Compare node positions
175	!	full.nodes.reduced <- full.nodes[full.codes %in% reduced.codes]
176	!	cor.nodes <- cor.test(full.nodes.reduced, reduced.nodes)
177	!	lower.CI.cor.nodes <- round(as.numeric(cor.nodes$conf.int[1]), digits = 4)
178	!	upper.CI.cor.nodes <- round(as.numeric(cor.nodes$conf.int[2]), digits = 4)
179
180	!	result = list(
181	!	"num.removed.codes" = length(full.codes) - length(reduced.codes),
182	!	"goodness.of.fit" = cor.set.gof,
183	!	"goodness.of.fit.ci_lower" = lower.CI.cor.set.gof,
184	!	"goodness.of.fit.ci_upper" = upper.CI.cor.set.gof,
185	!	"points" = as.numeric(cor.points$estimate),
186	!	"points.ci_lower" = lower.CI.cor.points,
187	!	"points.ci_upper" = upper.CI.cor.points,
188	!	"nodes" = as.numeric(cor.nodes$estimate),
189	!	"nodes.ci_lower" = lower.CI.cor.nodes,
190	!	"nodes.ci_upper" = upper.CI.cor.nodes
191	!	,"codes.used" = list(reduced.codes)
192		)
193
194	!	attr(result, "codes.removed") <- as.character(reduced.codes);
195	!	attr(result, "codes.maintained") <- as.character(full.codes %in% reduced.codes);
196	!	if (keep_sets == TRUE) {
197	!	attr(result, "ena.set") <- reduced.set;
198		}
199		# class(result) <- c("pria_result", class(result))
200	!	return(result)
201		}
202		#' Compare to models by distance
203		#'
204		#' @param full.set TBD
205		#' @param reduced.set TBD
206		#' @param FUN TBD
207		#' @param keep_sets TBD
208		#'
209		#' @return TBD
210		#' @export
211		compare.two.models.by.dist = function(full.set, reduced.set, keep_sets = TRUE,n_dim=2){
212		### Arguments:
213		# dimension.type: whether you would like to run comparison for MR dimension ("MR") or regular first SVD dimension.
214		# full.set (or model 0 set): you need to pass an ena.set object. It could be mean rotation or regular ENA.
215		# reduced.set (or model alternative set): you need to pass an ena.set object. It could be mean rotation or regular ENA.
216
217		### Return:
218		# upper or lower boundary of C.I.s for three criteria (1. goodneess of fit, 2. location of points, 3. location of nodes)
219	!	full.points <- as.matrix(full.set$points)
220	!	reduced.points <- as.matrix(reduced.set$points)[,1:n_dim]
221	!	if(n_dim>0)
222		{
223	!	full.points=full.points[,1:min(n_dim,ncol(full.points))]
224	!	reduced.points=reduced.points[,1:min(n_dim,ncol(reduced.points))]
225		}
226
227	!	full.codes <- full.set$rotation$codes
228	!	reduced.codes <- reduced.set$rotation$codes
229	!	dist_corr_score <-dist_corr(full.points,reduced.points)
230
231	!	result = list(
232	!	"num.removed.codes" = length(full.codes) - length(reduced.codes),
233	!	"score" = dist_corr_score,
234	!	"codes.used" = list(reduced.codes)
235		)
236
237	!	attr(result, "codes.removed") <- as.character(reduced.codes);
238	!	attr(result, "codes.maintained") <- as.character(full.codes %in% reduced.codes);
239	!	if (keep_sets == TRUE) {
240	!	attr(result, "ena.set") <- reduced.set;
241		}
242		# class(result) <- c("pria_result", class(result))
243	!	return(result)
244		}
245
246		dist_corr <- function(A, B, sample_size = 1000000, method = "pearson") {
247	!	m <- nrow(A)
248	!	total_possible_pairs <- choose(m, 2)
249
250	!	if (total_possible_pairs <= sample_size) {
251		# CASE 1: Small m - Compute exact correlation
252		#message("Computing exact correlation for all pairs...")
253	!	d_A <- dist(A)
254	!	d_B <- dist(B)
255	!	return(cor(d_A, d_B, method = method))
256
257		} else {
258		# CASE 2: Large m - Compute sampled correlation
259		#message(paste("Sampling", sample_size, "pairs for estimation..."))
260
261		# Generate random indices
262	!	idx1 <- sample.int(m, sample_size, replace = TRUE)
263	!	idx2 <- sample.int(m, sample_size, replace = TRUE)
264
265		# Filter out identity pairs (where idx1 == idx2)
266	!	diff_mask <- idx1 != idx2
267	!	idx1 <- idx1[diff_mask]
268	!	idx2 <- idx2[diff_mask]
269
270		# Memory-efficient distance calculation (Vectorized)
271		# Distance = sqrt(sum of squared differences)
272	!	dist_A <- sqrt(rowSums((A[idx1, , drop=FALSE] - A[idx2, , drop=FALSE])^2))
273	!	dist_B <- sqrt(rowSums((B[idx1, , drop=FALSE] - B[idx2, , drop=FALSE])^2))
274
275	!	return(cor(dist_A, dist_B, method = method))
276		}
277		}
278		#' Conmpare
279		#'
280		#' @param A matrix, m by n1
281		#' @param B matrix, m by n2
282		#' @param n_dims number of first dimensions
283		#'
284		#' @return space similarity in terms of weighted correlation
285		#' @export
286		compute_weighted_average_correlation <- function(A, B, n_dims = NULL) {
287		# A: Original Space (m x n1)
288		# B: Reduced Space (m x n2)
289		# k: Number of dimensions to compare (defaults to all available)
290
291		# 1. Determine how many dimensions we can actually compare
292	!	k = n_dims
293	!	max_k <- min(ncol(A), ncol(B))
294	!	if (is.null(k)) {
295	!	k <- max_k
296		} else {
297	!	k <- min(k, max_k)
298		}
299
300		# 2. Extract the first k dimensions
301	!	A_k <- A[, 1:k, drop = FALSE]
302	!	B_k <- B[, 1:k, drop = FALSE]
303
304		# 3. Calculate weights based on variance (Singular Values squared)
305		# This ensures Dim 1 is more important than Dim 2
306	!	weights <- apply(A_k, 2, var)
307
308		# 4. Compute 1-to-1 correlation for each dimension
309	!	r2_values <- numeric(k)
310	!	for (j in 1:k) {
311		# We use cor()^2 to get R-squared
312		# use = "complete.obs" handles potential NAs
313	!	r2_values[j] <- cor(A_k[, j], B_k[, j], use = "complete.obs")^2
314		}
315
316		# 5. Return the Weighted Average
317	!	weighted_r2 <- sum(r2_values * weights) / sum(weights)
318
319	!	return(list(
320	!	similarity = weighted_r2,
321	!	dimension_r2 = r2_values,
322	!	weights = weights / sum(weights) # Normalized weights for reporting
323		))
324		}

1		#' @noRd
2		pria_backward_greedy <- function(x, remove.num, min_gof, correlation_fun,use_dist_corr=TRUE) {
3	!	full_codes <- x$rotation$codes
4	!	current_codes <- full_codes
5	!	best_valid_set <- full_codes
6	!	best_valid_score <- get_score_internal(full_codes, x, correlation_fun, min_gof,use_dist_corr=use_dist_corr)
7
8	!	for (step in 1:remove.num) {
9	!	if (length(current_codes) <= 3) break
10	!	iteration_scores <- sapply(seq_along(current_codes), function(i) {
11	!	get_score_internal(current_codes[-i], x, correlation_fun, min_gof,use_dist_corr=use_dist_corr)
12		})
13
14	!	if (max(iteration_scores) <= 0) break
15
16	!	best_idx <- which.max(iteration_scores)
17	!	current_codes <- current_codes[-best_idx]
18	!	best_valid_set <- current_codes
19	!	best_valid_score <- max(iteration_scores)
20		}
21	!	list(best_code_set = best_valid_set, final_performance = best_valid_score)
22		}
23
24		#' @noRd
25		pria_backward_beam<- function(x, remove.num, width, min_gof, correlation_fun,use_dist_corr=TRUE) {
26		# 1. INITIALIZATION
27	!	current_codes <- x$rotation$codes
28	!	beams <- list(current_codes)
29
30		# Start with the full set
31	!	best_overall_set <- current_codes
32	!	best_overall_score <- get_score_internal(current_codes, x, correlation_fun, min_gof,use_dist_corr=use_dist_corr)
33
34	!	for (step in 1:remove.num) {
35		# Generate candidates (using lapply for speed/efficiency)
36	!	new_candidates <- list()
37	!	for (b in seq_along(beams)) {
38	!	curr <- beams[[b]]
39	!	if(length(curr) <= 3) next
40		# Add all possible 'minus-one' subsets
41	!	new_candidates <- c(new_candidates, lapply(seq_along(curr), function(i) curr[-i]))
42		}
43
44	!	if (length(new_candidates) == 0) break
45
46		# 2. UNIQUE & PRE-FILTERING (Improvement 1)
47		# Sort and remove duplicates BEFORE scoring to save time
48	!	candidate_strings <- vapply(new_candidates, function(s) paste(sort(s), collapse="\|"), character(1))
49	!	unique_idx <- which(!duplicated(candidate_strings))
50	!	new_candidates <- new_candidates[unique_idx]
51
52		# 3. EFFICIENT SCORING (Improvement 2)
53		# Using vapply instead of growing a vector with c()
54	!	candidate_scores <- vapply(new_candidates, function(s) {
55	!	get_score_internal(s, x, correlation_fun, min_gof,use_dist_corr=use_dist_corr)
56	!	}, numeric(1))
57
58		# Identify valid candidates (score > 0)
59	!	valid_idx <- which(candidate_scores > 0)
60
61		# If no valid candidates exist, we cannot remove any more codes
62	!	if (length(valid_idx) == 0) break
63
64		# Filter to only valid ones
65	!	new_candidates <- new_candidates[valid_idx]
66	!	candidate_scores <- candidate_scores[valid_idx]
67
68		# 4. PRUNING THE BEAM
69		# Even if scores are decreasing, we pick the 'width' best remaining scores
70	!	top_idx <- order(candidate_scores, decreasing = TRUE)[1:min(width, length(new_candidates))]
71	!	beams <- new_candidates[top_idx]
72
73		# 5. ALWAYS UPDATE (The "Greedy-Style" Fix)
74		# Since you want to remove as many as possible, we update the best set
75		# to be the best performer of THIS step, regardless of the previous step's score.
76	!	best_overall_set <- beams[[1]]
77	!	best_overall_score <- candidate_scores[top_idx[1]]
78		}
79
80	!	list(best_code_set = best_overall_set, final_performance = best_overall_score)
81		}
82
83		#' @noRd
84		pria_backward_pairwise <- function(x, remove.num, min_gof, correlation_fun,use_dist_corr=TRUE) {
85		# 1. INITIALIZATION
86	!	current_codes <- x$rotation$codes
87	!	best_score <- get_score_internal(current_codes, x, correlation_fun, min_gof,use_dist_corr=use_dist_corr)
88
89		# Calculate how many pairs we can remove
90		# e.g., if remove.num is 5, floor(5/2) = 2 pairs (4 codes total)
91	!	max_pair_steps <- floor(remove.num / 2)
92
93		# 2. PAIRWISE REMOVAL PHASE
94	!	if (max_pair_steps > 0) {
95	!	for (s in 1:max_pair_steps) {
96		# Safety: Need at least 5 codes to remove 2 and leave 3
97	!	if (length(current_codes) <= 4) break
98
99	!	combs <- combn(current_codes, 2)
100
101		# Evaluate removing each pair
102	!	iteration_scores <- apply(combs, 2, function(pair) {
103	!	get_score_internal(setdiff(current_codes, pair), x, correlation_fun, min_gof,use_dist_corr=use_dist_corr)
104		})
105
106		# If no pairs can be removed safely, stop the pairwise phase
107	!	if (all(iteration_scores <= 0)) break
108
109		# Update the set with the best pair removal found
110	!	best_idx <- which.max(iteration_scores)
111	!	current_codes <- setdiff(current_codes, combs[, best_idx])
112	!	best_score <- iteration_scores[best_idx]
113		}
114		}
115
116		# 3. GREEDY CLEANUP PHASE
117		# Calculate how many codes we have actually removed so far
118	!	removed_so_far <- length(x$rotation$codes) - length(current_codes)
119
120		# If we haven't reached remove.num (usually happens if remove.num is odd),
121		# try to remove the remaining required codes one by one.
122	!	while (removed_so_far < remove.num) {
123	!	if (length(current_codes) <= 3) break
124
125		# Single-item evaluation (Greedy)
126	!	step_scores <- vapply(seq_along(current_codes), function(i) {
127	!	get_score_internal(current_codes[-i], x, correlation_fun, min_gof)
128	!	}, numeric(1))
129
130		# If no further single removals are valid, stop completely
131	!	if (all(step_scores <= 0)) break
132
133		# Remove the best single item
134	!	best_step_idx <- which.max(step_scores)
135	!	current_codes <- current_codes[-best_step_idx]
136	!	best_score <- step_scores[best_step_idx]
137
138		# Update counter
139	!	removed_so_far <- length(x$rotation$codes) - length(current_codes)
140		}
141
142		# 4. FINAL OUTPUT
143	!	list(best_code_set = current_codes, final_performance = best_score)
144		}
145
146		#' @noRd
147		pria_brute_force <- function(x, remove.num, min_gof, correlation_fun, code_weights = NULL, lower_n = 12,use_dist_corr=TRUE) {
148		# 1. INITIALIZATION & FILTERING
149	!	full_codes <- x$rotation$codes
150	!	total_n <- length(full_codes)
151
152		# Logic for selecting candidate codes for removal
153	!	if (!is.null(code_weights)) {
154		# Ensure lower_n isn't larger than the total codes available
155	!	actual_lower_n <- min(lower_n, total_n)
156
157		# Sort codes by weight and pick the 'actual_lower_n' smallest
158		# Assumes code_weights is a data.frame with columns 'code' and 'weight'
159	!	sorted_weights <- code_weights[order(code_weights$weight), ]
160	!	candidate_to_remove <- sorted_weights$code[1:actual_lower_n]
161
162		# The 'fixed' codes that we definitely want to keep
163	!	fixed_codes <- setdiff(full_codes, candidate_to_remove)
164
165		# We need to remove 'remove.num' from the 'candidate_to_remove' pool
166		# Therefore, we need to choose (actual_lower_n - remove.num) to KEEP from that pool
167	!	pool_to_keep_count <- actual_lower_n - remove.num
168
169	!	if (pool_to_keep_count < 0) {
170	!	stop("remove.num cannot be greater than the number of candidate codes (lower_n).")
171		}
172
173	!	search_pool <- candidate_to_remove
174	!	target_from_pool <- pool_to_keep_count
175
176		} else {
177		# Fallback to original behavior if no weights provided
178	!	fixed_codes <- c()
179	!	search_pool <- full_codes
180	!	target_from_pool <- total_n - remove.num
181		}
182
183		# 2. COMBINATORIAL SAFETY CHECK
184	!	num_combs <- choose(length(search_pool), target_from_pool)
185
186	!	if (num_combs > 1e6) {
187	!	stop(paste("Too many combinations:", round(num_combs),
188	!	". Even with filtering, brute force is too slow."))
189		}
190
191		# 3. GENERATE SUBSETS FROM THE CANDIDATE POOL
192	!	pool_combinations <- combn(search_pool, target_from_pool, simplify = FALSE)
193
194		# 4. EVALUATE ALL SETS (Combining Fixed + Combinations)
195	!	scores <- vapply(pool_combinations, function(pool_set) {
196		# Merge the guaranteed codes with this specific combination
197	!	current_set <- c(fixed_codes, pool_set)
198	!	get_score_internal(current_set, x, correlation_fun, min_gof,use_dist_corr=use_dist_corr)
199	!	}, numeric(1))
200
201		# 5. FIND THE OPTIMUM
202	!	best_idx <- which.max(scores)
203
204	!	if (scores[best_idx] <= 0) {
205	!	warning("No valid code sets (score > 0) were found.")
206		}
207
208	!	list(
209	!	best_code_set = c(fixed_codes, pool_combinations[[best_idx]]),
210	!	final_performance = scores[best_idx],
211	!	combinations_tested = num_combs
212		)
213		}