function [bestmse, method] = compare(images)
% compares the expand methods on a training set
% returns a table mapping cell numbers to optimal method in that cell

% used for hybrid method training


[row col nimages] = size(images);

a = 0.4;
w1 = [1/4-a/2; 1/4; a; 1/4; 1/4-a/2];
W1 = w1*w1';
a = 0.6;
w2 = [1/4-a/2; 1/4; a; 1/4; 1/4-a/2];
W2 = w2*w2';
load NLcoeffs;

nmethods = 10;

maxpart = 400;
method = ['burt4';'burt6'; 'ideal'; 'haar '; 'MMF  '; 'spli1'; 'spli3'; 'opt_l'; 'opt_c';'optNL';];
   
   
mse = zeros([4, maxpart, nmethods]);
for index = 1:nimages
   image = images(:,:, index);
   simage = reduceFunc(image, 'opt_NL');
   
   rimage = zeros([row col nmethods]);
   if (1)
      rimage(:,:,1) = expandFunc(simage, 'burt', W1); disp('burt  a=0.4');
      rimage(:,:,2) = expandFunc(simage, 'burt', W2); disp('burt  a=0.6'); 
      rimage(:,:,3) = expandFunc(simage, 'ideal'); disp('ideal');
      rimage(:,:,4) = expandFunc(simage, 'haar'); disp('haar');
      rimage(:,:,5) = expandFunc(simage, 'MMF', 3); disp('MMF');
      rimage(:,:,6) = expandFunc(simage, 'spline', 1); disp('linear spline method');
      rimage(:,:,7) = expandFunc(simage, 'spline', 3); disp('cubic spline method');
      rimage(:,:,8) = expandFunc(simage, 'opt_lin', 15); disp('optimal linear');
      rimage(:,:,9) = expandFunc(simage, 'opt_cub', 15); disp('optimal cubic');
      rimage(:,:,10) = expandFunc(simage, 'opt_NL', NLcoeffs); disp('optimal non-linear');
   end
   
   error = zeros([row*col nmethods]);
   for i=1:nmethods
      error(:,i) = reshape((rimage(:,:,i)-image).^2, [row*col 1]);
      
   end
      
   part = partition(simage);

   for j = 0:max(max(part))
      [indx indy] = find([part == j]);
      for k = 0:1
         for l = 0:1
            mask = (2*(indx-1)+1+k)+col*(2*(indy-1)+l);
            cellsize = length(mask);
            disp(sprintf('subimage (%d,%d) cell %d : %d pixels', k, l, j, cellsize));
            if (cellsize > 0)
               for i = 1:nmethods
                  mse(2*k+l+1, j+1, i) = mse(2*k+l+1, j+1, i) + mean(error(mask, i)) / nimages;
               end
           end
        end
     end
  end
end

[bestmse, method] = min(mse, [], 3);


     
      