c
c program implementing Friedman-Rafsky two-sample graphics.
c
c Friedman J. H. and Rafsky L. C. (1981). Graphics for the multivariate
c    two-sample problem (with discussion). J. Amer. Statist. Assoc., 76,
c    277-295 (June).
c
c version of 7/17/78.
c
c coded and copyright (c) by:
c
c                           Jerome. H. Friedman
c                         Department of Statistics
c                                  and
c                    Stanford Linear Accelerator Center
c                           Stanford University
c
c
c call tsplot(p,n,data,deg,mst,order,plane)
c
c integer p,n
c
c real data(p,n),deg(n,2),plane(2,n)
c
c integer mst(n,5),order(n,2)
c
c   p = dimension (number of columns of data matrix).
c   n = sample cardinality (number of rows of data matrix) .
c   data = pooled sample data matrix.
c   deg = scratch storage for tsplot. upon return
c         deg(1:n,1) contain the degree of each node (point)
c         in the mst.
c   mst = scratch storage for mstcpr. upon return
c         the integer pairs (i,mst(i,1)), i=1,n-1
c         label the node pairs defining edges of the pooled
c         sample mst.
c
c output quanities:
c
c   order(1:n,1) = rank vector for 'standard' (linear)
c                  multivariate sequence.
c   order(1:n,2) = rank vector for 'radial' multivariate sequence.
c
c   plane(1:2,1:n) = planar locations of observations resulting
c                      from multivariate planing procedure.
c
c   note: the indicies 1:n above correspond to the rows of the
c         data matrix, data(1:p,1:n).
c
c labeled commons:
c
c common /dump/ dodump
c
c   logical dodump
c
c   dodump = print flag for dumping intermediate results during
c            planing. used for debugging purposes.
c          = .true., dumping enabled.
c          = .false., dumping disabled.
c            (default value = .false.)
c
c
      block data
      common /dump/ dodump
      logical dodump
      data dodump /.false./
      end
      subroutine prim (m,n,points,mst,ui,ji,nit)
      integer m,n
      real points(m,n),ui(n)
      integer mst(n),ji(n),nit(n)
      nitp=n-1
      kp=n
      do 10 i=1,nitp
         nit(i)=i
         ui(i)=1.0e73
 10      continue
 20   do 40 i=1,nitp
         ni=nit(i)
         d=0.0
         do 30 i001=1,m
            d=d+(points(i001,ni)-points(i001,kp))**2
 30         continue
         if (ui(i).le.d) go to 40
         ui(i)=d
         ji(i)=kp
 40      continue
      uk=ui(1)
      do 50 i=1,nitp
         if (ui(i).gt.uk) go to 50
         uk=ui(i)
         k=i
 50      continue
      kp=nit(k)
      mst(kp)=ji(k)
      ui(k)=ui(nitp)
      nit(k)=nit(nitp)
      ji(k)=ji(nitp)
      nitp=nitp-1
      if (nitp.ne.0) go to 20
      return
      end
      subroutine tsplot (d,n,x,deg,istor,order,x2)
      integer d,n
      real x(d,n),deg(n,2),x2(2,1)
      integer istor(1),order(n,1)
      integer ttn,ftn,node
      nm1=n-1
      ttn=n+n
      ftn=ttn+ttn
      call prim (d,n,x,istor,deg,istor(ttn+1),istor(ftn+1))
      do 10 j=1,nm1
         deg(j,1)=1.0
 10      continue
      deg(n,1)=0.0
      do 20 j=1,nm1
         deg(istor(j),1)=deg(istor(j),1)+1.0
 20      continue
      istor(ftn+1)=0
      do 30 j=2,n
         istor(ftn+j)=istor(ftn+j-1)+deg(j-1,1)
 30      continue
      do 40 j=1,nm1
         istor(ftn+j)=istor(ftn+j)+1
         istor(ttn+istor(ftn+j))=istor(j)
         node=istor(j)
         istor(ftn+node)=istor(ftn+node)+1
         istor(ttn+istor(ftn+node))=j
 40      continue
      istor(ftn+1)=0
      do 50 j=2,n
         istor(ftn+j)=istor(ftn+j-1)+deg(j-1,1)
 50      continue
      call map (x,d,n,deg,istor(ftn+1),istor(ttn+1),istor,deg(1,2),order
     1,x2)
      call restor (istor,n,deg,istor(ftn+1),istor(ttn+1),deg(1,2))
      return
      end
      subroutine map (x,d,n,deg,pntr,auux,stack,depth,order,x2)
      integer d,n,maxdph,enddia,center,diaend
      real x(d,1),deg(1),depth(1),x2(2,1)
      integer pntr(1),auux(1),stack(2,1),order(n,1)
      maxdph=n+1
      call maxtag (1,deg,pntr,auux,stack,n,depth)
      enddia=diaend(1,deg,pntr,auux,depth,maxdph)
      call maxtag (enddia,deg,pntr,auux,stack,n,depth)
      call line (enddia,deg,pntr,auux,depth,x,d,n,stack,order)
      maxdph=depth(enddia)/2.0+1.0
      center=diaend(enddia,deg,pntr,auux,depth,maxdph)
      call radlin (x,d,n,center,deg,pntr,auux,stack,depth,order(1,2))
      call planit (x,d,n,center,deg,pntr,auux,stack,depth,x2)
      return
      end
      subroutine line (root,deg,pntr,auux,depth,x,d,n,stack,order)
      integer root,d,n,node,stkptr,parent,next,son,it
      real depth(1),deg(1),x(d,1)
      integer pntr(1),auux(1),stack(1),order(1)
      node=root
      parent=0
      stkptr=parent
      deg(node)=-deg(node)
      it=1
      order(1)=node
 10   jdeg=abs(deg(node))
      next=0
      dphmin=9.9e73
      do 20 i=1,jdeg
         son=auux(pntr(node)+i)
         if (deg(son).le.0) go to 20
         if (depth(son).ge.dphmin) go to 20
         dphmin=depth(son)
         next=son
 20      continue
      if (next.ne.0) go to 30
      if (stkptr.eq.0) go to 40
      node=stack(stkptr)
      parent=0
      if (stkptr.gt.1) parent=stack(stkptr-1)
      stkptr=stkptr-1
      go to 10
 30   stkptr=stkptr+1
      parent=node
      stack(stkptr)=parent
      node=next
      it=it+1
      order(it)=node
      deg(node)=-deg(node)
      go to 10
 40   do 50 i=1,n
         deg(i)=abs(deg(i))
 50      continue
      return
      end
      subroutine radlin (x,d,n,start,deg,pntr,auux,stack,depth,order)
      integer d,start,n,node,stkptr,parent,next,istart
      real x(d,n),depth(1),deg(1)
      integer pntr(1),auux(1),stack(n,1),order(1)
      node=start
      parent=0
      istart=1
      stkptr=0
 10   jdeg=abs(deg(node))
      i=istart
      go to 30
 20   i=i+1
 30   if ((i).gt.(jdeg)) go to 40
      next=auux(pntr(node)+i)
      if (next.ne.parent) go to 50
      go to 20
 40   depth(node)=stkptr
      if (stkptr.eq.0) go to 60
      node=stack(stkptr,1)
      parent=0
      if (stkptr.gt.1) parent=stack(stkptr-1,1)
      istart=stack(stkptr,2)
      stkptr=stkptr-1
      go to 10
 50   stkptr=stkptr+1
      stack(stkptr,1)=node
      stack(stkptr,2)=i+1
      parent=node
      node=next
      istart=1
      go to 10
 60   do 70 i=1,n
         order(i)=i
 70      continue
      call sort (depth,order,1,n)
      m=0
 80   m=m+1
      m0=m
 90   if (depth(m).ne.depth(m+1)) go to 100
      m=m+1
      if (m.ge.n) go to 100
      go to 90
 100  if (m.le.m0) go to 130
      do 120 j=m0,m
         kj=order(j)
         depth(j)=0
         do 110 i=1,d
            depth(j)=depth(j)+(x(i,kj)-x(i,start))**2
 110        continue
 120     continue
      call sort (depth,order,m0,m)
 130  if (m.ge.n) go to 140
      go to 80
 140  return
      end
      subroutine planit (x,d,n,start,deg,pntr,auux,stack,depth,x2)
      integer d,start,n,node,stkptr,parent,next,istart,propts,frthst
      integer grandp,other,princp
      real x(d,n),depth(1),deg(1),x2(2,1)
      integer pntr(1),auux(1),stack(n,1)
      common /dump/ dodump
      logical dodump
      data eps /1.0e-70/
      node=start
      parent=0
      istart=1
      stkptr=0
 10   jdeg=abs(deg(node))
      i=istart
      go to 30
 20   i=i+1
 30   if ((i).gt.(jdeg)) go to 40
      next=auux(pntr(node)+i)
      if (next.ne.parent) go to 50
      go to 20
 40   depth(node)=stkptr
      x2(1,node)=parent
      if (stkptr.eq.0) go to 60
      node=stack(stkptr,1)
      parent=0
      if (stkptr.gt.1) parent=stack(stkptr-1,1)
      istart=stack(stkptr,2)
      stkptr=stkptr-1
      go to 10
 50   stkptr=stkptr+1
      stack(stkptr,1)=node
      stack(stkptr,2)=i+1
      parent=node
      node=next
      istart=1
      go to 10
 60   do 70 i=1,n
         stack(i,1)=x2(1,i)
         stack(i,2)=i
 70      continue
      call sort (depth,stack(1,2),1,n)
      m=0
      x2(1,start)=0.0
      x2(2,start)=x2(1,start)
      deg(start)=-deg(start)
 80   m=m+1
      m0=m
 90   if (depth(m).ne.depth(m+1)) go to 100
      m=m+1
      if (m.ge.n) go to 100
      go to 90
 100  if (m.eq.1) go to 80
      if (m.le.m0) go to 130
      do 120 j=m0,m
         node=stack(j,2)
         parent=stack(node,1)
         depth(j)=0.0
         do 110 i01=1,d
            depth(j)=depth(j)+(x(i01,node)-x(i01,parent))**2
 110        continue
         depth(j)=-depth(j)
 120     continue
      call sort (depth,stack(1,2),m0,m)
 130  k=m0
 140  if (deg(stack(k,2)).gt.0.0) go to 150
      k=k+1
      if (k.gt.m) go to 150
      go to 140
 150  if (k.gt.m) go to 330
      propts=stack(k,2)
      parent=stack(propts,1)
      grandp=stack(parent,1)
      if (grandp.ne.0) go to 170
      x2(1,propts)=0.0
      do 160 i01=1,d
         x2(1,propts)=x2(1,propts)+(x(i01,propts)-x(i01,start))**2
 160     continue
      x2(1,propts)=sqrt(x2(1,propts))
      x2(2,propts)=0.0
      deg(propts)=-deg(propts)
      other=propts
      go to 240
 170  dstmax=0.0
      do 190 j=1,m
         if (deg(stack(j,2)).gt.0.0) go to 190
         node=stack(j,2)
         if (node.eq.parent.or.node.eq.grandp) go to 190
         if ((x2(1,node)-x2(1,parent))**2+(x2(2,node)-x2(2,parent))**2.l
     1t.eps) go to 190
         dist=0.0
         do 180 i01=1,d
            dist=dist+(x(i01,propts)-x(i01,node))**2
 180        continue
         if (dist.le.dstmax) go to 190
         dstmax=dist
         frthst=node
 190     continue
      dstpar=0.0
      do 200 i01=1,d
         dstpar=dstpar+(x(i01,propts)-x(i01,parent))**2
 200     continue
      dstpar=sqrt(dstpar)
      dstgpr=0.0
      do 210 i01=1,d
         dstgpr=dstgpr+(x(i01,propts)-x(i01,grandp))**2
 210     continue
      dstgpr=sqrt(dstgpr)
      if (.not.(dodump)) go to 230
      write (6,220) propts,frthst,parent,grandp
 220  format (   30h propts,frthst,parent,grandp = 4i5)
 230  call plot2 (x2(1,frthst),x2(1,parent),x2(1,grandp),sqrt(dstmax),ds
     1tpar,dstgpr,x2(1,propts))
      deg(propts)=-deg(propts)
 240  if ((x2(1,propts)-x2(1,parent))**2+(x2(2,propts)-x2(2,parent))**2.
     1ge.eps) go to 250
      princp=frthst
      go to 260
 250  princp=propts
 260  jdeg=abs(deg(parent))
      do 320 i=1,jdeg
         node=auux(pntr(parent)+i)
         if (node.eq.grandp.or.node.eq.propts) go to 320
         dstprn=0.0
         do 270 i01=1,d
            dstprn=dstprn+(x(i01,node)-x(i01,princp))**2
 270        continue
         dstprn=sqrt(dstprn)
         dstpar=0.0
         do 280 i01=1,d
            dstpar=dstpar+(x(i01,node)-x(i01,parent))**2
 280        continue
         dstpar=sqrt(dstpar)
         if (grandp.ne.0) other=grandp
         dstoth=0.0
         do 290 i01=1,d
            dstoth=dstoth+(x(i01,node)-x(i01,other))**2
 290        continue
         dstoth=sqrt(dstoth)
         if (.not.(dodump)) go to 310
         write (6,300) node,princp,parent,other
 300  format (   27h node,princp,parent,other = 4i5)
 310     call plot2 (x2(1,princp),x2(1,parent),x2(1,other),dstprn,dstpar
     1,dstoth,x2(1,node))
         deg(node)=-deg(node)
         other=node
 320     continue
      m0=k+1
      if (m0.gt.m) go to 330
      go to 130
 330  if (m.ge.n) go to 340
      go to 80
 340  return
      end
      subroutine plot2 (x1,x2,xc,r1,r2,rc,xn)
      real xa(2),xb(2),xc(2),xna(2),xnb(2),xn(2),ra,rb,rc
      real dsq,d,d1,s,tmp,x1(2),x2(2)
      common /dump/ dodump
      logical dodump
      if (r1.ge.r2) go to 10
      ra=r1
      rb=r2
      xa(1)=x1(1)
      xa(2)=x1(2)
      xb(1)=x2(1)
      xb(2)=x2(2)
      go to 20
 10   ra=r2
      rb=r1
      xa(1)=x2(1)
      xa(2)=x2(2)
      xb(1)=x1(1)
      xb(2)=x1(2)
 20   dsq=(xa(1)-xb(1))**2+(xa(2)-xb(2))**2
      d=sqrt(dsq)
      if ((d.lt.ra+rb).and.(d.gt.rb-ra)) go to 40
      d1=ra/d
      if (d.le.rb-ra) d1=-d1
      xn(1)=(xb(1)-xa(1))*d1+xa(1)
      xn(2)=(xb(2)-xa(2))*d1+xa(2)
      if (.not.(dodump)) go to 80
      write (6,30)
 30   format (   13h approx soln.)
      go to 80
 40   d1=(dsq+ra**2-rb**2)*0.5/d
      s=sqrt(ra**2-d1**2)/d
      d1=d1/d
      xna(1)=(xb(2)-xa(2))*s
      xnb(1)=-xna(1)
      xna(2)=(xa(1)-xb(1))*s
      xnb(2)=-xna(2)
      tmp=xa(1)+(xb(1)-xa(1))*d1
      xna(1)=xna(1)+tmp
      xnb(1)=xnb(1)+tmp
      tmp=xa(2)+(xb(2)-xa(2))*d1
      xna(2)=xna(2)+tmp
      xnb(2)=xnb(2)+tmp
      if (abs(sqrt((xc(1)-xna(1))**2+(xc(2)-xna(2))**2)-rc).gt.abs(sqrt(
     1(xc(1)-xnb(1))**2+(xc(2)-xnb(2))**2)-rc)) go to 50
      xn(1)=xna(1)
      xn(2)=xna(2)
      go to 60
 50   xn(1)=xnb(1)
      xn(2)=xnb(2)
 60   if (.not.(dodump)) go to 80
      write (6,70)
 70   format (   12h exact soln.)
 80   if (.not.(dodump)) go to 110
      dst1=sqrt((x1(1)-xn(1))**2+(x1(2)-xn(2))**2)
      dst2=sqrt((x2(1)-xn(1))**2+(x2(2)-xn(2))**2)
      write (6,90) r1,dst1,r2,dst2
 90   format (   11h 1st dist = 2g13.5,    11h 2nd dist =2g13.5)
      write (6,100) x1(1),x1(2),x2(1),x2(2),xn(1),xn(2)
 100  format (   15h plane points =6g13.5)
 110  return
      end
      subroutine maxtag (start,deg,pntr,auux,stack,n,count)
      integer start,n,node,stkptr,parent,next,istart
      real count(1),deg(1),maxcnt
      integer pntr(1),auux(1),stack(2,1)
      do 10 i=1,n
         count(i)=0.0
 10      continue
      node=start
      parent=0
      istart=1
      stkptr=0
 20   jdeg=abs(deg(node))
      i=istart
      go to 40
 30   i=i+1
 40   if ((i).gt.(jdeg)) go to 50
      next=auux(pntr(node)+i)
      if (next.ne.parent) go to 90
      go to 30
 50   maxcnt=0.0
      i=1
      go to 70
 60   i=i+1
 70   if ((i).gt.(jdeg)) go to 80
      next=auux(pntr(node)+i)
      if (next.ne.parent) maxcnt=amax1(maxcnt,count(next))
      go to 60
 80   count(node)=maxcnt+1.0
      if (stkptr.eq.0) return
      node=stack(1,stkptr)
      parent=0
      if (stkptr.gt.1) parent=stack(1,stkptr-1)
      istart=stack(2,stkptr)
      stkptr=stkptr-1
      go to 20
 90   stkptr=stkptr+1
      stack(1,stkptr)=node
      stack(2,stkptr)=i+1
      parent=node
      node=next
      istart=1
      go to 20
      end
      integer function diaend(start,deg,pntr,auux,count,maxdph)
      integer start,node,parent,next,nexts,depth,maxdph
      real count(1),deg(1),maxcnt
      integer pntr(1),auux(1)
      node=start
      parent=0
      depth=0
 10   maxcnt=0.0
      depth=depth+1
      jdeg=abs(deg(node))
      i=1
      go to 30
 20   i=i+1
 30   if ((i).gt.(jdeg)) go to 40
      next=auux(pntr(node)+i)
      if (next.eq.parent) go to 20
      if (count(next).le.maxcnt) go to 20
      maxcnt=count(next)
      nexts=next
      go to 20
 40   if ((maxcnt.ne.0.0).and.(depth.lt.maxdph)) go to 50
      diaend=node
      return
 50   parent=node
      node=nexts
      go to 10
      end
      subroutine restor (istor,start,deg,pntr,auux,stack)
      integer start,node,stkptr,parent,next,istart
      real deg(1)
      integer istor(1),pntr(1),auux(1),stack(2,1)
      node=start
      parent=0
      istart=1
      stkptr=0
 10   jdeg=abs(deg(node))
      i=istart
      go to 30
 20   i=i+1
 30   if ((i).gt.(jdeg)) go to 40
      next=auux(pntr(node)+i)
      if (next.ne.parent) go to 50
      go to 20
 40   istor(node)=parent
      if (stkptr.eq.0) return
      node=stack(1,stkptr)
      parent=0
      if (stkptr.gt.1) parent=stack(1,stkptr-1)
      istart=stack(2,stkptr)
      stkptr=stkptr-1
      go to 10
 50   stkptr=stkptr+1
      stack(1,stkptr)=node
      stack(2,stkptr)=i+1
      parent=node
      node=next
      istart=1
      go to 10
      end
      subroutine sort (v,a,ii,jj)
