/*****************************************************************************
   Major portions of this software are copyrighted by the Medical College
   of Wisconsin, 1994-2000, and are released under the Gnu General Public
   License, Version 2.  See the file README.Copyright for details.
******************************************************************************/

/*
   This file implements the Simplex algorithm for non-linear optimization.
   Some routines are adapted from:  A Jump Start Course in C++ Programming

  File:     simplex.c
  Author:   B. Douglas Ward
  Date:     23 May 1997

  Mod:      Changed random number generator function from rand to drand48.
            29 August 1997
*/




/*---------------------------------------------------------------------------*/
/*
  Routine to generate a uniform U(0,1) random variate.
*/

float uniform ()
{
  return ( (float)drand48() );
}


/*---------------------------------------------------------------------------*/
/*
  Routine to generate a normal N(0,1) random variate.
*/

void normal (float * n1, float * n2)
{
  float u1, u2;
  float r;


  u1 = 0.0;
  while (u1 <= 0.0)
    {
      u1 = uniform();
    }
  u2 = uniform();

  r   = sqrt(-2.0*log(u1));
  *n1 = r * cos(2.0*PI*u2);
  *n2 = r * sin(2.0*PI*u2);
}


/*---------------------------------------------------------------------------*/
/*
  Routine to generate a U(a,b) random variate.
*/

float get_random_value(float a, float b)
{
  float fval;

  fval = a + uniform() * (b-a);

  return (fval);
}


/*---------------------------------------------------------------------------*/
/*
  Allocate memory for simplex algorithm.
*/

void allocate_arrays
(
  int dimension,              /* dimension of parameter space */
  float *** simplex,          /* the simplex itself */
  float ** centroid,          /* center of mass of the simplex */
  float ** response,          /* error sum of squares at each vertex */
  float ** step_size,         /* controls random placement of new vertex */
  float ** test1,             /* test vertex */
  float ** test2              /* test vertex */
)

{
  int i;

  *centroid = (float *) malloc (sizeof(float) * dimension);
  *step_size = (float *) malloc (sizeof(float) * dimension);
  *test1 = (float *) malloc (sizeof(float) * dimension);
  *test2 = (float *) malloc (sizeof(float) * dimension);

  *response = (float *) malloc (sizeof(float) * (dimension+1));
  *simplex = (float **) malloc (sizeof(float *) * (dimension+1));

  for (i = 0;  i < dimension+1;  i++)
    (*simplex)[i] = (float *) malloc (sizeof(float) * dimension);

}


/*---------------------------------------------------------------------------*/
/*
  Set up initial values for the simplex vertices.
*/

void initialize_simplex
(
  int dimension,          /* dimension of the full model */
  vfp nmodel,             /* pointer to noise model */
  vfp smodel,             /* pointer to signal model */
  int r,                  /* number of parameters in the noise model */
  int p,                  /* number of parameters in the signal model */
  int nabs,               /* use absolute constraints for noise parameters */
  float * min_nconstr,    /* minimum parameter constraints for noise model */
  float * max_nconstr,    /* maximum parameter constraints for noise model */
  float * min_sconstr,    /* minimum parameter constraints for signal model */
  float * max_sconstr,    /* maximum parameter constraints for signal model */
  float * par_rdcd,       /* estimated parameters for the reduced model */
  float * parameters,     /* starting point */
  float ** simplex,       /* the simplex itself */
  float * response,       /* sse at each vertex of the simplex */
  float * step_size,      /* amount of allowed variation at each parameter */
  int ts_length,          /* length of time series array */
  float ** x_array,       /* independent variable matrix */
  float * ts_array        /* observed time series */
)

