% Written by Carmen García-Comas
%for the article 
%'Automatic Image Processing to Determine the Community Size Structure of Riverine Macroinvertebrates'
% Rosa Gurí, Ignasi Arranz, Marc Ordeix, Carmen García-Comas
% published in JoVE 2022


clear all;close all;
let={'a','b','c','d','e','f','g','h','i','j'};

T=readtable('Id_from_sorted_vignettes_20220408_1215'); % output table from the validated zooscan
samples=unique(T(:,2));
G=[1,1,1,1,0,1,0,1,1,0,1,0,1,1,0,1,0,1,0,1,0,1,1,0,1,0,0,0,1,0,1,1,0,1,1,1,0,0,1,1,0,1,1,0,1]'; %large fraction
F=[0,0,0,0,1,0,1,0,0,1,0,1,0,0,1,0,1,0,1,0,1,0,0,1,0,1,1,1,0,1,0,0,1,0,0,0,1,1,0,0,1,0,0,1,0]'; %short fraction
samples=table2array(samples);

N=NaN(length(samples),3);Ntot=NaN(length(samples),1);
Nauto=N;
Recall=N;Precision=N;ellipMauto=N;ellipMmanual=N;ELLIPVOLauto={};ELLIPVOLmanual={};ellipSDauto=N;ellipSDmanual=N;
category1=unique(table2cell(T(:,14)));
category2=unique(table2cell(T(:,15)));

category=unique([category1;category2]);

N2=NaN(length(samples),length(category));
Recall2=N2;Precision2=N2;ellipMauto2=N2;ellipMmanual2=N2;ELLIPVOLauto2={};ELLIPVOLmanual2={};ellipSDauto2=N2;ellipSDmanual2=N2;

cat={'macroinv','other','out'};

 % to convert pixel to um
    
    Resol=2400;
    microm=25.4/Resol*1000; %pixel to um
    converse= microm*1E-3; %um to mm
     
% size classes

ellip=[];
for i=1:length(samples)
M=table2array(T(:,9));m=table2array(T(:,10));
   M=M*converse;m=m*converse;
   ellip= [ellip;(4/3)*pi.*(M./2).*((m./2).^2)]; 
end
mIn=nanmin(ellip); 
MIn=nanmax(ellip);
bins(1)=mIn;
for i=2:30
      bins(i)=bins(i-1)*2;
end

  bins=bins(bins<=MIn); 
  mbins=(bins(1:end-1)+bins(2:end))./2;
  
for i=1:length(samples)
    
   f=strmatch(samples{i},table2cell(T(:,2)),'exact');
   
   auto=table2cell(T(f,14));
   manual=table2cell(T(f,15)); 
   M=table2array(T(f,9));m=table2array(T(f,10));
   M=M*converse;m=m*converse;
   ellip= (4/3)*pi.*(M./2).*((m./2).^2); %formula for the ellipsoidal volume
   
   for j=1:length(cat)
     
       f1=strmatch(cat{j},auto);
       f2=strmatch(cat{j},manual);
       
      Nauto(i,j)=length(f1);
      N(i,j)=length(f2);Ntot(i)=length(auto);
      
       a=ismember(f1,f2);
       if isempty(a)
           TP=0;
       else
     
       TP=sum(a);
       end
       
       b=ismember(f2,f1);
       if isempty(b)
           FN=0;
       else
      
       FN=length(find(b==0));
       end
       
       c=find(a==0);
       if isempty(c)
           FP=0;
       else
      
       FP=length(c);
       end
       Recall(i,j)=TP/(TP+FN);
       Precision(i,j)=TP/(TP+FP);   
       
       % size and size indicators
       if isempty(f)==0
       ELLIPVOLauto{i,j}=ellip(f1);
       ellipMauto(i,j)=nanmean(ellip(f1));
       if length(f1)>1
       ellipSDauto(i,j)=size_diversity_kde(log2(ellip(f1)));
       end
       end
       if isempty(f2)==0
       ELLIPVOLmanual{i,j}=ellip(f2);
       ellipMmanual(i,j)=nanmean(ellip(f2));
       if length(f2)>1
       ellipSDmanual(i,j)=size_diversity_kde(log2(ellip(f2)));
       end
       end
   end
    
  
   for j=1:length(category)
     
       
       f1=strmatch(category{j},auto,'exact');
       
       f2=strmatch(category{j},manual,'exact');
       
     
      N2(i,j)=length(f2);
      
       a=ismember(f1,f2);
       if isempty(a)
           TP=0;
       else
      
       TP=sum(a);%TP (classified by the machine and you as that category)
       end
      
       b=ismember(f2,f1);
       if isempty(b)
           FN=0;
       else
       
       FN=length(find(b==0)); %FN (classified by the machine as something else)
       end
       
       c=find(a==0);
       if isempty(c)
           FP=0;
       else
      
       FP=length(c);
       end 
       Recall2(i,j)=TP/(TP+FN);
       Precision2(i,j)=TP/(TP+FP);   
       
       % size and size indicators
       
       if isempty(f)==0
       ELLIPVOLauto2{i,j}=ellip(f1);
       ellipMauto2(i,j)=nanmean(ellip(f1));
       if length(f1)>1
       ellipSDauto2(i,j)=size_diversity_kde(log2(ellip(f1)));
       end
       end
       if isempty(f2)==0
       ELLIPVOLmanual2{i,j}=ellip(f2);
       ellipMmanual2(i,j)=nanmean(ellip(f2));
       if length(f2)>1
       ellipSDmanual2(i,j)=size_diversity_kde(log2(ellip(f2)));
       end
       end
   end
   
  
