1## Copyright (C) 2012, 2016, 2020 Moreno Marzolla
2##
3## This file is part of the queueing toolbox.
4##
5## The queueing toolbox is free software: you can redistribute it and/or
6## modify it under the terms of the GNU General Public License as
7## published by the Free Software Foundation, either version 3 of the
8## License, or (at your option) any later version.
9##
10## The queueing toolbox is distributed in the hope that it will be
11## useful, but WITHOUT ANY WARRANTY; without even the implied warranty
12## of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
13## General Public License for more details.
14##
15## You should have received a copy of the GNU General Public License
16## along with the queueing toolbox. If not, see <http://www.gnu.org/licenses/>.
17
18## -*- texinfo -*-
19##
20## @deftypefn {Function File} {@var{V} =} qncsvisits (@var{P})
21## @deftypefnx {Function File} {@var{V} =} qncsvisits (@var{P}, @var{r})
22##
23## Compute the mean number of visits to the service centers of a
24## single class, closed network with @math{K} service centers.
25##
26## @strong{INPUTS}
27##
28## @table @code
29##
30## @item @var{P}(i,j)
31## probability that a request which completed service at center
32## @math{i} is routed to center @math{j} (@math{K \times K} matrix).
33## For closed networks it must hold that @code{sum(@var{P},2)==1}. The
34## routing graph must be strongly connected, meaning that each node
35## must be reachable from every other node.
36##
37## @item @var{r}
38## Index of the reference station, @math{r \in @{1, @dots{}, K@}};
39## Default @code{@var{r}=1}. The traffic equations are solved by
40## imposing the condition @code{@var{V}(r) = 1}. A request returning to
41## the reference station completes its activity cycle.
42##
43## @end table
44##
45## @strong{OUTPUTS}
46##
47## @table @code
48##
49## @item @var{V}(k)
50## average number of visits to service center @math{k}, assuming
51## @math{r} as the reference station.
52##
53## @end table
54##
55## @end deftypefn
56
57## Author: Moreno Marzolla <moreno.marzolla(at)unibo.it>
58## Web: http://www.moreno.marzolla.name/
59
60function V = qncsvisits( P, r )
61
62  if ( nargin < 1 || nargin > 2 )
63    print_usage();
64  endif
65
66  issquare(P) || ...
67      error("P must be a square matrix");
68
69  [res err] = dtmcchkP(P);
70  (res>0) || ...
71      error( "invalid transition probability matrix P" );
72
73  K = rows(P);
74
75  if ( nargin < 2 )
76    r = 1;
77  else
78    isscalar(r) || ...
79	error("r must be a scalar");
80
81    (r>=1 && r<=K) || ...
82	error("r must be an integer in the range [1, %d]",K);
83  endif
84
85  V = zeros(size(P));
86  A = P-eye(K);
87  b = zeros(1,K);
88  A(:,r) = 0; A(r,r) = 1;
89  b(r) = 1;
90  V = b/A;
91  ## Make sure that no negative values appear (sometimes, numerical
92  ## errors produce tiny negative values instead of zeros)
93  V = max(0,V);
94endfunction
95%!test
96%! P = [-1 0; 0 0];
97%! fail( "qncsvisits(P)", "invalid" );
98%! P = [1 0; 0.5 0];
99%! fail( "qncsvisits(P)", "invalid" );
100%! P = [1 2 3; 1 2 3];
101%! fail( "qncsvisits(P)", "square" );
102%! P = [0 1; 1 0];
103%! fail( "qncsvisits(P,0)", "range" );
104%! fail( "qncsvisits(P,3)", "range" );
105
106%!test
107%!
108%! ## Closed, single class network
109%!
110%! P = [0 0.3 0.7; 1 0 0; 1 0 0];
111%! V = qncsvisits(P);
112%! assert( V*P,V,1e-5 );
113%! assert( V, [1 0.3 0.7], 1e-5 );
114
115%!test
116%!
117%! ## Test tolerance of the qncsvisits() function.
118%! ## This test builds transition probability matrices and tries
119%! ## to compute the visit counts on them.
120%!
121%! for k=[5, 10, 20, 50]
122%!   P = reshape(1:k^2, k, k);
123%!   P = P ./ repmat(sum(P,2),1,k);
124%!   V = qncsvisits(P);
125%!   assert( V*P, V, 1e-5 );
126%! endfor
127
128%!demo
129%! P = [0 0.3 0.7; ...
130%!      1 0   0  ; ...
131%!      1 0   0  ];
132%! V = qncsvisits(P)
133
134