*Program FUZZSTEPJA.SAS for fuzzy-set ordination of binary data
 with a step-across routine to minimize "curl-over" distortion;
/*fuzzin is the input data, organized as sites in rows and
   species presence-absence 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 'paula.data';
filename fuzzout 'stepja.data';
filename gradient 'pelev.data';
OPTIONS LINESIZE=72 PAGESIZE=55;
TITLE 'Fuzzy Set Analysis--Binary Data';
title2 'Jaccard"s Coefficient';
DATA A;
  INFILE fuzzin;
  *Replace the 2944 in x2944 with the product of sites x species;
  INPUT x1-x2944;
  *Replace the values for sites and species with the
   appropriate values;
  sites=46; species=64;
  *put input values into data array & also change x3200 to correct
    value;
  ARRAY dta{46,64} x1-x2944;  /*replace 46 with # of sites,*/
  *input gradient scores;     /*64 with # of species,2944 with sitesxspecies*/
  infile gradient;
  *Change 46 in grad46 to # of gradient scores (sites);
  input grad1-grad46;
  array grad{46} grad1-grad46; /*replace 46 with # of sites*/
  *calculate Jaccard similarity indices;
  *possible range is 0 to 1;
  array s{46,46};              /*replace 46 with # of sites*/
  array d{46};                 /*replace 46 with # of sites*/
  array flags{46};             /*replace 46 with # of sites*/
  do x=1 to sites;
    do y = 1 to sites;
      a=0;b=0;c=0;
      do i=1 to species;
       if dta{x,i}=1 and dta{y,i}=1 then a=a+1;
       if dta{y,i}=1 and dta{x,i}=0 then b=b+1;
       if dta{x,i}=1 and dta{y,i}=0 then c=c+1;
      end;
      sdenom = a+b+c;
      if sdenom=0 then sdenom =1;
      s{x,y}=a/sdenom;
    end;
  end;
  *step-across routine--uses Dijkstra's algorithm;
  array dist{46,46};    /*replace 46 with # of sites,*/
  do i = 1 to sites;
    do j = 1 to sites;
      *first calculate distance, the complement of similarity;
      dist{i,j} = 1 - s{i,j};
      if dist{i,j} = 1 then dist {i,j} = 1000;/*sites with no */
    end;                                     /*species in common*/
  end;
  do i = 1 to sites;
    do j = 1 to sites;
      if dist{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}+dist{l,k});
        if d{k}<mini then do;
          mini=d{k};
          flag=k;
        end;
skip: end;
      if flag=j then do;
         dist{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 - dist{i,j};
    end;
  end;
  *relativize s to be between 0 and 1;
  amax=1;amin=0;
  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{46} ma1-ma46;    /*replace 46 with # of sites*/
  array mb{46};             /*replace 46 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{46};    /*replace 46 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{46};    /*replace 46 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{46} me1-me46;    /*replace 46 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-ma46;
proc print;
  var me1-me46;
quit;