/* 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 */
/*ParseOptimize parse the optimize command line */
/*GetOptColDen read observed column densities & errors for optimizer */
/*GetOptLineInt parse observed line intensites for optimization routines */
#include "cddefines.h"
#include "trace.h"
#include "optimize.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 );

void ParseOptimize(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 )
		{
			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 lcolumn 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 )
		{
			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;
			}
			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( lgMatch("SUBP",chCard) )
	{
		/* use subplex to optimize parameters */
		strcpy( optimize.chOptRtn, "SUBP" );

	}
	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
	{
		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 chColmnLab */
		cap4((char*)optimize.chColmnLab[optimize.ncobs] , chCard );

		/* now get the ion stage, this should be 1 for atom, up to element
		 * number plus one */
		i = 5;
		optimize.ionam[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.ionam[optimize.ncobs] < 0 )
		{
			fprintf( ioQQQ, " %s\n", chCard );
			fprintf( ioQQQ, " An ionization stage of%4ld does not make sense.  Sorry.\n", 
			  optimize.ionam[optimize.ncobs] );
			puts( "[Stop in GetOptColDen]" );
			cdEXIT(EXIT_FAILURE);
		}

		optimize.amcol[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.amcerr[optimize.ncobs] = (float)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
		if( optimize.amcerr[optimize.ncobs] <= 0.0 )
		{
			/* this is the relative error allowed */
			optimize.amcerr[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.amcerr[optimize.ncobs] = -optimize.amcerr[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.chColmnLab[i], optimize.ionam[i], optimize.amcol[i], 
			  optimize.amcerr[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.chAmLab[optimize.nlobs], chCard , 4 );
		/* null terminate the label*/
		optimize.chAmLab[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.amint[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.amint[optimize.nlobs] <= 0. )
		{
			fprintf( ioQQQ, " An observed intensity of %.2e is not allowed.  Sorry.\n", 
			  optimize.amint[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.amierr[optimize.nlobs] = (float)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
		/* most often will use the default */
		if( optimize.amierr[optimize.nlobs] <= 0.0 )
		{
			/* this is the relative error allowed */
			optimize.amierr[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.amierr[optimize.nlobs] = -optimize.amierr[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.chAmLab[i] );
			prt_wl( ioQQQ, optimize.wavelength[i] );

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

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

