			    SUBROUTINE FITEM
 
	   A GENERAL STRATEGY FOR NON-LINEAR PARAMETER OPTIMISATION
		USING VALUES OF FUZZY FUNCTIONS WITH NOISY DATA

	Three kinds of minimisation search are used in a dynamic variety
	of combinations.  Each kind looks in a different way, and differs
	in its rate of approach towards and convergence on a minimum.


	Given a set of data points and a model function with parameters which
 are to be chosen so that the function will represent the data as closely as
 possible, first construct a"goodness-of-fit" function FUNVAL which connects
 the model function with the data.  Then choose a starting set of parameter
 values PAR, each with a corresponding STEP value which will make a noticeable
 change to FUNVAL.  FITEM will iterate changes of PAR to get a better fit
 (defined to be a minimum of FUNVAL) and will also change each STEP so it
 makes the same change in FUNVAL.  The "precision" of the parameters obtained
 is indicated by the size of their steps in conjunction with reported values for
 the direction cosines between pairs of parameters. 

	Ideally the parameters would be independent, and a minimisation
 technique adopted which was based on a knowledge of the function behaviour
 i.e. the shape of the multi-dimensional surface described by the parameters
 and the data.  In practice, the degree of independence of the parameters and
 the characteristics of the surface are largely unknown.  FITEM provides a
 means of finding out about these properties, and, based on the knowledge
 gained, you can dictate the minimisation strategy to be used.  Alternatively
 you can let FITEM explore the surface and adopt whichever method it deems most
 suitable at any stage.  

	Efficiency ?  This depends on the relative burdens of calculation and
 re-compilation to try alternative procedures.  A free-running FITEM probably
 takes more function evaluations than would an optimised strategy - but it
 has probably obtained the desired result in less time and with less effort
 than are spent in discovering the "best" procedure.  Its argument list is
 structured in such a way that "difficult", "ill-conditioned", or "intractable"
 problems can be investigated from many viewpoints without re-compilation.


 BASIC SEARCH METHODS
 --------------------
	Three kinds of search are used in various combinations.  Each kind
 looks in a different way, and they differ in their rates of approach towards
 and convergence on a minimum when the multi-dimensional surface described by
 the parameters and their function values in relation to the data points is
 wrinkled, i.e. has contorted valleys and local minima.

	Line searches look for a minimum in one direction, using a parabolic
 approximation.  Since a zero slope cannot distinguish between a minimum and
 maximum, the search begins by looking for a "bracket" range, i.e. 3 points such
 that the midpoint is lower than both ends.  In a DESCENT cycle (based on
 Powells's Method of Steepest Descent) a sequence of 2-dimensional sections of
 the surface are examined - the initial search directions are those
 corresponding to single parameters, then a VECTOR is constructed from all the
 parameter changes and a search made in its direction.  Each search ends when
 the position of the parabola minimum is within a small fraction of the midpoint
 of the bracket range - the fraction being chosen so that the function value is
 within a requested fraction of the lowest function value already found.

	The SIMPLEX method () looks at many points on the surface
 simultaneously.  For N parameters it requires N+1 distinct starting points,
 each point being a set of N parameter values.  It makes no assumptions about
 the degree of polynomial which approximates to the shape of the surface.  It
 proceeds by "shrinking" i.e. replacing the highest point by a new-found point
 which is lower, and by "stretching" if it finds a new point lower than any
 found before.  This search terminates when it has replaced the original N+1
 points by ones which all give function values within the requested fractional
 significance.

	The CONTOUR method, instead of looking towards the low point, looks
 outward from it by displacing pairs of parameters to their bracket limits.
 It is unusual for the individual parameter cross sections of the surface to be
 orthogonal, and pair displacements enable correlations between parameters to
 be examined.  If the surface rises everywhere the minimum of a matching second
 degree surface is used to predict a fruitful direction to move; if a pair-
 function falls a Descent or Simplex search is used to follow the downward path.


 CHARACTERISTICS OF THE METHODS
 ------------------------------
	DESCENT goes downhill to a minimum, first looking at PARs one at a
 time, then along a VECTOR made up of all their changes.  It is fast far from
 a minimum but gets "tired",i.e. can't find a stationary set of parameters,
 when some PAR directions are not orthogonal.  The PAR set obtained can depend
 critically on the particular downhill path which has been taken.  The way
 down can be altered by presenting the PARs in a different order and by
 changing the signs of their STEPs; this code delays the onset of "tiring" by
 beginning successive search cycles at a different PAR.

	SIMPLEX moves cautiously by shrinking its highest point toward the
 average of all other points, and accepting a new lowest point when it finds
 one regardless of ups and downs between the low points.  It requires more
 function evaluations than Descent when the parameters are reasonably
 independent, but is faster and gets more consistent results when parameters
 are correlated.  It can wriggle down long, narrow and twisting valleys, but
 gets "tired" because it can't very easily look over ridges to see a lower
 point.  PAR and STEP values at the search start may produce jumps between
 alternative local minima more unexpectedly than with Descent - this may or
 may not be a blessing.

	CONTOUR can get out of a local minimum by looking around at a higher
 level rather than towards the lowest point already found, but if it doesn't
 happen to find a gap in the contour it also tires.  The contour level can be
 changed to explore the surface under different "magnifications".

	In essence, these methods prescribe different ways of sampling the
 surface and use different rules to decide on a subsequent action.  Even
 though a particular search may not lower the function value very much, it
 is often able to dislodge the PARs enough to revitalise one of the other
 types of search.  This effect is enhanced by subsidiary adaptations of the
 main procedures:
	TRICKLE makes a line search in the net direction moved in several
 cycles,
	SWAY looks "sideways" after each vector search has located its minimum
 by removing the largest component from the vector just used and looking in
 that direction
	STEP renormalisation at the end of each cycle changes the size of the
 bracket range for the next cycle.


CALLING THE ROUTINE
-------------------
original:
	SUBROUTINE FITEM(NPAR,PAR,STEP,W,KEEP,SIZE,F,FUNVAL,LOP,MUM,NOUT)
revised:
	SUBROUTINE FIT77(NPAR,PAR,STEP,D,NC,ASIZ,F,FUNVAL,LOP,MM,NOUT)
 
c	FITEM expanded to include simplex and trickle searches, keeping the
c		argument list consistent with established usage.  Internal
c		arrays D,NC replace the old workspace ones W,KEEP
c	Options for search kinds in signs of ASIZ,LOP,MM,NOUT
c	Language is Fortran77 consistent with Vax and MicroSoft

	implicit real*4(A-H,O-Z),integer*4(I-N)
cold 1800 2100 programs need:    integer*2
 
c Arguments: (no constants - any might be altered by FITEM)
c	dimension PAR,STEP to at least NPAR:  first NPAR STEPs are scanned
c		for non-zero values which identify the pars to be optimised
c		- nonzero step count is NP (limited to maxnp)
c	FUNCTION FUNVAL(PAR,NP) is EXTERNAL in calling program
c		It must not alter NP or any PAR(J) with non-zero STEP(J)
c		A -ve value will cause an error return from FITEM 
c	ASIZ - abs(ASIZ) is SIZE used to adjust STEPs to give funval values
c		at a contour level 1+SIZE times lowest funval found so far
c		(SIZE 1 to .01 typically larger than 100/10**LOP)
c		SIZE lt .1**LOP will change ASIZ so check after return
c			ASIZ -ve allows "sideways" look after vector search
c	F   - AT CALL, is a satisfactory low FUNVAL ge 0 at which the search
c		may be terminated
c	    - thereafter is the lowest FUNVAL found (F can be used, but
c		not changed, by FUNVAL)
c	    - AT RETURN, is the lowest FUNVAL found, for returned PAR,
c		and STEPs correspond to SIZE
c	LOP - abs(LOP) le 6 is the number of significant digits for F
c		Over 6 is reduced - check LOP after return
c			LOP -ve suppresses SIMPLEX cycles
c	MM  = 1 for no progress reports from fitem
c	    = 2 for reports on console unit at each cycle end
c	    = 3 for reports as for 2 and at every funval call
c			MM -ve suppresses DESCENT cycles
c	NOUT - AT CALL, abs value is the limit number of cycles,
c			NOUT -ve suppresses CONTOUR cycles
c	     - AT RETURN is one of the following completion codes:
c	   1  no par has changed by more than its step*small, where
c		 small = sqrt(0.1**LOP/SIZE) = 0.1 for SIZE=100*0.1**LOP
c	   2  no significant decrease in funval in the last cycle
c	   3  funval is below the limit set for F at call, or funval=0
c	   4  cycle limit reached
c	   5  all steps zero, no search made, funval has been called
c	   6  -ve funval found, error, search abandoned, and
c			f,par are set for lowest funval found so far
c	D	is unused dummy argument	(was  W   in old FITEM)
c	NC	dimension (2)			(was KEEP in old FITEM)
c	returns NC(1) = number of function calls
c		NC(2) = number of improvement cycles done


c RESTRICTIONS AND LIMITATIONS:

c   1	any of the FITEM arguments can be in COMMON if your FUNVAL needs
c	to refer to them.  The only values which may be changed by FUNVAL
c	are those PAR(J) which have STEP(J)=0.

c   2	FITEM will not accept negative FUNVAL (for historical reasons) -  it
c	operates by locating the minimum of a second degree surface which
c	has been matched to the actual surface.  This approximation is usually
c	good for FUNVALs constructed to look like a sum-of-squared-residuals
c	based on differences between data and model function values.

c   3	PAR changes by +-STEP should increase F by factor 1+SIZE.  When
c	searching to bracket a minimum the SCALE on STEP is doubled when F
c	fails to increase appreciably (and a CONTOUR cycle will revert to
c	DESCENT).  This reduces loss of significant digits for predicting the
c	position of the minimum and ensures a non-singular matrix inversion
c	for the CONTOUR search. In a single cycle, the range expansion is
c	limited to a factor 100 - this shows as a corresponding increase in
c	STEP for the next cycle. If a parameter has no effect on FUNVAL too
c	many cycles will increase STEP until "real overflow" occurs - it is
c	left to the user to correct this situation since FITEM has no prior
c	knowledge of how large or small a PAR can be handled by FUNVAL

c   4   The number of PARs with non-zero STEPs is limited by the dimensions
c	of the internal working arrays.  These are set by the statement
c	PARAMETER (MAXNP=30)


 USING FITEM IS AN ART
 ---------------------
	If you don't want to use FITEM merely as a black box, points to bear
 in mind include the following:

	Has it converged on a local or a global minimum - the only way you
 know there is a better minimum than the one you have found is by finding the
 better one !

	Starting PAR set - some can be chosen by insight, others by solving
 your model function at selected data points.

	A general strategy uses at least 3 successive calls to FITEM
 - 1st one lets the initial guesses for PARs and STEPs settle down. It may be
 limited to a few iterations or only vary those PARs in which you have least
 confidence,
 - 2nd one may let all non-linear PARS be free and run to completion,
 - 3rd should displace the results from the 2nd and see how well they are
 regained.

	Let FITEM teach you - use MUM = 2 or 3 - about the behaviour of your
 model function, and how sensitive FUNVAL is to PARs which may interact with
 only part of the range of data. In return, you will be able to use FITEM more
 effectively by using sub-ranges of the data and making the STEP zero for any
 PAR which has little effect in each data range being examined.

	PARs which are strongly correlated, cosine=1 or -1, far from a
 minimum may be less correlated at a minimum.  If the array of "cosines" has
 a row or column of off-diagonal elements larger than 1 it is a good indication
 that there is a "hole" in the contour nearby.  The route downhill can be
 altered in many ways, thus
 - vary SIZE e.g. between 1 and 0.01,
 - vary LOP  e.g. between 2 and 4
 - change the sign of STEP to begin searching in opposite directions
 - begin with different initial PARs and watch the path FITEM takes
 - suppress one or more of the kinds of search (Contour alone is not possible)

	SIZE and LOP reflect properties of your function values - statistical
 scatter in data and propagation of truncation errors in calculation - they do
 not determine the precision of the results. Making one or both of them smaller
 is more likely to get you stuck in a false or local minimum.  Too large a LOP
 in particular can get Simplex stuck - if it detects inability to shrink to the
 implied precision it calls on Contour or Descent to help it out.

	PARs which are amenable to linear least squares treatment need not be
 varied by FITEM - it is usually quicker to calculate them by a linear method
 each time FUNVAL is called. The choice of parameters to be handled this way
 may be made by having STEP in COMMON.

	CONSTRAINED OPTIMISATION can be handled by using COMMON variables.
 If F and STEP are COMMON then a constraint can impose a "graduated penalty"
 when it is violated e.g. increase the current lowest F by adding an amount
 proportional to the square of the violation in units of the current STEP.
 Beware of creating false minima and of preventing PARs from coming close to
 an invalid bound e.g. using the magnitude or the square of a PAR which should
 not be negative can make a cusp at PAR=zero and FITEM may oscillate between
 the pair of equivalent minima

	FUNVAL and COMMON
 FITEM assumes a definition real*4 FUNCTION FUNVAL(PAR,NP).  It is often
 convenient to have PAR in COMMON in which case the function definition
 statement and COMMON must use different names for PAR.

------------------------------------------------------------------------
