R/maxent.R

In rexpokit: R Wrappers for EXPOKIT; Other Matrix Functions

# maxent
# Copy of the FD::maxent function (to avoid the dependency on package FD)
# 
# Authored by: Bill Shipley \email{bill.shipley@usherbrooke.ca}
# http://pages.usherbrooke.ca/jshipley/recherche/
# 
# Ported to FD by Etienne Laliberte.
# Copied to rexpokit by Nick Matzke (just in order to avoid the dependency on package "FD")
# 
# 
maxent <- function (constr, states, prior, tol = 1e-07, lambda = FALSE) 
	{
	# check input
	if (is.vector(constr))
		{
		means.names <- names(constr)
		constr <- matrix(constr, 1, length(constr))
		dimnames(constr) <- list("set1", means.names)
		}

	if (is.data.frame(constr)) 
			constr <- as.matrix(constr)
	if (!is.numeric(constr)) 
			stop("constr must only contain numeric values\\n")
	if (!is.numeric(tol)) 
			stop("tol must be numeric\\n")
	if (!is.logical(lambda)) 
			stop("lambda must be logical\\n")

	if (is.vector(states))
		{
		s.names <- names(states)
		states <- matrix(states, nrow = 1, ncol = length(states))
		dimnames(states) <- list("constraint", s.names)
		}

	if (is.data.frame(states)) 
			states <- as.matrix(states)
	if (!is.numeric(states)) 
			stop("states must only contain numeric values\\n")
	if (any(states <= 0)) 
			stop("states cannot contain zero or negative values\\n")
	if (dim(states)[2] == 1 && dim(states)[1] > 1) 
			states <- t(states)

	s.names <- dimnames(states)[[2]]
	c.names <- dimnames(states)[[1]]
	set.names <- dimnames(constr)[[1]]
	n.states <- dim(states)[2]
	n.traits <- dim(states)[1]
	n.sets <- dim(constr)[1]

	if (n.traits != dim(constr)[2]) 
			stop("number of constraints in constr should be equal to number of constraints in states\\n")

	if (missing(prior))
		{
		prior <- matrix(1/n.states, n.sets, n.states)
		dimnames(prior) <- list(set.names, s.names)
		}

	if (is.vector(prior))
		{
		if (length(prior) != n.states)
			{
			stop("number of states in prior should be equal to number in states\\n")
			}
		if (n.sets == 1)
			{
			prior <- matrix(prior, 1, length(prior))
			dimnames(prior) <- list("set1", s.names)
			}
		else
			{
			prior <- matrix(rep(prior, n.sets), n.sets, length(prior), byrow = T)
			dimnames(prior) <- list(set.names, s.names)
			}
		}

	if (is.data.frame(prior)) 
			prior <- as.matrix(prior)
	if (!is.numeric(prior)) 
			stop("prior must only contain numeric values\\n")
	if (dim(prior)[2] == 1 && dim(prior)[1] > 1) 
			prior <- t(prior)
	if (dim(prior)[2] != n.states) 
			stop("number of states in prior should be equal to number in states\\n")
	if (dim(prior)[1] > 1 && dim(prior)[1] != n.sets) 
			stop("number of rows in prior should be 1 or equal to number of rows in constr\\n")
	if (dim(prior)[1] == 1)
		{
		prior <- matrix(rep(prior[1, ], n.sets), n.sets, ncol(prior), byrow = T)
		}
	if (any(prior < 0 | prior > 1))
		{
		stop("prior must contain probabilities between 0 and 1\\n")
		}
	prior <- t(apply(prior, 1, function(x) x/sum(x)))
	dimnames(prior) <- list(set.names, s.names)
	if (any(is.na(constr)) || any(is.na(states)) || any(is.na(prior)))
		{
		stop("no NA's allowed\\n")
		}
	
	
	# create storage objects
	allprobs <- matrix(NA, n.sets, n.states)
	dimnames(allprobs) <- list(set.names, s.names)
	moments <- matrix(NA, n.sets, n.traits)
	dimnames(moments) <- list(set.names, c.names)
	entropy <- rep(NA, n.sets)
	names(entropy) <- set.names
	iter <- rep(NA, n.sets)
	names(iter) <- set.names
	
	if (lambda)
		{
		lambdas <- matrix(NA, n.sets, n.traits + 2)
		dimnames(lambdas) <- list(set.names, c("intercept", c.names, "prior"))
		}
	
	
	# FORTRAN loop
	for (i in 1:n.sets)
		{
		itscale <- .Fortran("itscale5", as.double(t(states)), as.integer(n.states), as.integer(n.traits), as.double(constr[i, ]), as.double(prior[i, ]), prob = double(n.states), entropy = double(1), niter = integer(1), as.double(tol), moments = double(n.traits), PACKAGE = "rexpokit")
		allprobs[i, ] <- itscale$prob
		moments[i, ] <- itscale$moments
		entropy[i] <- itscale$entropy
		iter[i] <- itscale$niter
		if (lambda)
			{
			lambdas[i, ] <- coef(lm(log(itscale$prob) ~ t(states) + log(prior[i, ])))
			}
		}
		
	# BUILD OUTPUT
	res <- list()
	if (n.sets == 1)
		{
		res$prob <- allprobs[1, ]
		res$moments <- moments[1, ]
		names(entropy) <- NULL
		res$entropy <- entropy
		names(iter) <- NULL
		res$iter <- iter
		if (lambda) 
				res$lambda <- lambdas[1, ]
		res$constr <- constr[1, ]
		res$states <- states
		res$prior <- prior[1, ]
		} else {
		res$prob <- allprobs
		res$moments <- moments
		res$entropy <- entropy
		res$iter <- iter
		if (lambda) 
				res$lambda <- lambdas
		res$constr <- constr
		res$states <- states
		res$prior <- prior
		}
	return(res)
}