/* 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 */
/*optimize_do main driver for optimization runs*/
#include "cddefines.h"
#define	NPLXMX	(LIMPAR*(LIMPAR+6)+1)
#include "input.h"
#include "called.h"
#include "prt.h"
#include "punch.h"
#include "optimize.h"

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

#	ifdef DEBUG_FUN
	fputs( "<+>optimize_do()\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 */

	toler = (float)log10(1. + optimize.OptGlobalErr);

	if( strcmp(optimize.chOptRtn,"POWE") == 0 )
	{
		float xi[LIMPAR][LIMPAR];

		/*  optimize_powell's method */
		for( j=0; j < optimize.nvary; j++ )
		{
			ptem[j] = optimize.vparm[0][j];
			for( i=0; i < optimize.nvary; i++ )
			{
				xi[j][i] = 0.0;
			}
			xi[j][j] = optimize.vincr[j];
		}
		optimize_powell(ptem,(float*)xi,&optimize.nvary,LIMPAR,toler,&iteram,&fret);
		for( j=0; j < optimize.nvary; j++ )
		{
			optimize.vparm[0][j] = ptem[j];
		}

	}

	else if( strcmp(optimize.chOptRtn,"PHYM") == 0 )
	{
		/*  Phymir method */
		for( j=0; j < optimize.nvary; j++ )
		{
			ptem[j] = optimize.vparm[0][j];
			delta[j] = optimize.vincr[j];
		}
		/* >>chng 06 jan 02, fix uninitialized var problem detected by valgrind/purify, PvH */
		for( j=optimize.nvary; j < LIMPAR; j++ )
		{
			ptem[j] = -FLT_MAX;
			delta[j] = -FLT_MAX;
		}
		optimize_phymir(ptem,delta,optimize.nvary,&fret,toler);
		for( j=0; j < optimize.nvary; j++ )
		{
			optimize.vparm[0][j] = ptem[j];
		}

	}

	else if( strcmp(optimize.chOptRtn,"SUBP") == 0 )
	{
		fprintf( ioQQQ, " Begin optimization with SUBPLEX\n" );
		need = 2*optimize.nvary + optimize.nvary*(optimize.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 optimize_do]" );
			cdEXIT(EXIT_FAILURE);
		}
		for( j=0; j < optimize.nvary; j++ )
		{
			ptem[j] = optimize.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 08, remove first arg, optimize_func, since not used in routines */
		optimize_subplex(
			/* the number of parameters to vary */
			optimize.nvary,
			/* the relative error, single number */
			toler,
			/* maximum number of function evaluations before giving up */
			optimize.nIterOptim,
			/* mode of operation, we simply set to zero */
			mode,
			/* the initial changes in the guessed best coefficients, typically 0.2 to 1  */
			optimize.vincr,
			/* a vector of nvary initial parameters that are the starting guesses for the parameters */
			ptem,
			/* a float, this is simply ignored */
			&fx,
			/* another parameter that is simply ignored, a long int */
			&nfe,
			/* a float that is NPLXMX long, used for working space by the routine */
			/* an array that is 20*26 + 1 elements long, used for working space */
			worke,
			/* a long int that is NPLXMX long, used for working space by the routine */
			/* an array that is 20*26 + 1 elements long, used for working space */
			iworke,
			/* a long int - says what happened, if -1 then exceeded nIterOptim iteration */
			&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 < optimize.nvary; j++ )
		{
			optimize.vparm[0][j] = ptem[j];
		}

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

	else if( strcmp(optimize.chOptRtn,"AMOE") == 0 )
	{
		float p[LIMPAR][LIMPAR+1],y[LIMPAR+1];

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

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

		/*  initialize y, the value to be minimized */
		for( i=0; i < optimize.nvary+1; i++ )
		{
			for( j=0; j < optimize.nvary; j++ )
			{
				optimize.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)optimize_func(param);
		}

		/* optimize_amoeba minimizes the function optimize_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 */
		optimize_amoeba((float*)p,y,LIMPAR+1,optimize.nvary,toler,optimize_func,&iteram);
		/* print summary of results */
		ival = 1;
		ymn = FLT_MAX;
		/* keep pointer to minimum value */
		for( i=0; i < optimize.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 < optimize.nvary; j++ )
		{
			optimize.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", optimize.nOptimiz );

	for( i=0; i < optimize.nvary; i++ )
	{
		optimize.vparm[0][i] = (float)MIN2(optimize.vparm[0][i],optimize.varang[i][1]);
		optimize.vparm[0][i] = (float)MAX2(optimize.vparm[0][i],optimize.varang[i][0]);
		param[i] = optimize.vparm[0][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 me3\n");
			puts( "[Stop in optimize_do]" );
			cdEXIT(EXIT_FAILURE);
		}


		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", 
		  optimize.varmin[i], optimize.varmax[i], optimize.varang[i][0], 
		  optimize.varang[i][1] );
	}

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

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

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

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

	/* flag that punch files should be opened this time */
	punch.lgOpenUnits = TRUE;

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


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

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

