/* This file is part of Cloudy and is copyright (C)1978-2006 by Gary J. Ferland
 * For conditions of distribution and use see copyright notice in license.txt */
/*ParseOptimize parse the optimize command line */
/*GetOptColDen read observed column densities & errors for optimizer */
/*GetOptLineInt parse observed line intensites for optimization routines */
/*GetOptTemp read observed temperatures & errors for optimizer */
#include "cddefines.h"
#include "trace.h"
#include "optimize.h"
#include "grid.h"
#include "input.h"
#include "prt.h"
#include "parse.h"
#include "lines_service.h"
#ifdef DEFERR
#	undef DEFERR
#endif
#define	DEFERR	0.05f

/*GetOptLineInt parse observed line intensites for optimization routines */
static void GetOptLineInt(char *chCard );

/*GetOptColDen read observed column densities & errors for optimizer */
static void GetOptColDen(char *chCard );

/*GetOptTemp read observed temperatures & errors for optimizer */
static void GetOptTemp(char *chCard );

void ParseOptimize(
	/* command line, which was changed to all caps in main parsing routine */
	char *chCard)
{
	int lgEOL;
	long int i;
#	ifdef __unix
	long int dum;
#	endif

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

	/* this must come first so that contents of filename do not trigger wrong command */
	if( lgMatch("FILE",chCard) )
	{
		/* option to send final set of parameters to an input file 
		 * get name within double quotes, 
		 * and set to blanks in chCard and OrgCard */
		/* file handle set NULL in optimize_do when first set up,
		 * this routine called many times during optimization process,
		 * only want to open file one time */
		if( optimize.ioOptim == NULL )
		{
			/* chCard is all caps at this point.  GetQuote will work with 
			 * original version of command line, which preserves case of
			 * characters.  Also removes string between quotes */
			GetQuote( chOptimFileName , chCard , TRUE );

			/* open the file */
			optimize.ioOptim = fopen( chOptimFileName , "w" ) ;
			if( optimize.ioOptim == NULL )
			{
				fprintf(ioQQQ," error opening file %s\n", chOptimFileName );
				cdEXIT(EXIT_FAILURE);
			}
		}

	}
	else if( lgMatch("AMOE",chCard) )
	{
		/* use optimize_amoeba to optimize parameters */
		strcpy( optimize.chOptRtn, "AMOE" );
	}

	else if( lgMatch("COLU",chCard) )
	{
		/* optimize column density */
		optimize.lgOptCol = TRUE;

		/* read column densities to match */
		GetOptColDen(chCard);
	}

	else if( lgMatch("CONT",chCard) )
	{
		/* set flag saying that optimization should start from continue file */
		optimize.lgOptCont = TRUE;
	}

	else if( lgMatch("INCR",chCard) )
	{
		/* scan off increments for the previously selected paramter */
		if( optimize.nparm > 0 )
		{
			/* also called during optimization process, ignore then */
			i = 5;
			optimize.OptIncrm[optimize.nparm-1] = 
				(float)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
		}
	}

	else if( lgMatch("LUMI",chCard) || lgMatch("INTE",chCard) )
	{
		/* scan off intensity or luminosity of normalization line */
		i = 5;
		optimize.optint = (float)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
		optimize.optier = (float)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
		if( lgEOL )
		{
			optimize.optier = DEFERR;
		}

		/* set flag to say that intentity or luminosity of line set */
		optimize.lgOptLum = TRUE;
	}

	else if( lgMatch("ITER",chCard) )
	{
		/* scan off number of iterations */
		i = 5;
		optimize.nIterOptim = (long)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
	}

	else if( lgMatch("LINE",chCard) )
	{
		/* read lines to match */
		GetOptLineInt(chCard);

		/* set flag saying that something has been set */
		optimize.lgOptLin = TRUE;
	}

	else if( lgMatch("PHYM",chCard) )
	{
		/* use PHYMIR to optimize parameters */
		strcpy( optimize.chOptRtn, "PHYM" );
#		ifdef __unix
		optimize.lgParallel = ! lgMatch("SEQU",chCard);
		if( optimize.lgParallel ) {
			i = 5;
			dum = (long)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
			/* default has already been set in cdinit ! */
			if( ! lgEOL ) {
				optimize.maxCPU = dum;
			}
		}
		else {
			optimize.maxCPU = 1;
		}
#		else
		optimize.lgParallel = FALSE;
		optimize.maxCPU = 1;
#		endif
	}

	else if( lgMatch("POWE",chCard) )
	{
		/* use optimize_powell to optimize parameters */
		strcpy( optimize.chOptRtn, "POWE" );
	}

	else if( lgMatch("RANG",chCard) )
	{
		/* scan off range for the previously selected variable */
		if( optimize.nparm > 0 )
		{
			int firstOneReal = FALSE;

			i = 5;
			optimize.varang[optimize.nparm-1][0] = 
				(float)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
			if( lgEOL )
			{
				optimize.varang[optimize.nparm-1][0] = -1e38f;
			}
			else
			{
				firstOneReal = TRUE;
			}

			optimize.varang[optimize.nparm-1][1] = 
				(float)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
			if( lgEOL )
			{
				optimize.varang[optimize.nparm-1][1] = 1e38f;
			}
			else if( firstOneReal == TRUE )
			{
				/* If we get two numbers on the line, we'll increment the number
				 * of parameters with a real range set, but only if the second is
				 * greater than the first.  */
				if( optimize.varang[optimize.nparm-1][1] > optimize.varang[optimize.nparm-1][0] )
					++optimize.nRangeSet;
			}
		}
	}

	else if( lgMatch("SUBP",chCard) )
	{
		/* use subplex to optimize parameters */
		strcpy( optimize.chOptRtn, "SUBP" );
	}

	/* match a temperature */
	else if( lgMatch("TEMP",chCard) )
	{
		/* read temperatures to match */
		GetOptTemp(chCard);

		/* set flag saying that optimzie temps has been set */
		optimize.lgOptTemp = TRUE;
	}

	else if( lgMatch("TOLE",chCard) )
	{
		/* scan off tolerance of fit, sum of residuals must be smaller than this
		 * default is 0.10 */
		i = 5;
		optimize.OptGlobalErr = (float)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
	}

	else if( lgMatch("TRAC",chCard) )
	{
		if( lgMatch("STAR",chCard) )
		{
			/* trace start iteration number
			 * turn on trace printout starting on nth call to cloudy */
			i = 5;
			optimize.nTrOpt = (long)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
			if( lgEOL )
			{
				fprintf( ioQQQ, " optimize trace start command:\n" );
				fprintf( ioQQQ, " The iteration number must appear.\n" );
				puts( "[Stop in ParseOptimize]" );
				cdEXIT(EXIT_FAILURE);
			}
			optimize.lgTrOpt = TRUE;
		}
		else if( lgMatch("FLOW",chCard) )
		{
			/* trace flow
			 * follow logical flow within code */
			optimize.lgOptimFlow = TRUE;
		}
		else
		{
			fprintf( ioQQQ, " optimize trace flow command:\n" );
			fprintf( ioQQQ, " One of the sub keys START or FLOW must appear.\n" );
			puts( "[Stop in ParseOptimize]" );
			cdEXIT(EXIT_FAILURE);
		}
	}
	
	else if( lgMatch("XSPE",chCard) )
	{
		/* use subplex to optimize parameters */
		strcpy( optimize.chOptRtn, "XSPE" );
		grid.lgGrid = TRUE;
	}

	else
	{
		fprintf( ioQQQ, " Unrecognized keyword, consult HAZY.\n" );
		puts( "[Stop in ParseOptimize]" );
		cdEXIT(EXIT_FAILURE);
	}

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

/*GetOptColDen read observed column densities & errors for optimizer */
static void GetOptColDen(char *chCard )
{
	char chCap[INPUT_LINE_LENGTH];
	int lgEOF, 
	  lgEOL;
	long int i;

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

	/* read observed column densities & errors */
	optimize.ncobs = 0;

	/* get first line */
	input_readarray(chCard,&lgEOF);
	if( lgEOF )
	{
		fprintf( ioQQQ, " Hit EOF while reading column density list; use END to end list.\n" );
		puts( "[Stop in GetOptColDen]" );
		cdEXIT(EXIT_FAILURE);
	}

	/* now copy this line to chCap, then convert to caps */
	strcpy( chCap, chCard );
	caps(chCap);

	while( !lgEOF )
	{
		if( optimize.ncobs > NCOLLM )
		{
			fprintf( ioQQQ, " Too many column densities have been entered; the limit is%4ld.  Increase variable NCOLLM.\n", 
			  NCOLLM );
			puts( "[Stop in GetOptColDen]" );
			cdEXIT(EXIT_FAILURE);
		}

		/* order on line is element label (col 1-4), ionization stage, column density, err */
		/* copy cap'd version of first 4 char of chCard to chColDen_label */
		cap4((char*)optimize.chColDen_label[optimize.ncobs] , chCard );

		/* now get the ion stage, this should be 1 for atom, up to element
		 * number plus one */
		i = 5;
		optimize.ion_ColDen[optimize.ncobs] = (long int)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
		if( lgEOL )
		{
			fprintf( ioQQQ, " %s\n", chCard );
			fprintf( ioQQQ, " The ionization stage MUST appear on this line.  Sorry.\n" );
			puts( "[Stop in GetOptColDen]" );
			cdEXIT(EXIT_FAILURE);
		}

		/* the ion must be 1 or greater unless requesting a special,
		 * like a molecule or excited state population, in which
		 * case ion = 0
		 * can't check on upper limit yet since have not resolved element name */
		if( optimize.ion_ColDen[optimize.ncobs] < 0 )
		{
			fprintf( ioQQQ, " %s\n", chCard );
			fprintf( ioQQQ, " An ionization stage of%4ld does not make sense.  Sorry.\n", 
			  optimize.ion_ColDen[optimize.ncobs] );
			puts( "[Stop in GetOptColDen]" );
			cdEXIT(EXIT_FAILURE);
		}

		optimize.ColDen_Obs[optimize.ncobs] = (float)pow(10.,FFmtRead(chCard,&i,
		  INPUT_LINE_LENGTH,&lgEOL));
		if( lgEOL )
		{
			fprintf( ioQQQ, " %80.80s\n", chCard );
			fprintf( ioQQQ, " An observed column density MUST be entered.  Sorry.\n" );
			puts( "[Stop in GetOptColDen]" );
			cdEXIT(EXIT_FAILURE);
		}

		optimize.chColDen_error[optimize.ncobs] = (float)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
		if( optimize.chColDen_error[optimize.ncobs] <= 0.0 )
		{
			/* this is the relative error allowed */
			optimize.chColDen_error[optimize.ncobs] = (float)DEFERR;
		}

		/* check if number is a limit - if '<' appears on the line then it is */
		if( strchr( chCard , '<' ) != NULL )
		{
			/* value is an upper limit only, use negative error to flag */
			optimize.chColDen_error[optimize.ncobs] = -optimize.chColDen_error[optimize.ncobs];
		}

		input_readarray(chCard,&lgEOF);
		strcpy( chCap, chCard );
		caps(chCap);
		if( lgEOF )
		{
			fprintf( ioQQQ, " Hit EOF while reading column density list; use END to end list.\n" );
			puts( "[Stop in GetOptColDen]" );
			cdEXIT(EXIT_FAILURE);
		}

		if( strncmp( chCap , "END" , 3) == 0 )
		{
			lgEOF = TRUE;
		}

		/* now increment the number of columns we have entered */
		optimize.ncobs += 1;
	}

	if( trace.lgTrace && optimize.lgTrOpt )
	{
		fprintf( ioQQQ, "%4ld colunms were entered, they were;\n", 
		  optimize.ncobs );
		for( i=0; i < optimize.ncobs; i++ )
		{
			fprintf( ioQQQ, " %4.4s ion=%5ld%10.2e%10.2e\n", 
			  optimize.chColDen_label[i], optimize.ion_ColDen[i], optimize.ColDen_Obs[i], 
			  optimize.chColDen_error[i] );
		}
	}

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

/*GetOptLineInt parse observed line intensites for optimization routines */
static void GetOptLineInt(char *chCard )
{
	char chCap[INPUT_LINE_LENGTH];

	int lgEOF, 
	  lgEOL;
	long int i;

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

	/* read observed line fluxes & errors */
	optimize.nlobs = 0;

	input_readarray(chCard,&lgEOF);
	if( lgEOF )
	{
		fprintf( ioQQQ, " Hit EOF while reading line list; use END to end list.\n" );
		puts( "[Stop in GetOptLineInt]" );
		cdEXIT(EXIT_FAILURE);
	}

	strcpy( chCap, chCard );
	caps(chCap);

	while( !lgEOF )
	{
		if( optimize.nlobs >= NOBSLM )
		{
			fprintf( ioQQQ, 
				" Too many lines have been entered; the limit is %ld.  Increase variable NOBSLM.\n", 
			  NOBSLM );
			puts( "[Stop in GetOptLineInt]" );
			cdEXIT(EXIT_FAILURE);
		}

		/* order on line is label (col 1-4), wavelength, flux, error */
		strncpy( optimize.chLineLabel[optimize.nlobs], chCard , 4 );
		/* null terminate the label*/
		optimize.chLineLabel[optimize.nlobs][4] = 0;

		i = 5;
		/* next get the wavelength */
		optimize.wavelength[optimize.nlobs] = (float)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);

		/* now convert wavelength to angstroms */
		/* line was entered, look for possible micron or cm label */
		if( input.chCARDCAPS[i-1] == 'M' )
		{
			/* microns */
			optimize.wavelength[optimize.nlobs] *= 1e4f;
		}
		else if( input.chCARDCAPS[i-1] == 'C' )
		{
			/* microns */
			optimize.wavelength[optimize.nlobs] *= 1e8f;
		}

		/* get the error assocated with 4 significant figures */
		optimize.errorwave[optimize.nlobs] = 
			WavlenErrorGet( optimize.wavelength[optimize.nlobs] );
	
		/* next get the observed intensity */
		optimize.xLineInt_Obs[optimize.nlobs] = (float)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
		if( lgEOL )
		{
			fprintf( ioQQQ, " %s\n", chCard );
			fprintf( ioQQQ, " The wavelength and relative intensity MUST be entered on this line.  Sorry.\n" );
			fprintf( ioQQQ, " The command line is the following:\n %s\n", chCard );
			puts( "[Stop in GetOptLineInt]" );
			cdEXIT(EXIT_FAILURE);
		}

		if( optimize.xLineInt_Obs[optimize.nlobs] <= 0. )
		{
			fprintf( ioQQQ, " An observed intensity of %.2e is not allowed.  Sorry.\n", 
			  optimize.xLineInt_Obs[optimize.nlobs] );
			fprintf( ioQQQ, " The command line is the following:\n %s\n", chCard );
			puts( "[Stop in GetOptLineInt]" );
			cdEXIT(EXIT_FAILURE);
		}

		/* finally the optional error */
		optimize.xLineInt_error[optimize.nlobs] = (float)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
		/* most often will use the default */
		if( optimize.xLineInt_error[optimize.nlobs] <= 0.0 )
		{
			/* this is the relative error allowed */
			optimize.xLineInt_error[optimize.nlobs] = (float)DEFERR;
		}

		/* check if number is a limit - if '<' appears on the line then it is */
		if( strchr( chCard , '<' ) != NULL )
		{
			/* value is an upper limit only, use negative error to flag */
			optimize.xLineInt_error[optimize.nlobs] = -optimize.xLineInt_error[optimize.nlobs];
		}

		/* get next line */
		input_readarray(chCard,&lgEOF);
		if( lgEOF )
		{
			fprintf( ioQQQ, " Hit EOF while reading line list for optimize command; use END to end list.\n" );
			puts( "[Stop in GetOptLineInt]" );
			cdEXIT(EXIT_FAILURE);
		}

		strcpy( chCap, chCard );
		caps(chCap);
		if( strncmp( chCap ,"END" , 3 ) == 0 )
			lgEOF = TRUE;

		/* finally increment the number of observed lines */
		++optimize.nlobs;
	}

	if( trace.lgTrace && trace.lgTrOptm )
	{
		fprintf( ioQQQ, "%4ld lines were entered, they were;\n", 
		  optimize.nlobs );

		for( i=0; i < optimize.nlobs; i++ )
		{
			fprintf( ioQQQ, " %4.4s ", optimize.chLineLabel[i] );
			prt_wl( ioQQQ, optimize.wavelength[i] );

			fprintf( ioQQQ, " %10.2e%10.2e\n", 
				optimize.xLineInt_Obs[i], 
				optimize.xLineInt_error[i] );
		}
	}

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

/*GetOptTemp parse observed line intensites for optimization routines */
static void GetOptTemp(char *chCard )
{
	char chCap[INPUT_LINE_LENGTH];

	int lgEOF, 
	  lgEOL;
	long int i;

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

	/* read observed line fluxes & errors - first set total number of observe temps */
	optimize.nTempObs = 0;

	input_readarray(chCard,&lgEOF);
	if( lgEOF )
	{
		fprintf( ioQQQ, " Hit EOF while reading line list; use END to end list.\n" );
		puts( "[Stop in GetOptTemp]" );
		cdEXIT(EXIT_FAILURE);
	}

	/* make line caps so we can parse it */
	strcpy( chCap, chCard );
	caps(chCap);

	while( !lgEOF )
	{
		if( optimize.nTempObs >= NOBSLM )
		{
			fprintf( ioQQQ, 
				" Too many temperatures have been entered; the limit is %ld.  Increase variable NOBSLM.\n", 
			  NOBSLM );
			puts( "[Stop in GetOptTemp]" );
			cdEXIT(EXIT_FAILURE);
		}

		/* order on line is label (col 1-4), ion, temperature, error */
		strncpy( optimize.chTempLab[optimize.nTempObs], chCard , 4 );
		/* null terminate the label*/
		optimize.chTempLab[optimize.nTempObs][4] = 0;

		i = 5;
		/* next get the ion stage */
		optimize.ionTemp[optimize.nTempObs] = (long)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
	
		/* next get the observed temperature */
		optimize.temp_obs[optimize.nTempObs] = (float)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
		if( lgEOL )
		{
			fprintf( ioQQQ, " %s\n", chCard );
			fprintf( ioQQQ, " The ion stage and temperature MUST be entered on this line.  Sorry.\n" );
			fprintf( ioQQQ, " The command line is the following:\n %s\n", chCard );
			puts( "[Stop in GetOptTemp]" );
			cdEXIT(EXIT_FAILURE);
		}

		/* temperatures less than or equal to 10 are logs */
		if( optimize.temp_obs[optimize.nTempObs] <= 10. )
		{
			optimize.temp_obs[optimize.nTempObs] = (float)pow( 10. , (double)optimize.temp_obs[optimize.nTempObs] );
		}

		/* finally the optional error */
		optimize.temp_error[optimize.nTempObs] = (float)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
		/* most often will use the default */
		if( optimize.temp_error[optimize.nTempObs] <= 0.0 )
		{
			/* this is the relative error allowed */
			optimize.temp_error[optimize.nTempObs] = (float)DEFERR;
		}

		/* check if number is a limit - if '<' appears on the line then it is */
		if( strchr( chCard , '<' ) != NULL )
		{
			/* value is an upper limit only, use negative error to flag */
			optimize.temp_error[optimize.nTempObs] = -optimize.temp_error[optimize.nTempObs];
		}

		/* check for radius or volume for how to weight the mean temp 
		 * this will be the default */
		strcpy( optimize.chTempWeight[optimize.nTempObs] , "radius" );
		/* >>chng 05 dec 29, from chCard to chCap, unlike much of code in this file,
		 * chCard has been read in by this routine and contains the original form of 
		 * the command line.  It was converted to caps and stored in chCap above.
		 * As it was written only VOLUME was matched, would not match volumne.
		 * Bug caught and corrected by Bohdan Melekh */
		/*if( lgMatch( "VOLUME" , chCard ) )*/
		if( lgMatch( "VOLUME" , chCap ) )
		{
			strcpy( optimize.chTempWeight[optimize.nTempObs] , "volume" );
		}

		/* get next line */
		input_readarray(chCard,&lgEOF);
		if( lgEOF )
		{
			fprintf( ioQQQ, " Hit EOF while reading line list for optimize command; use END to end list.\n" );
			puts( "[Stop in GetOptTemp]" );
			cdEXIT(EXIT_FAILURE);
		}

		strcpy( chCap, chCard );
		caps(chCap);
		if( strncmp( chCap ,"END" , 3 ) == 0 )
			lgEOF = TRUE;

		/* finally increment the number of observed lines */
		++optimize.nTempObs;
	}

	if( trace.lgTrace && trace.lgTrOptm )
	{
		fprintf( ioQQQ, "%4ld temperatures were entered, they were;\n", 
		  optimize.nTempObs );

		for( i=0; i < optimize.nTempObs; i++ )
		{
			fprintf( ioQQQ, " %4.4s ", optimize.chTempLab[i] );
			fprintf( ioQQQ, " %li " , optimize.ionTemp[i] );

			fprintf( ioQQQ, " %.2e %.2e\n", 
				optimize.temp_obs[i], 
				optimize.temp_error[i] );
		}
	}

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

#undef	DEFERR

