function [yhat,coefs] = lwreg1(data,query,h,order)

% LWREG1  --  Locally weighted regression in one dimension.
%
%  YHAT = LWREG1(DATA,QUERY,H)
%  YHAT = LWREG1(DATA,QUERY,H,ORDER)
%  [YHAT,COEFS] = LWREG1(...)
%
%  DATA is a matrix with two columns: X and Y.
%  X is the predictor variable and Y is the target variable.
%  QUERY is a column vector of X values for which prediction is requested.
%  H is a number specifying the bandwidth of the kernel.
%  ORDER is either 1 for linear (=default) or 2 for quadratic regression.
%
%  Algorithm: For each query point Q:
%   1. Find the set of "neighbors" -- data points less than H units away.
%   2. Calculate the weight of each neighbor using a quadratic kernel:
%        dist = |q - neighbors| ; wght = max(0,1-(dist/h).^2) ;
%        (The quadratic kernel is reportedly optimal in a MSE sense.)
%   3. Find regression coefficients COEFS that minimize the weighted
%      sum of squares.
%   4. Use COEFS to predict the target value YHAT at the query point.
%
%  The function returns a column vector YHAT of the same length as QUERY.
%  A matrix of COEFS may also be requested. It has one line for each
%  query point. The first coefficient is the constant term, the second 
%  is the linear term, and the third is the quadratic term (if order=2).
%
%  Example:
%    x = sort(rand(100,1)*20 - 10) ;
%    y = sin(x) + randn(100,1)*0.2 ;
%    plot(x,y,'b.',x,sin(x),'b--') ;
%    query = [-10:.1:10]' ;
%    yhat = lwreg1([x,y],query,2,2) ;
%    hold on ; plot(query,yhat,'r-') ; hold off ;
%
%  Reference:  Atkeson, C.G., Moore, A.W., & Schaal, S. (1997).
%              Locally weighted learning. AI Review, 11, 11-73.
%
%  See also CVERR_LWR1, LOCALREG, KDE.

% Original coding by Alex Petrov, Ohio State University
% $Revision: 1.0 $  $Date: 2000/11/28 20:20 $
%
% Part of the utils toolbox version 1.1 for MATLAB version 5 and up.
% http://alexpetrov.com/softw/utils/
% Copyright (c) Alexander Petrov 1999-2006, http://alexpetrov.com
% Please read the LICENSE and NO WARRANTY statement in ../utils_license.m

%--  Decipher input arguments
if (nargin < 4)
  order = 1 ;       % default is quadratic
end

x = data(:,1) ;
y = data(:,2) ;
N_data = length(x) ;
N_query = length(query) ;

%--  Prepare left-hand-side values for the regressions
if (order==1)    % linear regression
  A = [ones(N_data,1) , x ] ;
  Q = [ones(N_query,1) , query] ;
else             % quadratic regression
  A = [ones(N_data,1) , x , x.^2] ;
  Q = [ones(N_query,1) , query , query.^2] ;
end

yhat = zeros(N_query,1) ;
coefs = zeros(1+order,N_query) ;   % will be transposed later

%--  Main loop
for k = 1:N_query
  q = query(k) ;
  %- Find the indices of all neighbors
  dist = abs(x-q) ;
  neighbors = find(dist<h) ;
  switch (length(neighbors))
  case 0
    warning('No neighbors found for this query. Try wider bandwidth.') ;
  case 1
    coefs(1,k) = y(neighbors) ;     % constant ; coefs(2:end,k) remain zero
  otherwise    % the typical case -- enough data to do regression
    %- Calculate the weights using the quadratic kernel
    wght = 1 - (dist(neighbors)/h).^2 ;
    wght = sqrt(wght) ;       % will be squared back inside the error term
    %- Filter and weight the data points
    X = A(neighbors,:) .* repmat(wght,1,1+order) ;
    Y = y(neighbors) .* wght ;
    %- Do the regression and store the coefficients for this query point
    coefs(:,k) = X\Y ;
  end
end 

%--  Calculate and export values
yhat = sum(Q' .* coefs)' ;
coefs = coefs' ;

%-- Return YHAT and possibly COEFS
%%%%% End of LWREG1.M
