%  *Homework, due tuesday, October 2*: Write a matlab function, following
%  this lecture, that given any size array of points x and y, finds the best
%  fit line for these points, *and* reports the squared error for each point
%  and the sum of squared errors (of all points).  Your function should
%  start with something similar to:
%
%  function [errors, sse] = fitLine(x,y);
%
%  Finally, when writing up your assignment, you should not look at
%  anyone elses source code.  Help between students is limited to general
%  questions and questions about specific error messages.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Lecture 4
%
% Grading policies for this class:
%
%  The grade for this class will be determined based on a numerical
% threshold for your score, with the common gradingscale.
%
%  [90-100] = A
%  [80-90): = B
%  [70-80): = C
%  [60-70): = D
%  [ 0-60): = F
%
%  The score will be made up of the weighted combination of in-class
%  quizes, labs, and homeworks.
%
%  This weight will be approximately 40% quizes, 30% homeworks, 30%
%    labs.
%
%  The quizes will start to be every week.  They will not have neither
%  partial credit nor extra credit.  However, for every question that
%  appears on quiz i, there will be a "replacement" question on quizes
%  i+1 and i+2.  Answering this "replacement" question correctly will
%  give you full credit on the original quiz question.
% 
% Last week:  Functions, if-then construct and for loops.
%
% This week.  
% Topic 1, More throughts on Gauss, and practice with functions and
%          variables in functions.
%
% Topic 2, Polynomial fitting and plotting.
%
%  *Homework*:  Write a matlab function, following this lecture, that given
%  any size array of points x and y, finds the best fit line for these
%  points, *and* reports the squared error for each point and the sum of
%  squared errors (of all points).  Your function should start with
%  something similar to:
%
%  function [errors, sse] = fitLine(x,y);
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Last time we talked about Gauss, and how he works with functions,
%   for-loops and if-then statements.
%
%
% Let's get started today with a warm up.  What does the following code
% do? 

array = [];
for ix = 1:40000
  array = [array,ix];
end

% And now, if we think about Gauss, what steps does he go through in
% executing this line?
%  Suppose that size(a) = [1 n]
%  a = [a, i]
%  The right hand side evaluates to a 1 x (n+1) array
%  Gauss asks the operating system for room to write it
%  He writes it down, and then labels it "a"
%  He crosses out the previous value of "a"

%  So how many steps does Gauss have to do?
%
%  array = [];
%  for i = 1:40000
%    array = [array,i];
%  end
%  First time: 1 number must be copied
%  Second time: 2 numbers must be copied
%  40000th time: 40,000 numbers must be copied
%  Total: 1 + 2 + 3 + ... + 40,000 ~ .8 billion steps!
%
%
% much faster would be:
array = zeros(1,40000);
for ix = 1:40000
  array(ix) = ix;
end
% (what does Gauss do here?)
% and much, much faster is the line:
array = 1:40000
%
% (which is really the same, but is a built-in matlab function.
%
%  Another program control structure that Gauss knows is a "while"
%  loop.  This is like a for loop, but for the case when you may not know how
%  many times the loop needs to be called.
%
% the general format is 
%  WHILE (CONDITION)
%   ... do stuff
%  end
%
% Pop Quiz: What is the log(n)?
%
function powerTwo(n)
% find a number k with 2^k > n
k = 0
while (n > 1)
  n = n / 2;
  k = k+1;
end

%
% Now, a reprise and extenstion of the last lab.
% define:
x = (1:10)';          % make a column vector
% and 
y = 2*x + rand(10,1);
%
%
% given collection of points x,y, let's fit a line to them...
%
% y = ax + b;
%
% we can write this in matlab notation to remember which variables are
% really vectors:
%
%   y(:) = a * x(:) + b;
%
% which is actually many equations
%    y(i) = a * x(i) + b;
%
% How can we write this in matrix form?
% 
%    y(i) = a * x(i) + b * matrixOfOnes(i);
 
%or: 

%[x(1)  1][a]     [y(1)]
%[x(2)  1][b]  =  [y(2)]
%
%[x(n)  1]        [y(n)]

%this has the form:
%M v = y
%
% How do I make M from x(:)?
%
% What are the unknowns?
%
% How to solve for this?
%
% (can I invert M?)
%
% (can I make it so that I can invert M)

% Solution 1:

M = [x, ones(10,1)];

%%% Classic trick, how to solve Mv = y?

%Mv = y
%M'Mv = M'y
%inv(M'M)M'Mv = inv(M'M)M'y
%v = inv(M'M)M'y
v = inv(M'*M)*M'*y   % should all be *, not .*
%
v = pinv(M)*y
%
% or... matlab can do all this for you...
%Mv = y
%M\Mv = M\y
%v = M\y

v = M\y   % which is a little better...

% so how can we plot these values?

% and what does this minimize?