{
  int i, j;
  float minval, maxval;


  /*----- copy parameter vector into first vertex of simplex -----*/
  for (i = 0;  i < dimension;  i++)
    simplex[0][i] = parameters[i];


  /*----- set up initial step sizes -----*/
  for (i = 0;  i < r;  i++)
    step_size[i] = 0.1 * (max_nconstr[i] - min_nconstr[i]);
  for (i = r;  i < dimension;  i++)
    step_size[i] = 0.1 * (max_sconstr[i-r] - min_sconstr[i-r]);


  /*----- choose random vectors for remaining vertices -----*/
  for (i = 1;  i < dimension+1;  i++)
    {
      /*----- choose noise model parameters -----*/
      for (j = 0;  j < r;  j++)
	{
	  minval = simplex[0][j] - step_size[j];
	  if (nabs)   /*--- absolute noise parameter constraints ---*/
	    {
	      if (minval < min_nconstr[j])
		minval = min_nconstr[j];
	    }
	  else        /*--- relative noise parameter constraints ---*/
	    {
	      if (minval < min_nconstr[j] + par_rdcd[j])
		minval = min_nconstr[j] + par_rdcd[j];
	    }

	  maxval = simplex[0][j] + step_size[j];
	  if (nabs)   /*--- absolute noise parameter constraints ---*/
	    {
	      if (maxval > max_nconstr[j])
		maxval = max_nconstr[j];
	    }
	  else        /*--- relative noise parameter constraints ---*/
	    {
	      if (maxval > max_nconstr[j] + par_rdcd[j])
		maxval = max_nconstr[j] + par_rdcd[j];
	    }
	  simplex[i][j] = get_random_value (minval, maxval);
	}


      /*----- choose signal model parameters -----*/
      for (j = r;  j < dimension;  j++)
	{
	  minval = simplex[0][j] - step_size[j];
	  if (minval < min_sconstr[j-r])
	    minval = min_sconstr[j-r];
	  maxval = simplex[0][j] + step_size[j];
	  if (maxval > max_sconstr[j-r])
	    maxval = max_sconstr[j-r];
	  simplex[i][j] = get_random_value (minval, maxval);
	}
    }


  /*----- calculate and save sse for each vertex of simplex -----*/
  for (i = 0;  i < dimension+1;  i++)
    response[i] = calc_sse(nmodel, smodel, r, p, nabs,
			   min_nconstr, max_nconstr, min_sconstr, max_sconstr,
			   par_rdcd, simplex[i], ts_length, x_array, ts_array);

}


/*---------------------------------------------------------------------------*/
/*
  Evaluate the vertices of the simplex.  Find indices of the best, worst, and
  next-to-worst vertices.
*/

void eval_vertices
(
 int dimension,            /* dimension of parameter space */
 float * response,         /* error sum of squares at each vertex */
 int * worst,              /* index of worst vertex in simplex */
 int * next,               /* index of next-to-worst vertex in simplex */
 int * best                /* index of best vertex in simplex */
)

{
  int i;

  /* initialize values */
  *worst = 0;
  *best = 0;

  /* find the best and worst */
  for (i = 1;  i < dimension+1;  i++)
    {
      if (response[i] > response[*worst])
	*worst = i;
      if (response[i] < response[*best])
	*best = i;
    }

  /* find the next worst index */
  if (*worst == 0)
    *next = 1;
  else
    *next = 0;

  for (i = 0;  i < dimension+1;  i++)
    if ((i != *worst) && (response[i] > response[*next]))
      *next = i;
}


/*---------------------------------------------------------------------------*/
/*
  Routine to restart the simplex algorithm by putting random vectors into
  the simplex (and keeping the best previous vertex).
*/

void restart
(
  int dimension,          /* dimension of the full model */
  vfp nmodel,             /* pointer to noise model */
  vfp smodel,             /* pointer to signal model */
  int r,                  /* number of parameters in the noise model */
  int p,                  /* number of parameters in the signal model */
  int nabs,               /* use absolute constraints for noise parameters */
  float * min_nconstr,    /* minimum parameter constraints for noise model */
  float * max_nconstr,    /* maximum parameter constraints for noise model */
  float * min_sconstr,    /* minimum parameter constraints for signal model */
  float * max_sconstr,    /* maximum parameter constraints for signal model */
  float * par_rdcd,       /* estimated parameters for the reduced model */
  float ** simplex,       /* the simplex itself */
  float * response,       /* sse at each vertex of the simplex */
  float * step_size,      /* amount of allowed variation at each parameter */
  int ts_length,          /* length of time series array */
  float ** x_array,       /* independent variable matrix */
  float * ts_array        /* observed time series */
)

