[Octave-cvsupdate] SF.net SVN: octave:[9979] trunk/octave-forge/main/queueing

[<mma...@us...>]

Revision: 9979
          http://octave.svn.sourceforge.net/octave/?rev=9979&view=rev
Author:   mmarzolla
Date:     2012-03-20 16:20:06 +0000 (Tue, 20 Mar 2012)
Log Message:
-----------
code restructuring

Modified Paths:
--------------
    trunk/octave-forge/main/queueing/doc/references.txi
    trunk/octave-forge/main/queueing/inst/ctmc.m
    trunk/octave-forge/main/queueing/inst/ctmc_check_Q.m
    trunk/octave-forge/main/queueing/inst/ctmc_fpt.m
    trunk/octave-forge/main/queueing/inst/ctmc_mtta.m
    trunk/octave-forge/main/queueing/inst/dtmc.m

Modified: trunk/octave-forge/main/queueing/doc/references.txi
===================================================================
--- trunk/octave-forge/main/queueing/doc/references.txi	2012-03-20 14:03:24 UTC (rev 9978)
+++ trunk/octave-forge/main/queueing/doc/references.txi	2012-03-20 16:20:06 UTC (rev 9979)
@@ -108,6 +108,11 @@
 Solution of Queueing Network Models}, @uref{http://www.cs.purdue.edu/research/technical_reports/1980/TR%2080-355.pdf, Technical Report CSD-TR-355},
 Department of Computer Sciences, Purdue University, feb 15, 1982.
 
+@items [Tij03]
+H. C. Tijms, @cite{A first course in stochastic models},
+John Wiley and Sons, 2003, ISBN 0471498807, ISBN 9780471498803,
+DOI @uref{http://dx.doi.org/10.1002/047001363X, 10.1002/047001363X}
+
 @item [ZaWo81]
 Zahorjan, J. and Wong, E. @cite{The solution of separable queueing
 network models using mean value analysis}. SIGMETRICS

Modified: trunk/octave-forge/main/queueing/inst/ctmc.m
===================================================================
--- trunk/octave-forge/main/queueing/inst/ctmc.m	2012-03-20 14:03:24 UTC (rev 9978)
+++ trunk/octave-forge/main/queueing/inst/ctmc.m	2012-03-20 16:20:06 UTC (rev 9979)
@@ -48,7 +48,7 @@
 ##
 ## @item p0
 ## @code{@var{p0}(i)} is the probability that the system
-## is in state @math{i} at time 0 .
+## is in state @math{i} at time 0.
 ##
 ## @end table
 ##
@@ -74,9 +74,7 @@
 
 function q = ctmc( Q, t, p0 )
 
-  persistent epsilon = 10*eps;
-
-  if ( nargin < 1 || nargin > 3 )
+  if ( nargin != 1 && nargin != 3 )
     print_usage();
   endif
 
@@ -85,22 +83,17 @@
   ( N>0 ) || \
       usage(err);
 
-  if ( nargin > 1 ) 
+  if ( nargin == 1 )
+    q = __ctmc_steady_state( Q );
+  else
     ( isscalar(t) && t>=0 ) || \
         usage("t must be nonnegative");
-  endif
 
-  if ( nargin > 2 )
-    ( isvector(p0) && length(p0) == N && all(p0>=0) && abs(sum(p0)-1.0)<epsilon ) || \
+    ( isvector(p0) && length(p0) == N && all(p0>=0) && abs(sum(p0)-1.0)<N*eps ) || \
         usage( "p0 must be a probability vector" );   
+
     p0 = p0(:)'; # make p0 a row vector
-  else
-    p0 = ones(1,N) / N;
-  endif
 
-  if ( nargin == 1 )
-    q = __ctmc_steady_state( Q );
-  else
     q = __ctmc_transient(Q, t, p0 );
   endif
 
@@ -198,6 +191,32 @@
 %! q = ctmc(Q, 100, [1 0]);
 %! assert( qlim, q, 1e-5 );
 
+## Example on p. 172 of [Tij03]
+%!test
+%! ll = 0.1;
+%! mu = 100;
+%! eta = 5;
+%! Q = zeros(9,9);
+%! ## 6--1, 7=sleep2 8=sleep1 9=crash
+%! Q(6,5) = 6*ll;
+%! Q(5,4) = 5*ll;
+%! Q(4,3) = 4*ll;
+%! Q(3,2) = 3*ll;
+%! Q(2,1) = 2*ll;
+%! Q(2,7) = mu;
+%! Q(1,9) = ll;
+%! Q(1,8) = mu;
+%! Q(8,9) = ll;
+%! Q(7,8) = 2*ll;
+%! Q(7,6) = eta;
+%! Q(8,6) = eta;
+%! Q -= diag(sum(Q,2));
+%! q0 = zeros(1,9); q0(6) = 1;
+%! q = ctmc(Q,10,q0);
+%! assert( q(9), 0.000504, 1e-6 );
+%! q = ctmc(Q,2,q0);
+%! assert( q, [3.83e-7 1.938e-4 0.0654032 0.2216998 0.4016008 0.3079701 0.0030271 0.0000998 5e-6], 1e-5 );
+
 %!demo
 %! Q = [ -1  1; \
 %!        1 -1  ];

Modified: trunk/octave-forge/main/queueing/inst/ctmc_check_Q.m
===================================================================
--- trunk/octave-forge/main/queueing/inst/ctmc_check_Q.m	2012-03-20 14:03:24 UTC (rev 9978)
+++ trunk/octave-forge/main/queueing/inst/ctmc_check_Q.m	2012-03-20 16:20:06 UTC (rev 9979)
@@ -33,7 +33,7 @@
 
 function [result err] = ctmc_check_Q( Q )
 
-  persistent epsilon = 10*eps;
+  persistent epsilon = 100*eps;
 
   if ( nargin != 1 )
     print_usage();
@@ -46,7 +46,7 @@
     return;
   endif
   
-  if (any(Q(~logical(eye(size(Q))))<0) || \ # there is any negavite non-diagonal element
+  if (any(Q(~logical(eye(size(Q))))<0) || \ # there is any negative non-diagonal element
       norm( sum(Q,2), "inf" ) > epsilon )
     err = "Q is not an infinitesimal generator matrix";
     return;

Modified: trunk/octave-forge/main/queueing/inst/ctmc_fpt.m
===================================================================
--- trunk/octave-forge/main/queueing/inst/ctmc_fpt.m	2012-03-20 14:03:24 UTC (rev 9978)
+++ trunk/octave-forge/main/queueing/inst/ctmc_fpt.m	2012-03-20 16:20:06 UTC (rev 9979)
@@ -79,16 +79,15 @@
       error(err);
 
   if ( nargin == 1 ) 
-    M = zeros(N,N);
+    result = zeros(N,N);
     for j=1:N
       QQ = Q;
       QQ(j,:) = 0; # make state j absorbing
       for i=1:N
 	p0 = zeros(1,N); p0(i) = 1;
-	M(i,j) = ctmc_mtta(QQ,p0);
+	result(i,j) = ctmc_mtta(QQ,p0);
       endfor
     endfor
-    result = M;
   else
     (isscalar(i) && i>=1 && j<=N) || usage("i must be an integer in the range 1..%d", N);
     (isvector(j) && all(j>=1) && all(j<=N)) || usage("j must be an integer or vector with elements in 1..%d", N);
@@ -112,13 +111,3 @@
 %! Q -= diag(sum(Q,2));
 %! M = ctmc_fpt(Q);
 %! assert( all(diag(M) < 10*eps) );
-
-%!xtest
-%! Q = unifrnd(0.1,0.9,10,10);
-%! Q -= diag(sum(Q,2));
-%! m = ctmc_fpt(Q,1,3);
-
-%!xtest
-%! Q = unifrnd(0.1,0.9,10,10);
-%! Q -= diag(sum(Q,2));
-%! m = ctmc_fpt(Q,1,[3 5 6]);

Modified: trunk/octave-forge/main/queueing/inst/ctmc_mtta.m
===================================================================
--- trunk/octave-forge/main/queueing/inst/ctmc_mtta.m	2012-03-20 14:03:24 UTC (rev 9978)
+++ trunk/octave-forge/main/queueing/inst/ctmc_mtta.m	2012-03-20 16:20:06 UTC (rev 9979)
@@ -127,3 +127,5 @@
 %! plot(initial_state,t,"+");
 %! xlabel("Initial state");
 %! ylabel("MTTA");
+
+

Modified: trunk/octave-forge/main/queueing/inst/dtmc.m
===================================================================
--- trunk/octave-forge/main/queueing/inst/dtmc.m	2012-03-20 14:03:24 UTC (rev 9978)
+++ trunk/octave-forge/main/queueing/inst/dtmc.m	2012-03-20 16:20:06 UTC (rev 9979)
@@ -46,7 +46,8 @@
 ## @var{P} must be equal to its dimension.
 ##
 ## @item n
-## Step at which to compute the transient probability
+## Number of transitions after which compute the state occupancy probabilities
+## (@math{n=0, 1, @enddots{}})
 ##
 ## @item p0
 ## @code{@var{p0}(i)} is the probability that at step 0 the system
@@ -63,7 +64,7 @@
 ## @code{@var{p}(i)} is the steady-state probability that the system is
 ## in state @math{i}. @var{p} satisfies the equations @math{p = p{\bf P}} and @math{\sum_{i=1}^N p_i = 1}. If this function is invoked
 ## with three arguments, @code{@var{p}(i)} is the marginal probability
-## that the system is in state @math{i} at step @var{n},
+## that the system is in state @math{i} after @var{n} transitions,
 ## given the initial probabilities @code{@var{p0}(i)} that the initial state is
 ## @math{i}.
 ##
@@ -74,11 +75,9 @@
 ## Author: Moreno Marzolla <marzolla(at)cs.unibo.it>
 ## Web: http://www.moreno.marzolla.name/
 
-function q = dtmc( P, n, p0 )
+function p = dtmc( P, n, p0 )
 
-  persistent epsilon = 10*eps;
-
-  if ( nargin < 1 || nargin > 3 )
+  if ( nargin != 1 && nargin != 3 )
     print_usage();
   endif
 
@@ -86,29 +85,24 @@
   
   ( N>0 ) || \
       usage( err );
-  
-  if ( nargin > 1 )
+
+  if ( nargin == 1 )
+    p = __dtmc_steady_state( P );  
+  else
     ( isscalar(n) && n>=0 ) || \
 	usage( "n must be >=0" );
-  endif
 
-  if ( nargin > 2 )
-    ( isvector(p0) && length(p0) == N && all(p0>=0) && abs(sum(p0)-1.0)<epsilon ) || \
+    ( isvector(p0) && length(p0) == N && all(p0>=0) && abs(sum(p0)-1.0)<N*eps ) || \
         usage( "p0 must be a probability vector" );   
-    p0 = p0(:)'; # make q0 a row vector
-  else
-    p0 = ones(1,N) / N;
-  endif
 
-  if ( nargin == 1 )
-    q = __dtmc_steady_state( P );
-  else
-    q = __dtmc_transient(P, n, p0);
+    p0 = p0(:)'; # make p0 a row vector
+
+    p = __dtmc_transient(P, n, p0);
   endif
 endfunction
 
 ## Helper function, compute steady-state probability
-function q = __dtmc_steady_state( P )
+function p = __dtmc_steady_state( P )
   N = rows(P);
   A = P-eye(N);
   A(:,N) = 1; # add normalization condition
@@ -116,25 +110,25 @@
       warning( "dtmc(): P is reducible" );
 
   b = [ zeros(1,N-1) 1 ];
-  q = b/A;
+  p = b/A;
 endfunction
 
 ## Helper function, compute transient probability
-function q = __dtmc_transient( P, n, p0 )
-  q = p0*P^n;
+function p = __dtmc_transient( P, n, p0 )
+  p = p0*P^n;
 endfunction
 
 %!test
 %! P = [0.75 0.25; 0.5 0.5];
-%! q = dtmc(P);
-%! assert( q*P, q, 1e-5 );
-%! assert( q, [0.6666 0.3333], 1e-4 );
+%! p = dtmc(P);
+%! assert( p*P, p, 1e-5 );
+%! assert( p, [0.6666 0.3333], 1e-4 );
 
 %!test
 %! #Example 2.11 p. 44 Bolch et al.
 %! P = [0.5 0.5; 0.5 0.5];
-%! q = dtmc(P);
-%! assert( q, [0.5 0.5], 1e-3 );
+%! p = dtmc(P);
+%! assert( p, [0.5 0.5], 1e-3 );
 
 %!test
 %! fail("dtmc( [1 1 1; 1 1 1] )", "square");

This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.