%%%%%%%%%%%%%  Polynomial Least Squares  with Mouse Input   %%%%%%%%%%%
%
% MOUSELSQ
% Usage: MOUSELSQ
%
%
% Demonstrates least squares interpolation for polynomials of several degrees  
%              - the data and interpolating polynomials are drawn
%              - points are entered via the mouse
%              
%              the window and the interpolating degrees can be changed 
%              by editing

clf; echo off;

%% set up regression parameters
degs = [1 2 3 4 5 6 7 8 9 10];   % degrees of regression polynomials
a  = 0;                          % left interpolation endpoint
b  = 10;                         % right interpolation endpoint
M  = max(degs);
                  
%% set up graphics screen
minx = a;  maxx = b; 
miny = 0;  maxy = 10;
win = [minx maxx miny maxy];
clf
axis(win)                       % set up sceen window
hold on
plot(minx-1,miny-1);             % this draws the axes 
title('Regression for several degrees')
xlabel('Select points, Press Enter or right mouse button to quit'); 

%% get data points
X = [];
Y = [];
btn = 1;                              % set up while loop
hold on
[x,y,btn] =ginputsab(1);              % get mouse coordinates
while (btn==1)                        % 1 = left 2 = right or Enter
   X = [X x];                         
   Y = [Y y];
   plot(x,y,'k*');
   [x,y,btn] =ginputsab(1);
   end;
hold off

%% set up regression matrix for all polys 

X = X(:); Y = Y(:);               % flatten x and y if neccessary 
A = [];                           
for j = 0:M A = [ X.^j A]; end    % calculate big interpolating matrix one
                                  % column at a time in ascending powers
                                  % x.^j exponentiates elementwise
                                   

% display results
clf
cols = 'ymcrgb';
axis(win);
plot(X,Y,'k*')                   % plot data points with white stars
title('Regression for several degrees')
xlabel('Press a key for the next curve') 
pause
hold on

for n = degs                      % for all degrees in degs
    B = A(:,M-n+1:M+1);           % pick off submatrix 
    c = B\Y;                      % calculate regression coefficients
                                  % by solving solves A'*Y=(A'*A)*c
    clc                           % clear screen and show values

    gcode = [cols(1+rem(n-1,6)) '-'];
                                  % different color for each curve
    px = a:(b-a)/100:b;           % x-values for plot (100 intervals)
    pp = polyval(c,px);           % y-values for polynomial approximation plot
    plot(px,pp,gcode)             % plot regression curve in color gcol
    pause
    end                          
ylabel('done, Press a key'); 
hold off