/* This file is part of Cloudy and is copyright (C)1978-2013 by Gary J. Ferland and
 * others.  For conditions of distribution and use see copyright notice in license.txt */
/*optimize_phymir Peter van Hoof's optimization routine */

#include "cddefines.h"
#include "version.h"
#include "optimize.h"
#if defined(__unix) || defined(__APPLE__)
#include <cstddef>
#include <sys/types.h>
#include <sys/wait.h>
#include <unistd.h>
#else
#define pid_t int
#define fork() TotalInsanityAsStub<pid_t>()
#define wait(X) TotalInsanityAsStub<pid_t>()
#endif
#include "mpi_utilities.h"


/** this is the name of the file that optimize_phymir automatically creates containing
 * information to continue the optimization at a later time.  */
const char* STATEFILE = "continue.pmr";
const char* STATEFILE_BACKUP = "continue.bak";

/* The optimization algorithm Phymir and its subsidiary routines have been
 * written by Peter van Hoof. */

inline void wr_block(const void*,size_t,const char*);
inline void rd_block(void*,size_t,const char*);


template<class X, class Y, int NP, int NSTR>
void phymir_state<X,Y,NP,NSTR>::p_clear1()
{
	DEBUG_ENTRY( "p_clear1()" );

	// first zero out the entire struct so that even the padding gets initialized
	// this prevents complaints about writing uninitialized bytes by valgrind
	memset( this, 0, sizeof(*this) );

	p_xmax = numeric_limits<X>::max();
	p_ymax = numeric_limits<Y>::max() / Y(2.);
	for( int i=0; i < 2*NP+1; ++i )
	{
		for( int j=0; j < NP; ++j )
			p_xp[i][j] = -p_xmax;
		p_yp[i] = -p_ymax;
	}
	for( int i=0; i < NP; ++i )
	{
		p_absmin[i] = -p_xmax;
		p_absmax[i] = p_xmax;
		p_varmin[i] = p_xmax;
		p_varmax[i] = -p_xmax;
		for( int j=0; j < NP; ++j )
			p_a2[i][j] = -p_xmax;
		p_c1[i] = -p_xmax;
		p_c2[i] = -p_xmax;
		p_xc[i] = -p_xmax;
		p_xcold[i] = -p_xmax;
	}
	p_vers = VRSNEW;
	p_toler = -p_xmax;
	p_dmax = X(0.);
	p_dold = X(0.);
	p_ymin = p_ymax;
	p_dim = int32(NP);
	p_sdim = int32(NSTR);
	p_nvar = int32(0);
	p_noptim = int32(0);
	p_maxiter = int32(0);
	p_jmin = int32(-1);
	p_maxcpu = int32(1);
	p_curcpu = int32(0);
	p_mode = PHYMIR_ILL;
	// p_chState, and p_chStr[123] have already been initialized with memset above

	// calculate the size of the struct from the start upto, but not including p_size
	// this section of the struct will be written in binary mode to a state file
	// cast pointers to size_t so that we get the size in bytes for fwrite / fread
	// make p_size an uint32 so that the width of the field is platform independent
	p_size = static_cast<uint32>(reinterpret_cast<size_t>(&p_size) - reinterpret_cast<size_t>(this));
	p_func = NULL;
}

template<class X, class Y, int NP, int NSTR>
void phymir_state<X,Y,NP,NSTR>::p_wr_state(const char *fnam) const
{
	DEBUG_ENTRY( "p_wr_state()" );

	if( cpu.i().lgMaster() && strlen(fnam) > 0 )
	{
		FILE *fdes = open_data( fnam, "wb", AS_LOCAL_ONLY_TRY );
		bool lgErr = ( fdes == NULL );
		lgErr = lgErr || ( fwrite( &p_size, sizeof(uint32), 1, fdes ) != 1 );
		lgErr = lgErr || ( fwrite( this, static_cast<size_t>(p_size), 1, fdes ) != 1 );
		lgErr = lgErr || ( fclose(fdes) != 0 );
		if( lgErr )
		{
			printf( "p_wr_state: error writing file: %s\n", fnam );
			remove(fnam);
		}
	}
	return;
}

