# This is a shell archive.  Remove anything before this line,
# then unpack it by saving it in a file and typing "sh file".
#
# Wrapped by ssyx.ucsc.edu!mmcohen on Tue Nov 21 00:05:04 PST 1989
# Contents:  README fx lram1.f lram1.fit lram1.lst lrd r s.f stepit.doc steps.f
 
echo x - README
sed 's/^@//' > "README" <<'@//E*O*F README//'
1) compile steps.f into steps.o and s.f into s.o
   and put in lib directory (e.g. /usr/lib/local). see file fx
   (which you will want to modify for compile flags, which may be diff
   for you, depending on machine. this was for a sparc2)
2) put fx and r file in your bin directory or
   wherever else you put your executables.
3) s.f contains stepit routine
   steps.f contains cvarious additional routines to
         make fitting easy
   lram1.f is typical program which fits flmp to some data
         - massaro & cohen 1983 exp2  phonlogical constraints exp
   lram1.fit is output file (should rename this and check that you get
           same results) has obs and predicted data
   lram1.lst is listing of rmsds, parameters, and correls
   stepit.doc is short document on how stepit works.
4) to run lram1, say "r lram1". This should compile lram1.f,
   link it with s.o and steps.o, rename a.out to vlram1, and run it.
   can also say "r lram1&". on vax, might want to add to r file
   to delete lram1.o which is otherwise leftover.
5) please write me at mmcohen@ssyx.ucsc.edu if you have any problems
   or questions.
@//E*O*F README//
chmod u=rw,g=r,o=r README
 
echo x - fx
sed 's/^@//' > "fx" <<'@//E*O*F fx//'
f77 $1 /usr/lib/local/s.o /usr/lib/local/steps.o
@//E*O*F fx//
chmod u=rwx,g=,o= fx
 
echo x - lram1.f
sed 's/^@//' > "lram1.f" <<'@//E*O*F lram1.f//'
c	lram1.f
c
	COMMON /FUN/ pre(7,4),obs(7,4)
	real d(28,7)
	external funk,pp
	call sas('lrd','lram1.lst','lram1.fit')
	write(7,"(' *** fuzzy model 1 of lr data ***')")
	read(9,"(7f8.5)")d
	call runstp(funk,11,28,7,d,obs,pre,pp)
	stop
	end
	subroutine funk
	common /FUN/  pre(7,4),obs(7,4)
	common /STEP/ nv,ntrace,matrix,chisq,mask(100),x(100)
	real g(7),c(4)
	equivalence (x(1),g),(x(8),c)
	do 500 j=1,7
	do 500 k=1,4
	t=g(j)*c(k)
	tn=(1.-g(j))*(1.-c(k))
500	pre(j,k)=t/(t+tn)
	call cchi(obs,pre,28)
	RETURN
	entry pp(is,r)
	write(7,117)is,r,g,c
117	format(/' Sub  'i2' rmsd= 'f7.4/' g  = '
&	7f7.4/' c  = '4f7.4)
	return
	end
@//E*O*F lram1.f//
chmod u=rw,g=r,o=r lram1.f
 
