# torus_2011.txt
# Kyle E. Harms, 1/14/2011

torusallspp=function(habmat, allabund, species.list, plotdim=c(1000,500), gridsize=20)
 {
  plotdimqx=plotdim[1]/gridsize  # Calculates no. of x-axis quadrats of plot.
  plotdimqy=plotdim[2]/gridsize  # Calculates no. of y-axis quadrats of plot.
  num.habs=length(unique(as.vector(habmat)))  # Determines tot. no. of habitat types.
  sppno=length(species.list)  # Determines tot. no. of spp. for which to conduct habitat analysis.
  GrLsEq=matrix(0,sppno,num.habs*4)  # Creates empty matrix for output.
  rownames(GrLsEq)=species.list  # Names rows of output matrix.
  for(i in 1:num.habs)  # Creates names for columns of output matrix.
   {
    if(i==1)
      cols=c(paste("N.Hab.",i,sep=""), paste("Gr.Hab.",i,sep=""), paste("Ls.Hab.",i,sep=""), paste("Eq.Hab.",i,sep=""))
    if(i>1)
      cols=c(cols, paste("N.Hab.",i,sep=""), paste("Gr.Hab.",i,sep=""), paste("Ls.Hab.",i,sep=""), paste("Eq.Hab.",i,sep=""))  
   }
  colnames(GrLsEq)=cols  # Names columns of output matrix.
  for(j in 1:sppno)  # Opens "for loop" through all spp.
   {
    cat("Currently on species = ", species.list[j], "\n")  
       # cat() prints current sp. mnem. to screen.

 # CALCULATIONS FOR OBSERVED RELATIVE DENSITIES ON THE TRUE HABITAT MAP
  
    spmat=matrix(allabund[species.list[j],],plotdimqy,plotdimqx,byrow=F)  # Fills a matrix, with no. rows = plotdimqy and no. columns = plotdimqx, with the indiv. counts per quadrat of one species.
    tot=allabund[sppno+1]  # Creates empty vector for tot. no. stems per quadrat.
    totmat=matrix(tot,plotdimqy,plotdimqx,byrow=F)  # Converts "tot" to a matrix, with no. rows = plotdimqy and no. columns = plotdimqx. 
    spstcnthab=numeric()  # Creates empty vector for stem counts per sp. per habitat.
    totstcnthab=numeric()  # Creates empty vector for tot. stem counts per habitat.
    for(i in 1:num.habs)  
     {
      totstcnthab[i]=sum(totmat[habmat==i])  # Determines tot. no. stems per habitat of the true map.
      spstcnthab[i]=sum(spmat[habmat==i])  # Determines tot. no. stems for focal sp. per habitat of the true map. 	
     }
    spprophab=spstcnthab/totstcnthab  # Calculates observed relative stem density of focal sp. per habitat of the true map.

 # CALCULATIONS FOR RELATIVE DENSITIES ON THE TORUS-BASED HABITAT MAPS
   
    for(x in 0:(plotdimqx-1))  # Opens "for loop" through all 20-m translations along x-axis.
      { 
      for (y in 0:(plotdimqy-1))  # Opens "for loop" through all 20-m translations along y-axis.
       { 
        newhab=matrix(0,plotdimqy,plotdimqx)  # Creates empty matrix of quadrats' habitat designations.
        
        # The following "if" statements create the x,y torus-translation of the habitat map.

        if(y==0 & x==0)          
           newhab=habmat

        if(y==0 & x>0)
           newhab=habmat[c(1:plotdimqy),c((plotdimqx-x+1):plotdimqx,1:(plotdimqx-x))] 
  
        if(x==0 & y>0)
           newhab=habmat[c((plotdimqy-y+1):plotdimqy,1:(plotdimqy-y)),c(1:plotdimqx)]

        if(x>0 & y>0)
            newhab=habmat[c((plotdimqy-y+1):plotdimqy,1:(plotdimqy-y)),c((plotdimqx-x+1):plotdimqx,1:(plotdimqx-x))] 

        Torspstcnthab=numeric()  # Creates empty vector for stem counts per sp. per habitat in torus-based maps.
        Tortotstcnthab=numeric()  # Creates empty vector for tot. stem counts per habitat in torus-based maps.
        for(i in 1:num.habs)
         {
          Tortotstcnthab[i]=sum(totmat[newhab==i])  # Determines tot. no. stems per habitat of the focal torus-based map.
          Torspstcnthab[i]=sum(spmat[newhab==i])  # Determines tot. no. stems for focal sp. per habitat of the focal torus-based map.
         }
        Torspprophab=Torspstcnthab/Tortotstcnthab  # Calculates relative stem density of focal sp. per habitat of the focal torus-based map.  

        for(i in 1:num.habs)
         {

          if(spprophab[i]>Torspprophab[i])  # If rel. dens. of focal sp. in focal habitat of true map is greater than rel. dens. of focal sp. in focal habitat of torus-based map, then add one to "greater than" count. 		
            GrLsEq[j,(4*i)-2]=GrLsEq[j,(4*i)-2]+1  

          if(spprophab[i]<Torspprophab[i])  # If rel. dens. of focal sp. in focal habitat of true map is less than rel. dens. of focal sp. in focal habitat of torus-based map, then add one to "less than" count. 
            GrLsEq[j,(4*i)-1]=GrLsEq[j,(4*i)-1]+1 
          	            
          if(spprophab[i]==Torspprophab[i])  # If rel. dens. of focal sp. in focal habitat of true map is equal to rel. dens. of focal sp. in focal habitat of torus-based map, then add one to "equal to" count.
            GrLsEq[j,4*i]=GrLsEq[j,4*i]+1 
  
         }
       }  # Closes "for loop" through all 20-m translations along y-axis.
     }    # Closes "for loop" through all 20-m translations along x-axis.
    for(i in 1:num.habs)
     {
      GrLsEq[j,(4*i)-3]=spstcnthab[i]
     } 
   }  # Closes "for loop" through all spp.                   
  return(GrLsEq)
 }