{
  const float STEP_FACTOR = 0.9;
  int i, j;
  int worst, next, best;
  float minval, maxval;


  /* find the current best vertex */
  eval_vertices (dimension, response, &worst, &next, &best);


  /* set the first vertex to the current best */
  for (i = 0; i < dimension;  i++)
    simplex[0][i] = simplex[best][i];


  /* decrease step size */
  for (i = 0;  i < dimension;  i++)
    step_size[i] *= STEP_FACTOR;


  /* set up remaining vertices of simplex using new step size */
  for (i = 1;  i < dimension+1;  i++)
    {
      /*----- choose noise model parameters -----*/
      for (j = 0;  j < r;  j++)
	{
	  minval = simplex[0][j] - step_size[j];
	  if (nabs)   /*--- absolute noise parameter constraints ---*/
	    {
	      if (minval < min_nconstr[j])
		minval = min_nconstr[j];
	    }
	  else        /*--- relative noise parameter constraints ---*/
	    {
	      if (minval < min_nconstr[j] + par_rdcd[j])
		minval = min_nconstr[j] + par_rdcd[j];
	    }

	  maxval = simplex[0][j] + step_size[j];
	  if (nabs)
	    {
	      if (maxval > max_nconstr[j])
		maxval = max_nconstr[j];
	    }
	  else
	    {
	      if (maxval > max_nconstr[j] + par_rdcd[j])
		maxval = max_nconstr[j] + par_rdcd[j];
	    }

	  simplex[i][j] = get_random_value (minval, maxval);
	}


      /*----- choose signal model parameters -----*/
      for (j = r;  j < dimension;  j++)
	{
	  minval = simplex[0][j] - step_size[j];
	  if (minval < min_sconstr[j-r])
	    minval = min_sconstr[j-r];
	  maxval = simplex[0][j] + step_size[j];
	  if (maxval > max_sconstr[j-r])
	    maxval = max_sconstr[j-r];
	  simplex[i][j] = get_random_value (minval, maxval);
	}
    }


  /* initialize response for each vector */
  for (i = 0;  i < dimension+1;  i++)
    response[i] = calc_sse (nmodel, smodel, r, p, nabs,
			   min_nconstr, max_nconstr, min_sconstr, max_sconstr,
			   par_rdcd, simplex[i], ts_length, x_array, ts_array);
}


/*---------------------------------------------------------------------------*/
/*
  Calculate the centroid of the simplex, ignoring the worst vertex.
*/

void calc_centroid
(
  int dimension,         /* dimension of parameter space */
  float ** simplex,      /* the simplex itself */
  int worst,             /* index of worst vertex in simplex */
  float * centroid       /* center of mass of the simplex minus worst vertex */
)

{
  int i, j;

  for (i = 0;  i < dimension;  i++)
    {
      centroid[i] = 0.0;

      /* add each vertex, except the worst */
      for (j = 0;  j < dimension+1;  j++)
	if (j != worst)
	  centroid[i] += simplex[j][i];
    }

  /* divide by the number of vertices */
  for (i = 0;  i < dimension;  i++)
    centroid[i] /= dimension;
}


/*---------------------------------------------------------------------------*/
/*
  Calculate the reflection of the worst vertex about the centroid.
*/

void calc_reflection
(
  int dimension,               /* dimension of parameter space */
  float ** simplex,            /* the simplex itself */
  float * centroid,            /* center of mass of the simplex */
  int worst,                   /* index of worst vertex in simplex */
  float coef,                  /* expansion or contraction factor */
  float * vertex               /* new vertex */
)

{
  int i;

  for (i = 0;  i < dimension;  i++)
    vertex[i] = centroid[i] + coef*(centroid[i] - simplex[worst][i]);
}