end

Recall(N==0)=NaN;
Precision(N==0)=NaN;
Recall2(N2==0)=NaN;
Precision2(N2==0)=NaN;

MacroinvAuto=NaN(1000,length(samples));InvAuto=MacroinvAuto;
MacroinvMan=MacroinvAuto;InvMan=MacroinvAuto;
NBSSauto=NaN(length(mbins),size(ELLIPVOLauto,1));ASSauto=NBSSauto;
NBSSmanual=NaN(length(mbins),size(ELLIPVOLauto,1));ASSmanual=NBSSmanual;
MODEauto=NaN(size(ELLIPVOLauto,1),1);MODEmanual=MODEauto;

for i=1:size(ELLIPVOLauto,1)
a=ELLIPVOLauto{i,1};
MacroinvAuto(1:length(a),i)=a;
b=ELLIPVOLauto{i,2};a=[a;b];
InvAuto(1:length(a),i)=a;
    for j=1:length(bins)-1
         if j<=length(bins)-1
         x=find(((a>=bins(j)).*(a<bins(j+1)))==1);   % a >= limite inferior y < limiter superior de la clase de talla 
         NBSSauto(j,i)=sum(a(x));
         ASSauto(j,i)=length(x);
         end
    end
    
     f=find(NBSSauto(:,i)./mbins'==max(NBSSauto(:,i)./mbins'));
     MODEauto(i)=f;
     
a=ELLIPVOLmanual{i,1};
MacroinvMan(1:length(a),i)=a;
b=ELLIPVOLmanual{i,2};a=[a;b];
InvMan(1:length(a),i)=a;

    for j=1:length(bins)-1
         if j<=length(bins)-1
         x=find(((a>=bins(j)).*(a<bins(j+1)))==1);   % ev >= limite inferior y < limiter superior de la clase de talla 
         NBSSmanual(i,j)=sum(a(x));
         ASSmanual(i,j)=length(x);
         end
    end
    
     f=find(NBSSauto(:,i)./mbins'==max(NBSSauto(:,i)./mbins'));
     MODEauto(i)=f;
end

%%%%%%%%%%% FIGURES %%%%%%%%%%%%%%%

% FIGURE 5: Recall & 1-Precision of the selected learningset

close all
figure
set(gcf, 'PaperUnits', 'centimeters','papersize',[10 8])
set(gcf,'paperposition',[0 0 10 8])
a=boxplot(Recall,'notch','on','color','k');
set(a,'linewidth',1)
hold on
a=boxplot(1-Precision,'notch','on','color','r');
set(a,'linewidth',1)
set(gca,'linewidth',1,'xtick',[1,2,3],'xticklabel',{'Macro inv.','Other inv.','Debris'},'ytick',[0:0.25:1],'FontName','arial','fontsize',10,'linewidth',1);
grid on
a=xlabel('Categories');
set(a,'FontName','arial','fontsize',11);
a=ylabel('Predictive value');
set(a,'FontName','arial','fontsize',11);
a=text(0.65,1.1,'Recall');
set(a,'FontName','arial','fontsize',11,'color','k');
a=text(0.65,0.4,'1-Precision');
set(a,'FontName','arial','fontsize',11,'color','r');

