global kxv;

global tr_period;
global w;
global Wtot;
global MM_max;
global MM;
global Cmiss;
global Cfalse;
global missed_cost;
global false_alarm_rate;
global Wmin;
global PKZ;
global PKZ1;

global Winpath_index;
global Winel_index;
global Winperm_index;
global WinMM;

global WINDOW_DEBUG;

% Edited by Jessica Bradley 3/14/2000
% This function is the same as PKK.m except that it subsitutes tr_period for
% N, Winpath_index for path_index and Winel_index for el_index. This allows
% the program only to look at the specific window and estimate the number of
% tracks in it.
%
% Also, WinMM is the subset of MM for the current window
%
% ------------------------------------------------------------------------

% compute -log(P(k|Z)) for k=0..K using results of K-path Viterbi
%
% Wmin will be used as a scale factor to prevent underflow in P
%
% Note that for unreachable final states, exp(-HUGE_VAL) == 0,
% so these states do not affect the probability calculations here.
%
% P will include all final cost results for k=K, and all costs for
% unique paths from the K-path L-best for 0<k<K
%
% P1 will include just the single lowest costs.
%
P  = zeros( K, 1);
P1 = zeros( K, 1);

% P(K) and P1(K) come directly from K-path Viterbi results
%
W = sort(w(:));	% all costs, not just the L-best
Wmin = W(1);
P(K) = sum( exp( -(W - Wmin)));
P1(K) = exp( -(W(1) - Wmin));

K1 = [1:K];

% for 0<k<K, extract the unique tracks from the L-best to estimate P and P1
%
for k = 1:K-1

  k

  cc = combinations(K1,k);
  ncc = size(cc,2);		% number of combinations

  x = zeros(k*tr_period,ncc*L);
  c = zeros(ncc*L,1);

  % L-best paths taken k at a time, and associated Kalman costs
  %
  for i = 1:L
    [p t] = SortCols(pathsLocal(i));
    tc = Wtot(K1+(Winpath_index(i,tr_period)-1)*K,Winel_index(i,tr_period));
    tc = tc(t);
    offset = (i-1)*ncc;
    for j = 1:ncc
      q = j + offset;
      t = cc(:,j);
      r = p(:,t);
      x(:,q) = r(:);
      c(q) = sum(tc(t));
    end
  end

  if( WINDOW_DEBUG)
    cc, x, c
  end

  % indices and costs for the unique sets of tracks
  %
  u = Unique(x);
  uc = c(u);
  nc = size(uc,1);

  if( WINDOW_DEBUG)
    u, x(:,u)
  end

  for j = 1:nc
    xj = reshape( x(:,u(j)), tr_period, k);
    for i = 1:tr_period
      n = sum(xj(i,:) <= MM_max);			% number of measurements used for tracks
      m = WinMM(i) - n;					% number of false detections
% i, j, k, n, m
      uc(j) = uc(j) + Cmiss(n+1,k) + Cfalse(m+1);	% apriori costs
      if( i > 1)					% missed target costs only for i>1
        uc(j) = uc(j) + m*missed_cost;
      end
    end
  end

  % W = sort(uc);
  [W track_index] = sort(uc);

  if( WINDOW_DEBUG)
    W, track_index
  end

  xj = reshape( x(:,u(track_index(1))), tr_period, k);

  xj

  P(k) = sum( exp( -(W - Wmin)));
  P1(k) = exp( -(W(1) - Wmin));

end


% compute P(0|Z) manually, Viterbi not needed
%
% if( tr_period > 1)
%   W0 = missed_cost * sum(WinMM(2:tr_period));
% else
%   W0 = 0
% end
%
W0 = missed_cost * sum(WinMM);

for i=1:tr_period
  W0 = W0 - log( p_poisson( WinMM(i), false_alarm_rate) );
end

P0=exp(-(W0-Wmin));


% Compute sum over tau' and normalize.
%
% The normalization is done only for convenience and does not affect the
% analysis of which P(k|Z) is maximum.
%
PKZ = [ P0; P];

PKZ = -log(PKZ/sum(PKZ));

disp( 'K    PKZ');

[ [0:K]' PKZ]


% Compute PKZ1 using just the single lowest costs
%
PKZ1 = [P0; P1];

PKZ1 = -log(PKZ1/sum(PKZ1));

% disp( 'K    PKZ1');
% 
% [ [0:K]' PKZ1]


% try something crazy here:
%
% kxv=zeros(size(kxv));