/*---------------------------------------------------------------------------*/
/*
  Replace a vertex of the simplex.
*/

void replace
(
  int dimension,              /* dimension of parameter space */
  float ** simplex,           /* the simplex itself */
  float * response,           /* error sum of squares at each vertex */
  int index,                  /* index of vertex to be replaced */
  float * vertex,             /* new vertex */
  float resp                  /* error sum of squares at new vertex */
)

{
  int i;

  for (i = 0;  i < dimension;  i++)
    simplex[index][i] = vertex[i];

  response[index] = resp;
}


/*---------------------------------------------------------------------------*/
/*
  Calculate the goodness of fit.  This is measured by the variation in
  responses at the different vertices relative to the average response.
*/

float calc_good_fit
(
  int dimension,                   /* dimension of parameter space */
  float * response                 /* error sum of squares at each vertex */
)

{
  int i;

  float avg, sd, tmp;

  /*----- average the response values -----*/
  avg = 0.0;
  for (i = 0;  i < dimension+1;  i++)
    avg += response[i];
  avg /= dimension+1;

  /*----- compute standard deviation of response -----*/
  sd = 0.0;
  for (i = 0;  i < dimension+1;  i++)
    {
      tmp = response[i] - avg;
      sd += tmp*tmp;
    }

  sd /= dimension;

  return (sqrt(sd) / avg);
}


/*---------------------------------------------------------------------------*/
/*
  Release memory required for simplex optimization.
*/

void deallocate_arrays
(
  int dimension,              /* dimension of parameter space */
  float *** simplex,          /* the simplex itself */
  float ** centroid,          /* center of mass of the simplex */
  float ** response,          /* error sum of squares at each vertex */
  float ** step_size,         /* controls random placement of new vertex */
  float ** test1,             /* test vertex */
  float ** test2              /* test vertex */
)

{
  int iv;                     /* vertex index */


  free (*centroid);    *centroid = NULL;
  free (*response);    *response = NULL;
  free (*step_size);   *step_size = NULL;
  free (*test1);       *test1 = NULL;
  free (*test2);       *test2 = NULL;

  for (iv = 0;  iv < dimension+1;  iv++)
    {
      free ((*simplex)[iv]);
      (*simplex)[iv] = NULL;
    }

  free (*simplex);     *simplex = NULL;

}


/*---------------------------------------------------------------------------*/
/*
  Implementation of the (non-linear) simplex optimization algorithm.
*/

void simplex_optimization
(
  vfp nmodel,             /* pointer to noise model */
  vfp smodel,             /* pointer to signal model */
  int r,                  /* number of parameters in the noise model */
  int p,                  /* number of parameters in the signal model */
  float * min_nconstr,    /* minimum parameter constraints for noise model */
  float * max_nconstr,    /* maximum parameter constraints for noise model */
  float * min_sconstr,    /* minimum parameter constraints for signal model */
  float * max_sconstr,    /* maximum parameter constraints for signal model */
  int nabs,               /* use absolute constraints for noise parameters */
  int ts_length,          /* length of time series array */
  float ** x_array,       /* independent variable matrix */
  float * ts_array,       /* observed time series */
  float * par_rdcd,       /* estimated parameters for the reduced model */
  float * parameters,     /* estimated parameters */
  float * sse             /* error sum of squares */
)

