% cost_init() - initialize time-varying Gauss or Gauss-Markov ISI FIR and G
%
function cost_init()

global KNOWN;
global DET;
global GAUSS;
global GM;
global ML;
global KALEH;

global GAUSS_H;
global GM_H;

global trans;
global states;
global FIR;
global FIR_save;
global FIR_length;
global Memory;

global G;
global avg_h;

global GG;
global qq;
global alpha;
global alpha2;

global N;
global B;
global Rh;
global inv_Rh;
global inv_Rh_d;
global var_noise;
global z;

% disp('cost_init');

% Note: inv_Rh here is inv(C) in the paper


% known channel or direct ML
%
if( KNOWN | ML)

  % FIR isn't used in direct ML, but trans is needed
  %
  FIR = FIR_save(:,1);
  trans = FSM_trans( states, FIR);
  return;

end


% Kalet-Vallet
%
if( KALEH)

  % FIR was set by kaleh_iter
  %
  trans = FSM_trans( states, FIR);
  return;

end


% EM initalization, no symbol estimates
%
if( isempty(B) )

  if( DET)

    % [0 ... 0 1 0 ... 0]
    %
    % FIR = zeros(FIR_length,1);
    % FIR( (size(FIR,1)+1)/2) = 1;

    %% using average avg_h
    %
    % if( GM_H)
    %   FIR = avg_h/(1-alpha);    % Gauss-Markov
    % else
    %   FIR = (sum(avg_h')/N)';   % Gaussian
    % end

    %% using FIR_save(:,1)
    %
     FIR = FIR_save(:,1);
    
  elseif( GAUSS)

    % disp( 'cost_init: isempty(B) and GAUSS');

    % reasonable initialization
    %
     FIR = avg_h(:,1);
     G = Rh;

    % junk initialization, just for testing
    %
    % original, unstable for high SNR:
    %
    % t = 2*(rand(FIR_length,1) > 0.5) - 1; % random [-1,+1]
    % G = inv( t*(t'/var_noise) + inv_Rh);
    % FIR = G*( t*(z(1)/var_noise) + inv_Rh*avg_h(:,1));
    %
    % new stable way, 26 June 2000
    %
    % t = 2*(rand(FIR_length,1) > 0.5) - 1; % random [-1,+1]
    % s2 = var_noise + t'*(Rh*t);
    % G = Rh - (Rh*t)*((t'*Rh)/s2);
    % FIR = avg_h(:,1) + Rh*(t*((z(1)-t'*avg_h(:,1))/s2));

  elseif( GM)

    %% using regular formulas with B=0 is equivalent to:
    %%
    % FIR = avg_h/(1-alpha);
    % G = Rh/(1-alpha2);

    %% [0 ... 0 1 0 ... 0]
    %%
    % FIR = zeros(FIR_length,1);
    % FIR( (size(FIR,1)+1)/2) = 1;
    % G = Rh/(1-alpha2);

    %% using FIR_save(:,1)
    %
     FIR = FIR_save(:,1);
     G = Rh/(1-alpha2);

  else

    error( 'cost_init type flag not set.');

  end

% have symbol estimates in B
%
else 

  if( DET)

    % with get_A() fix, 9 April 2000, B can be singular,
    % so use pinv()
    %
    % old:
    %       FIR = B\(z');
    %
    FIR = pinv(B)*(z');

  elseif( GAUSS)

    % disp( 'cost_init: ~isempty(B) and GAUSS');

    % original, unstable for high SNR:
    %
    % t = B(1,:)';
    % G = inv( t*(t'/var_noise) + inv_Rh);
    % FIR = G*( t*(z(1)/var_noise) + inv_Rh*avg_h(:,1));

    % new stable way, 24 June 2000
    %
    t = B(1,:)';
    s2 = var_noise + t'*(Rh*t);
    G = Rh - (Rh*t)*((t'*Rh)/s2);
    FIR = avg_h(:,1) + Rh*(t*((z(1)-t'*avg_h(:,1))/s2));

  elseif( GM)

    GG = zeros( FIR_length, N*FIR_length);
    qq = zeros( FIR_length, N);
    Fv = [1:FIR_length];

    for k=N:-1:2

      t = B(k,:)';
      G = t*(t'/var_noise) + inv_Rh;
      q = t*(z(k)/var_noise) + inv_Rh_d;
      if( k < N)
        CGk1 = inv_Rh * GG(:,Fv+k*FIR_length);
        G = G + alpha2*(inv_Rh - CGk1*inv_Rh);
        q = q + alpha*(CGk1*qq(:,k+1) - inv_Rh_d);
      end
      G = inv(G);
      GG(:,Fv+(k-1)*FIR_length) = G;
      qq(:,k) = q;
    end

    if(1) % known g1/G1
      %% using FIR_save(:,1)
      %
       FIR = FIR_save(:,1);
       G = Rh/(1-alpha2);
    else
      % k = 1
      %
      t = B(1,:)';
      CGk1 = inv_Rh * G;
      q = t*(z(1)/var_noise) + inv_Rh_d + alpha*(CGk1*q);		% RHS G and q are from k=2
      G = t*(t'/var_noise) + inv_Rh - alpha2*(CGk1*inv_Rh);
      G = inv(G);
      FIR = G*q;
    end

  end

end

trans = FSM_trans( states, FIR);
