swap(matrix,numswaps);

%-----------------------------------------------------------------------
%THIS FUNCTION CARRIES OUT A SWAP ALGORITHM AS DESCRIBED IN 
%STONE AND ROBERTS (1990), OECOLOGIA 85:74-79. THE DATA MATRIX MUST BE
%ORGANIZED WITH SITES (ISLANDS) AS ROWS, AND SPECIES AS COLUMNS. THE 
%FUNCTION ONLY WORKS WITH BINARY MATRICIES (BUT CAN BE EASILY MODIFIED
%TO CONVERT ABUNDANCE MATRICIES). 
%
%"Swapping" randomly shuffles a binary matrix while keeping the row and 
%column sums the same.
%
%  To call the function, send the data matrix and the number of swaps 
%	you wish to have carried out:  swap(matrix,1000). The 
%	function returns the shuffled matrix. For each swap, the function
%	randomly identifies two rows and columns such that:
%
%			[...0...1..;
%			 ..........;
%			 ...1...0..]     
%
%	and then switches (swaps) the 0's and 1's within the columns.
%
%Written by Brice X. Semmens ([e-mail]), 03/26/03
%Comments and modifications/improvements welcome. 
%-----------------------------------------------------------------------


function [matrix] = swap(matrix,numswaps);

sz = size(matrix); %row count in first cell, cols in second
swaptic=0;

while swaptic<numswaps %number of swaps to be carried out
   clear x1 & x2 & ledger & zero_ones & one_zeroes;
   
   x1 = ceil(sz(1,1) * rand); %randomly selects a row from data matrix
   x2 = ceil(sz(1,1) * rand); %randomly selects a row from data matrix
    
   %FIND 01 and 10 SPECIES
   ledger = matrix(x1,:)-matrix(x2,:);
       
   %IF THERE IS AT LEAST ONE OF EACH (01 and 10)...
   if max(ledger)==1 & min(ledger)==-1
      swaptic=swaptic+1;
      
      %INDEX 01 AND 10 CASES
      zero_ones = find(ledger==-1); %references all spps present in x1, not in x2
    	one_zeroes = find(ledger==1); %references all spps present in x2, not x1
  
  		%RANDOMLY SELECT A 01 AND 10 CASE
        pickZeroOnes=zero_ones(ceil(size(zero_ones,2)*rand));
        pickOneZeroes=one_zeroes(ceil(size(one_zeroes,2)*rand));
        
      %MODIFY MATRIX BY SWAPING CASES (turns 01 into 10, and vice versa)
      matrix(x1,pickZeroOnes)=1;
      matrix(x2,pickZeroOnes)=0;
      matrix(x1,pickOneZeroes)=0;
      matrix(x2,pickOneZeroes)=1;
   end
end