{
  const int MAX_ITERATIONS = 50;         /* maximum number of iterations */
  const int MAX_RESTARTS = 5;            /* maximum number of restarts */
  const float EXPANSION_COEF = 2.0;      /* expansion coefficient */
  const float REFLECTION_COEF = 1.0;     /* reflection coefficient */
  const float CONTRACTION_COEF = 0.5;    /* contraction coefficient */
  const float TOLERANCE = 1.0e-4;        /* solution convergence tolerance */

  float ** simplex   = NULL;    /* the simplex itself */
  float * centroid   = NULL;    /* center of mass of the simplex */
  float * response   = NULL;    /* error sum of squares at each vertex */
  float * step_size  = NULL;    /* controls random placement of new vertex */
  float * test1      = NULL;    /* test vertex */
  float * test2      = NULL;    /* test vertex */
  float resp1, resp2;           /* error sum of squares for test vertex */
  int i;                        /* vertex index */
  int worst;                    /* index of worst vertex in simplex */
  int best;                     /* index of best vertex in simplex */
  int next;                     /* index of next-to-worst vertex in simplex */
  int num_iter;                 /* number of simplex algorithm iterations */
  int num_restarts;             /* number of restarts of simplex algorithm */
  int done;                     /* boolean for search finished */
  float fit;                    /* array of fitted time series values */
  int dimension;                /* dimension of parameter space */


  /*----- dimension of parameter space -----*/
  dimension = r + p;

  /*----- allocate memory -----*/
  allocate_arrays (dimension, &simplex, &centroid, &response, &step_size,
		   &test1, &test2);

  /*----- initialization for simplex algorithm -----*/
  initialize_simplex (dimension, nmodel, smodel, r, p, nabs,
		      min_nconstr, max_nconstr, min_sconstr, max_sconstr,
		      par_rdcd, parameters, simplex, response, step_size,
		      ts_length, x_array, ts_array);

  /* start loop to do simplex optimization */
  num_iter = 0;
  num_restarts = 0;
  done = 0;

  while (!done)
    {
      /*----- find the worst vertex and compute centroid of remaining simplex,
	discarding the worst vertex -----*/
      eval_vertices (dimension, response, &worst, &next, &best);
      calc_centroid (dimension, simplex, worst, centroid);

      /*----- reflect the worst point through the centroid -----*/
      calc_reflection (dimension, simplex, centroid, worst,
		       REFLECTION_COEF, test1);
      resp1 = calc_sse (nmodel, smodel, r, p, nabs, min_nconstr, max_nconstr,
			min_sconstr, max_sconstr, par_rdcd, test1,
			ts_length, x_array, ts_array);


      /*----- test the reflection against the best vertex and expand it if the
	reflection is better.  if not, keep the reflection -----*/
      if (resp1 < response[best])
	{
	  /*----- try expanding -----*/
	  calc_reflection (dimension, simplex, centroid, worst,
			   EXPANSION_COEF, test2);

	  resp2 = calc_sse (nmodel, smodel, r, p, nabs,
			    min_nconstr, max_nconstr,
			    min_sconstr, max_sconstr, par_rdcd, test2,
			    ts_length, x_array, ts_array);

	  if (resp2 <= resp1)    /* keep expansion */
	    replace (dimension, simplex, response, worst, test2, resp2);
	  else                   /* keep reflection */
	    replace (dimension, simplex, response, worst, test1, resp1);
	}

      else if (resp1 < response[next])
	{
	  /*----- new response is between the best and next worst
	    so keep reflection -----*/
	  replace (dimension, simplex, response, worst, test1, resp1);
	}
          else
	{
	  /*----- try contraction -----*/
	  if (resp1 >= response[worst])
	    calc_reflection (dimension, simplex, centroid, worst,
			     -CONTRACTION_COEF, test2);
	  else
	    calc_reflection (dimension, simplex, centroid, worst,
			     CONTRACTION_COEF, test2);

	  resp2 = calc_sse (nmodel, smodel, r, p, nabs,
			    min_nconstr, max_nconstr,
			    min_sconstr, max_sconstr, par_rdcd, test2,
			    ts_length, x_array, ts_array);
	
	  /*---- test the contracted response against the worst response ----*/
	  if (resp2 > response[worst])
	    {
	      /*----- new contracted response is worse, so decrease step size
		and restart -----*/
	      num_iter = 0;
	      num_restarts += 1;
	      restart (dimension, nmodel, smodel, r, p, nabs,
		       min_nconstr, max_nconstr, min_sconstr, max_sconstr,
		       par_rdcd, simplex, response, step_size,
		       ts_length, x_array, ts_array);
	    }
	  else       /*----- keep contraction -----*/
	    replace (dimension, simplex, response, worst, test2, resp2);
	}

      /*----- test to determine when to stop.
	first, check the number of iterations -----*/
      num_iter += 1;    /*----- increment iteration counter -----*/
      if (num_iter >= MAX_ITERATIONS)
	{
	  /*----- restart with smaller steps -----*/
	  num_iter = 0;
	  num_restarts += 1;
	  restart (dimension, nmodel, smodel, r, p, nabs,
		   min_nconstr, max_nconstr, min_sconstr, max_sconstr,
		   par_rdcd, simplex, response, step_size,
		   ts_length, x_array, ts_array);
	}

      /*----- limit the number of restarts -----*/
      if (num_restarts == MAX_RESTARTS)  done = 1;

      /*----- compare relative standard deviation of vertex responses
	against a defined tolerance limit -----*/
      fit = calc_good_fit (dimension, response);
      if (fit <= TOLERANCE)  done = 1;

      /*----- if done, copy the best solution to the output array -----*/
      if (done)
	{
	  eval_vertices (dimension, response, &worst, &next, &best);
	  for (i = 0;  i < dimension;  i++)
	    parameters[i] = simplex[best][i];
	  *sse = response[best];
	}

    }  /*----- while (!done) -----*/

  deallocate_arrays (dimension, &simplex, &centroid, &response, &step_size,
		     &test1, &test2);

}

