function [u, v, lastexpy] = newton_gaussian(x, nb_axis, nb_iterations, ...
                                    startu, startv,eps_stop,epsilon_cond)

% newton_gaussian computes the Newton algorithm to factorize the matrix of 
% of the noisy patches in the Gaussian domain (mainly it perform a PCA 
% without enforcing the constraint of orthogonal axis).
%                 
%   INPUT:
%     x                      : noisy patches 
%  	  nb_axis                : number of axis/components choosen 
%  	  nb_iterations		     : number of iterations to perform the Newton
%                              step
%  	  startu                 : initialization of the coefficient in the
%                              matrix factorization / pca representation
%     startv                 : initialization of the axis/components in the
%                              matrix factorization / pca representation
%     eps_stop               : relative stopping  criterion for the Newton
%                              step
%     epsilon_cond           : small conditioning number to be sure the
%                              Hessian is well conditionned when doing 
%                              the Newton step
%
%   OUTPUT:
%     u                      : coefficient of the patches in the PCA-domain
%    
%     v                      : axis/component/dictionary element to 
%                              represent the patches in the PCA-domain
%
%     lastexpy               : representation of the patches using only the
%                              nb_axis elements
%
% [u, v, lastexpy] = newton_gaussian(x, nb_axis, nb_iterations, startu,...
%                                            startv,eps_stop,epsilon_cond)
%
%   Compute the Newton algorithm to factorize the matrix of noisy patches,
%   using nb_axis, a maxim amount of nb_iterations for the main loop and
%   a stoping criterion eps_stop, for controlling the relative change. 
%   startu and starv are initial guess for the factorization. 
%   
%   Copyright (C) 2012 NL-PCA project
%   Joseph Salmon, Charles-Alban Deledalle, Rebecca Willet, Zachary Harmany
%
%   See The GNU Public License (GPL)

%---------------------------------------------------------------------
%
%   This file is part of NL-PCA.
%
%   NL-PCA is free software: you can redistribute it and/or modify
%   it under the terms of the GNU General Public License as
%   published by the Free Software Foundation, either version 3 of
%   the License, or (at your option) any later version.
%
%   NL-PCA is distributed in the hope that it will be useful,
%   but WITHOUT ANY WARRANTY; without even the implied warranty of
%   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%   GNU General Public License for more details.
%
%   You should have received a copy of the GNU General Public
%   License along with NL-PCA.  If not, see
%   <http://www.gnu.org/licenses/>.
%   Joseph Salmon, Charles-Alban Deledalle, Rebecca Willet, Zachary Harmany
%
%   See The GNU Public License (GPL)

numelx=numel(x);
[n, d] = size(x);
u = startu;
v = startv;

v(1,:)=ones(1,d);
v = v./repmat(sqrt(sum( v.^2, 2 )),[1,d]);  
condi=epsilon_cond*eye(nb_axis);

lastexpy=u*v;
lasterr=x-lastexpy;

for iter = 1:nb_iterations
  for i = 1:n     
      err=lasterr(i,:);   
      covar = (v * v' +condi);
      u(i,:)=u(i,:)+ (err*v')/covar;
  end
  
  lastexpy2 = (u*v);
  lasterr= x-lastexpy2;
  
  for i = 1:d
      err=lasterr(:,i);       
      covar = (u' * u + condi);  
      v(:,i) =v(:,i)+ covar \ (u' * err); 
  end 
  
  normalization=sqrt(sum( v.^2, 2 ));
  v = v./repmat(normalization,[1,d]); 
  lastexpy2=u*v;
  curnorm = sum(sum(lastexpy2)) / numelx;
  change = sum(sum(abs(lastexpy2-lastexpy))) / numelx;
  ratio= change / curnorm;
  % verbose
  % fprintf('%3d: %g %g %g \n', iter, curnorm, change,ratio);
  if (i > 1 && (ratio < eps_stop))
  break;
  end
  lastexpy=lastexpy2;
  lasterr=x-lastexpy2; 
end
