%
% FIGURE4.M - generate figure 4 from JoZ paper
%

% start from a clean sheet

clear all

% load data

load mammals.dat 
load dinosaurs.dat 

% train models

y_mammal   = log10(1000*mammals(:,end));
x_mammal   = log10(mammals(:,1) + mammals(:,2));
x_dinosaur = log10(dinosaurs);

% get on with it

w = [x_mammal(1:end-1) ones(size(x_mammal(1:end-1)))]\y_mammal(1:end-1);


X          = [0.5:0.01:+3]';
Y          = [X          ones(size(X))]*w;
Y_mammal   = [x_mammal   ones(size(x_mammal))]*w;
Y_dinosaur = [x_dinosaur ones(size(x_dinosaur))]*w;
S2         = var(y_mammal(1:end-1) - Y_mammal(1:end-1));
U          = Y+2*sqrt(S2);
L          = Y-2*sqrt(S2);

X          = 10.^X;
x_mammal   = 10.^x_mammal;
x_dinosaur = 10.^x_dinosaur;
Y          = 10.^Y;
y_mammal   = 10.^y_mammal;
Y_mammal   = 10.^Y_mammal;
Y_dinosaur = 10.^Y_dinosaur;
U          = 10.^U;
L          = 10.^L;
b          = w(1);
a          = 10.^w(2);

figure(1);
clf
set(axes, 'FontSize', 14);
plot(x_mammal, y_mammal, 'r+', ...
     X, Y, 'r-',...
     X, U, 'g--', ...
     X, L, 'g--');
axis([0 1000 -5000000 20000000]);
xlabel('Humerus + Femur (mm)');
ylabel('Body mass (g)');
legend('mammals', '\mu', '\mu\pm 2\sigma', 'Location', 'NorthWest');
text(100,7500000,sprintf('%f X^{%f}',a,b), 'FontSize', 14);
line([650 860], [0e+6 5.6e+6], 'Color', 'black');
text(500,-1500000, sprintf('elephant (%.1f metric tons)', Y_mammal(end)/1e+6), 'FontSize', 14);
line([300 575], [7e+6 3.5e+6], 'Color', 'black');
drawnow;
print -depsc2 fig4a.eps

y_mammal   = 1000*mammals(:,end);
x_mammal   = mammals(:,1) + mammals(:,2);
x_dinosaur = dinosaurs;

% get on with it

opts = optimset('Display', 'iter');
w    = fminunc(@loss, log2([2 2]), opts, x_mammal(1:end-1), y_mammal(1:end-1));
a    = 2.^w(1);
b    = 2.^w(2);

X          = [0.1:0.1:1000]';
Y          = a*X.^b;
Y_mammal   = a*x_mammal.^b;
Y_dinosaur = a*x_dinosaur.^b;
S2         = var(y_mammal(1:end-1) - Y_mammal(1:end-1));
U          = Y+2*sqrt(S2);
L          = Y-2*sqrt(S2);

figure(2);
clf
set(axes, 'FontSize', 14);
plot(x_mammal, y_mammal, 'r+', ...
     X, Y, 'r-',...
     X, U, 'g--', ...
     X, L, 'g--');
xlabel('Humerus + Femur (mm)');
ylabel('Body mass (g)');
text(200,35000000,sprintf('%f X^{%f}',a,b), 'FontSize', 14);
line([400 700], [3.3e+7 1.9e+7], 'Color', 'black');
text(500,-5000000, sprintf('elephant (%.1f metric tons)', Y_mammal(end)/1e+6), 'FontSize', 14);
line([700 850], [-0.1e+7 0.5e+7], 'Color', 'black');
legend('mammals', '\mu', '\mu \pm 2\sigma', 'Location', 'NorthWest');
drawnow
print -depsc2 fig4b.eps

% all done...

