%%Mixing Polygon Simulation (3 isotopes) (v.1.0; 15/2/2013) 
%J. A. Smith et. al. 2013; The University of New South Wales
%Created using MATLAB v7.10.0.499(2010a)
%For updates and instructions, visit http://www.famer.unsw.edu.au/downloads.html

clear all;
close all;

%USER INPUT
mixture = xlsread('Mixture_example_3i.xls'); %NOTE: use .xls files, column headers only
sources = xlsread('Sources_example_3i.xls'); %NOTE: every source must have some error
TEF = xlsread('TEF_example_3i.xls');
its = 1500;  %specify the number of iterations

%RUN the simulation
[mrows,mcols] = size(mixture);
[srows,scols] = size(sources);
p = zeros(its,mrows);
par_vals = zeros(its,srows*6+3);
[prows,pcols] = size(par_vals);
c_p = pcols - 2;
c_i = pcols - 1;

for i = 1:its
    disp(['Run Number ',num2str(i),''])
    vC = zeros(srows,1);
    vN = zeros(srows,1);
    vS = zeros(srows,1);
    fC = zeros(srows,1);
    fN = zeros(srows,1);
    fS = zeros(srows,1);
    for j = 1:srows
        vC(j) = normrnd(sources(j,1),sources(j,2)); %isotope 1 (e.g. d13C)
        vN(j) = normrnd(sources(j,3),sources(j,4)); %isotope 2 (e.g. d15N)
        vS(j) = normrnd(sources(j,5),sources(j,6)); %isotope 3 (e.g. d34S)
        fC(j) = normrnd(TEF(j,1),TEF(j,2)); %enrichment for iso 1
        fN(j) = normrnd(TEF(j,3),TEF(j,4)); %enrichment for iso 2
        fS(j) = normrnd(TEF(j,5),TEF(j,6)); %enrichment for iso 3
    end
    VC = vC + fC;
    VN = vN + fN;
    VS = vS + fS;
    V = [VC, VN, VS];
    [hulln, vol] = convhulln(V); %gives indices of the hull (triangulation)
    hulln_i = unique(hulln);
    indices = V(hulln_i,:);
    indices_d = delaunayn(indices); %computes simplices using only the convex hull's vertices
    P = tsearchn(indices,indices_d,mixture);
    P = P';
    P(P > 0) = 1;
    P(isnan(P)) = 0;
    p(i,:) = P;
    
    vals = [vC;vN;vS;fC;fN;fS;0;0;0]';
    par_vals(i,:) = vals;
    par_vals(i,c_p) = vol;
    par_vals(i,c_i) = i;
    par_vals(i,pcols) = var(par_vals(1:i,c_p));
end

%FIGURE 1: variance of iterated polyhedron volumes; evaluates number of iterations
Variance = par_vals(:,pcols);
Iteration = par_vals(:,c_i);
figure
plot(Iteration,Variance)
xlabel('Iterations','FontSize',12)
ylabel('Variance polyhedron volume','FontSize',12)
 
%FIGURE 2: proportion of iterations that each consumer was inside mixing polygon
Probabilities = sum(p)/its
figure
bar(Probabilities); grid;
xlabel('Consumer','FontSize',12)
ylabel('Probability consumer is within mixing polyhedron','FontSize',12)
set(gca,'FontName','Ariel','FontSize',12,'LineWidth',1)

pp = [p ; Probabilities];
xlswrite('Consumer_Probabilities.xls',p,1,'A1')
Lets = {'Poly_Vol' 'Iteration' 'Variance'};
Headers = [strcat({'Source'},int2str(repmat(1:srows,1,6).')).' Lets];
xlswrite('Parameter_values.xls', Headers, 1, 'A1') 
xlswrite('Parameter_values.xls', par_vals, 1, 'A2') %stores isotope, TEF, volume, iteration, and variance data