template<class X, class Y, int NP, int NSTR>
void phymir_state<X,Y,NP,NSTR>::p_rd_state(const char *fnam)
{
	DEBUG_ENTRY( "p_rd_state()" );

	FILE *fdes = open_data( fnam, "rb", AS_LOCAL_ONLY );
	uint32 wrsize;
	bool lgErr = ( fread( &wrsize, sizeof(uint32), 1, fdes ) != 1 );
	lgErr = lgErr || ( p_size != wrsize );
	lgErr = lgErr || ( fread( this, static_cast<size_t>(p_size), 1, fdes ) != 1 );
	lgErr = lgErr || ( fclose(fdes) != 0 );
	if( lgErr )
	{
		printf( "p_rd_state: error reading file: %s\n", fnam );
		cdEXIT(EXIT_FAILURE);
	}
	return;
}

template<class X, class Y, int NP, int NSTR>
Y phymir_state<X,Y,NP,NSTR>::p_execute_job(const X x[],
					   int jj,
					   int runNr)
{
	DEBUG_ENTRY( "p_execute_job()" );

	pid_t pid;
	switch( p_mode )
	{
	case PHYMIR_SEQ:
		return p_lgLimitExceeded(x) ? p_ymax : p_func(x,runNr);
	case PHYMIR_FORK:
		p_curcpu++;
		if( p_curcpu > p_maxcpu )
		{
			/* too many processes, wait for one to finish */
			(void)wait(NULL);
			p_curcpu--;
		}
		/* flush all open files */
		fflush(NULL);
		pid = fork();
		if( pid < 0 )
		{
			fprintf( ioQQQ, "creating the child process failed\n" );
			cdEXIT(EXIT_FAILURE);
		}
		else if( pid == 0 )
		{
			/* this is child process so execute */
			p_execute_job_parallel( x, jj, runNr );
			/* at this point ioQQQ points to the main output of the parent process,
			 * the child should not close that in cdEXIT(), so wipe the handle here */
			ioQQQ = NULL;
			cdEXIT(EXIT_SUCCESS);
		}
		// the real y-value is not available yet...
		return p_ymax;
	case PHYMIR_MPI:
		if( (jj%cpu.i().nCPU()) == cpu.i().nRANK() )
			p_execute_job_parallel( x, jj, runNr );
		// the real y-value is not available yet...
		return p_ymax;
	default:
		TotalInsanity();
	}
}

// p_execute_job_parallel MUST be const, otherwise changes by child processes may be lost!
template<class X, class Y, int NP, int NSTR>
void phymir_state<X,Y,NP,NSTR>::p_execute_job_parallel(const X x[],
						       int jj,
						       int runNr) const
{
	DEBUG_ENTRY( "p_execute_job_parallel()" );

	char fnam1[20], fnam2[20];
	sprintf(fnam1,"yval_%d",jj);
	sprintf(fnam2,"output_%d",jj);
	/* each child should have its own output file */
	FILE *ioQQQ_old = ioQQQ;
	ioQQQ = open_data( fnam2, "w", AS_LOCAL_ONLY );
	// fail-safe in case p_func crashes without being caught...
	Y yval = p_ymax;
	wr_block(&yval,sizeof(yval),fnam1);
	if( !p_lgLimitExceeded(x) )
	{
		try
		{
			yval = p_func(x,runNr);
		}
		catch( ... )
		{
			// make sure output is complete since files may not have been closed properly...
			fflush(NULL);
			yval = p_ymax;
		}
		wr_block(&yval,sizeof(yval),fnam1);
	}
	fclose( ioQQQ );
	ioQQQ = ioQQQ_old;
}

