% SNIP OFF anything above this line.
% 
% In directory: ~/www/public_html/numanal/matlab/
%
% Hi Dad,
% 	I'm sending you the following matlab file, in the hopes that you can
% test it out and see how you like it.  It is one part of the demo that I'm
% going to do.
%
% In this part, I'm looking at the decomposition of Abe, and I'm also making
% an Abe function (which interpolates him), using spline-fit singular vectors.
%
%							Andy
%
colormap('default')
abe =[ 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8
 8 8 8 8 8 8 8 8 8 8 14 15 25 65 40 40 45 45 42 35 34 33 31 27 10 8 8 8 8 8 8 8 8 8 8 8
 8 8 8 8 8 8 8 8 8 8 28 30 50 65 80 80 90 90 83 71 68 67 63 54 20 8 8 8 8 8 8 8 8 8 8 8
 8 8 8 8 8 8 8 8 8 40 55 70 80 88 94 96 97 95 94 92 89 82 72 65 60 46 8 8 8 8 8 8 8 8 8 8
 8 8 8 8 8 8 8 8 30 50 75 87 94 98 100 100 100 98 97 94 91 80 75 70 65 60 40 8 8 8 8 8 8 8 8 8
 8 8 8 8 8 8 30 54 58 85 100 104 105 106 107 107 106 103 102 99 95 88 80 75 70 60 45 40 8 8 8 8 8 8 8 8
 8 8 8 8 8 8 35 45 56 70 82 93 98 103 106 108 107 108 108 106 103 100 94 78 74 63 50 40 8 8 8 8 8 8 8 8
 8 8 8 8 37 48 53 68 88 96 100 103 105 109 108 107 106 105 104 103 99 94 88 77 71 63 54 40 33 26 8 8 8 8 8 8
 8 8 8 8 26 45 60 77 89 95 96 99 101 103 106 107 108 107 106 105 104 100 93 86 78 68 57 50 47 46 44 40 8 8 8 8
 8 8 8 8 40 40 60 70 80 88 86 86 104 106 106 107 106 106 104 102 99 94 89 80 79 69 65 60 55 49 44 40 35 26 8 8
 8 8 8 8 30 45 47 58 66 73 82 82 82 86 87 92 95 99 100 100 97 94 89 83 76 70 65 63 62 56 48 43 37 31 26 8
 8 8 8 8 30 42 56 58 60 69 77 80 83 88 92 94 95 96 98 100 99 98 96 95 93 87 83 69 67 66 56 48 43 36 29 20
 8 8 8 10 32 38 50 55 60 65 75 85 90 93 95 97 100 102 103 104 104 104 103 102 100 97 92 85 77 69 68 67 52 35 26 8
 8 8 8 32 36 48 54 61 66 76 84 92 96 99 103 104 105 106 107 107 106 106 104 103 95 89 84 77 69 68 69 67 41 26 20 8
 8 8 28 29 42 50 56 63 74 82 89 97 104 105 107 108 108 108 108 107 106 104 102 99 96 90 85 77 70 68 69 68 45 27 18 8
 8 8 29 29 46 54 60 65 77 87 103 105 107 108 109 110 110 110 109 109 108 107 105 103 99 94 89 80 74 70 74 70 44 28 19 8
 8 8 14 30 46 53 58 68 77 90 94 98 102 105 107 108 109 109 109 108 107 106 105 103 101 99 95 87 74 68 69 64 42 28 10 8
 8 18 26 32 47 54 60 80 91 102 107 108 109 112 113 113 113 112 112 111 110 110 109 106 104 97 93 86 75 71 60 50 44 30 8 8
 8 14 28 36 40 49 57 59 85 95 101 106 108 111 113 113 113 113 113 113 112 111 110 107 106 100 95 84 75 67 58 50 44 38 25 8
 8 12 29 36 40 47 54 61 86 99 104 106 108 110 112 113 114 114 114 114 114 112 111 108 107 102 95 86 78 68 59 49 42 36 26 8
 23 21 23 29 41 48 56 68 85 105 116 116 115 114 114 114 114 114 114 114 114 115 116 116 115 114 111 108 84 66 58 50 40 35 30 20
 20 20 20 20 29 40 49 62 68 100 110 115 116 115 116 116 116 114 113 113 112 112 113 113 112 109 106 101 97 68 62 57 40 30 20 20
 20 20 20 20 30 40 54 66 103 109 111 113 114 115 116 115 115 117 116 115 114 113 112 112 111 109 107 103 98 68 64 59 50 24 20 20
 20 20 20 30 34 50 60 69 95 103 105 107 108 107 105 108 110 113 115 112 105 103 104 105 106 104 102 98 95 63 60 50 22 19 15 10
 20 20 20 30 30 27 56 60 63 92 101 104 105 104 104 99 100 113 117 116 104 102 105 105 106 105 94 90 80 63 57 30 24 20 12 10
 20 20 22 23 30 28 40 60 83 90 97 101 103 103 100 101 107 122 123 118 103 98 100 101 100 95 89 83 78 61 30 18 21 19 10 10
 8 8 20 25 26 30 50 64 80 95 97 99 100 100 100 103 113 123 125 118 104 99 99 99 98 95 90 85 80 55 40 20 18 10 8 8
 8 8 22 25 28 28 60 64 88 94 98 102 102 103 103 104 114 130 131 123 110 100 102 103 100 98 94 88 82 56 40 18 17 12 8 8
 8 8 23 24 26 27 58 65 90 97 99 100 101 102 103 107 120 131 131 125 110 103 101 101 100 97 93 88 83 56 45 20 18 14 8 8
 8 8 8 28 31 45 58 72 90 96 99 100 101 104 107 113 130 136 135 132 125 106 104 102 99 97 95 92 88 55 40 25 20 20 8 8
 8 8 8 24 29 55 60 76 88 95 98 100 103 104 104 118 131 137 136 132 125 104 102 100 97 95 93 88 80 54 40 20 18 8 8 8
 8 8 8 24 31 39 60 72 73 90 97 105 106 106 106 110 125 140 139 129 118 105 104 102 100 97 94 88 64 56 35 21 17 8 8 8
 8 8 8 8 23 41 70 70 74 85 97 102 103 104 105 107 109 111 112 109 107 102 100 97 98 95 89 83 64 53 29 22 8 8 8 8
 8 8 8 8 8 8 50 67 71 86 95 102 103 104 106 111 113 113 113 112 109 106 102 99 98 95 89 78 64 51 26 8 8 8 8 8
 8 8 8 8 8 8 40 74 77 87 94 98 104 107 109 112 113 114 114 112 108 105 101 98 97 90 84 69 60 50 20 8 8 8 8 8
 8 8 8 8 8 55 70 71 80 83 90 98 99 100 101 108 109 110 108 106 104 99 97 95 92 84 77 63 57 48 8 8 8 8 8 8
 8 8 8 8 8 61 62 72 78 82 85 97 98 101 107 109 111 111 111 110 106 100 96 95 94 82 76 67 58 49 8 8 8 8 8 8
 8 8 8 8 8 62 62 66 73 82 92 95 96 100 102 105 106 106 105 104 102 98 98 94 90 84 75 68 60 50 8 8 8 8 8 8
 8 8 8 8 8 55 62 75 78 90 96 98 99 101 102 104 108 108 107 104 100 96 94 92 89 83 73 69 64 55 8 8 8 8 8 8
 8 8 8 8 8 8 66 72 79 93 95 99 99 101 102 105 106 106 104 101 97 94 92 90 85 77 67 61 54 52 8 8 8 8 8 8
 8 8 8 8 8 8 8 67 89 98 99 105 110 110 110 112 110 110 108 107 106 105 101 96 90 81 86 77 70 50 8 8 8 8 8 8
 8 8 8 8 8 8 8 80 88 97 100 107 110 112 113 114 113 112 113 113 111 105 105 101 94 83 76 60 57 8 8 8 8 8 8 8
 8 8 8 8 8 8 8 8 63 85 101 106 109 112 113 114 114 114 114 113 112 110 108 102 99 90 80 64 62 8 8 8 8 8 8 8
 8 8 8 8 8 8 8 8 8 80 90 99 104 107 110 110 110 109 108 109 108 105 102 98 94 84 70 60 8 8 8 8 8 8 8 8
 8 8 8 8 8 8 8 8 8 80 91 100 103 107 108 109 108 107 107 106 104 101 99 95 90 80 66 8 8 8 8 8 8 8 8 8
 8 8 8 8 8 8 8 8 25 35 77 95 102 110 112 112 112 111 110 110 110 107 102 97 92 50 26 8 8 8 8 8 8 8 8 8
 8 8 8 8 8 8 8 8 8 20 34 40 90 96 100 101 106 106 107 104 102 99 97 85 40 30 30 8 8 8 8 8 8 8 8 8
 8 8 8 8 8 8 8 8 8 23 33 40 45 53 58 61 63 62 60 57 54 49 44 38 30 22 8 8 8 8 8 8 8 8 8 8
 8 8 8 8 8 8 8 8 8 15 27 35 41 48 54 57 58 57 56 54 50 45 41 33 23 8 8 8 8 8 8 8 8 8 8 8 ];

