# ==============================================================================
# STOICHCALC - R-Routines for Solving Stoichiometric Equations
# ==============================================================================
#
# Version 0.96                 Peter Reichert, May 03, 2005 (reichert@eawag.ch)
# ============
#
# Literature: - Peter Reichert, The mathematics of stoichiometric calculations,
#               submitted to Water Research.
#             - Peter Reichert, STOICHCALC - R routines for solving 
#               stoichiometric equations, in preparation.
#             - http://www.stoichcalc.eawag.ch
#
# ==============================================================================
# Core Functions
# ==============================================================================


# Definition of functions:
# ========================

# Compose composition matrix from list of substance composition vectors:
# ----------------------------------------------------------------------

calc.comp.matrix <- function(subst.comp)
{
   elements <- character(0)
   for ( i in 1:length(subst.comp) )
   {
      elements <- c(elements,names(subst.comp[[i]]))
   }
   elements <- unique(elements)
   alpha <- matrix(data=0,nrow=length(elements),ncol=length(subst.comp))
   colnames(alpha) <- names(subst.comp)
   rownames(alpha) <- elements
   for ( i in 1:length(subst.comp) )
   {
      for ( k in 1:length(subst.comp[[i]]) )
      {
         alpha[names(subst.comp[[i]])[k],i] <- subst.comp[[i]][k]
      }
   }
   return(alpha)
}


# Calculate basis of stoichiometry space that is compatible
# with mass balances and constraints:
# ---------------------------------------------------------

calc.stoich.basis <- function(alpha,constraints=list())
{
   # transpose composition matrix:
   a <- t(alpha)
   
   # add constraints as additional columns of t(alpha):
   if ( length(constraints) > 0 )
   {
      for ( i in 1:length(constraints) )
      {
         a <- cbind(a,rep(0,nrow(a)))
         for ( j in 1:length(constraints[[i]]) )
         {
            a[names(constraints[[i]])[j],ncol(a)] <- constraints[[i]][j]
         }
      }
   }
   
   # extend t(alpha) by additional zero columns if necessary:
   if ( nrow(a) > ncol(a) )
   {
      a <- cbind(a,matrix(0,nrow=nrow(a),ncol=nrow(a)-ncol(a)))
   }
   
   # perform singluar value decomposition:
   svd.res <- svd(a)
   
   # extract basis of null space:
   ut <- t(svd.res$u)
   d  <- svd.res$d
   d.max <- max(d)
   fact  <- 1e-10
   res <- matrix(nrow=0,ncol=ncol(ut))
   colnames(res) <- colnames(alpha)
   for ( i in 1:length(d) )
   {
      if ( d[i] < fact*d.max )
      {
         res <- rbind(res,ut[i,])
      }
   }
   
   # return extracted basis:
   return(res)
}


# Calculate stoichiometry of a process from involved substances and constraints:
# ------------------------------------------------------------------------------

calc.stoich.coef <- function(alpha,name,subst,subst.norm,nu.norm=1,
                             constraints=list())
{
   substances <- unique(c(subst.norm,subst))
   alpha.reduced <- alpha[,substances]
   if ( is.vector(alpha.reduced) ) 
   {
      alpha.reduced <- t(as.matrix(alpha.reduced))
   }
   nu.basis <- calc.stoich.basis(alpha.reduced,constraints)
   if ( nrow(nu.basis) == 0 ) stop("no process possible.")
   if ( nrow(nu.basis) >  1 ) stop(paste("process not unique,",
                                   nrow(nu.basis)-1,"constraints needed."))
   nu <- matrix(data=0,nrow=1,ncol=ncol(alpha))
   colnames(nu) <- colnames(alpha)
   nu[1,colnames(nu.basis)] <- nu.basis[1,]/nu.basis[1,subst.norm]*nu.norm
   rownames(nu) <- name
   return(nu)
}