function [K] = svm_calc_kernel(x1,x2,parms)
%
% Calculate kernel between two sets of samples
%
%      kertype - on of {'none','poly','rbf_','givn'}
%      dim     -
%      gamma   - rbf width 
%      const   - for sigmoid or poly kernels 
%      scalar  - for sigmoid or poly kernels 
% 
  
 % Init 
  [parms,kertype]   = take_from_struct(parms,'kertype');
  [parms,dim]       = take_from_struct(parms,'dim', 2);
  [parms,const]     = take_from_struct(parms,'const', 0);  
  [parms,scalar]    = take_from_struct(parms,'scalar', 1);  
  [parms,gamma]     = take_from_struct(parms,'gamma', 1);
  [parms,fill_mode] = take_from_struct(parms,'fill_mode');  
  [parms,geoms1]    = take_from_struct(parms,'geoms1', []);  % geom scalars
  [parms,geoms2]    = take_from_struct(parms,'geoms2', []);  % geom scalars
  [n1,d1] = size(x1);
  [n2,d2] = size(x2);
  if(d1~=d2), error('Inconsistent dimensions');end;

  x1 = full(x1);
  x2 = full(x2);  
  
  % Handle missing values (NaN)
  switch(fill_mode)  
   case 'avgw',     
    % Skip NaN during Kernel calculation, correct kernel by fixing
    % the samples    
    xx1=~isnan(x1);
    xx2=~isnan(x2);    
    pk = mean([xx1;xx2]);
    pksqrt = sqrt(pk);
    
    if strcmp(kertype(1:4),'poly')
      pksqrt = pksqrt.^(1.0/dim);
    elseif strcmp(kertype(1:4),'line')
	% Do nothing 
    else
      error('unsuported kernel');
    end
        
    for i=1:n1
      x1(i,:) = x1(i,:) ./pksqrt;
    end
    for i=1:n2
      x2(i,:) = x2(i,:) ./pksqrt;
    end    
    
   case 'geom',
    % Skip NaNs during Kernel calculation, fix by geoms later
    if(isempty(geoms1)),geoms1 = ones(n1,1);end;
    if(isempty(geoms2)),geoms2 = ones(n2,1);end;    
    
    if(length(geoms1) ~= size(x1,1)), 
      size(geoms1),error('Inconsistent dimensions');end;
    if(length(geoms2) ~= size(x2,1)), 
      size(geoms2),error('Inconsistent dimensions');end;
      
   case 'zero',    
    x1(isnan(x1))=0;
    x2(isnan(x2))=0;    
    
   case 'none',
    % Do nothing (skip NaN's in kernel function)    
    
   otherwise
    % Do nothing (filled before) 
    
  end
  
  % Calc Kernel
  fprintf('Calculate kernel (%d by %d) ',n1,n2);
  kertype(end+1:4) = '_';
  switch(kertype)
   case 'line', 
    fprintf('Linear kernel %s \n',kertype);
    K=ker_line(x1,x2);

   case 'rbf_', 
    % fprintf('     gamma=%d\n',gamma);
    K=ker_rbf(x1,x2,gamma); 

   case 'poly', 
    fprintf('Polynomial kernel:');
    fprintf(' dim=%d,',dim);
    fprintf(' const=%d,',const);
    fprintf(' scalar=%d\n',scalar);
    K=ker_poly(x1,x2,dim,const,scalar);
    
   case 'sigm', 
    fprintf('Sigmoidal kernel:');    
    fprintf(' const=%d,',const); 
    fprintf(' scalar=%d\n',scalar);
    K=ker_sigmoid(x1,x2,const,scalar);

   otherwise 
    error('Unknown kernel type "%s"\n',kertype);
  end    
  
  % Fix kernel 
  switch(fill_mode)     
   case 'geom', K = fix_kernel_geom(K,geoms1,geoms2);
   otherwise,  fprintf('No Kernel fixing\n');
  end  
    
  return
end
% ===================== poly kernel ====
function K=ker_poly(x1,x2,dim,const,scalar)
%
% In kernel poly, settgin the NaNs to zero is eqiuvlent to skipping
% them, but faster to calculate
%
  
  x1(isnan(x1))=0;
  x2(isnan(x2))=0;      
  K = x1*x2';
  K = (scalar*K+const).^dim;
  
  return
end
% ===================== linear kernel =
function K=ker_line(x1,x2)
  
  x1(isnan(x1))=0;
  x2(isnan(x2))=0;      
  K = x1*x2';
    
  return
end
% ===================== rbf kernel ====
function K=ker_rbf(x1,x2,gamma)
  n1 = size(x1,1);  
  n2 = size(x2,1);    
  K = zeros(n1,n2);

  for i1=1:n1, inds1{i1} = find(~isnan(x1(i1,:))); end  %#ok
  for i2=1:n2, inds2{i2} = find(~isnan(x2(i2,:))); end  %#ok

  for i=1:n1
    fprintf('.');
    if(mod(i,10)==0), fprintf('%d',i);end;    
    for j=1:n2
      conct_inds = sort([inds1{i},inds2{j}]);
      inds = conct_inds(diff(conct_inds)==0);
      d = (x1(i,inds)-x2(j,inds));
      K(i,j) = exp(-gamma*sqrt(d*d'));
    end
  end
  
  return
end
% ===========================================
function K = fix_kernel_geom(K,geoms1,geoms2)
%
  
  % Check argin 
  [n1,n2] = size(K);  
  if(length(geoms1) ~= n1), error('Inconsistent dimensions');end;
  if(length(geoms2) ~= n2), error('Inconsistent dimensions');end;

  factor = geoms1(:) * geoms2(:)';
  K = K ./factor;
  
  return
end

  