clf
set(gca,'visible','off')
text(.2,.5,'In the beginning, there was Abraham Lincoln...')
drawnow
pause

mesh(abe);
axis([0 40 0 50 -50 150]);
az = -135; el = 45;
view([az,el]);
title('Abe Lincoln as a Matrix (36 x 49)');
axis('equal')
drawnow
pause

xx=linspace(1,36,36)';
yy=linspace(1,49,49)';
[u,s,v]=svd(abe);
s = diag(s);
plot(xx,s,'y+')
title('Singular Values of Abe');
pause

plot(xx,s,'y+')
axis([0 37 0 400])
title('Closeup on Smaller Singular Values of Abe');
pause

clf
set(gca,'visible','off')
text(.4,.7,'For our first trick,')
text(.1,.6,'we decompose Abe, then look at some singular vectors.')
text(.1,.5,'The singular vectors are linearly interpolated for later')
text(.1,.4,'use in the "Abe Function" stage.')
drawnow
pause


n=18;
M = moviein(n);

for j=1:n
	plot(xx,v(:,j),'y-',xx,v(:,j),'g*')
	axis([0 37 -1 1]);
	txt=int2str(j);
	text(18,-.9,txt);
	txt=sprintf('Singular Value: %8.2f',s(j));
	text(11,.9,txt);
	title('Singular Vectors in the Shorter Dimension');
	M(:,j) = getframe;