print ('Gurietal_JoVE_fig5','-dpdf','-r300')

% FIGURE 7: individual size distributions

close all

figure
set(gcf, 'PaperUnits', 'centimeters','papersize',[15 11])
set(gcf, 'PaperPosition', [0 0 15 11]);
pos=find(G==1);
col=jet(length(pos));
subplot(2,2,1)
for i=1:length(pos)
    x=MacroinvAuto(:,pos(i));[f,xi] = ksdensity(log2(x)); plot(xi,f,'-','color',col(i,:),'linewidth',1);
    hold on
end
lim=[-15,25];
set(gca,'ylim',[0,0.8],'xlim',lim,'ytick',[0:0.2:1],'linewidth',1);
a=xlabel('log_2(mm^3)');
set(a,'FontName','arial','fontsize',11,'color','k');
a=ylabel('Density');
set(a,'FontName','arial','fontsize',11,'color','k');
grid on
a=title('Automatic coarse');
set(a,'FontName','arial','fontsize',12,'color','k');
a=text(lim(1),0.9,let{1});
set(a,'FontName','arial','fontsize',12);
subplot(2,2,2)
for i=1:length(pos)
    x=MacroinvMan(:,pos(i));[f,xi] = ksdensity(log2(x)); plot(xi,f,'-','color',col(i,:),'linewidth',1);
    hold on
end
set(gca,'ylim',[0,0.8],'xlim',lim,'ytick',[0:0.2:1],'linewidth',1);
a=xlabel('log_2(mm^3)');
set(a,'FontName','arial','fontsize',11,'color','k');
a=ylabel('Density');
set(a,'FontName','arial','fontsize',11,'color','k');
grid on
a=title('Validated coarse');
set(a,'FontName','arial','fontsize',12,'color','k');
a=text(lim(1),0.9,let{2});
set(a,'FontName','arial','fontsize',12);
pos=find(F==1);
col=jet(length(pos));
subplot(2,2,3)
lim=[-11,11];
for i=1:length(pos)
    x=MacroinvAuto(:,pos(i));[f,xi] = ksdensity(log2(x)); plot(xi,f,'-','color',col(i,:),'linewidth',1);
    hold on
end
set(gca,'ylim',[0,0.8],'xlim',lim,'ytick',[0:0.2:1],'linewidth',1);
a=xlabel('log_2(mm^3)');
set(a,'FontName','arial','fontsize',11,'color','k');
a=ylabel('Density');
set(a,'FontName','arial','fontsize',11,'color','k');
grid on
a=title('Automatic fine');
set(a,'FontName','arial','fontsize',12,'color','k');
a=text(lim(1),0.9,let{3});
set(a,'FontName','arial','fontsize',12);
subplot(2,2,4)
for i=1:length(pos)
    x=MacroinvMan(:,pos(i));[f,xi] = ksdensity(log2(x)); plot(xi,f,'-','color',col(i,:),'linewidth',1);
    hold on
end
set(gca,'ylim',[0,0.8],'xlim',lim,'ytick',[0:0.2:1],'linewidth',1);
a=xlabel('log_2(mm^3)');
set(a,'FontName','arial','fontsize',11,'color','k');
a=ylabel('Density');
set(a,'FontName','arial','fontsize',11,'color','k');
grid on
a=title('Validated fine');
set(a,'FontName','arial','fontsize',12,'color','k');
a=text(lim(1),0.9,let{4});
set(a,'FontName','arial','fontsize',12);
%print ('fig2','-dtiff','-r300')
print ('Gurietal_JoVE_fig7','-dpdf','-r300')