/*----------------------------------------------------------------------------*/
/******************************************************************************/
/***** Powell's NEWUOA Method (instead of Nelder-Mead) - cf. powell_int.c *****/
/******************************************************************************/
/*----------------------------------------------------------------------------*/

static int N_newuoa = 0 ;   /* indicates if NEWUOA method is to be used */

static vfp N_nmodel , N_smodel ;
static int N_r , N_p , N_nabs, N_ts_length ;
static float *N_min_nconstr , *N_max_nconstr ;
static float *N_min_sconstr , *N_max_sconstr ;
static float **N_x_array ;
static float *N_ts_array , *N_par_rdcd ;
static float *N_pbot , *N_psiz ;

static float *N_pv ;

#if defined(SOLARIS) && !defined(floorf)
#define floorf floor   /* is Solaris lame, or what? */
#endif

/* Macro to periodically reduce a float variable into the range 0..1:
   for example: PRED01(1.2) == 0.8, PRED01(1.8) == 0.2, et cetera;
   graphically
               PRED01(x)|
                        | /\      /\      /\      /\      /\
                        |/  \    /  \    /  \    /  \    /
                        |    \  /    \  /    \  /    \  /
                        |     \/      \/      \/      \/
                        +------------------------------------> x
                          -3  -2  -1   0  +1  +2  +3  +4  +5
*/

#undef  PRED01
#define PRED01(x) fabsf( (x) - 2.0f*floorf(0.5f*((x)+1.0f)) )

/* double precision version of the above */

#undef  DRED01
#define DRED01(x) fabs ( (x) - 2.0 *floor (0.5 *((x)+1.0 )) )

/*----------------------------------------------------------------------------*/

double newfunc( int np , double *pv )  /* parameters are scaled to [0,1] */
{
   double val ; int ii ; float x,y ;

   for( ii=0 ; ii < np ; ii++ ){
     x = (float)pv[ii] ;
     if( x < 0.0f || x > 1.0f ) x = PRED01(x); /* reduce to [0,1] range */
     N_pv[ii] = N_pbot[ii] + N_psiz[ii] * x ;  /* scale to true value  */
   }

   /* compute sum of squares between model fit and data */

   val = (double) calc_sse( N_nmodel, N_smodel, N_r, N_p, N_nabs,
                            N_min_nconstr, N_max_nconstr,
                            N_min_sconstr, N_max_sconstr,
                            N_par_rdcd, N_pv ,
                            N_ts_length, N_x_array, N_ts_array ) ;
   return val ;
}

/*----------------------------------------------------------------------------*/

static double N_rstart=0.04 ;
static double N_rend  =0.0005 ;
static int    N_maxit =9999 ;