template<class X, class Y, int NP, int NSTR>
void phymir_state<X,Y,NP,NSTR>::p_barrier(int jlo,
					  int jhi)
{
	DEBUG_ENTRY( "p_barrier()" );

	switch( p_mode )
	{
	case PHYMIR_SEQ:
		// nothing to do...
		break;
	case PHYMIR_FORK:
		/* wait for child processes to terminate */
		while( p_curcpu > 0 )
		{
			(void)wait(NULL);
			p_curcpu--;
		}
		p_process_output( jlo, jhi );
		break;
	case PHYMIR_MPI:
		// wait for all function evaluations to finish so that output is complete
		MPI::COMM_WORLD.Barrier();
		p_process_output( jlo, jhi );
		// y values are known now, so broadcast to other ranks
		MPI::COMM_WORLD.Bcast( &p_yp[jlo], jhi-jlo+1, MPI::type(p_yp), 0 );
		break;
	default:
		TotalInsanity();
	}
}

template<class X, class Y, int NP, int NSTR>
void phymir_state<X,Y,NP,NSTR>::p_process_output(int jlo,
						 int jhi)
{
	DEBUG_ENTRY( "p_process_output()" );

	if( cpu.i().lgMaster() )
	{
		char fnam[20];
		for( int jj=jlo; jj <= jhi; jj++ )
		{
			sprintf(fnam,"yval_%d",jj);
			rd_block(&p_yp[jj],sizeof(p_yp[jj]),fnam);
			remove(fnam);
			sprintf(fnam,"output_%d",jj);
			append_file(ioQQQ,fnam);
			remove(fnam);
		}
	}
}

template<class X, class Y, int NP, int NSTR>
void phymir_state<X,Y,NP,NSTR>::p_evaluate_hyperblock()
{
	DEBUG_ENTRY( "p_evaluate_hyperblock()" );

	int jlo = 1, jhi = 0;
	for( int j=0; j < p_nvar; j++ )
	{
		X sgn = X(1.);
		for( int jj=2*j+1; jj <= 2*j+2; jj++ )
		{
			sgn = -sgn;
			for( int i=0; i < p_nvar; i++ )
			{
				p_xp[jj][i] = p_xc[i] + sgn*p_dmax*p_c2[j]*p_a2[j][i];
				p_varmax[i] = max(p_varmax[i],p_xp[jj][i]);
				p_varmin[i] = min(p_varmin[i],p_xp[jj][i]);
			}
			if( !lgMaxIterExceeded() )
			{
				// p_execute_job will not return the correct chi^2 in parallel mode
				// only after p_barrier() has finished will p_yp[jj] contain the correct value
				p_yp[jj] = p_execute_job( p_xp[jj], jj, p_noptim++ );
				jhi = jj;
			}
		}
	}
	/* wait for jobs to terminate */
	p_barrier( jlo, jhi );
}

