/*func actual function called during evaluation of optimization run */
#include <string.h>
#include <float.h>
#include "cddefines.h"
#include "zero.h"
#include "func.h"
#include "varypar.h"
#include "linesave.h"
#include "called.h"
#include "norm.h"
#include "fourpi.h"
#include "input.h"
#include "cap4.h"
#include "cloudy.h"
#include "setinp.h"
#include "cddrive.h"

double func(float param[])
{

#define MAXCAT 10

	char chCAPLB[5], 
	  chFind[5];

	int lgHIT, 
	  lgLimOK, 
	  lgOK;

	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( "<+>func()\n", debug_fp );
#	endif
	/* limpar defined in VaryPar.com */

	/* 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( VaryPar.lgOptimFlow )
	{
		fprintf( ioQQQ, " trace, func variables" );
		for( i=1; i <= VaryPar.nvary; i++ )
		{
			fprintf( ioQQQ, "%10.2e", param[i-1] );
		}
		fprintf( ioQQQ, "\n" );
	}

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

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

	/* call routine to pack /kardsv/ variable with appropriate
	 * CLOUDY input lines given the array of variable parameters p(i) */
	setinp(&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", 
		  VaryPar.nOptimiz );

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

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

	cloudy();

	/* 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( norm.ipNormWavL < 0 )
	{
		fprintf( ioQQQ, " Normalization line pointer is bad.  What has gone wrong?\n" );
		puts( "[Stop in func]" );
		exit(1);
	}
	snorm = LineSv[norm.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 func]" );
		exit(1);
	}

	/* 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( VaryPar.lgOptLin )
	{
		/* set pointers to all optimized lines if first call */
		if( !lgLinSet )
		{
			lgLinSet = TRUE;
			for( i=0; i < VaryPar.nlobs; i++ )
			{
				cap4(chFind , (char*)chVaryPar.chAmLab[i]);
				lgHIT = FALSE;
				j = 1;

				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-1].chALab);
					if( LineSv[j-1].lin == VaryPar.linam[i] && 
					  strcmp(chCAPLB,chFind) == 0 )
					{
						/*       match, so set pointer */
						ipobs[i] = j;
						lgHIT = TRUE;
					}
					j += 1;
				}

				/* we did not find the line */
				if( !lgHIT )
				{
					fprintf( ioQQQ, " Optimizer could not find line with label %4.4s and wavelength%5ld. Sorry.\n", 
					  chVaryPar.chAmLab[i], VaryPar.linam[i] );
					puts( "[Stop in func]" );
					exit(1);
				}
			}
		}

		for( i=0; i < VaryPar.nlobs; i++ )
		{
			/* and find corresponding model value by straight search */
			nfound += 1;
			scld = (double)LineSv[ipobs[i]-1].sumlin/(double)snorm*norm.ScaleNormLine;
			chi1 = chi2_func(scld,(double)VaryPar.amint[i],(double)VaryPar.amierr[i]);
			cat = 0;
			nobs_cat[cat]++;
			chi2_cat[cat] += chi1;
			fprintf( ioQQQ, " %4.4s%6ld%12.5f%12.5f%12.5f%12.2e Relative intensity\n", 
			  LineSv[ipobs[i]-1].chALab, LineSv[ipobs[i]-1].lin, 
			  scld, VaryPar.amint[i], VaryPar.amierr[i], chi1 );
		}
	}

	/* option to optimize column densities */
	if( VaryPar.lgOptCol )
	{
		for( i=0; i < VaryPar.ncobs; i++ )
		{
			lgOK = cdColm((char*)chVaryPar.chColmnLab[i],VaryPar.ionam[i],
			  &theocl);
			if( !lgOK )
			{
				puts( "[Stop in func]" );
				exit(1);
			}
			nfound += 1;
			chi1 = chi2_func(theocl,(double)VaryPar.amcol[i],(double)VaryPar.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", 
			  chVaryPar.chColmnLab[i], VaryPar.ionam[i], theocl, 
			  VaryPar.amcol[i], VaryPar.amcerr[i], chi1 );
		}
	}

	/* option to optimize line flux */
	if( VaryPar.lgOptLum )
	{
		nfound += 1;
		if( LineSv[norm.ipNormWavL].sumlin > 0.f )
		{
			predin = log10(LineSv[norm.ipNormWavL].sumlin) + fourpi.pirsq;
			help = pow(10.,predin-(double)VaryPar.optint);
			chi1 = chi2_func(help,1.,(double)VaryPar.optier);
		}
		else
		{
			predin = -999.99999;
			chi1 = (double)FLT_MAX;
		}
		cat = 2;
		nobs_cat[cat]++;
		chi2_cat[cat] += chi1;
		fprintf( ioQQQ, " %4.4s%6ld%12.5f%12.5f%12.5f%12.2e Line intensity\n", 
		  LineSv[norm.ipNormWavL].chALab, 
		  LineSv[norm.ipNormWavL].lin, 
		  predin, 
		  VaryPar.optint, 
		  VaryPar.optier, 
		  chi1 );
	}

	if( nfound < VaryPar.nlobs )
	{
		fprintf( ioQQQ, " Observables were missed in optimization.\n" );
		puts( "[Stop in func]" );
		exit(1);
	}

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

	/* 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", VaryPar.nOptimiz, chisq );

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

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

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

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

			else
			{
				fprintf(ioQQQ,"the number of variable options on this line makes no sense to me3\n");
				puts( "[Stop in func]" );
				exit(1);
			}

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

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


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


double chi2_func(double ymodl,
	double ymeas,
	double yerr)
{
	double chi2_func_v;

#	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]" );
		exit(1);
	}

	if( yerr > 0. )
	{
		if( ymodl > 0. )
			chi2_func_v = MIN2(POW2((ymodl-ymeas)/(MIN2(ymodl,ymeas)*yerr)),(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 )
			chi2_func_v = MIN2(POW2((ymodl-ymeas)/(ymeas*yerr)),(double)FLT_MAX);
		else
			chi2_func_v = 0.;
	}
	else
	{
		fprintf( ioQQQ, "chi2_func: relative error is zero, this should not happen\n" );
		puts( "[Stop]" );
		exit(1);
	}

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