% Script 2 for week 1

% First a little puzzle, then its solution (with documentation along the way.)

% "A person who can, within a year, solve x^2-92*y^2=1 is a mathematician." 
%  Quote from Bramagupta, The Opening of the Universe, 628 AD.

% Suppose you want to find some integer solutions (which is what Bramagupta meant) to Bramagupta's diophantine equation x^2-92*y^2=1. Here we find solutions.

% First create a BIG vector of integers.  This little problem is an advertisement for vector methods of the kind built into MATLAB, and the philosophy of problem solving they encourage.

for n=1:10^7
x(n)=n;
end

% So integer j is in the jth coordinate of x for coordinates up to ten million. The y values of any solution will have to be just a LITTLE smaller than x/sqrt(92).

y=floor(x/sqrt(92));
z=x.^2-92*y.^2-1;

% Focus attention on the x values (which here ARE the indices) we want. Keep an eye on the workspace.

x=find(z==0)

% Now find the associated y values 

y=floor(x/sqrt(92))

% Above, we have solved Bramagupta's problem in 7 lines of MATLAB code.

% check:

x.^2-92*y.^2

% Interesting note. If you let x go to 1 billion you get no more solutions. If you do try that, be prepared to wait a few minutes for the output.

% Now, we go on to "more about graphing." What is the difference between the following two methods of creating a graph? (Answer: One uses vector methods, the other doesn't.) 

clearvars
x=0:0.1:1
y=sin(x);
figure(20)
plot(x,y,x,y,'*r','markersize',12)
clearvars


for i=1:11
t(i)=(i-1)*0.1;
z(i)=t(i)^2;
end
figure(21)
plot(t,z,t,z,'*r','markersize',12)


% We now look at how to read and write comma-separated files from variables which can be created/read by excel or any db program.


D=[pi,pi^2,pi^3;0,1,2];

format long
D
format short
D
format long

dlmwrite('report.txt',D);

% Open and examine the textfile report.txt, which will be in MATLAB's Current Folder.

jj=dlmread('report.txt')
jj==D

format short
D
jj

% MATLAB thinks the second rows of jj are the same as the second rows of D, but thinks the first rows differ.

% To get more precision try the following. You can specify the column delimiter and the number of digits (though the last few digits in report.txt entries, in this case, will be meaningless.) MATLAB will now think D and jj are indistinguishable.

dlmwrite('report.txt', D, 'delimiter', '\t', 'precision', 20)

jj=dlmread('report.txt')
jj==D

format long

D
jj

format short

% Understanding some things about the internal representation of numbers in MATLAB (or ANY finite precision arithmetic tool) is important and the following bears further thought. 


x=0.3-0.1;
y=0.5-0.3;

x==y

% Why do you think MATLAB thinks they are different? (Note: The BINARY representations of decimal 0.1 is 0.000110011001100110011... with 0011 repeating forever. But the binary representation of decimal 0.5 is exactly 0.1.)

x-y
eps(x)
eps(y)
eps(x-y)


% Look at the file report.txt after doing these next commands. 
% What do you think is happening? (How many digits does MATLAB preserve?)


d=12345678912345678912345678912345678912345678912345
f=12345678912345678912000000000000000000000000000000
eps(d)

d+1e33==d
d+3e33==d

g=f+1
h=d-f
dlmwrite('report.txt', [d;f;g;h], 'precision', 60)

% So the key idea is the following. MATLAB (or any finite precision tool) has an eps for each number which varies depending on its size but for MATLAB represents about 16 decimal places to the right of the first significant digit. MATLAB cannot distinguish numbers that differ by less than the eps of the bigger number.


% Here is a key idea: 
% Suppose you have a million different numbers OF WIDELY VARYING MAGNITUDES whose sum SHOULD, based on some physical theory, represent the value of a measurement. In order to compare actual measurement to theory you want to add these numbers up. Given the option, how should you carry out the addition?

% Should you add them from largest to smallest, or from smallest to largest?

format short
clear
a=10^-5*ones(1,10^8+1);
b=a;
a(1)=10^11;
b(10^8+1)=10^11;

max(a)
max(b)

min(a)
min(b)


% a and b are vectors containing the same 10^8+1 numbers but in reverse order. Now let's add them up from left to right in each case.


sum(a)
sum(b)
sum(a)-sum(b)


% Think about this! A large fraction of the things we do here correspond to adding up many numbers with widely varying "epsilons"!!! It could possibly make a significant difference depending on how you add them!

% On to a new topic.... fun with matrices!

A=[3 5;8 7;2 7]

transpose(A)

A'

% A' is conjugate-transpose, the same as transpose(A) unless there are complex number entries.

hh=A(:,1)

kk=A(2,:)

A

fliplr(A)

flipud(A)

A

A([1 3],:)


% Read and write files (and save an entire MATLAB environment) for internal MATLAB use. If you want to shut down operations for the day but come back with EVERYTHING just like now, use this.

save('hhhhh')

clearvars

load('hhhhh')

save('ggggg','A')

clearvars A

load('ggggg')

clearvars

% More on graphing!!!

t=0:0.01:1;

figure(42)
plot(t,t, '-..r','LineWidth',2)
hold on
plot(t,t.^2,'g','LineWidth',2)
legend('y=x','y=x^2')
legend('location','northwest')
xlabel('Time (Sec)')
ylabel('Distance (m)')
title('Graphing Exercise','FontSize',20,'FontWeight','bold')
hold off


figure(43)
hold on
subplot(1,2,1)
plot(t,t, '-.r','LineWidth',2)

% The decimals in a text command represent fraction of total graph size 
% from bottom left at which text is placed.

text(.35,0.3,'\leftarrow  eek','FontSize',16,'FontWeight','bold')
xlabel('Time (years)')
ylabel('Distance (yards)')
title('Y=X','FontSize',20,'FontWeight','bold')
subplot(1,2,2)
plot(t,t.^2, '--g','LineWidth',2)
text(.3,.8,'wow','FontSize',16,'FontWeight','bold')
xlabel('Time (sec)')
ylabel('Distance (m)')
title('Y=X^2','FontSize',20,'FontWeight','bold')
hold off


A=[3 5 7; 9 3 6; 3 2 1; 8 8 9]
t=[6 7 8 9]
figure(45)
plot(A)
figure(46)
plot(t,A)


% Here is a (much more complex and varied) example illustrating a variety of things you might want to do. First, input the data to be plotted - I found this example somewhere on the interweb.

TBdata = [1990 4889 16.4; 1991 5273 17.4; 1992 5382 17.4; 1993 5173 16.5;
          1994 4860 15.4; 1995 4675 14.7; 1996 4313 13.5; 1997 4059 12.5;
          1998 3855 11.7; 1999 3608 10.8; 2000 3297  9.7; 2001 3332  9.6;
          2002 3169  9.0; 2003 3227  9.0; 2004 2989  8.2; 2005 2903  7.9;
          2006 2779  7.4; 2007 2725  7.2];

measles = [38556 24472 14556 18060 19549 8122 28541 7880 3283 4135 7953 1884]';
mumps = [20178 23536 34561 37395 36072 32237 18597 9408 6005 6268 8963 13882]';
chickenPox = [37140 32169 37533 39103 33244 23269 16737 5411 3435 6052 12825 23332]';
years = TBdata(:, 1);
cases = TBdata(:, 2);
rate  = TBdata(:, 3);

% Create the pie chart in position 1 of a 2x2 grid

figure(47)

subplot(2, 2, 1)
pie([sum(measles) sum(mumps) sum(chickenPox)], {'Measles', 'Mumps', 'Chicken Pox'})
title('Childhood Diseases')

% Create the bar chart in position 2 of a 2x2 grid
subplot(2, 2, 2)
bar(1:12, [measles/1000 mumps/1000 chickenPox/1000], 0.5, 'stack')
xlabel('Month')
ylabel('Cases (in thousands)')
title('Childhood Diseases')
axis([0 13 0 100])
set(gca, 'XTick', 1:12)

% Create the stem chart in position 3 of a 2x2 grid
subplot(2, 2, 3)
stem(years, cases)
xlabel('Years')
ylabel('Cases')
title('Tuberculosis Cases')
axis([1988 2009 0 6000])

% Create the line plot in position 4 of a 2x2 grid
subplot(2, 2, 4)
plot(years, rate)
xlabel('Years')
ylabel('Infection Rate')
title('Tuberculosis Cases')
axis([1988 2009 5 20])


% Below we find the time taken to solve a typical matrix nxn equation A*x=b using A\b and inv(A)*b, compared using tic and toc. 

% Here we create a cell array of random matrices to work on. A cell array is like a matrix but its entries can be ANYTHING: strings, numbers, logicals or, in this case, matrices. 

% And cell arrays can have any number of dimensions, not just 2 like a matrix has.

% First we create empty 1x1000 cell arrays. It is more efficient to do this up-front rather than to let MATLAB deduce what you must be doing when you start naming cell entries. Also, if somewhere else you have already created a LONGER cell array with this name, this will chop it down to size. (You can do the same thing with a matrix M by starting with M=zeros(m,n) to initialize M as a matrix of the specified size with all zero entries.)

% We then populate these empty cell arrays with randomly generated vectors and matrices. Then we solve the equation A*x=b for 1000 different matrix dimmensions in two different ways and compare the time it takes as the size of the matrices increases. This will take a minute or two.

% What do the two graphs indicate, what are its important features?

A = cell(1,1000);
b = cell(1,1000);
for n = 1:1000;
    A{n} = rand(n,n);
    b{n} = rand(n,1);
    end
    
t = zeros(1,1000);
for n = 1:1000
    tic;
    x = A{n}\b{n};
    t(n) = toc;
end

s = zeros(1,1000);
for n = 1:1000
    tic;
    x = inv(A{n})*b{n};
    s(n) = toc;
end


figure(8)
plot(s)
hold on
plot(t)
hold off
legend('Solve by Inverse Matrix','Solve by Backslash Operator')
title('Solve Matrix Equation, Two Methods')