Mercurial > forge

changeset 9625:c2e1990614a2 octave-forge

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
New signature for ctmc_exps() and ctmc_taexps()
author mmarzolla
date Sat, 10 Mar 2012 20:27:25 +0000
parents dbea4e0ca957
children c4728d463424
files main/queueing/ChangeLog main/queueing/NEWS main/queueing/doc/queueing.html main/queueing/doc/queueing.pdf main/queueing/inst/ctmc_exps.m main/queueing/inst/ctmc_mtta.m main/queueing/inst/ctmc_taexps.m
diffstat 7 files changed, 113 insertions(+), 105 deletions(-) [+]
line wrap: on
line diff
--- a/main/queueing/ChangeLog	Sat Mar 10 19:46:27 2012 +0000
+++ b/main/queueing/ChangeLog	Sat Mar 10 20:27:25 2012 +0000
@@ -6,6 +6,9 @@
 	* 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.
+	* Function ctmc_exps() and ctmc_taexps() accept a scalar as second
+	  argument; the old syntax is still supported, but may be deprecated
+	  in future releases.
 	* Miscellaneous fixes/improvements to the documentation
 
 2012-02-04 Moreno Marzolla <marzolla@cs.unibo.it>
--- a/main/queueing/NEWS	Sat Mar 10 19:46:27 2012 +0000
+++ b/main/queueing/NEWS	Sat Mar 10 20:27:25 2012 +0000
@@ -3,6 +3,10 @@
 ** Function ctmc_exps() can now compute the expected sojourn time
    until absorption for absorming CTMC
 
+** Functions ctmc_exps() and ctmc_taexps() now accept a scalar as
+   the second argument (time). The old syntax is still supported,
+   but may be deprecated in the future.
+
 Summary of important user-visible changes for queueing-1.0.0
 ------------------------------------------------------------------------------
 
--- a/main/queueing/doc/queueing.html	Sat Mar 10 19:46:27 2012 +0000
+++ b/main/queueing/doc/queueing.html	Sat Mar 10 20:27:25 2012 +0000
@@ -1100,15 +1100,15 @@
    <p><a name="doc_002dctmc_005fexps"></a>
 
 <div class="defun">
-&mdash; Function File: <var>L</var> = <b>ctmc_exps</b> (<var>Q, tt, p </var>)<var><a name="index-ctmc_005fexps-20"></a></var><br>
+&mdash; Function File: <var>L</var> = <b>ctmc_exps</b> (<var>Q, t, p </var>)<var><a name="index-ctmc_005fexps-20"></a></var><br>
 &mdash; Function File: <var>L</var> = <b>ctmc_exps</b> (<var>Q, p</var>)<var><a name="index-ctmc_005fexps-21"></a></var><br>
 <blockquote>
         <p><a name="index-Markov-chain_002c-continuous-time-22"></a><a name="index-Expected-sojourn-time-23"></a>
-With three arguments, compute the expected time <var>L</var><code>(t,j)</code>
-spent in each state j during the time interval
-<code>[0,</code><var>tt</var><code>(t))</code>, assuming that at time 0 the state occupancy
-probability was <var>p</var>. With two arguments, compute the expected
-time <var>L</var><code>(j)</code> spent in each state j until absorption.
+With three arguments, compute the expected times <var>L</var><code>(i)</code>
+spent in each state i during the time interval
+[0,t], assuming that the state occupancy probabilities
+at time 0 are <var>p</var>. With two arguments, compute the expected time
+<var>L</var><code>(i)</code> spent in each state i until absorption.
 
         <p><strong>INPUTS</strong>
 
@@ -1118,10 +1118,7 @@
 &le; i \neq j &le; N. The matrix <var>Q</var> must also satisfy the
 condition \sum_j=1^N Q_ij = 0.
 
-          <br><dt><var>tt</var><dd>This parameter is a vector used for numerical integration. The first
-element <var>tt</var><code>(1)</code> must be 0, and the last element
-<var>tt</var><code>(end)</code> must be the upper bound of the interval
-[0,t) of interest (<var>tt</var><code>(end) == t</code>).
+          <br><dt><var>t</var><dd>Time
 
           <br><dt><var>p</var><dd>Initial occupancy probability vector; <var>p</var><code>(i)</code> is the
 probability the system is in state i at time 0, i = 1,
@@ -1132,14 +1129,12 @@
         <p><strong>OUTPUTS</strong>
 
           <dl>
