R/eliminate.rep.R

In rapidphylo: Rapidly Estimates Phylogeny from Large Allele Frequency Data Using Root Distances Method

#' Eliminating repetitions in split set indicator matrix
#'
#' This function is used in \code{zeroRDD_to_splits()}. This function eliminates the repeated splits in the split set indicator matrix, a matrix with P columns, where each row represents a split. In a given row, a k-way split is represented by assigning each set of taxa that is descendent of a split as an identifying number from 1 to k, and numbering the position corresponding to each taxa that is part of such a set by that identifying number. The other positions are set to zero.
#'
#' @param split_set_indicator Generated by \code{check.equivalence()} function
#'
#' @return A matrix \code{new_split_set_indicator}, a matrix with P columns, where each row represents a split, but eliminating repeating splits. In a given row, a k-way split is represented by assigning each set of taxa that is descendent of a split as an identifying number from 1 to k, and numbering the position corresponding to each taxa that is part of such a set by that identifying number. The other positions are set to zero.
#' @noRd
#'
#' @examples NA
eliminate.rep<-function(split_set_indicator){

  P<-ncol(split_set_indicator)
  n_split_set<-nrow(split_set_indicator)
  increment <- P-1
  max_split <- increment
  next_split_set_indicator <- matrix(nrow=max_split,ncol=P)
  next_n_split_set <- 0
  merge_record <- NULL

  for (i_split in 1:n_split_set)
  {

    j_split <- 0

    repeat
    {
      j_split <- j_split + 1

      if (j_split > next_n_split_set)
      {


        next_n_split_set <- next_n_split_set + 1


        if (next_n_split_set > max_split)
        {
          max_split <- max_split + increment
          next_array <- matrix(nrow=max_split,ncol=P)
          next_array[(1:(max_split-increment)),] <- next_split_set_indicator
          next_split_set_indicator <- next_array
          rm(next_array)
        }



        next_split_set_indicator[next_n_split_set,] <- split_set_indicator[i_split,]

        break
      }

      if (  ((length( intersect( which(next_split_set_indicator[j_split,] == 1),which(split_set_indicator[i_split,] == 1) ) ) > 0) && (length( intersect( which(next_split_set_indicator[j_split,] == 2),which(split_set_indicator[i_split,] == 2) ) ) > 0) )
            &&    ((length( intersect( which(next_split_set_indicator[j_split,] == 1),which(split_set_indicator[i_split,] == 2) ) ) == 0) && (length( intersect( which(next_split_set_indicator[j_split,] == 2),which(split_set_indicator[i_split,] == 1) ) ) == 0) )  )
      {


        next_split_set_indicator[j_split,(union( which(next_split_set_indicator[j_split,] == 1),which(split_set_indicator[i_split,] == 1) ))] <- 1
        next_split_set_indicator[j_split,(union( which(next_split_set_indicator[j_split,] == 2),which(split_set_indicator[i_split,] == 2) ))] <- 2


        break
      }
      else if (  ((length( intersect( which(next_split_set_indicator[j_split,] == 1),which(split_set_indicator[i_split,] == 2) ) ) > 0) && (length( intersect( which(next_split_set_indicator[j_split,] == 2),which(split_set_indicator[i_split,] == 1) ) ) > 0) )
                 &&    ((length( intersect( which(next_split_set_indicator[j_split,] == 1),which(split_set_indicator[i_split,] == 1) ) ) == 0) && (length( intersect( which(next_split_set_indicator[j_split,] == 2),which(split_set_indicator[i_split,] == 2) ) ) == 0) )  )
      {


        next_split_set_indicator[j_split,(union( which(next_split_set_indicator[j_split,] == 1),which(split_set_indicator[i_split,] == 2) ))] <- 1
        next_split_set_indicator[j_split,(union( which(next_split_set_indicator[j_split,] == 2),which(split_set_indicator[i_split,] == 1) ))] <- 2


        break
      }
    }
  }

  new_split_set_indicator <- matrix(nrow=max_split,ncol=P)


  new_n_split_set <- 0


  for (i_split in 1:next_n_split_set)
  {

    j_split <- 0

    repeat
    {
      j_split <- j_split + 1

      if (j_split > new_n_split_set)
      {


        new_n_split_set <- new_n_split_set + 1


        if (new_n_split_set > max_split)
        {
          max_split <- max_split + increment
          new_array <- matrix(nrow=max_split,ncol=P)
          new_array[(1:(max_split-increment)),] <- new_split_set_indicator
          new_split_set_indicator <- new_array
          rm(new_array)
        }



        new_split_set_indicator[new_n_split_set,] <- next_split_set_indicator[i_split,]

        break
      }

      if (setequal(which(new_split_set_indicator[j_split,] == 1),which(next_split_set_indicator[i_split,] == 1)))
      {
        new_split_set_indicator[j_split,(union( which(new_split_set_indicator[j_split,] == 2),which(next_split_set_indicator[i_split,] == 2) ))] <- 2

        break
      }
      else if (setequal(which(new_split_set_indicator[j_split,] == 2),which(next_split_set_indicator[i_split,] == 2)))
      {
        new_split_set_indicator[j_split,(union( which(new_split_set_indicator[j_split,] == 1),which(next_split_set_indicator[i_split,] == 1) ))] <- 1

        break
      }
      else if (setequal(which(new_split_set_indicator[j_split,] == 1),which(next_split_set_indicator[i_split,] == 2)))
      {
        new_split_set_indicator[j_split,(union( which(new_split_set_indicator[j_split,] == 2),which(next_split_set_indicator[i_split,] == 1) ))] <- 2

        break
      }
      else if (setequal(which(new_split_set_indicator[j_split,] == 2),which(next_split_set_indicator[i_split,] == 1)))
      {
        new_split_set_indicator[j_split,(union( which(new_split_set_indicator[j_split,] == 1),which(next_split_set_indicator[i_split,] == 2) ))] <- 1

        break
      }

    }
  }

  return(new_split_set_indicator)
}