clear all
global el n k tk N r0 w0 a0 h0 b0 d0 m0 c0 A0 H0 B0 D0 M0 C0 AA0 HH0 BB0 DD0 MM0 CC0 p10 p20 p30
global r w a h b d m c A H B D M C AA HH BB DD MM CC p1 p2 p3
%return
parnames={'el' 'n' 'k' 'tk' 'N' 'r0' 'w0' 'a0' 'h0' 'b0' 'm0' 'd0' 'c0'...
    'A0' 'H0' 'B0' 'D0' 'M0' 'C0' 'AA0' 'HH0' 'BB0' 'DD0' 'MM0' 'CC0' 'q10' 'q20' 'p10' 'p20' 'p30'};
%method={'??'};

% Specify the dimensions of the food web.
el=4; % Number of resources
n=100;% Max. possible number of plants
k=200;  % Max. possible number of herbivores
tk=80;%tk=50;  % Max. possible number of carnivores ('top killers').
N=60; % Number of simulation runs. 
%
% Open the file 'parameterZeros' to specify the means of parameter values.
%
parameterZeros;


% Specify the varying parameters for which the foodweb data of model "BioDiversityModel4" is collected for
% "N" many runs.

partialr0=[1 2 1e+3]; %'partialr0(1)'=1 iff the first 'partialr0(2)' is enriched with the other fixed at 'partialr0(3)'.
chngpars={'r0', 'MM0', 'CC0'}; % Name the varying parameters for the data. 
logscales={'yes','no', 'no'}; % Using log scale ('yes') on the varying parameters. 
chngparranges=[1e-2 1e+5; MM0 2*MM0; CC0 2*CC0]; % Specify the ranges of the varying parameters.
                                           % data generated for
                                           % 'partialr0(1)=0' with [1e-1 1e+7; ... ];
chngparpoints=[60; 1; 1]; % Specify numbers of regularly partitioned points for the varying parameters.
randparpts={'no','no','no'}; % Specify if the discrete parameter points should be randomized 
                               % around the regular partition ('yes' or 'no').
popthrhd=[1e-4, 1e-2, 1e-0, 1e+1, 1e+2]; %Specify any population size levels to detect. 

dataname='Data.mat'; % Save the data to 'dataname'.

Tini=0;  % Time initial point for the numercial solutions.
Tend=250; % Time end point for the numercial solutions.
Ttrn=.9*Tend; % Time transient for the numercial solutions.

options=odeset('RelTol',1e-5,'AbsTol',1e-8);

%%
%%%% For the ode solver of this problem, the default relative and absolute
%%%% tolerances are (1e-3 and 1e-6, respectively) 
% It takes about 6 minutes to go N=1 run for r0pts=60. Or 15 hrs for N=50.
%



%%%%%%%%%%%%
%%%%%%%%%%%%  Run the whole script when you reach this point. 
%%%%%%%%%%%%

par0values=[el n k tk N r0 w0 a0 h0 b0 d0 m0 c0 A0 H0 B0 D0 M0 C0 AA0 HH0 BB0 DD0 MM0 CC0 q10 q20 p10 p20 p30];

startclk=clock;
OutData=[];
for i0=1:N
    
[el1 n1 k1 tk1 N1 r1 w1 a1 h1 b1 d1 m1 c1 A1 H1 B1 D1 M1 C1 AA1 HH1 BB1 DD1 MM1 CC1 q11 q21 p11 p21 p31]=UnitParValues(par0values);

[el n k tk N r0 w0 a0 h0 b0 d0 m0 c0 A0 H0 B0 D0 M0 C0 AA0 HH0 BB0 DD0 MM0 CC0 q10 q20 p10 p20 p30]=Par0Values(par0values);
par0values=[el n k tk N r0 w0 a0 h0 b0 d0 m0 c0 A0 H0 B0 D0 M0 C0 AA0 HH0 BB0 DD0 MM0 CC0 q10 q20 p10 p20 p30];

%%%%%%%%%
%%%%%%%%% Discretize parameters to get par1pts and par2pts. 
%%%%%%%%%

DiscretizeParameters;

%%%%%%%%%
%%%%%%%%%
%%%%%%%%%
xinit=[(r0>0)*(ones(1,el)) .1*min(popthrhd)*(ones(1,n+k+tk))];
        
