function f = fil_fac(s,reg_param,method,s1,V1)
%FIL_FAC Filter factors for some regularization methods.
%
% f = fil_fac(s,reg_param,method)
% f = fil_fac(sm,reg_param,method)  ,  sm = [sigma,mu]
% f = fil_fac(s1,k,'ttls',s1,V1)
%
% Computes all the filter factors corresponding to the
% singular values in s and the regularization parameter
% reg_param, for the following methods:
%    method = 'dsvd' : damped SVD
%    method = 'tsvd' : truncated SVD
%    method = 'Tikh' : Tikhonov regularization
%    method = 'ttls' : truncated TLS.
% If sm = [sigma,mu] is specified, then the filter factors
% for the corresponding generalized methods are computed.
%
% If method = 'ttls' then the singular values s1 and the
% right singular matrix V1 of [A,b] must also be supplied.
%
% If method is not specified, 'Tikh' is default.

% Per Christian Hansen, UNI-C, 06/23/93.

% Initialization.
[p,ps] = size(s); lr = length(reg_param);
if (nargin==2), method = 'Tikh'; end
f = zeros(p,lr);

% Check input data.
if (min(reg_param) <= 0)
  error('Regularization parameter must be positive')
end
if (method ~= 'Tikh' & min(reg_param) > p)
  error('Truncation parameter too large')
end

% Compute the filter factors.
for j=1:lr
  if (method(1:2)=='cg' | method(1:2)=='nu' | method(1:2)=='ls')
    error('Filter factors for iterative methods are not supported')
  elseif (method(1:4)=='dsvd')
    if (ps==1)
      f(:,j) = s./(s + reg_param(j));
    else
      f(:,j) = s(:,1)./(s(:,1) + reg_param(j)*s(:,2));
    end
  elseif (method(1:4)=='Tikh' | method(1:4)=='tikh')
    if (ps==1)
      f(:,j) = (s.^2)./(s.^2 + reg_param(j)^2);
    else
      f(:,j) = (s(:,1).^2)./(s(:,1).^2 + reg_param(j)^2*s(:,2).^2);
    end
  elseif (method(1:4)=='tsvd' | method(1:4)=='tgsv')
    if (ps==1)
      f(:,j) = [ones(reg_param(j),1);zeros(p-reg_param(j),1)];
    else
      f(:,j) = [zeros(p-reg_param(j),1);ones(reg_param(j),1)];
    end
  elseif (method(1:4)=='ttls')
    if (ps==1)
      coef = ((V1(p+1,:).^2)')/norm(V1(p+1,reg_param(j)+1:p+1))^2;
      for i=1:p
        k = reg_param(j);
        f(i,j) = s(i)^2*...
          sum( coef(k+1:p+1)./(s(i)+s1(k+1:p+1))./(s(i)-s1(k+1:p+1)) );
        if (f(i,j) < 0), f(i,j) = eps; end
        if (i > 1)
          if (f(i-1,j) <= eps & f(i,j) > f(i-1,j)), f(i,j) = f(i-1,j); end
        end
      end
    else
      error('The SVD of [A,b] must be supplied')
    end
  elseif (method(1:4)=='mtsv')
    error('Filter factors for MTSVD are not supported')
  else
    error('Illegal method')
  end
end
