R/hyper2.R

`hyper2` <-  function(L=list(), d=0, pnames=NA){
  if(length(d)==1){d <- rep(d,length(L))}  
  stopifnot(is_valid_hyper2(L,d,pnames))
  L <- lapply(L,as.integer)
  out <- identityL(L,d)
  out$pnames <- pnames
  class(out) <- 'hyper2'  # This is the only class assignment in the package
  return(out)
}

## Following three functions are the only accessor methods in the package
setGeneric("brackets",function(x){standardGeneric("brackets")})
setGeneric("powers"  ,function(x){standardGeneric("powers"  )})
setGeneric("pnames"  ,function(x){standardGeneric("pnames"  )})
`brackets` <- function(H){UseMethod("brackets")}
`powers`   <- function(H){UseMethod("powers"  )}
`pnames`   <- function(H){UseMethod("pnames"  )}

`brackets.hyper2` <- function(H){ H$brackets }
`powers.hyper2` <- function(H){if(is_constant(H)){
                                 return(0)
                               } else {
                                 return(H$powers)
                               }
}
  `pnames.hyper2` <- function(H){ H$pnames   }
## accessor methods end

## Following function is the only setter method in the package
setGeneric("pnames<-",function(x,value){standardGeneric("pnames<-")})
`pnames<-` <- function(x,value){UseMethod("pnames<-")}
`pnames<-.hyper2` <- function(x,value){
  hyper2(brackets(x),powers(x),pnames=value)
}
## setter methods end

`is_constant` <- function(H){ length(brackets(H))==0 }


`is_valid_hyper2` <- function(L,d,pnames){
  stopifnot(is.list(L))
  stopifnot(is.vector(d))
  stopifnot(is.numeric(d))
  stopifnot(length(L) == length(d))
  stopifnot(all(unlist(lapply(L,function(x){all(x==round(x))}))))
  if(!identical(pnames,NA)){
    if(length(L)==0){
      hsize <- 0
    } else {
      hsize <- max(c(L,recursive=TRUE))
    }
    stopifnot(length(pnames) >= hsize)
  }
  return(TRUE)
}

`size` <- function(H){
  if(is.null(H)){return(0)}
  
  if(identical(pnames(H),NA)){
    return(max(c(brackets(H),recursive=TRUE)))
  } else {
    return(length(pnames(H)))
  }
}

`as.hyper2` <- function(L,d,pnames){
  if(is.matrix(L)){
    L <- as.list(as.data.frame(t(L)))
  } else if(is.list(L)){
    ignore <- 0
  } else {
    stop("first argument must be a list or a matrix")
  }
  return(hyper2(L,d,pnames))
}

`.print.helper` <- function(pnames,vec){
  if(length(vec)==1){
    out <- pnames[vec]
  } else {
    out <- paste("(",pnames[vec[1]],sep="")
    
    for(i in vec[-1]){
      out <- paste(out," + ",pnames[i],sep="")
    }
    out <- paste(out,")",sep="")
  }
  return(out)
}
  
`print.hyper2` <- function(x,...){
  b <- brackets(x)
  powers <- powers(x)
  if(length(b)==0){  # not is.null(b)
    n <- 1
    if(identical(pnames(x),NA)){
      pn <- "1"
    } else {  # empty but with specified pnames
      pn <- pnames(x)
    }
    n <- length(pn)
    b <- list(seq_len(n))
  } else {  # non-empty
    pn <- pnames(x)
    n <- max(c(b,recursive=TRUE))   # b non-empty; n>0
    if(identical(pnames(x),NA)){
      pn <- paste("p",seq_len(n),sep="")
    } else {
      pn <- pnames(x)
    }
  }

  out <- ""
  for(i in seq_along(b)){
    jj <- unlist(b[i])
    pp <- powers[i]
    if(pp==1){
      out <- paste(out, .print.helper(pn, jj))
    } else {
      out <- paste(out, .print.helper(pn, jj), "^",pp,sep="")
    }
    if(i < length(b)){out <- paste(out," * ",sep="")}
  }
  
  out <- paste(out,"\n",sep="")
  for(i in strwrap(out)){
    cat(noquote(i))
    cat("\n")
  }
  return(invisible(x))
}

