--- a/main/queueing/doc/queueingnetworks.txi	Wed Dec 05 22:49:14 2012 +0000
+++ b/main/queueing/doc/queueingnetworks.txi	Thu Dec 06 13:06:43 2012 +0000
@@ -235,7 +235,7 @@
 @iftex
 @tex
 $$\left\{\eqalign{V_j & = \sum_{i=1}^K V_i P_{i, j} \quad j=1, \ldots, K \cr
-                  V_1 & = 1}\right. $$
+                  V_r & = 1 \quad \hbox{for a selected reference station $r$}}\right. $$
 @end tex
 @end iftex
 @ifnottex
@@ -248,16 +248,33 @@
 |       /___
 |        i=1
 |
-| V_1 = 1
+| V_r = 1                     for a selected reference station r
 \
 @end example
 @end ifnottex
 
-@GETHELP{qnvisits}
+Note that the set of traffic equations @math{V_j = \sum_{i=1}^K V_i
+P_{i, j}} alone can only be solved up to a multiplicative constant; to
+get a unique solution we need an additional constraint @math{V_r =
+1}. This is equivalent to define station @math{r} as the
+@emph{reference station} (the default is @math{r=1},
+@xref{doc-qncsvisits}). The choice of the reference station for closed
+networks is extremely important, since it determines the performance
+results of the whole network.  First, the network throughput is set to
+the throughput of the reference station. Second, when a job returns to
+the reference station, it is assumed that the job completed its
+activity cycle, and enters a new activity cycle. Therefore, the system
+response time is also measured with respect to the reference station.
 
+@GETHELP{qncsvisits}
+
+@GETHELP{qnosvisits}
+
+@c
+@c
+@c
 @noindent @strong{EXAMPLE}
 
-
 @float Figure,fig:qn_closed_single
 @center @image{qn_closed_single,3in}
 @caption{Closed network with a single class of requests}
@@ -292,10 +309,11 @@
 @end example
 @end ifnottex
 
-The visit ratios @math{V} can be computed with:
+The visit ratios @math{V} using station 1 as the reference
+station can be computed with:
 
 @example
-@GETDEMO{qnvisits,1}
+@GETDEMO{qncsvisits,1}
    @result{} V = 1.00000   0.30000   0.70000
 @end example
 
@@ -337,12 +355,12 @@
 @math{V} as follows:
 
 @example
-@GETDEMO{qnvisits,2}
+@GETDEMO{qnosvisits,1}
    @result{} V = 5.0000   1.5000   2.5000
 @end example
 
-Function @command{qnvisits} expects a vector with @math{K} elements as a
-second parameter, for open networks only. The vector contains the
+Function @command{qnosvisits} expects a vector with @math{K} elements
+as a second parameter, for open networks only. The vector contains the
 arrival rates at each individual node; since in our example external
 arrivals exist only for node @math{S_1} with rate @math{\lambda =
 1.2}, the second parameter is @code{[1.2, 0, 0]}.
@@ -829,17 +847,25 @@
 
 @iftex
 @tex
