R/toporanga.R

In toporanga: Topological Sort-Based Hierarchy Inference

# Load the native code
#' @useDynLib toporanga, .registration=TRUE
NULL

#' Example EcoHAB data
#'
#' Mice dominance events recorded by the Eco-HAB system.
#' Covers a group of 10 mice and 6575 events recorded over a span of almost 5 days, corresponding to followings in the system.
#' Contains three columns: \code{winner}, ID of winning mouse, \code{loser}, ID of a losing mice, \code{t} time of event occurrence, in seconds since the first event.
#' @usage data(EcoHAB)
"EcoHAB"


#' A catch-all function to access package functionality in one call
#'
#' The general pipeline of toporanga is to tally the domination events, optionally modify the tally by diffusing event weights and/or re-scoring reciprocal weights, arrange agents into an order represented by an acyclic graph, finally to reduce it to an order or a dominance score.
#' This function collects all of this into a single call.
#'
#' @param winning vector of IDs of winning agents.
#' @param opposing vector of IDs of opposing agents.
#' @param weight optional weight of the event; set to 1 for each event if not given.
#' @param agents optional vector of IDs of all agents; allows for ghost agents that are in the set but are not involved in any event.
#' @param diffuse if \code{TRUE}, weight diffusion is applied. See \code{diffuse} for more details.
#' @param epochs optional driver for applying analysis on subsets of events. See \code{cumulative_epochs} and \code{window_epochs} for more information.
#' If given, changes the output to an aggregated version of whatever was selected with \code{output}. 
#' @param reciprocals a vector of commands applied to the tally.
#' \code{"keep"} does nothing;
#' \code{"marginise"} reduces weights into weight margins, see \code{marginise},
#' \code{"zero"} zeroes opposing weights, see \code{zero_opposition},
#' \code{"normalise"} normalises weights within conflicts, see \code{normalise};
#' \code{"diffuse"} applies \code{diffuse}, allows to manage order in which diffuse is mixed with other modifiers; cannot be mixed with \code{diffuse=TRUE}.
#' Normally, diffuse is applied first.
#' @param output specifies what to return from the function.
#' \code{"parameters"} returns the some agent parameters, see \code{parameters};
#' \code{"graph"} a \code{toporanga_graph} object;
#' \code{"order_sample"} an order of agents; in case multiple is possible, a single random sample;
#' \code{"subordinate_count"} a count of subordinate agents for each agent, which can be interpreted as a dominance score. 
#' @param ... ignored.
#' @returns A result of the caclulation in a form specified by the \code{output} argument, aggregated over epochs is \code{epochs} argument is given.
#' @examples data(EcoHAB)
#' toporanga(EcoHAB$winner,EcoHAB$loser)
#' @export
toporanga<-function(winning,opposing,weight,
  ...,
  agents,
  diffuse=FALSE,
  reciprocals=c(),
  output=c("parameters","graph","order_sample","subordinate_count"),
  epochs
  ){

  reciprocals<-match.arg(
    reciprocals,
    choices=c("keep","marginise","zero","normalise","diffuse"),
    several.ok=TRUE
  )
  output<-match.arg(output,several.ok=FALSE)

  if(diffuse && "diffuse"%in%reciprocals) stop("Diffuse given and diffuse in reciprocals")

  if(!missing(epochs)){
    winning<-as.character(winning)
    opposing<-as.character(opposing)
    stopifnot(all(!is.na(winning)))
    stopifnot(all(!is.na(opposing)))
    stopifnot(length(winning)==length(opposing))
    m<-length(winning)
    if(missing(weight)){
      weight<-rep(1,m)
    }else{
      stopifnot(all(is.finite(weight)))
      stopifnot(length(weight)==length(winning))
    }
    if(missing(agents))
      agents<-sort(unique(c(unique(winning),unique(opposing))))
    
    epochs(m)->es
    es$n->n_epochs
    es$t->t
    lapply(1:n_epochs,function(e){
      epo<-epochs(m,e)
      mask<-epo$mask
      toporanga(
        winning[mask],opposing[mask],weight[mask],
        agents=agents,diffuse=diffuse,reciprocals=reciprocals,output=output
      )->ans
      if(output%in%c("parameters","graph")){
        ans$Epoch<-e
        for(k in names(epo$meta)) ans[[k]]<-epo$meta[[k]]
        ans
      }else ans
    })->ans
    return(switch(
      output,
      "parameters"=do.call(rbind,ans),
      "graph"=ans,
      "order_sample"=data.frame(do.call(rbind,lapply(ans,'[',agents))),
      "subordinate_count"=data.frame(do.call(rbind,lapply(ans,'[',agents))),
    ))
  }

  tally_log(winning,opposing,weight=weight,agents=agents)->tally
  for(e in reciprocals)
    tally<-switch(
      e,
      "keep"=tally,
      "marginise"=marginise(tally),
      "zero"=zero_opposition(tally),
      "normalise"=normalise(tally),
      "diffuse"=diffuse(tally)
    )

  arrange(tally)->graph

  parameter_pull<-function(x,p){
    parameters(x)->g
    stats::setNames(g[,p],g$Agent)
  }

  ans<-switch(
    output,
    "parameters"=parameters(graph),
    "graph"=graph,
    "order_sample"=parameter_pull(graph,"Order"),
    "subordinate_count"=parameter_pull(graph,"Subordinates")
  )

  ans
}

