/* 1 -1 {xr}=search(f,x0[,dx[,x1...xN]]), search for minimum of f(xr)

Uses the Nelder-Mead simplex algorithm to determine values for the N elements
of range xr which minimize the formula defined by cell or symbol f.

x0 specifies the initial value for xr.  Optionally, dx specifies step sizes
used to create the initial simplex, and x1...xN specifies the initial simplex
vertices.  If x1...xN are specified then dx is only used for restarts.  If dx
is not specified then a default step size of +10% in each direction is used.
x0, dx, and x1...xN may be expression lists or ranges.

xr must be a range, even in the one-dimensional case, e.g. {a0:a0}=search(...)
must be used instead instead of a0=search(...).  This allows the values of xr to
be accessed through the left search Node pointer so they can be set during the
simplex function evaluations.  Also, xr will not contain the search formula,
it will only be marked as depending on it, so xr can be initialized to values
to use for x0.
*/

//----------------------------------------

#define SEARCH_xSS sspp

// set xSS pointers
//
static void rf_search_xSS_callback( Data *d)
{
  *d->SEARCH_xSS++ = d->cp;

 // have to mark the cells EVALUATED so their indirect dependency
 // on the search() formula will not be flagged as cyclic
 //
  d->cp->state = EVALUATED;
}
 
//----------------------------------------

#define SEARCH_X     dp
#define SEARCH_STEP  dp2
#define SEARCH_COUNT i[0]
#define SEARCH_N     i[1]
#define SEARCH_2N    i[2]
#define SEARCH_MAX   i[3]

// set x0, dx, x1...xN values
//
static void rf_search_x0_callback( Data *d)
{
  ++d->SEARCH_COUNT;

  if( d->SEARCH_COUNT <= d->SEARCH_N) // getting x0
    *d->SEARCH_X++ = d->list_val;
  else if( d->SEARCH_COUNT <= d->SEARCH_2N) // getting dx
    *d->SEARCH_STEP++ = d->list_val;
  else if( d->SEARCH_COUNT <= d->SEARCH_MAX) // getting x1...xN
    *d->SEARCH_X++ = d->list_val;
  else
    fprintf( out, "search: too many argument values\n");
}
 
//----------------------------------------

double rf_search(const Node *e, const Cell *c)
{
  Data d; Node *xr = Left(e);

  if( debug) fprintf( out, "search: start\n");

  if( !depend)
    fprintf( out, "search: warning: dependency evaluation not set\n");

  if( check_range( xr, "search: return value")) return 0;

  Node *f = Right(e);

  if( !f || (f->type != CELL && f->type != ID))
  {
    fprintf( out, "search: first argument must be a cell or symbol\n");
    return 0;
  }

  if( f->type == CELL && check_cell( &f->u.c)) return 0;

  Node *x0 = f->next;

  if( !x0)
  {
    fprintf( out, "search: missing x0 argument\n");
    return 0;
  }

  int N = Rsize( xr);

  if( debug) fprintf( out, "search: creating simplex, N = %d\n", N);

  Simplex *s = simplex_create( N);

 // set simplex reference cell and function formula
 //
  if( f->type == CELL)
  {
    s->cSS = &f->u.c;
    s->fSS = ss[f->u.c.row][f->u.c.col].n;
  }
  else
  {
    s->cSS = &cell_00;
    s->fSS = f->u.s->n;
  }

 // initialize xSS with the N cell pointers for range xr
 //
  d.SEARCH_xSS = s->xSS;

  ss_traverse_range( &xr->u.r.start, &xr->u.r.end, rf_search_xSS_callback, &d);

 // initialize x0 and optionally dx, x1...xN
 //
  d.SEARCH_X = s->x[0];
  d.SEARCH_STEP = s->x[s->S];
  d.SEARCH_N = N;
  d.SEARCH_2N = 2*N;
  d.SEARCH_MAX = N*(2+N);
  d.SEARCH_COUNT = 0;

  if( debug) fprintf( out, "search: traversing argument list\n");

  ss_traverse_list( x0, c, rf_search_x0_callback, &d);

  if( d.SEARCH_COUNT != d.SEARCH_N  &&
      d.SEARCH_COUNT != d.SEARCH_2N &&
      d.SEARCH_COUNT != d.SEARCH_MAX)
  {
    fprintf( out, "search: wrong number of argument values\n");
    simplex_free( s);
    return 0;
  }

  int init_x = 0; // set if we have x1...xN

  if( d.SEARCH_COUNT == d.SEARCH_MAX) init_x = 1;
  else if( d.SEARCH_COUNT == d.SEARCH_N) // use default dx, 10% x0
  {
    double *x0 = s->x[0], *dx = s->x[s->S];
    for( int j = 0; j < N; ++j) dx[j] = (x0[j] == 0) ? 0.1 : 0.1*x0[j];
  }

  double r = simplex( s, init_x);

 // set final f value
 //
  if( f->type == CELL)
  {
    ss[f->u.c.row][f->u.c.col].val = r;
    ss[f->u.c.row][f->u.c.col].state = EVALUATED;
  }
  else
  {
    f->u.s->val = r;
    f->u.s->state = EVALUATED;
  }

  simplex_free( s);

  return r;
}