template<class X, class Y, int NP, int NSTR>
void phymir_state<X,Y,NP,NSTR>::p_setup_next_hyperblock()
{
	DEBUG_ENTRY( "p_setup_next_hyperblock()" );

	const Y F0 = Y(1.4142136);
	const X F1 = X(0.7071068);
	const X F2 = X(0.1);

	/* find best model */
	for( int jj=1; jj <= 2*p_nvar; jj++ )
	{
		if( p_yp[jj] < p_ymin )
		{
			p_ymin = p_yp[jj];
			p_jmin = jj;
		}
	}
	bool lgNewCnt = p_jmin > 0;

	/* determine minimum and relative uncertainties */
	bool lgNegd2 = false;
	X xnrm = X(0.);
	X xhlp[NP];
	for( int i=0; i < p_nvar; i++ )
	{
		Y d1 = p_yp[2*i+2] - p_yp[2*i+1];
		Y d2 = Y(0.5)*p_yp[2*i+2] - p_yp[0] + Y(0.5)*p_yp[2*i+1];
		if( d2 <= Y(0.) )
			lgNegd2 = true;
		xhlp[i] = -p_dmax*p_c2[i]*(static_cast<X>(max(min((Y(0.25)*d1)/max(d2,Y(1.e-10)),F0),-F0)) -
					   p_delta(2*i+1,p_jmin) + p_delta(2*i+2,p_jmin));
		xnrm += pow2(xhlp[i]);
	}
	xnrm = static_cast<X>(sqrt(xnrm));
	/* set up new transformation matrix */
	int imax = 0;
	X amax = X(0.);
	for( int j=0; j < p_nvar; j++ )
	{
		for( int i=0; i < p_nvar; i++ )
		{
			if( xnrm > X(0.) )
			{
				if( j == 0 )
				{
					p_a2[0][i] *= xhlp[0];
				}
				else
				{
					p_a2[0][i] += xhlp[j]*p_a2[j][i];
					p_a2[j][i] = p_delta(i,j);
					if( j == p_nvar-1 && abs(p_a2[0][i]) > amax )
					{
						imax = i;
						amax = abs(p_a2[0][i]);
					}
				}
			}
			else
			{
				p_a2[j][i] = p_delta(i,j);
			}
		}
	}
	/* this is to assure maximum linear independence of the base vectors */
	if( imax > 0 )
	{
		p_a2[imax][0] = X(1.);
		p_a2[imax][imax] = X(0.);
	}
	/* apply Gram-Schmidt procedure to get orthogonal base */
	p_phygrm( p_a2, p_nvar );

	/* set up hyperblock dimensions in new base and choose new center */
	for( int i=0; i < p_nvar; i++ )
	{
		p_c2[i] = X(0.);
		for( int j=0; j < p_nvar; j++ )
		{
			p_c2[i] += pow2(p_a2[i][j]/p_c1[j]);
		}
		p_c2[i] = static_cast<X>(1./sqrt(p_c2[i]));
		p_xc[i] = p_xp[p_jmin][i];
		p_xp[0][i] = p_xc[i];
	}
	p_yp[0] = p_yp[p_jmin];
	p_jmin = 0;
	/* choose size of next hyperblock */
	X r1, r2;
	if( lgNegd2 )
	{
		/* this means that the hyperblock either is bigger than the scale
		 * on which the function varies, or is so small that we see noise.
		 * in both cases make the hyperblock smaller. */
		r1 = F1;
		r2 = F1;
	}
	else
	{
		r1 = F2;
		if( lgNewCnt )
		{
			/* when we make progress, p_dmax may become bigger */
			r2 = sqrt(1.f/F1);
		}
		else
		{
			/* when there is no progress force p_dmax smaller, otherwise there
			 * may never be an end */
			r2 = sqrt(F1);
		}
	}
	p_dmax = min(max(xnrm/p_c2[0],p_dmax*r1),p_dmax*r2);
	/* fail-safe against excessive behaviour */
	p_dmax = MIN2(p_dmax,p_dold);
}

template<class X, class Y, int NP, int NSTR>
void phymir_state<X,Y,NP,NSTR>::p_reset_hyperblock()
{
	DEBUG_ENTRY( "p_reset_hyperblock()" );

	if( !lgConvergedRestart() )
	{
		/* reset hyperblock so that we can restart the optimization */
		for( int i=0; i < p_nvar; i++ )
		{
			p_xcold[i] = p_xc[i];
			p_c2[i] = p_c1[i];
		}
		p_dmax = p_dold;
		p_reset_transformation_matrix();
	}
}

template<class X, class Y, int NP, int NSTR>
void phymir_state<X,Y,NP,NSTR>::p_phygrm(X a[][NP],
					 int n)
{
	DEBUG_ENTRY( "p_phygrm()" );

	/* apply modified Gram-Schmidt procedure to a */
	for( int i=0; i < n; i++ )
	{
		X ip = X(0.);
		for( int k=0; k < n; k++ )
			ip += pow2(a[i][k]);
		ip = sqrt(ip);
		for( int k=0; k < n; k++ )
			a[i][k] /= ip;
		for( int j=i+1; j < n; j++ )
		{
			X ip = X(0.);
			for( int k=0; k < n; k++ )
				ip += a[i][k]*a[j][k];
			for( int k=0; k < n; k++ )
				a[j][k] -= ip*a[i][k];
		}
	}
	return;
}