`hyper2_add` <- function(e1,e2){
  out <- addL(brackets(e1),powers(e1),brackets(e2),powers(e2))

  if(identical(pnames(e1),NA) & identical(pnames(e2),NA)){
    jj <- NA
  } else if(identical(pnames(e1),NA) & !identical(pnames(e2),NA)){
    jj <- pnames(e2)
  } else if(identical(pnames(e2),NA) & !identical(pnames(e1),NA)){
    jj <- pnames(e1)
  } else if(!identical(pnames(e2),NA) & !identical(pnames(e1),NA)){
    if(size(e2)>size(e1)){
      jj <- pnames(e2)
    } else {
      jj <- pnames(e1)
    }
  } else {
    stop("this cannot happen")
  }

  return(hyper2(out[[1]],out[[2]],pnames=jj))
}


`lhyper2` <- function(H,p,log=TRUE){
  if(is.matrix(p)){
    return(apply(p,1,function(o){.lhyper2_single(H,p=o, log=log)}))
  } else {
    return(.lhyper2_single(H,p,log=log))
  }
}

`.lhyper2_single` <- function(H,p,log=TRUE){
  stopifnot(length(p) == size(H)-1)
  stopifnot(all(p>=0))
  stopifnot(sum(p)<=1)

  out <- evaluate(brackets(H), powers(H), probs=fillup(p))
  if(log){
    return(out)
  } else {
    return(exp(out))
  }
}

`hyper2_equal` <- function(e1,e2){
  equality(
      brackets(e1),powers(e1),
      brackets(e2),powers(e2)
  )
}

`hyper2_prod` <- function(H,n){
  stopifnot(is.numeric(n))
  stopifnot(length(n)==1)
  return(hyper2(brackets(H),powers(H)*n,pnames=pnames(H)))
}

`Ops.hyper2` <-
  function (e1, e2 = NULL) 
{
  f <- function(...){stop("odd---neither argument has class hyper2?")}
  unary <- nargs() == 1
  lclass <- inherits(e1,"hyper2")
  rclass <- !unary && inherits(e2,"hyper2")
  
  if(unary){
    stop("Unary operator '", .Generic, "' is not implemented for hyper2 objects")
    }

  if (!is.element(.Generic, c("+", "==", "!=", "*", "^" ))){
    stop("Binary operator '", .Generic, "' is not implemented for hyper2 objects")
  }
  
  if(lclass && rclass){  
    if (.Generic == "+"){
      return(hyper2_add(e1,e2))
    } else if (.Generic == "==") {
      return(hyper2_equal(e1, e2))
    } else if (.Generic == "!=") {
      return(!hyper2_equal(e1, e2))
    } else {
      stop("H1 ", .Generic, " H2 not defined")
    }
  } else {  # H * n
    stopifnot(is.element(.Generic, c("*")))
    if(lclass && !rclass){
      return(hyper2_prod(e1,e2))
    } else if (!lclass && rclass){
      return(hyper2_prod(e2,e1))
    } else {
    stop("method not defined for hyper2")
    }
  }
}
  
`character_to_number` <- function(char,pnames){ # char = c("a","b","D")
    stopifnot(all(char %in% pnames))
    unlist(apply(outer(char,pnames, "=="),1,which))
    }

`char2num` <- character_to_number

`[.hyper2` <- function(x, ...){
    dots <- list(...)
    first <- dots[[1]]

    if(nargs() > 2){  # something like H[3,4]
      wanted <- list(c(dots,recursive=TRUE))
    } else if(is.list(first)){
      wanted <- dots[[1]]
    } else if(is.matrix(first)){
      wanted <- as.list(as.data.frame(t(first)))
    } else if(is.vector(first)){
      wanted <- list(first)
    } else {
      wanted <- dots
    }

    if(any(unlist(lapply(wanted, is.character)))){
        wanted <- lapply(wanted, character_to_number, pnames=pnames(x))
    }
    
    out <- accessor(x[[1]],x[[2]],wanted)
    return(hyper2(out[[1]],out[[2]],pnames=pnames(x)))
}