% FIGURE 6: Nauto vs Nval EllipVol auto vs EllipVol val coarse & fine samples
figure
set(gcf, 'PaperUnits', 'centimeters','papersize',[25 10])
set(gcf, 'PaperPosition', [0 0 25 10]);
pos=find(G==1);
subplot(1,2,1)
plot(Nauto(G==1,1),N(G==1,1),'ko','markerfacecolor','k')
hold on
plot(Nauto(F==1,1),N(F==1,1),'ko','markerfacecolor',[0.9,0.9,0.9])
hold on
plot([0,max([xlim,ylim])],[0,max([xlim,ylim])],'k--')
a=xlabel('automatic N ind. counted');
set(a,'FontName','arial','fontsize',11,'color','k');
a=ylabel('validated N ind. counted');
set(a,'FontName','arial','fontsize',11,'color','k');
grid on
% P = polyfit(Nauto(G==1,1),N(G==1,1),1)
% P =
%     0.9565   -0.0745
a=text(50,720,'\mu=0.85 (r=0.99)');
set(a,'FontName','arial','fontsize',11);
a=text(50,630,'coarse \mu=0.75 (r=0.92)');
set(a,'FontName','arial','fontsize',11);
a=text(50,540,'fine \mu=0.84 (r=0.98)');
set(a,'FontName','arial','fontsize',11);
 a=legend('coarse','fine','location','northwest');
 set(a,'FontName','arial','fontsize',11);
set(gca,'ytick',[0:200:max([xlim,ylim])],'xtick',[0:200:max([xlim,ylim])],'FontName','arial','fontsize',11);
a=text(0,1050,let{1});
set(a,'FontName','arial','fontsize',12);
subplot(1,2,2)
plot(ellipMauto(G==1,1),ellipMmanual(G==1,1),'ko','markerfacecolor','k')
hold on
plot(ellipMauto(F==1,1),ellipMauto(F==1,1),'ko','markerfacecolor',[0.9,0.9,0.9])
hold on
plot([0,max([xlim,ylim])],[0,max([xlim,ylim])],'k--')
a=xlabel('automatic mean ellip. vol. (mm^3)');
set(a,'FontName','arial','fontsize',11,'color','k');
a=ylabel('validated mean ellip. vol. (mm^3)');
set(a,'FontName','arial','fontsize',11,'color','k');
set(gca,'ytick',[0:150:max([xlim,ylim])],'xtick',[0:150:max([xlim,ylim])],'FontName','arial','fontsize',11);
grid on
%corr(ellipMauto(F==1,1),ellipMmanual(F==1,1))
%P = polyfit(ellipMauto(F==1,1),ellipMmanual(F==1,1),1)
% ans =
% 
%     0.998
% 
% corr(ellipMauto(G==1,1),ellipMmanual(G==1,1))
% 
% ans =
% 
%     0.913
% 
% corr(ellipMauto(:,1),ellipMmanual(:,1))
% 
% ans =
% 
%     0.987
a=text(45,520,'\mu=1.27 (r=0.99)');
set(a,'FontName','arial','fontsize',11);
a=text(45,450,'coarse \mu=2.24 (r=0.91)');
set(a,'FontName','arial','fontsize',11);
a=text(45,380,'fine \mu=1 (r=1)');
set(a,'FontName','arial','fontsize',11);
a=text(0,733,let{2});
set(a,'FontName','arial','fontsize',12);

print ('Gurietal_JoVEfig6','-dpdf','-r300')

% FIGURE 9: size spectra ind size distributions & size diversity estimates 

T2=readtable('Id_from_sorted_vignettes_20220509_1705');
samples2=unique(T2(:,2));
samples2=table2array(samples2);
fraction=[8;1;16;16;1;16;1;16;1;8;1];%from abundancies_2 excel
cat=unique(table2cell(T2(:,15)));
%concatenate size vectors
N3=NaN(length(samples2),length(cat));ELLIPVOL3={};ellipM3=N3;
for i=1:length(samples2)
    
   f=strmatch(samples2{i},table2cell(T2(:,2)),'exact');

   manual=table2cell(T2(f,15)); 
   
   M=table2array(T2(f,9));m=table2array(T2(f,10));
   M=M*converse;m=m*converse;
   ellip= (4/3)*pi.*(M./2).*((m./2).^2); %formula for the ellipsoidal volume

   for j=1:length(cat)
      f2=strmatch(cat{j},manual,'exact');
      N3(i,j)=length(f2);
      % size and size indicators
      if isempty(f2)==0
           ELLIPVOL3{i,j}=repmat(ellip(f2),fraction(i),1);
           ellipM3(i,j)=nanmean(ellip(f2));
           
      end
   end
   f2=strmatch(samples2{i},samples,'exact');
   match(i)=f2;
   ellipvolInv{i,1}=repmat(ELLIPVOLauto{f2,1},fraction(i),1);
   ellipvolInv{i,2}=repmat(ELLIPVOLmanual{f2,1},fraction(i),1);
    Ninv(i,1)=N(f2,1);