echo x - lram1.fit
sed 's/^@//' > "lram1.fit" <<'@//E*O*F lram1.fit//'
    0.02500    0.05000    0.00000    0.20000    0.47500    0.97500
    0.97500    0.20000    0.40000    0.55000    0.70000    0.87500
    0.97500    0.97500    0.00000    0.02500    0.00000    0.00000
    0.17500    0.95000    1.00000    0.00000    0.00000    0.10000
    0.30000    0.70000    1.00000    1.00000    0.01455    0.03608
    0.06255    0.15187    0.48699    0.98838    0.99887    0.19486
    0.38019    0.52233    0.74583    0.93960    0.99928    0.99993
    0.00310    0.00782    0.01386    0.03634    0.16661    0.94713
    0.99467    0.03129    0.07567    0.12736    0.28143    0.67494
    0.99465    0.99948
    0.47500    0.57500    0.62500    0.72500    0.95000    0.97500
    0.97500    0.27500    0.52500    0.77500    0.70000    0.92500
    1.00000    1.00000    0.00000    0.05000    0.00000    0.05000
    0.32500    0.85000    0.90000    0.05000    0.17500    0.37500
    0.57500    0.82500    0.97500    0.95000    0.37648    0.55146
    0.71244    0.78938    0.95429    0.99616    0.99748    0.33706
    0.50867    0.67599    0.75938    0.94618    0.99544    0.99701
    0.01287    0.02586    0.05078    0.07487    0.31075    0.84839
    0.89523    0.13025    0.23369    0.38062    0.48176    0.83815
    0.98468    0.98991
    0.07500    0.47500    0.87500    0.97500    1.00000    1.00000
    1.00000    0.62500    0.97500    0.97500    1.00000    1.00000
    0.97500    1.00000    0.00000    0.00000    0.00000    0.30000
    0.67500    0.95000    1.00000    0.60000    0.67500    0.75000
    0.90000    0.75000    0.85000    0.80000    0.18806    0.44845
    0.80572    0.98397    0.99668    0.99963    0.99996    0.67229
    0.87808    0.97350    0.99816    0.99962    0.99996    1.00000
    0.00160    0.00558    0.02781    0.29748    0.67411    0.94964
    0.99492    0.48927    0.77080    0.94491    0.99608    0.99919
    0.99991    0.99999
    0.07500    0.00000    0.00000    0.27500    0.75000    0.90000
    0.82500    0.00000    0.05000    0.32500    0.70000    0.90000
    0.97500    1.00000    0.00000    0.00000    0.02500    0.02500
    0.17500    0.52500    0.72500    0.02500    0.00000    0.07500
    0.30000    0.92500    0.95000    0.92500    0.00167    0.00584
    0.05080    0.22142    0.75150    0.91934    0.96031    0.01378
    0.04686    0.30941    0.70423    0.96200    0.98963    0.99509
    0.00016    0.00055    0.00502    0.02612    0.22193    0.51807
    0.69531    0.00298    0.01040    0.08739    0.33725    0.84402
    0.95326    0.97742
    0.00000    0.07500    0.32500    0.70000    0.97500    1.00000
    0.97500    0.00000    0.00000    0.00000    0.22500    0.70000
    0.95000    1.00000    0.00000    0.00000    0.00000    0.02500
    0.72500    0.97500    1.00000    0.02500    0.05000    0.15000
    0.42500    0.87500    0.95000    1.00000    0.00718    0.08196
    0.32489    0.69774    0.97221    0.99716    0.99994    0.00058
    0.00710    0.03713    0.15609    0.73704    0.96564    0.99929
    0.00043    0.00534    0.02813    0.12192    0.67784    0.95475
    0.99906    0.00228    0.02742    0.13191    0.42160    0.91698
    0.99105    0.99982
    0.00000    0.20000    0.37500    0.85000    1.00000    1.00000
    1.00000    0.37500    0.47500    0.65000    0.57500    0.70000
    0.72500    0.80000    0.00000    0.02500    0.05000    0.12500
    0.72500    1.00000    1.00000    0.00000    0.07500    0.32500
    0.85000    0.95000    1.00000    0.97500    0.11004    0.19195
    0.37327    0.79570    0.98285    0.99984    0.99983    0.28688
    0.43595    0.65960    0.92686    0.99467    0.99995    0.99995
    0.00540    0.01033    0.02550    0.14610    0.71574    0.99628
    0.99604    0.08892    0.15790    0.31978    0.75456    0.97837
    0.99979    0.99978
    0.02500    0.10000    0.52500    0.90000    0.97500    0.92500
    0.90000    0.90000    1.00000    1.00000    1.00000    1.00000
    1.00000    1.00000    0.00000    0.00000    0.00000    0.02500
    0.27500    1.00000    1.00000    0.00000    0.00000    0.10000
    0.55000    1.00000    1.00000    1.00000    0.01547    0.10077
    0.52089    0.91332    0.99800    0.99999    0.99999    0.90160
    0.98492    0.99842    0.99984    1.00000    1.00000    1.00000
    0.00001    0.00009    0.00083    0.00799    0.27623    0.99243
    0.99243    0.00179    0.01265    0.11053    0.54633    0.98277
    0.99995    0.99995
    0.09643    0.21071    0.38929    0.66071    0.87500    0.96786
    0.95000    0.33929    0.48929    0.61071    0.70000    0.87143
    0.94286    0.96786    0.00000    0.01429    0.01071    0.07857
    0.43929    0.89286    0.94643    0.10000    0.13929    0.26786
    0.55714    0.86071    0.96071    0.95000    0.14787    0.24017
    0.37388    0.61526    0.89662    0.98926    0.99427    0.30567
    0.44503    0.60238    0.80226    0.95653    0.99574    0.99773
    0.01427    0.02570    0.04746    0.11773    0.41987    0.88483
    0.93535    0.10478    0.17572    0.28712    0.51890    0.85401
    0.98415    0.99152
@//E*O*F lram1.fit//
chmod u=rw,g=r,o=r lram1.fit
 
echo x - lram1.lst
sed 's/^@//' > "lram1.lst" <<'@//E*O*F lram1.lst//'
 *** fuzzy model 1 of lr data ***

 Sub   1 rmsd=  0.0307
 g  =  0.0177 0.0437 0.0753 0.1794 0.5369 0.9905 0.9991
 c  =  0.4502 0.9307 0.1471 0.6417
 obs~ pre: r= 0.997 r2= 0.995 t=    70.64 df=   26

 Sub   2 rmsd=  0.0484
 g  =  0.1430 0.2536 0.4064 0.5088 0.8523 0.9862 0.9909
 c  =  0.7835 0.7529 0.0725 0.4730
 obs~ pre: r= 0.991 r2= 0.981 t=    37.11 df=   26

 Sub   3 rmsd=  0.0896
 g  =  0.0754 0.2227 0.5937 0.9558 0.9906 0.9990 0.9999
 c  =  0.7395 0.9617 0.0192 0.9215
 obs~ pre: r= 0.970 r2= 0.940 t=    20.24 df=   26

 Sub   4 rmsd=  0.0419
 g  =  0.0031 0.0107 0.0896 0.3434 0.8476 0.9545 0.9780
 c  =  0.3522 0.8199 0.0488 0.4932
 obs~ pre: r= 0.995 r2= 0.990 t=    50.08 df=   26

 Sub   5 rmsd=  0.0306
 g  =  0.0016 0.0200 0.0991 0.3453 0.8888 0.9877 0.9998
 c  =  0.8140 0.2596 0.2084 0.5802
 obs~ pre: r= 0.998 r2= 0.995 t=    72.32 df=   26

 Sub   6 rmsd=  0.1162
 g  =  0.0606 0.1103 0.2372 0.6704 0.9677 0.9997 0.9997
 c  =  0.6570 0.8617 0.0776 0.6019
 obs~ pre: r= 0.960 r2= 0.922 t=    17.53 df=   26

 Sub   7 rmsd=  0.0251
 g  =  0.0009 0.0065 0.0596 0.3806 0.9668 0.9999 0.9999
 c  =  0.9449 0.9999 0.0129 0.6622
 obs~ pre: r= 0.998 r2= 0.997 t=    92.89 df=   26

 Sub   8 rmsd=  0.0390
 g  =  0.0957 0.1617 0.2670 0.4939 0.8411 0.9825 0.9906
 c  =  0.6211 0.8061 0.1203 0.5250
 obs~ pre: r= 0.995 r2= 0.990 t=    51.52 df=   26

 Sub  ** rmsd=  0.0546
 g  =  0.0432 0.0954 0.2230 0.4834 0.8644 0.9882 0.9953
 c  =  0.6773 0.7981 0.0838 0.6248