end
movie(M,-2,1)
pause
% for j=1:n
% 	plot(yy,u(:,j),'y-',yy,u(:,j),'g*')
% 	axis([0 50 -1 1]);
% 	txt=int2str(j);
% 	text(25,-.9,txt);
% 	txt=sprintf('Singular Value: %8.2f',s(j));
% 	text(15,.9,txt);
%	title('Singular Vectors in the Longer Dimension')
% 	M(:,j) = getframe;
% end

clf
set(gca,'visible','off')
text(.2,.5,'For our next trick, we recompose Abe...')
drawnow
pause

z = zeros(49,36);
k = 1;
l = 0;
for j=1:36
	z = z + s(j)*u(:,j)*v(:,j)';
	if (j == k) 
		k = k + 2;
		l = l + 1;
		mesh(z);
		axis([0 40 0 50 -50 150]);
		axis('equal');
		txt=int2str(j);
		text(36,40,100,txt);
		view([az,el]);
		title('Abe Recomposed (as We Add Singular Vectors)');
		M(:,l) = getframe;
	end
end
movie(M,-2,1)
pause

clf
set(gca,'visible','off')
text(.1,.5,'Now we are going to spline-fit Abe to double his mesh...')
drawnow
pause

xxx=linspace(1,36,71)';
yyy=linspace(1,49,97)';
z = zeros(97,71);
for j=1:36
	spx = interp1(xx,v(:,j),xxx);
	spy = interp1(yy,u(:,j),yyy);
	z = z + s(j)*spy*spx';
end
mesh(z)
axis('equal')
az = -135; el = 45;
title('Abe Lincoln as a Matrix (71 x 97)');
view([az,el]);
spinmap
pause

% z = abe;
% az = -135; el = 45;
% k = 0;
% for j=36:-1:1
% 	mesh(pinv(z))
% 	axis([0 40 0 50 -.1 .1]);
% 	txt=int2str(j);
% 	text(36,40,100,txt);
% 	view([az,el]);
% 	k = k+1;
% 	M(:,k) = getframe;
% 	z = z - s(j)*u(:,j)*v(:,j)';
% end
% movie(M,-3,4)

mesh(abe/2);
% I added the /2 to scale abe so that the "axis('equal')" gives the same
% image.
title('Abe Lincoln as a Matrix (36 x 49)');axis('equal');
view([az,el]);spinmap
pause
mesh(z);title('Abe Lincoln as a Matrix (71 x 97)');axis('equal');
view([az,el]);spinmap
pause


clf
set(gca,'visible','off')
text(.5,.5,'The End!')
drawnow
pause