c
c     puts into a the permutation vector which sorts v into
c     increasing order.  only elements from ii to jj are considered.
c     arrays iu(k) and il(k) permit sorting up to 2**(k+1)-1 elements
c
c     this is a modification of cacm algorithm #347 by r. c. singleton,
c     which is a modified hoare quicksort.
c
      dimension a(jj),v(1),iu(20),il(20)
      integer t,tt
      integer a
      real v
      m=1
      i=ii
      j=jj
 10   if (i.ge.j) go to 80
 20   k=i
      ij=(j+i)/2
      t=a(ij)
      vt=v(ij)
      if (v(i).le.vt) go to 30
      a(ij)=a(i)
      a(i)=t
      t=a(ij)
      v(ij)=v(i)
      v(i)=vt
      vt=v(ij)
 30   l=j
      if (v(j).ge.vt) go to 50
      a(ij)=a(j)
      a(j)=t
      t=a(ij)
      v(ij)=v(j)
      v(j)=vt
      vt=v(ij)
      if (v(i).le.vt) go to 50
      a(ij)=a(i)
      a(i)=t
      t=a(ij)
      v(ij)=v(i)
      v(i)=vt
      vt=v(ij)
      go to 50
 40   a(l)=a(k)
      a(k)=tt
      v(l)=v(k)
      v(k)=vtt
 50   l=l-1
      if (v(l).gt.vt) go to 50
      tt=a(l)
      vtt=v(l)
 60   k=k+1
      if (v(k).lt.vt) go to 60
      if (k.le.l) go to 40
      if (l-i.le.j-k) go to 70
      il(m)=i
      iu(m)=l
      i=k
      m=m+1
      go to 90
 70   il(m)=k
      iu(m)=j
      j=l
      m=m+1
      go to 90
 80   m=m-1
      if (m.eq.0) return
      i=il(m)
      j=iu(m)
 90   if (j-i.ge.11) go to 20
      if (i.eq.ii) go to 10
      i=i-1
 100  i=i+1
      if (i.eq.j) go to 80
      t=a(i+1)
      vt=v(i+1)
      if (v(i).le.vt) go to 100
      k=i
 110  a(k+1)=a(k)
      v(k+1)=v(k)
      k=k-1
      if (vt.lt.v(k)) go to 110
      a(k+1)=t
      v(k+1)=vt
      go to 100
      end