@//E*O*F lram1.lst//
chmod u=rw,g=r,o=r lram1.lst
 
echo x - lrd
sed 's/^@//' > "lrd" <<'@//E*O*F lrd//'
 0.02500 0.05000 0.      0.20000 0.47500 0.97500 0.97500
 0.20000 0.40000 0.55000 0.70000 0.87500 0.97500 0.97500
 0.      0.02500 0.      0.      0.17500 0.95000 1.00000
 0.      0.      0.10000 0.30000 0.70000 1.00000 1.00000
 0.47500 0.57500 0.62500 0.72500 0.95000 0.97500 0.97500
 0.27500 0.52500 0.77500 0.70000 0.92500 1.00000 1.00000
 0.      0.05000 0.      0.05000 0.32500 0.85000 0.90000
 0.05000 0.17500 0.37500 0.57500 0.82500 0.97500 0.95000
 0.07500 0.47500 0.87500 0.97500 1.00000 1.00000 1.00000
 0.62500 0.97500 0.97500 1.00000 1.00000 0.97500 1.00000
 0.      0.      0.      0.30000 0.67500 0.95000 1.00000
 0.60000 0.67500 0.75000 0.90000 0.75000 0.85000 0.80000
 0.07500 0.      0.      0.27500 0.75000 0.90000 0.82500
 0.      0.05000 0.32500 0.70000 0.90000 0.97500 1.00000
 0.      0.      0.02500 0.02500 0.17500 0.52500 0.72500
 0.02500 0.      0.07500 0.30000 0.92500 0.95000 0.92500
 0.      0.07500 0.32500 0.70000 0.97500 1.00000 0.97500
 0.      0.      0.      0.22500 0.70000 0.95000 1.00000
 0.      0.      0.      0.02500 0.72500 0.97500 1.00000
 0.02500 0.05000 0.15000 0.42500 0.87500 0.95000 1.00000
 0.      0.20000 0.37500 0.85000 1.00000 1.00000 1.00000
 0.37500 0.47500 0.65000 0.57500 0.70000 0.72500 0.80000
 0.      0.02500 0.05000 0.12500 0.72500 1.00000 1.00000
 0.      0.07500 0.32500 0.85000 0.95000 1.00000 0.97500
 0.02500 0.10000 0.52500 0.90000 0.97500 0.92500 0.90000
 0.90000 1.00000 1.00000 1.00000 1.00000 1.00000 1.00000
 0.      0.      0.      0.02500 0.27500 1.00000 1.00000
 0.      0.      0.10000 0.55000 1.00000 1.00000 1.00000
@//E*O*F lrd//
chmod u=rw,g=r,o=r lrd
 
echo x - r
sed 's/^@//' > "r" <<'@//E*O*F r//'
fx $1.f; mv a.out v$1; v$1; rm v$1
@//E*O*F r//
chmod u=rwx,g=,o= r
 
echo x - s.f
sed 's/^@//' > "s.f" <<'@//E*O*F s.f//'

C  STEPIT 5.15 ...  STEPIT WITH /JVARY/.  SEPTEMBER 2, 1966.
      SUBROUTINE STEPIT(FUNK)
C  AVAILABLE FROM....        QUANTUM CHEMISTRY PROGRAM EXCHANGE
C                            I.U. CHEMISTRY DEPT., BLOOMINGTON, INDIANA.
      DIMENSION VEC(100),TRIAL(100),XSAVE(100),CHI(100),DX(100)
      DIMENSION OLDVEC(100),SALVO(100),XOSC(100,15),CHIOSC(15)
      DIMENSION VECFLG(100)
      COMMON /STEP/ NV,NTRACE,MATRIX,CHISQ,MASK(100),X(100),XMAX(100),
     *    XMIN(100),DELTAX(100),DELMIN(100),ERR(100,15)
      COMMON/FASTER/ JVARY
      COMMON/TEST/ NFB,CRIT
	EXTERNAL FUNK
	NVMAX = 100
      MOSQUE=15
      RATIO = 10.0
      COLIN=0.99
      NCOMP=5
      ACK=2.0
      SIGNIF=2.E8
      HUGE = 1.E35
      IF(NFB.EQ.0)NFB=5000
C
   40 IF(NV)290,290,50
   50 NACTIV=0
      DO 150 I=1,NV
      IF(MASK(I))150,60,150
   60 IF(SIGNIF*ABS (DELTAX(I))-ABS (X(I)))70,70,100
   70 IF(X(I))90,80,90
   80 DELTAX(I)=0.01
      GO TO 100
   90 DELTAX(I)=0.01*X(I)
  100 IF(DELMIN(I))120,110,120
  110 DELMIN(I)=DELTAX(I)/SIGNIF
  120 IF(XMAX(I)-XMIN(I))130,130,140
  130 XMAX(I)=HUGE
      XMIN(I)=-HUGE
  140 NACTIV=NACTIV+1
      X(I)=AMAX1(XMIN(I),AMIN1(XMAX(I),X(I)))
  150 CONTINUE
      COMPAR=0.0
      IF(NACTIV-1)160,190,180
  160 DO 170 J=1,NV
  170 MASK(J)=0
      GO TO 50
  180 A=NACTIV
      SUB=2.0/(A-1.0)
      P=2.0*(1.0/SQRT (A)/(1.0-0.5**SUB)-1.0)
      COMPAR=AMIN1(.999,ABS((1.0-(1.0-COLIN)**SUB)*(1.0+P*(1.-COLIN))))
  190 IF(NTRACE)280,200,200
  200 WRITE(6,210)
  210 FORMAT('\f'
&' ENTER SUBROUTINE STEPIT V5.15 - COPYRIGHT 1965 JP CHANDLER'
& /,' PHYSICS DEPT., INDIANA UNIVERSITY.',/,' INITIAL VALUES...'/)
      WRITE(6,220) (MASK(J),J = 1,NV)
      WRITE(6,230) (X(J),J=1,NV)
      WRITE(6,240) (XMAX(J),J=1,NV)
      WRITE(6,250) (XMIN(J) ,J=1,NV)
      WRITE(6,260) (DELTAX(J),J=1,NV)
      WRITE(6,270) (DELMIN(J),J=1,NV)
