%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                            %
%  EP - Metodos Numericos de Algebra Linear                  % 
%  Professor Marcelo Queiroz                                 %
%  Willian Gigliotti, nUSP 4963454                           %
%                                                            %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Funcao Plot_elipse
% Ela, como diz o nome plota uma elipse, porem de centro na origem, 
% pois sao apenas essas que eu uso no EP.

% para plotar a elipse eu apenas a parametrizo a elipse e plot 1001 pontos dela!

% o parametro num deve ser 1 ou 0, e diz se eu devo ou nao imprimir os vetores da elipse
% o parametro s_plot eh a strig que vai junto com a funcao plot, podendo entao dar nomes para as elipses

function Plot_elipse(v1,v2,num,s_plot)
	rx = sqrt(v1(1)^2 + v1(2)^2);
	ry = sqrt(v2(1)^2 + v2(2)^2);
	if(rx==0)
		rx = ry;
		v1 = v2;
		ry=0;
		v2 = [0,0];
	endif
	if(rx>0) 

	% essa T esta ligada com a inclinacao na elipse;
	T(1,1) = v1(1)/rx;
	T(1,2) = -v1(2)/rx;
	T(2,1) = v1(2)/rx;
	T(2,2) = v1(1)/rx;

	t = 0:2*pi/(1001):2*pi;
	for a=1:1001
		F(1) = cos(t(a))*rx;
		F(2) = sin(t(a))*ry;
		% uso T para inclina-la, provavelmente existe uma maneira melhor, mas...
		F = T*F;
		
		x(a) = F(1);
		y(a) = F(2);
	endfor
	
	plot(x, y, s_plot)
	if(num!=0)
		hold on;

		arrowx(1) = 0;
		arrowx(2) = v1(1);
		arrowy(1) = 0;
		arrowy(2) = v1(2);
		plot(arrowx, arrowy,"-@;Vetor 1;");

		arrowx(2) = v2(1);
		arrowy(2) = v2(2);
		plot(arrowx, arrowy, "-*;Vetor 2;");
		hold off;

	endif
	else
		x=0;
		y=0;
		plot(x,y);
		pause
	endif;
	
endfunction
% funcao que plota concluida!!!

% recebe At*A e acha os seus autovalores
function auto_valores = autovalor(B)
	% B recebe A transposta vezes A;
	% a,b e c abaixo recebem os valores do coeficientes do polinomio
	% caracteristico de B, ax^2 + bx + c , det(B-LI)=0 , L autovalor
	a = 1;
	b = -B(1,1)-B(2,2);
	c = B(1,1)*B(2,2) -B(2,1)*B(1,2);
	% autovalores 1 e 2
	avt(1) = (-b + sqrt(b*b -4*a*c))/(2*a);
	avt(2) = (-b - sqrt(b*b -4*a*c))/(2*a);

	% eu retorno os valores em ordem crescente, facilitando na decomposicao espectral mais tarde
	if(avt(1) < avt(2))
		aux = avt(2);
		avt(2) = avt(1);
		avt(1) = aux;
	endif
	auto_valores = avt;
	
endfunction

% funcao que auxilia a proxima funcao (auto_vetor), nela entra a matriz M = A - LI e acho um vetor V, associado a L
function avet = aux_avet(M)
	% como estamos em R2, temos 2 autovalores e cada auto-valor tem pelo menos 
	% um auto vetor associado a ele, entao, a matriz M tem suas linhas multiplas 
	% uma da outra, pois elas nao podem ser LI, pra resolver o sistemas escolho uma delas!
	% escolho a nao nula
	if(M(1,1)==0 && M(1,2)==0) 
		indice = 2;
	else
		indice = 1;
	endif
	v(1) = M(indice,2);
	v(2) = -M(indice,1);
	v = v/(sqrt(v(1)^2 + v(2)^2));
	avet = v;
endfunction


function auto_vetor = autovetor(B,x)
	% observacao , como sempre AtA tem decomposicao espectral
	% entao, se x1 = x2, o meu espaco vai ser R2
	I = [1,0;0,1];
	if(x(1) == x(2) )
		v = [1,0;0,1];
	else
		M = B - x(1)*I;
		avet = aux_avet(M);
		v(:,1) = avet;
		M = B - x(2)*I;
		avet = aux_avet(M);
		v(:,2) = avet;
	endif
	auto_vetor = v;
endfunction

function normas(A, norma2)
	% achando norma 1 da matriz A
	max1 = abs(A(1,1)) + abs(A(2,1));
	max2 = abs(A(1,2)) + abs(A(2,2));

	norma1 = max(max1, max2);
	
	printf("Norma1 = %f\n", norma1 );
	
	printf("Norma2 = %f\n", norma2);

	% achando norma inf de A
	max1 = abs(A(1,1)) + abs(A(1,2));
	max2 = abs(A(2,1)) + abs(A(2,2));

	normainf = max(max1, max2);

	printf("Norma inf = %f\n" ,normainf);

	% achando Norma F
	AtA = A'*A;
	normaf = sqrt( AtA(1,1) + AtA(2,2));
	printf("Norma F = %f\n", normaf);

	printf("Aperte qualquer Enter pra continuar\n");
	pause;
	
endfunction
	
		