-<dt><var>L</var><dd>If this function is called with three arguments, <var>L</var> is a matrix
-of size <code>[length(</code><var>tt</var><code>), N]</code> where <var>L</var><code>(t,j)</code> is the
-expected time spent in state j during the interval
-<code>[0,</code><var>tt</var><code>(t)]</code>. If this function is called with two
-arguments, <var>L</var> is a vector with N elements where
-<var>L</var><code>(j)</code> is the expected time spent in state j until
-absorption, if j is a transient state. If j
-is an absorbing state, <var>L</var><code>(j) = 0</code>.
+<dt><var>L</var><dd>If this function is called with three arguments, <var>L</var><code>(i)</code> is
+the expected time spent in state j during the interval
+[0,t]. If this function is called with two arguments
+<var>L</var><code>(i)</code> is the expected time spent in state i until
+absorption (if i is a transient state), or zero
+(if <var>i</var> is an absorbing state).
 
         </dl>
 
@@ -1185,13 +1180,13 @@
 <p><a name="doc_002dctmc_005ftaexps"></a>
 
 <div class="defun">
-&mdash; Function File: <var>M</var> = <b>ctmc_taexps</b> (<var>Q, tt, p</var>)<var><a name="index-ctmc_005ftaexps-24"></a></var><br>
+&mdash; Function File: <var>M</var> = <b>ctmc_taexps</b> (<var>Q, t, p</var>)<var><a name="index-ctmc_005ftaexps-24"></a></var><br>
 <blockquote>
         <p><a name="index-Markov-chain_002c-continuous-time-25"></a><a name="index-Time_002dalveraged-sojourn-time-26"></a>
-Compute the <em>time-averaged sojourn time</em> <var>M</var><code>(t,j)</code>,
-defined as the fraction of the time interval <code>[0,</code><var>tt</var><code>(t))</code> spent in
-state j, assuming that at time 0 the state occupancy
-probability was <var>p</var>.
+Compute the <em>time-averaged sojourn time</em> <var>M</var><code>(i)</code>,
+defined as the fraction of the time interval [0,t] spent in
+state i, assuming that the state occupancy probabilities at
+time 0 are <var>p</var>.
 
         <p><strong>INPUTS</strong>
 
@@ -1199,15 +1194,9 @@
 <dt><var>Q</var><dd>Infinitesimal generator matrix. <var>Q</var><code>(i,j)</code> is the transition
 rate from state i to state j,
 1 &le; i \neq j &le; N. The
-matrix <var>Q</var> must also satisfy the condition <code>sum(</code><var>Q</var><code>,2) == 0</code>
-
-          <br><dt><var>tt</var><dd>This parameter is a vector used for numerical integration of the
-sujourn time. The first element <var>tt</var><code>(1)</code> must be slightly
-larger than 0, and the
-last element <var>tt</var><code>(end)</code> must be the upper limit of the
-interval [0,t) of interest (<var>tt</var><code>(end) == t</code>). 
-This vector is used by the ODE solver to compute the solution
-<var>M</var>.
+matrix <var>Q</var> must also satisfy the condition \sum_j=1^N Q_ij = 0
+
+          <br><dt><var>t</var><dd>Time
 
           <br><dt><var>p</var><dd><var>p</var><code>(i)</code> is the probability that, at time 0, the system was in
 state i, for all i = 1, <small class="dots">...</small>, N
@@ -1217,10 +1206,9 @@
         <p><strong>OUTPUTS</strong>
 
           <dl>
-<dt><var>M</var><dd><var>M</var><code>(t,j)</code> is the expected fraction of time spent in state
-j during the interval [0,tt(t)) assuming that the state
-occupancy probability at time zero was <var>p</var>. <code>1 &le;
-</code><var>t</var><code> &le; length(</code><var>tt</var><code>)</code>
+<dt><var>M</var><dd><var>M</var><code>(i)</code> is the expected fraction of time spent in state
+i during the interval [0,t] assuming that the state
+occupancy probability at time zero is <var>p</var>.
 
         </dl>
 
@@ -1279,7 +1267,7 @@
         <p><a name="index-Markov-chain_002c-continuous-time-28"></a><a name="index-Mean-time-to-absorption-29"></a>
 Compute the Mean-Time to Absorption (MTTA) of the CTMC described by
 the infinitesimal generator matrix <var>Q</var>, starting from initial
-occupancy probability <var>p</var>. If there are no absorbing states, this
+occupancy probabilities <var>p</var>. If there are no absorbing states, this
 function fails with an error.
 
         <p><strong>INPUTS</strong>
Binary file main/queueing/doc/queueing.pdf has changed
--- a/main/queueing/inst/ctmc_exps.m	Sat Mar 10 19:46:27 2012 +0000
+++ b/main/queueing/inst/ctmc_exps.m	Sat Mar 10 20:27:25 2012 +0000
@@ -17,17 +17,17 @@
 
 ## -*- texinfo -*-
 ##
-## @deftypefn {Function File} {@var{L} =} ctmc_exps (@var{Q}, @var{tt}, @var{p} )
+## @deftypefn {Function File} {@var{L} =} ctmc_exps (@var{Q}, @var{t}, @var{p} )
 ## @deftypefnx {Function File} {@var{L} =} ctmc_exps (@var{Q}, @var{p})
 ##
 ## @cindex Markov chain, continuous time
 ## @cindex Expected sojourn time
 ##
