%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% DESCRIPTION: Compute chi2 and its gradient.
%
% VERSION:     $Id: chi2.m,v 1.1 2011/03/10 16:55:34 congedo Exp $
%
% HISTORY:     10-03-2011 G. Congedo
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [chi2,g] = chi2(p,data,models,Dmodels,lb,ub)

  weights = 1;
  Nfreeparams = numel(p);
  Nmdls = numel(models);
  Ndata = numel(data(:,1));
  Np = numel(p);
  
  mdldata = zeros(Ndata,Nmdls);
  for kk=1:Nmdls
    mdldata(:,kk) = models{kk}(p);
  end  
  res = (mdldata-data).*weights;
  if all(p>=lb & p<=ub)
    chi2 = res'*res;
    chi2 = sum(diag(chi2));
  else
    chi2 = 10e50;
  end
  
  if nargout > 1 % gradient required
    grad = cell(Nmdls,1);
    g = zeros(Nmdls,Nfreeparams);
    for kk=1:Nmdls
      grad{kk} = zeros(Ndata,Np);
      for ii=1:Np        
        grad{kk}(:,ii) = Dmodels{kk}{ii}(p);
        g(kk,ii) = 2.*res(:,kk)'*grad{kk}(:,ii);
      end
    end
    if Nmdls>1
      g = sum(g);
    end
  end
  
end