`.assign_lowlevel`<- function(x,index,value){ #H[index] <- value

    stopifnot(class(x) == 'hyper2')
    
    if(is.list(index)){
        ignore <- 3  # 'index' is supposed to be a list
    } else if(is.matrix(index)){
        index <- as.list(as.data.frame(t(index)))
    } else if(is.vector(index)){
        index <- list(index)
    } else {
        stop("replacement index must be a list, a matrix, or a vector")
    }

    if(any(unlist(lapply(index, is.character)))){
        index <- lapply(index, character_to_number, pnames=pnames(x))
    }

    stopifnot(is.numeric(value)) # coercion to integer is done in C
    stopifnot(is.vector(value))
    if(length(value)==1){
        value <- rep(value, length(index))
    }
    return(assigner(brackets(x),powers(x),index,value))
}

`.overwrite_lowlevel` <- function(x,value){
  stopifnot(class(x)     == 'hyper2')
  stopifnot(class(value) == 'hyper2')

  overwrite(brackets(x), powers(x), brackets(value), powers(value))
}

`[<-.hyper2` <- function(x, index, ..., value){
  if(missing(index)){  # A[] <- B
    out <- .overwrite_lowlevel(x,value)
  } else {
    out <- .assign_lowlevel(x,index,value)
  }
    return(hyper2(out[[1]],out[[2]],pnames=pnames(x)))
}

`gradient` <- function(H,probs){
  stopifnot(length(probs) == size(H)-1)
  differentiate(brackets(H), powers(H), fillup(probs), size(H))$grad_comp
}

`fillup` <- function(x){
  if(is.matrix(x)){
    return(cbind(x,1-rowSums(x)))
  } else {  # assumed to be a vector
    return(c(x,1-sum(x)))
  }
}

`indep` <- function(x){
   if(is.matrix(x)){
     return(x[,-ncol(x)])
  } else {  # assumed to be a vector
    return(x[-length(x)])
  }
}

`maxp` <- function(H, startp=NULL, give=FALSE, ...){
    SMALL <- 1e-6
    
    n <- size(H)
    if(is.null(startp)){
        startp <- rep(1/n,n-1)
        }

    objective <- function(p){ -lhyper2(H,p) }
    gradfun   <- function(p){ -(gradient(H,p))} #NB minus signs
    
    UI <- rbind(diag(nrow=n-1),-1)
    CI <- c(rep(SMALL,n-1),-1+SMALL)
    
    out <- constrOptim(
        theta = startp,
        f     = objective,
        grad  = gradfun,
        ui    = UI,
        ci    = CI,
        ...)
        
    
    if(give){
      return(out)
    } else {
      jj <- fillup(out$par)
      if(!identical(pnames(H),NA)){names(jj) <- pnames(H)}
      return(jj)
    }
}

`sum.hyper2` <- function(x, ..., na.rm=FALSE){
  if(nargs()==1){
    return(x)
  } else if (nargs()==2){
    return(hyper2_add(x, ...))
  } else {
    return(hyper2_add(x, Recall(...)))
  }
}

`head.hyper2` <- function(x,...){ x[head(brackets(x),...)] }

`order_likelihood` <- function(M,times=1){
  M <- rbind(M)  # deals with vectors
  times <- cbind(seq_len(nrow(M)),times)[,2]
  if(is.null(colnames(M))){
    cn <- NA
  } else {
    cn <- colnames(M)
  }
  out <- hyper2(pnames=cn)
  for(i in seq_len(nrow(M))){
    v <- rev(M[i,,drop=TRUE])
    jj1 <- hyper2(as.list(v),times[i])   # numerators
    jj2 <- hyper2(sapply(seq_along(v),function(i){v[seq_len(i)]}),-times[i]) # denominators
    out <- out + jj1 + jj2
  }  # i loop closes
  return(out)
}