end
n=0;
posF=[1,3,6,8,10]; n=0;
for i=1:length(posF)
    n=n+1;
    if posF(i)==3
        indsizeAuto{n,:}= [ellipvolInv{i,1}; ellipvolInv{i+1,1}; ellipvolInv{i+2,1}];
        indsizeMan{n,:}=  [ellipvolInv{i,2}; ellipvolInv{i+1,2}; ellipvolInv{i+2,2}];
        indsizeManNoM{n,:}= [ELLIPVOL3{i,1};ELLIPVOL3{i+1,1};ELLIPVOL3{i+2,1}];
    else
        indsizeAuto{n,:}= [ellipvolInv{i,1}; ellipvolInv{i+1,1}];
        indsizeMan{n,:}=  [ellipvolInv{i,2}; ellipvolInv{i+1,2}];
        indsizeManNoM{n,:}= [ELLIPVOL3{i,1};ELLIPVOL3{i+1,1}];
    end
end
Nind=NaN(5,3);mS=Nind;
SD=NaN(5,3);

figure
set(gcf, 'PaperUnits', 'centimeters','papersize',[20 30])
set(gcf, 'PaperPosition', [0 0 20 30]);
Sd1={'2.98';'2.65';'2.63';'2.62';'2.82'};
Sd2={'3';'2.62';'2.60';'2.56';'2.87'};
Sd3={'3';'2.60';'2.59';'2.57';'2.72'};
for i=1:5
    subplot(3,2,i)
 [f,xi] = ksdensity(log2(indsizeAuto{i})); plot(xi,f,'-','color',[1,0,0],'linewidth',1);
hold on
[f,xi] = ksdensity(log2(indsizeMan{i})); plot(xi,f,'-','color',[0,0,0.5],'linewidth',1);
hold on
[f,xi] = ksdensity(log2(indsizeManNoM{i})); plot(xi,f,'-','color',[0,0.7,0],'linewidth',1);

a=xlabel('log_2(mm^3)');
set(a,'FontName','arial','fontsize',12,'color','k');
a=ylabel('Density');
set(a,'FontName','arial','fontsize',12,'color','k');
if i==1
   a=legend('Automatic','Validated','Corrected') ;
   set(a,'FontName','arial','fontsize',10);
end
grid on

