/* This file is part of Cloudy and is copyright (C) 1978-2004 by Gary J. Ferland.
 * For conditions of distribution and use, see copyright notice in license.txt */
/*optimize_func actual function called during evaluation of optimization run */
#include "cddefines.h"
#include "zero.h"
#include "lines.h"
#include "prt.h"
#include "called.h"
#include "radius.h"
#include "input.h"
#include "cloudy.h"
#include "cddrive.h"
#include "optimize.h"
/* used below to derive chi2 */
double chi2_func(double,double,double);

double optimize_func(float param[])
{

#define MAXCAT 10

	char chCAPLB[5], 
	  chFind[5];

	int lgBAD,
		lgHIT, 
	  lgLimOK;

	long int cat,
	  i, 
	  j, 
	  nfound, 
	  nobs_cat[MAXCAT],
	  np;

	static long int ipobs[NOBSLM];

	double chi1, 
	  chi2_cat[MAXCAT],
	  chisq, 
	  func_v, 
	  help,
	  predin,
	  scld,
	  snorm,
	  theocl;

	static char name_cat[MAXCAT][13] = 
	{
	  "rel flux    ",
	  "column dens ",
	  "abs flux    "
	};

	static int lgLinSet = FALSE;

#	ifdef DEBUG_FUN
	fputs( "<+>optimize_func()\n", debug_fp );
#	endif
	/* limpar defined in optimize.h */

	/* This routine is called by AMOEBA (or POWELL) with values of the
	 * variable parameters for CLOUDY in the array p(i). It returns
	 * the value FUNC = SUM (obs-model)**2/sig**2 for the lines
	 * specified in the observational data file, values held in the
	 * common blocks /OBSLIN/ & /OBSINT/
	 * replacement input strings for CLOUDY READR held in /chCardSav/
	 * parameter information for setting chCardSav held in /parmv/
	 * additional variables
	 * Gary's variables
	 */

	/* initialize the code for this run */
	/*cdInit();*/
	/* zero out lots of variables */
	zero();

	if( optimize.lgOptimFlow )
	{
		fprintf( ioQQQ, " trace, optimize_func variables" );
		for( i=1; i <= optimize.nvary; i++ )
		{
			fprintf( ioQQQ, "%10.2e", param[i-1] );
		}
		fprintf( ioQQQ, "\n" );
	}

	for( i=0; i < optimize.nvary; i++ )
	{
		optimize.vparm[0][i] = param[i];
	}

	/* allways increment nOptimiz, even if parameters are out of bounds,
	 * this prevents optimizer to get stuck in infinite loop */
	optimize.nOptimiz += 1;

	/* call routine to pack /kardsv/ variable with appropriate
	 * CLOUDY input lines given the array of variable parameters p(i) */
	optimize_input(&lgLimOK);

	if( !lgLimOK )
	{
		/* these parameters are not within limits of parameter search
		 * >>chng 96 apr 26, as per Peter van Hoof comment */
		fprintf( ioQQQ, " Iteration%4ld not within range.\n", 
		  optimize.nOptimiz );

		/* this is error; very bad since not within range of parameters */
		func_v = (double)FLT_MAX;
		
#		ifdef DEBUG_FUN
		fputs( " <->optimize_func()\n", debug_fp );
#		endif
		return( func_v );
	}

	for( i=0; i < optimize.nvary; i++ )
	{
		optimize.varmax[i] = (float)MAX2(optimize.varmax[i],optimize.vpused[i]);
		optimize.varmin[i] = (float)MIN2(optimize.varmin[i],optimize.vpused[i]);
	}

	lgBAD = cloudy();
	if( lgBAD )
	{
		fprintf( ioQQQ, " Cloudy returned error condition - what happened?\n" );
	}

	/* Line fluxes now in commn blocks, so extract and
	 * compare with observations */
	chisq = 0.0;
	for( i=0; i < MAXCAT; i++ )
	{
		nobs_cat[i] = 0;
		chi2_cat[i] = 0.0;
	}

	if( LineSave.ipNormWavL < 0 )
	{
		fprintf( ioQQQ, 
			" Normalization line array index is bad.  What has gone wrong?\n" );
		puts( "[Stop in optimize_func]" );
		cdEXIT(EXIT_FAILURE);
	}
	snorm = LineSv[LineSave.ipNormWavL].sumlin;

	if( snorm == 0. )
	{
		fprintf( ioQQQ, " Normalization line has zero intensity.  What has gone wrong?\n" );
		fprintf( ioQQQ, " Is spectrum normalized to a species that does not exist?\n" );
		puts( "[Stop in optimize_func]" );
		cdEXIT(EXIT_FAILURE);
	}

	/* first print all warnings */
	cdWarnings(ioQQQ);
	fprintf( ioQQQ, "  ID              Model    Observed       error      chi**2     Type\n" );
	/* cycle through the observational values */
	nfound = 0;

	/* first is to optimize relative emission line spectrum */
	if( optimize.lgOptLin )
	{
		/* set pointers to all optimized lines if first call */
		if( !lgLinSet )
		{
			lgLinSet = TRUE;
			for( i=0; i < optimize.nlobs; i++ )
			{
				cap4(chFind , (char*)optimize.chAmLab[i]);
				lgHIT = FALSE;
				j = 0;

				while( !lgHIT && j < LineSave.nsum )
				{
					/* check wavelength for a match
					 * change label to all caps to be like input chAmLab */
					cap4(chCAPLB , (char*)LineSv[j].chALab);
					/*if( fabs(LineSv[j].wavelength - optimize.wavelength[i])/SDIV(optimize.wavelength[i])< optimize.errorwave[i] && 
					  strcmp(chCAPLB,chFind) == 0 )*/
					if( fabs(LineSv[j].wavelength - optimize.wavelength[i]) < optimize.errorwave[i] && 
					  strcmp(chCAPLB,chFind) == 0 )
					{
						/* match, so set array index */
						ipobs[i] = j;
						lgHIT = TRUE;
					}
					++j;
				}

				/* we did not find the line */
				if( !lgHIT )
				{
					fprintf( ioQQQ, " Optimizer could not find line with label %4.4s and wavelength ", 
					  optimize.chAmLab[i] );
					prt_wl( ioQQQ, optimize.wavelength[i] );
					fprintf( ioQQQ, " . Sorry.\n");
					puts( "[Stop in optimize_func]" );
					cdEXIT(EXIT_FAILURE);
				}
			}
		}

		for( i=0; i < optimize.nlobs; i++ )
		{
			/* and find corresponding model value by straight search */
			nfound += 1;
			scld = (double)LineSv[ipobs[i]].sumlin/(double)snorm*LineSave.ScaleNormLine;
			chi1 = chi2_func(scld,(double)optimize.amint[i],(double)optimize.amierr[i]);
			cat = 0;
			nobs_cat[cat]++;
			chi2_cat[cat] += chi1;

			fprintf( ioQQQ, " %4.4s", 
			  LineSv[ipobs[i]].chALab);

			prt_wl( ioQQQ,LineSv[ipobs[i]].wavelength);

			fprintf( ioQQQ, "%12.5f%12.5f%12.5f%12.2e Relative intensity\n", 
			  scld, 
			  optimize.amint[i], 
			  optimize.amierr[i], 
			  chi1 );
		}
	}

	/* option to optimize column densities */
	if( optimize.lgOptCol )
	{
		for( i=0; i < optimize.ncobs; i++ )
		{
			if( cdColm((char*)optimize.chColmnLab[i],optimize.ionam[i], &theocl) )
			{
				puts( "[Stop in optimize_func]" );
				cdEXIT(EXIT_FAILURE);
			}
			nfound += 1;
			chi1 = chi2_func(theocl,(double)optimize.amcol[i],(double)optimize.amcerr[i]);
			cat = 1;
			nobs_cat[cat]++;
			chi2_cat[cat] += chi1;

			fprintf( ioQQQ, " %4.4s%6ld%12.4e%12.4e%12.5f%12.2e Column density\n", 
			  optimize.chColmnLab[i], optimize.ionam[i], theocl, 
			  optimize.amcol[i], optimize.amcerr[i], chi1 );
		}
	}

	/* option to optimize line flux */
	if( optimize.lgOptLum )
	{
		nfound += 1;
		if( LineSv[LineSave.ipNormWavL].sumlin > 0.f )
		{
			predin = log10(LineSv[LineSave.ipNormWavL].sumlin) + radius.Conv2PrtInten;
			help = pow(10.,predin-(double)optimize.optint);
			chi1 = chi2_func(help,1.,(double)optimize.optier);
		}
		else
		{
			predin = -999.99999;
			chi1 = (double)FLT_MAX;
		}
		cat = 2;
		nobs_cat[cat]++;
		chi2_cat[cat] += chi1;

		fprintf( ioQQQ, " %4.4s%6f%12.5f%12.5f%12.5f%12.2e Line intensity\n", 
		  LineSv[LineSave.ipNormWavL].chALab, 
		  LineSv[LineSave.ipNormWavL].wavelength, 
		  predin, 
		  optimize.optint, 
		  optimize.optier, 
		  chi1 );
	}

	if( nfound < optimize.nlobs )
	{
		fprintf( ioQQQ, " Observables were missed in optimization.\n" );
		puts( "[Stop in optimize_func]" );
		cdEXIT(EXIT_FAILURE);
	}

	if( nfound <= 0 )
	{
		fprintf( ioQQQ, " WARNING; no line matches found\n" );
		puts( "[Stop in optimize_func]" );
		cdEXIT(EXIT_FAILURE);
	}

	/* write out chisquared for this iteration */
	fprintf( ioQQQ, "\n" );
	for( i=0; i < MAXCAT; i++ )
	{
		if( nobs_cat[i] > 0 )
		{
			chisq += chi2_cat[i]/nobs_cat[i];
			fprintf( ioQQQ, " Category %s #obs.%3ld  Total Chi**2%11.3e  Average Chi**2%11.3e\n",
			  name_cat[i],nobs_cat[i],chi2_cat[i],chi2_cat[i]/nobs_cat[i] );
		}
	}

	fprintf( ioQQQ, "\n Iteration%4ld Chisq=%13.5e\n", optimize.nOptimiz, chisq );

	/* only print this if output has been turned on */
	if( called.lgTalk )
	{
		fprintf( ioQQQ, "\n" );
		for( i=0; i < optimize.nvary; i++ )
		{
			optimize.vparm[0][i] = (float)MIN2(optimize.vparm[0][i],optimize.varang[i][1]);
			optimize.vparm[0][i] = (float)MAX2(optimize.vparm[0][i],optimize.varang[i][0]);
			param[i] = optimize.vparm[0][i];
			np = optimize.nvfpnt[i];

			/* now generate the actual command with parameter,
			 * there will be from 1 to 3 numbers on the line */
			if( optimize.nvarxt[i] == 1 )
			{
				/* case with 1 parameter */
				sprintf( input.chCardSav[np] , optimize.chVarFmt[i], optimize.vparm[0][i] );
			}

			else if( optimize.nvarxt[i] == 2 )
			{
				/* case with 2 parameter */
				sprintf( input.chCardSav[np] , optimize.chVarFmt[i], optimize.vparm[0][i], optimize.vparm[1][i]);
			}

			else if( optimize.nvarxt[i] == 3 )
			{
				/* case with 3 parameter */
				sprintf( input.chCardSav[np] , optimize.chVarFmt[i], 
					optimize.vparm[0][i], optimize.vparm[1][i] , optimize.vparm[2][i] );
			}

			else if( optimize.nvarxt[i] == 4 )
			{
				/* case with 4 parameter */
				sprintf( input.chCardSav[np] , optimize.chVarFmt[i], 
					optimize.vparm[0][i], optimize.vparm[1][i] , optimize.vparm[2][i], optimize.vparm[3][i] );
			}

			else if( optimize.nvarxt[i] == 5 )
			{
				/* case with 5 parameter */
				sprintf( input.chCardSav[np] , optimize.chVarFmt[i], 
					optimize.vparm[0][i], optimize.vparm[1][i] , optimize.vparm[2][i],
					optimize.vparm[3][i] , optimize.vparm[4][i]);
			}

			else
			{
				fprintf(ioQQQ,"The number of variable options on this line makes no sense to me4\n");
				puts( "[Stop in optimize_func]" );
				cdEXIT(EXIT_FAILURE);
			}

			fprintf( ioQQQ, " Optimal command: %s\n", 
			  input.chCardSav[np] );
		}
	}

	func_v = MIN2(chisq,(double)FLT_MAX);

#	ifdef DEBUG_FUN
	fputs( " <->optimize_func()\n", debug_fp );
#	endif
	return( func_v );
}

