function outputs = ...
    classify_with_absent(trn_data, tst_data, trn_lbls, tst_lbls, parms)
%
% classify_with_absent features. Supports mulktiple fill modes
% 
% Inputs: 
%   trn_data - Each row is a sample, Missing values are denoted with NaNs.
%   trn_lbls - column vector
%   parms - options (structure)
%
% (C) Gal Chechik

  % Check inputs
  n_trn_samples = size(trn_data,1);
  n_tst_samples = size(tst_data,1);
  if n_trn_samples ~= size(trn_lbls,1), error('Inconsistent input size');end;
  if n_tst_samples ~= size(tst_lbls,1), error('Inconsistent input size');end;

  % Init 
  parms.dummy=0;
  default_gstr = '0.001,0.01,0.1,1,10,100,1000,10000,100000';
  default_cs_list =  10.^[-2:3];
  [parms,do_force]   = take_from_struct(parms,'do_force', 0);
  [parms,seed]       = take_from_struct(parms,'seed', 1);
  [parms,fill_mode]  = take_from_struct(parms,'fill_mode', 'geom');
  [parms,dig1]       = take_from_struct(parms,'dig1', 5);
  [parms,dig2]       = take_from_struct(parms,'dig2', 6);
  [parms,X]          = take_from_struct(parms,'X', 1);
  [parms,t]          = take_from_struct(parms,'t', 1);
  [parms,gstr]       = take_from_struct(parms,'gstr', default_gstr);
  [parms,n_neighbors]= take_from_struct(parms,'n_neighbors', 30);
  [parms,dbname]     = take_from_struct(parms,'dbname', 'mnist');
  [parms,cs_list]    = take_from_struct(parms,'cs_list',default_cs_list);
  [parms,cgstep]     = take_from_struct(parms,'cgstep',0.1);
  [parms,n_restarts] = take_from_struct(parms,'n_restarts',1);


  % Load results from disk (or calculate)
  rslt_filename = set_filenames('rslt',parms);
  vars = {'outputs'};
  [do_calc, outputs] = cond_load(rslt_filename, do_force, vars{1:end});
  if do_calc<1, return;  end   

  
  ncs = length(cs_list);
  rand('state',seed); randn('state',seed);
  outputs.dummy=[];

  % Fill missing values as a preprocessing phase (if needed)  
  [trn_data,tst_data] = preprocess_fill_missing(fill_mode,parms,trn_data, ...
						tst_data);
  trn_data = full(trn_data);
  tst_data = full(tst_data);
   

% Prepare labels and file names 
if(isempty(dig1)==0 & isempty(dig2)==0)
  trn_lbls(find(trn_lbls==dig1))=1;
  tst_lbls(find(tst_lbls==dig1))=1;
  trn_lbls(find(trn_lbls==dig2))= -1;
  tst_lbls(find(tst_lbls==dig2))= -1;
end

% Add constant and split to trn/tst
trn_data = [ones(size(trn_data,1),1) trn_data(:,:)];
tst_data = [ones(size(tst_data,1),1) tst_data(:,:)];

n_trn_actual = size(trn_data,1);
n_trn1 = ceil(n_trn_actual*0.8);

trn_data1 = trn_data(1       :n_trn1,:);
trn_data2 = trn_data(n_trn1+1:end   ,:);
trn_lbls1 = trn_lbls (1       :n_trn1,:);
trn_lbls2 = trn_lbls (n_trn1+1:end   ,:);

% Train svm using cplex
x_trn = trn_data1;
y_trn = trn_lbls1;
x_tst = trn_data2;
y_tst = trn_lbls2;
x_vld = tst_data;
y_vld = tst_lbls;


% Calculate kernels
K_trn = load_kernel(x_trn,x_trn,parms);
K_tst = load_kernel(x_trn,x_tst,parms);
K_vld = load_kernel(x_trn,x_vld,parms);

x_zro=x_trn;x_zro(find(isnan(x_zro)))=0;
% Train classifier (can loop over C values)
if((fill_mode(1:4) == 'geom') || ...
   (fill_mode(1:4) == 'soge') || ...          
   (fill_mode(1:4) == 'geop') || ...       
   (fill_mode(1:4) == 'invg') || ...    
   (fill_mode(1:4) == 'bfge'))
  n_iters = 7;
elseif(...
    fill_mode(1:4) == 'prgr'||...
    fill_mode(1:4) == 'qcqp')  
  n_iters = 20;
else
  n_iters = 1;
end

outputs.w         = zeros(ncs,n_iters*n_restarts,size(x_trn,2));
outputs.geoms_trn = zeros(ncs,n_iters*n_restarts,length(y_trn));    
outputs.geoms_tst = zeros(ncs,n_iters*n_restarts,length(y_tst));    
outputs.geoms_vld = zeros(ncs,n_iters*n_restarts,length(y_vld));    
outputs.alphas    = zeros(ncs,n_iters*n_restarts,length(y_trn));

% Loop over CS list
  fprintf('n_iters=%d n_restarts=%d\n',n_iters,n_restarts);  
