c This main program reads in the data arrays necessary to find the
C Ivakhnenko polynomial (ip).  Subroutine GMDH finds the IP and stores
c it in the arrays TREE, ITREE, and ITER.  These arrays are then used
c by subsroutine COMP to evaluate the IP at the original data points.
c The user can easily modify this program by either adding a 
c preprocessor subroutine to transform the data before GMDH or by
c adding another read statement to evaluate the IP at other points.
c
c Input data is m = number of variables
c               n = number of data points
c              nt = number of data points in the training set (the first
c                   nt data points are assumed to be in the training set)
c           niter = number of levels GMDH performs before stopping
c                   (if niter = 0 GMDH decides itself)
c           nprnt = 0 if input is not printed
c                 = 1 if input is printed
c              pi = fractional increase in the number of variables
c                   at each iteration (should be between 0 and 1)
c               x = n by m array of independent variables
c               y = n array of dep variables
c
c Order of input: cards 1-10 hollerith info to be printed
c card 11  m, n, nt, pi    format (3i5,f10.2)
c cards 12-  the i'th card contains  y(i), (x(i,j),j=1,m)  i = 1, n
c
	dimension x(100,30), y(100), itree(10,75), tree(10,75,6)
	dimension zz(30), ev(100,30), ysave(100), itr(256), data(80)
c
c       read no. of variables, no. of observations, no. of observations
c       in the training set, rms value
	write (3,4)
	do 60 i = 1, 10
	   read (1,42) data
  42       format (80a1)
	   write (3,45) data
  45       format (10x,80a1)
  60    continue
        read (1,1) m,n,nt,niter,nprnt,pi
   1    format (5i5,5f10.2)
 	do 2 i = 1, n
