# File rvsao/Emvel/emmin.x
# June 14, 1994
# By Doug Mink, Harvard-Smithsonian Center for Astrophysics
# After John Tonry ( MINI Revision 2.0, 11/17/82)

# Copyright(c) 1994 Smithsonian Astrophysical Observatory
# You may do anything you like with this file except remove this copyright.
# The Smithsonian Astrophysical Observatory makes no representations about
# the suitability of this software for any purpose.  It is provided "as is"
# without express or implied warranty.

define	MAX_PARMS	16
define	MAX_PARMS2	256

#  Black box minimization and fitting program.

procedure emmin (y,ycont,npar,a,debug,covar,ierr)

real	y[ARB]		# Data vector including data to be fit
real	ycont[ARB]	# Continuum vector including region to be fit
int	npar		# Number of parameters to be varied in searching
			#  for a minimum.
double	a[ARB]		# This vector serves three purposes:
			#  it passes the initial estimate for the parameters
			#  it passes arguments to the function to be minimized
			#  it returns the minimum parameter values
bool	debug		# Print diagnostic messages if TRUE
double	covar[ARB]	# RETURNED diagonal of covariance matrix 
int	ierr		# RETURNED error condition:
			#  0 for no error
			#  1 for maximum iteration exceeded
			#  2 for singular matrix
			#  3 for parameter value exceeded

#     descriptions of some of the variables:
#     a - argument for the function
#     a0 - current guess for the minimum
#     da - vector from a0 to new guess for the minimum
#     df - first derivatives of the function
#     d2f - second derivatives of the function
#     lambda - governs mix of gradient and analytic searches
#     iter - maximum number of iterations
#     qfrac - maximum fractional change for successful exit
#
#     The calling program should be as follows (eg):
#********************************************************
#     real rdata[ARB]
#     double a[4]
#     double covar[4]
#     ... (initialize a to the guess for the minimum)
#     call emmin (rdata,4,a,debug,covar,ier)
#     ...
#     procedure chi(a)
#     double a[4], chi
#     ... (define the function)
#     procedure dchi(a,df)
#     double a[4], df[4]
#     ... (dummy or compute the derivatives)
#     procedure d2chi(a,cov)
#     double a[4], cov[4,4]
#     ... (dummy or compute the derivatives)
#************************************************************

double	df[MAX_PARMS],a0[MAX_PARMS],da[MAX_PARMS]
double	d2f[MAX_PARMS,MAX_PARMS],cov[MAX_PARMS,MAX_PARMS]
double	fnow, fthen, qfrac, dfrac, places, vary, det
double	err, base, blambda
int	lambda,itermax,j,k, ii, iter

begin

#  quit if npar is greater than the space allocated
	if (npar > MAX_PARMS) {
	    call printf ("EMMIN: too many parameters, maximum = MAX_PARMS\n")
	    return
	    }
	if (debug) {
	    call printf ("EMMIN: %d-parameter fit\n")
		call pargi (npar)
	    }

#  define a few parameters
	lambda = -3
	base = 10.d0
#	itermax = 20
	itermax = 200
	qfrac = 1.d-6
	dfrac = .02d0
	places = 1.d-7
	vary = 1.d-5

# Initialize a0
	call aclrd (a0,MAX_PARMS)
	call amovd (a, a0, npar)
	call aclrd (cov,MAX_PARMS2)
	call aclrd (d2f,MAX_PARMS2)

	iter = 0
	call gauss (y,ycont,a,debug,fnow)
	fthen = fnow
	if (debug) call vprint (npar,a0,fthen,fnow,iter,lambda)

# Begin iteration to find minimum
	do iter = 1, itermax {
	    fthen = fnow

	# Compute the function derivatives.
	    call amovd (a0, a, npar)

	# Compute 1st and approximate 2nd derivatives
	    call dgauss (y,ycont,a,df,d2f)

	# Begin loop to find a direction along which function decreases
	    do ii = 1, 15 {

	# Get weight matrix
		blambda = 1.d0 + base**lambda
		do j = 1, npar {
		    do k = 1, j-1 {
			cov[j,k] = d2f[j,k]
			cov[k,j] = d2f[j,k]
			}
		    cov[j,j] = d2f[j,j] * blambda
		    da[j] = df[j]
		    }
	# Matrix solution using Gauss-Jordan elimination
		call jgauss (cov,da,npar,MAX_PARMS,ierr)
		if (ierr > 0)
		    return

	# Get new function value
		do j = 1, npar {
		    a[j] = a0[j] - da[j]
		    }
		call gauss (y,ycont,a,debug,fnow)

	# Test for whether the function has decreased
	# If so, adopt the new point and decrement lambda
	# else, increment lambda, and get a new weight matrix

		if (fnow < fthen) break
		lambda = lambda + 1
		if (debug) call vprint(npar,a,fthen,fnow,iter,lambda)
		}

	# Normal exit, the function at a0 + da is less than at a0

	    call amovd (a,a0,npar)
	    lambda = lambda - 1

	# Print the current status and test to see whether the function
	# has varied fractionally less than qfrac.

	    if (debug) call vprint(npar,a0,fthen,fnow,iter,lambda)
	    if (dabs ((fthen - fnow)/fnow) < qfrac) break
	    }

# This is the final computation of the covariance matrix
# Quit if no minimum was found in the allowed number of iterations

	if (iter > itermax) {
	    if (debug)
		call printf ("EMMIN: maximum iteration exceeded\n")
	    ierr = 1
	    }
	else
	    ierr = 0

# Finally, compute the covariance matrix

	do j = 1, npar {
	    do k = 1, npar {
	        cov[j,k] = d2f[j,k]
		}
	    }
	call jgauss (cov,da,npar,MAX_PARMS,ierr)
	do j = 1, npar {
	    err = sqrt (abs (cov[j,j]))
	    if (cov[j,j] < 0) err = -err
	    cov[j,j] = err
	    }
	do j = 2, npar {
	    do k = 1, j-1 {
	        cov[j,k] = cov[j,k] / (cov[j,j]*cov[k,k])
	        cov[k,j] = cov[j,k]
		}
	    }
	do j = 1, npar {
	    covar[j] = cov[j,j]
	    }

	return
end


# simple subroutine to print the current parameters

procedure vprint (n,a,fthen,fnow,niter,lambda)

int	n
double	a[ARB]
double	fthen,fnow
int	niter
int	lambda

int	i

begin
	call printf (" a(i) = ")
	do i = 1, n {
	    call printf ("%10.4f ")
		call pargd (a[i])
	    }

	call printf ("\n f = %15.7g -> %15.7g   iter = %d   lambda = %d\n")
	    call pargd (fthen)
	    call pargd (fnow)
	    call pargi (niter)
	    call pargi (lambda)
	return
end
# Oct 31 1991	New program

# Apr 22 1994	Drop length of spectrum from GAUSS, DGAUSS, and D2GAUS calls
# Apr 28 1994	Parameterize maximum number of parameters
# Apr 29 1994	Make sure temporary arrays are freed
# May  2 1994	Drop derivative type as argument; set it internally
# May  5 1994	Drop debugging output per iteration
# May 24 1994	Add parameter limits
# May 27 1994	Return with error flag set if parameter limits are exceeded
# Jun 13 1994	Drop parameter limits; pass continuum through
# Jun 13 1994	Drop options for analytical and numerical 2nd derivatives
# Jun 14 1994	Use Gauss-Jordan instead of simple matrix inversion