template<class X, class Y, int NP, int NSTR>
bool phymir_state<X,Y,NP,NSTR>::p_lgLimitExceeded(const X x[]) const
{
	DEBUG_ENTRY( "p_lgLimitExceeded()" );

	for( int i=0; i < p_nvar; i++ )
	{
		if( x[i] < p_absmin[i] )
			return true;
		if( x[i] > p_absmax[i] )
			return true;
	}
	return false;
}

template<class X, class Y, int NP, int NSTR>
void phymir_state<X,Y,NP,NSTR>::p_reset_transformation_matrix()
{
	DEBUG_ENTRY( "p_reset_transformation_matrix()" );

	/* initialize transformation matrix to unity */
	for( int i=0; i < p_nvar; i++ )
		for( int j=0; j < p_nvar; j++ )
			p_a2[j][i] = p_delta(i,j);
}

template<class X, class Y, int NP, int NSTR>
void phymir_state<X,Y,NP,NSTR>::init_minmax(const X pmin[], // pmin[nvar]: minimum values for params
					    const X pmax[], // pmax[nvar]: maximum values for params
					    int nvar)       // will not be initialized yet
{
	DEBUG_ENTRY( "init_minmax()" );

	ASSERT( !lgInitialized() );

	for( int i=0; i < nvar; i++ )
	{
		p_absmin[i] = pmin[i];
		p_absmax[i] = pmax[i];
	}
}

template<class X, class Y, int NP, int NSTR>
void phymir_state<X,Y,NP,NSTR>::init_strings(const char* date,      // date string for inclusion in state file
					     const char* version,   // version string for inclusion in state file
					     const char* host_name) // host name for inclusion in state file
{
	DEBUG_ENTRY( "init_strings()" );

	// use NSTR-1 so that last 0 byte is not overwritten...
	if( date != NULL )
		strncpy( p_chStr1, date, NSTR-1 );
	if( version != NULL )
		strncpy( p_chStr2, version, NSTR-1 );
	if( host_name != NULL )
		strncpy( p_chStr3, host_name, NSTR-1 );
}

template<class X, class Y, int NP, int NSTR>
void phymir_state<X,Y,NP,NSTR>::init_state_file_name(const char* fnam) // name of the state file that will be written
{
	DEBUG_ENTRY( "init_state_file_name()" );

	// use NSTR-1 so that last 0 byte is not overwritten...
	strncpy( p_chState, fnam, NSTR-1 );
}

template<class X, class Y, int NP, int NSTR>
void phymir_state<X,Y,NP,NSTR>::initial_run(Y (*func)(const X[],int),  // function to be optimized
					    int nvar,                  // number of parameters to be optimized
					    const X start[],           // start[n]: initial estimates for params
					    const X del[],             // del[n]: initial stepsizes for params
					    X toler,                   // absolute tolerance on parameters
					    int maxiter,               // maximum number of iterations
					    phymir_mode mode,          // execution mode: sequential, fork, mpi
					    int maxcpu)                // maximum no. of CPUs to use
{
	DEBUG_ENTRY( "initial_run()" );

	ASSERT( nvar > 0 && nvar <= NP );

	p_func = func;
	p_nvar = nvar;
	p_toler = toler;
	p_maxiter = maxiter;
	p_mode = mode;
	p_maxcpu = maxcpu;
	p_noptim = 0;

	/* initialize hyperblock dimensions and center */
	p_dmax = X(0.);
	for( int i=0; i < p_nvar; i++ )
		p_dmax = max(p_dmax,abs(del[i]));

	p_dold = p_dmax;
	for( int i=0; i < p_nvar; i++ )
	{
		p_xc[i] = start[i];
		p_xcold[i] = p_xc[i] + X(10.)*p_toler;
		p_c1[i] = abs(del[i])/p_dmax;
		p_c2[i] = p_c1[i];
		p_xp[0][i] = p_xc[i];
		p_varmax[i] = max(p_varmax[i],p_xc[i]);
		p_varmin[i] = min(p_varmin[i],p_xc[i]);
	}

	// p_execute_job will not return the correct chi^2 in parallel mode
	// only after p_barrier() has finished will p_yp[0] contain the correct value
	p_yp[0] = p_execute_job( p_xc, 0, p_noptim++ );
	p_barrier( 0, 0 );

	p_ymin = p_yp[0];
	p_jmin = 0;

	p_reset_transformation_matrix();

	p_wr_state( p_chState );
}

