R/check.equivalence.R

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

#' Checking equivalence of a set of splits, and reporting each equivalent class of splits as a single split
#'
#' \code{check.equivalence()} is used in function \code{zeroRDD_to_splits()}. It is an internal function. This function is used to check the equivalence of each split, and then report each equivalent class of splits as a single split in a matrix \code{split_set_indicator}. \code{split_set_indicator} is a matrix with \eqn{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 array_zero_ID an array where zero valued RD differences have a value zero, and the rest have a value of 1. Note that the array is P X P X P, and for each i, array_zero_ID[i,,] is symmetric. This array could be generated by function `BMaf_to_zeroRDD`
#'
#' @return  A matrix \code{split_set_indicator}, which indicates each participation of each taxon (by column) in each split (by row).
#' @noRd
#'
#' @examples NA
#'
check.equivalence<-function(array_zero_ID){

  P<-dim(array_zero_ID)[1]
  increment <- 2*P
  max_split <- increment
  split_set_indicator <- matrix(nrow=max_split,ncol=P)
  n_split_set <- 0
  # This will be incremented as we find new splits
  for (i in 1:P)
    for (j in 1:P)
      for (k in 1:P)
      {

        if ( ((i != j) && (i != k)) && ((j < k) && (array_zero_ID[i,j,k] == 0)))
        {

          i_split <- 0
          repeat
          {

            ## i_split is counting the number of already-identified splits compared
            i_split <- i_split + 1

            if (i_split > n_split_set)
            {
              n_split_set <- i_split

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

              this_split_set_indicator <- integer(P)
              this_split_set_indicator[i] <- 1
              this_split_set_indicator[j] <- 2
              this_split_set_indicator[k] <- 2

              split_set_indicator[n_split_set,] <- this_split_set_indicator

              break
            }
            ## that is, if no other equivalent split is found, define (i,j,k) as a new split; then get out

            if ( (split_set_indicator[i_split,i] > 0) && (split_set_indicator[i_split,k]==0) && (( (split_set_indicator[i_split,j] > 0)  &&  (split_set_indicator[i_split,j] !=  split_set_indicator[i_split,i]) ) ) )
            {
              split_set_indicator[i_split,k] <- split_set_indicator[i_split,j]

              break
            }
            else if ( (split_set_indicator[i_split,i] > 0) && (split_set_indicator[i_split,j]==0) && ( ( (split_set_indicator[i_split,k] > 0)  &&  (split_set_indicator[i_split,k] !=  split_set_indicator[i_split,i]) ) ) )
            {
              split_set_indicator[i_split,j] <- split_set_indicator[i_split,k]

              break
            }

          }
        }

      }

  split_set_indicator <- matrix(ncol=P,split_set_indicator[(1:n_split_set),])

  return(split_set_indicator)

}