function idx=HW13Asol(X,Y)
% Input:  X is n x p and Y is n x m.  Output the indices found by the OLS
% algorithm.

[n,p]=size(X);  [n,m]=size(Y);

%Normalize X and Y

%SOLUTION:
d1=sqrt(sum(X.*X)); X=X./repmat(d1,n,1);
d2=sqrt(sum(Y.*Y)); Y=Y./repmat(d2,n,1);

%Alternative:  X=normc(X);  Y=normc(Y);  % normc is only included in the
%neural network toolbox however, so would not be available in Matlab only.

for i=1:p-1  
    
  %Compute cos^2(theta) between the p columns of X and m columns of Y
  %(Should be a p x m matrix for the first time through).
  C=(X'*Y).^2;

  %Get the best row norm of C.  You might find Matlab's max function useful
  %  (You should look up the help file for it)
  d3=sum(C'.*C');
  [vals,winidx]=max(d3);
  
  
  %Put the winning column as xwin and delete the corresponding column from
  % the matrix X.
  xwin=X(:,winidx);
  X(:,winidx)=[];
  idx(i)=winidx;
  
  %Deflate the remaining columns of X.  You may use a loop if you are more
  %comfortable with that:
  
  [nn,pp]=size(X);
  for j=1:pp
    X(:,j)=X(:,j)-(X(:,j)'*xwin)*xwin;
  end
  
  %Re-normalize the columns of X:
  d1=sqrt(sum(X.*X)); 
  X=X./repmat(d1,nn,1);
  
end