/* 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 */
#include "cddefines.h"
#include "conv.h"
#include "input.h"
#include "called.h"
#include "version.h"
#include "prt.h"
#include "trace.h"
#include "grains.h"
#include "parse.h"
#include "punch.h"
#include "optimize.h"
/* #include "cddrive.h" */
#include "grid.h"
/* #include "rfield.h" */
/* #include "cloudy.h"
#include "lines.h" */

/* called by cdDrive, this returns 0 if things went ok, 1 for disaster */
int grid_do(void)
{
	char chLine[INPUT_LINE_LENGTH], 
	  chNote[8];
	long int i, 
	  ii, 
	  j, 
	  nint;
	float ptem[LIMPAR], 
	  delta[LIMPAR], 
	  toler;
/* 	double _e0[210]; */
/* 	float (*const p)[21] = (float(*)[21])_e0; */
/* 	float (*const xi)[20] = (float(*)[20])_e0; */

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

	/* main driver for optimization runs
	 * Drives cloudy to grid 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++ )
	{
		optimize.OptIncrm[i] = 0.;
		optimize.varang[i][0] = -FLT_MAX;
		optimize.varang[i][1] = FLT_MAX;
		/* this should be overwritten by format of vary line */
		strcpy( optimize.chVarFmt[i], "error - no optimizer line image was set" );
	}

	/* this moved here from zerologic */
	optimize.nIterOptim = 20;
	optimize.ioOptim = NULL;
	optimize.OptGlobalErr = 0.10f;
	optimize.nlobs = 0;
	optimize.nTempObs = 0;
	optimize.ncobs = 0;
	optimize.nRangeSet = 0;
	strcpy( optimize.chOptRtn, "SUBP" );

	/* flags says what is to be matched */
	optimize.lgOptLin = FALSE;
	optimize.lgOptLum = FALSE;
	optimize.lgOptCol = FALSE;
	optimize.lgOptTemp = FALSE;

	/* trace flag for optim process */
	optimize.lgTrOpt = FALSE;

	optimize.lgOptimFlow = FALSE;
	optimize.optint = 0.;
	optimize.optier = 0.;
#	ifdef __unix
	optimize.lgParallel = TRUE;
#	else
	optimize.lgParallel = FALSE;
#	endif
	optimize.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 */
	prt.lgFaintOn = FALSE;
	conv.LimFail = 1000;

	/* flag that no punch files should be opened until the final model */
	punch.lgOpenUnits = FALSE;

	/* >>chng 01 sep 11, removed call to OldStyleGrainBlockDataInit, now done in GrainZero() PvH */
	/* OldStyleGrainBlockDataInit(); */

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

	/* >>chng 00 aug 09, return memory allocated for grains, they are not used, PvH */
	ReturnGrainBins();

	optimize.nvary = optimize.nparm;

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

	/* check that more than 1 observed intensiteis or column densites were entereed */
	if( ((optimize.nlobs + nint + optimize.nTempObs + optimize.ncobs) < 1 ) && ( grid.lgGrid == FALSE ) )
	{
		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 grid_do]" );
		cdEXIT(EXIT_FAILURE);
	}

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

	if( ( strcmp(optimize.chOptRtn,"XSPE") == 0 ) && ( optimize.nRangeSet != optimize.nvary ) )
	{
		fprintf( ioQQQ, " For XSpec option, every parameter to be varied must have a range specified,\n" );
		fprintf( ioQQQ, " but %ld parameter(s) did not.",  optimize.nvary - optimize.nRangeSet );
		puts( "[Stop in grid_do]" );
		cdEXIT(EXIT_FAILURE);
	}

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

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

			else if( optimize.nvarxt[i] == 3 )
			{
				/* case with 3 parameter */
				sprintf( chLine , 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( chLine , 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( chLine , 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 me1\n");
				puts( "[Stop in grid_do]" );
				cdEXIT(EXIT_FAILURE);
			}

			/* 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( optimize.OptIncrm[i] != 0. )
		{
			optimize.vincr[i] = optimize.OptIncrm[i];
		}
	}

	/* say who we are */
	if( strcmp(optimize.chOptRtn,"XSPE") == 0 )
	{
		fprintf( ioQQQ, "                                                           Grid  Driver\n" );
	}
	else
	{
		fprintf( ioQQQ, "                                                      Optimization  Driver\n" );
	}
	fprintf( ioQQQ, "                                                        Cloudy %s\n\n",version.chVersion);
	fprintf( ioQQQ, "                       **************************************%7.7s**************************************\n", version.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 < optimize.nvary; j++ )
		{
			/* input.nSave is on C array counting, rest are on fortran */
			if( i == optimize.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( optimize.lgOptimFlow )
	{
		for( j=0; j < optimize.nvary; j++ )
		{
			i = optimize.nvfpnt[j];
			fprintf( ioQQQ, " trace:%80.80s\n", input.chCardSav[i]);
			fprintf( ioQQQ, "%80.80s\n", optimize.chVarFmt[j]);
			fprintf( ioQQQ, " number of variables on line:%4ld\n", 
			  optimize.nvarxt[j] );
			fprintf( ioQQQ, " Values:" );
			for( ii=1; ii <= optimize.nvarxt[j]; ii++ )
			{
				fprintf( ioQQQ, "%10.2e", optimize.vparm[ii-1][j] );
			}
			fprintf( ioQQQ, "\n" );
		}
	}

	fprintf( ioQQQ, " Up to%5ld iterations will be performed,\n", 
	  optimize.nIterOptim );
	fprintf( ioQQQ, " and the final version of the input file will be written to the file %s\n", 
	  chOptimFileName );

	if( strcmp(optimize.chOptRtn,"AMOE") == 0 )
	{
		fprintf( ioQQQ, " The amoeba method will be used.\n" );
	}

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

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

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

	else if( strcmp(optimize.chOptRtn,"XSPE") == 0 )
	{
		fprintf( ioQQQ, " Producing Xspec output.\n" );
	}

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

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

	for( i=0; i < optimize.nvary; i++ )
	{
		optimize.varmax[i] = -FLT_MAX;
		optimize.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( optimize.nvarxt[i] == 1 )
		{
			/* case with 1 parameter */
			sprintf( chLine , optimize.chVarFmt[i], optimize.vparm[0][i] );
		}

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

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

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

	/* this will be number of times gridr calls cloudy */
	optimize.nOptimiz = 0;
	toler = (float)log10(1. + optimize.OptGlobalErr);
	grid.numParamValues = 3;

	if( strcmp(optimize.chOptRtn,"XSPE") == 0 )
	{
		for( j=0; j < optimize.nvary; j++ )
		{
			/* ptem[j] = optimize.vparm[0][j]; */
			ptem[j] = optimize.varang[j][0]; 
			/* delta[j] = optimize.vincr[j]; */
			delta[j] = ( optimize.varang[j][1] - optimize.varang[j][0] )/
				((float)grid.numParamValues - 1.f);
		}
		gridXspec(ptem,delta,optimize.nvary,grid.numParamValues);
		for( j=0; j < optimize.nvary; j++ )
		{
			optimize.vparm[0][j] = ptem[j];
		}

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

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

		called.lgTalk = TRUE;
		called.lgTalkIsOK = TRUE;
		prt.lgFaintOn = TRUE;
		/* flag that punch files should be opened this time */
		punch.lgOpenUnits = TRUE;
	}
	else
	{
		optimize_do();
	}

#	ifdef DEBUG_FUN
	fputs( " <->grid_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;
	}
}

