*Program FUZZSTEPPS.SAS for fuzzy-set ordination of abundance data
 with a step-across routine to minimize "curl-over" distortion;
/*fuzzin is the input data, organized as sites in rows and
   species abundances in columns
 fuzzout is the output data, containing the gradient scores
   (e.g. elevations) and the fuzzy ordination scores (e.g.
   apparent elevations)
 gradient is the gradient score of each site (e.g. elevation)*/
*;
filename fuzzin 'rba.data';
filename fuzzout 'ascutstepps.data';
filename gradient 'elev.data';
OPTIONS LINESIZE=72 PAGESIZE=55;
TITLE 'Fuzzy Set Analysis--Abundance Data';
title2 'Percent Similarity';
title3 'Step-Across Routine';
*It is often necessary to replace DATA A with DATA _NULL_, otherwise you
 may get a a "too many variables" error statement;
*However, the PROC PRINT statements won't work with DATA _NULL_;
DATA A;
  INFILE fuzzin;
  *Replace the 2000 in x2000 with the product of sites x species;
  INPUT x1-x2000;
  *Replace the values for sites and species with the appropriate values;
  sites=100; species=20;
  *put input values into data array & also change x2000 to correct value;
  ARRAY dta{100,20} x1-x2000;     /*replace 100 with # of sites, */
  *input gradient scores;         /*20 with # of species         */
  infile gradient;
  input grad1-grad100;
  array grad{100} grad1-grad100;  /*replace 100 with # of sites*/
  *relativize grad to be between 0 and 1 if using raw numbers;
  amax=2000;amin=2000; /*chose initial amax and amin to be in middle*/
  do i=1 to sites;     /* of gradient                               */
    if grad{i}<amin then amin=grad{i};
    if grad{i}>amax then amax=grad{i};
  end;
  do i=1 to sites;
    grad{i}=(grad{i}-amin)/(amax-amin);
  end;
  *calculate percent similarity indices;
  *possible range is 0 to 1;
  array s{100,100};    /*replace 100 with # of sites*/
  array d{100};        /*replace 100 with # of sites*/
  array flags{100};    /*replace 100 with # of sites*/
  do x=1 to sites;
    do y=1 to sites;
      t=0;
      do i=1 to species;
        t=t+min(dta{x,i},dta{y,i});
      end;
      s{x,y}=t;
    end;
  end;
  *step-across routine--uses Dijkstra's algorithm;
  array dt{100,100};   /*replace 100 with # of sites*/
  do i = 1 to sites;
    do j = 1 to sites;
      *first calculate distance, the complement of similarity;
      dt{i,j} = 1 - s{i,j};
      if dt{i,j} = 1 then dt {i,j} = 1000; /*sites with no species in common*/
    end;
  end;
 *go to b;
  do i = 1 to sites;
    do j = 1 to sites;
      if dt{i,j} <1000 then go to skip2;
      do x =1 to sites;
        d{x} = 1000;
        flags{x} = 0;
      end;
      l=i;
      d{i} = 0;
      mini=1000;
      flag=0;
      flags{i} =1;
z:    count = 0;
      do k = 1 to sites;
        m=k;
        if flags{k} = 1 then go to skip;
        d{k} = min(d{k}, d{l}+dt{l,k});
        if d{k}<mini then do;
          mini=d{k};
          flag=k;
        end;
skip: end;
      if flag=j then do;
         dt{i,j}=d{j};
         go to skip2;
      end;
      flags{flag}=1;
      l=flag;
      mini=1000;
      go to z;
skip2:end;
  end;
*recalculated similarities from new distances;
  do i = 1 to sites;
    do j = 1 to sites;
      s{i,j} = 1 -dt{i,j};
    end;
  end;
*relativize s to be between 0 and 1;
  amax=0.5;amin=0.5;
  do i=1 to sites;
    do j=1 to sites;
    if s{i,j}<amin then amin=s{i,j};
    if s{i,j}>amax then amax=s{i,j};
    end;
  end;
  do i=1 to sites;
    do j=1 to sites;
      s{i,j}=(s{i,j}-amin)/(amax-amin);
    end;
  end;
  *set up arrays for high (ma) and low (ma) elevation plots;
  array ma{100} ma1-ma100; /*replace 100 with # of sites*/
  array mb{100}; /*replace 100 with # of sites*/
  do i=1 to sites;
    ma{i}=grad{i};
    mb{i}=1-ma{i};
  end;
 *set up array for plots similar to high-elevation plots;
  array mc{100}; /*replace 100 with # of sites*/
  do x=1 to sites;
    sumn=0;sumd=0;
    do y=1 to sites;
      if y=x then sumn=sumn;
        else do;
          sumn=sumn+s{x,y}*ma{y};
        end;
      if y=x then sumd=sumd;
        else sumd=sumd+ma{y};
    end;
    mc{x}=sumn/sumd;
  end;
  *set up array for plots similar to low-elevation plots;
  array md{100}; /*replace 100 with # of sites*/
  do x=1 to sites;
    sumn=0;sumd=0;
    do y=1 to sites;
      if y=x then sumn=sumn;
        else do;
          sumn=sumn+s{x,y}*mb{y};
        end;
      if y=x then sumd=sumd;
        else sumd=sumd+mb{y};
    end;
    md{x}=sumn/sumd;
  end;
  *set up anticommutative difference of c & d--plots similar to high
   elevation while not similar to low elevation plots;
  array me{100} me1-me100;  /*replace 100 with # of sites*/
  do x=1 to sites;
    me{x}=(1+(1-md{x})*(1-md{x})-(1-mc{x})*(1-mc{x}))/2;
  end;
  *relativize me to be between 0 and 1;
  amax=0.50;amin=0.50;
  do i=1 to sites;
    if me{i}<amin then amin= me{i};
    if me{i}>amax then amax= me{i};
  end;
  do i=1 to sites;
    me{i}=(me{i}-amin)/(amax-amin);
  end;
  *put actual and apparent elevations into a file;
  do i=1 to sites;
    FILE fuzzout;
    PUT ma{i} 6.4 me{i} 7.4;
  end;
*print actual gradient scores (e.g. elevations) & fuzzy ordination scores
 (e.g. apparent elevations);
PROC PRINT;
  var ma1-ma100; /*replace 100 with # of sites*/
PROC PRINT;
  var me1-me100; /*replace 100 with # of sites*/
quit;