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

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

Revision: 10146
          http://octave.svn.sourceforge.net/octave/?rev=10146&view=rev
Author:   mmarzolla
Date:     2012-04-05 14:32:41 +0000 (Thu, 05 Apr 2012)
Log Message:
-----------
Bug fixes

Modified Paths:
--------------
    trunk/octave-forge/main/queueing/inst/ctmc.m
    trunk/octave-forge/main/queueing/inst/ctmc_taexps.m
    trunk/octave-forge/main/queueing/inst/dtmc.m
    trunk/octave-forge/main/queueing/inst/dtmc_fpt.m
    trunk/octave-forge/main/queueing/inst/qnmmm.m

Modified: trunk/octave-forge/main/queueing/inst/ctmc.m
===================================================================
--- trunk/octave-forge/main/queueing/inst/ctmc.m	2012-04-04 12:47:24 UTC (rev 10145)
+++ trunk/octave-forge/main/queueing/inst/ctmc.m	2012-04-05 14:32:41 UTC (rev 10146)
@@ -25,11 +25,14 @@
 ## @cindex Markov chain, state occupancy probabilities
 ## @cindex Stationary probabilities
 ##
+## Compute stationary or transient state occupancy probabilities
+## for a continuous-time Markov chain.
+##
 ## With a single argument, compute the stationary state occupancy
 ## probability vector @var{p}(1), @dots{}, @var{p}(N) for a
-## Continuous-Time Markov Chain with infinitesimal generator matrix
-## @var{Q} of size  @math{N \times N}. With three arguments, compute the
-## state occupancy probabilities @var{p}(1), @dots{}, @var{p}(N) at time
+## continuous-time Markov chain with @math{N \times N} infinitesimal
+## generator matrix @var{Q}. With three arguments, compute the state
+## occupancy probabilities @var{p}(1), @dots{}, @var{p}(N) at time
 ## @var{t}, given initial state occupancy probabilities @var{p0} at time
 ## 0.
 ##
@@ -74,6 +77,8 @@
 
 function q = ctmc( Q, t, p0 )
 
+  persistent epsilon = 10*eps;
+
   if ( nargin != 1 && nargin != 3 )
     print_usage();
   endif
@@ -83,52 +88,43 @@
   ( N>0 ) || \
       usage(err);
 
-  if ( nargin == 1 )
-    q = __ctmc_steady_state( Q );
-  else
+  if ( nargin == 1 ) # steady-state analysis
+
+    ## non zero columns
+    nonzero=find( any(abs(Q)>epsilon,1 ) );
+    if ( length(nonzero) == 0 )
+      error( "Q is the zero matrix" );
+    endif
+    
+    normcol = nonzero(1); # normalization condition column
+    
+    ## force probability of unvisited states to zero
+    for i=find( all(abs(Q)<epsilon,1) )
+      Q(i,i) = 1;
+    endfor
+    
+    ## assert( rank(Q) == N-1 );
+    
+    Q(:,normcol) = 1; # add normalization condition
+    b = zeros(1,N); b(normcol)=1;
+    q = b/Q; # qQ = b;
+    
+  else # transient analysis
+
     ( isscalar(t) && t>=0 ) || \