i1cnt=0;            
for i1=par1pts%i1=chngparranges(1,1):inc(1):chngparranges(1,2)
i1cnt=i1cnt+1;
id1=strcmp(chngpars(1), parnames);
par0values(id1)=i1;
i2cnt=0;
    for i2=par2pts%i2=chngparranges(2,1):inc(2):chngparranges(2,2)
        i2cnt=i2cnt+1;
        id2=strcmp(chngpars(2), parnames);
        par0values(id2)=i2;
        i3cnt=0;
        for i3=par3pts
            i3cnt=i3cnt+1;
            id3=strcmp(chngpars(3), parnames);
            par0values(id3)=i3;
            
            progress=[i0 i1cnt i2cnt i3cnt]
yy=[];
zz=[];
pp=[];
    
    [el n k tk N r0 w0 a0 h0 b0 d0 m0 c0 A0 H0 B0 D0 M0 C0 AA0 HH0 BB0 DD0 MM0 CC0 q10 q20 p10 p20 p30]=Par0Values(par0values);
    par0values=[el n k tk N r0 w0 a0 h0 b0 d0 m0 c0 A0 H0 B0 D0 M0 C0 AA0 HH0 BB0 DD0 MM0 CC0 q10 q20 p10 p20 p30];

    [el n k tk N r w a h b d m c A H B D M C AA HH BB DD MM CC q1 q2 p1 p2 p3]=ParValues(par0values,...
        el1, n1, k1, tk1, N1, r1, w1, a1, h1, b1, d1, m1, c1, A1, H1, B1, D1, M1, C1, AA1, HH1, BB1, DD1, MM1, CC1, q11, q21, p11, p21, p31);

    if partialr0(1)==1
        r(partialr0(2)+1:end)=partialr0(3)*r1(partialr0(2)+1:end);
    end
    
[t,x]=ode15s('BioDiversityModel4',[0,Tend],xinit,options);


%  xinit=x(end,:);

