\documentclass[article,nojss]{jss}
\DeclareGraphicsExtensions{.pdf,.eps}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Add-on packages and fonts
\usepackage{mathptmx}
\usepackage{float}
\usepackage{graphicx}
\usepackage{amsmath}
%\usepackage{amssymb}
\usepackage{array} % tabel commands
\usepackage{xspace}
\usepackage{verbatim}
\usepackage[english]{babel}
\usepackage{times}
\usepackage{textcomp}
%\usepackage{mathptmx}
%\usepackage{helvet}
\usepackage[T1]{fontenc}
\usepackage[latin1]{inputenc}
\renewcommand{\ttdefault}{lmtt}
\renewcommand{\familydefault}{\rmdefault}
\usepackage[T1]{fontenc}
\usepackage[latin1]{inputenc}
\usepackage{geometry}


\newcommand{\noun}[1]{\textsc{#1}}
%% Bold symbol macro for standard LaTeX users
\providecommand{\boldsymbol}[1]{\mbox{\boldmath $#1$}}

%% Because html converters don't know tabularnewline
\providecommand{\tabularnewline}{\\}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% User specified LaTeX commands.
\newcommand{\li}{\textbf{\textsf{LIM }}}
\newcommand{\ls}{\textbf{\textsf{limSolve }}}
\newcommand{\footnoteremember}[2]{
  \footnote{#2}
  \newcounter{#1}
  \setcounter{#1}{\value{footnote}}
}
\newcommand{\footnoterecall}[1]{
  \footnotemark[\value{#1}]
}

\title{Package \li, implementing linear inverse models in \proglang{R} }
\Plaintitle{Package LIM, implementing linear inverse models in R}

\Keywords{Linear inverse models, flux balance analysis, linear programming,
text files, \proglang{R}}

\Plainkeywords{Linear inverse models, flux balance analysis, linear programming,
text files, R}


\author{Karline Soetaert and
Dick van oevelen\\
Centre for Estuarine and Marine Ecology\\
Netherlands Institute of Ecology\\
The Netherlands}

\Plainauthor{Karline Soetaert, Dick van oevelen}

\Abstract{We present \proglang{R} package \li \citep{LIM} which is designed
  for reading and solving linear inverse models (LIM). The model
  problem is formulated in text files in a way that is natural and
  comprehensible. \li then converts this input into the required
  linear equality and inequality conditions, which can be solved
  either by least squares or by linear programming techniques. By
  letting an algorithm formulate the mathematics, it becomes very
  simple to reformulate the model in case a parameter value changes,
  or a component is added or removed.

  Three different types of problems are supported: flow networks,
  reaction networks and other (operations research) problems. The
  first two cases are based on mass balances of the components.

  We give three examples, a food web example, a biogechemical
  reaction example and a blending example.
  
  If you use this package, please cite as: \citep{oevelen09}.

}

%% The address of (at least) one author should be given
%% in the following format:
\Address{
  Karline Soetaert\\
  Centre for Estuarine and Marine Ecology (CEME)\\
  Netherlands Institute of Ecology (NIOO)\\
  4401 NT Yerseke, Netherlands
  E-mail: \email{k.soetaert@nioo.knaw.nl}\\
  URL: \url{http://www.nioo.knaw.nl/users/ksoetaert}\\
   \\
  Dick van Oevelen\\
  Centre for Estuarine and Marine Ecology (CEME)\\
  Netherlands Institute of Ecology (NIOO)\\
  4401 NT Yerseke, Netherlands
  E-mail: \email{d.vanoevelen@nioo.knaw.nl}\\
  URL: \url{http://www.nioo.knaw.nl/users/dvanoevelen}\\
}


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% R/Sweave specific LaTeX commands.
%% need no \usepackage{Sweave}
%\VignetteIndexEntry{LIM: implementing linear inverse models in R}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Begin of the document
\begin{document}

\maketitle


\section{Introduction}
  In many disciplines, mathematical formulation of problems lead to a
  combination of linear equalities that are
  supplemented with linear inequality constraints. Such linear
  equations arise for instance:
  \begin{itemize}
    \item by considerations that certain quantities have to be
      positive, that the summed values should not exceed a certain value
      (i.e. summed fractions or probabilities should remain smaller or
      equal to 1), etc.
    \item In curve fitting problems, inequality
      constraints may arise by requirements of monoticity,
      nonnegativity, convexity, while in piecewise linear fitting,
      equality conditions result from the need to guarantee continuity
      and smoothness of the curves.
    \item In biochemical applications, the linear equalities arise because of
      linear conservation relationships such as the conservation of
      mass, charge, etc.., while inequalities ensure that mass remains a
      positive quantity.
  \end{itemize}

\section{Linear Inverse Models}

  Mathematically, linear inverse problems can be written in matrix
  notation as: \footnote{notations: vectors and matrices are in
  \textbf{bold}; scalars in normal font. Vectors are indicated with a
  small letter; matrices with capital letter. }
  \begin{align*}
    \label{eq:1}
      \mathbf{A}\cdot \mathbf{x}\simeq\mathbf{b}   \qquad  (1)\\
      \mathbf{E}\cdot \mathbf{x}=\mathbf{f}     \qquad (2)\\
      \mathbf{G}\cdot \mathbf{x}\geq\mathbf{h}    \qquad  (3)
  \end{align*}

  These are three sets of linear equations: equalities that have to be
  met as closely as possible (1), equalities that have to be met exactly
  (2) andinequalities (3).

  Often the problem originally only contains the latter two types of
  equations (2-3), and the approximate equalities are added to single
  out one solution.

  Quadratic and linear programming methods are the main mathematical
  techniques to solve for the vector x in this type of models. In
  \proglang{R}, these are made available through package \ls
  \citep{limSolve}.

  Depending on the active set of equalities (2) and constraints (3), the
  system may either be underdetermined, even determined, or
  overdetermined. Solving these problems requires different mathematical
  techniques.
  \begin {itemize}
    \item If the model is even determined, there is only one solution
      that satisfies the equations exactly. This solution can be singled
      out by matrix inversion (e.g. the \code{solve} function, in case there
      are no inequalities) or using the least squares method \emph{lsei}
      from package \ls.
    \item If the model is overdetermined, there is only one solution in the
      least squares sense; this solution is singled out by function
      \emph{lsei} (\emph{l}east \emph{s}quares with \emph{e}qualities and
      \emph{i}nequalities). This function also returns the parameter
      covariance matrix, which gives indication on the confidence interval
      and relationship among the estimated unknowns (elements in \textbf{x}).
    \item If the model is underdetermined, there exist an infinite amount of
      solutions. To solve such models, there are several options:
      \begin{itemize}
        \item \emph{ldei} - finds the "least distance" solution, i.e. the
          solution with minimal sum of squared unknowns.
        \item \emph{lsei}- minimises some other set of linear functions
          ($\mathbf{A}\cdot \mathbf{x} \simeq \mathbf{b}$) in a least
          squares sense
        \item \emph{linp} - finds the solution where \textbf{one} linear
          function (i.e. the sum of unknowns) is either minimized (a "cost"
          function) or maximized (a "profit" function)
        \item \emph{xranges} -  finds the possible ranges ([min,max])
          for each unknown.
        \item \emph{xsample} -  randomly samples the solution space using
          a Markov chain. This method returns the marginal probability
          density function for each unknown. \citep{xsample}
      \end{itemize}
  \end{itemize}
  All these functions are also available from package \li.
\section{Three types of LIM}

  One of the main remaining challenges in LIM models constitutes the
  setup of this type of problems. Especially when many unknowns have to
  be simultaneously estimated and the problem contains many equality
  and inequality constraints, the construction of the matrix equations
  may be quite complicated and error-prone.
  In addition to providing methods of solution,
  R-package \li has been designed to facilitate problem implementation.

  Depending on how the problem is formulated and which are the
  unknowns, \li distinguishes three types of Linear Inverse Models
  (Figure \ref{fig:gen}).

  \begin{itemize}
    \item flow networks. Here the problem consists of a number of
      compartments, connected by flows. Solving the model then constitutes
      of deriving the values of the flows between the compartments.
    \item reaction networks. The problem consists of a number of
      compartments that are involved in reactions. The LIM will estimate
      the reaction rates.
    \item other. \li can also solve problems often occurring in
    operational research, e.g. to find the optimal allocation of resources,
    optimal diet composition etc....
  \end{itemize}
  We give examples of these three types below.
\begin{figure}[h]
\centering
\includegraphics[width=0.8\textwidth]{FigGeneral}
\caption{Three types of Linear Inverse Models that can be created
  and solved with \proglang{R} package \li. A. Flow networks. B.
  Reaction networks, C. Other. In type (A) and (B), a mass balance
  of components is generated. This is not the case for type
  C.}\label{fig:gen}
\end{figure}

\subsection{Flow network problems}

  Flow networks are represented as a set of nodes (compartments), which are
  connected by arrows (flows).
  The arrows generally have a direction, i.e. the flows are positive.
  Thus
  \[
    A \rightarrow B
  \]
  denotes a flow directed from A to B, while
  \[A
    \leftrightarrow B
  \]
  denotes a flow that can proceed in both directions.

  There can only be one flow from A to B (but there can also be a flow
  from B to A).


  Solving the LIM-problem consists of finding the values of the flows.

  After solution, several indices and food web properties can be estimated,
  using functions from package \pkg{NetIndices} \citep{NetIndices, NetIndicesArt}

\subsubsection{Example: a simple food-web}
  Organisms eat and are eaten; they use part of their food for
  biomass production and reproduction, part is expelled as faeces or
  respired. Other (so-called autotrophic) organisms
  produce biomass from light energy and inorganic compounds, whilst
  dead matter (detritus) may be consumed by animals and bacteria.

  When the mass balances of several groups of organisms (and dead
  matter) are considered together, we obtain a food web model. In
  this type of LIM, the unknowns are the food web \emph{flows} that connect
  the components (organisms and dead matter).

  Assume a simple food web comprising a plant, detritus and an
  animal that eats both the plant and detritus.
  For simplicity we assume that the system is
  in a climax situation, i.e. the masses, which are expressed in
  $moles~C~m^{-2}$ are invariant in time. There are eight flows that
  connect the components with each other and with the outside world.
  \footnote{Since the foodweb is a subsystem of a larger system, we need
  to distinguish between model compartments, i.e. compartments whose dynamics
  are fully described in the model and external compartments, whose dynamics is
  coupled to processes occurring outside the model realm. The difference is
  essential: \li will create mass balance equations for model compartments only.
  In the example, there is no balance for $CO_2$}
  The mass balance equation for the three components and with the
  rate of change = 0, is given by:
  \begin{eqnarray*}
    \frac{\mathrm{dPLANT}}{dt} = 0 &=& \mathrm{net~primary~production}
    - \mathrm{grazing~ on~ plant} - \mathrm{plant~
    mortality}\\
    \frac{\mathrm{dANIMAL}}{dt} = 0 &=& \mathrm{grazing~on~plant} +
    \mathrm{grazing~on~detritus}- \mathrm{animal~respiration} \\&-&
    \mathrm{animal~mortality} - \mathrm{faeces~production}\\
    \frac{\mathrm{dDETRITUS}}{dt} = 0 &=& \mathrm{plant~mortality} +
    \mathrm{animal~mortality} + \mathrm{faeces~production} \\&-&
    \mathrm{grazing~on~detritus} - \mathrm{detritus~mineralisation}
  \end{eqnarray*}

These mass balances can be written in a more general way, and
using shorthand notation for the flows, as:
\begin{eqnarray}
 0 = 1 \cdot NPP - 1 \cdot Pgraz - 1 \cdot Pmort + 0 \cdot Dgraz + 0 \cdot Aresp + 0 \cdot Amort + 0 \cdot Faeces + 0 \cdot Detmin  \\
 0 = 0 \cdot NPP + 1 \cdot Pgraz + 0 \cdot Pmort + 1 \cdot Dgraz - 1 \cdot Aresp - 1 \cdot Amort - 1 \cdot Faeces + 0 \cdot Detmin  \\
 0 = 0 \cdot NPP + 0 \cdot Pgraz + 1 \cdot Pmort - 1 \cdot Dgraz + 0 \cdot Aresp + 1 \cdot Amort + 1 \cdot Faeces - 1 \cdot Detmin
 \end{eqnarray}

These equations relate, on the left hand side, the zero rates of
changes to a sum of products, where each product is composed of
the flows and a coefficient. The coefficient indicates if and how
much these flows contribute to the rate of change.

Now assume that net primary production and the total grazing rate
(Grazing) of the animal has been measured (30 $mmol~C~m^{-2}~ d^{-1}$ and
10 $mmol~C~m^{-2}~d^{-1}$ respectively). Thus, we can add two extra
equations:
\begin{eqnarray}
NPP = 30\\
Pgraz + Dgraz=10
\end{eqnarray}

In matrix notation, we obtain
\[
\left[ {\begin{array}{*{20}c}
   1 & { - 1} & { - 1} & 0 & 0 & 0 & 0 & 0  \\
   0 & 1 & 0 & 1 & { - 1} & { - 1} & { - 1} & 0  \\
   0 & 0 & 1 & { - 1} & 0 & 1 & 1 & { - 1}  \\
   1 & 0 & 0 & 0 & 0 & 0 & 0 & 0  \\
   0 & 1 & 0 & 1 & 0 & 0 & 0 & 0  \\
\end{array}} \right] \cdot \left[ {\begin{array}{*{20}c}
   {NPP} \\{Pgraz}\\{Pmort}\\{Dgraz}  \\{Aresp} \\{Amort}\\{Faeces}\\{Detmin}\\
\end{array}} \right] = \left[ {\begin{array}{*{20}c}
   0  \\   0  \\   0  \\   {30}  \\   {10}  \\
\end{array}} \right]
\]



  The feeding, defaecation and respiration flows are not
  independent of one another. Firstly, organisms cannot produce more faeces
  than the amount of food they ingest.
  Thus it is customary in foodweb modelling, to assume that
  faeces production lies in between some range of food ingested.
  For our example we assume that in between 30 and 60\% of total food ingested is
  defaecated (the food is not of high quality).

  Secondly, organisms respire carbohydrates to provide the energy for growth.
  Thus, of the fraction of the food that is assimilated (i.e. not defaecated),
  part will be used to create new biomass, the other part will provide
  the energy to do so (this is referred to as the cost of growth).
  Here we assume that 30\% of the assimilated food is respired.
  As total animal respiration also includes basal respiration (for the
  animal's maintenance), we impose that the animal respiration has to be
  larger than - or equal - to this amount:

  \begin{eqnarray}
    0.3\cdot{Pgraz} + 0.3\cdot{Dgraz} <= Faeces\\
    0.6\cdot{Pgraz} + 0.6\cdot{Dgraz} >= Faeces\\
    0.3\cdot({Pgraz} + {Dgraz} - {Faeces}) <=Aresp
  \end{eqnarray}

  Adding to that the requirement that the flows have to be positive
  gives the following set of inequality conditions:
  \[
    \left[ {
      \begin{array}{*{20}c}
         0 & -0.3 & 0 & -0.3 & 1 & 0 & 0.3 & 0  \\
         0 &  0.6 & 0 &  0.6 & 0 & 0 &  -1 & 0  \\
         0 & -0.3 & 0 & -0.3 & 0 & 0 &   1 & 0  \\
         1 & 0 & 0 & 0 & 0 & 0 & 0 & 0  \\
         0 & 1 & 0 & 0 & 0 & 0 & 0 & 0  \\
         0 & 0 & 1 & 0 & 0 & 0 & 0 & 0  \\
         0 & 0 & 0 & 1 & 0 & 0 & 0 & 0  \\
         0 & 0 & 0 & 0 & 1 & 0 & 0 & 0  \\
         0 & 0 & 0 & 0 & 0 & 1 & 0 & 0  \\
         0 & 0 & 0 & 0 & 0 & 0 & 1 & 0  \\
         0 & 0 & 0 & 0 & 0 & 0 & 0 & 1  \\
      \end{array}}
    \right]
    \cdot
    \left[ {
      \begin{array}{*{20}c}
        {NPP}\\{Pgraz}\\{Pmort}\\{Dgraz}\\{Aresp}\\{Amort}\\{Faeces}\\{Detmin }\\
      \end{array}}
    \right]
    >  =
    \left[ {
      \begin{array}{*{20}c}
        0  \\   0  \\   0  \\   0  \\   0  \\   0  \\   0  \\   0  \\
      \end{array}}
    \right]
  \]


  This model comprises 5 equations and 11 inequalities; there are
  8 unknown flows.
  We will outline below how this particular problem can be implemented
  and solved in package \li.
\subsection{reaction problems}

  These are LIM problems which are written as a set of reactions that
  connect the dynamics of several constituents. For instance, in the
  reaction
  \[
    A+B \rightarrow C
  \]
  C is produced while A and B are consumed in a stoichiometric ratio of 1 to 1.
  Some reactions can occur in two directions, e.g.
  \[
    A+B \leftrightarrow C
  \]
  In contrast to previous ("flow") network problems, where only one link
  between two compartments was allowed, in reaction problems there may exist
  many links between the constituents.
  Solving the LIM amounts to finding values for the reaction rates.

\subsubsection{The core metabolism of E.coli}
The LIM software can be used for performing flux balance analysis 
(e.g. \citep{Edwards02}).

See vignette ("LIMecoli") \citep{limecoli} for an example of how to do that.


\subsubsection{Example: chemical reactions.}
  In the natural environment, the cycles of many constituents are
  linked via chemical reactions that produce and consume them.
  We take the biogeochemical cycling of carbon (C), nitrogen (N) and
  oxygen (O) in a marine sediment as an example.

  Organic matter ($(CH_2O)_{106}(NH_3)_{16}(H_3PO_4)$) is mineralized (respired),
  using a series of oxidants: oxygen ($O_2$), nitrate ($HNO_3$) and
  some other, undefined oxidant ($XO$). The reduced byproducts of
  this mineralization process, ammonium ($NH_3$), and an undefined
  reduced substance ($X$) can be re-oxidized by a reaction with oxygen.
  All dissolved substances are exchanged with the water column.
  $N_2$, produced by the reaction of organic matter with nitrate,
  does not react in the sediment.

  The mineralisation reactions can be written as:

  \begin{align*}
    &r1 & (CH_2O)_{106}(NH_3)_{16}(H_3PO_4)+106O_2&\rightarrow106CO_2+16NH_3+H_3PO_4+106H_2O\\
    &r2 & (CH_2O)_{106}(NH_3)_{16}(H_3PO_4)+84.8HNO_3&\rightarrow106CO_2+42.4N_2+16NH_3+H_3PO_4+148.4H_2O\\
    &r3 & (CH_2O)_{106}(NH_3)_{16}(H_3PO_4)+106O_2X&\rightarrow106CO_2+106X+16NH_3+H_3PO_4+106H_2O\\
  \end{align*}
  for the oxic mineralisation, denitrification and anoxic mineralisation
  respectively.

  The secondary reactions (nitrification and reoxidation of other
  reduced substances):
  \begin{align*}
    &r4& NH_3+2O_2&\rightarrow HNO_3+H_2O \\
    &r5& X+O_2&\rightarrow O_2X \\
  \end{align*}
  and the exchange with the bottom water:
  \begin{align*}
    &r6& OMBW &\rightarrow(CH_2O)_{106}(NH_3)_{16}(H_3PO_4) \\
    &r7& O_2 &\leftrightarrow O_2BW          \\
    &r8& HNO_3 &\leftrightarrow HNO_3BW      \\
    &r9& NH_3 &\leftrightarrow NH_3BW        \\
    &r10& O_2X &\leftrightarrow O_2XBW       \\
    &r11& H_3PO_4 &\leftrightarrow H_3PO_4BW \\
    &r12& CO_2 &\leftrightarrow  CO_2BW      \\
  \end{align*}
  Note that the deposition of organic matter (r6) is directed \emph{into} the
  sediment, while the direction of the other fluxes can go either into or out of
  the sediment.

  In this LIM, the rates of the mineralisation reactions, of the
  secondary reactions and the exchange reactions with the bottom
  water are the unknowns (r1-r12). Based on these reactions, and following
  the law of conservation of mass, we can write a mass
  balance reaction for the following 7 constituents:
  $(CH_2O)_{106}(NH_3)_{16}(H_3PO_4)$,
  $O_2$, $CO_2$, $NH_3$, $H_3PO_4$, $HNO_3$, while for the others (e.g. $O_2BW$)
  only part of the reactions are specified. Hence these are considered to be
  outside the domain of the model (externals).

  We give only the mass balance reactions for $O_2$ and $HNO_3$:
  \begin{eqnarray*}
    \frac{dO_2}{dt}&=&0=-106\cdot r_1 -2 \cdot r_4 - r_5 - r_7\\
    \frac{dHNO_3}{dt}&=&0=-84.8\cdot r_1 + \cdot r_4 - r_8\\
    \cdots
  \end{eqnarray*}

  As the exchange of dissolved substance across the sediment-water interface
  can go either way, directed into or out of the sediment, they can be
  either positive or negative. Only the rates of unidirectional reactions
  need be positive, and the following inequalities hold:
  \begin{eqnarray*}
    r_1>=0\\
    r_2>=0\\
    r_3>=0\\
    r_4>=0\\
    r_5>=0\\
    r_6>=0\\
  \end{eqnarray*}

  In this particular example, the oxygen, nitrate, and ammonium fluxes
  have been estimated; they are -15 (influx), 1 (efflux) and 2
  $mmol~m^{-2}~d^{-1}$ respectively. These measurements lead to the equations:

  \begin{eqnarray*}
    r7=-15\\
    r8=1\\
    r9=2\\
  \end{eqnarray*}

  Thus there are 10 equations (7 mass balances, 3 measurements) and 12 unknowns.
  In addition, there are 6 inequality conditions\footnote{Note a difference
  with the flow networks, where the coefficients were either -1, 0, 1. Here
  the coefficients reflect the stoichiometry of the reaction and can differ
  from these numbers}.

\subsection{other problems}
  It is also possible to use \li for specifying more general (linear)
  operational research problems that do not classify as network problems.

  These problems often try to find the most efficient, or least costly,
  way of achieving something.
  They are often solved with linear programming techniques that optimize
  some function (cost or profit) given a set of linear constraints.

\subsubsection{blending problems}
  This example is borrowed from \ls and comes from the website of
  J E Beasley (find it on the web).

  A manufacturer produces a feeding mix for animals. The feed mix
  contains two nutritive ingredients and one ingredient (filler) to
  provide bulk. One kg of feed mix must contain a minimum quantity
  of each of four nutrients as below:
  \begin{table}[H]
  \center
    \begin{tabular}{l|llll}\hline
      Nutrient & A&  B&  C&  D\\
      gram  &80&  50&  25&  5\\ \hline
    \end{tabular}
  \end{table}


  The ingredients have the following nutrient values and cost:
  \begin{table}[H]
  \center
    \begin{tabular}{l|llll|l}\hline
     (gram/kg) & A  &B  &C  &D  &Cost/kg\\\hline
     Ingredient 1 & 100  &50  &40  &10  &40\\
     Ingredient 2 & 200  &150  &10  &-  &60\\
     Filler & -  &-  &-  &-  &0\\
     \hline
   \end{tabular}
  \end{table}

  The problem is to find the composition of the feeding mix that
  minimises the production costs subject to the constraints above.
  Stated otherwise: what is the optimal amount of ingredients in one
  kg of feeding mix?

  Mathematically this can be estimated by solving a linear
  programming problem, where the equalities ensure that the sum of
  the three fractions equals 1, and the inequalities enforce the
  nutritional constraints; the quantity to be minimized is the cost
  function.
  \begin{eqnarray*}
    &&\min(x_1\cdot 40+x_2 \cdot 60) \\
    &&x_i \geq 0\\
    &&x_1+x_2+x_3=1\\ and\\
    &&100\cdot x_1+200\cdot x_2\geq80\\
    &&50\cdot x_1+150\cdot x_2\geq50\\
    &&40 \cdot x_1+10 \cdot x_2\geq25\\
    &&10 \cdot x_1 \geq 5
  \end{eqnarray*}

\section{Specifying a Linear Inverse Model in R-package LIM }
  The previous examples were quite simple, and the resulting
  matrices of small or moderate size. Nevertheless, it is easy to make mistakes.
  Moreover, once the matrices are constructed, it may be quite a challenge
  to update them after adding or removing constituents.
  Also, based on the resulting set of linear equations
  it is not straightforward to infer the underlying model assumptions.

  In general, a linear inverse model is first formulated verbally, after which
  the verbal description of the problem is translated into an equivalent
  mathematical description.

  Typically the equations are specified on
  aggregated unknowns, i.e. unknowns that are themselves linear
  combinations of other unknowns. For instance, in the food web model
  example, the faeces production (the flow from the animal to detritus)
  is specified as a part of the amount of food ingested.
  Ingested food is itself the sum of the flow from the plant to
  the animal and from detritus to the animal.

  Model input in \li is close to these verbal statements. Thus to
  implement the food web model we first define a variable called
  "Ingestion" that consists of the sum of the two feeding flows and then
  define the defaecation constraints on this variable. When the LIM input
  is parsed, the constraints will be rewritten as a function of the unknowns.

  Apart from this more natural input, there are many other benefits of
  using \li.  For instance, for the \emph{flow network} and
  \emph{reaction network} type of problems, \li generates the mass balances
  for each component, based on the flows or reactions that were defined.
  This facilitates adding or removing flows or constituents.
  Finally, solving the model will also generate estimates of all
  defined variables.

  We now document the input for each of the above introduced problems.

\subsection{food web problem}
\begin{verbatim}
=====================================================
Header of the file - ignored
file: foodweb.lim

Solve the model in R with:
require(LIM)
lim <- Setup("foodweb.lim")
Ldei(lim)
Xranges(lim)
=====================================================

## EXTERNAL
  CO2
  EXP   ! export
## END EXTERNAL

## COMPONENT
  Pl    ! plant
  AN    ! Animal
  Det   ! Detritus
## END COMPONENT

## Flows
  NPP   : CO2 -> Pl
  Pgraz : Pl -> An
  Pmort : Pl -> Det
  Dgraz : Det -> An
  Aresp : An ->CO2
  Amort : An ->EXP
  Faeces: An ->Det
  Detmin: Det -> CO2
## END Flows


## PARAMETERS
  minFaeces  = 0.3
  maxFaeces  = 0.6
  growthCost = 0.3
## END PARAMETERS

## VARIABLES
  Ingestion     = Pgraz + Dgraz
  Assimilation  = Ingestion - Faeces
  GrowthResp    = Assimilation*growthCost
## END VARIABLES

## Equalities
  Faeces = 30
  Det -> CO2 = 10
## End equalities

## Inequalities
  growthcost : Aresp > GrowthResp
  defaecation: Faeces = [minFaeces,maxFaeces]*Ingestion
## End inequalities

\end{verbatim}

  Note the use of sections (\#\# SECTIONNAME ... \#\# END SECTIONNNAME) to
  declare items;
  the sections "COMPONENT" and "EXTERNAL" define the names;
  a mass balance equation is only generated for components, not for externals.
  A name is declared as "name:  ", an exclamation mark ("!") demarcates
  the start of a comment.

  Although more lengthy, this problem formulation is much more elegant,
  more flexible, less error-prone, and easier to understand than the
  resulting matrices themselves.

  Based on this input file, the matrices are generated using \li function
  \emph{Setup} and put in a list (see below).
  The resulting LIM input can then be solved with
  \emph{Lsei(lim, parsimonious=TRUE)} or with \emph{Ldei},
  which will generate the simplest -parsimonious- solution,
  with \emph{Xranges} which will estimate ranges of unknowns,
  or with \emph{Xsample} which will generate the conditional
  probability distribution of each flow.

  In the table below is what we obtained from running the following R-code:
\begin{verbatim}
require(LIM)
web.lim <- Setup("foodweb.lim")
pars <- Ldei(web.lim)
webranges<- Xranges(web.lim)

data.frame(webranges,parsimonious=pars$X)
\end{verbatim}

\begin{table}[H]
\center
\begin{tabular}{l|ll|l}\hline
& min& max& parsimonious\\ \hline
NPP   &  30&  80&    30.000000\\
Pgraz &   0&  80&    29.554950\\
Pmort &   0&  80&     0.445050\\
Dgraz &  20& 100&    20.445050\\
Aresp &   6&  70&     9.489658\\
Amort &   0&  49&    10.510342\\
Faeces&  30&  30&    30.000000\\
Detmin&  10&  10&    10.000000\\
\hline
\end{tabular}
\end{table}

  Based on these results it is simple to create a plot which depicts the
  parsimonious solution and the ranges (see Figure \ref{fig:range}):
\begin{verbatim}
xlim <- range(webranges)
dotchart(x=pars$X,labels=rownames(webranges),xlim=xlim,
         main="Food web",pch=16)
cc <- 1:nrow(webranges)
segments(x0=webranges[,1],y0=cc,x1=webranges[,2],y1=cc)
\end{verbatim}


\begin{figure}[H]
\centering
\includegraphics[width=0.3\textwidth]{FigWeb}
\caption{Ranges and parsimonious solution of foodweb example -
seee text for \proglang{R}-code}\label{fig:range}
\end{figure}

\subsection{chemical reaction problem}
  The input of the chemical reaction problem is:
\begin{verbatim}
=====================================================
Header of the file - ignored
file reaction.lim
0-dimensional sediment coupled C, N, O, P model
run with:
require(LIM)
reaction.lim <- Setup("reaction.lim")
X  <- Ldei(reaction.lim)
xr <- Xranges(reaction.lim, ispos = FALSE)
=====================================================
### COMPONENTS
OM
O2
CO2
NH3
H3PO4
HNO3
X
### END COMPONENTS

### EXTERNAL
H2O
N2
O2X
OMBW
O2BW
HNO3BW
NH3BW
XBW
H3PO4BW
CO2BW
### END EXTERNAL

### REACTIONS
r1: OM + 106*O2 -> 106*CO2 + 16*NH3 + H3PO4 + 106*H2O
r2: OM + 84.8*HNO3 -> 106*CO2 + 42*N2+16*NH3 + H3PO4 + 148.4*H2O
r3: OM + 106*O2X -> 106*CO2 + 106*X + 16*NH3 + H3PO4 + 106*H2O

r4: NH3+2*O2 ->HNO3 +H2O
r5: X+ O2 ->O2X

r6: OMBW -> OM
r7: O2 <-> O2BW
r8: HNO3 <-> HNO3BW
r9: NH3 <-> NH3BW
r10: X <-> XBW
r11: H3PO4 <-> H3PO4BW
r12: CO2 <-> CO2BW
### END REACTIONS

### EQUATIONS
r7 = -15
r8 = 1
r9 = 2
### END EQUATIONS

### INEQUALITY
 r1>0
 r2>0
 r3>0
 r4>0
 r5>0
 r6>0
### END INEQUALITY


\end{verbatim}

Results are in the following table:

\begin{table}[H]
\center
\begin{tabular}{l|ll|l}\hline
& min& max& parsimonious\\ \hline
r1 &   0.00000000&   0.12264151&   0.1226364\\
r2 &   0.00000000&   0.07665094&   0.00000000\\
r3 &   0.06485849&   0.51709906&   0.06486363\\
r4 &   1.00000000&   7.50000000&   1.00000000\\
r5 &   0.00000000&  13.00000000&   0.00005450\\
r6 &   0.18750000&   0.59375000&   0.18750000\\
r7 & -15.00000000& -15.00000000& -15.00000000\\
r8 &   1.00000000&   1.00000000&   1.00000000\\
r9 &   2.00000000&   2.00000000&   2.00000000\\
r10&   6.87500000&  54.81250000&   6.87500000\\
r11&   0.18750000&   0.59375000&   0.18750000\\
r12&  19.87500000&  62.93750000&  19.87500000\\
\hline
\end{tabular}
\end{table}

The marginal probability distribution of all reaction rates can
be generated by \code{Xsample} and then simply plotted using
\proglang{R}-function \code{pairs}. This is done in the
\proglang{R}-script below. Before creating the pairs plot, we
first remove the rates that were given a fixed value. On the
diagonal of the pairs plot, we plot a histogram; we define this
function first (it is copied from one of the examples in the pairs
help file). We plot only the lower part of the pairs plot (i.e.
set \code{upper.panel = NULL}).
\begin{verbatim}
xs <- Xsample(reaction.lim, jmp = 10, iter = 500)
panel.hist <- function(x, ...)  {
    usr <- par("usr")
    on.exit(par(usr))
    par(usr = c(usr[1:2], 0, 1.5) )  #redefine y-axis; x-axis stays the same
    h <- hist(x, plot = FALSE)
    breaks <- h$breaks; nB <- length(breaks)
    y <- h$counts
    y <- y/max(y)
    rect(breaks[-nB], 0, breaks[-1], y, col = "grey", ...)
}
xs <- xs[ ,-(7:9)]                    #remove constant rates
pairs(xs, upper.panel = NULL, diag.panel = panel.hist,
      pch = ".", cex = 2, main = "Reaction network")
\end{verbatim}

  The results are in Figure \ref{fig:mcmc}

\begin{figure}[h]
\centering
\includegraphics[width=0.8\textwidth]{FigReaction.pdf}
\caption{MCMC sample of the reaction network}\label{fig:mcmc}
\end{figure}

\subsection{blending problems}
  Finally we give the input for the blending problem.
\begin{verbatim}
=====================================================
Header of the file - ignored
Blending problem file blending.lim
run with:
require(LIM)
blend.lim <- Setup("blending.lim")
lp <- Linp(blend.lim)
xr <- Xranges(blend.lim, ispos = TRUE)
xs <- Xsample(blend.lim)
=====================================================

### COMPONENT
 X1   ! Part ingredient 1
 X2   ! Part ingredient 2
 X3   ! Part ingredient 3 = filler
### END COMPONENTS

## PARAMETERS
! Minimal nutrient requirements in feeding mix
 NutA = 80
 NutB = 50
 NutC = 25
 NutD = 5

! Cost of ingredients 1,2,3
 Cost1 = 40
 Cost2 = 60
 Cost3 =  0

! nutrient contents in X1 and X2
 N1_A =100
 N1_B =50
 N1_C =40
 N1_D =10

 N2_A =200
 N2_B =150
 N2_C =10
 N2_D =0

 N3_A =0
 N3_B =0
 N3_C =0
 N3_D =0
## END PARAMETERS

## COST
 Cost1*X1 + Cost2*X2 + Cost3*X3
## END COST

### EQUATIONS
 X1 + X2 + X3 = 1
### END EQUATIONS

### INEQUALITY
 X1>0
 X2>0
 X3>0

 N1_A*X1 + N2_A*X2 + N3_A*X3 >NutA
 N1_B*X1 + N2_B*X2 + N3_B*X3 >NutB
 N1_C*X1 + N2_C*X2 + N3_C*X3 >NutC
 N1_D*X1 + N2_D*X2 + N3_D*X3 >NutD

### END INEQUALITY
\end{verbatim}

  The following code generates multiple solutions (small dots) and plots these
  together with the minimal cost solution (large red dots) (see figure
  \ref{fig:blending}).
  Note that the MCMC-generated matrix is extended first with the
  parsimonious results.
\begin{verbatim}
XS <- rbind(lp$X,xs)

xsplot <- function (x,y,...) {
points(x,y,pch=".")
points(x[1],y[1],pch=16,cex=2,col="red")
}

pairs(XS,upper.panel=NULL,lower.panel=xsplot,main="blending")
\end{verbatim}

\begin{figure}[H]
\centering
\includegraphics[width=0.8\textwidth]{Figblending.pdf}
\caption{MCMC sample of the blending problem}\label{fig:blending}
\end{figure}

\subsection{A simple linear programming problem}
  It is also possible to write simple linear programming problems as in the
  following example, from Vanderbei (2006), example 2.11:

\begin{verbatim}
======================================================
File linprog.lim
Simple linear programming example
======================================================

## EQUALITIES
-x12             + x23 + x24       = 0
     - x13       - x23       + x34 = 0
## END EQUALITIES

## INEQUALITIES
 x12 + x13 + x14                   > 1
             x14       + x24 + x34 < 1
## END INEQUALITIES

## MINIMISATION
x12 + 8*x13 + 9*x14 + 2*x23 + 7*x24 + 3*x34
## END MINIMISATION
\end{verbatim}

  In this simple example, the components ("decision variables" in linear
  programming jargon) are not explicitly declared. Their names are inferred
  from the minimization function.

  This model is solved as follows:
\begin{verbatim}
Linp("linprog.lim")
\end{verbatim}

which gives:
\begin{verbatim}
$residualNorm
[1] 1.110223e-16

$solutionNorm
[1] 6

$X
      X12 X13 X14 X23 X24 X34
cost1   1   0   0   1   0   1
\end{verbatim}

\section{Structure of the LIM input file}
  Based on the examples above, we now discuss the general structure of the
  LIM input files.

  The structure of the LIM declaration file has to obey the following
  \textbf{rules}:

  \begin{itemize}
    \item Declarations are case-INsensitive: flows, Flows, FLOWS is
        all the same.
    \item The declaration file is divided into several sections, each
       contained between '\#\# section name' and '\#\# END section name'.
       Only the text embraced by "\#\#" and "\#\# END" couples is considered
       by the \li parser. The number of \#s does not matter. Only the first
       four characters of the section names are considered, e.g. to designate
       the parameter section, we can write \#\# PARAM or \#\#PARAMETERS.
       The declaration sections allowed are summarised in table \ref{tb:dec}.
    \item text inbetween the declaration sections is ignored (and can be used
      to write comments). In the foodweb example for instance, all text
      positioned in front of "\#\# EXTERNAL" will be ignored.
    \item An input file can contain declarations for externals, components,
      flows, parameters, variables and also defines the additional equalities
      (i.e. not the mass balances) and inequalities, costs, and profits
      (see below).
    \item Any line that starts with "!" or any blank line is ignored.
      The exclamation mark can also be used to discard part of an input line
      (i.e. everything past "!" is ignored).
    \item Simple calculations are allowed, i.e. addition and multiplication.
      The use of brackets for a calculation is not allowed.
    \item Continuation of a line is allowed via the use of the "\&" sign,
      at the end of the line.
    \item flows can also be given a name  using ("name :"). Although this is
      not mandatory, it may make the equations more readable.
    \item equalities and inequalities can also be given a name. This is only
      used for output.
  \end {itemize}

  A number of \textbf{shorthand} notations are available:
  \begin{itemize}
    \item If the LIM is a flow network, then \emph{FLOWfrom(x)} is shorthand
      for the sum of all flows directed out of component x, while
     \emph{FLOWto(x)} is shorthand for all flows directed into component x.
       In the foodweb model example, we wrote:
\begin{verbatim}
## VARIABLES
  Ingestion     = Pgraz + Dgraz
## VARIABLES
\end{verbatim}
This could have been written as:
\begin{verbatim}
## VARIABLES
  Ingestion     = Flowto(An)
## VARIABLES
\end{verbatim}
    \item In the inequality section, using \begin{verbatim}[]\end{verbatim}
      assigns in one statement lower and upper bounds. In the foodweb
      example for instance, we wrote:
\begin{verbatim}
  defaecation: Faeces = [minFaeces,maxFaeces]*Ingestion
\end{verbatim} which is the same as:
\begin{verbatim}
  defaecation1: Faeces > minFaeces*Ingestion
  defaecation2: Faeces < maxFaeces*Ingestion
\end{verbatim}
\end{itemize}
\clearpage


\begin{table}[H]
\caption{Nomenclature for LIM elements; sometimes several names are allowed for one
  type of element; the parser only considers the first four characters.}\label{tb:dec}
\centering
\begin{tabular}{p{.28\textwidth}p{.72\textwidth}}\\
  Name & Description\\
\hline
  COMP, STATE, STOCKS, DECI, UNKN   &  compartments (or components) between
    which flows are defined.
    If neither \emph{FLOWS} nor \emph{REACTIONS} are defined,
    then the compartments constitute the unknowns to be estimated (so-called
    "decision variables", see the blending example).
    If \emph{FLOWS} or \emph{REACTIONS}
    are present, then there will be one mass balance generated for each
    compartment; this distinguishes them from \emph{EXTERNALS}. If the
    components are not explicitly specified, they will be generated from the
    \emph{FLOWS}, \emph{REACTIONS} or \emph{COST} or \emph{PROFIT} section.
    To avoid errors due to typing mistakes, it is recommended to explicitly
    define components; in this case the parser can check if all items used
    in the flow or reaction section actually exist.
    Compartments may be given a value.\\
  EXTERNALS       &  compartments that represent the external world. There
    is no mass balance generated for \emph{EXTERNALS}\\
  PARAMETERS      &  Parameters (with their values) that have constant
    values during one model application but whose value can be altered for
    other applications. They can be changed in monte carlo runs (see below).
    Parameters may be calculated based on other parameters (that have been
    declared in front of the derived parameter).\\
  FLOWS           &  Flows between two components, written either as
    Flow(Source,Sink) or Source $\rightarrow$ Sink, and where source and
    sink are components.
    If this section is defined, then the \emph{FLOWS} will be considered
    the unknowns that have to be estimated.
    This declaration section forms the basis of a set of mass balance
    equations, one for each component.
    Cannot co-occur together with \emph{REACTIONS} declarations\\
  REACTION        &  A reaction occurring between components, e.g.
    A $\rightarrow$ B + 2*C. Reactions that can occur in two directions are
    denoted as in A $\leftrightarrow$ B+C.
    For each unidirectional reaction, the rates are positive; thus an
    inequality condition (rate>0) is imposed. Reactions occurring in two
    directions need not be positive.
    If this section is defined, then the \emph{reaction rates} will be
    considered the unknowns that have to be estimated.
    The reactions form the basis of a set of mass balance equations,
    one for each component.
    A \emph{REACTION} section can not co-occur with a \emph{FLOWS}
    declaration section.\\
  VARIABLES       & A linear expression involving the unknowns, parameters
    or other variables. Variables are derived quantities, useful to make
    the declaration of e.g. inequality constraints easier and more readable.
    Their values are estimated by function \code{Variables}; the ranges of
    variables can be solved with function \code{Varranges}\\
  COST, MINIMUM   &  One or more linear expression(s) that should be
    minimized \\
  PROFIT, MAXIMUM,&  One or more linear expression(s) that should be
    maximized \\
  RATES           &  Only present when the problem is a flow or reaction
    network and the net rate of change is not 0. It should give the net rate
    of change for a compartment. If not specified, steady-state is assumed and
    \emph{RATES} gets the value 0. \\
  EQUALITIES      &  Relationships between unknowns, or measured values that
    are assumed to be exactly known  \\
  INEQUALITIES    &  Relationships between unknowns, or measured values that
    are assumed to be known only with certain bounds. \\

  \hline
  \end{tabular}
  \end{table}

\section{setting up the linear inverse model}
  The LIM input is used to create the matrices and vectors that constitute
  the lsei problem (least squares with equalities and inequalities).

  This is done in two steps.
  \begin{itemize}
    \item Function \emph{Read} performs the first step, which creates the
      \emph{liminput}, a list that defines all elements of the LIM as a
      function of the other elements.
    \item Function \emph{Setup.lim} or \emph{Setup} performs the second
      step. Based on the \emph{liminput}, all terms are written as a function
      of the unknowns only and the matrices and vectors \textbf{A},
      \textbf{b}, \textbf{G} and \textbf{h} are created.
  \end{itemize}

  It is also possible to create the LIM matrices and vectors directly from
  an input file. This is done by calling function \emph{Setup} which takes \
  as input a file name.

  Thus:
\begin{verbatim}
lim <- Setup("linprog.lim")
\end{verbatim}
does the same as:
\begin{verbatim}
liminput <- Read("linprog.lim")
lim <- Setup(liminput)
\end{verbatim}

  Splitting problem generation in two steps is convenient when several
  runs need to be performed with different parameter values, e.g. for
  performing a monte carlo analysis. Thus, the values of the parameters
  can be directly altered in the \emph{liminput} list, after which
  \emph{Setup} will recreate the corresponding matrices and vectors.

  In the next sections we take a closer look at how setting up the LIM
  is achieved. This rather technical information can be skipped.
\subsection{creating a liminput list}
  Here the names of all section elements, and their calculations are saved
  as a list of type "liminput". All elements are considered to result from
  linear calculations, which are saved as a data frame. In this data frame,
  one line denotes a product, while subsequent lines belonging to the same
  calculation are sums. A product can be composed of the following items
  (columns): constant values (column "val"), parameters (up to 4,
  columns "par1",..."par4"), variables ("var"), flows ("flow"),
  components ("comp"), externals and reactions.

  Except for the constants, all items are denoted with their number.

  Consider the following part of the liminput generated by reading the
  simple linear programming example.

  The data frame captures the calculation of the inequalities. They were
  defined as:
\begin{verbatim}
## INEQUALITIES
 x12 + x13 + x14                   > 1
             x14       + x24 + x34 < 1
## END INEQUALITIES
\end{verbatim}
  and are parsed into the following data.frame:
\begin{verbatim}
$constraints
   name nr val par1 par2 par3 par4 var flow  comp external reaction
1 ineq1  1   1   NA   NA   NA   NA  NA   NA     1       NA       NA
2 ineq1  1   1   NA   NA   NA   NA  NA   NA     2       NA       NA
3 ineq1  1   1   NA   NA   NA   NA  NA   NA     3       NA       NA
4 ineq1  1  -1   NA   NA   NA   NA  NA   NA    NA       NA       NA
5 ineq2  2  -1   NA   NA   NA   NA  NA   NA     3       NA       NA
6 ineq2  2  -1   NA   NA   NA   NA  NA   NA     5       NA       NA
7 ineq2  2  -1   NA   NA   NA   NA  NA   NA     6       NA       NA
8 ineq2  2   1   NA   NA   NA   NA  NA   NA    NA       NA       NA
\end{verbatim}
  which should be understood as follows:

  First of all, there are two inequalities, numbered 1 and 2 (column "nr");
  both are the sum of 4 terms (there are 4 lines in the data frame for each
  ineqauality). On line 2, only column "val" and "comp" are not a NA. This
  term should be read as 1*comp[2].

  Inequality 1 can thus be reconstructed as:
  $1 * comp[1] + 1 * comp[2] + 1 * comp[3] - 1 > 0$


  The liminput data.frame contains the following elements:
  \begin{itemize}
    \item "file" - the name of the input file
    \item "pars" - a data.frame with the parameters
    \item "comp" - a data.frame with the components (state variables)
    \item "rate" - a data.frame with rates of change
    \item "extern" - a data.frame with externals
    \item "flows" - a data.frame with flows  (these have a different
      -simplified structure)
    \item "vars" - a data.frame with variables
    \item "cost" - a data.frame with cost function
    \item "profit" - a data.frame with profit function
    \item "equations" - a data.frame with equality conditions
    \item "constraints" - a data.frame with inequality conditions
    \item "reactions" - a data.frame with reaction
    \item "posreac" - a vector of logical elements: TRUE if corresponding
      reaction is positive (i.e. unidirectional reaction, $\rightarrow$) ,
      FALSE otherwise (i.e. two-ways reaction, $\leftrightarrow$)
    \item "marker" - a data.frame with  markers (see below)
    \item "parnames" - a vector with parameter names
    \item "varnames" - a vector with variable names
    \item "compnames" - a vector with component (state variable) names
    \item "externnames" - a vector with names of externals
    \item "Type" - one of "flow", "reaction",or "simple"
  \end{itemize}

\subsection{creating LIM matrices and vectors}

  Based on the \emph{liminput} list, function \emph{Setup} rewrites all
  terms as a function of the unknowns only.
  It creates an instance of class \emph{lim}, a list that contains, amongst
  other things the matrices and vectors \textbf{A}, \textbf{b}, \textbf{G}
  and \textbf{h}.

  The following elements are in type \emph{lim}
  \begin{itemize}
    \item "file" - The name of the input file
    \item "NUnknowns" - the number of unknowns
    \item "NEquations" - the number of equations inputted (the "true"
      equality conditions, i.e. excluding the mass balances for flow and
       reaction networks)
    \item "NConstraints" - the number of constraints inputted (the "true"
      inequality conditions, i.e. excluding the positivity constraints which
      are assumed for flow networks)
    \item "NComponents" - the number of components or state variables
    \item "NExternal" - the number of externals
    \item "NVariables" - the number of variables
    \item "A" - the matrix \textbf{A}, containing the coefficients of the
      equalities. If the problem is a flow or reaction network, then the
      first \emph{NComponent} equations are the mass balances, the last
      \emph{NEquations} rows correspond to the inputted equalities.
    \item "B" - the vector \textbf{b}, containing the right hand side of
      the equalities
    \item "G" - the matrix \textbf{G}, containing the coefficients of the
      inequalities. If the problem is a flow or reaction network, then the
      first \emph{NConstraints} rows correspond to the inputted inequalities,
      while the last rows correspond to the imposed positivity constraints.
      For flow networks, there are \emph{NComponents} positivity constraints ;
      for reaction networks the number of positivity constraints are less
      than or equal to this amount.
    \item "H" - the vector \textbf{h}, containing the right hand side of
      the inequalities
    \item "Cost" - the cost vector, contains the coefficients of the
      cost function
    \item "Profit" - the profit vector, contains the coefficients of the
      profit function
    \item "Flowmatrix" - if a flow network: a matrix whose i,jth value
      denotes the flow number from i to j
    \item "VarA" - variable matrix; contains the coefficients of the variables.
       \footnote{Variables are linear functions of the unknowns defined as:
       VarX*x-VarB}
    \item "VarB" - variable vector; contains the right hand side of the
      variable declarations
    \item "Parameters" - names and values of all parameters
    \item "Components" - names and values of all components
    \item "Externals" - names and values of all externals
    \item "rates"  - names and values of all rates
    \item "markers" - names and values of all markers (see below)
    \item "Variables"  - names of the variables
    \item "eqnames" - names of the equalities
    \item "ineqnames" - names of the inequalities
    \item "costnames" - names of the cost functions
    \item "profitnames" - names of the profit functions
    \item "Unknowns" - names of the unknowns
    \item "ispos" - true if the unknowns have to be positive, i.e. for
      certain flows in a flow network
  \end{itemize}

\subsection{Solving LIM problems}

  During an inverse solution, two norms are calculated:
  \begin{itemize}
    \item the residual to the equations E * X = F. This is called the
       residual norm.
       A model that can be solved has a residual norm $\sim 0$.
    \item the value of the function that has been minimised or maximised:
       MIN(f(Flows)) or MAX (f(Flows). This is the solution norm.
  \end{itemize}

\section{Examples}
  Package \li contains many example input files. They are present in three
  different subdirectories of the packages \url{examples} subdirectory.

  \begin{itemize}
    \item directory \url{Foodweb} contains the following food web examples:
    \begin{itemize}
      \item \code{BrouageMudflat.input} \citep{Leguerrier03}
      \item \code{CaliforniaSediment.input} \citep{Eldridge93}
      \item \code{Everglades.input} \citep{Diffendorfer01}
      \item \code{foodweb.lim}, the simple food web in this vignette
      \item \code{RigaAutumn.input},\code{RigaSummer.input},
        \code{RigaSpring.input} \citep{donali99}
      \item \code{ScheldtIntertidal.input} \citep{Oevelen06}
      \item \code{Takapoto.input} \citep{Niquil98}
    \end{itemize}
    \item directory \url{Reactions} contains
    \begin{itemize}
      \item \code{E_coli.lim} \citep{Edwards02}
      \item \code{fba_simple.lim}
      \item \code{reaction.lim} the example from this vignette
    \end{itemize}
    \item directory \url{LinearProg} contains
    \begin{itemize}
      \item \code{alloymixture.input}
      \item \code{blending.input}
      \item \code{diet.input}
      \item \code{FoodManufacture.001}
      \item \code{Greenberg1.input}
      \item \code{linprog.lim} the example from this vignette
      \item \code{Machines.input}
      \item \code{manpower.001}
      \item \code{refinery.001}
      \item \code{transportation.input}
      \item \code{Vanderbei1.input}, \code{Vanderbei2.input},
        \code{Vanderbei3.input}, \code{Vanderbei4.input},
        \code{Vanderbei5.input}
        5 examples from the book of \citep{Vanderbei00}
    \end{itemize}
  \end{itemize}
  Additional foodweb examples can be found in package \pkg{ToxLim} \citep{ToxLim}
\section{Extensions}
\subsection{Markers}
  For some flow network applications it may be useful to solve not
  only the default mass balance that is setup based on the section "\#\#
  COMP", but also an additional mass balances based e.g. on stable
  isotope data. Stable isotope data are collected in many food web
  studies to decipher food sources or trophic position of a consumer.
  Mathematically, isotope data are evaluated with so-called linear
  mixing models. Verbally, a linear mixing model assumes that the
  isotope value of a consumer is a flow-weighted average of its
  resources. Mathematically this is represented by:

  \begin{equation}
    X_j = \frac {\sum_i F_{i\rightarrow j} \left( X_i + \Delta_j\right)}
              {\sum_i F_{i\rightarrow j}}
    \label{eqs:IsoMod}
  \end{equation}

  in which $X_j$ is the isotope composition of the consumer,
  $F_{i\rightarrow j}$ is the flow from resource $i$ to consumer $j$
  and $\Delta_j$ is trophic fractionation (e.g. $^{15}$N fractionates
  with $\sim$3 \textperthousand per trophic level).

  Stable isotope data can be inputted as data equations in a \li.
  Suppose that in the food web example we have $\delta^{15}$N isotope
  data for plants (Pl\_iso), detritus (Det\_iso) and the consumer AN
  (AN\_iso). Further, we assume a common fractionation factor of
  3.4\textperthousand for $\delta^{15}$ (Frac\_iso). The resulting
  isotope equation can be included in the \#\# EQUALITIES section as:

\begin{verbatim}
## EQUALITIES
 ...
 AN_iso*Ingestion = Pl_iso*Pl->AN + Frac_iso*Pl->AN +
Det_iso*Det->AN + Frac_iso*Det->AN
 ...
## END EQUALITIES
\end{verbatim}

  This notation works fine and is easily manageable for this simple
  example. When however, food webs interactions are more complex, the
  \li package offers a short-hand notation for the above. We can
  rewrite equation \ref{eqs:IsoMod} as

  \begin{equation}
    \sum_i F_{i\rightarrow j} \left(X_i - X_j\right) = -\Delta_j \sum_i
    F_{i\rightarrow j}
  \end{equation}

  in which the left-hand side of the equation is the mass balance of
  the isotope and the right-hand side accounts for trophic
  fractionation processes.

  In the file we first include the parameter value for Frac\_iso in
  the \#\# PARAMETER section and define an additional section "\#\#
  MARKER" in which the marker values (in this case $\delta^{15}$N
  data) are assigned to the components of the food web.

\begin{verbatim}
## PARAMETERS
 ...
Frac_iso = 3.4
 ...
## END PARAMETERS
\end{verbatim}

\begin{verbatim}
## MARKER
 ...
 Pl  =   10
 Det =   5
 AN  =   11
 ...
## END MARKER
\end{verbatim}

  Subsequently, we define the additional mass balance for the isotope
  in the equalities.

\begin{verbatim}
## EQUALITIES
 ...
 Massbalance(AN) = -flowto(AN)*Frac_iso
 ...
## END EQUALITIES
\end{verbatim}

  In the above formulation, the left-hand side of the equality gives
  the mass balance notation as derived above and the right-hand side
  takes care of the trophic fractionation. This short-hand notation
  has several advantages. First of all, the notation in the equalities
  becomes must shorter (especially in large flow networks). Secondly,
  when one decides to remove a particular food source from a consumer
  in the \#\# FLOWS section, this change in the food web structure is
  also propagated to the isotope balance, reducing the change of
  inconsistencies between different mass balances. Currently, it is
  only possible to define one additional marker in the \li package
  using the short-hand notation.

\subsection{Monte Carlo runs}
  Sometimes, a certain model needs to be solved with different parameter
  values.

  Consider the following flow network:
\begin{verbatim}


================================================================================
file simple.input
================================================================================
##Parameter
 f1 = 1
## end parameter

## FLOWS
A             ->   B
A             ->   C
B             ->   C
C             ->   A
B             ->   D
D             ->   A

## END FLOWS

## EQUATIONS
0.5* A->B = A->C

C->A = D->A
C->A = f1
## END EQUATIONS

\end{verbatim}

  where there are four compartments, and 6 flows. The value of the flow
  from C to A is specified with a parameter (f1).

  To solve this model, we write:
\begin{verbatim}
Ldei("simple.input")
\end{verbatim}
  which outputs:
\begin{verbatim}
$X
     A->B      A->C      B->C      C->A      B->D      D->A
1.3333333 0.6666667 0.3333333 1.0000000 1.0000000 1.0000000

$unconstrained.solution
[1] 1.3333333 0.6666667 0.3333333 1.0000000 1.0000000 1.0000000

$residualNorm
[1] 1.887379e-15

$solutionNorm
[1] 5.333333

$IsError
[1] FALSE

$type
[1] "ldei"
\end{verbatim}

  Now we want solve the model successively for increasing values of the flow f1.

  Here is how to do this: we first create the \code{liminput} structure by just
  reading the input file (\code{Read("simple.input")}.

  This structure has a list item called "pars" which looks like:
\begin{verbatim}
$pars
     name nr val par1 par2 par3 par4 var flow comp external reaction
left   f1  1   1   NA   NA   NA   NA  NA   NA   NA       NA       NA
\end{verbatim}
  (it specifies the name \code{f1} and the value (\code{val}) = 1.

  We then loop over all required values of f1 (\code{for ( pars in
  seq(0,1,by=0.2))}), each time setting the value of the parameter in
  the liminput structure (\code{ls$pars$val <- pars}) and then solving the
  model with this input.
\begin{verbatim}
ls <- Read("simple.input")

Res <- NULL
for ( pars in seq(0,1,by=0.2))
{
 ls$pars$val <- pars
 Res<- rbind(Res,c(pars,Ldei(Setup(ls))$X))
}
Res
\end{verbatim}

  The results look like:
\begin{verbatim}
              A->B      A->C       B->C C->A B->D D->A
[1,] 0.0 0.0000000 0.0000000 0.00000000  0.0  0.0  0.0
[2,] 0.2 0.2666667 0.1333333 0.06666667  0.2  0.2  0.2
[3,] 0.4 0.5333333 0.2666667 0.13333333  0.4  0.4  0.4
[4,] 0.6 0.8000000 0.4000000 0.20000000  0.6  0.6  0.6
[5,] 0.8 1.0666667 0.5333333 0.26666667  0.8  0.8  0.8
[6,] 1.0 1.3333333 0.6666667 0.33333333  1.0  1.0  1.0
\end{verbatim}

  \footnote{For the die-hard who has actually reached this part of the vignette.
  R makes a vignette only from files that have extension "rnw" and that are then
  processed by R-function Sweave \citep{Leisch02}. Sweave interprets the
  R-code and makes a tex file.
  When I started writing, I did not yet know how to Sweave.
  By renaming the "tex" file as "rnw", R is tricked to believing it
  is a true Sweave file (which it is not) and thus makes a vignette.
  This means that the "R-code" that you can read is not interpreted}

  A list of all functions in \li is in table (1).

  A list of useful functions in other packages is table (2)
\clearpage
\begin{table*}[t]
\caption{Table 1. Summary of the functions in package \ls}\label{tb:ls}
\centering
\begin{tabular}{p{.15\textwidth}p{.75\textwidth}}\\
 Function & Description\\
\hline
\code{Flowmatrix}    &  Generates a flow matrix from a LIM problem\\
\code{Plotranges}    &  Plots minimum and maximum (ranges) and a central value of a LIM problem\\
\code{PrintMat}      &  Print the matrices of a LIM problem\\
\code{Read}          &  Reads a LIM input file and creates a liminput list\\
\code{Setup}         &  Composes a LIM problem from either a liminput list or from a file\\
\code{Ldei}          &  Solves a LIM problem using ldei (Least distance programming with equalities and inequalities)\\
\code{Linp}          &  Solves a LIM problem using Linear programming \\
\code{Lsei}          &  Solves a LIM problem using lsei (Least squares with equality and inequality conditions)\\
\code{Xranges}       &  Calculates ranges of unknowns\\
\code{Varranges}     &  Calculates ranges of variables (linear combinations of unknonws)\\
\code{Xsample}       &  Randomly samples a LIM problem for the unknowns\\
\code{Varsample}     &  Randomly samples a LIM problem for the inverse variables\\
\hline
\end{tabular}
\end{table*}
 \clearpage
\begin{table*}[t]
\caption{Table 2. Useful functions from other packages: 
\pkg{diagram} \citep{diagram}, \pkg{ToxLim} \citep{ToxLim} 
and \pkg{NetIndices} \citep{NetIndices, NetIndicesArt}}
\centering
\begin{tabular}{p{.15\textwidth}p{.15\textwidth}p{.65\textwidth}}\hline
\rule[-3mm]{0mm}{8mm}       Function  & Package        &Description\\
\hline \hline
\code{plotweb}        & \code{diagram}   & Plots a web, based on a flowmatrix \footnoteremember{flowmatrix}{created using \li function \code{Flowmatrix}}, thickness of arrow =value of flow  \\
\code{plotmat}        & \code{diagram}   & Visualises the transpose of a flowmatrix \footnoterecall{flowmatrix} as labeled boxes connected by arrows (created using \li function \code{Flowmatrix}  \\ \hline
\code{LimOmega}       & \code{ToxLim}    & Predicts internal concentrations of hydrophobic chemicals in aquatic organisms, based on a LIM  \\
\code{AscInd}         & \code{NetIndices} & Based on a flowmatrix\footnoterecall{flowmatrix}, estimates the ascendency network indices\\
\code{Dependency}     & \code{NetIndices} & Based on a flowmatrix\footnoterecall{flowmatrix}, estimates the dependency network indices\\
\code{EffInd}         & \code{NetIndices} & Based on a flowmatrix\footnoterecall{flowmatrix}, estimates the effective connectivity, flows, nodes, roles network indices\\
\code{EnvInd}         & \code{NetIndices} & Based on a flowmatrix\footnoterecall{flowmatrix}, estimates the indices of homogenization, synergism,...\\
\code{GenInd}         & \code{NetIndices} & Based on a flowmatrix\footnoterecall{flowmatrix}, estimates the general network indices\\
\code{PathInd}        & \code{NetIndices} & Based on a flowmatrix\footnoterecall{flowmatrix}, estimates the direct and indirecgt pathways in a network\\
\code{UncInd}         & \code{NetIndices} & Based on a flowmatrix\footnoterecall{flowmatrix}, estimates the statistical, realised and conditional uncertainty,...\\
\code{TrophInd}       & \code{NetIndices} & Based on a flowmatrix\footnoterecall{flowmatrix}, estimates the trophic level and omnivory indices\\

\hline
\end{tabular}
\end{table*}

\clearpage
\bibliography{vignettes}

\end{document}