-        usage("t must be nonnegative");
+        usage("t must be a scalar >= 0");
 
     ( 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
 
-    q = __ctmc_transient(Q, t, p0 );
-  endif
+    q = p0*expm(Q*t);
 
-endfunction
-
-## Helper function, compute steady state probability
-function q = __ctmc_steady_state( Q )
-  persistent epsilon = 10*eps;
-  N = rows(Q);
-
-  ## non zero columns
-  nonzero=find( any(abs(Q)>epsilon,1 ) );
-  if ( length(nonzero) == 0 )
-    error( "Q is the zero matrix" );
   endif
 
-  normcol = nonzero(1); # normalization condition column
-
-  ## force probability of unvisited states to zero
-  for i=find( all(abs(Q)<epsilon,1) )
-    Q(i,i) = 1;
-  endfor
-
-  ## assert( rank(Q) == N-1 );
-
-  Q(:,normcol) = 1; # add normalization condition
-  b = zeros(1,N); b(normcol)=1;
-  q = b/Q; # qQ = b;
 endfunction
 
-## Helper function, compute transient probability
-function q = __ctmc_transient( Q, t, p0 )
-  q = p0*expm(Q*t);
-endfunction
-
 %!test
 %! Q = [-1 1 0 0; 2 -3 1 0; 0 2 -3 1; 0 0 2 -2];
 %! q = ctmc(Q);

Modified: trunk/octave-forge/main/queueing/inst/ctmc_taexps.m
===================================================================
--- trunk/octave-forge/main/queueing/inst/ctmc_taexps.m	2012-04-04 12:47:24 UTC (rev 10145)
+++ trunk/octave-forge/main/queueing/inst/ctmc_taexps.m	2012-04-05 14:32:41 UTC (rev 10146)
@@ -75,7 +75,7 @@
   endif
 
   L = ctmc_exps(Q,varargin{:});
-  M = L ./ sum(L);
+  M = L ./ repmat(sum(L,2),1,columns(L));
 
 #{
   if ( nargin == 3 ) 
@@ -111,6 +111,7 @@
 %! for i=1:length(t)
 %!   M(i,:) = ctmc_taexps(Q,t(i),p);
 %! endfor
+%! clf;
 %! plot(t, M(:,1), ";State 1;", "linewidth", 2, \
 %!      t, M(:,2), ";State 2;", "linewidth", 2, \
 %!      t, M(:,3), ";State 3;", "linewidth", 2, \
@@ -141,7 +142,7 @@
 %!          0     0         0      d -d];
 %! p = ctmc(Q);
 %! printf("System availability: %f\n",p(1)+p(4));
-%! TT = linspace(1e-5,1*day,101);
+%! TT = linspace(0,1*day,101);
 %! PP = ctmc_taexps(Q,TT,[1 0 0 0 0]);
 %! A = At = Abart = zeros(size(TT));
 %! A(:) = p(1) + p(4); # steady-state availability
@@ -151,6 +152,7 @@
 %!   At(n) = p(1) + p(4); # instantaneous availability
 %!   Abart(n) = PP(n,1) + PP(n,4); # interval base availability
 %! endfor
+%! clf;
 %! semilogy(TT,A,";Steady-state;", \
 %!      TT,At,";Instantaneous;", \
 %!      TT,Abart,";Interval base;");

Modified: trunk/octave-forge/main/queueing/inst/dtmc.m
===================================================================
--- trunk/octave-forge/main/queueing/inst/dtmc.m	2012-04-04 12:47:24 UTC (rev 10145)
+++ trunk/octave-forge/main/queueing/inst/dtmc.m	2012-04-05 14:32:41 UTC (rev 10146)
@@ -86,9 +86,15 @@
   ( N>0 ) || \
       usage( err );
 
-  if ( nargin == 1 )
-    p = __dtmc_steady_state( P );  
-  else
+  if ( nargin == 1 ) # steady-state analysis
+    A = P-eye(N);
+    A(:,N) = 1; # add normalization condition
+    rank( A ) == N || \
+	warning( "dtmc(): P is reducible" );
+    
+    b = [ zeros(1,N-1) 1 ];
+    p = b/A;
+  else # transient analysis
     ( isscalar(n) && n>=0 ) || \
 	usage( "n must be >=0" );
 
@@ -97,27 +103,10 @@
 
     p0 = p0(:)'; # make p0 a row vector
 
-    p = __dtmc_transient(P, n, p0);
+    p = p0*P^n;
   endif
 endfunction
 
-## Helper function, compute steady-state probability
-function p = __dtmc_steady_state( P )
-  N = rows(P);
-  A = P-eye(N);
-  A(:,N) = 1; # add normalization condition
-  rank( A ) == N || \
-      warning( "dtmc(): P is reducible" );
-
-  b = [ zeros(1,N-1) 1 ];
-  p = b/A;
-endfunction
-
-## Helper function, compute transient probability
-function p = __dtmc_transient( P, n, p0 )
-  p = p0*P^n;
-endfunction
-
 %!test
 %! P = [0.75 0.25; 0.5 0.5];
 %! p = dtmc(P);

Modified: trunk/octave-forge/main/queueing/inst/dtmc_fpt.m
===================================================================
--- trunk/octave-forge/main/queueing/inst/dtmc_fpt.m	2012-04-04 12:47:24 UTC (rev 10145)
+++ trunk/octave-forge/main/queueing/inst/dtmc_fpt.m	2012-04-05 14:32:41 UTC (rev 10146)
@@ -140,15 +140,3 @@
 %! M = dtmc_fpt(P);
 %! assert( M(1:9 != 5,5)', [6 5 6 5 5 6 5 6], 10*eps );
 
-%!demo
-%! P = [ 0.0 0.9 0.1; \
-%!       0.1 0.0 0.9; \
-%!       0.9 0.1 0.0 ];
-%! M = dtmc_fpt(P);
-%! w = dtmc(P);
-%! N = rows(P);
-%! W = repmat(w,N,1);
-%! Z = inv(eye(N)-P+W);
-%! M1 = (repmat(diag(Z),1,N) - Z) ./ repmat(w',1,N);
-%! assert(M, M1);
-

Modified: trunk/octave-forge/main/queueing/inst/qnmmm.m
===================================================================
--- trunk/octave-forge/main/queueing/inst/qnmmm.m	2012-04-04 12:47:24 UTC (rev 10145)
+++ trunk/octave-forge/main/queueing/inst/qnmmm.m	2012-04-05 14:32:41 UTC (rev 10146)
@@ -141,7 +141,7 @@
   R = Q ./ X;
 endfunction
 %!demo
-%! disp("This is figure 6.4 on p. 220 Bolch et al.");
+%! # This is figure 6.4 on p. 220 Bolch et al.
 %! rho = 0.9;
 %! ntics = 21;
 %! lambda = 0.9;
@@ -153,6 +153,6 @@
 %! axis([0,ntics,0,25]);
 %! legend("Jobs in the system","Queue Length","location","northwest");
 %! xlabel("Number of servers (m)");
-%! title("\lambda = 0.9, \mu = 0.9");
+%! title("\\lambda = 0.9, \\mu = 0.9");
 
 

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