% cost_update( col, v, vl) - update FIR and G for Markov time-varying ISI
%
function cost_update( col, v, vl)

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

global FIR_M;
global FIR_M_t;
global FIR_M_p_;

global G;
global DD;
global qq;
global alpha;
global alpha2;

global B;
global avg_h;
global inv_Rh;
global inv_Rh_d;
global var_noise;
global z;


% disp('cost_update'), col

if( col < N)

  c1 = col + 1;

  if( isempty(B) ) % called from x4_em_init

    FIR = FIR_save(:,c1);
    trans = FSM_trans( states, FIR);

  elseif( ~isempty(G))

    if( size(avg_h,2) == 1)	% Gauss-Markov

      Fv = [1:FIR_length_p];

      % update current G and q
      %
      if( length(var_noise) > 1)
        vn = var_noise(col);
      else
        vn = var_noise;
      end

      t = B(col,:)';
      D = FIR_M_t*t*((t'/vn)*FIR_M) + inv_Rh;
      q = FIR_M_t*t*((z(col) - t' * FIR_M_p_)/vn) + inv_Rh_d;

      if( col > 1)
        CDk1 = inv_Rh * DD(:,Fv+(col-2)*FIR_length_p);
        D = D + alpha2*(inv_Rh - CDk1*inv_Rh);
        q = q - alpha*inv_Rh_d + alpha*(CDk1*qq(:,col-1));
      end

      D = inv(D);
      DD(:,Fv+(col-1)*FIR_length_p) = D;
      qq(:,col) = q;

      % set next G and q
      %
      if( length(var_noise) > 1)
        vn = var_noise(c1);
      else
        vn = var_noise;
      end

      t = B(c1,:)';
      CDk1 = inv_Rh * D;
      q = FIR_M_t*t*((z(c1) - t' * FIR_M_p_)/vn) + inv_Rh_d + alpha*(CDk1*q); % RHS D and q are from previous section
      D = FIR_M_t*t*((t'/vn)*FIR_M) + inv_Rh - alpha2*(CDk1*inv_Rh);

      if( c1 < N)
        CDk1 = inv_Rh * DD(:,Fv+c1*FIR_length_p);
        D = D + alpha2*(inv_Rh - CDk1*inv_Rh);
        q = q + alpha*(CDk1*qq(:,c1+1) - inv_Rh_d);
      end

      D = inv(D);
      d = D*q;
      G = FIR_M * D * FIR_M_t;
      FIR = FIR_M * d + FIR_M_p_;

    else  % time-varying Gaussian

      if( length(var_noise) > 1)
        vn = var_noise(c1);
      else
        vn = var_noise;
      end

      t = B(c1,:)';
      D = inv( FIR_M_t*t*((t'/vn)*FIR_M) + inv_Rh);
      G = FIR_M * D * FIR_M_t;
      d = D * (FIR_M_t*t*(z(c1) - t' * FIR_M_p_)/vn + inv_Rh*avg_h(:,c1));
      FIR = FIR_M * d + FIR_M_p_;

    end % size(avg_h) == 1

    trans = FSM_trans( states, FIR);

  end % ~isempty(G)

end % col < N

% FIR