SD(i,1)=size_diversity_kde(log2(indsizeAuto{i}));
SD(i,2)=size_diversity_kde(log2(indsizeMan{i}));
SD(i,3)=size_diversity_kde(log2(indsizeManNoM{i}));
Nind(i,1)=length(indsizeAuto{i});
Nind(i,2)=length(indsizeMan{i});
Nind(i,3)=length(indsizeManNoM{i});
mS(i,1)=nanmean(indsizeAuto{i});
mS(i,2)=nanmean(indsizeMan{i});
mS(i,3)=nanmean(indsizeManNoM{i});
set(gca,'ylim',[0,0.33],'ytick',[0:0.1:0.3],'xlim',[-10,11],'xtick',[-10:5:10],'FontName','arial','fontsize',12,'linewidth',1);
a=text(-10,0.35,let{i});
set(a,'FontName','arial','fontsize',12);
a=text(-8,0.29,['Sd=',Sd1{i}],'color','r');
set(a,'FontName','arial','fontsize',12);
a=text(-8,0.25,['Sd=',Sd2{i}],'color',[0,0,0.5]);
set(a,'FontName','arial','fontsize',12);
a=text(-8,0.21,['Sd=',Sd3{i}],'color',[0,0.7,0]);
set(a,'FontName','arial','fontsize',12);
end
subplot(3,2,6)
plot(SD(:,2),SD(:,3),'ko','markerfacecolor','k');
hold on
plot([min([xlim,ylim]),max([xlim,ylim])],[min([xlim,ylim]),max([xlim,ylim])],'k--')
a=xlabel('Validated size diversity');
set(a,'FontName','arial','fontsize',12,'color','k');
a=ylabel('Corrected size diversity');
set(a,'FontName','arial','fontsize',12,'color','k');
grid on
%P = polyfit(SD(:,2),SD(:,3),1))
a=text(2.55,2.9,'\mu=0.85 (r=0.94)');
set(a,'FontName','arial','fontsize',12);
set(gca,'xlim',[2.5,3.05],'ylim',[2.5,3.05],'xtick',[2.5:0.1:3],'ytick',[2.5:0.1:3],'FontName','arial','fontsize',12,'linewidth',1);
a=text(2.5,3.08,let{6});
set(a,'FontName','arial','fontsize',12);
print ('Gurietal_JoVEfig9','-dpdf','-r300')
close all

%Figure 8: N validated v N separated mean size
figure
set(gcf, 'PaperUnits', 'centimeters','papersize',[25 10])
set(gcf, 'PaperPosition', [0 0 25 10]);
pos=find(G==1);
subplot(1,2,1)
plot(Nind(:,2),Nind(:,3),'ko','markerfacecolor','k')
hold on
set(gca,'ytick',[0:10^4:max([xlim,ylim])],'xtick',[0:10^4:max([xlim,ylim])],'ylim',[min([xlim,ylim]),max([xlim,ylim])],'xlim',[min([xlim,ylim]),max([xlim,ylim])],'FontName','arial','fontsize',11);
plot([0,max([xlim,ylim])],[0,max([xlim,ylim])],'k--')
a=xlabel('validated N individuals frame^-^1');
set(a,'FontName','arial','fontsize',11,'color','k');
a=ylabel('separated N individuals frame^-^1');
set(a,'FontName','arial','fontsize',11,'color','k');
grid on
a=text(5000,32000,'\mu=1.38');
set(a,'FontName','arial','fontsize',11,'color','k');
a=text(5000,28000,' r=1');
set(a,'FontName','arial','fontsize',11,'color','k');
set(gca,'ytick',[0:10^4:max([xlim,ylim])],'xtick',[0:10^4:max([xlim,ylim])],'ylim',[min([xlim,ylim]),max([xlim,ylim])],'xlim',[min([xlim,ylim]),max([xlim,ylim])],'FontName','arial','fontsize',11);
a=text(-2000,37000,let{1});
set(a,'FontName','arial','fontsize',12);
subplot(1,2,2)
plot(mS(:,2),mS(:,3),'ko','markerfacecolor','k')
hold on
plot([0,max([xlim,ylim])],[0,max([xlim,ylim])],'k--')
a=xlabel('validated mean ellip. vol. (mm^3)');
set(a,'FontName','arial','fontsize',11,'color','k');
a=ylabel('separated mean ellip. vol. (mm^3)');
set(a,'FontName','arial','fontsize',11,'color','k');
set(gca,'ytick',[0:2:max([xlim,ylim])],'xtick',[0:2:max([xlim,ylim])],'FontName','arial','fontsize',11);
grid on
% P = polyfit(mS(:,2),mS(:,3),1)
% 
% P =
% 
%     0.9565   -0.0745
% 
% corr(mS(:,2),mS(:,3))
% 
% ans =
% 
%     0.9977
a=text(1,7.2,'\mu=0.96');
set(a,'FontName','arial','fontsize',11,'color','k');
a=text(1,6.4,' r=1');
set(a,'FontName','arial','fontsize',11,'color','k');
a=text(-0.5,8.5,let{2});
set(a,'FontName','arial','fontsize',12);
%print ('fig4','-dtiff','-r300')
print ('Gurietal_JoVEfig8','-dpdf','-r300')

%create size classes
ellip=[];
for i=1:length(indsizeMan)
    ellip=[ellip;indsizeMan{i}];