c          read dep variable, m indep variables
	   read (1,3) y(i), (x(i,j),j=1,m)
   2    ysave(i) = y(i)
   3    format (8f10.2)
   4    format (1h1)
	write (3,43)
  43    format (' --------------------------------------------------------------
     x--------------------------')
	write (3,5) m,n
   5    format (' number of indep var = ',i4/' number of obs = ',i4)
	write (3,6) nt
   6    format (' number of var in training set = ', i4)
	write (3,7) pi
   7    format (' fractional increase in variables = ', f8.3//)
	write (3, 43)
	if (nprnt.eq.0.) goto 32
	write (3, 30)
  30    format (' print input data'//)
	write (3,8)
   8    format (3x,' obs   y     (x(i),i=1,m)'/)
	do 9 i = 1, n
   9       write (3,10) i,y(i),(x(i,j),j=1,m)
  10    format (i6,10f7.2)
        write (3,4)
c       save the data array for later
  32    do 11 i = 1, n
	   do 11 j = 1, m
  11          ev(i,j) = x(i,j)
c       call gmdh subroutine
	call gmdh(n,m,nt,pi,x,y,itree,tree,iter,niter)
	write (3,40)
  40    format (' case no.   observed value      estimate          error
     x            percent error'/)
	do 12 i = 1, n
	   do 13 j = 1, m
  13          zz(j) = ev(i,j)
c          call subroutine to evaluate the Ivakhnenko polynomial
	   call comp(zz,yy,itree,tree,iter,itr)
	   er = ysave(i) - yy
	   perer = 100.*er/ysave(i)
  12    write (3,14) i,ysave(i),yy,er,perer
  14    format (i5,3x,4e18.8)
        write (3,43)
        write (3,4)
c       The Ivakhnenko polynomial is printed only if it is a simple
c       quadratic
        if (iter.gt.1) stop
        write (3,70) (tree(1,1,j), j=1,6)
  70    format(//' Ivakhnenko polynomial '//
     x   ' y = a + b*u + c*v + d*u*u + e*v*v + f*u*v   '/
     x   ' a = ',e12.4/' b = ',e12.4/' c = ',e12.4/
     x   ' d = ',e12.4/' e = ',e12.4/' f = ',e12.4//)
        write (3,71) itr(2),itr(3)
  71    format (' u = x(',i2,' )      v = x(',i2,' )')
        write (3,4)
        stop
        end
c
c-----------------------------------------------------------------------------
	subroutine gmdh(n,m,nt,pi,x,y,itree,tree,iter,niter)
	dimension x(100,30), y(100), tree(10,75,6), itree(10,75)
	dimension poly(6,75), xtx(6,7), xty(6), index(435), work(100,75)
	dimension d(435), ind(435), zzz(6), ma(20), wrk(100)
	dimension xwork(100), ywork(100)
	rms = pi
	ntp1 = nt+ 1
	nc = n - nt
	niter = niter + 1
	mm = m
	iter = 1
	dmin = 1.0e20
  28    l = 1
	mm1 = m-1
c
c       begin loop to compute new variables
c
	do 4 i = 1, mm1
	ip1 = i + 1
	   do 4 j = ip1, m
	      do 5 ii = 1, 6
	   	 xty(ii) = 0.0
	         do 5 jj = 1, 6
   5             xtx(ii,jj) = 0.
              xtx(1,1) = nt
              do 7 k = 1, nt
                 xtx(1,2) = xtx(1,2) + x(k,i)
                 xtx(1,3) = xtx(1,3) + x(k,j)
                 xtx(1,4) = xtx(1,4) + x(k,i)**2
                 xtx(1,5) = xtx(1,5) + x(k,j)**2
                 xtx(1,6) = xtx(1,6) + x(k,i)*x(k,j)
                 xtx(2,2) = xtx(2,2) + x(k,i)**2
                 xtx(2,3) = xtx(2,3) + x(k,i)*x(k,j)
                 xtx(2,4) = xtx(2,4) + x(k,i)**3
                 xtx(2,5) = xtx(2,5) + x(k,i)*x(k,j)**2
                 xtx(2,6) = xtx(2,6) + x(k,i)*x(k,i)*x(k,j)
                 xtx(3,3) = xtx(3,3) + x(k,j)**2
                 xtx(3,4) = xtx(3,4) + x(k,i)*x(k,i)*x(k,j)
                 xtx(3,5) = xtx(3,5) + x(k,j)**3
                 xtx(3,6) = xtx(3,6) + x(k,i)*x(k,j)**2
                 xtx(4,4) = xtx(4,4) + x(k,i)**4
                 xtx(4,5) = xtx(4,5) + (x(k,i)*x(k,j))**2
                 xtx(4,6) = xtx(4,6) + x(k,j)*x(k,i)**3
                 xtx(5,5) = xtx(5,5) + x(k,j)**4
                 xtx(5,6) = xtx(5,6) + x(k,i)*x(k,j)**3
                 xtx(6,6) = xtx(6,6) + (x(k,i)*x(k,j))**2
                 xty(1) = xty(1) + y(k)
                 xty(2) = xty(2) + x(k,i)*y(k)
                 xty(3) = xty(3) + x(k,j)*y(k)
                 xty(4) = xty(4) + x(k,i)*x(k,i)*y(k)
                 xty(5) = xty(5) + x(k,j)*x(k,j)*y(k)
  7              xty(6) = xty(6) + x(k,i)*x(k,j)*y(k)
              do 45 ii = 2, 6
                 im1 = ii-1
                 do 45 jj = 1, im1
  45          xtx(ii,jj) = xtx(jj,ii)
c
c	      solve regression equation for variables i and j
c
	      call sys(xtx,zzz,xty,6,iflag)
              if (iflag.eq.1) goto  41
              do 40 iii = 1, 6
  40          poly(iii,l) = zzz(iii)
              do 15 k = 1, n
                 ww = poly(1,l) + poly(2,l)*x(k,i) + poly(3,l)*x(k,j)
		 ww = ww + poly(4,l)*x(k,i)**2 + poly(5,l)*x(k,j)**2
		 ww = ww + poly(6,l)*x(k,i)*x(k,j)
  15          work(k,l) = ww
              ind(l) = 100*(i+10) + (j+10)
              if (l.eq.75) goto 385
              l = l+1
  41       continue
   4  	continue
c
c       completed construction of m*(m-1)/2 new variables
c
	l = l-1
 385    do 120 i = 1, nc
 120    ywork(i) = y(nt+i)
        do 122 j = 1,l
           do 124 i = 1, nc
 124       xwork(i) = work(nt+i,j)
c
c       compute the goodness of fit statistics
c
	call stat(nc,ywork,xwork,st)
	d(j) = st
 122	index(j) = j
c
c	sort the values of the statistics from low to high
c
	call sort(d,l,index)
	m = m + rms*m
	if (m.gt.l) m = l
  53 	continue
c	the largest number of var is set to 75
	if (m.gt.75) m = 75
	if (m.lt.mm) m = mm
	do 22 j = 1, m
	   itree(iter,j) = ind(index(j))
	   do 23 k = 1, 6
  23       tree(iter,j,k) = poly(k,index(j))
  22    continue
c
c       test for convergence of gmdh algorithm
c
	if (niter.eq.1) goto 55
	if (iter.eq.niter) goto 60
	goto 56
  55    test = d(1) - dmin + .0005
	if (test.ge.0.) goto 60
  56    write (3,300)
 300    format(' ----------------------------------------------------------
     x-------------------------')
	write  (3,200) iter
 200 	format (' Level number = ',i3)
	write  (3,203) m,d(1)
 203 	format (' No. variables saved = ',i3/' rmin value (summed ',
     x	'over checking set) = ',e12.5)
	dmin = d(1)
	ma(iter) = m
	iter = iter + 1
c  	max number of iteration is I1
	do 26 i = 1, n
	   do 26 j = 1, m
  26	x(i,j) = work(i,index(j))
	sum = 0.
	do 190 i = 1, nt
 190	sum = sum + y(i)
	sum = sum/nt
	sum1 = 0.
	sum2 = 0.
	do 191 i = 1, nt
	   sum1 = sum1 + (sum - x(i,1))**2
 191	sum2 = sum2 + (y(i) - sum)**2
	sum = sum1/sum2
	write (3,92) sum
	write (3,300)
	goto 28
  60	if (iter.eq.1) goto 65
	iter = iter - 1
  65	m = mm
	write  (3,130) iter
 130	format ( /' GMDH converged after ',i3,' generation(s)')
	sum = 0.
	do 90 i = 1, nt
  90	sum = sum + y(i)
	sum = sum/nt
	sum1 = 0.
	sum2 = 0.
	do 91 i = 1, nt
	   sum1 = sum1 + (sum-x(i,1))**2
  91	sum2 = sum2 + (y(i) - sum)**2
	sum = sum1/sum2
	write  (3, 92) sum
  92	format (' Multiple correlation (summed over training set) = '
     x	,f6.2)
	write (3,300)
	return
	end
c---------------------------------------------------------------------------
	subroutine sort(a,la,ir)
	integer la,ir(la)
	real a(la)
	integer iu(21),il(21),i,m,j,k,ij,it,l,itt
	real t,tt,r
	if (la.le.0) return
	do 5 i = 1, la
	   if (a(i).lt.0) a(i) = -a(i)
  5	continue
	m = 1
	i = 1
	j = la
	r = .375
  10	if (i.eq.j) goto 55
  15	if (r.gt..5898437) goto 20
	r = r+3.90625e-2
	goto 25
  20	r = r - .21875
  25 	k = i
c
	ij = i + (j-i)*r
	t = a(ij)
	it = ir(ij)
c
	if (a(i).le.t) goto 30
	a(ij) = a(i)
	a(i) = t
	t = a(ij)
	ir(ij) = ir(i)
	ir(i) = it
	it = ir(ij)
  30	l = j
c
	if (a(j).ge.t) goto 40
	a(ij) = a(j)
	a(j) = t
	t = a(ij)
	ir(ij) = ir(j)
	ir(j) = it
	it = ir(ij)
c
	if (a(i).le.t) goto 40
	a(ij) = a(i)
	a(i) = t
	t = a(ij)
	ir(ij) = ir(i)
	ir(i) = it
	it = ir(ij)
	goto 40
 35	if (a(l).eq.a(k)) goto 40
	tt = a(l)
	a(l) = a(k)
	a(k) = tt
	itt = ir(l)
	ir(l) = ir(k)
	ir(k) = itt
c
  40	l = l-1
	if (a(l).gt.t) goto 40
c
  45	k = k+1
	if (a(k).lt.t) goto 45
c
	if (k.le.l) goto 35
c
	if (l-i.le.j-k) goto 50
	 il(m) = i
	iu(m) = l
	i = k
	m = m+1
	goto 60
  50	il(m) = k
	iu(m) = j
	j = l
	m = m + 1
	goto 60
c
  55 	m = m -1
	 if(m.eq.0) return
	i = il(m)
	j = iu(m)
  60	if (j-i.ge.11) goto 25
	if (i.eq.1) goto 10
	i = i-1
  65	i = i+1
	if (i.eq.j) goto 55
	t = a(i+1)
	it = ir(i+1)
	if (a(i).le.t) goto 65
	k = i
  70	a(k+1) = a(k)
	ir(k+1) = ir(k)
	k=k-1
	if (t.lt.a(k)) goto 70
	a(k+1) = t
	ir(k+1) = it
	goto 65
	end
c--------------------------------------------------------------------------
	subroutine sys(a,x,b,n,iflag)
	dimension a(6,7),x(6),index(6),b(6)
	np1 = n + 1
	do 1 i = 1, n
	   index(i) = i
  1	a(i,np1) = b(i)
	do 5 k = 1, n
	   if (k.eq.n) goto 2
	   call inter(a,n,k,index)
  2	   kp1 = k + 1
	   do 4 j = kp1, np1
	      if (abs(a(k,k)).lt..000001) goto 8
	      a(k,j) = a(k,j)/a(k,k)
	      do 4 i = 1, n
	         if (k-i) 3,4,3
  3	         a(i,j) = a(i,j) -a(i,k)*a(k,j)
  4	      continue
  5	continue
	do 7 i = 1,n
	   do 6 j = 1, n
	      if (index(j).ne.i) goto 6
	      x(i) = a(j,np1)
	      goto 7
  6	   continue
  7	continue
	iflag = 0
	return
  8	iflag = 1
	return 
	end
c---------------------------------------------------------------------------
	subroutine inter(a,n,k,index)
	dimension a(6,7),index(6)
	np1 = n + 1
	nr = k
	nc = k
	ab = abs(a(k,k))
	do 2 i = k, n
	   do 2 j = k, n
	      if (abs(a(i,j))-ab) 2,2,1
  1	      ab = abs(a(i,j))
  	      nr = i
	      nc = j
  2	continue
	do 3 j = k, np1
	   de = a(nr,j)
	   a(nr,j) = a(k,j)
  3	a(k,j) = de
	do 4 i = 1, n
	   de = a(i,nc)
	   a(i,nc) = a(i,k)
  4	a(i,k) = de
	is = index(nc)
	index(nc) = index(k)
	index(k) = is
	return
	end
c-------------------------------------------------------------------------
	subroutine comp(x,y,itree,tree,iter,itr)
	dimension x(30),itree(10,75),tree(10,75,6),work(256)
	dimension itr(256)
	it = iter
	itr(1) = 1
	i = 1
  6 	l = 0
	nn = 2**(i-1)
 	n1 = 2**i
	nz = 2**(i+1)-1
	j = n1
  4	jj = itr(nn+l)
	xx = itree(iter,jj)
	itr(j) = itree(iter,jj)/100-10
	iz = itree(iter,jj)/100
	itr(j+1) = xx -100*iz - 10
	j = j + 2
	if (j.gt.nz) goto 3
	l = l + 1
	goto 4
  3	if (iter.eq.1) goto 5
	iter = iter -1
	i = i + 1
	goto 6
  5	iter = it
	nz = 2**(iter-1)
	nzz = nz
	n1 = 2**iter
	do 8 j = 1, nzz
	   jj1 = itr(nz)
	   jj2 = itr(n1)
	   jj3 = itr(n1+1)
	   wk = tree(1,jj1,1) + tree(1,jj1,2)*x(jj2) + 
     x                                           tree(1,jj1,3)*x(jj3)
	   wk = wk + tree(1,jj1,4)*x(jj2)**2 + tree(1,jj1,5)*x(jj3)**2
	   wk = wk + tree(1,jj1,6)*x(jj2)*x(jj3)
	   work(j) = wk
	   nz = nz + 1
  8	n1 = n1 + 2
	iter = iter -1
	if (iter.ne.0) goto 12
	y = work(1)
	iter = it
 	return
  12    i = 2
  18    nz = 2**(iter-1)
  	n1 = 2**iter
	nzz = nz
	n11 = n1
	do 13 j = 1, nzz
	   jj = 2*j-1
	   jj1 = itr(nz)
	   jj2 = itr(n1)
	   jj3 = itr(n1+1)
	   wk = tree(i,jj1,1) + tree(i,jj1,2)*work(jj)
	   wk = wk + tree(i,jj1,3)*work(jj+1) + tree(i,jj1,4)*work(jj)**2
	   wk = wk + tree(i,jj1,5)*work(jj+1)**2
	   wk = wk + tree(i,jj1,6)*work(jj)*work(jj+1)
	   work(n11+j) = wk
	   nz = nz + 1
  13    n1 = n1 + 2
	iter = iter - 1
        if (iter.eq.0) goto 15
	do 14 l = 1, nzz
  14	work(l) = work(n11+l)
	i = i + 1
	goto 18
  15	y = work(3)
	iter = it
	return
        end
c---------------------------------------------------------------------------
	subroutine stat(n,y,x,rms)
	dimension y(50),x(50)
	ss = 0.
	sss = 0.
	do 5 i = 1, n
	   ss = ss + (y(i)-x(i))**2
  5	sss = sss + y(i)**2
	rms = ss/sss
	return
	end