/* ============================================================================== */
double chi2_func(double ymodl,
	double ymeas,
	double yerr)
{
	double chi2_func_v,
		temp ;

#	ifdef DEBUG_FUN
	fputs( "<+>chi2_func()\n", debug_fp );
#	endif

	/* compute chi**2 by comparing model quantity ymodl with a measured
	 * quantity ymeas with relative error yerr (negative means upper limit)
	 */

	if( ymeas <= 0. )
	{
		fprintf( ioQQQ, "chi2_func: non-positive observed quantity, this should not happen\n" );
		puts( "[Stop]" );
		cdEXIT(EXIT_FAILURE);
	}

	if( yerr > 0. )
	{
		if( ymodl > 0. )
		{
			temp = POW2((ymodl-ymeas)/(MIN2(ymodl,ymeas)*yerr));
			chi2_func_v = MIN2( temp , (double)FLT_MAX );
		}
		else
			chi2_func_v = (double)FLT_MAX;
	}
	else if( yerr < 0. )
	{
		/* value quoted is an upper limit, so add to chisq
		 * only if limit exceeded, otherwise return zero.
		 */
		if( ymodl > ymeas )
		{
			temp = POW2((ymodl-ymeas)/(ymeas*yerr));
			chi2_func_v = MIN2(temp,(double)FLT_MAX);
		}
		else
			chi2_func_v = 0.;
	}
	else
	{
		fprintf( ioQQQ, "chi2_func: relative error is zero, this should not happen\n" );
		puts( "[Stop]" );
		cdEXIT(EXIT_FAILURE);
	}

#	ifdef DEBUG_FUN
	fputs( " <->chi2_func()\n", debug_fp );
#	endif
	return chi2_func_v;
}