end   
    
mIn=quantile( ellip ,0.01 );%nanmin(ellip); % doing this, the size classes are too thin at the beginning and too large at the end and thus the spectra have an inverse size spectrum (positive size class)
MIn=quantile( ellip ,0.99 )%nanmax(ellip);
bins(1)=mIn;
for i=2:30
      bins(i)=bins(i-1)*2;
end
mbins=(bins(1:end-1)+bins(2:end))./2;
wbins=(bins(2:end)-bins(1:end-1));
NBSSa=NaN(length(mbins),size(indsizeAuto,1));NBSSm=NBSSa;NBSSma=NBSSa;
ASSa=NBSSa;ASSm=NBSSa;ASSma=NBSSa;
MODEa=NaN(size(indsizeAuto,1),1);MODEm=MODEa;MODEms=MODEa;
for i=1:size(indsizeAuto,1)
a=indsizeAuto{i};b=indsizeMan{i};c=indsizeManNoM{i};
    for j=1:length(bins)-1
         if j<=length(bins)-1
         x=find(((a>=bins(j)).*(a<bins(j+1)))==1);   % a >= low limit & < upper limit of the size class 
         NBSSa(j,i)=sum(a(x));
         ASSa(j,i)=length(x);
         x=find(((b>=bins(j)).*(b<bins(j+1)))==1);   
         NBSSm(j,i)=sum(b(x));
         ASSm(j,i)=length(x);
         x=find(((c>=bins(j)).*(c<bins(j+1)))==1); 
         NBSSms(j,i)=sum(c(x));
         ASSms(j,i)=length(x);
         end
    end
    
     f=find(NBSSa(:,i)./wbins'==max(NBSSa(:,i)./wbins'));
     MODEa(i)=f;
     f=find(NBSSm(:,i)./wbins'==max(NBSSm(:,i)./wbins'));
     MODEm(i)=f;
     f=find(NBSSms(:,i)./wbins'==max(NBSSms(:,i)./wbins'));
     MODEms(i)=f;

end

NBSSa=NBSSa./wbins';
NBSSm=NBSSm./wbins';
NBSSms=NBSSms./wbins';

moda=mode([MODEa;MODEm;MODEms]);

%    -0.8794   -0.8833   -0.9037
%    -1.2431   -1.2002   -1.2425
%    -1.0684   -1.0684   -1.1242
%    -1.2314   -1.2276   -1.2992
%    -0.9213   -0.9164   -0.9102


%     0.9345    0.9335    0.9313
%     0.9366    0.9336    0.9579
%     0.9608    0.9279    0.9477
%     0.9601    0.9174    0.9004
%     0.8617    0.8624    0.8226
   
Mu1={'-0.88','-1.24','-1.07','-1.23','-0.92'};
Mu2={'-0.88','-1.20','-1.07','-1.23','-0.92'};
Mu3={'-0.90','-1.24','-1.12','-1.30','-0.91'};

R_1={'0.93','0.94','0.96','0.96','0.86'};
R_2={'0.93','0.94','0.93','0.92','0.86'};
R_3={'0.93','0.96','0.95','0.90','0.82'};

% Figure 10: NBSS & Slope

figure
set(gcf, 'PaperUnits', 'centimeters','papersize',[20 30])
set(gcf, 'PaperPosition', [0 0 20 30]);
for i=1:5
    subplot(3,2,i)
plot(log2(mbins), log2(NBSSa(:,i)),'X','color',[1,0,0],'markersize',10,'linewidth',1);
hold on
plot(log2(mbins), log2(NBSSm(:,i)),'^','color',[0,0,0.5],'markersize',10,'linewidth',1);
hold on
plot(log2(mbins), log2(NBSSms(:,i)),'o','color',[0,0.7,0],'markersize',10,'linewidth',1);
hold on
plot([log2(mbins(moda)),log2(mbins(moda))],[0,15],'k--','linewidth',1);
a=xlabel('log_2(mm^3)');
set(a,'FontName','arial','fontsize',12,'color','k');
a=ylabel('log_2(frame^-^1)');
set(a,'FontName','arial','fontsize',12,'color','k');
if i==1
   a=legend('Automatic','Validated','Corrected','location','southwest') ;
   set(a,'FontName','arial','fontsize',10);
