function curve = parabolicregress(X,Y,Z)
% PARABOLICREGRESS Fit a parabolic regression to data.
%   CURVE = PARABOLICREGRESS(X,Y) returns a vector of coefficients [c b a]
%   which correspond to the parabolic function ax^2 + bx + c.  For n data
%   points, X is a 1-by-n vector of predictors and Y is a 1-by-n vector 
%   of observations.
%
%   CURVE = PARABOLICREGRESS(X,Y,Z) where Z is a vector of weights,
%   performs a weighted least squares regression.
%
%   See also LINEARREGRESS, REGRESS, LSCOV

% Copyright (c) 2011
% Jerod Weinman
% jerod@acm.org
  
% $Id$
 
% Jerod Weinman, Zach Butler, Dugan Knoll, Alexander Turner, Yirong Jing
% jerod@acm.org
% Copyright 2011, 2012, 2013, 2014
 
%    This file is part of PIONEER WordNorm.
%
%    PIONEER WordNorm 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.
%
%    PIONEER WordNorm 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 PIONEER WordNorm.  If not, see <http://www.gnu.org/licenses/>.
  
% Subtract off the averages of the values
    x_mu = mean(X);
    X = X - x_mu;
    y_mu = mean(Y);
    Y = Y - y_mu;
    
    % Set up the X values for the regression  Y = aX^2 + bX + c
    Xfunc = [X.^2 X ones([length(X) 1])];
    
    if nargin<3
      abc = regress(Y, Xfunc);
    else
      abc = lscov(Xfunc,Y,Z);
    end;

    % We found the regression with the origin being centered around the
    % middle of the data points.  Now we want to readjust the origin to be
    % where it was initially.  To do this, we will complete the square and
    % subtract off the mu values.
    
    % Since we cannot divide by zero, determine if it is quadratic
    if ~(abc(1)==0)
    % ax^2 + bx + c = a(x - h)^2 + k  where h = -b/2a  and k = c - b^2 / 4a
        h = - abc(2) /( 2*abc(1) );
        k = abc(3) - abc(2)^2 / (4*abc(1));
    % Adjust the origin
        h = h + x_mu;
        k = k + y_mu;
        h = -h;
    
    % Put it into our quadratic curve
        curve(3) = abc(1);
        curve(2) = abc(1) * 2 * h;
        curve(1) = abc(1) * h^2 + k;
    else
        curve(3) = 0;
        curve(2) = abc(2);
        curve(1) = abc(3) + abc(2)*x_mu - y_mu;
    end
    
end
