% KVIS file kvis.m  Time to complete k visits to a set of states
% Version of 1/2/96 Revised 7/31/97 for version 4.2 and 5.1
% Calculates repeatedly the time to complete
% visits to k of the states in a prescribed set.
% Default is visit to all the target states.
% Assumes the procedure chainset has been run.
r  = input('Enter the number of repetitions  ');
disp('The target state set is:')
disp(e)
ks = input('Enter the number of target states to visit  ');
if isempty(ks)
  ks = ne;
end
if ks > ne
  ks = ne;
end
st = input('Enter the initial state  ');
if ~isempty(st)
  s1 = MS(st==states); % Initial state number
else
  s1 = 1;
end
disp(' ')
clear T                    % Trajectory in state numbers (reset)
R0 = zeros(1,ms);          % Indicator for target state numbers
R0(E) = ones(1,ne);        % reset
S = zeros(1,r);            % Terminal transitions for each rep (reset)
for k = 1:r
  R = R0;
  s  = s1;
  if R(s) == 1
    R(s) = 0;
  end
  i  = 1;
  T(1) = s;
  while sum(R) > ne - ks
    u = rand(1,1);
    s = ((A(s,:) < u)&(u <= B(s,:)))*MS';
    if R(s) == 1
      R(s) = 0;
    end
    i = i+1;
    T(i) = s;
  end
  S(k) = i-1;
end
if r == 1
  disp(['The time for completion is ',int2str(i-1),])
  N = 0:i-1;
  TR = [N;states(T)]';
  disp('To view the trajectory of states, call for TR')
else
  [t,f]  = csort(S,ones(1,r));
  p  = f/r;
  D  = [t;f]';
  AV = dot(t,p);
  SD = sqrt(dot(t.^2,p) - AV^2);
  MN = min(t);
  MX = max(t);
  disp(['The average completion time is ',num2str(AV),])
  disp(['The standard deviation is ',num2str(SD),])
  disp(['The minimum completion time is ',int2str(MN),])
  disp(['The maximum completion time is ',int2str(MX),])
  disp(' ')
  disp('To view a detailed count, call for D.')
  disp('The first column shows the various completion times;')
  disp('the second column shows the numbers of trials yielding those times')
end