%
% FIGURE1 - generate figure 1 (a)-(d) for JoZ paper
%

% start from a clean sheet

clear all

% load data

load mammals.dat 
load dinosaurs.dat 

% first the conventional allometric models

y_mammal   = log10(1000*mammals(:,end));
x_mammal   = log10(mammals(:,1) + mammals(:,2));
x_dinosaur = log10(dinosaurs);
w          = [x_mammal ones(size(x_mammal))]\y_mammal;

% generate predictions

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 - Y_mammal);
U          = Y+2*sqrt(S2);
L          = Y-2*sqrt(S2);
b          = w(1);
a          = w(2);

% plot some graphs

figure(1);
clf;
set(axes, 'FontSize', 14);
plot(x_mammal, y_mammal, 'r+', ...
     X, Y, 'r-',...
     X, U, 'g--', ...
     X, L, 'g--');
xlabel('log_{10}(Humerus + Femur)');
ylabel('log10 body mass');
text(1.5,1,sprintf('log Y = %.4f + %.4f log X',a,b), 'FontSize', 14);
legend('mammals', '\mu', '\mu\pm 2\sigma', 'Location', 'NorthWest');
line([1.45 1.15], [1.1 2.00], 'Color', 'black');
drawnow;
print -depsc2 fig1a.eps

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(2);
clf
set(axes, 'FontSize', 14);
plot(x_mammal, y_mammal, 'r+', ...
     X, Y, 'r-',...
     X, U, 'g--', ...
     X, L, 'g--');
ax = axis;
axis(ax - [0 0 5000000 0]);
text(200,1e+7,sprintf('%f X^{%f}',a,b), 'FontSize', 14);
line([380 600], [0.95e+7 0.4e+7], 'Color', 'black');
text(500,-2.5e+6, sprintf('elephant (%.1f metric tons)',Y_mammal(end)/1e+6), 'FontSize', 14);
line([650 850], [-0.1e+7 0.5e+7], 'Color', 'black');
xlabel('Humerus + Femur (mm)');
ylabel('Body mass (g)');
legend('mammals', '\mu', '\mu\pm 2\sigma', 'Location', 'NorthWest');
drawnow;
print -depsc2 fig1b.eps

% now for non-linear regression models 

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, y_mammal);
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 - Y_mammal);
U          = Y+2*sqrt(S2);
L          = Y-2*sqrt(S2);

figure(3);
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,5e+6,sprintf('%.4f X^{%.4f}',a,b), 'FontSize', 14);
line([375 575], [4.5e+6 2.6e+6], 'Color', 'black');
text(500,-0.5e+6, sprintf('elephant (%.1f metric tons)',Y_mammal(end)/1e+6), 'FontSize', 14);
line([675 865], [0e+6 5.6e+6], 'Color', 'black');
legend('mammals', '\mu', '\mu \pm 2\sigma', 'Location', 'NorthWest');
drawnow;
print -depsc2 fig1d.eps

X          = log10(X);
x_mammal   = log10(x_mammal);
x_dinosaur = log10(x_dinosaur);
Y          = log10(Y);
y_mammal   = log10(y_mammal);
Y_dinosaur = log10(Y_dinosaur);
U          = log10(U);
L(L<0)     = -10;
L(L>0)     = log10(L(L>0));

figure(4);
clf;
set(axes, 'FontSize', 14);
plot(x_mammal, y_mammal, 'r+', ...
     X, Y, 'r-',...
     X, U, 'g--', ...
     X, L, 'g--');
axis([0.5 3 1 8]);
xlabel('log_{10}(Humerus + Femur)');
ylabel('log_{10} body mass');
text(1.5,1.5,sprintf('log Y = %.4f + %.4f log X',a,b), 'FontSize', 14);
line([1.45 1], [1.6 2.5], 'Color', 'black');
legend('mammals', '\mu', '\mu \pm 2\sigma', 'Location', 'NorthWest');
drawnow;
print -depsc2 fig1c.eps

% all done...

