R/EquivalenceClass.R

In MCARtest: Optimal Nonparametric Testing of Missing Completely at Random

#' Simplifies H-representation by exploiting symmetry
#'
#' The marginal polytope and related objects have many symmetries. By relabelling the levels
#' of discrete variables we transform facets into other facets. This function reduces a list of
#' halfspace normals to its equivalence classes.
#'
#' @param bS A binary matrix specifying the set of observation patterns. Each row encodes a single pattern.
#' @param M A vector of positive integers giving the alphabet sizes of the discrete variables.
#' @param Hrep An H-representation generated by \code{MargPolyHrep}, \code{ConsMinkSumHrep} or \code{InconsMinkSumHrep}.
#'
#' @return A list of representative halfspace normals.
#' @export
#'
#' @importFrom rcdd q2d
#' @importFrom gtools permutations
#' @importFrom Epi in.span
#'
#' @examples
#' bS=matrix(c(1,1,0, 1,0,1, 0,1,1),byrow=TRUE,ncol=3) # Our canonical 3d example
#' Hrep=MargPolyHrep(bS,c(2,2,2))
#' EquivalenceClass(bS,c(2,2,2),Hrep)
#'


EquivalenceClass=function(bS,M,Hrep){
  d=length(M); v=ncol(Hrep)-2; cardS=nrow(bS)
  l=Hrep[,1]; LinearConstr=Hrep[l=="1",]; A2=q2d(t(LinearConstr[,-c(1,2)]))
  IneqConstr=Hrep[l=="0",]; NIneq=nrow(IneqConstr); A1t=q2d(IneqConstr[,-c(1,2)])
  Classes=A1t[1,]; new=rep(1,NIneq)
  for(nineq in 2:NIneq){
    fS=A1t[nineq,]
    for(nperm in 1:prod(factorial(M))){
      nperms=as.vector(arrayInd(nperm,factorial(M))); Perms=matrix(rep(0,d*max(M)),nrow=d)
      for(j in 1:d){
        Perms[j,1:M[j]]=permutations(M[j],M[j])[nperms[j],]
        #Perms[[j]]=as(permvector, "pMatrix")
      }
      fSPerm=fS-fS
      for(S in 1:cardS){
        MS=M[bS[S,]==1]; inds=which(bS[S,]==1)
        for(i in 1:prod(MS)){
          xS=as.vector(arrayInd(i,MS)); xSPerm=xS
          for(j in 1:length(xS)){
            xSPerm[j]=Perms[inds[j],xS[j]]
          }
          fSPerm[row_index(bS,M,S,xSPerm)]=fS[row_index(bS,M,S,xS)]
        }
      }
      if(sum(apply(as.matrix(Classes), 1, function(x, want) isTRUE(in.span(A2, x-want)), fSPerm))>0){
        new[nineq]=0
      }
    }
    if(new[nineq]==1) Classes=rbind(Classes,fS)
    # In this case fS is not equivalent to anything above it.
  }
  return(rbind(LinearConstr,IneqConstr[new==1,]))
}