%program SSLpopStats1.m
% SSLpopStats1 Load the SSLpopStats1.txt and conduct analyses on these data.

% Tomo Eguchi
% 25 February 2002

% Read the data file.
data = dlmread('d:\tomosFolders\classes\multiVariateSpring2002\SSLpopStats1.txt','\t');

% Separate them into different stats:
locations = data(:, 1);
indexNon0NP = find(abs(data(:,2))>eps);    % index for non-pups
rNPearlyAll = data(indexNon0NP, 2);        % Early instantaneous rates
rNPlateAll = data(indexNon0NP, 3);         % Late instantaneous rates

[indexNon0] = find(abs(data(:, 4)) > eps); % index for 18 locations with pup counts
rNPearly18 = data(indexNon0, 2);           % Early instantaneous rates for non-pups at 18 locs.
rNPlate18 = data(indexNon0, 3);            % Late instantaneous rates for non-pups at 18 locs.

pupData = data(indexNon0, 4:end);          % Pup data
pupLocations = data(indexNon0, 1);         % pup count locations

rPupsEarly18 = pupData(:, 1);              % Early instantaneous rates for pups.
rPupsLate18 = pupData(:, 2);               % Late instantaneous rates for pups.
NPpupsRatioEarly = pupData(:, 3);          % Early non-pup vs. pup ratio
NPpupsRatioMiddle = pupData(:, 4);         % Middle non-pup vs. pup ratio
NPpupsRatioLate = pupData(:, 5);           % Late non-pup vs. pup ratio

clf;
figure(1);                                 % replication of figure 3
plot(rNPearlyAll, rNPlateAll, 'o');
hold on;
for i = 1:length(rNPearlyAll),
   string = num2str(i);
   text(rNPearlyAll(i), rNPlateAll(i), string);
end

PCA1data = [rNPearlyAll, rNPlateAll];              % Do the PCA
[scores1, U1, lambda1, loadings1] = pca(PCA1data, 0);
varExplained1 = sum(lambda1)./(sum(sum(lambda1)));
Xbar = mean([rNPearlyAll, rNPlateAll]);
slopes = U1(:,2)./U1(:,1);                         % find the slopes of axes
refline(slopes(1), Xbar(2) - slopes(1).*Xbar(1));  % plot them
refline(slopes(2), Xbar(2) - slopes(2).*Xbar(1));
axis([-0.5, 0.15, -0.5, 0.15]);
vertline(0.0, 'r:');
horizline(0.0, 'r:');
refline(1, 0);
xlabel('Early Change');
ylabel('Late Change')

hold off;

figure(2);                                % plot scores
plot(scores1(:,1), scores1(:,2), 'o');
for i = 1:length(scores1(:,1)),
   string = num2str(i);
   text(scores1(i,1), scores1(i,2), string);
end
xlabel('PC 1');
ylabel('PC 2');
vertline(0.0, 'r:');
horizline(0.0, 'r:');

PCA2data = [rPupsEarly18, rPupsLate18, NPpupsRatioEarly, NPpupsRatioMiddle, NPpupsRatioLate];
[scores2, U2, lambda2, loadings2] = pca(PCA2data, 0);
varExplained2 = sum(lambda2)./(sum(sum(lambda2)));

figure(3);                                % replication of figure 4.
hold on;
plot(-scores2(:,1), -scores2(:,2), 'o');
for i = 1:length(scores2(:,1)),
   string = num2str(pupLocations(i));
   text(-scores2(i,1), -scores2(i,2), string);
end
xlabel('PC 1');
ylabel('PC 2');
vertline(0.0, 'r:');
horizline(0.0, 'r:');
hold off;

PCA3data = [rNPearly18, rNPlate18];
[scores3, U3, lambda3, loadings3] = pca(PCA3data, 0);
varExplained3 = sum(lambda3)./(sum(sum(lambda3)));

% combine Chernabura and Atkins (take means?)
index20 = find(pupLocations == 20);
CA = PCA2data(index20, :);
CAnew = mean(CA);
PCA4data = [PCA2data(1:(min(index20) - 1), :); CAnew; PCA2data((max(index20)+1):end, :)];
[scores4, U4, lambda4, loadings4] = pca(PCA4data, 0);
varExplained4 = sum(lambda4)./(sum(sum(lambda4)));

% load the fishery data
%fishData = load('d:\tomoFolders\classes\multiVariateSpring2002\SSLfishery.txt');