template<class X, class Y, int NP, int NSTR>
void phymir_state<X,Y,NP,NSTR>::continue_from_state(Y (*func)(const X[],int), // function to be optimized	
						    int nvar,                 // number of parameters to be optimized
						    const char* fnam,         // file name of the state file
						    X toler,                  // absolute tolerance on parameters
						    int maxiter,              // maximum number of iterations
						    phymir_mode mode,         // execution mode: sequential, fork, mpi
						    int maxcpu)               // maximum no. of CPUs to use
{
	DEBUG_ENTRY( "continue_from_state()" );

	p_rd_state( fnam );
	// sanity checks
	if( !fp_equal( p_vers, VRSNEW ) )
	{
		printf( "optimize continue - file has incompatible version, sorry\n" );
		cdEXIT(EXIT_FAILURE);
	}
	if( p_dim != NP )
	{
		printf( "optimize continue - arrays have wrong dimension, sorry\n" );
		cdEXIT(EXIT_FAILURE);
	}
	if( p_sdim != NSTR )
	{
		printf( "optimize continue - strings have wrong length, sorry\n" );
		cdEXIT(EXIT_FAILURE);
	}
	if( p_nvar != nvar )
	{
		printf( "optimize continue - wrong number of free parameters, sorry\n" );
		cdEXIT(EXIT_FAILURE);
	}
	// this pointer was not part of the state file, it would be useless...
	p_func = func;
	// Option to override the tolerance set in the state file. This is useful
	// if you want to refine an already finished run to a smaller tolerance.
	p_toler = toler;
	// you ran out of iterations, but wanted to continue...
	p_maxiter = maxiter;
	// it is OK to continue in a different mode
	p_mode = mode;
	p_maxcpu = maxcpu;
}

template<class X, class Y, int NP, int NSTR>
void phymir_state<X,Y,NP,NSTR>::optimize()
{
	DEBUG_ENTRY( "optimize()" );

	ASSERT( lgInitialized() );

	while( !lgConverged() )
	{
		p_evaluate_hyperblock();
		if( lgMaxIterExceeded() )
			break;
		p_setup_next_hyperblock();
		p_wr_state( p_chState );
	}
}

template<class X, class Y, int NP, int NSTR>
void phymir_state<X,Y,NP,NSTR>::optimize_with_restart()
{
	DEBUG_ENTRY( "optimize_with_restart()" );

	ASSERT( lgInitialized() );

	while( !lgConvergedRestart() )
	{
		optimize();
		if( lgMaxIterExceeded() )
			break;
		p_reset_hyperblock();
	}
}

template<class X, class Y, int NP, int NSTR>
bool phymir_state<X,Y,NP,NSTR>::lgConvergedRestart() const
{
	DEBUG_ENTRY( "lgConvergedRestart()" );

	if( lgConverged() )
	{
		X dist = X(0.);
		for( int i=0; i < p_nvar; i++ )
			dist += pow2(p_xc[i] - p_xcold[i]);
		dist = static_cast<X>(sqrt(dist));
		return ( dist <= p_toler );
	}
	return false;
}