-## 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.
+## With three arguments, compute the expected times @code{@var{L}(i)}
+## spent in each state @math{i} during the time interval
+## @math{[0,t]}, assuming that the state occupancy probabilities
+## at time 0 are @var{p}. With two arguments, compute the expected time
+## @code{@var{L}(i)} spent in each state @math{i} until absorption.
 ##
 ## @strong{INPUTS}
 ##
@@ -39,11 +39,8 @@
 ## @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
-## element @code{@var{tt}(1)} must be 0, and the last element
-## @code{@var{tt}(end)} must be the upper bound of the interval
-## @math{[0,t)} of interest (@code{@var{tt}(end) == @math{t}}).
+## @item t
+## Time
 ##
 ## @item p
 ## Initial occupancy probability vector; @code{@var{p}(i)} is the
@@ -57,14 +54,12 @@
 ## @table @var
 ##
 ## @item L
-## 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}.
+## If this function is called with three arguments, @code{@var{L}(i)} is
+## the expected time spent in state @math{j} during the interval
+## @math{[0,t]}. If this function is called with two arguments
+## @code{@var{L}(i)} is the expected time spent in state @math{i} until
+## absorption (if @math{i} is a transient state), or zero
+## (if @var{i} is an absorbing state).
 ##
 ## @end table
 ##
@@ -85,12 +80,12 @@
       usage( "Q must be a square matrix" );
 
   ( norm( sum(Q,2), "inf" ) < epsilon ) || \
-      error( "Q is not an infinitesimal generator matrix" );
+      usage( "Q is not an infinitesimal generator matrix" );
 
   if ( nargin == 2 )
     p = varargin{1};
   else
-    tt = varargin{1};
+    t = varargin{1};
     p = varargin{2};
   endif
 
@@ -98,15 +93,27 @@
       usage( "p must be a probability vector" );
 
   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 );
+    if ( isscalar(t) )
+      (t >= 0 ) || \
+	  usage( "t must be >= 0" );
+      ## F(x) are the transient state occupancy probabilities at time x.
+      ## It is known that F(x) = p*expm(Q*x) (see function ctmc()).
+      F = @(x) (p*expm(Q*x));
+      L = quadv(F,0,t);
+    else
+      ## FIXME: deprecate this?
+      ( isvector(t) && abs(t(1)) < epsilon ) || \
+	  usage( "t must be a vector, and t(1) must be 0.0" );
+      t = t(:)'; # 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)), t );
+    endif
   else
 #{
+    ## This code is left for information only
+
     ## 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
@@ -138,6 +145,10 @@
     L(nzrows) = LN;
   endif
 endfunction
+%!test
+%! Q = [-1 1; 1 -1];
+%! L = ctmc_exps(Q,10,[1 0]);
+%! L = ctmc_exps(Q,linspace(0,10,100),[1 0]);
 
 %!demo
 %! lambda = 0.5;
@@ -146,18 +157,16 @@
 %! death = zeros(1,N-1);
 %! Q = diag(birth,1)+diag(death,-1);
 %! Q -= diag(sum(Q,2));
-%! tt = linspace(0,10,100);
+%! t = 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);
+%! L = zeros(length(t),N);
+%! for i=1:length(t)
+%!   L(i,:) = ctmc_exps(Q,t(i),p0);
+%! endfor
+%! plot( t, L(:,1), ";State 1;", "linewidth", 2, \
+%!       t, L(:,2), ";State 2;", "linewidth", 2, \
+%!       t, L(:,3), ";State 3;", "linewidth", 2, \
+%!       t, L(:,4), ";State 4;", "linewidth", 2 );
 %! legend("location","northwest");
 %! xlabel("Time");
 %! ylabel("Expected sojourn time");
--- a/main/queueing/inst/ctmc_mtta.m	Sat Mar 10 19:46:27 2012 +0000
+++ b/main/queueing/inst/ctmc_mtta.m	Sat Mar 10 20:27:25 2012 +0000
@@ -24,7 +24,7 @@
 ##
 ## 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
+## occupancy probabilities @var{p}. If there are no absorbing states, this
 ## function fails with an error.
 ##
 ## @strong{INPUTS}
@@ -71,8 +71,8 @@
   
   N = rows(Q);
 
-  all( abs( sum(Q,2) ) < epsilon )  || \
-      usage( "Q is not an infinitesimal generator matrix" );
+  ( norm( sum(Q,2), "inf" ) < epsilon ) || \
+      usage( "Q must be an infinitesimal generator matrix" );
 
   ( isvector(p) && length(p) == N && all(p>=0) && abs(sum(p)-1.0)<epsilon ) || \
       usage( "p must be a probability vector" );
