

#testujemy czy \ rozwiazwiaze LZNK
#czy metoda z rozkladem qr rozwiaze LZNK 
#czy metoda z ukladem r nomalnych A'Ax=A'b rozwize LZNK
N=20;
for k=1:2,
 A=randn(N,10);
 x=randn(10,1);
 b=A*x;
 [q,r]=qr(A);
 norm(q*q'-eye(N),1)
 norm(q'*q-eye(N),1)
 norm(A-q*r,1)
 y=A\b;
 g=q'*b;
 z=r(1:10,:)\g(1:10);
 norm(g(11:N)) #czy zero?
 w=(A'*A)\(A'*b);
 [norm(x-y), norm(x-z),norm(x-w)]
 N=N*10;
endfor

#porownajmy czas \ qr i uklad r normlanych
N=20;
for k=1:2,
 A=randn(N,10);
 b=randn(N,1);
 [q,r]=qr(A);
 tic;y=A\b;t(k,1)=toc;
 tic; g=q'*b;z=r(1:10,:)\g(1:10);t(k,2)=toc;
 tic;w=(A'*A)\(A'*b);t(k,3)=toc;
 N=N*10;
endfor
t

#krzywa ax^2+by^2-1=0  do pktow x_k y_k
N=10;
for k=1:1,
 x=linspace(0,1,N+1);x=x';
 #pkty z elipsy
 y1=sqrt(1-x.*x)*0.5; 
 #pkty blisko elipsy
 y2=y1+1e-1*(2*rand(N+1,1)-1);
 A1=[x.*x,y1.*y1]; 
 A2=A1;A2(:,2)=y2.*y2;
 b=ones(N+1,1);
 w1=A1\b
 w2=A2\b
 s=linspace(1e-3,1-1e-3,100);
 if(N<20)
  plot(x,y2,"+",s,sqrt((1-w2(1)*s.*s)/w2(2)));
  pause(1);
 endif
 N*=2;
endfor

#to samo ale na hiperboli
disp("szukamy hiperboli");
 N=20;x=linspace(1,3,N);x=x';
 y1=sqrt((1+x.*x)/2);
#pkty blisko 
 y2=y1+1e-1*(2*rand(N,1)-1);
 A1=[x.*x,y1.*y1]; 
 A2=A1;A2(:,2)=y2.*y2;
 b=ones(N,1);
 w1=A1\b
 w2=A2\b
 s=linspace(-1,3,100);
 pause(1);
 plot(x,y2,"+",s,sqrt((1-w2(1)*s.*s)/w2(2)),s,sqrt((1-w1(1)*s.*s)/w1(2)),";y^2-x^2=1;")
 pause(1);
norm(A2*w2-b)
norm(A2*w1-b) #czy wieksze?


#testy Householdera

function Y=HV(X,w,snw=w'*w)
#mnozymy X m x n przez 
#macierz Householdera H mxm z wektorem House. w  m x 1
#Input
#X macierz m x n
#w wektor Householdera m x 1
#snw=norm(w,2) (opcja)
#OUTPUT
#Y=H*X - m x n -

Y=X- (2/snw)*w*(w'*X);

endfunction


function Y=BadHV(X,w,snw=w'*w)
#bad implementation of H*X (H Householder matrix)
#wrong brackets placement - cost O(m^2) +O(m*n)
Y=X- (2*(w*w')*X)/snw;
endfunction

function Y=BadHV1(X,w,snw=w'*w)
#bad implementation of H*X (H Householder matrix)
#wrong brackets placement - cost O(m^2) +O(m*n)
H=eye(length(w))- 2*(w*w')/snw;
Y=H*X;
endfunction


#testujemy czy dziala

m=5; n=3;
X=rand(m,n);
w=rand(m,1);
Y=HV(X,w);
Y1=BadHV(X,w);
Y2=BadHV1(X,w);
norm(Y-Y1,1)
norm(Y-Y2,1)
#same results...
#if H=H^{-1}
X1=HV(Y,w);
norm(X1-X,1)
#if Frobenius norms are equal 
[norm(X,'fro') ,norm(Y,'fro')]

H=eye(length(w))- 2*(w*w')/(w'*w);
#if H==H'
norm(H-H',1)
#If H*H=I
norm(H*H-eye(size(H)),1)

#we zero elements in the 1st col of random X 
X=rand(m,n);
w=X(:,1);w(1)-=norm(X(:,1),2);
Y=HV(X,w)
norm(Y(2:m,1),2) #==0?

#we test times
N=10;M=2*N;
for k=1:4,
 X=rand(M,N);
 w=rand(M,1);
 hvt(k,1)=M;hvt(k,2)=N;
 tic;Y=HV(X,w);hvt(k,3)=toc;
 tic;Y1=BadHV(X,w);hvt(k,4)=toc;
 tic;Y2=BadHV1(X,w); hvt(k,5)=toc;
 N*=5;M*=5;
endfor
hvt
hvt(:,4)./hvt(:,3)
#ponad 100 razy wolniejszy