#`addrank` <- function(H,ranks){
#     .Defunct(new = 'order_likelihood',msg='use H <- H+order_likelihood(...)')
#     ranks <- rbind(ranks)
# 
#     for(r in seq_len(nrow(ranks))){
#         rank <- rev(ranks[r,])
#         for(i in seq_along(rank)){
#             H[rank[i]] <- powers(H[rank[i]]) + 1
#             H[rank[seq_len(i)]] <- powers(H[rank[seq_len(i)]])-1
#         }
#     }
#     return(H)
# }

`.rrank_single` <- function(p){
    ell <- length(p)
    r <- integer(ell)
    i <- 1   # 'i' means placing.
    while(any(p>0)){  
        m <- sample(ell,1,prob=p) #randomly chosen competitor
        r[i] <- m   # place 'i' is competitor 'm'  
        p[m] <- 0   # competitor 'm' can't place again.
        i <- i+1    # consider the next place
    }
    return(r)
}

`rrank` <- function(n=1, p, pnames=NULL, fill=FALSE){
    if(fill){ p <- fillup(p) }
    out <- t(replicate(n, .rrank_single(p)))
    colnames(out) <- pnames
    class(out) <- "rrank"
    return(drop(out))
}

`print.rrank` <- function(x, ...){
  if(is.null(colnames(x))){
    cn <- seq_len(ncol(x))
  } else {
    cn <- colnames(x)
  }
  thing <- noquote(matrix(cn[x],ncol=ncol(x)))
  colnames(thing) <- paste("c",seq_len(ncol(x)),sep="")
  print(thing)
  return(NA) # sic: this is confusing!
}

`.allorders` <- function(x){
  out <- perms(length(x))
  matrix(x[out],nrow=length(x))
}

`pair_grid` <- function(a,b){  ## sort of a  generalized expand.grid()
  rbind(
      kronecker(a,t(rep(1L,ncol(b)))),
      kronecker(t(rep(1L,ncol(a))),b))
}

`mult_grid` <- function(L){
  if(length(L) == 0){
    return(1)
  } else if(length(L)==1){
    return(L)
  } else {
    return(Recall(c(list(pair_grid(L[[1]],L[[2]])),L[-(1:2)])))
  }
}

`general_grouped_order_likelihood` <- function(H, ...){
  ## ggol(1:3,4,5:9) means group(1:3) came first, 4 in the middle, 5:9
  ## last
    dotargs <- list(...)
    if(any(unlist(lapply(dotargs, is.character)))){
        dotargs <- lapply(dotargs, character_to_number, pnames=pnames(H))
    }

    if(any(table(c(dotargs,recursive=TRUE))>1)){
        stop('repeated competitor')
    }

    out <- mult_grid(lapply(dotargs, .allorders))[[1]]
    out <- apply(out,2,function(rank){H+order_likelihood(rank)})
  return(out)
}

`ggol` <- general_grouped_order_likelihood

`choose_losers` <- function(H,all_players,losers){
  stopifnot(losers %in% all_players)
  winners <- all_players[!all_players %in% losers]
  ggol(H,winners,losers)
}

`choose_winners` <- function(H,all_players, winners){
  stopifnot(all(winners %in% all_players))
  losers <- all_players[!all_players %in% winners]
  ggol(H,winners,losers)
}

`like_single_list` <- function(p,Lsub){;  # eg like_single(p,Lc)
  sum(unlist(lapply(Lsub,lhyper2,p=p,log=FALSE)))
  ## sum, because it can be any one of the orders specified in the
  ## elements of Lsub
}

`like_series` <- function(p,L,log=TRUE){
    out <- prod(unlist(lapply(L, function(Lsub){like_single_list(p,Lsub)})))
    if(log){ out <- log(out) }
    return(out)
}

Try the hyper2 package in your browser

Any scripts or data that you put into this service are public.

hyper2 documentation built on July 6, 2017, 9:02 a.m.