220   FORMAT(/,' MASK   = ', 5(I6,6X)/(4X,5I12))
  230 FORMAT(/,' X      = ', 5E12.4,/,' ', 9X,5E12.4)
  240 FORMAT(/,' XMAX   = ', 5E12.4,/,' ', 9X,5E12.4)
  250 FORMAT(/,' XMIN   = ', 5E12.4,/,' ', 9X,5E12.4)
  260 FORMAT(/,' DELTAX = ', 5E12.4,/,' ', 9X,5E12.4)
  270 FORMAT(/,' DELMIN = ', 5E12.4,/,' ', 9X,5E12.4)
  280 JVARY=0
      CALL FUNK
      NF=1
      JOCK=1
  290 IF(NTRACE)320,300,300
  300 WRITE(6,310)NV,NACTIV,MATRIX,NCOMP,RATIO,ACK,COLIN,COMPAR,CHISQ
  310 FORMAT(//' 'I3,' VARIABLES',I3,' ACTIVE.'10X'MATRIX ='I4,10X
&	'NCOMP ='
&	I2//' RATIO ='F5.1,6X'ACK ='F5.1,6X'COLIN ='F6.3,6X'COMPAR ='
&	F6.3///' CHISQ ='E15.8)
320   IF(NV.GT.0)GO TO 330
      WRITE(6,321)
321   FORMAT(' NV=0 !!!!!!!!!!')
      RETURN
  330 IF(NTRACE)360,360,340
  340 WRITE(6,350)
  350 FORMAT(//38(' *')//10X,' TRACE MAP OF THE MINIMIZATION PROCESS'//)
C
  360 DO 370 I=1,NV
      DX(I)=DELTAX(I)
      VEC(I)=0.
      VECFLG(I) = 0.0
C     DO 370 J=1,NV
      DO 370 J=1,MOSQUE
  370 ERR(I,J)=0.
      CHIOLD=CHISQ
      NOSC=0
C
  380 NCIRC=0
      NZIP=0
C
C  MAIN DO LOOP FOR CYCLING THROUGH THE VARIABLES.
C  FIRST TRIAL STEP WITH EACH VARIABLE IS SEPARATE.
C
  390 NACK=0
      DO 1350 I=1,NV
      OLDVEC(I)=VEC(I)
      VEC(I)=0.0
      VECFLG(I) = 0.0
      TRIAL(I)=0.0
      IF(MASK(I))400,410,400
  400 VEC(I)= 0.0
C 400 VEC(I)=-0.0
      VECFLG(I) = 1.0
      GO TO 1350
  410 NACK=NACK+1
      SAVE=X(I)
      IF(SIGNIF*ABS (DX(I))-ABS (X(I)))580,580,420
  420 X(I)=SAVE+DX(I)
      JVARY=0
      IF(JOCK)440,440,430
  430 JOCK=0
      JVARY=I
  440 NFLAG=1
C                    CHECK LIMITS
      IF(X(I)-XMIN(I))460,450,450
  450 IF(X(I)-XMAX(I))470,470,460
  460 NFLAG=NFLAG+3
      GO TO 490
  470 CALL FUNK
      NF=NF+1
      JVARY=I
      CHIME=CHISQ
      IF(CHISQ-CHIOLD)620,480,490
C                    BETTER, NO DIFF, WORSE
  480 NFLAG=NFLAG+1
  490 X(I)=SAVE-DX(I)
      IF(X(I)-XMIN(I))590,500,500
C                    CHECK LIMITS...
  500 IF(X(I)-XMAX(I))510,510,590
  510 CALL FUNK
      NF=NF+1
      JVARY=I
      IF(CHISQ-CHIOLD)610,520,530
  520 NFLAG=NFLAG+1
  530 IF(NFLAG-3)540,580,590
  540 IF((CHISQ-CHIME)*(CHIME-2.*CHIOLD+CHISQ))550,590,550
  550 TRIAL(I)=.5*DX(I)*(CHISQ-CHIME)/(CHIME-2.*CHIOLD+CHISQ)
      VEC(I)=TRIAL(I)/ABS (DX(I))
      VECFLG(I) = 0.0
      X(I)=SAVE+TRIAL(I)
      CALL FUNK
      NF=NF+1
      IF(CHISQ-CHIOLD)560,570,570
  560 CHIOLD=CHISQ
      JOCK=1
      GO TO 600
  570 TRIAL(I)=0.0
      VECFLG(I) = 0.0
      VEC(I)=0.0
      GO TO 590
  580 VEC(I)= 0.0
C 580 VEC(I)=-0.0
      VECFLG(I) = 1.0
C                    RESTORE...
  590 X(I)=SAVE
  600 NCIRC=NCIRC+1
      IF(NCIRC-NACTIV)690,1430,1430
  610 DX(I)=-DX(I)
C
C          A LOWER VALUE HAS BEEN FOUND. HENCE THIS VARIABLE WILL CHANGE
C
  620 NCIRC=0
      DEL=DX(I)
  630 CHIME=CHIOLD
      CHIOLD=CHISQ
      VEC(I)=VEC(I)+DEL/ABS (DX(I))
      TRIAL(I)=TRIAL(I)+DEL
      VECFLG(I) = 0.0
      DEL=ACK*DEL
      SAVE=X(I)
      X(I)=SAVE+DEL
      IF(X(I)-XMIN(I))680,640,640
  640 IF(X(I)-XMAX(I))650,650,680
  650 CALL FUNK
      NF=NF+1
      IF(CHISQ-CHIOLD)630,655,655
C                    PERFORM PARABOLIC INTERPOLATION...
655   IF(ACK*CHIME-(ACK+1.0)*CHIOLD+CHISQ)660,680,660
  660 CINDER=(0.5/ACK)*(ACK**2*CHIME-(ACK**2-1.0)*CHIOLD-CHISQ)/
     *   (ACK*CHIME-(ACK+1.0)*CHIOLD+CHISQ)
      X(I)=SAVE+CINDER*DEL
      CALL FUNK
      NF=NF+1
      IF(CHISQ-CHIOLD)670,680,680
  670 CHIOLD=CHISQ
      TRIAL(I)=TRIAL(I)+CINDER*DEL
      VEC(I)=VEC(I)+CINDER*DEL/ABS (DX(I))
      VECFLG(I) = 0.0
      GO TO 690
  680 X(I)=SAVE
  690 IF(NZIP-1)1340,700,700
  700 IF(ABS (VEC(I))-ACK)750,710,710
  710 DX(I)=ACK*ABS (DX(I))
      VEC(I)=VEC(I)/ACK
      OLDVEC(I)=OLDVEC(I)/ACK
      DO 720 J=1,MOSQUE
  720 ERR(I,J)=ERR(I,J)/ACK
      IF(NTRACE)750,750,730
  730 WRITE(6,740) I,DX(I)
  740 FORMAT(' STEP SIZE'I3,' INCREASED TO 'E12.5)
  750 SUMO=0.0
      SUMV=0.0
      DO 760 J=1,NV
      SUMO=SUMO+OLDVEC(J)**2
  760 SUMV=SUMV+VEC(J)**2
      IF(SUMO*SUMV)1340,1340,770
  770 SUMO=SQRT (SUMO)
      SUMV=SQRT (SUMV)
      COSINE=0.0
      DO 780 J=1,NV
  780 COSINE=COSINE+(OLDVEC(J)/SUMO)*(VEC(J)/SUMV)
      IF(NZIP-1)1340,790,800
  790 IF(NACK-NACTIV)1340,820,820
  800 IF(NACK-NACTIV)820,810,810
  810 IF(NZIP-NCOMP)820,830,830
  820 IF(COSINE-COMPAR)1340,830,830
C
C  SIMON SAYS, TAKE AS MANY GIANT STEPS AS POSSIBLE...
C
  830 IF(NTRACE)860,860,840
  840 WRITE(6,850) CHIOLD,NF,(VEC(J),J=1,I)
  850 FORMAT(/' CHISQ ='E15.8,'   NF = ',I6,
& /' # OF STEPS ='5F9.2/(13X5F9.2))
  860 NGIANT=0
      NTRY=0
      NRETRY=0
C                    POINTER FOR OSC CHECK
      KL=1
      NOSC=NOSC+1
      IF(NOSC-MOSQUE)890,890,870
  870 NOSC=MOSQUE
C                    PUSH DOWN OSC INFO QUEUE AND LOOSE BOTTOM...
      DO 880 K=2,MOSQUE
      CHIOSC(K-1)=CHIOSC(K)
      DO 880 J=1,NV
      XOSC(J,K-1)=XOSC(J,K)
  880 ERR(J,K-1)=ERR(J,K)
  890 DO 900 J=1,NV
      XOSC(J,NOSC)=X(J)
  900 ERR(J,NOSC)=VEC(J)/SUMV
      CHIOSC(NOSC)=CHIOLD
      IF(NOSC-3)960,910,910
910   COXCOM=0.0
C
C SEARCH FOR A PREVIOUS SUCCESSFUL GIANT STEP IN A DIRECTION MORE
C  NEARLY PARALLEL TO THE DIRECTION OF THE PROPOSED STEP THAN WAS THE
C  IMMEDIATELY PREVIOUS ONE.
C
      DO 920 J=1,NV
  920 COXCOM=COXCOM+ERR(J,NOSC)*ERR(J,NOSC-1)
      NAH=NOSC-2
  930 NTRY=0
      DO 950 K=KL,NAH
      NRETRY=NAH-K
      COSINE=0.0
      DO 940 J=1,NV
  940 COSINE=COSINE+ERR(J,NOSC)*ERR(J,K)
      IF(COSINE-COXCOM)950,950,970
  950 CONTINUE
  960 CHIBAK=CHI(I)
      GO TO 1020
  970 NTRY=1
C                    OSC DETECTED...
      KL=K+1
      IF(NTRACE)1000,1000,980
  980 NT=NOSC-K
      WRITE(6,990) NT
  990 FORMAT(/1X,8('*'),5X,' POSSIBLE OSCILLATION WITH PERIOD ',I2,
&         ' DETECTED.')
 1000 DO 1010 J=1,NV
      SALVO(J)=TRIAL(J)
1010  TRIAL(J)=(X(J)-XOSC(J,K))/ACK
      CHIBAK=CHIOLD+(CHIOSC(K)-CHIOLD)/ACK
C
 1020 DO 1040 J=1,NV
      XSAVE(J)=X(J)
      TRIAL(J)=ACK*TRIAL(J)
      IF(MASK(J))1040,1030,1040
1030  X(J)=AMAX1(AMIN1(X(J)+TRIAL(J),XMAX(J)),XMIN(J))
 1040 CONTINUE
      JOCK=0
      JVARY=0
      CALL FUNK
      NF=NF+1
      IF(CHISQ-CHIOLD)1050,1080,1080
 1050 CHIBAK=CHIOLD
      CHIOLD=CHISQ
      NGIANT=NGIANT+1
      IF(NTRACE)1020,1020,1060
 1060 WRITE(6,1070) CHISQ,NF,(X(J),J=1,NV)
 1070 FORMAT(/' CHISQ ='E15.8,'   NF = ',I6,
&    /' X(I)....'/(5(1XE12.5)))
      GO TO 1020
C
 1080 IF(NRETRY)1095,1095,1090
 1090 IF(NGIANT)1150,1150,1095
C                    PERFORM PARABOLIC INTERPOLATION
1095  IF(ACK*CHIBAK-(ACK+1.0)*CHIOLD+CHISQ)1100,1150,1100
 1100 CINDER=(0.5/ACK)*(ACK**2*CHIBAK-(ACK**2-1.0)*CHIOLD-CHISQ)/
&   (ACK*CHIBAK-(ACK+1.0)*CHIOLD+CHISQ)
      DO 1120 J=1,NV
      IF(MASK(J))1120,1110,1120
1110  X(J)=AMAX1(AMIN1(XSAVE(J)+CINDER*TRIAL(J),XMAX(J)),XMIN(J))
 1120 CONTINUE
      JOCK=0
      JVARY=0
      CALL FUNK
      NF=NF+1
      IF(CHISQ-CHIOLD)1280,1130,1130
 1130 IF(NGIANT)1170,1140,1170
 1140 IF(NTRY)1150,1170,1150
 1150 DO 1160 J=1,NV
      TRIAL(J)=SALVO(J)
 1160 X(J)=XSAVE(J)
      GO TO 1190
 1170 DO 1180 J=1,NV
      TRIAL(J)=TRIAL(J)/ACK
 1180 X(J)=XSAVE(J)
 1190 IF(NTRACE)1240,1240,1200
 1200 WRITE(6,1210) CHIOLD,NGIANT
 1210 FORMAT(/' CHISQ ='E15.8,'  AFTER'I3,' GIANT STEPS.')
      WRITE(6,1220) (X(J),J=1,NV)
 1220 FORMAT(' X(I)....'/(5(1XE12.5)))
      WRITE(6,1230)
 1230 FORMAT(/)
 1240 IF(NGIANT)1250,1250,1310
 1250 IF(NRETRY)1260,1260,930
 1260 IF(NTRY)1270,1330,1270
 1270 NTRY=0
      GO TO 960
C
 1280 CHIOLD=CHISQ
      JOCK=1
      IF(NTRACE)1310,1310,1290
 1290 STEPS=FLOAT (NGIANT)+CINDER
      WRITE(6,1300) CHIOLD,STEPS
 1300 FORMAT(/' CHISQ ='E15.8,' AFTER'F6.1,' GIANT STEPS. ')
      WRITE(6,1220) (X(J),J = 1,NV)
      WRITE(6,1230)
 1310 IF(NTRY)1320,380,1320
 1320 NOSC=0
      GO TO 380
 1330 NOSC= MAX0 (NOSC-1,0)
 1340 CHI(I)=CHIOLD
 1350 CONTINUE
C
C  ANOTHER CYCLE THROUGH THE VARIABLES HAS BEEN COMPLETED.
C  PRINT ANOTHER LINE OF TRACES.
C
      IF(NTRACE)1380,1380,1360
 1360 WRITE(6,850) CHIOLD,NF, (VEC(J),J=1,NV)
 1380 IF(NZIP)1420,1390,1420
 1390 IF(NTRACE)1420,1420,1400
 1400 WRITE(6,1220) (X(J),J=1,NV)
      WRITE(6,1410)
 1410 FORMAT(' ')
C
 1420 NZIP=NZIP+1
      GO TO 390
C
C  A MINIMUM HAS BEEN FOUND.  PRINT THE REMAINING TRACES.
C
 1430 IF(NTRACE)1450,1450,1440
 1440 WRITE(6,850) CHIOLD,NF, (VEC(J),J=1,NV)
 1450 IF(NTRACE)1480,1480,1460
 1460 WRITE(6,1220) (X(J),J=1,NV)
      WRITE(6,1230)
C
C  DECREASE THE SIZE OF THE STEPS FOR ALL VARIABLES.
C
 1480 NOSC=0
      NGATE=1
      DO 1520 J=1,NV
C     IF(MASK(J))1520,1490,1520
      IF(MASK(J))1520,1481,1520
 1481 IF(VEC(J)) 1490,1483,1490
 1483 IF(VECFLG(J)) 1500,1500,1520
1490  IF(AMAX1(VEC(J),SIGN(1.0,VEC(J))))1500,1520,1500
 1500 IF(ABS (DX(J))-ABS (DELMIN(J)))1520,1520,1510
 1510 NGATE=0
 1520 DX(J)=DX(J)/RATIO
 1530 IF(NTRACE) 1590, 1590,1540
 1540 WRITE(6,1550) (DX(J),J=1,NV)
 1550 FORMAT(38(' *')//' STEP SIZES REDUCED TO....'//(5(1XE12.5)))
      WRITE(6,1560)
 1560 FORMAT(//)
1590  IF(NGATE.EQ.1)GO TO 1600
      IF(NF.LT.NFB)GO TO 380
 1602 WRITE(6,1604)
 1604 FORMAT( //' NF .GT. NFB '/)
 1600 CHISQ=CHIOLD
	if(ntrace)2190,1610,1610
 1610 WRITE(6,1620) NF
 1620 FORMAT(// 1X,I5,' FUNCTION COMPUTATIONS ')
C
 2190 JVARY=0
      CALL FUNK
      IF(NTRACE)2230,2200,2200
 2200 WRITE(6,2210) (X(J),J=1,NV)
 2210 FORMAT(///10X' FINAL VALUES OF X(I)....'//(5(E15.8)))
      WRITE(6,2220)CHISQ
 2220 FORMAT(// ' FINAL VALUE OF CHISQ =  'E15.8//)
 2230 RETURN
	END

@//E*O*F s.f//
chmod u=r,g=,o= s.f
 
echo x - stepit.doc
sed 's/^@//' > "stepit.doc" <<'@//E*O*F stepit.doc//'
                             SUBROUTINE STEPIT

'STEPIT' is a Fortran II subroutine for finding local minima of almost
any real function (called CHISQ) of several variables (called X(I)).
The X(I) might, for example, be parameters appearing nonlinearly in an
expression which is to be fitted to data by minimizing chi-square.
Other uses include the solution of systems of transcendental
equations by minimizing the sum of squares of the residuals, the 
solution of nearly redundant systems of linear equations, Hartree-
Fock wave function calculations, fitting of bubble-chamber events, etc.,
i.e. any problem for which some measure of the poorness of a trial
solution is a continuous and reasonably well behaved real function of
parameters in the solution.  See the reference by Hook and Jeeves below.
STEPIT uses only function values (no derivatives).  

The function to be minimized must be smooth.  That is, it must be a
continuous function defined on a continuous domain, and all of its first
derivatives must be continuous.  For example, the function 
ABSF(X(1))+ABSF(X(2)) (an inverted pyramid) is not admissible because
its first derivatives are not continuous at the edges of the pyramid.

Author...  J.P. Chandler  I.U.  Physics Dept.  Bloomington, Ind.
All criticisms and suggestions will be welcomed.              
Subroutine Stepit is a copyrighted program (Copyright 1965 J. P. Chandler).     

Stepit is available from the Quantum Chemistry Exchange (QCPE),
Dept. of Chemistry, Indiana University, Bloomington, Indiana.
It may be used by anyone so long as they include an acknowledgement
in any publication of work in which it is used.
Another minimization program (mentioned below), well suited to problems
with complicated constraints, is also available.

BIBLIOGRAPHY....

Direct Search Solution of Numerical and Statistical Problems,
    Robert Hook and T.A. Jeeves, J, Assoc. For Computing Machinery 8, 212.

Non-Linear Estimation,
    IBM Mathematics and Applications Dept., Program Manual MA-3.

A Simplex Method For Function Minimization,
    J.A. Nedler and R. Mead, The Computer Journal 7, 308.

Notes On Statistics For Physicists,
    Jay Orear, Lawrence Radiation Lab. Report UCRL-8417 (U. of California).

Analysis of Experiments in Particle Physics,
    Frank T. Solmitz, Annual Review of Nuclear Science 14, 375.    


USAGE....

The call statement is	CALL STEPIT(FUNK)
FUNK may be any name of a function computation subroutine, of which
there may be any number. An 'external card'
is required in the program calling STEPIT for each subroutine.

All program linkage is via labelled common.
The common statement follows....

	COMMON/STEP/ NV,NTRACE,MATRIX,CHISQ,MASK(100),X(100),XMAX(100)
&	XMIN(100),DELTAX(100),DELMIN(100),ERR(100,15)

CHISQ is the value of the function to be minimized.  The calculation of 
CHISQ is the sole purpose of the subroutine named in the call to Stepit.

NV is the number of variables, X(I).  With this common statement NV
may be as large as 100.
To increase the maximum number of variables it is only necessary to
change NVMAX (an internal parameter of Stepit) and the dimension 
and common statements, provided storage does not become a problem.
(Instructions are given below for handling larger cases.)

Nonzero NTRACE causes a trace map of the minimization process to be 
printed. Initial and final values are printed in any case.

A nonzero mask(I) causes X(I) to be fixed.

X and DELTAX contain the initial values of the 
variables and the step sizes, and the DELMIN the step size cutoffs.
parameters are roughly the same.                                                Cutting of step size continues until all of the step sizes are less
than the corresponding DELMIN(I), for those DELMIN which are nonzero.
Note well that the DELMIN(I) do not, repeat *not*, give the accuracy 
of the final values of the X(I).  The DELMIN(I) should be much 
smaller than the desired accuracy of minimization, and may all be set
equal to zero if maximum accuracy is desired and if no parameter has
its optimum value at exactly zero.

XMAX(I) and XMIN(I) are limits on X(I).  Generally they should be used
only for excluding regions in which the final function is complex, 
discontinuous, or otherwise uncomputable (and not merely unreasonable).
Their use can prevent the location of minima even in the acceptable
region, if the valley leaves but re-enters the region, because the function
is by definition positive infinite outside the boundaries and this can 
introduce new minima.
If limits are used and the final position of the minimum is against one 
of them, any errors calculated will probably be meaningless.
If equal or reversed, XMAX and XMIN will be ignored.

All reference below to error calculation assumes that CHISQ actually is
chi-square (or, in general, twice the negative of the natural log of the
likelihood function).  It is further assumed that the problem is
heteroscedastic, that is, that the errors in the data points are known
and are used in the calculation of the chi-square.  It is assumed that users
running homoscedastic problems (unknown errors, assumed equal) know how to 
scale the errors correctly.  This scaling cannot be done in Stepit because
the number of degrees of freedom is not available internally.
For definitions of chi-square and the likelihood function, see the 
references by Orear or Solmitz.


also see:

Chandler, J.P. Subroutine STEPIT - Finds local minima of
	a smooth function of several parameters.,
	_Behavioral Science_, 1969, _14_, 81-82.
@//E*O*F stepit.doc//
chmod u=r,g=,o= stepit.doc
 
echo x - steps.f
sed 's/^@//' > "steps.f" <<'@//E*O*F steps.f//'
	subroutine steps(nva,xmn,xmx,dlx,dlm,xx,msk,ntr,nf)
	COMMON /STEP/ NV,NTRACE,MATRIX,CHISQ,MASK(100),X(100),XMAX(100),
&	XMIN(100),DELTAX(100),DELMIN(100),ERR(100,15)
	COMMON /TEST/ NFB,CRIT
	NV=nva
	if(nv.le.0) stop
	DO 700 J=1,NV
700	call stepp(j,xmn,xmx,dlx,dlm,xx,msk)
	NTRACE=ntr
	MATRIX=105
	NFB=nf
	return
	END
	subroutine stepp(j,xmn,xmx,dlx,dlm,xx,msk)
	COMMON /STEP/ NV,NTRACE,MATRIX,CHISQ,MASK(100),X(100),XMAX(100),
&	XMIN(100),DELTAX(100),DELMIN(100),ERR(100,15)
	XMIN(J)=xmn
	XMAX(J)=xmx
	DELTAX(J)=dlx
	DELMIN(J)=dlm
	x(j)=xx
	MASK(J)=msk
	return
	end
	subroutine assign(n,name)
	character name*15 
	open(n, file=name)
	rewind(n)
	return
	end
	subroutine sas(fd,fl,ff)
	character fd*15,fl*15,ff*15
	call assign(9,fd)
	call assign(7,fl)
	call assign(4,ff)
	return
	end	
	subroutine clear(x,n)
	dimension x(1)
	do 20 j=1,n
20	x(j)=0.
	return
	end
	subroutine cpvec(a,b,n)
	dimension a(1),b(1)
	do 20 j=1,n
20	b(j)=a(j)
	return
	end
	subroutine admean(x,m,s,n)
c	add 1 nth of the m things in x(...) to s(...)
	dimension x(1),s(1)
	xn=n
	do 20 j=1,m
20	s(j)=s(j)+x(j)/xn
	return
	end
	subroutine means(x,m,s,n)
c	scrunch x(...sub) to s(...)
c	m things/sub, n subs
	dimension x(1),s(1)
	call clear(s,m)
	do 20 j=1,n
	jj=(j-1)*m+1
20	call admean(x(jj),m,s,n)
	return
	end
	subroutine runstp(funk,nnv,nd,nsa,d,do,dp,ppar)
	COMMON /STEP/ NV,NTRACE,MATRIX,CHISQ,MASK(100),X(100),XMAX(100),
&	XMIN(100),DELTAX(100),DELMIN(100),ERR(100,15)
	save
	real d(1),do(1),dp(1),sx(100)
	external funk,ppar
	ns=nsa
	if(ns.lt.0)ns=-ns
	rr=0.
	call clear(sx,nnv)
	nn=ns+1
	do 800 is=1,nn
	if(is.eq.nn.and.(nn.eq.2.or.nsa.lt.0))go to 800
	if(is.lt.nn)call cpvec(d((is-1)*nd+1),do,nd)
	if(is.eq.nn)call means(d,nd,do,ns)
	call steps(nnv,.0001,.9999,.03,.0,.5,0,-1,30000)
	call stepit(funk)
	rmsd=sqrt(chisq/nd)
	if(is.lt.nn)call admean(x,nnv,sx,ns)
	if(is.lt.nn)rr=rr+rmsd/(1.*ns)
	write(4,114)(do(j),j=1,nd),(dp(j),j=1,nd)
114	format(6f11.5)
	call ppar(is,rmsd)
	call corp7(do,dp,nd,' obs',' pre')
800	continue
	if(nn.eq.2.or.nsa.lt.0)return
	call cpvec(sx,x,nnv)
	call ppar(10000,rr)
	return
	end
	subroutine cchi(do,dp,nd)
	real do(1),dp(1)
	COMMON /STEP/ NV,NTRACE,MATRIX,CHISQ,MASK(100),X(100),XMAX(100),
&	XMIN(100),DELTAX(100),DELMIN(100),ERR(100,15)
	chisq=0.
	do 900 j=1,nd
900	chisq=chisq+(do(j)-dp(j))**2
	return
	end
	subroutine cor(x,y,n,rxy,r2xy,ts,xm,ym,byx)
	real x(1),y(1)
	sxy=0.
	sx=0.
	sy=0.
	sx2=0.
	sy2=0.	
	do 500 j=1,n
	sxy=sxy+x(j)*y(j)
	sx=sx+x(j)
	sy=sy+y(j)
	sx2=sx2+x(j)*x(j)
500	sy2=sy2+y(j)*y(j)
	if(sx.eq.0.0.or.sy.eq.0.0)go to 600
	xn=n
	rxy=(xn*sxy-sx*sy)/sqrt((xn*sx2-sx*sx)*(xn*sy2-sy*sy))
	r2xy=rxy*rxy
	ts=1000.
	if(r2xy.lt..99999)ts=rxy*sqrt(xn-2.)/sqrt(1.-r2xy)
	byx=(xn*sxy-sx*sy)/(xn*sx2-sx*sx)
590	xm=sx/xn
	ym=sy/xn
	return
600	rxy=0.
	r2xy=0.
	ts=0.
	byx=0.
	go to 590
	end
	subroutine corp7(x,y,n,lx,ly)
	character*4 lx,ly
	call cor(x,y,n,rxy,r2xy,t,xm,ym,byx)
	ndf=n-2
	write(7,101)lx,ly,rxy,r2xy,t,ndf
101	format(a4'~'a4': r='f6.3' r2='f6.3' t='f9.2' df='i5)
	return
	end
	subroutine corp(x,y,n,lx,ly)
	character*4 lx,ly
	call cor(x,y,n,rxy,r2xy,t,xm,ym,byx)
	ndf=n-2
	write(6,101)lx,ly,rxy,r2xy,t,ndf
101	format(a4'~'a4': r='f6.3' r2='f6.3' t='f7.3' df='i5)
	return
	end
	function ts(x,y,n)
	real x(1),y(1)
	sd=0.
	s2=0.
	do 800 j=1,n
	z=x(j)-y(j)
	sd=sd+z
800	s2=s2+z*z
	ts=sd/sqrt((s2*n-sd*sd)/(n-1.))
	return
	end
@//E*O*F steps.f//
chmod u=r,g=,o= steps.f
 
exit 0