[tplot bbb]=find(t>Ttrn);
yy=min(x(tplot,(el+1):(el+n)),[],1);
zz=min(x(tplot,(el+n+1):(el+n+k)),[],1);
pp=min(x(tplot,(el+n+k+1):(el+n+k+tk)),[],1);

 %speciesnumbers=[sum(yy(yy>popthrhd)), sum((yy>popthrhd)), ...
 %sum((zz>popthrhd)), sum(pp>popthrhd)]

        %rsrmass=sum(x(end,1:el)); 
        %pltmass=sum(yy(yy>popthrhd));
        %hbvrmass=sum(zz(zz>popthrhd));
        %cnvrmass=sum(pp(pp>popthrhd));
        %pltnmbr=sum((yy>popthrhd));
        %hbvrnmbr=sum((zz>popthrhd));
        %cnvrnmbr=sum((pp>popthrhd));
        
        rsrmass=sum(x(end,1:el));
        
        Pltmass=[];
        Hbvrmass=[];
        Cnvrmass=[];
        Pltnmbr=[];
        Hbvrnmbr=[];
        Cnvrnmbr=[];
        for i4=1:length(popthrhd)
            pidx=(yy>popthrhd(i4));
            hidx=(zz>popthrhd(i4));
            cidx=(pp>popthrhd(i4));
            Pltmass=[Pltmass sum(yy(pidx))];
            Hbvrmass=[Hbvrmass sum(zz(hidx))];
            Cnvrmass=[Cnvrmass sum(pp(cidx))];
            Pltnmbr=[Pltnmbr sum(pidx)];
            Hbvrnmbr=[Hbvrnmbr sum(hidx)];
            Cnvrnmbr=[Cnvrnmbr sum(cidx)];
        end
        %pltmass=sum(yy(yy>popthrhd));
        %hbvrmass=sum(zz(zz>popthrhd));
        %cnvrmass=sum(pp(pp>popthrhd));
        %pltnmbr=sum((yy>popthrhd));
        %hbvrnmbr=sum((zz>popthrhd));
        %cnvrnmbr=sum((pp>popthrhd));
        xplot=x(tplot,:);
            %poprate=[];
            livingpoprate=[];
            endangeredspecies=[];
            competitivespecies=[];
            for i=1:length(tplot)
                derivativeofstate=BioDiversityModel4(tplot(i),xplot(i,:)',[]);
                %poprate=[poprate derivativeofstate((el+1):end)];
                %pcptr=PerCapitaRate4(tplot(i),xplot(i,:)',[]);
                %percapitagrate=[percapitagrate pcptr((el+1):end)];
                percapitagrate=PerCapitaRate4(tplot(i),xplot(i,:)',[]);
                %meanpercapitagrate=mean(percapitagrate');
                noncompetitivespecies=(percapitagrate((el+1):end)'<-1e-4);
                livingspecies=(xplot(i,(el+1):end)>1e-4);
                competitivespecies=[competitivespecies;(percapitagrate((el+1):end)'>1e-4).*(1-livingspecies)];
                livingpoprate=[livingpoprate derivativeofstate((el+1):end).*livingspecies'];
                %selectionindx=sum(abs((1-noncompetitivespecies)-livingspecies));%/(n+k+tk);
                endangeredspecies=[endangeredspecies; noncompetitivespecies.*livingspecies];
            end

            %meanplantpoprate=mean(mean(abs(poprate(1:n,:)')));
            %meanherbivorepoprate=mean(mean(abs(poprate((n+1):(n+k),:)')));
            %meancarnivorepoprate=mean(mean(abs(poprate((n+k+1):(n+k+tk),:)')));

            %maxplantpoprate=max(mean(abs(poprate(1:n,:)')));
            %maxherbivorepoprate=max(mean(abs(poprate((n+1):(n+k),:)')));
            %maxcarnivorepoprate=max(mean(abs(poprate((n+k+1):(n+k+tk),:)')));

            endangeredspecies=mean(endangeredspecies,1);
            endangeredplants=sum(endangeredspecies(1:n));
            endangeredherbivores=sum(endangeredspecies((n+1):(n+k)));
            endangeredcarnivores=sum(endangeredspecies((n+k+1):(n+k+tk)));

            meancompetitivespecies=mean(competitivespecies,1);
            competitveplants=sum(meancompetitivespecies(1:n));
            competitveherbivores=sum(meancompetitivespecies((n+1):(n+k)));
            competitvecarnivores=sum(meancompetitivespecies((n+k+1):(n+k+tk)));

            meanlivingpoprate=mean(mean(abs(livingpoprate'),1));
            maxpoprate=max(mean(abs(livingpoprate'),1));
            
        OutData=[OutData; i0 i1 i2 i3 rsrmass Pltmass Hbvrmass Cnvrmass Pltnmbr Hbvrnmbr Cnvrnmbr endangeredplants endangeredherbivores endangeredcarnivores competitveplants competitveherbivores competitvecarnivores meanlivingpoprate maxpoprate];
            %meanplantpoprate meanherbivorepoprate meancarnivorepoprate ...
            %maxplantpoprate maxherbivorepoprate maxcarnivorepoprate ...
            %endangeredplants endangeredherbivores endangeredcarnivores meanpoprate maxpoprate];
        
     save(dataname,'chngpars','partialr0','parnames','par0values','popthrhd','OutData'); 
        end

    end
    
end

end

if N<5
    OutData(:,end-4:end)
end
endclk=clock;

startclk(4:6)

endclk(4:6)

save('Databkup.mat','chngpars','partialr0','parnames','par0values','popthrhd','OutData'); 
%

OutDatabkup=OutData;
        
return
        
        
        
        
        
figure(10)
hold off

[yvlu,indexy]=max(max(x(tplot,(el+1):(el+n)),[],1));
[zvlu,indexz]=max(max(x(tplot,(el+n+1):(el+n+k)),[],1));
[pvlu,indexp]=max(max(x(tplot,(el+n+k+1):(el+n+k+tk)),[],1));

plot(t,x(:,1)+x(:,2)+x(:,3)+x(:,4),'g-', t,x(:,el+indexy(1)),'b-', t,x(:,el+n+indexz(1)),'m-',...
    t,x(:,el+n+k+indexp(1)),'r-')

axis([0 Tend 0 max([yvlu,zvlu,pvlu])])