@@ -90,7 +90,7 @@
 
 %!test
 %! Q = [0 1 0; 1 0 1; 0 0 0 ];
-%! fail( "ctmc_mtta(Q,[1 0 0])", "not an infinitesimal");
+%! fail( "ctmc_mtta(Q,[1 0 0])", "must be an infinitesimal");
 
 %!test
 %! Q = [ 0 0.1 0 0; \
--- a/main/queueing/inst/ctmc_taexps.m	Sat Mar 10 19:46:27 2012 +0000
+++ b/main/queueing/inst/ctmc_taexps.m	Sat Mar 10 20:27:25 2012 +0000
@@ -17,15 +17,15 @@
 
 ## -*- texinfo -*-
 ##
-## @deftypefn {Function File} {@var{M} =} ctmc_taexps (@var{Q}, @var{tt}, @var{p})
+## @deftypefn {Function File} {@var{M} =} ctmc_taexps (@var{Q}, @var{t}, @var{p})
 ##
 ## @cindex Markov chain, continuous time
 ## @cindex Time-alveraged sojourn time
 ##
-## Compute the @emph{time-averaged sojourn time} @code{@var{M}(t,j)},
-## defined as the fraction of the time interval @code{[0,@var{tt}(t))} spent in
-## state @math{j}, assuming that at time 0 the state occupancy
-## probability was @var{p}.
+## Compute the @emph{time-averaged sojourn time} @code{@var{M}(i)},
+## defined as the fraction of the time interval @math{[0,t]} spent in
+## state @math{i}, assuming that the state occupancy probabilities at
+## time 0 are @var{p}.
 ##
 ## @strong{INPUTS}
 ##
@@ -35,16 +35,10 @@
 ## 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}
+## 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 of the
-## sujourn time. The first element @code{@var{tt}(1)} must be slightly
-## larger than 0, and the
-## last element @code{@var{tt}(end)} must be the upper limit of the
-## interval @math{[0,t)} of interest (@code{@var{tt}(end) == @math{t}}).
-## This vector is used by the ODE solver to compute the solution
-## @var{M}.
+## @item t
+## Time
 ##
 ## @item p
 ## @code{@var{p}(i)} is the probability that, at time 0, the system was in
@@ -57,10 +51,9 @@
 ## @table @var
 ##
 ## @item M
-## @code{@var{M}(t,j)} is the expected fraction of time spent in state
-## @math{j} during the interval @math{[0,tt(t))} assuming that the state
-## occupancy probability at time zero was @var{p}. @code{1 @leq{}
-## @var{t} @leq{} length(@var{tt})}
+## @code{@var{M}(i)} is the expected fraction of time spent in state
+## @math{i} during the interval @math{[0,t]} assuming that the state
+## occupancy probability at time zero is @var{p}.
 ##
 ## @end table
 ##
@@ -87,11 +80,19 @@
 
   ( isvector(p) && length(p) == N && all(p>=0) && abs(sum(p)-1.0)<epsilon ) || \
       usage( "p must be a probability vector" );
-  t = t(:)'; # make t a row vector
-  p = p(:)'; # make p a row vector
-  ff = @(x,t) (((x')*(Q-eye(N)/t).+p/t)');
-  fj = @(x,t) (Q-eye(N)/t);
-  M = lsode( {ff, fj}, zeros(size(p)), t );
+  if ( isscalar(t) )
+    (t >= 0) || \
+	usage( "t must be >= 0" );
+    F = @(x) (p*expm(Q*x));
+    M = quadv(F,0,t) / t;
+  else
+    ## FIXME: deprecate this?
+    t = t(:)'; # make t a row vector
+    p = p(:)'; # make p a row vector
+    ff = @(x,t) (((x')*(Q-eye(N)/t).+p/t)');
+    fj = @(x,t) (Q-eye(N)/t);
+    M = lsode( {ff, fj}, zeros(size(p)), t );
+  endif
 endfunction
 
 %!demo
@@ -101,9 +102,12 @@
 %! death = zeros(1,N-1);
 %! Q = diag(birth,1)+diag(death,-1);
 %! Q -= diag(sum(Q,2));
-%! t = linspace(1e-3,50,500);
+%! t = linspace(1e-5,30,100);
 %! p = zeros(1,N); p(1)=1;
-%! M = ctmc_taexps(Q,t,p);
+%! M = zeros(length(t),N);
+%! for i=1:length(t)
+%!   M(i,:) = ctmc_taexps(Q,t(i),p);
+%! endfor
 %! plot(t, M(:,1), ";State 1;", "linewidth", 2, \
 %!      t, M(:,2), ";State 2;", "linewidth", 2, \
 %!      t, M(:,3), ";State 3;", "linewidth", 2, \