void optimize_phymir(realnum xc[], 
		     const realnum del[], 
		     long int nvarPhymir, 
		     chi2_type *ymin, 
		     realnum toler)
{
	DEBUG_ENTRY( "optimize_phymir()" );

	if( nvarPhymir > LIMPAR )
	{
		fprintf( ioQQQ, "optimize_phymir: too many parameters are varied, increase LIMPAR\n" );
		cdEXIT(EXIT_FAILURE);
	}

	phymir_state<realnum,chi2_type,LIMPAR,STDLEN> phymir;

	// first check if a statefile already exists, back it up in that case
	(void)remove(STATEFILE_BACKUP);
	FILE *tmp = open_data( STATEFILE, "r", AS_LOCAL_ONLY_TRY );
	if( tmp != NULL )
	{
		fclose( tmp );
		// create backup copy of statefile
		FILE *dest = open_data( STATEFILE_BACKUP, "wb", AS_LOCAL_ONLY_TRY );
		if( dest != NULL )
		{
			append_file( dest, STATEFILE );
			fclose( dest );
		}
	}

	phymir_mode mode = optimize.lgParallel ? ( cpu.i().lgMPI() ? PHYMIR_MPI : PHYMIR_FORK ) : PHYMIR_SEQ;
	long nCPU = optimize.lgParallel ? ( cpu.i().lgMPI() ? cpu.i().nCPU() : optimize.useCPU ) : 1;
	if( optimize.lgOptCont )
	{
		phymir.continue_from_state( optimize_func, nvarPhymir, STATEFILE, toler,
					    optimize.nIterOptim, mode, nCPU );
	}
	else
	{
		phymir.init_state_file_name( STATEFILE );
		phymir.init_strings( t_version::Inst().chDate, t_version::Inst().chVersion,
				     cpu.i().host_name() );
		phymir.initial_run( optimize_func, nvarPhymir, xc, del, toler,
				    optimize.nIterOptim, mode, nCPU );
	}

	phymir.optimize_with_restart();

	if( phymir.lgMaxIterExceeded() )
	{
		fprintf( ioQQQ, " Optimizer exceeding maximum iterations.\n" );
		fprintf( ioQQQ, " This can be reset with the OPTIMIZE ITERATIONS command.\n" );
	}

	// updates to optimize.nOptimiz and optimize.varmin/max in child processes are lost, so repair here...
	optimize.nOptimiz = phymir.noptim();
	for( int i=0; i < nvarPhymir; i++ )
	{
		xc[i] = phymir.xval(i);
		optimize.varmax[i] = min(phymir.xmax(i),optimize.varang[i][1]);
		optimize.varmin[i] = max(phymir.xmin(i),optimize.varang[i][0]);
	}
	*ymin = phymir.yval();
	return;
}

/** wr_block: write <len> bytes of data from buffer <*ptr> into unformatted file <fnam> */
inline void wr_block(const void *ptr,
		     size_t len,
		     const char *fnam)
{
	DEBUG_ENTRY( "wr_block()" );

	FILE *fdes = open_data( fnam, "wb", AS_LOCAL_ONLY );
	if( fwrite(ptr,len,size_t(1),fdes) != 1 ) {
		printf( "error writing on file: %s\n",fnam );
		fclose(fdes);
		cdEXIT(EXIT_FAILURE);
	}
	fclose(fdes);
	return;
}

/** rd_block: read <len> bytes of data into buffer <*ptr> from unformatted file <fnam> */
inline void rd_block(void *ptr,
		     size_t len,
		     const char *fnam)
{
	DEBUG_ENTRY( "rd_block()" );

	FILE *fdes = open_data( fnam, "rb", AS_LOCAL_ONLY );
	if( fread(ptr,len,size_t(1),fdes) != 1 ) {
		printf( "error reading on file: %s\n",fnam );
		fclose(fdes);
		cdEXIT(EXIT_FAILURE);
	}
	fclose(fdes);
	return;
}