%%% This is the code used for the optimizations in the paper "A framework
%%% for mititaging tradeoffs between biodiversity conservation and
%%% agricultural production while minimizing species loss."  The code is
%%% written as a loop to limit the number of files.  In addition to the 5
%%% assembledges of species solved for in the paper we all include an assembledge
%%% (loop 1) that is all intermediate species, 
%%%(loop 7) 5 intermediate and
%%% 4 sharing species, (loop 8) five intermediate and 4 sparing species.
%%% The general setting in the opimization
%%% "'InitBarrierParam',.000000001,'MaxFunEvals',100000000,'MaxIter',100000000"
%%% works well for all assembledges.  However, a more precise answer for assembladge 5 can be obtained by increasing thses 
%%% parameters. sing a laptop computer running windows 7 home premium with the 64 bit operating systems, Intel Core i7-3520m CPU @
%%% 2.90GHZ and running Matlab R2013b, the full loop takes slightly over an
%%% hour. If values are changed to increase the precisions of assembladge
%%% 5, the time to solve takes much longer. The graph of assembledge 5
%%% shown in the manuscript took about 1 days time of processing. 
%%% Good initial values for each assemblege lead to much faster processing
%%% time. We thus seed the starting values differently for many of the
%%% assembledges. 
%%% 
tic

global a b0 b1 b2 maxbird bfmax a1 a2 a3 s maxbird maxyield a1 a2 a3 a4 a5 a6 a7 a8 a9 a10 e obi area c b g in in2 %%%Set global variabes to be call in all functions
makebspareshare9() %%% call function makebspareshare9 which contains the coeffcients of the species assembledges 
area()% set an open maxtix for total area
yy=zeros(1,10); %% create some matrixs that will be filled with results. 
h=ones(100,10);
q=ones(1);
  
profittab=ones(9,10); %% a matrix to store the results of yeild.
fname={'all inter','all share', 'all spare', '8 share 1 spare', '1 share 8 spare', 'all forms',' 5 inter 4 share', '5 inter 4 spare'};
%%names used for graphs
yname={'psallinter','psallshare', 'psallspare', 'ps8share1spare', 'ps1share8spare', 'psallforms','psintershare', 'psinterspare'};
%%% names used for graphs
dname={'psallinterprofit', 'psallshareprofit', 'psallspareprofit', 'ps8share1spareprofit', 'ps1share8spareprofit', 'psallformsprofit','psintershareprofit', 'ointerspareprofit'};
%%% names used for graphs
mname={'psallinteryield', 'psallshareyield', 'psllspareyield', 'ps8share1spareyield', 'ps1share8spareyield', 'psallformsyield','ointershareyield', 'psinterspareyield'};


a1=b; 
g1=g;
in1=in;
in22=in2;
x=zeros(100,1);%%% x is the yield of each parcel and goes from zero to 1 (we seed it with ones here).
for u=2:1:6 %%% this loops through the five species assembledges used in the paper
    %% a, in, g, in2 are coeffcients from makebspareshare9() and enters the density yeild curve 1-((x.^a(1))+(-in(1)*x.^g(1)))-in2(1))
%%where x is the yeild, a is a coeffcient, 
    a=a1(u*9-8:u*9,1);
    g=g1(u*9-8:u*9,1);
    in=in1(u*9-8:u*9,1);
    in2=in22(u*9-8:u*9,1);
    u
      x(1)=0;
for e=.1:.1:1;%%%% e is the min % of species on the landscape we simulate this from 10% to 100% (although there is zero production at 100% so no feasible solution exists. 
    e
if sum(x)~=100;   
if e==.1; %%% seeding the upper and lower bound for the first loop around the % of species populations
ub=ones(100,1); %%seed the lower bound
lb=zeros(100,1);%%%seed the upper bound
x0=ones(100,1)*.9;%%%seed the actual values
%X0(1:round((e+.1)*100),1)=x0(1:round((e+.1)*100)).*0;
else %%%Specific seed values for specific assemblages. 
    if u<3
x0=x;
    else
        if u==4
        x0=ones(100,1)*.9;%
        else
        if u==5
             x0=ones(100,1)*.5;
        else
     
x0=x*.99;
   end
        end
    end
end
ub=ones(100,1); 
lb=zeros(100,1);
%%% here we start the actual optimization usign the fmincon function in the
%%% optimization toolbox. 
options = optimoptions('fmincon','Algorithm','interior-point','InitBarrierParam',.000000001,'MaxFunEvals',1000000,'MaxIter',1000000);
[x, fval,exitflag] = ...
fmincon('profit',x0,[],[],[],[],lb,ub,'consconsletter1',options);
z=x;
yield1=sum(x); %% here we calculate the yield after solving for x.

%%% now we reset the starting values and run the optimization. We found the
%%% solutions to be a bit dependent on the starting guesses so we cycle
%%% through two starting points to be sure we find the best solution. 
x1=ones(100,1).*e;% Make a starting guess at the solution
ub=ones(100,1);
lb=zeros(100,1);
options = optimoptions('fmincon','Algorithm','interior-point','InitBarrierParam',.000000001,'MaxFunEvals',1000000,'MaxIter',1000000);
[x, fval,exitflag] = ...0
fmincon('profit',x1,[],[],[],[],lb,ub,'consconsletter1',options);
z2=x;
yield2=sum(x); %% calculate yield for the second solution

if yield1>yield2; %%% choose the solution with the higher yield. 
    x=z;
else
    x=z2;
end
if exitflag==-2||exitflag==0; % if the model does not solve we make sure that the program skips the rest of the loop in order to save time. 
    k=ones(100,1);
    x=k;
else
k=sort(x);
end
else
    k=ones(100,1);
    x=k;
end
if x==ones(100,1);
    yield=0;
else
yield=sum(x); % calculate the total yeild 
end



yy(1,round(e*10))=yield; %%% save yield
h(:,round(e*10))=k;%%% save solutions
end
profittab(u,:)=yy; 
p=gca; %%% graph optimal landscape:i.e. what is shown in figure 3. 
set(p,'XLim',[0 11])
set(p,'yLim',[0 1])
h(h==1)=nan;
 fig=figure(round(e*10+60));set(fig,'Position',[0 0 500 500]);
      bar(h');title ([fname{u}]);
      xlabel('species target %(x 10)')
      ylabel('yield')
        axis([0,10,0,1])
      myStyle.Format='jpeg';
      myStyle.Resolution = 700;
      hgexport(fig,[yname{u}],myStyle);
%%% graph yeild (not shown in paper)
     profitfig=figure(round(e*10+70));set( profitfig,'Position',[0 0 500 500]);
      bar(yy');title ([fname{u}]);
      xlabel('species target % (x 10)')
      ylabel('yield')
        axis([0,10,0,100])
      myStyle.Format='jpeg';
      myStyle.Resolution = 700;
      hgexport(profitfig,[dname{u}],myStyle);
   
      v=ones(100,9);
   for j=1:1:9  %%% Graph yield density curves for the assembledge
   for n=0.01:.01:1
    q=round(n*100);
    v(q,j)=1-((n.^a(j))+(-in(j)*(n.^g(j))))-in2(j);
   end
    end

    yieldfig=figure(round(e*10+80));set(yieldfig,'Position',[0 0 500 500]);
    plot(v);title ([fname{u}]);
    xlabel('Yield %')
    ylabel('Bird Density')
    hgexport(yieldfig,[mname{u}],myStyle);
  
end
save ssprofit9scale
toc