end
grid on
set(gca,'ylim',[0,15],'ytick',[0:5:15],'xlim',[-5,10],'xtick',[-5:5:10],'FontName','arial','fontsize',12,'linewidth',1);
a=log2(mbins(moda:end))';
b=log2(NBSSa(moda:end,i));
a=a(isinf(b)==0);b=b(isinf(b)==0);
 P = polyfit(a,b,1);
 
Slope(i,1)=P(1);
 R2(i,1)=corr(a,b)^2;
a=log2(mbins(moda:end))';
b=log2(NBSSm(moda:end,i));
a=a(isinf(b)==0);b=b(isinf(b)==0);
 P = polyfit(a,b,1);

Slope(i,2)=P(1);
R2(i,2)=corr(a,b)^2;
a=log2(mbins(moda:end))';
b=log2(NBSSms(moda:end,i));
a=a(isinf(b)==0);b=b(isinf(b)==0);
 P = polyfit(a,b,1);
 
Slope(i,3)=P(1);
R2(i,3)=corr(a,b)^2;

a=text(-7,17,let{i});
set(a,'FontName','arial','fontsize',12);


a=text(3.45,14,['\mu=',Mu1{i},'(r^2=',R_1{i},')'],'color','r');
set(a,'FontName','arial','fontsize',10);
a=text(3.45,12.5,['\mu=',Mu2{i},'(r^2=',R_2{i},')'],'color',[0,0,0.5]);
set(a,'FontName','arial','fontsize',10);
a=text(3.45,11,['\mu=',Mu3{i},'(r^2=',R_3{i},')'],'color',[0,0.7,0]);
set(a,'FontName','arial','fontsize',10);   
end
subplot(3,2,6)
plot(Slope(:,2),Slope(:,3),'ko','markerfacecolor','k');
hold on
plot([min([xlim,ylim]),max([xlim,ylim])],[min([xlim,ylim]),max([xlim,ylim])],'k--')
a=xlabel('Validated slope');
set(a,'FontName','arial','fontsize',12,'color','k');
a=ylabel('Corrected slope');
set(a,'FontName','arial','fontsize',12,'color','k');
grid on
%P = polyfit(Slope(:,2),Slope(:,3),1)
a=text(-1.25,-0.9,'\mu=1.15 (r=0.99)');
set(a,'FontName','arial','fontsize',12);
set(gca,'ylim',[-1.3,-0.8],'ytick',[-1.3:0.2:-0.8],'xlim',[-1.3,-0.8],'xtick',[-1.3:0.2:-0.8],'FontName','arial','fontsize',12,'linewidth',1);

a=text(-1.34,-0.75,let{6});
set(a,'FontName','arial','fontsize',12);
print ('Gurietal_JoVEfig10','-dpdf','-r300')


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Size diversity function:

%written by Ye Lin for Ye et al. 2016:
% Increasing zooplankton size diversity enhances the strength of top-down control on phytoplankton through diet niche partitioning
% Lin Ye, Chun-Yi Chang, Carmen García-Comas, Gwo-Ching Gong and Chih-hao Hsieh
% doi: 10.1111/1365-2656.12067 

function u=size_diversity_kde(data,n,MIN,MAX)
% Calculate the size diversity based on
% Quintana, X.D., Brucet S., et al. 2008. A nonparametric method for the
% measurement of size diversity with emphasis on data standardization.
% Limnology and Oceanography: Methods, 6, 75-86

if nargin<2 % if n is not supplied switch to the default
    [bandwidth,density,xmesh]=kde(data);
else
    [bandwidth,density,xmesh]=kde(data,n,MIN,MAX);
end



y=zeros(1,length(density));
for i=1:length(density)
    if density(i)>0
        y(i)=density(i)*log2(density(i));
    else
        y(i)=0;
    end
end

% Using trapezoidal integration
u=-1*trapz(xmesh,y);
% if u<0
%    % u=-1*trapz(flip(xmesh),flip(y)); %added by carmen to deal with spurious negative size diversity (?) integrate from right to left?
% % Using Monte Carlo integration
% end
end

