%The code is constructed for a 2 x 4 landscape with spatially dependent
%benefits.  The rows are labeled A and B in the paper and the columns are
%labeled 1 through 4 in the paper.  A letter-number combination, for example, A4, 
%gives the parcel's address on the map.  

%The C matrix gives conservation costs.

%The B1 matrix gives the conservation benefit (b) when no neighboring
%parcel is conserved.

%The B2 matrix gives the conservation benefit (b) when one neighboring
%parcel is conserved.

%The B3 matrix gives the conservation benefit (b) when two neighboring
%parcels are conserved.

%The B4 matrix gives the conservation benefit (b) when three neighboring
%parcels are conserved.

%To solve the spatailly-independent problem define B1 abd then set B1=B2=B3=B4.


C = [3 1 3 3; 1 1 1 3];  %User input
B1 = [6 5 4 2; 1 3 5 6]; %User input
B2 = [9 8 5 5; 2 6 8 9];  %User input
B3 = [11 10 7 7; 3 8 10 11];  %User input
B4 = [0 11 9 0; 0 9 11 0];  %User input

V = 0.25; %User input


%Find optimal landscape. (Find X-star-i) vector, net benefit, etc.  This
%calls the funciton 'findoptimal.m'
conserveoption=ones(2,4);
[NB,BSumFinal,Pattern,FinalB]=findoptimal(C,B1,B2,B3,B4,V,conserveoption);
OptNB=NB; OptBSum=BSumFinal; OptPattern=Pattern; OptB=FinalB; 

%OptPattern gives a value of '1' in a cell if the parcel is optimally conserved and
%a 0 otherwise.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%Find W(X-star-i) for each conserved i.
for j=1:2; for k=1:4;
    index=ones(2,4); index(j,k)=0; BV(j,k)=BSumFinal-sum(sum(C.*Pattern.*index));    
end; end;
W=BV.*Pattern;  % The matrix 'W' gives the values of "W(X-star-i)".
                % The value for each parcel is given at the parcel's location
                % on the landscape.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%Find X-star-~i for each conserved i
count = 0;
Wnoti=zeros(2,4);
for j=1:2; for k=1:4;
    if OptPattern(j,k)==1
       count = count + 1;
       conserveoption=ones(2,4);
       conserveoption(j,k)=0;
       [NB,BSumFinal,Pattern,FinalB]=findoptimal(C,B1,B2,B3,B4,V,conserveoption);
       OptNBnoti(count,1)=NB; OptBSumnoti(count,1)=BSumFinal; OptPatternnoti(((count-1)*2)+1:count*2,1:4)=Pattern; OptBnoti(((count-1)*2)+1:count*2,1:4)=FinalB;
       
      %Find W(X-star-~i) for each conserved i. 
      Wnoti(j,k)=BSumFinal-sum(sum(C.*Pattern));  % The matrix 'Wnoti' gives the values of "W(X-star-~i)".
                                                  % The value for each parcel is given at the parcel's location
                                                  % on the landscape.
    end
end; end


%Calculate Delta-W(i)
deltaWi = W - Wnoti;