static void set_newuoa_parm( double rs , double re , int mm )
{
   if( rs > re && re > 0.0 && mm > 9 ){
     N_rstart = rs ;
     N_rend   = re ;
     N_maxit  = mm ;
   } else {
     N_rstart = 0.04 ;
     N_rend   = 0.0005 ;
     N_maxit  = 9999 ;
   }
}

/*----------------------------------------------------------------------------*/
/*! Supposed to be a dropin replacement for simplex_optimization().
------------------------------------------------------------------------------*/

void newuoa_optimization
(
  vfp nmodel,             /* pointer to noise model */
  vfp smodel,             /* pointer to signal model */
  int r,                  /* number of parameters in the noise model */
  int p,                  /* number of parameters in the signal model */
  float * min_nconstr,    /* minimum parameter constraints for noise model */
  float * max_nconstr,    /* maximum parameter constraints for noise model */
  float * min_sconstr,    /* minimum parameter constraints for signal model */
  float * max_sconstr,    /* maximum parameter constraints for signal model */
  int nabs,               /* use absolute constraints for noise parameters */
  int ts_length,          /* length of time series array */
  float ** x_array,       /* independent variable matrix */
  float * ts_array,       /* observed time series */
  float * par_rdcd,       /* estimated parameters for the reduced model */
  float * parameters,     /* estimated parameters */
  float * sse             /* error sum of squares */
)
{
  double *dv ; int ii ;

  N_nmodel      = nmodel ;
  N_smodel      = smodel ;
  N_r           = r ;
  N_p           = p ;
  N_min_nconstr = min_nconstr ;
  N_max_nconstr = max_nconstr ;
  N_min_sconstr = min_sconstr ;
  N_max_sconstr = max_sconstr ;
  N_nabs        = nabs ;
  N_ts_length   = ts_length ;
  N_x_array     = x_array ;
  N_ts_array    = ts_array ;
  N_par_rdcd    = par_rdcd ;

  N_pv          = (float *) malloc(sizeof(float) *(r+p)) ;
  N_pbot        = (float *) malloc(sizeof(float) *(r+p)) ;
  N_psiz        = (float *) malloc(sizeof(float) *(r+p)) ;
  dv            = (double *)malloc(sizeof(double)*(r+p)) ;

  /* set N_pbot = min allowed values
         N_psiz = max-min = size of allowed values range */

  if( nabs ){
    for( ii=0 ; ii < r ; ii++ ){
      N_pbot[ii] = min_nconstr[ii] ;
      N_psiz[ii] = max_nconstr[ii] - min_nconstr[ii] ;
    }
  } else {
    for( ii=0 ; ii < r ; ii++ ){
      N_pbot[ii] = min_nconstr[ii] + par_rdcd[ii] ;
      N_psiz[ii] = max_nconstr[ii] - min_nconstr[ii] ;
    }
  }
  for( ii=0 ; ii < p ; ii++ ){
    N_pbot[ii+r] = min_sconstr[ii] ;
    N_psiz[ii+r] = max_sconstr[ii] - min_sconstr[ii] ;
  }

  /* scale all parameters into the range [0,1];
     NEWUOA will operate on these, since it is scale-free */

  for( ii=0 ; ii < r+p ; ii++ ){
    if(N_psiz[ii] == 0.0) dv[ii] = 0.0;
    else
       dv[ii] = (double) ((parameters[ii]-N_pbot[ii])/N_psiz[ii]);
    if( dv[ii] < 0.0 || dv[ii] > 0.0 ) dv[ii] = DRED01(dv[ii]); /* 03 Nov 2006 */
  }

  powell_newuoa( r+p , dv , N_rstart , N_rend , N_maxit , newfunc ) ;

  *sse = (float)newfunc( r+p , dv ) ;

  for( ii=0 ; ii < r+p ; ii++ ){
    if( dv[ii] < 0.0 || dv[ii] > 0.0 ) dv[ii] = DRED01(dv[ii]); /* 03 Nov 2006 */
    parameters[ii] = (float)( N_pbot[ii] + N_psiz[ii]*dv[ii] );
  }

  free((void *)dv)     ;
  free((void *)N_pbot) ;
  free((void *)N_psiz) ;
  free((void *)N_pv)   ;
  return ;
}

