%% Polynomial Data Fitting 
% Given a set of points in R^2, find the polynomial function that contains them


%% The Given Data
% The top row is the independant variable, "x", the bottom row "y"

%% The General Associated Polynomial
%
% For the given number of points ...
	X = -5+10*rand(2,5); 
	n = max(size(X));

%% The "Coefficient" Matrix
% The polynomial is a_0 + a_1*x + ... + a_m*x^m where m = n-1;
% This general form determines one equation for every given pair
% a(1) + a(2)*X(1,1) + a(2)*X(1,1)^2 = X(1,2)
	M = ones(n,n+1);
	for point = 1:n
		for i = 1:n;
			M(point,i) = X(1,point)^(i-1);
		end
		M(point,n+1) = X(2,point);	
	end


%% Solve the system
	A = rref(M);
	for i = 1:size(M,1); 
		a(i) = A(i,end);
	end

%% The polynomial
% For example, y = a(1) + a(2)*x + a(3)*x.^2;

	% domain of appropriate x-values
	r = max(abs(X(1,:)))+20; 
	x = -r:.01:r;

	% create the poly using the coef found in rref(A)
	s = [num2str(a(1))];
	for i = 2:length(a);
		s = [s,' + ',num2str(a(i)),'*x.^',num2str(i-1)];
	end

	% evaluate the poly over the domain to get y-values
	y = eval(s);
	
%% Plot the points and the curve
	%figure;
	plot(0,0); % Blank plot to clear everything
	l = 2; % enlarge view window
	ylim([min(X(2,:))-l max(X(2,:))+l]);
	xlim([min(X(1,:))-l max(X(1,:))+l]);
	hold on
	% display the poly
	title(s);
	% display points
	plot(X(1,:),X(2,:),'*b','MarkerSize',12);
	pause
	% display graph of poly
	plot(x,y,'r',X(1,:),X(2,:),'*b','MarkerSize',12);
	hold off;

