/*DoOptimize main driver for optimization runs*/
#include <float.h>
#include <string.h>
#include "cddefines.h"
#define	NPLXMX	200
#include "varypar.h"
#include "limfal.h"
#include "input.h"
#include "faint.h"
#include "called.h"
#include "date.h"
#include "trace.h"
#include "amoeba.h"
#include "phymir.h"
#include "func.h"
#include "powell.h"
#include "dooptimize.h"
#include "parse.h"
#include "subplx.h"

/* returns 0 if things went ok, 1 for disaster */
int DoOptimize(void)
{
	char chLine[81], 
	  chNote[8];
	long int i, 
	  iflag, 
	  ii, 
	  iteram, 
	  ival, 
	  iworke[NPLXMX], 
	  j, 
	  mode, 
	  need, 
	  nfe, 
	  nint, 
	  np;
	float fret, 
	  fx, 
	  param[LIMPAR], 
	  ptem[LIMPAR], 
	  delta[LIMPAR], 
	  toler,
	  worke[NPLXMX], 
	  y[21], 
	  ymn;
	double _e0[210];
	float (*const p)[21] = (float(*)[21])_e0;
	float (*const xi)[20] = (float(*)[20])_e0;

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

	/* main driver for optimization runs
	 * Drives cloudy to optimize variables;*/

	/* code originally written by R.F. Carswell, IOA Cambridge */

	/* equivalence the POWELL & AMOEBA workspaces */
	/* Stuff for Subplex */

	/* variables with optimizer */
	for( i=0; i < LIMPAR; i++ )
	{
		VaryPar.OptIncrm[i] = 0.;
		VaryPar.varang[i][0] = -FLT_MAX;
		VaryPar.varang[i][1] = FLT_MAX;
	}

	/* this moved here from zerologic */
	VaryPar.nIterOptim = 20;
	VaryPar.ioOptim = NULL;
	VaryPar.OptGlobalErr = 0.10f;
	VaryPar.nlobs = 0;
	VaryPar.ncobs = 0;
	strcpy( chVaryPar.chOptRtn, "SUBP" );
	VaryPar.lgOptLin = FALSE;
	VaryPar.lgOptLum = FALSE;
	VaryPar.lgOptCol = FALSE;
	VaryPar.lgTrOpt = FALSE;
	VaryPar.lgOptimFlow = FALSE;
	VaryPar.optint = 0.;
	VaryPar.optier = 0.;
#	ifdef __unix
	VaryPar.lgParallel = TRUE;
#	else
	VaryPar.lgParallel = FALSE;
#	endif
	VaryPar.lgOptCont = FALSE;

	called.lgTalk = FALSE;
	/* this flag is needed to turn print on to ahave effect */
	called.lgTalkIsOK = FALSE;

	/* necessary to do this to keep all lines in */
	faint.lgFaintOn = FALSE;
	limfal.LimFail = 1000;

	/* call READR the first time to scan off all variable options */
	ParseCommands();

	VaryPar.nvary = VaryPar.nparm;

	if( VaryPar.lgOptLum )
	{
		nint = 1;
	}
	else
	{
		nint = 0;
	}

	/* check that more than 1 observed intensiteis or column densites were entereed */
	if( (VaryPar.nlobs + nint + VaryPar.ncobs) < 1 )
	{
		fprintf( ioQQQ, " The input stream has vary commands, but\n" );
		fprintf( ioQQQ, " no observed quantities were entered.  Whats up?\n" );
		fprintf( ioQQQ, " Use the NO VARY command to input vary options but not try to perform this.\n" );
		puts( "[Stop in DoOptimize]" );
		exit(1);
	}

	/* check that the total number of parameters to vary is greater than 1 */
	if( VaryPar.nvary < 1 )
	{
		fprintf( ioQQQ, " No parameters to vary were entered. Whats up?\n" );
		puts( "[Stop in DoOptimize]" );
		exit(1);
	}

	/* lgTrOptm set with trace optimize command */
	if( trace.lgTrOptm )
	{
		for( i=0; i < VaryPar.nvary; i++ )
		{
			/*print the command format as debugging aid */
			fprintf( ioQQQ, "%s\n", chVaryPar.chVarFmt[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( chLine , chVaryPar.chVarFmt[i], VaryPar.vparm[0][i] );
			}

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

			else if( VaryPar.nvarxt[i] == 3 )
			{
				/* case with 1 parameter */
				sprintf( chLine , 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 me\n");
				puts( "[Stop in DoOptimize]" );
				exit(1);
			}

			/* print the resulting command line*/
			fprintf( ioQQQ, "%s\n", chLine );
		}
	}

	/* option to change default increments; if zero then leave as is */
	for( i=0; i < LIMPAR; i++ )
	{
		if( VaryPar.OptIncrm[i] != 0. )
		{
			VaryPar.vincr[i] = VaryPar.OptIncrm[i];
		}
	}

	/* say who we are */
	fprintf( ioQQQ, "                                                      Optimization  Driver\n" );
	fprintf( ioQQQ, "                                                         Cloudy %7.7s\n\n",date.chVersion);
	fprintf( ioQQQ, "                       **************************************%7.7s**************************************\n", date.chDate);
	fprintf( ioQQQ, "                       *                                                                                 *\n");


	/* now echo initial input quantities with flag for vary */
	/* first loop steps over all command lines entered */
	for( i=0; i <= input.nSave; i++ )
	{
		/* put space to start line, overwrite if vary found */
		strcpy( chNote, "       " );
		/* loop over all vary commands, see if this is one */
		for( j=0; j < VaryPar.nvary; j++ )
		{
			/* input.nSave is on C array counting, rest are on fortran */
			if( i == VaryPar.nvfpnt[j] )
			{
				/* this is a vary command, put keyword at start */
				strcpy( chNote, "VARY>>>" );
			}
		}

		fprintf( ioQQQ, "               %7.7s * ",chNote);
		j = 0;
		/* print actual command */
		while( input.chCardSav[i][j] !='\0' )
		{
			fprintf( ioQQQ, "%c", input.chCardSav[i][j] );
			++j;
		}
		/* flush out rest of space then * */
		while( j<80 )
		{
			fprintf( ioQQQ, "%c", ' ' );
			++j;
		}
		/* end the line */
		fprintf( ioQQQ, "*\n" );
	}
	fprintf( ioQQQ, "                       *                                                                                 *\n                       ***********************************************************************************\n\n\n" );

	/* option to trace logical flow within this sub */
	if( VaryPar.lgOptimFlow )
	{
		for( j=0; j < VaryPar.nvary; j++ )
		{
			i = VaryPar.nvfpnt[j];
			fprintf( ioQQQ, " trace:%80.80s\n", input.chCardSav[i]);
			fprintf( ioQQQ, "%80.80s\n", chVaryPar.chVarFmt[j]);
			fprintf( ioQQQ, " number of variables on line:%4ld\n", 
			  VaryPar.nvarxt[j] );
			fprintf( ioQQQ, " Values:" );
			for( ii=1; ii <= VaryPar.nvarxt[j]; ii++ )
			{
				fprintf( ioQQQ, "%10.2e", VaryPar.vparm[ii-1][j] );
			}
			fprintf( ioQQQ, "\n" );
		}
	}

	fprintf( ioQQQ, " Up to%5ld iterations will be performed,\n", 
	  VaryPar.nIterOptim );
	fprintf( ioQQQ, " and the final version of the input file will be written to the file %s\n", 
	  chOptimFileName );
	if( strcmp(chVaryPar.chOptRtn,"AMOE") == 0 )
	{
		fprintf( ioQQQ, " The amoeba method will be used.\n" );
	}

	else if( strcmp(chVaryPar.chOptRtn,"PHYM") == 0 )
	{
		fprintf( ioQQQ, " The phymir method will be used" );
		if( VaryPar.lgParallel ) {
			fprintf( ioQQQ, " in parallel mode.\n The maximum no. of CPU's to be used is %1d.\n",VaryPar.maxCPU );
		}
		else {
			fprintf( ioQQQ, " in sequential mode.\n" );
		}
	}

	else if( strcmp(chVaryPar.chOptRtn,"POWE") == 0 )
	{
		fprintf( ioQQQ, " Powells method will be used.\n" );
	}

	else if( strcmp(chVaryPar.chOptRtn,"SUBP") == 0 )
	{
		fprintf( ioQQQ, " Subplex method will be used.\n" );
	}

	else
	{
		fprintf( ioQQQ, " I dont understand what method to use.\n" );
		fprintf( ioQQQ, " Sorry.\n" );
		puts( "[Stop in DoOptimize]" );
		exit(1);
	}

	fprintf( ioQQQ, "\n%4ld parameters will be varied.  The first lines, and the increments are:\n", 
	  VaryPar.nvary );

	for( i=0; i < VaryPar.nvary; i++ )
	{
		VaryPar.varmax[i] = -FLT_MAX;
		VaryPar.varmin[i] = FLT_MAX;
		/*  write formatted to output using the format held in chVarFmt(np) */

		/* 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( chLine , chVaryPar.chVarFmt[i], VaryPar.vparm[0][i] );
		}

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

		else if( VaryPar.nvarxt[i] == 3 )
		{
			/* case with 1 parameter */
			sprintf( chLine , 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 me2\n");
			puts( "[Stop in DoOptimize]" );
			exit(1);
		}

		fprintf( ioQQQ, "\n %s\n", chLine );
		fprintf( ioQQQ, " Initial increment is%6.3f, the limits are%10.2e to %10.2e\n", 
		  VaryPar.vincr[i], VaryPar.varang[i][0], VaryPar.varang[i][1] );
	}

	/* this will be number of times optimizer calls cloudy */
	VaryPar.nOptimiz = 0;
	toler = (float)log10(1. + VaryPar.OptGlobalErr);

	if( strcmp(chVaryPar.chOptRtn,"POWE") == 0 )
	{
		/*  Powell's method */
		for( j=0; j < VaryPar.nvary; j++ )
		{
			ptem[j] = VaryPar.vparm[0][j];
			for( i=0; i < VaryPar.nvary; i++ )
			{
				xi[j][i] = 0.0;
			}
			xi[j][j] = VaryPar.vincr[j];
		}
		Powell(ptem,(float*)xi,&VaryPar.nvary,20,toler,&iteram,&fret);
		for( j=0; j < VaryPar.nvary; j++ )
		{
			VaryPar.vparm[0][j] = ptem[j];
		}

	}

	else if( strcmp(chVaryPar.chOptRtn,"PHYM") == 0 )
	{
		/*  Phymir method */
		for( j=0; j < VaryPar.nvary; j++ )
		{
			ptem[j] = VaryPar.vparm[0][j];
			delta[j] = VaryPar.vincr[j];
		}
		phymir(ptem,delta,VaryPar.nvary,&fret,toler);
		for( j=0; j < VaryPar.nvary; j++ )
		{
			VaryPar.vparm[0][j] = ptem[j];
		}

	}

	else if( strcmp(chVaryPar.chOptRtn,"SUBP") == 0 )
	{
		fprintf( ioQQQ, " Begin optimization with SUBPLEX\n" );
		need = 2*VaryPar.nvary + VaryPar.nvary*(VaryPar.nvary + 4) + 
		  1;
		if( need > NPLXMX )
		{
			fprintf( ioQQQ, " Increase size of NPLXMX in parameter statements to handle this many variables.\n" );
			fprintf( ioQQQ, " I need at least %5ld\n", need );
			puts( "[Stop in DoOptimize]" );
			exit(1);
		}
		for( j=0; j < VaryPar.nvary; j++ )
		{
			ptem[j] = VaryPar.vparm[0][j];
		}

		/* The subrouting SUBPLX input into cloudy 8/4/94.
		 * The program itself is very well commented.
		 * The mode must set to 0 for the default values.
		 * The switch iflag tells if the program terminated normally.   */
		mode = 0;

		/*  >>chng 97 dec 8, remove first arg, func, since not used in routines */
		subplx(VaryPar.nvary,toler,VaryPar.nIterOptim,
		  mode,VaryPar.vincr,ptem,&fx,&nfe,worke,iworke,&iflag);

		if( iflag == -1 )
		{
			fprintf( ioQQQ, " SUBPLEX exceeding maximum iterations.\n This can be reset with the OPTIMZE ITERATIONS command.\n" );
		}

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

		if( VaryPar.lgOptimFlow )
		{
			fprintf( ioQQQ, " trace return subplx:\n" );
			for( j=0; j < VaryPar.nvary; j++ )
			{
				fprintf( ioQQQ, " Values:" );
				for( ii=1; ii <= VaryPar.nvarxt[j]; ii++ )
				{
					fprintf( ioQQQ, "%10.2e", VaryPar.vparm[ii-1][j] );
				}
				fprintf( ioQQQ, "\n" );
			}
		}
	}

	else if( strcmp(chVaryPar.chOptRtn,"AMOE") == 0 )
	{
		/* AMOEBA
		 * set up the two-dimensional array of values for the parameters
		 * required by AMOEBA, NVARY+1 of them. */
		for( j=0; j < VaryPar.nvary; j++ )
		{
			for( i=0; i < (VaryPar.nvary + 1); i++ )
			{
				p[j][i] = VaryPar.vparm[0][j];
			}
		}

		/*  set the single values in each which are different from the
		 *  first guess */
		for( i=1; i < (VaryPar.nvary + 1); i++ )
		{
			j = i - 1;
			p[j][i] = VaryPar.vparm[0][j] + VaryPar.vincr[j];
		}

		/*  initialize y, the value to be minimized */
		for( i=0; i < (VaryPar.nvary + 1); i++ )
		{
			for( j=0; j < VaryPar.nvary; j++ )
			{
				VaryPar.vparm[0][j] = p[j][i];
				param[j] = p[j][i];
			}
			/* FUNC calls CLOUDY with the variable parameters in VPARM,
			 * returning the equivalent of chi**2 as FUNC */
			y[i] = (float)func(param);
		}

		/* Amoeba minimizes the function func(param)
		 * where param is the array of parameters
		 * dimensional array of parameters for the model. It is described
		 * on pages 402 - 406 of the book. [Press et al., Numerical
		 * Recipes (Fortran version), Cambridge U.P., 1992].  the current
		 * version has been totally recoded in c and IS NOT their version,
		 * although it is closely based on their logic */
		Amoeba((float*)p,y,21,VaryPar.nvary,toler,func,&iteram);
		/* print summary of results */
		ival = 1;
		ymn = 1.0e34f;
		/* keep pointer to minimum value */
		for( i=0; i < (VaryPar.nvary + 1); i++ )
		{
			if( y[i] < ymn )
			{
				ymn = y[i];
				ival = i;
			}
		}

		/* end set variable values to those which gave the minimum Y */
		for( j=0; j < VaryPar.nvary; j++ )
		{
			VaryPar.vparm[0][j] = p[j][ival];
		}
	}

	fprintf( ioQQQ, " **************************************************\n" );
	fprintf( ioQQQ, " **************************************************\n" );
	fprintf( ioQQQ, " **************************************************\n" );
	fprintf( ioQQQ, "\n Cloudy was called %4ld times.\n\n", VaryPar.nOptimiz );

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


		fprintf( ioQQQ, " Optimal command: %s\n", input.chCardSav[np]);
		fprintf( ioQQQ, "  Smallest value:%10.2e Largest value:%10.2e Allowed range %10.2e to %10.2e\n", 
		  VaryPar.varmin[i], VaryPar.varmax[i], VaryPar.varang[i][0], 
		  VaryPar.varang[i][1] );
	}

	called.lgTalk = TRUE;
	called.lgTalkIsOK = TRUE;
	faint.lgFaintOn = TRUE;

	/* though a page eject */
	fprintf( ioQQQ, "\f" );

	/* punch optimal parameters on unit ioOptim */
	if( VaryPar.ioOptim == NULL )
	{
		/* open default file name, optimal.in */
		VaryPar.ioOptim = fopen( chOptimFileName , "w" );
		if( VaryPar.ioOptim == NULL )
		{
			fprintf( ioQQQ," could not open optimal.in\n");
			exit(1);
		}
	}

	for( i=0; i <= input.nSave; i++ )
	{
		fprintf( VaryPar.ioOptim, "%s\n", input.chCardSav[i]);
	}
	fclose( VaryPar.ioOptim );

	/* recalculate values in cloudy for the best fit, and print out
	 * all the information */
	fret = (float)func(param);


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

	if( called.lgBusted )
	{
		/* busted set means there were serious problems somewhere */
		return 1;
	}
	else
	{
		/* return 0 if everything is ok */
		return 0;
	}
}