for(i_c = 1:length(cs_list))
  parms.C       = cs_list(i_c);
  parms.verbose = 0;
  fprintf('train CPLEX SVM C=%4.2f. Iter ',parms.C);

  % Loop over restarts
  for(i_restart=1:n_restarts)  
  
    i_iter=0;
    % Loop over iterations
    for(ii_iter=1:n_iters)
      i_iter=i_iter+1;
      
      % Opimize S(i)
      if exist('w','var')==0, w=[];end;
      [geoms_trn, geoms_tst, geoms_vld] = ...
        optimize_geoms(fill_mode, w, x_trn, y_trn, x_tst, y_tst, ...
		       x_vld, y_vld, ii_iter, i_restart);

      fixed_K_trn = (geoms_trn' * geoms_trn).* K_trn;
      fixed_K_tst = (geoms_trn' * geoms_tst).* K_tst;
      fixed_K_vld = (geoms_trn' * geoms_vld).* K_vld;	      
      [signed_alphas,b,f,w] = svmtrain_cplex(x_zro, y_trn, fixed_K_trn, parms);
      
      
      % Calc errors 	  
      [y_prd,y_mrg] = svmclassify(fixed_K_trn,signed_alphas,b);
            
      trn_acc  = sum (y_trn == y_prd) / length(y_prd);
      [trn_auc,trn_tp,trn_fp,trn_sens,trn_spec,acc,threshes]  = roc(y_trn,y_mrg);
      [m,i] = max(acc); empb = threshes(i);
      
      [y_prd,y_mrg] = svmclassify(fixed_K_tst,signed_alphas,b);
      tst_acc  = sum (y_tst == y_prd)/ length(y_prd);
      [tst_auc,tst_tp,tst_fp,tst_sens,tst_spec,acc,threshes]  = roc(y_tst,y_mrg);
            
      [y_prd,y_mrg] = svmclassify(fixed_K_vld,signed_alphas,b);
      vld_acc  = sum (y_vld == y_prd)/ length(y_prd);
      [vld_auc,vld_tp,vld_fp,vld_sens,vld_spec]  = roc(y_vld,y_mrg);      

      % Output resutls 
      disp([trn_acc,tst_acc,vld_acc]);
      disp([trn_auc,tst_auc,vld_auc]);
      
      % Save results 
      outputs.alphas(i_c,i_iter,:) = signed_alphas;         
      outputs.b(i_c,i_iter)   = b;
      outputs.f(i_c,i_iter)   = f;
      outputs.w(i_c,i_iter,:) = w(:)';
      
      outputs.geoms_trn(i_c,i_iter,1:length(geoms_trn)) = geoms_trn(:)';
      outputs.geoms_tst(i_c,i_iter,1:length(geoms_tst)) = geoms_tst(:)';
      outputs.geoms_vld(i_c,i_iter,1:length(geoms_vld)) = geoms_vld(:)';
      outputs.acc(i_c,i_iter,1:3)  = [trn_acc,tst_acc,vld_acc];
      outputs.auc(i_c,i_iter,1:3)  = [trn_auc,tst_auc,vld_auc];
      
      outputs.trn_tp(i_c,i_iter,1:length(trn_tp)) = trn_tp;
      outputs.trn_fp(i_c,i_iter,1:length(trn_tp)) = trn_fp;
      outputs.tst_tp(i_c,i_iter,1:length(tst_tp)) = tst_tp;
      outputs.tst_fp(i_c,i_iter,1:length(tst_tp)) = tst_fp;
      outputs.vld_tp(i_c,i_iter,1:length(vld_tp)) = vld_tp;
      outputs.vld_fp(i_c,i_iter,1:length(vld_tp)) = vld_fp;
      
      outputs.trn_sens(i_c,i_iter,1:length(trn_sens)) = trn_sens;
      outputs.trn_spec(i_c,i_iter,1:length(trn_spec)) = trn_spec;
      outputs.tst_sens(i_c,i_iter,1:length(tst_sens)) = tst_sens;
      outputs.tst_spec(i_c,i_iter,1:length(tst_spec)) = tst_spec;
      outputs.vld_sens(i_c,i_iter,1:length(vld_sens)) = vld_sens;
      outputs.vld_spec(i_c,i_iter,1:length(vld_spec)) = vld_spec;
      
      outputs.vld_y_mrg(i_c,i_iter,1:length(y_mrg)) = y_mrg;        
      outputs.vld_y_lbl(i_c,i_iter,1:length(y_mrg)) = y_vld;
      
      fprintf('%d ',i_iter);
    end
    fprintf('\n');
    
  end % over restarts 
  % fprintf('\n');    
  
end % Loop over CS list 

% Save resutls to tmp file
save(rslt_filename,'outputs');

return
end

% ====================================================================
function  [trn_data,tst_data] = preprocess_fill_missing(fill_mode, ...
						  parms,trn_data, tst_data)
%
% fill missing values as an unsupervised step, depending on fill_mode
%

  fill_mode(end+1:4) = '_';
  switch(fill_mode(1:4)) 
   case 'mean',[trn_data, tst_data] = load_filled_mean(parms, trn_data, tst_data);    
   case 'flag',[trn_data, tst_data] = load_flagged_data(parms, trn_data, tst_data);
   case 'blfl',[trn_data, tst_data] = load_blflagged_data(parms, trn_data, tst_data);
   case 'none',% do nothing  Kernel skips NaNs
   case 'zero',% do nothing  later fill with zeros;  
   case 'avgw',% do nothing, later use mean norm fix
   case 'geom',% do nothing, later use geometric based margin  
   case 'soge',% do nothing, later use soft-geometric based margin    
   case 'geop',% do nothing, later use geometric based margin 
   case 'qcqp',% do nothing, 
   otherwise,  error('Illegal fill mode');
  end
  return
end