/* This file is part of Cloudy and is copyright (C)1978-2013 by Gary J. Ferland and
 * others.  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 "init.h"
#include "lines.h"
#include "prt.h"
#include "called.h"
#include "predcont.h"
#include "radius.h"
#include "rfield.h"
#include "input.h"
#include "cloudy.h"
#include "cddrive.h"
#include "optimize.h"
#include "grid.h"
/* used below to derive chi2 */
STATIC double chi2_func(double,double,double);

chi2_type optimize_func(const realnum param[],
			int grid_index)
{
	const int MAXCAT = 6;

	static const char name_cat[MAXCAT][13] = 
	{
	  "rel flux    ",
	  "column dens ",
	  "abs flux    ",
	  "mean temp   ",
	  "ang diameter",
	  "photometry  "
	};

	char chFind[5];

	bool lgBAD,
	  lgLimOK;

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

	chi2_type chi1, 
	  chi2_cat[MAXCAT],
	  chisq, 
	  help,
	  predin,
	  scld,
	  snorm,
	  theocl,
	  temp_theory;

	DEBUG_ENTRY( "optimize_func()" );

	if( grid_index >= 0 )
		optimize.nOptimiz = grid_index;

	/* This routine is called by optimizer 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
	 */

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

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

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

	// nothing more to be done...
	if( strcmp(optimize.chOptRtn,"XSPE") == 0 )
		return 0.;

	/* zero out lots of variables */
	zero();

	for( i=0; i < optimize.nvary; i++ )
	{
		optimize.varmax[i] = max(optimize.varmax[i],min(param[i],optimize.varang[i][1]));
		optimize.varmin[i] = min(optimize.varmin[i],max(param[i],optimize.varang[i][0]));
	}

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

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

		/* this is error; very bad since not within range of parameters */
		return BIG_CHI2;
	}

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

	if( grid.lgGrid )
	{
		/* this is the function's return value */
		chisq = 0.;
	}
	else
	{
		/* this branch optimizing, not grid 
		/ * extract line fluxes 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" );
			cdEXIT(EXIT_FAILURE);
		}

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

		/* first print all warnings */
		cdWarnings(ioQQQ);

		/* print header before any of the individual chi2 values are printed */
		if( optimize.lgOptimize )
			fprintf( ioQQQ, "  ID               Model    Observed       error      chi**2     Type\n" );
		else
			ASSERT( grid.lgGrid );

		/* cycle through the observational values */
		nfound = 0;

		/* first is to optimize relative emission line spectrum */
		if( optimize.xLineInt_Obs.size() > 0 )
		{
			/* set pointers to all optimized lines if first call */
			if( optimize.ipobs.size() == 0 )
			{
				optimize.ipobs.resize( optimize.xLineInt_Obs.size() );
				bool lgHIT = true;
				for( i=0; i < long(optimize.xLineInt_Obs.size()); i++ )
				{
					chi2_type temp1, temp2;
					cap4( chFind, optimize.chLineLabel[i].c_str() );

					/* >> chng 06 may 04, use cdLine instead of ad hoc treatment.
					 * no need to complain, cdLine will do it automatically.  */
					/* this is an intensity, get the line, returns false if could not find it */
					j = cdLine( chFind, optimize.wavelength[i], &temp1, &temp2 );
					if( j <= 0 )
					{
						fprintf( ioQQQ, "\n" );
						lgHIT = false;
					}
					else
					{
						optimize.ipobs[i] = j;
					}
				}

				/* we did not find the line */
				if( !lgHIT )
				{
					fprintf( ioQQQ, "\n\n Optimizer could not find one or more lines.\n" );
					fprintf( ioQQQ, " Sorry.\n");
					cdEXIT(EXIT_FAILURE);
				}
			}

			for( i=0; i < 10; i++ )
			{
				optimize.SavGenericData[i] = 0.;
			}

			for( i=0; i < long(optimize.xLineInt_Obs.size()); i++ )
			{
				/* and find corresponding model value by straight search */
				nfound += 1;
				scld = (chi2_type)LineSv[optimize.ipobs[i]].SumLine[optimize.nEmergent[i]]/
					(chi2_type)snorm*LineSave.ScaleNormLine;
				chi1 = chi2_func(scld,(chi2_type)optimize.xLineInt_Obs[i],
						 (chi2_type)optimize.xLineInt_error[i]);
				cat = 0;
				nobs_cat[cat]++;
				chi2_cat[cat] += chi1;

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

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

				fprintf( ioQQQ, "%12.5f%12.5f%12.5f%12.2e Relative intensity", 
				  scld, 
				  optimize.xLineInt_Obs[i], 
				  optimize.xLineInt_error[i], 
				  chi1 );

				fprintf( ioQQQ, "\n" );

				if( i<10 )
				{
					optimize.SavGenericData[i] = chi1;
				}
			}
		}

		/* this is to optimize a mean temperature */
		for( i=0; i < long(optimize.temp_obs.size()); i++ )
		{
			if( cdTemp( optimize.chTempLab[i].c_str(), optimize.ionTemp[i],
				    &temp_theory, optimize.chTempWeight[i].c_str()) )
			{
				/* did not find column density */
				fprintf(ioQQQ," optimizer did not find column density %s %li \n",
					optimize.chTempLab[i].c_str(),optimize.ionTemp[i] );
				cdEXIT(EXIT_FAILURE);
			}
			nfound += 1;
			chi1 = chi2_func(temp_theory,(chi2_type)optimize.temp_obs[i],
					 (chi2_type)optimize.temp_error[i]);
			cat = 3;
			nobs_cat[cat]++;
			chi2_cat[cat] += chi1;

			fprintf( ioQQQ, " %4.4s%7ld%12.4e%12.4e%12.5f%12.2e Temperature\n",
				 optimize.chTempLab[i].c_str(), optimize.ionTemp[i], temp_theory,
				 optimize.temp_obs[i], optimize.temp_error[i], chi1 );
		}

		/* option to optimize column densities */
		for( i=0; i < long(optimize.ColDen_Obs.size()); i++ )
		{
			if( cdColm(optimize.chColDen_label[i].c_str(),optimize.ion_ColDen[i], &theocl) )
			{
				/* did not find column density */
				fprintf(ioQQQ," optimizer did not find column density %s %li \n",
					optimize.chColDen_label[i].c_str(), optimize.ion_ColDen[i] );
				cdEXIT(EXIT_FAILURE);
			}
			nfound++;
			chi1 = chi2_func(theocl,(chi2_type)optimize.ColDen_Obs[i],
					 (chi2_type)optimize.ColDen_error[i]);
			cat = 1;
			nobs_cat[cat]++;
			chi2_cat[cat] += chi1;

			fprintf( ioQQQ, " %4.4s%7ld%12.4e%12.4e%12.5f%12.2e Column density\n", 
				 optimize.chColDen_label[i].c_str(), optimize.ion_ColDen[i], theocl, 
				 optimize.ColDen_Obs[i], optimize.ColDen_error[i], chi1 );
		}

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

			fprintf( ioQQQ, " %4.4s ", 
				 LineSv[LineSave.ipNormWavL].chALab );

			prt_wl( ioQQQ, LineSv[LineSave.ipNormWavL].wavelength );

			fprintf( ioQQQ, "%12.5f%12.5f%12.5f%12.2e Line intensity\n", 
			  predin, 
			  optimize.optint, 
			  optimize.optier, 
			  chi1 );
		}

		/* option to optimize the absolute continuum flux */
		for( size_t k=0; k < optimize.ContIndex.size(); k++ )
		{
			nfound++;
			// there are 4 entries for each wavelength: nFnu, nInu, InwT, InwC
			long ind = t_PredCont::Inst().offset() + 4*optimize.ContIndex[k];
			chi2_type nFnu_model = 0.;
			if( LineSv[ind].SumLine[0] > SMALLFLOAT )
			{
				nFnu_model = chi2_type( log10(LineSv[ind].SumLine[0]) + radius.Conv2PrtInten );
				nFnu_model = pow(10.,nFnu_model);
			}
			Flux F_model( optimize.ContNFnu[k].E(), nFnu_model );

			chi1 = chi2_func(nFnu_model,optimize.ContNFnu[k].get("erg/s/cm2"),optimize.ContNFnuErr[k]);
			const char* catstr;
			// treat radio continuum flux as absolute flux so that it can be used
			// as a more accurate replacement of the normalization line intensity
			if( optimize.ContEner[k].mm() <= 1. )
			{
				cat = 5;
				catstr = "Photometry";
			}
			else
			{
				cat = 2;
				catstr = "Radio intensity";
			}
			nobs_cat[cat]++;
			chi2_cat[cat] += chi1;

			fprintf( ioQQQ, " %4.4s ", LineSv[ind].chALab);
			prt_wl( ioQQQ,LineSv[ind].wavelength);
			string unit = optimize.ContNFnu[k].uu();
			fprintf( ioQQQ, "%12.4g%12.4g%12.5f%12.2e %s [%s]\n", 
				 F_model.get(unit),
				 optimize.ContNFnu[k].get(unit),
				 optimize.ContNFnuErr[k], chi1,
				 catstr, unit.c_str() );
		}

		/* option to optimize angular diamater */
		if( optimize.lgOptDiam )
		{
			nfound++;
			chi2_type diam_model;
			// get diameter in cm
			if( rfield.lgUSphON )
				diam_model = 2.*rfield.rstrom; // ionization bounded -> use Stroemgren radius
			else
				diam_model = 2.*radius.Radius; // density bounded -> use outer radius
			// now convert to arcsec if necessary
			if( !optimize.lgDiamInCM && radius.distance > 0. )
				diam_model *= AS1RAD/radius.distance;

			chi1 = chi2_func(diam_model,optimize.optDiam,optimize.optDiamErr);
			cat = 4;
			nobs_cat[cat]++;
			chi2_cat[cat] += chi1;

			fprintf( ioQQQ, "            %12.4g%12.4g%12.5f%12.2e Angular diameter\n", 
				 diam_model, optimize.optDiam, optimize.optDiamErr, chi1 );
		}

		/* do not have to have line matches if doing grid. */
		if( nfound <= 0 && !grid.lgGrid )
		{
			fprintf( ioQQQ, " WARNING; no line matches found\n" );
			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] );
			}
		}
		if( nfound )
		{
			fprintf( ioQQQ, "\n Iteration%4ld Chisq=%13.5e\n", optimize.nOptimiz, chisq );
		}
	}

	/* increment nOptimiz, the grid / optimizer counter */
	++optimize.nOptimiz;

	/* only print this if output has been turned on */
	if( called.lgTalk )
	{
		fprintf( ioQQQ, "\n" );
		for( i=0; i < optimize.nvary; 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");
				cdEXIT(EXIT_FAILURE);
			}

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

	return min(chisq,BIG_CHI2);
}

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

	DEBUG_ENTRY( "chi2_func()" );

	/* 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" );
		cdEXIT(EXIT_FAILURE);
	}

	if( yerr > 0. )
	{
		if( ymodl > 0. )
		{
			temp = pow2((ymodl-ymeas)/(min(ymodl,ymeas)*yerr));
			chi2_func_v = min(temp,BIG_CHI2);
		}
		else
			chi2_func_v = BIG_CHI2;
	}
	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 = min(temp,BIG_CHI2);
		}
		else
			chi2_func_v = 0.;
	}
	else
	{
		fprintf( ioQQQ, "chi2_func: relative error is zero, this should not happen\n" );
		cdEXIT(EXIT_FAILURE);
	}
	return chi2_func_v;
}