-$$\left\{\eqalign{ V_{s, j} & = \sum_{r=1}^C \sum_{i=1}^K V_{r, i} P_{r, i, s, j}, \quad s=1, \ldots, K, j=1, \ldots, K \cr
-                   V_{s, 1} & = 1 \quad s=1, \ldots, K}\right. $$
+$$\left\{\eqalign{ V_{s, j} & = \sum_{r=1}^C \sum_{i=1}^K V_{r, i} P_{r, i, s, j}, \quad s=1, \ldots, C, j=1, \ldots, K \cr
+                   V_{s, r_s} & = 1 \quad s=1, \ldots, C}\right. $$
 @end tex
 @end iftex
 @ifnottex
 @group
-V_sj = sum_r sum_i V_ri P_risj    s=1,...,K, j=1,...,K
+V_sj = sum_r sum_i V_ri P_risj    s=1,...,C, j=1,...,K
+V_s r_s = 1                       s=1,...,C
 @end group
 @end ifnottex
 
-@noindent while for open networks:
+@noindent where @math{r_s} is the index of class @math{s}
+reference station; similarly to single class models, the traffic
+equation for closed multiclass networks can be solved up to
+multiplicative constants unless we choose one referenc station per
+class and set their visit ratios to 1. When multiple classes belong to
+the same chain, then only a single constraint per chain must be given.
+
+For open networks the traffic equations are as follows:
 
 @iftex
 @tex
@@ -859,6 +885,10 @@
 \sum_s \sum_j \lambda_{s, j}} is the overall external arrival rate to
 the whole system, then @math{P_{0, s, j} = \lambda_{s, j} / \lambda}.
 
+@GETHELP{qncmvisits}
+
+@GETHELP{qnomvisits}
+
 @c
 @c Open Networks
 @c
@@ -1269,8 +1299,8 @@
 @end table
 
 We can compute @math{V_k} from the routing probability matrix
-@math{P_{i, j}} using the @command{qnvisits} function
-@pxref{doc-qnvisits}.  We can analyze the network for a given
+@math{P_{i, j}} using the @command{qncsvisits} function
+@pxref{doc-qncsvisits}.  We can analyze the network for a given
 population size @math{N} (for example, @math{N=10}) as follows:
 
 @example
@@ -1278,7 +1308,7 @@
 @kbd{N = 10;}
 @kbd{S = [1 2 0.8];}
 @kbd{P = [0 0.3 0.7; 1 0 0; 1 0 0];}
-@kbd{V = qnvisits(P);}
+@kbd{V = qncsvisits(P);}
 @kbd{[U R Q X] = qnclosed( N, S, V )}
    @result{} U = 0.99139 0.59483 0.55518
    @result{} R = 7.4360  4.7531  1.7500
@@ -1361,14 +1391,14 @@
 arrivals at centers 2 and 3.
 
 Similarly to closed networks, we first need to compute the visit
-counts @math{V_k} to center @math{k}. Again, we use the
-@command{qnvisits} function as follows:
+counts @math{V_k} to center @math{k}. We use the
+@command{qnosvisits} function as follows:
 
 @example
 @group
 @kbd{P = [0 0.3 0.5; 1 0 0; 1 0 0];}
 @kbd{lambda = [0.15 0 0];}
-@kbd{V = qnvisits(P, lambda)}
+@kbd{V = qnosvisits(P, lambda)}
    @result{} V = 5.00000 1.50000 2.50000
 @end group
 @end example

--- a/main/queueing/inst/dtmcisir.m	Wed Dec 05 22:49:14 2012 +0000
+++ b/main/queueing/inst/dtmcisir.m	Thu Dec 06 13:06:43 2012 +0000
@@ -91,7 +91,7 @@
 %! assert( max(s), 1 );
 %! assert( min(s), 1 );
 
-## FIXME: (mosty) copied from qnvisits.m; use a better algorithm for SCC
+## FIXME: (mosty) copied from qncmvisits.m; use a better algorithm for SCC
 ## (e.g.,
 ## http://pmtksupport.googlecode.com/svn/trunk/gaimc1.0-graphAlgo/scomponents.m
 function s = __scc(G)
@@ -111,7 +111,7 @@
   endfor
 endfunction
 
-## FIXME: (mosty) copied from qnvisits.m
+## FIXME: (mosty) copied from qncmvisits.m
 function v = __dfs(G, s)
   assert( issquare(G) );
   N = rows(G);

--- a/main/queueing/inst/qnclosed.m	Wed Dec 05 22:49:14 2012 +0000
+++ b/main/queueing/inst/qnclosed.m	Thu Dec 06 13:06:43 2012 +0000
@@ -77,7 +77,7 @@
 %! m = ones(size(S));             # All centers are single-server
 %! Z = 2;                         # External delay
 %! N = 15;                        # Maximum population to consider
-%! V = qnvisits(P); # Compute number of visits from P
+%! V = qncsvisits(P);             # Compute number of visits
 %! X_bsb_lower = X_bsb_upper = X_ab_lower = X_ab_upper = X_mva = zeros(1,N);
 %! for n=1:N
 %!   [X_bsb_lower(n) X_bsb_upper(n)] = qncsbsb(n, S, V, m, Z);

--- a/main/queueing/inst/qnclosedgb.m	Wed Dec 05 22:49:14 2012 +0000
+++ b/main/queueing/inst/qnclosedgb.m	Thu Dec 06 13:06:43 2012 +0000
@@ -85,7 +85,7 @@
 %! P = [0 0.3 0.7; 1 0 0; 1 0 0];
 %! S = [1 0.6 0.2];
 %! m = ones(1,3);
-%! V = qnvisits(P);
+%! V = qncsvisits(P);
 %! Nmax = 20;
 %! tol = 1e-5; # compensate for numerical errors
 %! ## Test case with Z>0
@@ -103,7 +103,7 @@
 %! P = [0 0.3 0.7; 1 0 0; 1 0 0];
 %! S = [1 0.6 0.2];
 %! m = ones(1,3);
-%! V = qnvisits(P);
+%! V = qncsvisits(P);
 %! Nmax = 20;
 %! tol = 1e-5; # compensate for numerical errors
 %!

--- a/main/queueing/inst/qncmmva.m	Wed Dec 05 22:49:14 2012 +0000
+++ b/main/queueing/inst/qncmmva.m	Thu Dec 06 13:06:43 2012 +0000
@@ -99,7 +99,8 @@
 ## @item P
 ## @code{@var{P}(r,i,s,j)} is the probability that a class @math{r}
 ## job completing service at center @math{i} is routed to center @math{j}
-## as a class @math{s} job. @strong{If you pass this argument,
+## as a class @math{s} job; the reference station for all jobs
+## is considered station 1. @strong{If you pass argument @var{P},
 ## class switching is allowed}, but you can not specify any external delay
 ## (i.e., @var{Z} must be zero).
 ##
@@ -228,7 +229,7 @@
   U = R = Q = X = zeros(C,K);
 
   ## 1. Compute visit counts
-  [V ch] = qnvisits(P);
+  [V ch] = qncmvisits(P);
 
   ## 2. Identify chains
   nch = max(ch);
@@ -634,7 +635,7 @@
 %! N = [40 4];
 %! S = [ 5.0 .010 .035 .035 .035 .035; \
 %!      10.0 .100 .035 .035 .035 .035 ];
-%! V = qnvisits(P);
+%! V = qncmvisits(P);
 %! [U R Q X] = qncmmva(N, S, V, [-1 1 1 1 1 1]);
 %! # FIXME: The results below were computed with JMVA; the numbers
 %! # in the paper are different (wrong?!?)!!
@@ -680,7 +681,7 @@
 %!
 %! # Compute visits
 %!
-%! V = qnvisits(P);
+%! V = qncmvisits(P);
 %!
 %! # Define population and service times
 %! 
@@ -735,7 +736,6 @@
 %! P(1,2,1,1) = P(2,2,2,1) = 1;
 %! P(1,3,1,1) = P(2,3,2,1) = 1;
 %! N = [3 0];
-%! V = qnvisits(P);
 %! [U R Q X] = qncmmva(N, S, P);
 %! assert( R, [.015 .133 .163; .073 .133 .163], 1e-3 );
 %! assert( X, [12.609 8.826 2.522; 6.304 1.891 3.152], 1e-3 );

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/main/queueing/inst/qncmvisits.m	Thu Dec 06 13:06:43 2012 +0000
@@ -0,0 +1,349 @@
+## Copyright (C) 2008, 2009, 2010, 2011, 2012 Moreno Marzolla
+##
+## This file is part of the queueing toolbox.
+##
+## The queueing toolbox is free software: you can redistribute it and/or
+## modify it under the terms of the GNU General Public License as
+## published by the Free Software Foundation, either version 3 of the
+## License, or (at your option) any later version.
+##
+## The queueing toolbox is distributed in the hope that it will be
+## useful, but WITHOUT ANY WARRANTY; without even the implied warranty
+## of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+## General Public License for more details.
+##
+## You should have received a copy of the GNU General Public License
+## along with the queueing toolbox. If not, see <http://www.gnu.org/licenses/>.
+
+## -*- texinfo -*-
+##
+## @deftypefn {Function File} {[@var{V} @var{ch}] =} qncmvisits (@var{P})
+## @deftypefnx {Function File} {[@var{V} @var{ch}] =} qncmvisits (@var{P}, @var{r})
+##
+## Compute the average number of visits to the service centers of a closed multiclass network.
+##
+## @strong{INPUTS}
+##
+## @table @var
+##
+## @item P
+## @code{@var{P}(r,i,s,j)} is the probability that a
+## class @math{r} request which completed service at center @math{i} is
+## routed to center @math{j} as a class @math{s} request. Class switching
+## is allowed.
+##
+## @item r
+## @code{@var{r}(c)} is the index of class @math{c} reference station.
+## The class @math{c} visit count to server @code{@var{r}(c)} is
+## conventionally set to 1. The reference station serves two purposes:
+## its throughput is assumed to be the system throughput, and a job
+## returning to the reference station is assumed to have completed one
+## cycle. Default is station 1 for all classes.
+##
+## @end table
+##
+## @strong{OUTPUTS}
+##
+## @table @var
+##
+## @item V
+## @code{@var{V}(c,i)} is the number of visits of class @math{c}
+## requests at center @math{i}. Note that @code{@var{V}(c,@var{r}(c)) =
+## 1}.
+##
+## @item ch
+## @code{@var{ch}(c)} is the chain number that class @math{c} belongs
+## to. Different classes can belong to the same chain. Chains are
+## numbered sequentially starting from 1 (@math{1, 2, @dots{}}). The
+## total number of chains is @code{max(@var{ch})}.
+##
+## @end table
+## 
+## @end deftypefn
+
+## Author: Moreno Marzolla <marzolla(at)cs.unibo.it>
+## Web: http://www.moreno.marzolla.name/
+
+function [V chains] = qncmvisits( P, r )
+
+  if ( nargin < 1 || nargin > 2 )
+    print_usage();
+  endif
+
+  ndims(P) == 4 || \
+      error("P must be a 4-dimensional matrix");
+
+  [C, K, C2, K2] = size( P );
+  (K == K2 && C == C2) || \
+      error( "P must be a [C,K,C,K] matrix");
+
+  if ( nargin < 2)
+    r = ones(1,C);
+  else
+    isvector(r) && length(r) == C || \
+	error("r must be a vector with %d elements",C);
+    all( r>=1 && r<=K ) || \
+	error("elements of r are out of range [1,%d]",K);
+    r = r(:)';
+  endif
+
+  ## solve the traffic equations: V(s,j) = sum_r sum_i V(r,i) *
+  ## P(r,i,s,j), for all s,j V(s,r(s)) = 1 for all s.
+  A = reshape(P,[K*C K*C])-eye(K*C);
+  b = zeros(1,K*C);
+
+  CH = __scc(reshape(P,[C*K C*K])>0);
+  nCH = max(CH); # number of chains
+  CH = reshape(CH,C,K); # CH(c,k) is the chain that class c at center k belongs to
+
+  chains = zeros(1,C);
+  
+  for k=1:K
+    for c=1:C
+      if ( chains(c) == 0 )
+        chains(c) = CH(c,k);
+      else
+        ( CH(c,k) == 0 || chains(c) == CH(c,k) ) || \
+            error("Class %d belongs to different chains",c);
+      endif
+    endfor
+  endfor
+
+  constraints = zeros(1,nCH); # constraint(cc) = 1 iff we set a constraint for a class belonging to chain cc; we only set one constraint per chain
+
+  for c=1:C
+    cc = CH(c,r(c));
+    if ( cc == 0 || constraints(cc) == 0 ) 
+      ii = sub2ind([C K],c,r(c));
+      A(:,ii) = 0;
+      A(ii,ii) = 1;
+      if ( cc > 0 ) ## if r(c) is not an isolated node
+	constraints(cc) = 1;
+	b(ii) = 1;
+      endif
+    endif
+  endfor
+
+  V = reshape(b/A, C, K);
+  ## Make sure that no negative values appear (sometimes, numerical
+  ## errors produce tiny negative values instead of zeros)
+  V = max(0,V);
+endfunction
+
+%!test
+%!
+%! ## Closed, multiclass network
+%!
+%! C = 2; K = 3; 
+%! P = zeros(C,K,C,K);
+%! P(1,1,1,2) = 1;
+%! P(1,2,1,1) = 1;
+%! P(2,1,2,3) = 1;
+%! P(2,3,2,1) = 1;
+%! V = qncmvisits(P);
+%! for c=1:C
+%!   for k=1:K
+%!     assert(V(c,k), sum(sum(V .* P(:,:,c,k))), 1e-5);
+%!   endfor
+%! endfor
+
+%!test
+%!
+%! ## Test multiclass network. Example from Schwetman (figure 7, page 9 of
+%! ## http://docs.lib.purdue.edu/cgi/viewcontent.cgi?article=1258&context=cstech
+%! ## "Testing network-of-queues software, technical report CSD-TR 330,
+%! ## Purdue University).
+%!
+%! C = 2; K = 4;
+%! P = zeros(C,K,C,K);
+%! # class 1 routing
+%! P(1,1,1,1) = .05;
+%! P(1,1,1,2) = .45;
+%! P(1,1,1,3) = .5;
+%! P(1,2,1,1) = 1;
+%! P(1,3,1,1) = 1;
+%! # class 2 routing
+%! P(2,1,2,1) = .01;
+%! P(2,1,2,3) = .5;
+%! P(2,1,2,4) = .49;
+%! P(2,3,2,1) = 1;
+%! P(2,4,2,1) = 1;
+%! V = qncmvisits(P);
+%! for c=1:C
+%!   for i=1:K
+%!     assert(V(c,i), sum(sum(V .* P(:,:,c,i))), 1e-5);
+%!   endfor
+%! endfor
+
+%!test
+%!
+%! ## Network with class switching.
+%! ## This is example in figure 9 of
+%! ## Schwetman, "Implementing the Mean Value Analysis
+%! ## Algorithm fort the solution of Queueing Network Models", Technical
+%! ## Report CSD-TR-355, Department of Computer Science, Purdue Univrsity,
+%! ## Feb 15, 1982
+%! ## http://docs.lib.purdue.edu/cgi/viewcontent.cgi?article=1285&context=cstech
+%!
+%! C = 2; K = 3;
+%! S = [.01 .07 .10; \
+%!      .05 0.7 .10 ];
+%! P = zeros(C,K,C,K);
+%! P(1,1,1,2) = .7;
+%! P(1,1,1,3) = .2;
+%! P(1,1,2,1) = .1;
+%! P(2,1,2,2) = .3;
+%! P(2,1,2,3) = .5;
+%! P(2,1,1,1) = .2;
+%! P(1,2,1,1) = P(1,3,1,1) = 1;
+%! P(2,2,2,1) = P(2,3,2,1) = 1;
+%! N = [3 0];
+%! V = qncmvisits(P);
+%! VV = [10 7 2; 5 1.5 2.5]; # result given in Schwetman; our function computes something different, but that's ok since visit counts are actually ratios
+%! assert( V ./ repmat(V(:,1),1,K), VV ./ repmat(VV(:,1),1,K), 1e-5 );
+
+%!test
+%!
+%! ## two disjoint classes: must produce two disjoing chains
+%!
+%! C = 2; K = 3;
+%! P = zeros(C,K,C,K);
+%! P(1,1,1,2) = 1;
+%! P(1,2,1,1) = 1;
+%! P(2,1,2,3) = 1;
+%! P(2,3,2,1) = 1;
+%! [nc r] = qncmvisits(P);
+%! assert( r(1) != r(2) );
+
+%!test
+%!
+%! ## two classes, one chain
+%!
+%! C = 2; K = 3;
+%! P = zeros(C,K,C,K);
+%! P(1,1,1,2) = .5;
+%! P(1,2,2,1) = 1;
+%! P(2,1,2,3) = .5;
+%! P(2,3,1,1) = 1;
+%! [nc r] = qncmvisits(P);
+%! assert( r(1) == r(2) );
+
+%!test
+%! 
+%! ## a "Moebius strip". Note that this configuration is invalid, and
+%! ## therefore our algorithm must raise an error. This is because this
+%! ## network has two chains, but both chains contain both classes
+%!
+%! C = 2; K = 2;
+%! P = zeros(C,K,C,K);
+%! P(1,1,2,2) = 1;
+%! P(2,2,1,1) = 1;
+%! P(2,1,1,2) = 1;
+%! P(1,2,2,1) = 1;
+%! fail( "qncmvisits(P)", "different");
+
+%!test
+%!
+%! ## Network with two classes representing independent chains.
+%! ## This is example in figure 8 of
+%! ## Schwetman, "Implementing the Mean Value Analysis
+%! ## Algorithm fort the solution of Queueing Network Models", Technical
+%! ## Report CSD-TR-355, Department of Computer Science, Purdue Univrsity,
+%! ## Feb 15, 1982
+%! ## http://docs.lib.purdue.edu/cgi/viewcontent.cgi?article=1285&context=cstech
+%! 
+%! C = 2; K = 2;
+%! P = zeros(C,K,C,K);
+%! P(1,1,1,3) = P(1,3,1,1) = 1;
+%! P(2,2,2,3) = P(2,3,2,2) = 1;
+%! V = qncmvisits(P,[1,2]);
+%! assert( V, [1 0 1; 0 1 1], 1e-5 );
+
+%!test
+%! C = 2;
+%! K = 3;
+%! P = zeros(C,K,C,K);
+%! P(1,1,1,2) = 1;
+%! P(1,2,1,3) = 1;
+%! P(1,3,2,2) = 1;
+%! P(2,2,1,1) = 1;
+%! [V ch] = qncmvisits(P);
+%! assert( ch, [1 1] );
+
+## The following transition probability matrix is not well formed: note
+## that there is an outgoing transition from center 1, class 1 but not
+## incoming transition.
+%!test
+%! C = 2;
+%! K = 3;
+%! P = zeros(C,K,C,K);
+%! P(1,1,1,2) = 1;
+%! P(1,2,1,3) = 1;
+%! P(1,3,2,2) = 1;
+%! P(2,2,2,1) = 1;
+%! P(2,1,1,2) = 1;
+%! [V ch] = qncmvisits(P);
+%! assert( ch, [1 1] );
+
+%!demo
+%! P = [ 0 0.4 0.6 0; \
+%!       0.2 0 0.2 0.6; \
+%!       0 0 0 1; \
+%!       0 0 0 0 ];
+%! lambda = [0.1 0 0 0.3];
+%! V = qncmvisits(P,lambda);
+%! S = [2 1 2 1.8];
+%! m = [3 1 1 2];
+%! [U R Q X] = qnos( sum(lambda), S, V, m );
+
+## compute strongly connected components using Kosaraju's algorithm,
+## which requires two DFS visits. A better solution would be to use
+## Tarjan's algorithm.
+##
+## In this implementation, an isolated node without self loops will NOT
+## belong to any SCC. Although this is not formally correct from the
+## graph theoretic point of view, it is necessary to compute chains
+## correctly.
+function s = __scc(G)
+  assert(issquare(G));
+  N = rows(G);
+  GF = (G>0);
+  GB = (G'>0);
+  s = zeros(N,1);
+  c=1;
+  for n=1:N
+    if (s(n) == 0)
+      fw = __dfs(GF,n);
+      bw = __dfs(GB,n);
+      r = (fw & bw);
+      if (any(r))
+	s(r) = c++;
+      endif
+    endif
+  endfor
+endfunction
+
+## Executes a dfs visit on graph G, starting from source node s
+function v = __dfs(G, s)
+  assert( issquare(G) );
+  N = rows(G);
+  v = stack = zeros(1,N); ## v(i) == 1 iff node i has been visited
+  q = 1; # first empty slot in queue
+  stack(q++) = s;
+  ## Note: node s is NOT marked as visited; it will me marked as visited
+  ## only if we visit it again. This is necessary to ensure that
+  ## isolated nodes without self loops will not belong to any SCC.
+  while( q>1 )
+    n = stack(--q);
+    ## explore neighbors of n: all f in G(n,:) such that v(f) == 0
+    
+    ## The following instruction is equivalent to:
+    ##    for f=find(G(n,:))
+    ##      if ( v(f) == 0 )
+    for f = find ( G(n,:) & (v==0) )
+      stack(q++) = f;
+      v(f) = 1;
+    endfor
+  endwhile
+endfunction
+

--- a/main/queueing/inst/qncsgb.m	Wed Dec 05 22:49:14 2012 +0000
+++ b/main/queueing/inst/qncsgb.m	Thu Dec 06 13:06:43 2012 +0000
@@ -230,7 +230,7 @@
 %! P = [0 0.3 0.7; 1 0 0; 1 0 0];
 %! S = [1 0.6 0.2];
 %! m = ones(1,3);
-%! V = qnvisits(P);
+%! V = qncsvisits(P);
 %! Z = 2;
 %! Nmax = 20;
 %! tol = 1e-5; # compensate for numerical errors
@@ -248,7 +248,7 @@
 %!test
 %! P = [0 0.3 0.7; 1 0 0; 1 0 0];
 %! S = [1 0.6 0.2];
-%! V = qnvisits(P);
+%! V = qncsvisits(P);
 %! Nmax = 20;
 %! tol = 1e-5; # compensate for numerical errors
 %!

--- a/main/queueing/inst/qncsmvablo.m	Wed Dec 05 22:49:14 2012 +0000
+++ b/main/queueing/inst/qncsmvablo.m	Thu Dec 06 13:06:43 2012 +0000
@@ -125,7 +125,7 @@
   z = ones(1,N);
   lambda = 0;
   ## Computation of the visit counts
-  v = qnvisits(P);
+  v = qncsvisits(P);
   D = S .* v; # Service demand
   ## Main loop
   for k=1:K
@@ -193,7 +193,7 @@
 %! P = [0 0.7 0.3; 1 0 0; 1 0 0];
 %! K = 40;
 %! [U1 R1 Q1] = qncsmvablo( K, S, M, P );
-%! V = qnvisits(P);
+%! V = qncsvisits(P);
 %! [U2 R2 Q2] = qncsmva( K, S, V );
 %! assert( U1, U2, 1e-5 );
 %! assert( R1, R2, 1e-5 );

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/main/queueing/inst/qncsvisits.m	Thu Dec 06 13:06:43 2012 +0000
@@ -0,0 +1,131 @@
+## Copyright (C) 2012 Moreno Marzolla
+##
+## This file is part of the queueing toolbox.
+##
+## The queueing toolbox is free software: you can redistribute it and/or
+## modify it under the terms of the GNU General Public License as
+## published by the Free Software Foundation, either version 3 of the
+## License, or (at your option) any later version.
+##
+## The queueing toolbox is distributed in the hope that it will be
+## useful, but WITHOUT ANY WARRANTY; without even the implied warranty
+## of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+## General Public License for more details.
+##
+## You should have received a copy of the GNU General Public License
+## along with the queueing toolbox. If not, see <http://www.gnu.org/licenses/>.
+
+## -*- texinfo -*-
+##
+## @deftypefn {Function File} {@var{V} =} qncsvisits (@var{P})
+## @deftypefnx {Function File} {@var{V} =} qncsvisits (@var{P}, @var{r})
+##
+## Compute the mean number of visits to the service centers of a
+## single class, closed network.
+##
+## @strong{INPUTS}
+##
+## @table @var
+##
+## @item P
+## Routing probability matrix.
+## @code{@var{P}(i,j)} is the probability that a request which completed
+## service at center @math{i} is routed to center @math{j}. For closed
+## networks it must hold that @code{sum(@var{P},2)==1}. The routing
+## graph myst be strongly connected, meaning that it must be possible to
+## eventually reach each node starting from each node. 
+##
+## @item r
+## Index of the reference station; if omitted, it is assumed
+## @code{@var{r}=1}. The traffic equations are solved by imposing the
+## condition @code{@var{V}(r) = 1}. A request returning to the reference
+## station completes its activity cycle.
+##
+## @end table
+##
+## @strong{OUTPUTS}
+##
+## @table @var
+##
+## @item V
+## @code{@var{V}(i)} is the average number of visits to service center
+## @math{i}, assuming center @math{r} as the reference station.
+##
+## @end table
+## 
+## @end deftypefn
+
+## Author: Moreno Marzolla <marzolla(at)cs.unibo.it>
+## Web: http://www.moreno.marzolla.name/
+
+function V = qncsvisits( P, r )
+
+  if ( nargin < 1 || nargin > 2 )
+    print_usage();
+  endif
+
+  ndims(P) == 2 && issquare(P) || \
+      error("P must be a 2-dimensional square matrix");
+
+  [res err] = dtmcchkP(P);
+  (res>0) || \
+      error( "P is not a transition probability matrix for closed networks" );
+
+  N = rows(P);
+
+  if ( nargin < 2 )
+    r = 1;
+  else
+    issclarar(r) || \
+	error("r must be a scalar");
+
+    (r>=1 && r<=N) || \
+	error("r must be an integer in [%d,%d]",1,N);
+  endif
+
+  V = zeros(size(P));
+  A = P-eye(N);
+  b = zeros(1,N);
+  A(:,r) = 0; A(r,r) = 1;
+  b(r) = 1;
+  V = b/A;
+
+  ## Make sure that no negative values appear (sometimes, numerical
+  ## errors produce tiny negative values instead of zeros)
+  V = max(0,V);
+
+endfunction
+%!test
+%! P = [-1 0; 0 0];
+%! fail( "qncsvisits(P)", "not a transition probability" );
+%! P = [1 0; 0.5 0];
+%! fail( "qncsvisits(P)", "not a transition probability" );
+
+%!test
+%!
+%! ## Closed, single class network
+%!
+%! P = [0 0.3 0.7; 1 0 0; 1 0 0];
+%! V = qncsvisits(P);
+%! assert( V*P,V,1e-5 );
+%! assert( V, [1 0.3 0.7], 1e-5 );
+
+%!test
+%!
+%! ## Test tolerance of the qncsvisits() function. 
+%! ## This test builds transition probability matrices and tries
+%! ## to compute the visit counts on them. 
+%!
+%! for k=[5, 10, 20, 50]
+%!   P = reshape(1:k^2, k, k);
+%!   P = P ./ repmat(sum(P,2),1,k);
+%!   V = qncsvisits(P);
+%!   assert( V*P, V, 1e-5 );
+%! endfor
+
+%!demo
+%! P = [0 0.3 0.7; \
+%!      1 0   0  ; \
+%!      1 0   0  ];
+%! V = qncsvisits(P)
+

--- a/main/queueing/inst/qnjackson.m	Wed Dec 05 22:49:14 2012 +0000
+++ b/main/queueing/inst/qnjackson.m	Thu Dec 06 13:06:43 2012 +0000
@@ -145,7 +145,7 @@
   endif
 
   ## Compute the arrival rates using the traffic equation:
-  l = sum(lambda)*qnvisits( P, lambda );
+  l = sum(lambda)*qnosvisits( P, lambda );
   ## Check ergodicity
   for i=1:N
     if ( m(i)>0 && l(i)>=m(i)/S(i) )      

--- a/main/queueing/inst/qnmarkov.m	Wed Dec 05 22:49:14 2012 +0000
+++ b/main/queueing/inst/qnmarkov.m	Thu Dec 06 13:06:43 2012 +0000
@@ -328,7 +328,7 @@
 %! C = [6 6 6];
 %! S = [1 0.8 1.8];
 %! N = 6;
-%! [U1 R1 Q1 X1] = qnclosed( N, S, qnvisits(P) );
+%! [U1 R1 Q1 X1] = qnclosed( N, S, qncsvisits(P) );
 %! [U2 R2 Q2 X2] = qnmarkov( N, S, C, P );
 %! assert( U1, U2, 1e-6 );
 %! assert( R1, R2, 1e-6 );

--- a/main/queueing/inst/qnmvablo.m	Wed Dec 05 22:49:14 2012 +0000
+++ b/main/queueing/inst/qnmvablo.m	Thu Dec 06 13:06:43 2012 +0000
@@ -82,7 +82,7 @@
 %! P = [0 0.7 0.3; 1 0 0; 1 0 0];
 %! K = 40;
 %! [U1 R1 Q1] = qnmvablo( K, S, M, P );
-%! V = qnvisits(P);
+%! V = qncsvisits(P);
 %! [U2 R2 Q2] = qnclosedsinglemva( K, S, V );
 %! assert( U1, U2, 1e-5 );
 %! assert( R1, R2, 1e-5 );

--- a/main/queueing/inst/qnom.m	Wed Dec 05 22:49:14 2012 +0000
+++ b/main/queueing/inst/qnom.m	Thu Dec 06 13:06:43 2012 +0000
@@ -79,7 +79,7 @@
 ##
 ## @end table
 ##
-## @seealso{qnopen,qnos,qnvisits}
+## @seealso{qnopen,qnos,qnomvisits}
 ##
 ## @end deftypefn
 

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/main/queueing/inst/qnomvisits.m	Thu Dec 06 13:06:43 2012 +0000
@@ -0,0 +1,112 @@
+## Copyright (C) 2012 Moreno Marzolla
+##
+## This file is part of the queueing toolbox.
+##
+## The queueing toolbox is free software: you can redistribute it and/or
+## modify it under the terms of the GNU General Public License as
+## published by the Free Software Foundation, either version 3 of the
+## License, or (at your option) any later version.
+##
+## The queueing toolbox is distributed in the hope that it will be
+## useful, but WITHOUT ANY WARRANTY; without even the implied warranty
+## of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+## General Public License for more details.
+##
+## You should have received a copy of the GNU General Public License
+## along with the queueing toolbox. If not, see <http://www.gnu.org/licenses/>.
+
+## -*- texinfo -*-
+##
+## @deftypefn {Function File} {@var{V} =} qnomvisits (@var{P}, @var{lambda})
+##
+## Compute the average number of visits to the service centers of an open multiclass network.
+##
+## @strong{INPUTS}
+##
+## @table @var
+##
+## @item P
+## @code{@var{P}(r,i,s,j)} is the probability that a
+## class @math{r} request which completed service at center @math{i} is
+## routed to center @math{j} as a class @math{s} request. Class switching
+## is supported.
+##
+## @item lambda
+## @code{@var{lambda}(r,i)} is the arrival rate of class @math{r}
+## requests to center @math{i}.
+##
+## @end table
+##
+## @strong{OUTPUTS}
+##
+## @table @var
+##
+## @item V
+## @code{@var{V}(r,i)} is the number of visits of class @math{r}
+## requests at center @math{i}.
+##
+## @end table
+## 
+## @end deftypefn
+
+## Author: Moreno Marzolla <marzolla(at)cs.unibo.it>
+## Web: http://www.moreno.marzolla.name/
+
+function V = qnomvisits( P, lambda )
+
+  if ( nargin != 2 )
+    print_usage();
+  endif
+
+  ndims(P) == 4 || \
+      error("P must be a 4-dimensional matrix");
+
+  [C, K, C2, K2] = size( P );
+  (K == K2 && C == C2) || \
+      error( "P must be a [C,K,C,K] matrix");
+
+  [C,K] == size(lambda) || \
+      error( "lambda must be a %d x %d matrix", C, K );
+
+  all(lambda(:)>=0) || \
+      error(" lambda contains negative values" );
+
+  ## solve the traffic equations: V(s,j) = lambda(s,j) / lambda + sum_r
+  ## sum_i V(r,i) * P(r,i,s,j), for all s,j where lambda is defined as
+  ## sum_r sum_i lambda(r,i)
+  A = eye(K*C) - reshape(P,[K*C K*C]);
+  b = reshape(lambda / sum(lambda(:)), [1,K*C]);
+  V = reshape(b/A, [C, K]);
+
+  ## Make sure that no negative values appear (sometimes, numerical
+  ## errors produce tiny negative values instead of zeros)
+  V = max(0,V);
+endfunction
+%!test
+%! fail( "qnomvisits( zeros(3,3,3), [1 1 1] )", "matrix");
+
+%!test
+%! C = 2; K = 4;
+%! P = zeros(C,K,C,K);
+%! # class 1 routing
+%! P(1,1,1,1) = .05;
+%! P(1,1,1,2) = .45;
+%! P(1,1,1,3) = .5;
+%! P(1,2,1,1) = 0.1;
+%! P(1,3,1,1) = 0.2;
+%! # class 2 routing
+%! P(2,1,2,1) = .01;
+%! P(2,1,2,3) = .5;
+%! P(2,1,2,4) = .49;
+%! P(2,3,2,1) = 0.2;
+%! P(2,4,2,1) = 0.16;
+%! lambda = [0.1 0 0 0.1 ; 0 0 0.2 0.1];
+%! lambda_sum = sum(lambda(:));
+%! V = qnomvisits(P, lambda);
+%! assert( all(V(:)>=0) );
+%! for i=1:K
+%!   for c=1:C
+%!     assert(V(c,i), lambda(c,i) / lambda_sum + sum(sum(V .* P(:,:,c,i))), 1e-5);
+%!   endfor
+%! endfor
+

--- a/main/queueing/inst/qnopensingle.m	Wed Dec 05 22:49:14 2012 +0000
+++ b/main/queueing/inst/qnopensingle.m	Thu Dec 06 13:06:43 2012 +0000
@@ -82,7 +82,7 @@
 %!test
 %! lambda=[1];
 %! P=[0];
-%! V=qnvisits(P,lambda);
+%! V=qnosvisits(P,lambda);
 %! S=[0.25];
 %! [U1 R1 Q1 X1]=qnopensingle(sum(lambda),S,V); 
 %! [U2 R2 Q2 X2]=qnmm1(lambda(1),1/S(1));

--- a/main/queueing/inst/qnos.m	Wed Dec 05 22:49:14 2012 +0000
+++ b/main/queueing/inst/qnos.m	Thu Dec 06 13:06:43 2012 +0000
@@ -88,7 +88,7 @@
 ##
 ## @end table
 ##
-## @seealso{qnopen,qnclosed,qnvisits}
+## @seealso{qnopen,qnclosed,qnosvisits}
 ##
 ## @end deftypefn
 
@@ -189,7 +189,7 @@
 %!test
 %! lambda=[1];
 %! P=[0];
-%! V=qnvisits(P,lambda);
+%! V=qnosvisits(P,lambda);
 %! S=[0.25];
 %! [U1 R1 Q1 X1]=qnos(sum(lambda),S,V); 
 %! [U2 R2 Q2 X2]=qsmm1(lambda(1),1/S(1));

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/main/queueing/inst/qnosvisits.m	Thu Dec 06 13:06:43 2012 +0000
@@ -0,0 +1,117 @@
+## Copyright (C) 2012 Moreno Marzolla
+##
+## This file is part of the queueing toolbox.
+##
+## The queueing toolbox is free software: you can redistribute it and/or
+## modify it under the terms of the GNU General Public License as
+## published by the Free Software Foundation, either version 3 of the
+## License, or (at your option) any later version.
+##
+## The queueing toolbox is distributed in the hope that it will be
+## useful, but WITHOUT ANY WARRANTY; without even the implied warranty
+## of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+## General Public License for more details.
+##
+## You should have received a copy of the GNU General Public License
+## along with the queueing toolbox. If not, see <http://www.gnu.org/licenses/>.
+
+## -*- texinfo -*-
+##
+## @deftypefn {Function File} {@var{V} =} qnosvisits (@var{P}, @var{lambda})
+##
+## Compute the average number of visits to the service centers of a single 
+## class open Queueing Network with @math{N} service centers.
+##
+## @strong{INPUTS}
+##
+## @table @var
+##
+## @item P
+## @code{@var{P}(i,j)} is the probability that a request which completed
+## service at center @math{i} is routed to center @math{j}. 
+##
+## @item lambda
+## @code{@var{lambda}(i)} is the external arrival rate to
+## center @math{i}. 
+##
+## @end table
+##
+## @strong{OUTPUTS}
+##
+## @table @var
+##
+## @item V
+## @code{@var{V}(i)} is the average number of
+## visits to server @math{i}.
+##
+## @end table
+## 
+## @end deftypefn
+
+## Author: Moreno Marzolla <marzolla(at)cs.unibo.it>
+## Web: http://www.moreno.marzolla.name/
+
+function V = qnosvisits( P, lambda )
+
+  if ( nargin != 2 )
+    print_usage();
+  endif
+
+  ndims(P) == 2 && issquare(P) || \
+      error("P must be a square matrix");
+
+  N = size(P,1);
+  V = zeros(N,N);
+
+  ##
+  ## Open network
+  ##
+  all(P(:)>=0) && all( sum(P,2)<=1+1e-5 ) || \
+      error( "P is not a transition probability matrix for open networks" );
+  
+  ( isvector(lambda) && length(lambda) == N ) || \
+      error( "lambda must be a vector with %d elements", N );
+  all( lambda>= 0 ) || \
+      error( "lambda contains negative values" );
+
+  lambda = lambda(:)';
+  
+  A = eye(N)-P;
+  b = lambda / sum(lambda);
+  V = b/A;
+  ## Make sure that no negative values appear (sometimes, numerical
+  ## errors produce tiny negative values instead of zeros)
+  V = max(0,V);
+endfunction
+%!test
+%! P = [1 0; 0 1]; lambda=[0 -1];
+%! fail( "qnosvisits(P,lambda)", "contains negative" );
+
+%!test
+%!
+%! ## Open, single class network
+%!
+%! P = [0 0.2 0.5; 1 0 0; 1 0 0];
+%! lambda = [ 0.1 0.3 0.2 ];
+%! V = qnosvisits(P,lambda);
+%! assert( V*P+lambda/sum(lambda),V,1e-5 );
+
+%!demo
+%! p = 0.3;
+%! lambda = 1.2
+%! P = [0 0.3 0.5; \
+%!      1 0   0  ; \
+%!      1 0   0  ];
+%! V = qnosvisits(P,[1.2 0 0])
+
+%!demo
+%! P = [ 0 0.4 0.6 0; \
+%!       0.2 0 0.2 0.6; \
+%!       0 0 0 1; \
+%!       0 0 0 0 ];
+%! lambda = [0.1 0 0 0.3];
+%! V = qnosvisits(P,lambda);
+%! S = [2 1 2 1.8];
+%! m = [3 1 1 2];
+%! [U R Q X] = qnos( sum(lambda), S, V, m );
+

--- a/main/queueing/inst/qnsolve.m	Wed Dec 05 22:49:14 2012 +0000
+++ b/main/queueing/inst/qnsolve.m	Thu Dec 06 13:06:43 2012 +0000
@@ -754,7 +754,7 @@
 %!        qnmknode( "m/m/m-fcfs", 0.05 ) };
 %! P = [ 0 0.5 0.5 0; 1 0 0 0; 0.6 0 0 0; 1 0 0 0 ];
 %! lambda = [0 0 0 4];
-%! [U R Q X] = qnsolve("open", sum(lambda), QQ, qnvisits(P,lambda) );
+%! [U R Q X] = qnsolve("open", sum(lambda), QQ, qnosvisits(P,lambda) );
 %! assert( X, [20 10 10 4], 1e-4 );
 %! assert( U, [0.8 0.3 0.6 0.2], 1e-2 );
 %! assert( R, [0.2 0.043 0.15 0.0625], 1e-3 );

--- a/main/queueing/inst/qnvisits.m	Wed Dec 05 22:49:14 2012 +0000
+++ b/main/queueing/inst/qnvisits.m	Thu Dec 06 13:06:43 2012 +0000
@@ -43,8 +43,8 @@
 ## networks, @code{@var{lambda}(i)} is the external arrival rate to
 ## center @math{i}. For multiple class networks,
 ## @code{@var{lambda}(r,i)} is the arrival rate of class @math{r}
-## requests to center @math{i}. If this parameter is omitted, the
-## network is assumed to be closed.
+## requests to center @math{i}. If this function is called with a single
+## argument, the network is assumed to be closed.
 ##
 ## @end table
 ##
@@ -54,9 +54,10 @@
 ##
 ## @item V
 ## For single class networks, @code{@var{V}(i)} is the average number of
-## visits to server @math{i}. For multiple class networks,
-## @code{@var{V}(r,i)} is the class @math{r} visit ratio at center
-## @math{i}.
+## visits to server @math{i}, assuming center 1 as the reference station
+## (i.e., @code{@var{V}(1) = 1}). For multiple class networks,
+## @code{@var{V}(r,i)} is the number of visits of class @math{r}
+## requests at center @math{i}.
 ##
 ## @item ch
 ## (For closed networks only). @code{@var{ch}(c)} is the chain number
@@ -74,6 +75,13 @@
 
 function [V ch] = qnvisits( P, varargin )
 
+  persistent warned = false;
+  if (!warned)
+    warned = true;
+    warning("qn:deprecated-function",
+	    "qnvisits is deprecated. Please use one of qncsvisits, qnosvisits, qncmvisits, qnomvisits instead");
+  endif
+
   if ( nargin < 1 || nargin > 2 )
     print_usage();
   endif