mat.index <- function(aind, dims) {
# WHO: Jim Hughes
# WHAT: Returns the vector index of the aind'th element of
#       an array with dimension dims.
# WHEN: Dec 16,1991
# WHERE: Any machine running Splus
# WHY: 
# HOW: Source this file in Splus
# HOW MUCH: Not resource intensive
# EG:   x <- mat.index(c(2,3,1),c(4,3,2))
#       x
#       [1] 10
#
#
        delta <- exp(cumsum(log(dims)))
        indices <- aind[1]
        if(length(dims) > 1)
                for(j in 2:length(dims))
                        indices <- indices + (aind[j] - 1) * delta[j - 1]
        round(indices)
}
#####

count<-function(chain, states)
{
#
# INPUT: chain - time series of events
#        states - number of states in chain
#
# OUTPUT: the transition matrix for chain
#
# This function returns
# a transition matrix for chain (missing values allowed)
#
# check the number of states
        if(max(chain[!is.na(chain)]) > states) 
                        stop(paste("you need", max(chain),
                        "states for your\ndata"))       
# initialize the counts
        kount <- matrix(0, nrow = states, ncol = states)        
# start counting
        for(i in 2:length(chain))
                if (sum(is.na(chain[(i-1):i]))==0)
                kount[chain[i - 1], chain[i]] <- kount[chain[i - 1], chain[
                        i]] + 1
        return(kount)
}
#####

qcount <- function(chain,states) {
#
# INPUT: chain - time series of events
#        states - number of states in chain
#
# OUTPUT: the transition matrix for chain
#
# this function is a quick version of count. It returns
# a transition matrix for chain if
# chain does not have missing values
#
# check the number of states
        if(max(chain[!is.na(chain)]) > states) stop(paste("you need", max(
                        chain), "states for your\ndata"))
  n <- length(chain)
  temp1 <- chain[1:n-1]
  temp2 <- chain[2:n]
  table(temp1,temp2)
}
#####

gcount <- function(chain, states, iord, maxord)
{
#
# Generalization of the count function.
# Returns a transition matrix of counts which is the
# number or transitions from state (sub iord+1) at time t-iord to
# state (sub iord) at time t-iord+1 to ... to state (sub 1) at
# time t. Optionally, one can include maxord, the maximum order you
# intend to fit. In this case, only data that would be available
# for the maximum order is used in the computations.
#
#
# check the number of states
        if(max(chain[!is.na(chain)]) > states) stop(paste("you need", max(
                        chain), "states for your\ndata"))
# if maxord is missing, set it equal to iord
        if(missing(maxord)) maxord <- iord      
# initialize the counts
        dims <- rep(states, iord + 1)
        if(length(dims) == 1)
                kount <- matrix(0, 1, states)
        else kount <- array(0, dims)
        for(i in (maxord + 1):length(chain)) {
                aind <- chain[(i - iord):i]
                if(sum(is.na(chain[(i - maxord):i])) == 0) {
                        vind <- mat.index(aind, dims)
                        kount[vind] <- kount[vind] + 1
                }
        }
        return(kount)
}
#####

ccount <- function(chain, states, cond, cstate)
{
# conditional counts;
# return a transition matrix for each value of cond
#
# check the number of states
        if(max(chain[!is.na(chain)]) > states) 
                        stop(paste("you need", max(chain),
                        "states for your\ndata"))       
# initialize the counts
        kount <- array(0, c(states, states, cstate))        
# start counting
        for(i in 2:length(chain))
                if (sum(is.na(c(chain[(i-1):i],cond[i])))==0)
                kount[chain[i - 1], chain[i], cond[i]] <- 
                     kount[chain[i - 1], chain[i], cond[i]] + 1
        return(kount)
}
#####