--- a/main/queueing/ChangeLog	Sat Mar 10 12:44:08 2012 +0000
+++ b/main/queueing/ChangeLog	Sat Mar 10 16:00:45 2012 +0000
@@ -1,3 +1,13 @@
+2012-02-XX Moreno Marzolla <marzolla@cs.unibo.it>
+
+	* Version 1.0.X released
+	* Fixed bug in qnvisits() which made the function behave incorrectly
+	  under particular degenerate cases.
+	* Fixed bug in ctmc_exps() (wrong initial value in call to lsode)
+	* Function ctmc_exps() can now also compute the expected sojourn time
+	  until absorption for absorbing CTMCs.
+	* Miscellaneous fixes/improvements to the documentation
+
 2012-02-04 Moreno Marzolla <marzolla@cs.unibo.it>
 
 	* Version 1.0.0 released under the name "queueing" (initial

--- a/main/queueing/NEWS	Sat Mar 10 12:44:08 2012 +0000
+++ b/main/queueing/NEWS	Sat Mar 10 16:00:45 2012 +0000
@@ -1,3 +1,8 @@
+Summary of important user-visible changes for queueing-1.0.X
+
+** Function ctmc_exps() can now compute the expected sojourn time
+   until absorption for absorming CTMC
+
 Summary of important user-visible changes for queueing-1.0.0
 ------------------------------------------------------------------------------
 

--- a/main/queueing/inst/ctmc_exps.m	Sat Mar 10 12:44:08 2012 +0000
+++ b/main/queueing/inst/ctmc_exps.m	Sat Mar 10 16:00:45 2012 +0000
@@ -17,24 +17,27 @@
 
 ## -*- texinfo -*-
 ##
-## @deftypefn {Function File} {@var{L} =} ctmc_exps (@var{Q}, @var{tt}, @var{p})
+## @deftypefn {Function File} {@var{L} =} ctmc_exps (@var{Q}, @var{tt}, @var{p} )
+## @deftypefnx {Function File} {@var{L} =} ctmc_exps (@var{Q}, @var{p})
 ##
 ## @cindex Markov chain, continuous time
 ## @cindex Expected sojourn time
 ##
-## Compute the expected total time @code{@var{L}(t,j)} spent in state
-## @math{j} during the time interval @code{[0,@var{tt}(t))}, assuming
-## that at time 0 the state occupancy probability was @var{p}.
+## With three arguments, compute the expected time @code{@var{L}(t,j)}
+## spent in each state @math{j} during the time interval
+## @code{[0,@var{tt}(t))}, assuming that at time 0 the state occupancy
+## probability was @var{p}. With two arguments, compute the expected
+## time @code{@var{L}(j}} spent in each state @math{j} until absorption.
 ##
 ## @strong{INPUTS}
 ##
 ## @table @var
 ##
 ## @item Q
-## Infinitesimal generator matrix. @code{@var{Q}(i,j)} is the transition
-## rate from state @math{i} to state @math{j},
-## @math{1 @leq{} i \neq j @leq{} N}. The matrix @var{Q} must also satisfy the
-## condition @code{sum(@var{Q},2) == 0}
+## @math{N \times N} infinitesimal generator matrix. @code{@var{Q}(i,j)}
+## is the transition rate from state @math{i} to state @math{j}, @math{1
+## @leq{} i \neq j @leq{} N}. The matrix @var{Q} must also satisfy the
+## condition @math{\sum_{j=1}^N Q_{ij} = 0}.
 ##
 ## @item tt
 ## This parameter is a vector used for numerical integration. The first
@@ -43,8 +46,9 @@
 ## @math{[0,t)} of interest (@code{@var{tt}(end) == @math{t}}).
 ##
 ## @item p
-## @code{@var{p}(i)} is the probability that at time 0 the system was in
-## state @math{i}, for all @math{i = 1, @dots{}, N}
+## Initial occupancy probability vector; @code{@var{p}(i)} is the
+## probability the system is in state @math{i} at time 0, @math{i = 1,
+## @dots{}, N}
 ##
 ## @end table
 ##
@@ -53,8 +57,14 @@
 ## @table @var
 ##
 ## @item L
-## @code{@var{L}(t,j)} is the expected time spent in state @math{j}
-## during the interval @code{[0,@var{tt}(t))}. @code{1 @leq{} @var{t} @leq{} length(@var{tt})}
+## If this function is called with three arguments, @var{L} is a matrix
+## of size @code{[length(@var{tt}), N]} where @code{@var{L}(t,j)} is the
+## expected time spent in state @math{j} during the interval
+## @code{[0,@var{tt}(t)]}. If this function is called with two
+## arguments, @var{L} is a vector with @math{N} elements where
+## @code{@var{L}(j)} is the expected time spent in state @math{j} until
+## absorption, if @math{j} is a transient state. If @math{j}
+## is an absorbing state, @code{@var{L}(j) = 0}.
 ##
 ## @end table
 ##
@@ -63,11 +73,11 @@
 ## Author: Moreno Marzolla <marzolla(at)cs.unibo.it>
 ## Web: http://www.moreno.marzolla.name/
 
-function L = ctmc_exps( Q, tt, p )
+function L = ctmc_exps( Q, varargin )
 
   persistent epsilon = 10*eps;
 
-  if ( nargin != 3 )
+  if ( nargin < 2 || nargin > 3 )
     print_usage();
   endif
 
@@ -77,16 +87,56 @@
   ( norm( sum(Q,2), "inf" ) < epsilon ) || \
       error( "Q is not an infinitesimal generator matrix" );
 
+  if ( nargin == 2 )
+    p = varargin{1};
+  else
+    tt = varargin{1};
+    p = varargin{2};
+  endif
+
   ( isvector(p) && length(p) == size(Q,1) && all(p>=0) && abs(sum(p)-1.0)<epsilon ) || \
       usage( "p must be a probability vector" );
 
-  ( isvector(tt) && abs(tt(1)) < epsilon ) || \
-      usage( "tt must be a vector, and tt(1) must be 0.0" );
-  tt = tt(:)'; # make tt a row vector
-  p = p(:)'; # make p a row vector
-  ff = @(x,t) (x(:)'*Q+p);
-  fj = @(x,t) (Q);
-  L = lsode( {ff, fj}, p, tt );
+  if ( nargin == 3 ) 
+    ( isvector(tt) && abs(tt(1)) < epsilon ) || \
+	usage( "tt must be a vector, and tt(1) must be 0.0" );
+    tt = tt(:)'; # make tt a row vector
+    p = p(:)'; # make p a row vector
+    ff = @(x,t) (x(:)'*Q+p);
+    fj = @(x,t) (Q);
+    L = lsode( {ff, fj}, zeros(size(p)), tt );
+  else
+#{
+    ## F(t) are the transient state occupancy probabilities at time t.
+    ## It is known that F(t) = p*expm(Q*t) (see function ctmc()).
+    ## The expected times spent in each state until absorption can
+    ## be computed as the integral of F(t) from t=0 to t=inf
+    F = @(t) (p*expm(Q*t)); ## FIXME: this must be restricted to transient states ONLY!!!!
+
+    ## Since function quadv does not support infinite integration
+    ## limits, we define a new function G(u) = F(tan(pi/2*u)) such that
+    ## the integral of G(u) on [0,1] is the integral of F(t) on [0,
+    ## +inf].
+    G = @(u) (F(tan(pi/2*u))*pi/2*(1+tan(pi/2*u)**2));
+
+    L = quadv(G,0,1);
+#}
+    ## Find nonzero rows. Nonzero rows correspond to transient states,
+    ## while zero rows are absorbing states. If there are no zero rows,
+    ## then the Markov chain does not contain absorbing states and we
+    ## raise an error
+    N = rows(Q);
+    nzrows = find( any( abs(Q) > epsilon, 2 ) );
+    if ( length( nzrows ) == N )
+      error( "There are no absorbing states" );
+    endif
+    
+    QN = Q(nzrows,nzrows);
+    pN = p(nzrows);
+    LN = -pN*inv(QN);
+    L = zeros(1,N);
+    L(nzrows) = LN;
+  endif
 endfunction
 
 %!demo
@@ -99,10 +149,25 @@
 %! tt = linspace(0,10,100);
 %! p0 = zeros(1,N); p0(1)=1;
 %! L = ctmc_exps(Q,tt,p0);
+%! #L2 = 0*L;
+%! #for i=1:length(tt)
+%! #  L2(i,:) = quadv( @(t) (p0*expm(Q*t)) , 0, tt(i) );
+%! #endfor
 %! plot( tt, L(:,1), ";State 1;", "linewidth", 2, \
+%! #      tt, L2(:,1), "+;State 1 (quadv);", "markersize", 8, \
 %!       tt, L(:,2), ";State 2;", "linewidth", 2, \
 %!       tt, L(:,3), ";State 3;", "linewidth", 2, \
 %!       tt, L(:,4), ";State 4 (absorbing);", "linewidth", 2);
 %! legend("location","northwest");
 %! xlabel("Time");
 %! ylabel("Expected sojourn time");
+
+%!demo
+%! lambda = 0.5;
+%! N = 4;
+%! birth = lambda*linspace(1,N-1,N-1);
+%! death = 0*birth;
+%! Q = diag(birth,1)+diag(death,-1);
+%! Q -= diag(sum(Q,2));
+%! p0 = zeros(1,N); p0(1)=1;
+%! L = ctmc_exps(Q,p0)

--- a/main/queueing/inst/ctmc_mtta.m	Sat Mar 10 12:44:08 2012 +0000
+++ b/main/queueing/inst/ctmc_mtta.m	Sat Mar 10 16:00:45 2012 +0000
@@ -22,9 +22,10 @@
 ## @cindex Markov chain, continuous time
 ## @cindex Mean time to absorption
 ##
-## Compute the Mean-Time to Absorption (MTTA) starting from initial
-## occupancy probability @var{p} at time 0. If there are no absorbing
-## states, this function fails with an error.
+## Compute the Mean-Time to Absorption (MTTA) of the CTMC described by
+## the infinitesimal generator matrix @var{Q}, starting from initial
+## occupancy probability @var{p}. If there are no absorbing states, this
+## function fails with an error.
 ##
 ## @strong{INPUTS}
 ##
@@ -76,18 +77,7 @@
   ( isvector(p) && length(p) == N && all(p>=0) && abs(sum(p)-1.0)<epsilon ) || \
       usage( "p must be a probability vector" );
 
-  ## Find nonzero rows. Nonzero rows correspond to transient states,
-  ## while zero rows are absorbing states. If there are no zero rows,
-  ## then the Markov chain does not contain absorbing states and we
-  ## raise an error
-  nzrows = find( any( abs(Q) > epsilon, 2 ) );
-  if ( length( nzrows ) == N )
-    error( "There are no absorbing states" );
-  endif
-
-  QN = Q(nzrows,nzrows);
-  pN = p(nzrows);
-  L = -pN*inv(QN);
+  L = ctmc_exps(Q,p);
   t = sum(L);
 endfunction
 %!test
@@ -122,7 +112,7 @@
 %! death = [ 3 4 5 ] * mu;
 %! Q = diag(death,-1);
 %! Q -= diag(sum(Q,2));
-%! t = ctmc_mtta(Q,[0 0 0 1])
+%! [t L] = ctmc_mtta(Q,[0 0 0 1])
 
 %!demo
 %! N = 100;