#' Cumulative epoch generator
#'
#' Splits the event log into cumulative epochs, i.e., starting from the start up to selected number of final points.
#' Returns an object which has to be fed into \code{epochs} argument of \code{toporanga} function.
#' @param n number of epochs.
#' Silently capped to the number of events, if larger.
#' @param t optional event times, vector of a same length and order as the event log given to \code{toporanga} function.
#' When given, epochs will be generated based on equal-time breaks not equal-event-number breaks.
#' @note If time is not given, event log must be sorted; otherwise, epochs will not make any sense.
#' @returns A special function that can be passed to the \code{epochs} argument of \code{toporanga} function.
#' @export
cumulative_epochs<-function(n,t){
  if(missing(t)){
    t<-NULL
  }else{
    stopifnot(is.numeric(t))
    stopifnot(all(is.finite(t)))
    tr<-diff(range(t))
    tm<-min(t)
  }
  function(m,e){
    if(missing(e)){
      #Set-up
      if(!is.null(t)){
        if(length(t)!=m) stop("Time vector doesn't have the same length as event count ",m)
      }
      if(n>m){
        E<-parent.env(environment())
        assign("n",m,envir=E)
      }
      return(list(n=n,t=t))
    }else{
      if(is.null(t)){
        round((e-1)/(n-1)*(m-1))+1->end
        return(list(mask=1:end,meta=list(StartEvent=1,EndEvent=end)))
      }else{
        te<-tm+tr*e/n
        mask<-t<=te
        return(list(mask=mask,meta=list(StartTime=tm,EndTime=te,Events=sum(mask))))
      }
    }
  }
}

#' Moving window epoch generator
#'
#' Splits the event log into cumulative epochs, i.e., starting from the start up to selected number of final points.
#' Returns an object which has to be fed into \code{epochs} argument of \code{toporanga} function.
#' @param n number of epochs.
#' Silently capped to the number of events, if larger.
#' Inferred from \code{window} if missing.
#' @param window size of the moving window; if \code{t} is given, in the units of time, otherwise in event count.
#' Set to \code{range/n} if not given.
#' Larger windows are just made overlapping.
#' @param t optional event times, vector of a same length and order as the event log given to \code{toporanga} function.
#' When given, epochs will be generated based on equal-time breaks not equal-event-number breaks.
#' @note If time is not given, event log must be sorted; otherwise, epochs will not make any sense.
#' @returns A special function that can be passed to the \code{epochs} argument of \code{toporanga} function.
#' @export
window_epochs<-function(n,window,t){
  if(missing(t)){
    t<-NULL
  }else{
    stopifnot(is.numeric(t))
    stopifnot(all(is.finite(t)))
    tr<-range(t)
    tm<-tr[1]
    tr<-diff(tr)
  }
  if(missing(n) && missing(window)) stop("Need at least n or window, but both are missing")
  if(missing(n)) n<-NULL
  if(missing(window)){
    window<-NULL
  }else{
    stopifnot(is.numeric(window))
    stopifnot(window>0)
    if(is.null(t) && (window<1)) window<-1
  }
  
  function(m,e){
    E<-parent.env(environment())
    if(missing(e)){
      #Set-up
      if(!is.null(t)) if(length(t)!=m) stop("Time vector doesn't have the same length as event count ",m)
      if(is.null(window)){
        if(is.null(t)){
          assign("window",pmax(m/n,1),envir=E)
        }else{
          assign("window",diff(range(t))/n,envir=E)
        }
      }
      if(is.null(n)){
        if(is.null(t)){
          assign("n",round(m/window),envir=E)
        }else{
          assign("n",round(range(t)/window),envir=E)
        }
      }
      if(is.null(t) && window<=3) warning("Very small window ",window)
      if(n>m) assign("n",m,envir=E)
      return(list(n=n,t=t))
    }else{
      if(is.null(t)){
        pmax(1,(m-1-window)/(n-1)*(e-1)+1)->si
        pmin(si+window-1,m)->se
        return(list(
          mask=si:se,
          meta=list(StartEvent=si,EndEvent=se)
        ))
      }else{
        tf<-tm+(tr-window)/(n-1)*(e-1)
        te<-tf+window
        mask<-(t>=tf) & (t<=te)
        return(list(
          mask=mask,
          meta=list(StartTime=tf,EndTime=te,Events=sum(mask))
        ))
      }
    }
  }
  
}