File: qncsvisits.m

package info (click to toggle)
octave-queueing 1.2.8-3
  • links: PTS, VCS
  • area: main
  • in suites: forky, sid
  • size: 2,288 kB
  • sloc: makefile: 56
file content (132 lines) | stat: -rw-r--r-- 3,637 bytes parent folder | download | duplicates (2)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
## Copyright (C) 2012, 2016, 2020 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 with @math{K} service centers.
##
## @strong{INPUTS}
##
## @table @code
##
## @item @var{P}(i,j)
## probability that a request which completed service at center
## @math{i} is routed to center @math{j} (@math{K \times K} matrix).
## For closed networks it must hold that @code{sum(@var{P},2)==1}. The
## routing graph must be strongly connected, meaning that each node
## must be reachable from every other node.
##
## @item @var{r}
## Index of the reference station, @math{r \in @{1, @dots{}, K@}};
## Default @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 @code
##
## @item @var{V}(k)
## average number of visits to service center @math{k}, assuming
## @math{r} as the reference station.
##
## @end table
##
## @end deftypefn

## Author: Moreno Marzolla <moreno.marzolla(at)unibo.it>
## Web: http://www.moreno.marzolla.name/

function V = qncsvisits( P, r )

  if ( nargin < 1 || nargin > 2 )
    print_usage();
  endif

  issquare(P) || ...
      error("P must be a square matrix");

  [res err] = dtmcchkP(P);
  (res>0) || ...
      error( "invalid transition probability matrix P" );

  K = rows(P);

  if ( nargin < 2 )
    r = 1;
  else
    isscalar(r) || ...
	error("r must be a scalar");

    (r>=1 && r<=K) || ...
	error("r must be an integer in the range [1, %d]",K);
  endif

  V = zeros(size(P));
  A = P-eye(K);
  b = zeros(1,K);
  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)", "invalid" );
%! P = [1 0; 0.5 0];
%! fail( "qncsvisits(P)", "invalid" );
%! P = [1 2 3; 1 2 3];
%! fail( "qncsvisits(P)", "square" );
%! P = [0 1; 1 0];
%! fail( "qncsvisits(P,0)", "range" );
%! fail( "qncsvisits(P,3)", "range" );

%!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)