% Basicamente eu faco o seguinte, acho AtA, e seus auto valores e auto vetores pela funcoes auxiliares, que sao melhor explicadas nelas mesmas
% Mando Plotar num grafico o circulo Unitario e depois a transformacao dele

function decompoe2x2(A)

B = A'*A;
x = autovalor(B);
V = autovetor(B,x);
Vl = B*V;

printf("Decomposicao espectral de (At)(A): \n");
printf("AtA =\n");	
disp(V);
disp("*");
disp(diag(x));
disp("*");
disp(V');


printf("Aperte algum botao para ver o grafico.\n");
pause;
Plot_elipse( [ Vl(1,1),Vl(2,1) ] ,[Vl(1,2),Vl(2,2)], 0, "-;Elipse Imagem;");
hold on;
Plot_elipse( [ V(1,1),V(2,1) ] , [V(1,2), V(2,2)], 1, "-;Elipse Unitaria;");
hold off;
pause;

% na decomposicao de A em valores singulares eu uso que que 
% AtA = ((US(Vt))t)*(US(Vt))    * eh vezes e nao adjunto, no exemplo
%     = ((Vt)t)S(Ut)US(Vt))     Como U eh Unitaria
%       (VSIS(Vt))              onde I eh a Identidade
%       V(S^2)(Vt)             
	
%       conclusao:
%	V da Decomposicao em Valores Singulares de A eh
%       o V da Decomposicao Espectral de AtA, ja feito acima
%       E a matriz S de A=US(Vt) eh a raiz dos Autovalores de AtA
%       Como AV = US
%       Facilmente encontro US, depois encontro U, explicado Abaixo...
%       Depois eh so juntar os Resultados e Plotar no Grafico!

       

printf("Decomposicao em Valores Singulares de A: \n");
printf("A =\n");
S = diag(sqrt(x));
Ul = A*V;
%vou achar o U!
if(S(1,1)==0)
	% se S(1,1) for 0 entao o outro valor singular tb eh zero, portanto U = Identidade;
	U = [1 0; 0 1];
else
	U(:,1) = Ul(:,1)/(S(1,1));
	if(S(2,2)==0)
		% se S(2,2) for 0 entao, vou pegar um vetor ortogonal quaquer de U1
		U(1,2) = U(2,1);
		U(2,2) = -U(1,1);
	else
	U(:,2) = Ul(:,2)/(S(2,2));
	endif
endif

disp(U);
disp("*");
S = diag(sqrt(x));
disp(S);
disp("*");
disp(V');

printf("Aperta algum botao para ver o grafico. Os Vetores Singulares\n");
pause;
Plot_elipse( V(:,1), V(:,2),1, "-5;Elipse Unitaria de R2;" );
printf("Aperta algum botao para ver o Grafico. As Transformacoes de A nos valores Singulares\n");
pause;
Plot_elipse( Ul(:,1), Ul(:,2),1,"-5;Elipse Imagem da Unitaria por A;" );
pause;


normas(A, S(1,1));

% como S(2,2) eh o menor valor de S, e todos sao maiores ou iguais a zero, entao A vai ser inversivel se S(2,2) diferente de 0
% A^(-1) = (USV*)^(-1)
%  = V(S^(-1))U*
% Todo o resto a feito de Maneira Analoga A Matriz A

if(S(2,2)!=0)
	Sl = diag(diag(S).^(-1));

	% Apesar de Al = V*Sl*U estar correto, a matriz SL
	% não eh Sigma de Al, pois a ordem dos valores singulares estao invertidas
	% se x > y e x>0 e y>0 entao 1/x < 1/y
	% portanto posso escrever que a matriz Sln =
	% |0 1|        |0 1|       
	% |1 0| x Sl x |1 0|
	% onde Sln fica sendo a Sigma verdadeira da Inversa de A
	% agora eh so fazer V * [ (0,1), (1,0) ] 
	% pegar a Sl correga
	% e [ (0,1) , (1,0) ] * (Ut)
	% pronto, agora eu tenho a decomposicao de valores singulares da inversa de A
	M = [0 1; 1 0]
	V = V*M;
	Ut = M*U';
	auxi = Sl(1,1);
	Sl(1,1) = Sl(2,2);
	Sl(2,2) = auxi;
	
	Al = V*Sl*Ut;
	printf("A eh Inversivel e sua inversa: \n");
	disp(Al);
	
	printf("Decomposicao em Valores Singulares de A^(-1): \n");
	printf("A^(-1) =\n");

	disp(V);
	disp("*");
	disp(Sl);
	disp("*");
	disp(Ut);

	printf("Aperta algum botao para ver o grafico. Os Vetores Singulares\n");
	pause;
	Plot_elipse( Ut(:,1), Ut(:,2) , 1, "-5;Elipse Unitaria de A^(-1); ");
	printf("Aperta algum botao para ver o Grafico. As Transformacoes de A^(-1) nos valores Singulares\n");
	pause;
	Vl = Al*U;
	Plot_elipse( Vl(:,1), V(:,2), 1, ";Elipse Imagem da Inversa de A;" );
	pause;

	normas(Al, Sl(1,1));

else
	printf("A Matriz A nao eh inversivel\n");
endif
endfunction
% Matriz A para textes
% A = [0.01 -0.07; -0.4 0.2]
% decompoe2x2(A);


