function [Px P]=reflect(x,u)

% Implmement householder reflections:
% Relect x across orth(u).  In 2d, graph it


[n m]=size(x);

if n~=length(u)
    error('x andu must be the same length')
end

if m~=1 | n==1
    error('x and u must be column vectors')
end

% P is the Householder reflection operator
if nargout==2
    P=eye(n) - 2*u*u'/norm(u)^2;   % Inefficient to calculate P, then do P*x
end

% Reflect x across orth(u): calculate P*x
Px=x-2*u*(u'*x/norm(u)^2);

if n==2
    figure(1)
    clf
    plot([0 x(1)],[0,x(2)])
    hold on
    plot([-u(1) u(1)],[-u(2),u(2)],'g')
    e1=randn(2,1);
    orthu=(e1 - u*(u'*e1/norm(u)^2));
    orthu=orthu*norm(u)/norm(orthu);   % Make the same length as u for graphing purposes
    plot([-orthu(1) orthu(1)],[-orthu(2) orthu(2)],'r')
    plot([0 Px(1)],[0 Px(2)],'k')
    axis equal
    legend('x','u','orth(u)','x reflected across orth(u)')
end