/*----------------------------------------------------------------------------*/

#define WIN_SIM 1
#define WIN_POW 2
#define WIN_STP 3
static int opt_winner ;

/*----------------------------------------------------------------------------*/
/*! Chooses which optimization function(s) to call.
------------------------------------------------------------------------------*/

void generic_optimization
(
  vfp nmodel,             /* pointer to noise model */
  vfp smodel,             /* pointer to signal model */
  int r,                  /* number of parameters in the noise model */
  int p,                  /* number of parameters in the signal model */
  float * min_nconstr,    /* minimum parameter constraints for noise model */
  float * max_nconstr,    /* maximum parameter constraints for noise model */
  float * min_sconstr,    /* minimum parameter constraints for signal model */
  float * max_sconstr,    /* maximum parameter constraints for signal model */
  int nabs,               /* use absolute constraints for noise parameters */
  int ts_length,          /* length of time series array */
  float ** x_array,       /* independent variable matrix */
  float * ts_array,       /* observed time series */
  float * par_rdcd,       /* estimated parameters for the reduced model */
  float * parameters,     /* estimated parameters */
  float * sse             /* error sum of squares */
)
{
   float *powv , *simv , spow=1.e+33 , ssim=1.e+33 ;
   int dopow = (N_newuoa  > 0) ;
   int dosim = (N_newuoa == 2 || N_newuoa == 0) ;
   int stp=0 ;

   if( dopow && dosim ){
     powv = (float *)malloc(sizeof(float)*(r+p)) ;
     simv = (float *)malloc(sizeof(float)*(r+p)) ;
     memcpy(powv,parameters,sizeof(float)*(r+p)) ;
     memcpy(simv,parameters,sizeof(float)*(r+p)) ;
   } else {
     powv = simv = parameters ;
   }

   if( dosim ){                               /* Simplex from same start pt */
     simplex_optimization(nmodel, smodel, r, p,
                          min_nconstr, max_nconstr, min_sconstr, max_sconstr,
                          nabs, ts_length, x_array, ts_array, par_rdcd,
                          simv, &ssim );

     if( dopow ){
       float *qv = (float *)malloc(sizeof(float)*(r+p)) , qs=1.e+33 ;
       memcpy(qv,simv,sizeof(float)*(r+p)) ;
       set_newuoa_parm( 0.01 , 0.0009 , 666 ) ;   /* touchup with NEWUOA */
       newuoa_optimization(nmodel, smodel, r, p,
                           min_nconstr, max_nconstr, min_sconstr, max_sconstr,
                           nabs, ts_length, x_array, ts_array, par_rdcd,
                           qv, &qs );
       if( qs < ssim ){
         memcpy(simv,qv,sizeof(float)*(r+p)) ; ssim = qs ; stp = 1 ;
       }
       free((void*)qv) ;
     }
   }

   if( dopow ){                               /* NEWUOA from same start pt */
     set_newuoa_parm( 0.0 , 0.0 , 0 ) ;     /* default settings for NEWUOA */
     newuoa_optimization (nmodel, smodel, r, p,
                          min_nconstr, max_nconstr, min_sconstr, max_sconstr,
                          nabs, ts_length, x_array, ts_array, par_rdcd,
                          powv , &spow );
   }

   opt_winner = 0 ;
   if( dopow && dosim ){
     if( spow < ssim ) memcpy(parameters,powv,sizeof(float)*(r+p)) ;
     else              memcpy(parameters,simv,sizeof(float)*(r+p)) ;
     free((void *)simv); free((void *)powv);

     opt_winner = (spow < ssim) ? WIN_POW
                                : (stp) ? WIN_STP : WIN_SIM ;
   }
   *sse = (spow < ssim) ? spow : ssim ;
}
